MED Salome to dolfinx: Issues with XDMF

stefano.tronci · March 11, 2023, 4:22pm

Hi guys!

Some time ago I managed to write some code that can read MED files into FEniCS. Now I am refactoring that code to work with FEniCSx, but I am running into issues. This is what I have so far:

import typing
import pathlib
import meshio
import enum
import numpy as np
import tempfile
import mpi4py
import dolfinx


def load_raw(path: typing.Union[pathlib.Path, str]) -> meshio.Mesh:
    return meshio.read(path)


CellTagsT = typing.Dict[int, typing.List[str]]


def calculate_tags_offset(cell_tags: CellTagsT) -> int:
    tags_offset = 0
    min_tag = min(k for k in cell_tags.keys())
    if min_tag < 0:
        tags_offset = -min_tag

    return tags_offset


def get_cell_tags(mesh: meshio.Mesh) -> typing.Tuple[CellTagsT]:
    cell_tags = mesh.cell_tags
    tags_offset = calculate_tags_offset(cell_tags=cell_tags)
    return {
        k + tags_offset: cell_tags[k]
        for k in cell_tags
    }


class CellType(enum.Enum):
    HEXAHEDRON = 'hexahedron'
    INTERVAL = 'line'
    PRISM = 'prism'
    PYRAMID = 'pyramid'
    QUADRILATERAL = 'quad'
    TETRAHEDRON = 'tetra'
    TRIANGLE = 'triangle'


def convert_entity(
        comm: mpi4py.MPI,
        mesh: meshio.Mesh, 
        cell_type: CellType,
        tag_key: str,
        prune_mask: np.ndarray,
        coordinate_scaling: float
) -> typing.Tuple[dolfinx.mesh.Mesh, dolfinx.mesh.MeshTagsMetaClass]:
    cells = mesh.get_cells_type(cell_type=cell_type.value)
    cell_data = mesh.get_cell_data(name=tag_key, cell_type=cell_type.value)
    points = mesh.points[:, prune_mask] * coordinate_scaling
    tags_offset = calculate_tags_offset(cell_tags=mesh.cell_tags)

    entity_mesh = meshio.Mesh(
        points=points,
        cells={cell_type.value: cells},
        cell_data={'marker': [cell_data + tags_offset]}
    )

    with tempfile.TemporaryDirectory() as xdmf_dir:
        target_xdmf_dir = pathlib.Path(xdmf_dir)
        tmp_path = target_xdmf_dir.joinpath(cell_type.value + '_mesh.xdmf')
        # entity_mesh.write(tmp_path)
        meshio.write(tmp_path, entity_mesh)
        
        with dolfinx.io.XDMFFile(comm=comm, filename=tmp_path, file_mode='r') as f:
            out_mesh = f.read_mesh(name='Grid')
            out_tags = f.read_meshtags(out_mesh, name='Grid')
            out_mesh.topology.create_connectivity(
                *[mesh.topology.dim - n for n in range(mesh.topology.dim)]
            )

    return out_mesh, out_tags


convert_entity(
    comm=mpi4py.MPI.COMM_WORLD,
    mesh=load_raw('/home/admin/test.med'),
    cell_type=CellType.TETRAHEDRON,
    tag_key='cell_tags',
    prune_mask=np.zeros(3, dtype=bool),
    coordinate_scaling=1.0
)

A bit long, but I guess the only important bit is the convert_entity function. Basically, it is all based off:

However, when I use dolfinx.io.XDMFFile’s read_mesh method I get this error:

