New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fix drag and lift computations in TF problem #40
Comments
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 |
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
to
turtleFSI/turtleFSI/problems/TF_fsi.py Lines 232 to 233 in b4df373
to
|
Another thing about lift and drag computation is that normal vector, |
You can use Nanson's formula, ref: 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) |
Thank you @jorgensd ! I'll look into this when the time permits. |
Interior facet integrals needs a volume tag to have a consistent orientation:
turtleFSI/turtleFSI/problems/TF_fsi.py
Lines 230 to 233 in e1b5350
ref
https://fenicsproject.discourse.group/t/integrating-over-an-interior-surface/247/4?u=dokken
The text was updated successfully, but these errors were encountered: