MPI Error: 'Invalid rank, error stack' when reading mesh

I/O
, ,
1

Hi everyone,
I am trying to get my first 3D project to work. I have this simple script to just read a mesh I created from a .geo file in Gmsh 4.13.1 (.geo file at the very end):

from mpi4py import MPI
import dolfinx.io.gmshio
mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("3D-chamber4.msh", MPI.COMM_WORLD, 0, gdim=3)
print(mesh, cell_tags, facet_tags)

However, I keep running into this error that I have read a couple of posts about but couldn’t figure it out:

Info : Done reading ‘3D-chamber4.msh’
Invalid rank, error stack:
internal_Issend(61644): MPI_Issend(buf=0x5632a62ba6e1, count=1, MPI_BYTE, 1, 1, comm=0x84000005, request=0x5632a62b9024) failed
internal_Issend(61605): Invalid rank has value 1 but must be nonnegative and less than 1
Abort(1010434822) on node 0 (rank 0 in comm 416): application called MPI_Abort(comm=0x84000005, 1010434822) - process 0

Trying to run the .py with mpirun -n 2 python3 read-mesh.py did not help. From other posts I read that it is usually a problem with unassigned physical surfaces. However, if I comment out all the physical surfaces in the .geo file, the read-mesh.py suddenly works and the output is

Info : Done reading ‘3D-chamber4.msh’
<dolfinx.mesh.Mesh object at 0x7fef52a31fd0> <dolfinx.mesh.MeshTags object at 0x7fef52a31e80> <dolfinx.mesh.MeshTags object at 0x7fef52bfafd0>

I’ve also tried to add Recombine; in the extrude command, but that just led to mixed element types and a different error message. The current .msh file only has Tetrahedas for 3D elements.

SetFactory("OpenCASCADE");

// Mesh sizing
h=1;
a=0.01;
Mesh.CharacteristicLengthMin = a/2*h;
Mesh.CharacteristicLengthMax = a*h;//+

Circle(1) = {0, 0, 0, 0.035, 0, 2*Pi};
Circle(2) = {0, 0, 0, 0.01, 0, 2*Pi};
Circle(3) = {0, 0, 0, 0.0127, 0, 2*Pi};

// Chamber inlet annulus
Curve Loop(1) = {1};
Curve Loop(2) = {3};
Plane Surface(1) = {1, 2};

// Fuel inlet
Curve Loop(3) = {3};
Curve Loop(4) = {2};
Plane Surface(2) = {3, 4};
Physical Surface(12) = {2};

// Lox inlet
Curve Loop(5) = {2};
Plane Surface(3) = {5};
Physical Surface(11) = {3};

//Extrude
Extrude {0, 0, 0.2} {
  Curve{1}; Layers{10};
}

// Outlet
Curve Loop(7) = {5};
Plane Surface(5) = {7};
Physical Surface(13) = {5};

// Walls
Physical Surface(14) = {4,1};

// Volume
Surface Loop(1) = {4, 5, 1, 2, 3};
Volume(1) = {1};
Physical Volume(15) = {1};
2

This is due to the fact that you are not tagging the correct surfaces in your geo file.
This is a common mistake when using the GMSH geo format, rather than their C++ or Python API.

This can be illustrated with:

import meshio
import numpy as np

mesh = meshio.read("mesh.msh")

tetrahedrons = mesh.cells_dict["tetra"]
triangles = mesh.cells_dict["triangle"]

tetra_vertices = np.unique(tetrahedrons.flatten())
tri_vertices = np.unique(triangles.flatten())
print(np.isin(tri_vertices, tetra_vertices))

for your mesh, you will observe that many of the vertices in the triangles are not part of any tetrahedra.

You should use coherence to merge duplicate nodes: