Fix drag and lift computations in TF problem

[kvslab]

Fix drag and lift computations in TF problem

[jorgensd]

Here is an illustration of how it should be:

from dolfin import *
import numpy as np
parameters["ghost_mode"] = "shared_facet"

L = 0.2
H = 0.2
if MPI.comm_world.rank == 0:
    mesh = UnitSquareMesh(MPI.comm_self, 10, 10)

    class Domain(SubDomain):
        def inside(self, x, on_boundary):
            return abs(x[0] - 0.5) < L/2 + 100*DOLFIN_EPS and abs(x[1] - 0.5) < H/2 + 100*DOLFIN_EPS
    tdim = mesh.topology().dim()
    # Mark part of the domain
    cf = MeshFunction('size_t', mesh, mesh.topology().dim(), 1)

    Domain().mark(cf, 2)

    # Find interface between domains
    ff = MeshFunction("size_t", mesh, tdim-1, 1)
    mesh.init(tdim-1, tdim)
    f_to_c = mesh.topology()(tdim-1, tdim)
    for facet in facets(mesh):
        cells = f_to_c(facet.index())
        values = cf.array()[cells]
        if len(np.unique(values)) == 2:
            ff.array()[facet.index()] = 2
    with XDMFFile(MPI.comm_self, "ff.xdmf") as xdmf:
        xdmf.write(ff)
    with XDMFFile(MPI.comm_self, "cf.xdmf") as xdmf:
        xdmf.write(mesh)
        xdmf.write(cf)
MPI.comm_world.Barrier()
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("cf.xdmf") as xdmf:
    xdmf.read(mesh)
    xdmf.read(mvc)
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("ff.xdmf") as xdmf:
    xdmf.read(mvc)
ff = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# Intgral over a s

dS = Measure("dS", domain=mesh, subdomain_data=ff)
dx = Measure("dx", domain=mesh, subdomain_data=cf)

V = FunctionSpace(mesh, "DG", 0)
u = TrialFunction(V)
v = TestFunction(V)
ah = inner(u, v)*dx
Lh = inner(3, v)*dx(2) + inner(7, v)*dx(1)
uh = Function(V)
solve(ah == Lh, uh)

n = FacetNormal(mesh)
x = SpatialCoordinate(mesh)
integral = uh("+")*dS(2)
restriction = Constant(0)*dx
print("No restriction:", assemble(integral))
print("Restriction:", assemble(integral+restriction))
print("Exact:", 2*3*L + 2*3*H)

which yields:

No restriction: 3.9999999999999987
Restriction: 2.399999999999999
Exact: 2.4000000000000004

@keiyamamo

[jorgensd]

You can use Nanson's formula, ref:
https://en.wikipedia.org/wiki/Finite_strain_theory#Transformation_of_a_surface_and_volume_element

import dolfin

mesh = dolfin.UnitSquareMesh(10, 10)
n = dolfin.FacetNormal(mesh)

V = dolfin.VectorFunctionSpace(mesh, "DG", 2)
u = dolfin.TrialFunction(V)
v = dolfin.TestFunction(V)
a = dolfin.Constant(0)*dolfin.inner(u, v)*dolfin.dx+ dolfin.inner(u, v)*dolfin.ds
I = dolfin.Identity(len(u))
u = dolfin.Function(V, name="u")
u.interpolate(dolfin.Expression(("x[0]","sin(x[0])"), degree=1))
F = I+dolfin.grad(u)
J = dolfin.det(F)
n_new = J*dolfin.inv(F.T)*n
L = dolfin.inner(n_new, v)*dolfin.ds

A = dolfin.assemble(a, keep_diagonal=True)
A.ident_zeros()
b = dolfin.assemble(L)
nh = dolfin.Function(V, name="n")
dolfin.solve(A, nh.vector(), b)

with dolfin.XDMFFile("n_deformed.xdmf") as xdmf:
    dolfin.ALE.move(mesh, u)
    xdmf.write_checkpoint(u, "u", 0.0, append=False)
    xdmf.write_checkpoint(nh, "nh", 0.0, append=True)

[keiyamamo]

Hi @jorgensd

I looked into this but adding restriction actually made the results worse. I probably made some mistake in the implementation... while I’m looking into my implementation again, here is what I modified and please let me know if you find any mistakes.

turtleFSI/turtleFSI/monolithic.py

Line 108 in b4df373

dx = Measure("dx", subdomain_data=domains)

to

dx = Measure("dx", domain=mesh, subdomain_data=domains)

turtleFSI/turtleFSI/problems/TF_fsi.py

Lines 232 to 233 in b4df373

	Dr += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[0]*dS(5))
	Li += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[1]*dS(5))

to

Dr += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[0]*dS(5) + Constant(0)*dx)
Li += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[1]*dS(5) + Constant(0)*dx)
    

[keiyamamo]

Another thing about lift and drag computation is that normal vector, n, needs to be updated based on the flag deformation. Fixing this is not the top priority for now, but should be noted.

[keiyamamo]

Thank you @jorgensd ! I'll look into this when the time permits.