Fix drag and lift computations in TF problem

[jorgensd]

Interior facet integrals needs a volume tag to have a consistent orientation: https://github.com/KVSlab/turtleFSI/blob/e1b5350f5f105487bd77b811796f24e3b6075ea2/turtleFSI/problems/TF_fsi.py#L230-L233 ref https://fenicsproject.discourse.group/t/integrating-over-an-interior-surface/247/4?u=dokken

[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

[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.
https://github.com/KVSlab/turtleFSI/blob/b4df3739fabfa7b9592ba9511a8b7fa464e9860f/turtleFSI/monolithic.py#L108 to

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

https://github.com/KVSlab/turtleFSI/blob/b4df3739fabfa7b9592ba9511a8b7fa464e9860f/turtleFSI/problems/TF_fsi.py#L232-L233 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.

[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]

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