Hello,
I’ve been trying to add some force on an inner surface. However, the facets seems to have incoherent orientation.
According to Integrating over an interior surface - #4 by MiroK using different mesh domains on both sides of the surface should be sufficient to get a coherent orientation. Let me introduce a MWE. The mesh is generated with gmsh and consist of two blocks of same height, one inserted in the other. I have one tag for each: tag_silicone and tag_homogeneous. I will also tag each interface surface separately. Here is the code:
import gmsh
gmsh.initialize()
tag_frame = 0
gmsh.model.occ.addBox(-1, -1, -1, 10, 10, 3, tag=tag_frame)
gmsh.model.occ.addBox(0, 0, -1, 8, 8, 3, tag=tag_frame+1)
tag_silicone = 2
tag_homogeneous = 3
gmsh.model.occ.cut([(3, tag_frame)], [(3, tag_frame+1)],tag=tag_silicone)
gmsh.model.occ.addBox(0, 0, -1, 8, 8, 3, tag=tag_homogeneous)
gmsh.model.occ.fragment([(3,tag_silicone)],[(3,tag_homogeneous)])
gmsh.model.occ.synchronize()
inner_marker_1 = 10
inner_marker_2 = 11
inner_marker_3 = 12
inner_marker_4 = 13
inner_interfaces_1=[]
inner_interfaces_2=[]
inner_interfaces_3=[]
inner_interfaces_4=[]
for surface in gmsh.model.getEntities(dim=2):
com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
if (abs(com[2]-0.5)<1.e-6):
if (abs(com[0])<1.e-6): inner_interfaces_1.append(surface[1])
elif (abs(com[0]-8)<1.e-6): inner_interfaces_2.append(surface[1])
elif (abs(com[1])<1.e-6): inner_interfaces_3.append(surface[1])
elif (abs(com[1]-8)<1.e-6): inner_interfaces_4.append(surface[1])
gmsh.model.addPhysicalGroup(2, inner_interfaces_1, inner_marker_1)
gmsh.model.setPhysicalName(2, inner_marker_1, "neutral interface")
gmsh.model.addPhysicalGroup(2, inner_interfaces_2, inner_marker_2)
gmsh.model.setPhysicalName(2, inner_marker_2, "neutral interface")
gmsh.model.addPhysicalGroup(2, inner_interfaces_3, inner_marker_3)
gmsh.model.setPhysicalName(2, inner_marker_3, "neutral interface")
gmsh.model.addPhysicalGroup(2, inner_interfaces_4, inner_marker_4)
gmsh.model.setPhysicalName(2, inner_marker_4, "neutral interface")
volumes = gmsh.model.getEntities(dim=3)
silicon_marker = 11
homogeneous_marker = 12
gmsh.model.addPhysicalGroup(3, [tag_silicone], silicon_marker)
gmsh.model.setPhysicalName(3, silicon_marker, "Silicone")
gmsh.model.addPhysicalGroup(3, [tag_homogeneous], homogeneous_marker)
gmsh.model.setPhysicalName(3, homogeneous_marker, "HomogeneousMaterial")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.15)
gmsh.model.mesh.generate(3)
gmsh.write("Box.msh")
gmsh.finalize()
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
msh = meshio.read("Box.msh")
tetra_mesh = create_mesh(msh, "tetra")
triangle_mesh = create_mesh(msh, "triangle")
meshio.write("Box.xdmf", tetra_mesh)
meshio.write("BoxT.xdmf", triangle_mesh)
Now I will compute the surface of each interface surface of 4 different manners.
- Integrate 1dS over the corresponding domain
- Integrate n\cdot n_0 dS where n_0 is the reference normal, and n the computed normal
- Integrate n\cdot n_0 dS + 0 dx i.e. try the trick to add dx with subdomains to get coherent orientation of facets
- Integrate (n\cdot n_0)^2 dS
Note: In practice, I don’t need the surface, but I have to use the normal.
The resulting script is
from weakref import ref
import numpy as np
import ufl
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh
from dolfinx.io import XDMFFile
file_name = "Box.xdmf"
file_nameT = "BoxT.xdmf"
with XDMFFile(MPI.COMM_WORLD, file_name, "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(domain, name = "Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, file_nameT, "r") as xdmf:
ft = xdmf.read_meshtags(domain, name = "Grid")
V = fem.VectorFunctionSpace(domain, ("CG", 2))
silicon_marker = 11
homogeneous_marker = 12
silicon_tets = ct.indices[ct.values==silicon_marker]
silicon_dofs = fem.locate_dofs_topological(V, 3, silicon_tets)
homogeneous_tets = ct.indices[ct.values==homogeneous_marker]
homogeneous_dofs = fem.locate_dofs_topological(V, 3, homogeneous_tets)
marked_tets = np.hstack([silicon_tets, homogeneous_tets])
marked_values = np.hstack([np.full(len(silicon_tets), 1, dtype = np.int32), np.full(len(homogeneous_tets), 2, dtype = np.int32)])
sorted_tets = np.argsort(marked_tets)
tet_tag = mesh.meshtags(domain, 3, marked_tets[sorted_tets], marked_values[sorted_tets])
inner_marker_1 = 10
inner_marker_2 = 11
inner_marker_3 = 12
inner_marker_4 = 13
pressure_facets_1 = ft.indices[ft.values==inner_marker_1]
pressure_facets_2 = ft.indices[ft.values==inner_marker_2]
pressure_facets_3 = ft.indices[ft.values==inner_marker_3]
pressure_facets_4 = ft.indices[ft.values==inner_marker_4]
marked_facets = np.hstack([pressure_facets_1, pressure_facets_2, pressure_facets_3, pressure_facets_4])
marked_values = np.hstack([
np.full(len(pressure_facets_1), 0, dtype = np.int32),
np.full(len(pressure_facets_2), 1, dtype = np.int32),
np.full(len(pressure_facets_3), 2, dtype = np.int32),
np.full(len(pressure_facets_4), 3, dtype = np.int32),])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, 2, marked_facets[sorted_facets], marked_values[sorted_facets])
metadata = {"quadrature_degree": 4}
dS = ufl.Measure('dS', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata, subdomain_data=tet_tag)
u = fem.Function(V)
I = ufl.variable(ufl.Identity(3))
F = ufl.variable(I+ufl.grad(u))
J = ufl.variable(ufl.det(F))
N = ufl.FacetNormal(domain)
Finv = ufl.inv(F)
n = J*Finv.T*N
vec_x = fem.Constant(domain,PETSc.ScalarType((1.0,0,0)))
vec_y = fem.Constant(domain,PETSc.ScalarType((0.0,1.0,0)))
vecs=[-vec_x,vec_x,-vec_y,vec_y]
reference = 8.0*3.0
for i in range(0,4):
SurfaceForm_inner = fem.form(fem.Constant(domain, PETSc.ScalarType(1.))*dS(i))
surface_inner = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner), op=MPI.SUM)
SurfaceForm_inner_2 = fem.form(ufl.inner(n('-'),vecs[i])*dS(i))
surface_inner_2 = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner_2), op=MPI.SUM)
SurfaceForm_inner_3 = fem.form(ufl.inner(n('-'),vecs[i])*dS(i)+ fem.Constant(domain, PETSc.ScalarType(0.0))*dx)
surface_inner_3 = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner_3), op=MPI.SUM)
SurfaceForm_inner_4 = fem.form(ufl.inner(n('-'),vecs[i])**2*dS(i))
surface_inner_4 = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner_4), op=MPI.SUM)
print(f"Checking pressure on surface {i}: reference={reference}, {surface_inner} {surface_inner_2} {surface_inner_3} {surface_inner_4}")
Each output line should contains 4 output values, first and last should be coherent as they do not depends on the facets orientation. A correct last value should indicate the n is normal to the facet. If I understood the trick to get coherent orientation, third value should be correct.
However my output is
Checking pressure on surface 0:
reference=24.0, 23.999999999999897 23.80972755772505 23.80972755772505 23.999999999999897
Checking pressure on surface 1:
reference=24.0, 23.99999999999983 -18.501788312810014 -18.501788312810014 23.99999999999983
Checking pressure on surface 2:
reference=24.0, 23.99999999999986 23.736066937512884 23.736066937512884 23.99999999999986
Checking pressure on surface 3:
reference=24.0, 24.0 -18.570081136289787 -18.570081136289787 24.0
Any idea what am I doing wrong ?
Best