Assertion Error while reading mesh with gmshio.read_from_msh: What's wrong here?

Errors
1

Dear Community,

I got an AssertionError while executing this line of code:

msh, cell_tags, facet_tags = gmshio.read_from_msh("sinerot_1_1merged_final.msh", MPI.COMM_WORLD,0, gdim=3)

The weird thing is I already read many meshs successfully with the same code. I can’t explain why this one isn’t working. Am I missing something? Did I something wrong during mesh generation?

The full error message reads:

Warning : Gmsh has aleady been initialized
Info    : Reading 'sinerot_1_1merged_final.msh'...
Info    : 496 entities
Info    : 59854 nodes
Info    : 362949 elements
Info    : Done reading 'sinerot_1_1merged_final.msh'                       
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[24], line 1
----> 1 msh, cell_tags, facet_tags = gmshio.read_from_msh("sinerot_1_1merged_final.msh", MPI.COMM_WORLD,0, gdim=3)

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/io/gmshio.py:373, in read_from_msh(filename, comm, rank, gdim, partitioner)
    371 gmsh.model.add("Mesh from file")
    372 gmsh.merge(filename)
--> 373 msh = model_to_mesh(gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner)
    374 gmsh.finalize()
    375 return msh

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/io/gmshio.py:243, in model_to_mesh(model, comm, rank, gdim, partitioner, dtype)
    241 assert model is not None, "Gmsh model is None on rank responsible for mesh creation."
    242 # Get mesh geometry and mesh topology for each element
--> 243 x = extract_geometry(model)
    244 topologies = extract_topology_and_markers(model)
    246 # Extract Gmsh cell id, dimension of cell and number of nodes to
    247 # cell for each

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/io/gmshio.py:199, in extract_geometry(model, name)
    195 # In some cases, Gmsh does not return the points in the same
    196 # order as their unique node index. We therefore sort nodes in
    197 # geometry according to the unique index
    198 perm_sort = np.argsort(indices)
--> 199 assert np.all(indices[perm_sort] == np.arange(len(indices)))
    200 return points[perm_sort]

AssertionError:

Thank you for any help!

2

This usually means that there are gaps in the point ids given by Gmsh. You would have to provide your mesh generation script

3

Thanks for the advice. I found one gap in my point ids. Unfortunately, the error remains and I couldn’t find another one. Here’s the code for the mesh:

period = 1
resolution = 30

hoehe = 50/2*0.001
durchmesser = 102/2*0.001
k = (resolution*4)

zeiger_theta = np.linspace(0, period*2*np.pi, k)
theta = np.cos(zeiger_theta)
plane = hoehe*theta
x_data = np.linspace(0, durchmesser, k)

plane = np.concatenate((plane, -np.array([hoehe*1.2])))
plane = np.concatenate((plane, -np.array([hoehe*1.2])))

x_data = np.concatenate((x_data, np.array([x_data[len(x_data)-1]])))
x_data = np.concatenate((x_data, np.array([0])))

if gmsh.isInitialized(): 
    gmsh.finalize()
gmsh.initialize()

modelName = "sinerot_1_"+str(period)+"merged_final"

gmsh.model.add(modelName)

lc = 0.005

gmsh.option.setString("Geometry.OCCTargetUnit", "M")

gmsh_surface_tag = 1
gmsh_curve_tag = 1
surface_array = []
surface_array_2 = []

for i in range(1, len(plane)+1): 
    gmsh.model.occ.addPoint(x_data[i-1],plane[i-1],0, lc, i)

plane[-1:][0] = 1
plane[-2:][0] = 1

gmsh.model.occ.addPoint(x_data[-1:][0],plane[-1:][0],0, lc, len(plane)+1)
gmsh.model.occ.addPoint(x_data[-2:][0],plane[-2:][0],0, lc, len(plane)+2)

for i in range(1, len(plane)): 
    gmsh.model.occ.addLine(i, i+1 ,gmsh_surface_tag)
    surface_array.append(gmsh_surface_tag)
    gmsh_surface_tag+=1
    if i == len(plane)-1: 
        gmsh.model.occ.addLine(i+1, 1 ,gmsh_surface_tag)
        surface_array.append(gmsh_surface_tag)
        gmsh_surface_tag+=1

