When using gmshio to read hexahedral mesh, an error occurs

As per the title, I have confirmed that there is only one type of mesh (hexahedral), but I still encounter an error when reading the .msh file:

from dolfinx.io import gmshio
proc = MPI.COMM_WORLD.rank
mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
Traceback (most recent call last):
  File "/home/dyfluid/Desktop/data_3D/convert.py", line 10, in <module>
    mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
                                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py", line 373, in read_from_msh
    msh = model_to_mesh(gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py", line 244, in model_to_mesh
    topologies = extract_topology_and_markers(model)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py", line 138, in extract_topology_and_markers
    assert len(entity_types) == 1
           ^^^^^^^^^^^^^^^^^^^^^^
AssertionError

Here is my .geo file:

Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};

Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 0.3} {Surface{1};Layers {3};Recombine;}

Physical Surface("top_bottom", 4) = {1,3};
Physical Surface("cylinder", 5) = {2};
Physical Volume("v", 6) = {1};

Mesh 2;
RecombineMesh;
Mesh.SubdivisionAlgorithm = 1;

and the visualization results from the Gmsh interface are shown below. According to the log message, it’s evident that all the face meshes are quadrilateral:


I hope to receive guidance, thank you!

Very first step could be to edit /home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py, adding a new line before line 138 with a print(entity_types) to see what entity_types contains.

According to your suggestion, I got the output below, and I’m not quite sure what it means. Could you help me analyze it.

[3]
[3]
[3]
[]

Which version of gmsh are you using? The file you provided does not seem to be a valid geo file for my version (4.12.2).

