Hi! In the below MWE, I am creating a dented solid cylinder mesh in gmsh (version 4.11.1) using the extrude() command as given in the gmsh file given below. When importing it to FEniCS, it says there’s error in creating 2D surface mesh. However, I see no issue when I read the surface.xdmf file in Paraview!
Although there are several posts addressing errors reading 3D geometry imported from Gmsh in the community page, I found no suggestion that addresses the issue I am facing
SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.5, 0, 1.0};
//+
Point(2) = {0.0, -0.5, 0.0, 1.0};
//+
Point(3) = {0.0, 0.0, 0.5, 1.0};
//+
Point(4) = {0.0, 0.0, -0.5, 1.0};
//+
Point(5) = {0.0, 0.0, 0.0, 1.0};
//+
Circle(1) = {1, 5, 4};
//+
Circle(2) = {4, 5, 2};
//+
Circle(3) = {2, 5, 3};
//+
Circle(4) = {3, 5, 1};
//+
Point(6) = {50, 0.5, 0.0, 1.0};
//+
Point(7) = {50, -0.5, 0.0, 1.0};
//+
Point(8) = {50, 0.0, 0.5, 1.0};
//+
Point(9) = {50, 0.0, -0.5, 1.0};
//+
Point(10) = {50, 0.0, 0.0, 1.0};
//+
Circle(5) = {6, 10, 8};
//+
Circle(6) = {8, 10, 7};
//+
Circle(7) = {7, 10, 9};
//+
Circle(8) = {9, 10, 6};
//+
Point(11) = {70, 0.0, 0.0, 1.0};
//+
Point(12) = {70, 0.5, 0.0, 1.0};
//+
Point(13) = {70, -0.5, 0.0, 1.0};
//+
Point(14) = {70, 0.0, 0.5, 1.0};
//+
Point(15) = {70, 0.0, -0.5, 1.0};
//+
Circle(9) = {12, 11, 14};
//+
Circle(10) = {14, 11, 13};
//+
Circle(11) = {13, 11, 15};
//+
Circle(12) = {15, 11, 12};
//+
Point(16) = {120, 0.0, 0.0, 1.0};
//+
Point(17) = {120, 0.5, 0.0, 1.0};
//+
Point(18) = {120, -0.5, 0.0, 1.0};
//+
Point(19) = {120, 0.0, 0.5, 1.0};
//+
Point(20) = {120, 0.0, -0.5, 1.0};
//+
Circle(13) = {17, 16, 19};
//+
Circle(14) = {19, 16, 18};
//+
Circle(15) = {18, 16, 20};
//+
Circle(16) = {20, 16, 17};
//+
Point(21) = {60, 0.5-0.174551, 0.0, 1.0};
//+
Point(22) = {60.0, -0.5-0.174551, 0.0, 1.0};
//+
Point(23) = {60.0, 0.0-0.174551, 0.0, 1.0};
//+
Point(24) = {60.0, 0.0-0.174551, 0.5, 1.0};
//+
Point(25) = {60.0, 0.0-0.174551, -0.5, 1.0};
//+
Circle(17) = {21, 23, 24};
//+
Circle(18) = {24, 23, 22};
//+
Circle(19) = {22, 23, 25};
//+
Circle(20) = {25, 23, 21};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("leftsurface", 1) = {1, 4, 3, 2};
//+
Physical Surface("left", 6) = {1};
//+
Line(21) = {1, 6};
//+
Line(22) = {6, 21};
//+
Line(23) = {21, 12};
//+
Line(24) = {12, 17};
//+
Physical Curve("topline", 10) = {21, 22, 23, 24};
//+
Physical Curve("right", 5) = {13, 14, 15, 16};
//+
Curve Loop(2) = {13, 14, 15, 16};
//+
Plane Surface(2) = {2};
//+
Physical Surface("right", 7) = {2};
//+
Wire(3) = {21, 22, 23, 24};
Extrude { Surface{1}; } Using Wire {3}
//+
Physical Surface("boundary", 8) = {16, 12, 8, 4, 17, 13, 9, 5, 6, 18, 14, 10, 11, 15, 19, 7};
//+
Physical Volume("volume", 9) = {1};
Gmsh isn’t spitting out any error when reading it. Then I am meshing it with Automatic for 2D mesh and Delaunay for 3D mesh.
Then I import it to FEniCS with the below commands:
from dolfin import *
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)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
######################################################################################################################
msh = meshio.read("buckledcylinderl120w1.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")
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"topline = {assemble_edge(10)}")
print(f"leftcircle = {assemble_edge(1)}")
print(f"rightcircle = {assemble_edge(5)}")
print(f"leftsurface = {assemble(Constant(1.0)*dx_custom(6))}")
print(f"rightsurface = {assemble(Constant(1.0)*dx_custom(7))}")
print(f"boundarysurface = {assemble(Constant(1.0)*dx_custom(8))}")
This gives the following error:
Traceback (most recent call last):
File "/home/cylindermesh/mwe.py", line 35, in <module>
infile.read(mvc, "name_to_read")
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 never created meshes using the extrude() command before. But since Gmsh isn’t giving me any error when creating the 3D geometry and mesh, I guess the problem is with importing the mesh.
But I am not sure. Any help would be great!
Thanks in advance!