gmsh.model.occ.addCurveLoop(surface_array, gmsh_curve_tag)
gmsh.model.occ.addPlaneSurface([gmsh_curve_tag], gmsh_curve_tag)
gmsh_curve_tag+=1

gmsh.model.occ.addLine(123, 124 ,gmsh_surface_tag)
surface_array_2.append(gmsh_surface_tag)
gmsh_surface_tag+=1

gmsh.model.occ.addLine(124, 120 ,gmsh_surface_tag)
surface_array_2.append(gmsh_surface_tag)
gmsh_surface_tag+=1

for i in reversed(range(1, 120)): 
    surface_array_2.append(i)

gmsh.model.occ.addLine(1, 123 ,gmsh_surface_tag)
surface_array_2.append(gmsh_surface_tag)
gmsh_surface_tag+=1

gmsh.model.occ.addCurveLoop(surface_array_2, gmsh_curve_tag)
gmsh.model.occ.addPlaneSurface([gmsh_curve_tag], gmsh_curve_tag)
gmsh_curve_tag+=1

gmsh.model.occ.revolve([(2,1), (2, 2)], 0, 25, 0, 0, 1, 0, -2*np.pi)

gmsh.model.occ.remove([(2, 1), (2,2)])

gmsh.model.occ.fragment([(3, 1)], [(3, 2)])

gmsh.model.occ.synchronize()

volumes = gmsh.model.getEntities(dim=3)

gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], 1)
gmsh.model.setPhysicalName(volumes[0][0], 1, "Geometry volume")

gmsh.model.addPhysicalGroup(volumes[1][0], [volumes[1][1]], 2)
gmsh.model.setPhysicalName(volumes[1][0], 2, "Tube volume")

gmsh.model.addPhysicalGroup(2, [125], 1)
gmsh.model.setPhysicalName(2, 1, "Tube Boundary")

gmsh.model.addPhysicalGroup(2, [124], 2)
gmsh.model.setPhysicalName(2, 2, "Velocity Plain")

gmsh.model.addPhysicalGroup(2, [123], 3)
gmsh.model.setPhysicalName(2, 3, "Bottom Boundary")

gmsh.model.addPhysicalGroup(2, [122], 4)
gmsh.model.setPhysicalName(2, 4, "Geometry Boundary sides")

for i in range(3, 122): 
    gmsh.model.addPhysicalGroup(2, [i], i+2)
    gmsh.model.setPhysicalName(2, i+2, "Geometry Cone "+str(i-2))

gmsh.model.mesh.generate(3)

# GUI
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.write(modelName+".msh")
#gmsh.write(modelName+".stl")

gmsh.finalize()
4

I cannot reproduce the error using:

mesh, ct, ft = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)

but I can reproduce it with:

mesh, ct, ft = dolfinx.io.gmshio.read_from_msh(modelName+".msh", MPI.COMM_WORLD, 0, gdim=3)

indicating that something is weirdly stored in the msh file.

Adding:

gmsh.model.mesh.generate(3)


indices, points, _ = gmsh.model.mesh.getNodes()
points = points.reshape(-1, 3)
indices -= 1
perm_sort = np.argsort(indices)
assert np.all(indices[perm_sort] == np.arange(len(indices)))



gmsh.model.add("Mesh from file")
gmsh.merge(modelName+".msh")



indices2, points2, _ = gmsh.model.mesh.getNodes()
points2 = points2.reshape(-1, 3)
indices2 -= 1

perm_sort2 = np.argsort(indices2)
assert np.all(indices2[perm_sort2] == np.arange(len(indices2)))

you will see that points and points2 have different shapes.

It seems to relate to the addition of the second volume:

gmsh.model.addPhysicalGroup(volumes[1][0], [volumes[1][1]], 2)

but I cannot tell why that changes the gmsh/msh model. I think this is better asked in the GMSH forum (now that I’ve removed DOLFINx from the code to illustrate the inconsistency within gmsh.

5

Thank you! I’m going to ask this in the Gmsh forum but until then the model_to_mesh method is a good workaround.