Hello,
I have a 3D mesh containing an inner surface, and I want to add some force on it. By inner surface, I mean a surface between two blocks.
I know how to do that with external forces.
Let me give a MWE: first generate some example with gmsh, and get the xdmf files with meshio.
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 = 2
xmin_marker = 3
xmin_interfaces=[]
inner_interfaces=[]
for surface in gmsh.model.getEntities(dim=2):
com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
if (abs(com[0]+1)<1.e-6):
xmin_interfaces.append(surface[1])
elif (abs(com[2]-0.5)<1.e-6):
if (abs(com[0])<1.e-6): inner_interfaces.append(surface[1])
elif (abs(com[0]-8)<1.e-6): inner_interfaces.append(surface[1])
elif (abs(com[1])<1.e-6): inner_interfaces.append(surface[1])
elif (abs(com[1]-8)<1.e-6): inner_interfaces.append(surface[1])
print(f"Nb of left surfaces: {len(xmin_interfaces)}")
print(f"Nb of inner surfaces: {len(inner_interfaces)}")
gmsh.model.addPhysicalGroup(2, inner_interfaces, inner_marker)
gmsh.model.setPhysicalName(2, inner_marker, "neutral interface")
gmsh.model.addPhysicalGroup(2, xmin_interfaces, xmin_marker)
gmsh.model.setPhysicalName(2, xmin_marker, "xmin 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 if I open BoxT.xdmf with paraview, I can seem my interfaces, everything seems good.
So let’s load this with fenicsX:
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")
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))
inner_marker = 2
xmin_marker = 3
attach_facets = ft.indices[ft.values==xmin_marker]
attach_dofs = fem.locate_dofs_topological(V, 3, attach_facets)
pressure_facets = ft.indices[ft.values==inner_marker]
marked_facets = np.hstack([attach_facets, pressure_facets])
marked_values = np.hstack([
np.full(len(attach_facets), 1, dtype = np.int32),
np.full(len(pressure_facets), 2, 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)
SurfaceForm_inner = fem.form(fem.Constant(domain, PETSc.ScalarType(1.))*ds(2))
surface_inner = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner), op=MPI.SUM)
SurfaceForm_attach = fem.form(fem.Constant(domain, PETSc.ScalarType(1.))*ds(1))
surface_attach = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_attach), op=MPI.SUM)
print(f"Checking attach surface={surface_attach}")
print(f"Checking pressure surface={surface_inner}")
The output is:
Checking attach surface=30.00000000000009
Checking pressure surface=0.0
So something is wrong, the outer surface evaluates to correct value, but inner surface gets zero . I checked that pressure_facets was not empty. Just taking one of the inner side instead also leads to zero value.
One other strange thing is that changing mesh characteristic lengths
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.15)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5)
leads to PETSC ERROR: Caught signal number 11 SEGV.
Would someone have a hint how to correct this, so that I can work with inner surfaces?
Thanks