I created a 3D cylindrical mesh in Gmsh (Version 4.11.1). I exported it as “Version 4 ASCII”. I did not check the “Save all elements” and “Save parametric coordinates” options while exporting. I usually export all my 3D meshes like that and all of them worked good so far.
I created the cylindrical mesh by extruding an annulus as given in the following Gmsh script:
// Gmsh project created on Mon Feb 17 01:33:47 2025
SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.0, 0.0, 1.0};
//+
Point(2) = {0.0, 0.5, 0.0, 1.0};
//+
Point(3) = {0.0, -0.5, 0.0, 1.0};
//+
Point(4) = {0.0, 0.0, 0.5, 1.0};
//+
Point(5) = {0.0, 0.0, -0.5, 1.0};
//+
Point(6) = {0.0, 1.0, 0.0, 1.0};
//+
Point(7) = {0.0, -1.0, 0.0, 1.0};
//+
Point(8) = {0.0, 0.0, 1.0, 1.0};
//+
Point(9) = {0.0, 0.0, -1.0, 1.0};
//+
Circle(1) = {3, 1, 4};
//+
Circle(2) = {4, 1, 2};
//+
Circle(3) = {2, 1, 5};
//+
Circle(4) = {5, 1, 3};
//+
Circle(5) = {7, 1, 8};
//+
Circle(6) = {8, 1, 6};
//+
Circle(7) = {6, 1, 9};
//+
Circle(8) = {9, 1, 7};
//+
Curve Loop(1) = {6, 7, 8, 5};
//+
Curve Loop(2) = {2, 3, 4, 1};
//+
Plane Surface(1) = {1, 2};
//+
Extrude {10, 0, 0} {
Curve{7}; Curve{6}; Curve{5}; Curve{8}; Curve{4}; Curve{1}; Curve{2}; Curve{3}; Surface{1};
}
//+
Physical Curve("leftinner", 1) = {2, 1, 4, 3};
//+
Physical Curve("leftouter", 2) = {8, 5, 6, 7};
//+
Physical Curve("rightouter", 4) = {15, 16, 11, 13};
//+
Physical Curve("rightinner", 3) = {21, 19, 24, 23};
//+
Physical Surface("rightcellsurface", 5) = {18};
//+
Physical Surface("leftcellsurface", 6) = {1};
//+
Physical Surface("outercellsurface", 8) = {3, 2, 5, 4};
//+
Physical Surface("innercellsurface", 7) = {8, 7, 6, 9};
//+
Physical Volume("cellvolume", 9) = {1};
I imported it to FEniCS (legacy version and DOLFIN version: 2019.2.0.13.dev0), to calculate area, volume and length of curves, as follows:
from dolfin import *
import meshio
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
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)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
msh = meshio.read("cylindermwe.msh")
tetra_mesh = create_mesh(msh, "tetra", True)
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("surface.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("surface.xdmf") as infile:
infile.read(mvc, "name_to_read") # This line gives the error
sf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
def assemble_edge(tag):
mv = MeshView.create(mf, tag)
return assemble(Constant(1.0)*dx(domain=mv)) # Line measure
dx_custom = Measure("ds", domain=mesh, subdomain_data=sf) # Area measure
dv_custom = Measure("dx", domain=mesh, subdomain_data=cf) # Volume measure
print(f"Volume = {assemble(Constant(1.0)*dv_custom(9))}")
print(f"leftouterboundary = {assemble_edge(2)}")
print(f"leftinnerboundary = {assemble_edge(1)}")
print(f"leftcellsurface = {assemble(Constant(1.0)*dx_custom(6))}")
print(f"rightouterboundary = {assemble_edge(4)}")
print(f"rightinnerboundary = {assemble_edge(3)}")
print(f"rightcellsurface = {assemble(Constant(1.0)*dx_custom(5))}")
print(f"outercellsurface = {assemble(Constant(1.0)*dx_custom(8))}")
print(f"innercellsurface = {assemble(Constant(1.0)*dx_custom(7))}")
This generates the error
line 39, in <module>
infile.read(mvc, "name_to_read") # This line gives the error
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://fenicsproject.discourse.group/
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to find entity in map.
*** Reason: Error reading MeshValueCollection.
*** Where: This error was encountered inside HDF5File.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------
I have previously imported and read 3D meshes successfully using the same code. I think the problem comes from extruding the surfaces and lines. I checked some of the solutions in this forum and tried some of them. But none of them worked.
I am not sure if I am creating the mesh wrong or am I importing it incorrectly.
Please advise.
Thanks in advance!
