Problem with xdmf.read_mesh

General
1

Hi,

I am running a script to convert a .msh file to a .xdmf file, using the following part of the script:

if __name__ == "__main__":
    path_msh = Path("../meshes/patient9.msh")

    msh = meshio.read(path_msh)
    
    mesh_interior = create_mesh(msh, "tetra")
    mesh_facets = create_mesh(msh, "triangle")
    # store intermediate meshes in XDMF format
    
    path_tmp_tets = path_msh.with_name(path_msh.stem + "_tets.xdmf")
    path_tmp_tri = path_msh.with_name(path_msh.stem + "_tri.xdmf")
    
    meshio.write(path_tmp_tets, mesh_interior, file_format="xdmf")
    meshio.write(path_tmp_tri, mesh_facets, file_format="xdmf")
    
    # read back in with dolfinx to obtain the mesh in the correct format
    with io.XDMFFile(MPI.COMM_WORLD, path_tmp_tets, "r") as xdmf:
        print("test1")
        mesh = xdmf.read_mesh(name="Grid")
        print("test2")

However, my script runs indefinitely en gets stuck at the mesh = xdmf.read_mesh(name=“Grid”) line: it prints test1 but not test2. do you have any idea why this happens?

2

Please add the output of the plaintext xdmf file in 3x` encapsulation, ie.

```
Output of xdmf
```
3 
<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="37093 3" Format="HDF" Precision="8">patient9cgs_tets.h5:/data0</DataItem></Geometry><Topology TopologyType="Tetrahedron" NumberOfElements="185097" NodesPerElement="4"><DataItem DataType="Int" Dimensions="185097 4" Format="HDF" Precision="4">patient9cgs_tets.h5:/data1</DataItem></Topology><Attribute Name="markers" AttributeType="Scalar" Center="Cell"><DataItem DataType="Int" Dimensions="185097" Format="HDF" Precision="4">patient9cgs_tets.h5:/data2</DataItem></Attribute></Grid></Domain></Xdmf>
4

To me this mesh looks perfectly valid (as xdmf).
Are you running this in serial or parallel?

Could you provide the msh file that you are using?

5

I am running it in serial. The msh file is a rather large file: with 200.000 elements and 3700 nodes, so i cannot post it here, unfortunately. maybe a small part is useful?:

$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
37093
1 0.2698241770267487 -0.2756259739398956 -0.5569514036178589
2 0.2884699404239655 -0.2480147480964661 -0.5331186056137085
3 0.294191300868988 -0.2771624624729156 -0.5462889671325684
4 0.2622835636138916 -0.2926439344882965 -0.5369392037391663
5 0.2864086031913757 -0.2938072383403778 -0.5268594622612
6 0.3178472518920898 -0.2793562710285187 -0.5350016951560974
7 0.3099758923053741 -0.2957432270050049 -0.5159472823143005
8 0.2445784211158752 -0.2754203379154205 -0.5661513209342957
9 0.2649220824241638 -0.2454865276813507 -0.5443495512008667
6

A small part would not be useful, as I would have to try to reproduce your error message.

How is the mesh generated? Is is created with GMSH or is it just an intermediate format?

7

The mesh is created with GMSH. by importing as an stl file adding a volume as an physical entity, and the adding a volume mesh an exporting it as an msh file. I do now think of the fact that I exported it as a version 2 ASCII instead of the a version 4 ASCII, could that have anything to do with it?

8

what does dolfinx.io.gmshio.read_from_read_from_msh return?

9

for the following line:

mesh, cell_markers, facet_markers = gmshio.read_from_msh("../meshes/patient9cgs.msh", MPI.COMM_WORLD, gdim=3)

The following returned.

Info    : Reading '../meshes/patient9cgs.msh'...
Info    : 37093 nodes
Info    : 213865 elements
Info    : Done reading '../meshes/patient9cgs.msh'                         
Traceback (most recent call last):
  File "/home3/s3234266/nsx/examples/scripts/readmsh.py", line 4, in <module>
    mesh, cell_markers, facet_markers = gmshio.read_from_msh("../meshes/patient9cgs.msh", MPI.COMM_WORLD, gdim=3)
  File "/home3/s3234266/spack/var/spack/environments/fenicsx/.spack-env/view/lib/python3.10/site-packages/dolfinx/io/gmshio.py", line 311, in read_from_msh
    msh = model_to_mesh(gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner)
  File "/home3/s3234266/spack/var/spack/environments/fenicsx/.spack-env/view/lib/python3.10/site-packages/dolfinx/io/gmshio.py", line 221, in model_to_mesh
    cell_id = cell_information[perm_sort[-1]]["id"]
IndexError: index -1 is out of bounds for axis 0 with size 0
10

So this error message means that you are not using phyiscal entity grouping in Gmsh, Which is adviced to avoid Getting duplicate nodes and cells in your mesh.

11

ok, thanks!
I managed to get rid of the indexerror above by adding a volume and a surface as a physical entity. However, when running the original code I still get the problem that it gets stuck at the mesh = xdmf.read_mesh(name=“Grid”) line…
I’ve provided the .msh file that I now use to run the code:

12

I cannot reproduce your error with the following code:

import dolfinx
import meshio
from pathlib import Path
from mpi4py import MPI
import numpy
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)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh


path_msh = Path("patient9.msh")

msh = meshio.read(path_msh)

mesh_interior = create_mesh(msh, "tetra")
mesh_facets = create_mesh(msh, "triangle")
# store intermediate meshes in XDMF format

path_tmp_tets = path_msh.with_name(path_msh.stem + "_tets.xdmf")
path_tmp_tri = path_msh.with_name(path_msh.stem + "_tri.xdmf")

meshio.write(path_tmp_tets, mesh_interior, file_format="xdmf")
meshio.write(path_tmp_tri, mesh_facets, file_format="xdmf")

# read back in with dolfinx to obtain the mesh in the correct format
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, path_tmp_tets, "r") as xdmf:
    print("test1")
    mesh = xdmf.read_mesh(name="Grid")
    print("test2")

using docker

docker run -ti  -v $(pwd):/root/shared -w /root/shared --rm ghcr.io/fenics/dolfinx/dolfinx:nightly
 python3 -m pip install --no-binary=h5py h5py meshio