Traceback (most recent call last):
  File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher/../../debugpy/__main__.py", line 39, in <module>
    cli.main()
  File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher/../../debugpy/../debugpy/server/cli.py", line 430, in main
    run()
  File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher/../../debugpy/../debugpy/server/cli.py", line 284, in run_file
    runpy.run_path(target, run_name="__main__")
  File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/_vendored/pydevd/_pydevd_bundle/pydevd_runpy.py", line 321, in run_path
    return _run_module_code(code, init_globals, run_name,
  File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/_vendored/pydevd/_pydevd_bundle/pydevd_runpy.py", line 135, in _run_module_code
    _run_code(code, mod_globals, init_globals,
  File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/_vendored/pydevd/_pydevd_bundle/pydevd_runpy.py", line 124, in _run_code
    exec(code, run_globals)
  File "/home/admin/projects/acoupy/acoupy_meshutil/acoupy_meshutil/util.py", line 81, in <module>
    convert_entity(
  File "/home/admin/projects/acoupy/acoupy_meshutil/acoupy_meshutil/util.py", line 72, in convert_entity
    out_mesh = f.read_mesh(name='Grid')
  File "/usr/lib/python3.10/site-packages/dolfinx/io/utils.py", line 163, in read_mesh
    x = super().read_geometry_data(name, xpath)
RuntimeError: Cannot determine geometric dimension. GeometryType "" in XDMF file is unknown or unsupported

If I hack the XDMF file and I change it to <Geometry GeometryType="XYZ"> things get even worse?

python: /build/dolfinx/src/dolfinx-0.6.0/cpp/dolfinx/io/xdmf_mesh.cpp:264: std::pair<std::vector<double>, std::array<long unsigned int, 2> > dolfinx::io::xdmf_mesh::read_geometry_data(MPI_Comm, hid_t, const pugi::xml_node&): Assertion `gdims[1] == (int)gdim' failed.

Loguru caught a signal: SIGABRT
Stack trace:
53      0x5606f6449045 _start + 37
52      0x7fcf143ee84a __libc_start_main + 138
51      0x7fcf143ee790 /usr/lib/libc.so.6(+0x23790) [0x7fcf143ee790]
50      0x7fcf141e336b Py_BytesMain + 59
49      0x7fcf14096057 /usr/lib/libpython3.10.so.1.0(+0x96057) [0x7fcf14096057]
48      0x7fcf14216617 /usr/lib/libpython3.10.so.1.0(+0x216617) [0x7fcf14216617]
47      0x7fcf141562d9 _PyFunction_Vectorcall + 121
46      0x7fcf141454d6 _PyEval_EvalFrameDefault + 838
45      0x7fcf141562d9 _PyFunction_Vectorcall + 121
44      0x7fcf141454d6 _PyEval_EvalFrameDefault + 838
43      0x7fcf141564cf /usr/lib/libpython3.10.so.1.0(+0x1564cf) [0x7fcf141564cf]
42      0x7fcf141f7f1b /usr/lib/libpython3.10.so.1.0(+0x1f7f1b) [0x7fcf141f7f1b]
41      0x7fcf141f1b94 PyEval_EvalCode + 148
40      0x7fcf14144120 /usr/lib/libpython3.10.so.1.0(+0x144120) [0x7fcf14144120]
39      0x7fcf1414a03f _PyEval_EvalFrameDefault + 20143
38      0x7fcf141562d9 _PyFunction_Vectorcall + 121
37      0x7fcf141454d6 _PyEval_EvalFrameDefault + 838
36      0x7fcf141562d9 _PyFunction_Vectorcall + 121
35      0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
34      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
33      0x7fcf141562d9 _PyFunction_Vectorcall + 121
32      0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
31      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
30      0x7fcf141562d9 _PyFunction_Vectorcall + 121
29      0x7fcf1414c696 _PyEval_EvalFrameDefault + 29958
28      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
27      0x7fcf141562d9 _PyFunction_Vectorcall + 121
26      0x7fcf1414c696 _PyEval_EvalFrameDefault + 29958
25      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
24      0x7fcf141564cf /usr/lib/libpython3.10.so.1.0(+0x1564cf) [0x7fcf141564cf]
23      0x7fcf141f7f1b /usr/lib/libpython3.10.so.1.0(+0x1f7f1b) [0x7fcf141f7f1b]
22      0x7fcf141f1b94 PyEval_EvalCode + 148
21      0x7fcf14144120 /usr/lib/libpython3.10.so.1.0(+0x144120) [0x7fcf14144120]
20      0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
19      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
18      0x7fcf141562d9 _PyFunction_Vectorcall + 121
17      0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
16      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
15      0x7fcf14161446 /usr/lib/libpython3.10.so.1.0(+0x161446) [0x7fcf14161446]
14      0x7fcf1414b3e5 _PyEval_EvalFrameDefault + 25173
13      0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
12      0x7fcf141616a6 /usr/lib/libpython3.10.so.1.0(+0x1616a6) [0x7fcf141616a6]
11      0x7fcf1414f4eb _PyObject_MakeTpCall + 683
10      0x7fcf14155e21 /usr/lib/libpython3.10.so.1.0(+0x155e21) [0x7fcf14155e21]
9       0x7fcf01674129 /usr/lib/python3.10/site-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x74129) [0x7fcf01674129]
8       0x7fcf017984fb /usr/lib/python3.10/site-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x1984fb) [0x7fcf017984fb]
7       0x7fcf01537210 dolfinx::io::XDMFFile::read_geometry_data(std::string, std::string) const + 624
6       0x7fcf0154d62e dolfinx::io::xdmf_mesh::read_geometry_data(ompi_communicator_t*, long, pugi::xml_node const&) + 5598
5       0x7fcf143fc9f6 /usr/lib/libc.so.6(+0x319f6) [0x7fcf143fc9f6]
4       0x7fcf143ed45c /usr/lib/libc.so.6(+0x2245c) [0x7fcf143ed45c]
3       0x7fcf143ed53d abort + 215
2       0x7fcf14403ea8 gsignal + 24
1       0x7fcf144528ec /usr/lib/libc.so.6(+0x878ec) [0x7fcf144528ec]
0       0x7fcf14403f50 /usr/lib/libc.so.6(+0x38f50) [0x7fcf14403f50]
2023-03-11 16:07:34.629 (  44.973s) [main            ]                       :0     FATL| Signal: SIGABRT

It makes me think the mesh is the problem but it is similar to others I successfully loaded into FEnICS. The mesh is here.

Any idea?

dokken · March 11, 2023, 4:49pm

Your prune_mask is wrong, as you fill it with zeros, you actually remove all coordinates.
So if you instead do:

    points = mesh.points[:, ~prune_mask] * coordinate_scaling
    tags_offset = calculate_tags_offset(cell_tags=mesh.cell_tags)
    entity_mesh = meshio.Mesh(
        points=points,
        cells={cell_type.value: cells},
        cell_data={'marker': [cell_data + tags_offset]}
    )

    with tempfile.TemporaryDirectory() as xdmf_dir:
        target_xdmf_dir = pathlib.Path(xdmf_dir)
        tmp_path = target_xdmf_dir.joinpath(cell_type.value + '_mesh.xdmf')
        # entity_mesh.write(tmp_path)
        meshio.write(tmp_path, entity_mesh)

        with dolfinx.io.XDMFFile(comm=comm, filename=tmp_path, file_mode='r') as f:
            out_mesh = f.read_mesh(name='Grid')
            out_tags = f.read_meshtags(out_mesh, name='Grid')
    return out_mesh, out_tags

the code runs

stefano.tronci · March 11, 2023, 5:23pm

Ouch. Classic me coding a dumb bug. Thanks for your help! I see you removed the call to out_mesh.topology.create_connectivity. I put a couples of bugs there to. I meant to do

            out_mesh.topology.create_connectivity(
                *[out_mesh.topology.dim - n for n in range(out_mesh.topology.dim)]
            )

But that will produce

Traceback (most recent call last):
  File "/home/admin/projects/acoupy/acoupy_meshutil/acoupy_meshutil/util.py", line 80, in <module>
    convert_entity(
  File "/home/admin/projects/acoupy/acoupy_meshutil/acoupy_meshutil/util.py", line 73, in convert_entity
    out_mesh.topology.create_connectivity(
TypeError: create_connectivity(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.mesh.Topology, d0: int, d1: int) -> None

Invoked with: <dolfinx.cpp.mesh.Topology object at 0x7f994428c810>, 3, 2, 1

I get that create_connectivity is meant for 2D meshes?

dokken · March 11, 2023, 7:26pm

No, create_connectivity is used to make maps between different topological dimensions (cells, facets,edges,vertices).

You create connectivity between two of these at the time. For instance if you want to create the cell to facet connectivity, you would call mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)

If you want to create the vertex to cell connectivity you would call
mesh.topology.create_connectivity(0, mesh.topology.dim)

You would usually not call these functions except if you Need the specifically, or a dolfinx function throws and error saying you need to create any of them.

Ekrem_Ekici · March 24, 2023, 3:13pm

I used your code like this;

import typing
import pathlib
import meshio
import enum
import numpy as np
import tempfile
import mpi4py
import dolfinx


def load_raw(path: typing.Union[pathlib.Path, str]) -> meshio.Mesh:
    return meshio.read(path)


CellTagsT = typing.Dict[int, typing.List[str]]


def calculate_tags_offset(cell_tags: CellTagsT) -> int:
    tags_offset = 0
    min_tag = min(k for k in cell_tags.keys())
    if min_tag < 0:
        tags_offset = -min_tag

    return tags_offset


def get_cell_tags(mesh: meshio.Mesh) -> typing.Tuple[CellTagsT]:
    cell_tags = mesh.cell_tags
    tags_offset = calculate_tags_offset(cell_tags=cell_tags)
    return {
        k + tags_offset: cell_tags[k]
        for k in cell_tags
    }


class CellType(enum.Enum):
    HEXAHEDRON = 'hexahedron'
    INTERVAL = 'line'
    PRISM = 'prism'
    PYRAMID = 'pyramid'
    QUADRILATERAL = 'quad'
    TETRAHEDRON = 'tetra'
    TRIANGLE = 'triangle'


def convert_entity(
        comm: mpi4py.MPI,
        mesh: meshio.Mesh, 
        cell_type: CellType,
        tag_key: str,
        prune_mask: np.ndarray,
        coordinate_scaling: float
) -> typing.Tuple[dolfinx.mesh.Mesh, dolfinx.mesh.MeshTagsMetaClass]:
    cells = mesh.get_cells_type(cell_type=cell_type.value)
    cell_data = mesh.get_cell_data(name=tag_key, cell_type=cell_type.value)
    points = mesh.points[:, ~prune_mask] * coordinate_scaling
    tags_offset = calculate_tags_offset(cell_tags=mesh.cell_tags)
    entity_mesh = meshio.Mesh(
        points=points,
        cells={cell_type.value: cells},
        cell_data={'marker': [cell_data + tags_offset]}
    )

    with tempfile.TemporaryDirectory() as xdmf_dir:
        target_xdmf_dir = pathlib.Path(xdmf_dir)
        tmp_path = target_xdmf_dir.joinpath(cell_type.value + '_mesh.xdmf')
        # entity_mesh.write(tmp_path)
        meshio.write(tmp_path, entity_mesh)

        with dolfinx.io.XDMFFile(comm=comm, filename=tmp_path, file_mode='r') as f:
            out_mesh = f.read_mesh(name='Grid')
            out_tags = f.read_meshtags(out_mesh, name='Grid')
    return out_mesh, out_tags




volume_mesh, tag_mesh = convert_entity(
    comm=mpi4py.MPI.COMM_WORLD,
    mesh=load_raw('pipe.med'),
    cell_type=CellType.TETRAHEDRON,
    tag_key='cell_tags',
    prune_mask=np.zeros(3, dtype=bool),
    coordinate_scaling=1.0
)

name = "pipe"
xdmf_name = name + ".xdmf"
xdmf_tags_name = name + "_tags.xdmf"
meshio.write(xdmf_name, volume_mesh)
meshio.write(xdmf_tags_name, tag_mesh)

to read this med file. But it gives;

  File "mesh.py", line 80, in <module>
    mesh=load_raw('pipe.med'),
  File "mesh.py", line 12, in load_raw
    return meshio.read(path)
  File "/home/ee331/Dev/Venvs/v060/lib/python3.8/site-packages/meshio/_helpers.py", line 71, in read
    return _read_file(Path(filename), file_format)
  File "/home/ee331/Dev/Venvs/v060/lib/python3.8/site-packages/meshio/_helpers.py", line 103, in _read_file
    return reader_map[file_format](str(path))
  File "/home/ee331/Dev/Venvs/v060/lib/python3.8/site-packages/meshio/med/_med.py", line 114, in read
    mesh = Mesh(
  File "/home/ee331/Dev/Venvs/v060/lib/python3.8/site-packages/meshio/_mesh.py", line 181, in __init__
    if len(data[k]) != len(self.cells[k]):
TypeError: len() of unsized object

I didn’t understand why it is not detecting the cells?

dokken · March 24, 2023, 5:27pm

I would suggest you make an issue at the Meshio issue tracker as this seems meshio specific.

stefano.tronci · March 26, 2023, 3:11pm

Yes, I agree. All is needed to reproduce the issue is

>>> import meshio
>>> meshio.read('pipe.med')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.10/site-packages/meshio/_helpers.py", line 71, in read
    return _read_file(Path(filename), file_format)
  File "/usr/lib/python3.10/site-packages/meshio/_helpers.py", line 103, in _read_file
    return reader_map[file_format](str(path))
  File "/usr/lib/python3.10/site-packages/meshio/med/_med.py", line 114, in read
    mesh = Mesh(
  File "/usr/lib/python3.10/site-packages/meshio/_mesh.py", line 181, in __init__
    if len(data[k]) != len(self.cells[k]):
TypeError: len() of unsized object

I can open your mesh with gmsh so it is probably something to bring to meshio attention if you haven’t already.

Ekrem_Ekici · March 26, 2023, 4:06pm

I can work with gmsh without any problem. I just wanted to extend my code’s capability of importing different mesh formats. I will post this on issue tracker of meshio then. But it seems like meshio is not properly maintained.

dokken · March 26, 2023, 6:34pm

If you want a general way to read in meshes to DOLFINx, you could work with the raw data, as described at
https://jsdokken.com/dolfinx_docs/meshes.html