$ gmsh mesh.geo -2 -o mesh.msh
Info    : Running 'gmsh mesh.geo -2 -o mesh.msh' [Gmsh 4.12.2, 1 node, max. 1 thread]
Info    : Started on Tue Jul 30 09:20:53 2024
Info    : Reading 'mesh.geo'...
Error   : 'mesh.geo', line 1: Circle requires 3 points
Error   : Unknown curve 1
Error   : 'mesh.geo', line 4: Could not add plane surface
Warning : Skipping unknown surface 2 in physical surface 5
Info    : Meshing 1D...
Info    : Done meshing 1D (Wall 1.489e-05s, CPU 1.6e-05s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Warning : Mesh generation of surface 1 skipped: only 0 nodes on the boundary
Warning : Surface 1 consists of no elements
Info    : [ 50%] Meshing surface 3 (Extruded)
Info    : Done meshing 2D (Wall 0.000101788s, CPU 9.7e-05s)
Info    : 0 nodes 0 elements
Warning : ------------------------------
Warning : Mesh generation error summary
Warning :     2 warnings
Warning :     0 errors
Warning : Check the full log for details
Warning : ------------------------------
Info    : Recombining 2D mesh...
Info    : Blossom recombination completed (Wall 1.8739e-05s, CPU 1.8e-05s): 0 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = -nan, min Q = 1
Info    : Blossom recombination completed (Wall 1.529e-06s, CPU 1e-06s): 0 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = -nan, min Q = 1
Info    : Done recombining 2D mesh (Wall 4.3256e-05s, CPU 4.1e-05s)
Warning : Skipping unknown surface 2 in physical surface 5
Info    : Done reading 'mesh.geo'
Info    : Meshing 1D...
Info    : Done meshing 1D (Wall 3.665e-06s, CPU 4e-06s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Warning : Mesh generation of surface 1 skipped: only 0 nodes on the boundary
Warning : Surface 1 consists of no elements
Info    : [ 50%] Meshing surface 3 (Extruded)
Info    : Done meshing 2D (Wall 2.6817e-05s, CPU 2.5e-05s)
Info    : 0 nodes 0 elements
Warning : ------------------------------
Warning : Mesh generation error summary
Warning :     2 warnings
Warning :     0 errors
Warning : Check the full log for details
Warning : ------------------------------
Info    : Writing 'mesh.msh'...
Info    : Done writing 'mesh.msh'
Info    : Stopped on Tue Jul 30 09:20:53 2024 (From start: Wall 0.0534747s, CPU 0.076493s)
$ gmsh mesh.geo -3 -o mesh.msh
Info    : Running 'gmsh mesh.geo -3 -o mesh.msh' [Gmsh 4.12.2, 1 node, max. 1 thread]
Info    : Started on Tue Jul 30 09:21:34 2024
Info    : Reading 'mesh.geo'...
Error   : 'mesh.geo', line 1: Circle requires 3 points
Error   : Unknown curve 1
Error   : 'mesh.geo', line 4: Could not add plane surface
Warning : Skipping unknown surface 2 in physical surface 5
Info    : Meshing 1D...
Info    : Done meshing 1D (Wall 5.503e-06s, CPU 6e-06s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Warning : Mesh generation of surface 1 skipped: only 0 nodes on the boundary
Warning : Surface 1 consists of no elements
Info    : [ 50%] Meshing surface 3 (Extruded)
Info    : Done meshing 2D (Wall 4.4268e-05s, CPU 4.2e-05s)
Info    : 0 nodes 0 elements
Warning : ------------------------------
Warning : Mesh generation error summary
Warning :     2 warnings
Warning :     0 errors
Warning : Check the full log for details
Warning : ------------------------------
Info    : Recombining 2D mesh...
Info    : Blossom recombination completed (Wall 1.1907e-05s, CPU 1.1e-05s): 0 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = -nan, min Q = 1
Info    : Blossom recombination completed (Wall 6.78e-07s, CPU 1e-06s): 0 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = -nan, min Q = 1
Info    : Done recombining 2D mesh (Wall 2.446e-05s, CPU 2.3e-05s)
Warning : Skipping unknown surface 2 in physical surface 5
Info    : Done reading 'mesh.geo'
Info    : Meshing 1D...
Info    : Done meshing 1D (Wall 1.694e-06s, CPU 2e-06s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Warning : Mesh generation of surface 1 skipped: only 0 nodes on the boundary
Warning : Surface 1 consists of no elements
Info    : [ 50%] Meshing surface 3 (Extruded)
Info    : Done meshing 2D (Wall 1.2938e-05s, CPU 1.2e-05s)
Info    : Meshing 3D...
Info    : Meshing volume 1 (Extruded)
Error   : No elements in volume 1 
Info    : Done meshing 3D (Wall 3.6784e-05s, CPU 3.5e-05s)
Info    : 0 nodes 0 elements
Error   : ------------------------------
Error   : Mesh generation error summary
Error   :     2 warnings
Error   :     1 error
Error   : Check the full log for details
Error   : ------------------------------
Info    : Writing 'mesh.msh'...
Info    : Done writing 'mesh.msh'
Info    : Stopped on Tue Jul 30 09:21:34 2024 (From start: Wall 0.00126867s, CPU 0.078055s)

I am using version 4.10.5. I modified the msh file to eliminate the previous error. However, when I tried to read it again, MPI-related errors occurred:

from mpi4py import MPI
from dolfinx.io import gmshio

mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
Invalid rank, error stack:
internal_Issend(60788): MPI_Issend(buf=0x25f3e61, count=1, MPI_BYTE, 1, 1, comm=0x84000005, request=0x25f3e84) failed
internal_Issend(60749): Invalid rank has value 1 but must be nonnegative and less than 1
Abort(674890502) on node 0 (rank 0 in comm 416): application called MPI_Abort(comm=0x84000005, 674890502) - process 0

Below is my current .geo file:

// Gmsh project created on Tue Jul 30 10:39:54 2024
SetFactory("OpenCASCADE");

Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};

Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 0.3} {Surface{1};Layers {3};Recombine;}

Physical Surface("top_bottom", 4) = {1,3};
Physical Surface("cylinder", 5) = {2};
Physical Volume("v", 6) = {1};

Mesh.SubdivisionAlgorithm = 1;

Here is the visualization result in ParaView:
100de092ec876778ad21e4d31de2bbb

When I run this, I get

Info    : Reading 'mesh.msh'...
Info    : 9 entities
Info    : 260 nodes
Info    : 596 elements
Info    : Done reading 'mesh.msh'
Unknown cell type 6.
Traceback (most recent call last):
  File "/tmp/test.py", line 4, in <module>
    mesh, cell_markers, facet_markers = gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, gdim=3)
                                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/lib/python3.12/site-packages/dolfinx/io/gmshio.py", line 374, in read_from_msh
    msh = model_to_mesh(gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/lib/python3.12/site-packages/dolfinx/io/gmshio.py", line 291, in model_to_mesh
    ufl_domain = ufl_mesh(cell_id, gdim, dtype=dtype)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/lib/python3.12/site-packages/dolfinx/io/gmshio.py", line 75, in ufl_mesh
    raise e
  File "/lib/python3.12/site-packages/dolfinx/io/gmshio.py", line 72, in ufl_mesh
    shape, degree = _gmsh_to_cells[gmsh_cell]
                    ~~~~~~~~~~~~~~^^^^^^^^^^^
KeyError: 6

From Gmsh 4.13.1 you can see that 6 means 6-node prism, which not in the list at dolfinx/python/dolfinx/io/gmshio.py at 5d15d263f20a9d8ada051dd2ece0f35894d38d4a · FEniCS/dolfinx · GitHub

Please ask at Issues · gmsh / gmsh · GitLab how to change your geo file so that it uses one of the supported cell types (hex, I guess), and then report back the answer.

Sorry, this is my mistake: I just found that the .msh files generated from the same .geo file are different in different versions of Gmsh. When you use the new software to convert the .geo file I provided into a .msh file, the surface mesh shape is triangular, corresponding to the 3D mesh that you mentioned earlier with six vertices forming a prism. However, in the old version of the software, the same .geo file gave me the result as in my previous reply, which is a hexahedral mesh, but it produced an MPI error. Perhaps I should directly provide you with the .msh file (I tried just now, but I don’t know which cloud storage I should upload the file to), so we wouldn’t get different results due to the different versions of Gmsh.
So, my question is regarding the MPI-related error mentioned in my previous reply. Do you have any suggestions for this error? (I have searched for related posts, but there isn’t much content about this error, so I can’t find help)

Please sort out the gmsh part first with them. It’s unlikely that a version upgrade changes so dramatically (different cell types) the final mesh.

Sure, I will do my best to resolve the Gmsh issue. However, the same MPI problem still exists and has not been resolved, as mentioned in this post.
https://fenicsproject.discourse.group/t/read-meshtags-causing-invalid-rank-issues/14349?u=french_fries
I hope you can take a look at it.

As I said, please ask at Sign in · GitLab a relevant question could be why with a certain version you get a hex mesh, and with another a prism mesh. That is the gmsh support channel, which is surely more qualified than me to help in that respect. I won’t be replying any further until that part has been sorted.

The MPI error you link in that post is not necessarily related: MPI_Issend is a very low level function, which may be called in several places in the code, hence there is no guarantee that the two failures occur at exactly the same place.

I see, it looks like I indeed need to address the issue with the grid first. Thank you for your explanation.

Dear, I found that it is not a Gmsh version issue. In fact, the same geo file produces different results depending on how it is opened: for the mesh below, if I use gmsh test.geo -format msh2 -3 to convert it to an msh file, the result is as shown in Figure 1; however, if I open the geo file in the UI interface and save it as an msh file (which may require opening it twice), the result is as shown in Figure 2 (These are two different shapes of hexahedral meshes, perhaps you can reproduce this based on my description).

// Gmsh project created on Tue Jul 30 10:39:54 2024
SetFactory("OpenCASCADE");

Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};
Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 0.3} {Surface{1};Layers {3};Recombine;}
Physical Surface("top_bottom", 4) = {1,3};
Physical Surface("cylinder", 5) = {2};
Physical Volume("v", 6) = {1};

Mesh 2;
RecombineMesh;
Mesh.SubdivisionAlgorithm = 1;


The msh file generated in the latter case produces the following error when reading with gmshio.read_from_msh :

Traceback (most recent call last):
  File "/home/dyfluid/Desktop/data_3Dunit/convert.py", line 10, in <module>
    mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
                                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py", line 374, in read_from_msh
    msh = model_to_mesh(gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py", line 245, in model_to_mesh
    topologies = extract_topology_and_markers(model)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/gmshio.py", line 139, in extract_topology_and_markers
    assert len(entity_types) == 1
           ^^^^^^^^^^^^^^^^^^^^^^
AssertionError

As I’ve already stated, you must ask at the gmsh forum, since this is a gmsh issue. I am closing this.