I am facing another mesh problem, which I used to solve by using an old mesh format with gmsh (msh2), but this time this won’t fix the problem.
I have the following mesh.geo:
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1e-2, 5e-3, 1e-4};
//+
Box(2) = {2.5e-3, 2.5e-3, 1e-4, 5e-3, 1e-4, 1e-4};
//+
BooleanFragments{ Volume{1}; Delete; }{ Volume{2}; Delete; }
//+
Physical Volume("the_box") = {1};
//+
Physical Volume("thin box") = {2};
//+
Physical Surface("end_left") = {7};
//+
Physical Surface("end_right") = {8};
//+
Physical Surface("interface") = {11};
//+
MeshSize {19, 23, 21, 17, 18, 20, 24, 22} = 5e-4;
//+
MeshSize {11, 9, 10, 12, 15, 13, 16, 14} = 1e-4;
//+
Physical Surface("bottom_surface") = {17};
from which I produce a corresponding file.msh using the command gmsh file.geo -3.
Then I use a dolfinx docker container (and Python3) to produce files.xdmf like so:
import meshio
msh = meshio.read("the_mesh.msh")
cells = np.vstack(np.array([cells.data for cells in msh.cells if cells.type == "tetra"]))
tetra_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tetra_mesh = meshio.Mesh(points=msh.points, cells=[("tetra", cells)], cell_data={"name_to_read": [tetra_data]})
facet_cells = np.vstack([cells.data for cells in msh.cells if cells.type == "triangle"])
facet_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
facet_mesh = meshio.Mesh(points=msh.points,
cells=[("triangle", facet_cells)],
cell_data={"name_to_read": [facet_data]})
# Write mesh
meshio.xdmf.write("mesh_test0.xdmf", tetra_mesh)
meshio.xdmf.write("facet_mesh_test0.xdmf", facet_mesh)
Up to there, everything seems to work fine, the code executes and the files are produced.
Then I use dolfin to read the mesh and compute the volumes and surfaces to make sure everything is okay, like so:
from dolfin import *
mesh = Mesh()
# Read mesh.
with XDMFFile("mesh_test0.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facet_mesh_test0.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh_test0.xdmf") as infile:
infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = Measure("dx", domain=mesh,subdomain_data=cf)
# External surface
ds = Measure("ds", domain=mesh, subdomain_data=mf)
# Internal surface
dS = Measure("dS", domain=mesh, subdomain_data=mf)
subdomains = cf
boundary_parts = mf
# Volumes
# 1 the_box
# 2 long thin box
# Surfaces
# 3 left end long box.
# 4 right end long box.
# 5 Interface between the boxes.
# 6 Bottom surface big box.
dx_1 = Measure("dx", domain=mesh,subdomain_id=1)
dx_2 = Measure("dx", domain=mesh,subdomain_id=2)
ds_bottom_1 = Measure("ds", domain=mesh,subdomain_id=6)
print('Expected volume big box:', 1e-2 * 5e-3 * 1e-4)
print('Computed volume: ',assemble(Constant(1)*dx_1))
print('Expected volume thin box:', 5e-3 * 1e-4 * 1e-4)
print('Computed volume: ',assemble(Constant(1)*dx_2))
print('Expected interface area:', 5e-3* 1e-4)
print('Computed surface area: ',assemble(Constant(1)*dS))
print('Expected bottom big box surface: ', 1e-2 * 5e-3)
print('Computed bottom big box surface: ', assemble(Constant(1)*ds_bottom_1))
Here’s the output:
Expected volume big box 1: 5e-09
Computed volume 1: 0.0
Expected volume thin box: 5.000000000000001e-11
Computed volume 2: 0.0
Expected interface area: 5.000000000000001e-07
Computed surface area: 0.00015914247160792985
Expected bottom big box surface: 5e-05
Computed bottom big box surface: 0.0
I am not sure where the problem lies. Any idea to retrieve the correct volumes and surfaces is welcome!