MED format converter

LICENSE

Copyright 2026 Jørgen S. Dokken and Simula Research Laboratory

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

script.py

# READ MED files into DOLFINx with cell and facet tags
# Author: Jørgen S. Dokken and Simula Research Laboratory
# SPDX-License-Identifier: MIT

import h5py
from pathlib import Path
import dolfinx
import numpy as np
import basix.ufl
from mpi4py import MPI
from dolfinx.cpp.mesh import create_mesh, create_cell_partitioner


# MAP: type to (celltype, degree)
_salome_to_dolfinx = {
    "SE2": (dolfinx.mesh.CellType.interval, 1),
    "TR3": (dolfinx.mesh.CellType.triangle, 1),
    "QU4": (dolfinx.mesh.CellType.quadrilateral, 1),
    "TE4": (dolfinx.mesh.CellType.tetrahedron, 1),
    "PY5": (dolfinx.mesh.CellType.pyramid, 1),
    "HE8": (dolfinx.mesh.CellType.hexahedron, 1),
}

def compute_local_range(comm: MPI.Intracomm, N: np.int64):
    """
    Divide a set of `N` objects into `M` partitions, where `M` is
    the size of the MPI communicator `comm`.

    NOTE: If N is not divisible by the number of ranks, the first `r`
    processes gets an extra value

    Returns the local range of values
    """
    rank = comm.rank
    size = comm.size
    n = N // size
    r = N % size
    # First r processes has one extra value
    if rank < r:
        return [rank * (n + 1), (rank + 1) * (n + 1)]
    else:
        return [rank * n + r, (rank + 1) * n + r]


def read_med_mesh(filename: str | Path, comm: MPI.Intracomm):
    facet_data = None
    if comm.size == 1:
        h5_kwargs={}
    else:
        h5_kwargs= {"driver" : "mpio",
        "comm": MPI.COMM_WORLD}
    with h5py.File(filename, "r", **h5_kwargs) as h5file:
        if "ENS_MAA" not in h5file.keys():
            raise ValueError(
                f"Could not find a set of meshes (Ensemble de Maillages) in {filename}."
            )
        mesh_set = h5file["ENS_MAA"]
        meshes = mesh_set.keys()
        if len(meshes) > 1:
            raise RuntimeError(f"Multiple meshes found in {filename}")
        mesh = mesh_set[meshes.__iter__().__next__()]
        time_steps = mesh.keys()
        if len(time_steps) == 1:
            mesh = mesh[time_steps.__iter__().__next__()]

            # Nodes must be re-ordered according to NUM
            nodes = mesh["NOE"]
            attrs = mesh.parent.attrs
            if "DIM" in attrs:
                gdim = attrs["DIM"]
            elif "ESP" in attrs:
                gdim = attrs["ESP"]
            elif "ESPACE" in attrs:
                gdim = attrs["ESPACE"]
            else:
                gdim = 3  # Fallback assumption

            # Read all nodes on rank 0 as they must be reordered
            num_nodes = nodes["NUM"].size
            coords = nodes["COO"]
            in_dim = int(coords.size / num_nodes)
            assert coords.size % num_nodes == 0
            if comm.rank == 0:
                coordinates = coords[:].reshape(in_dim, -1).T
                node_ordering = nodes["NUM"][:]
                assert min(node_ordering) == 1
            else:
                coordinates = np.zeros((0, in_dim), dtype=np.float64)
                node_ordering = np.zeros(0, dtype=np.int64)

            # node_marker = nodes["FAM"]
            if len(node_ordering) != coordinates.shape[0]:
                raise RuntimeError("Not all coordinates have global index")
            coordinates = coordinates[node_ordering[:] - 1].astype(np.float64).copy()

            # Create a sorter to map Global ID -> Local Index (0, 1, 2...)
            # This allows us to look up where ID '105' is located in the array

            # Find cell types of highest topological dimension
            elements = mesh["MAI"]
            cell_types = []
            # Sort to have same accessing on each process to align
            # coordinate elements
            element_keys = sorted(list(elements.keys()))
            for el in element_keys:
                dx_type = _salome_to_dolfinx[el]
                cell_types.append((el, dx_type))

            def tdim(ct: tuple[str, tuple[dolfinx.mesh.CellType, int]]) -> int:
                return dolfinx.mesh.cell_dim(ct[1][0])

            max_tdim = tdim(max(cell_types, key=tdim))
            cells = filter(lambda el: tdim(el) == max_tdim, cell_types)
            facets = list(filter(lambda el: tdim(el) == max_tdim - 1, cell_types))

            topologies = []
            top_markers = []
            coord_els = []
            for salome_cell, (ct, degree) in cells:
                # Ignore NUM key, we do not use it
                top = elements[salome_cell]
                num_cells = top["NUM"].size
                local_cell_range = compute_local_range(comm, num_cells)

                # Extract data and permute (assume VTK ordering)
                num_dofs = int(top["NOD"].size / num_cells)
                assert np.isclose(top["NOD"].size % num_cells, 0)
                conn_data = np.zeros((local_cell_range[1]-local_cell_range[0], num_dofs), dtype=np.int64)
                for i in range(num_dofs):
                    conn_data[:,i] = top["NOD"][local_cell_range[0]+i*num_cells:local_cell_range[1]+i*num_cells]

                permutation = np.argsort(
                    dolfinx.cpp.io.perm_vtk(ct, conn_data.shape[1])
                )
                conn_flattened = (conn_data[:, permutation] - 1).flatten()

                # Store data
                top_markers.append(top["FAM"][slice(*local_cell_range)])
                topologies.append(conn_flattened)
                coord_els.append(
                    dolfinx.fem.coordinate_element(
                        basix.ufl.element(
                            "Lagrange",
                            dolfinx.mesh.to_string(ct),
                            degree,
                            shape=(int(gdim),),
                        ).basix_element
                    )._cpp_object
                )
            assert len(facets) < 2
            if len(facets) == 1:
                salome_facet, (ft, _) = facets[0]
                ftop = elements[salome_facet]
                num_facets = ftop["NUM"].size

                # Extract data and permute (assume VTK ordering)
                num_facet_dofs = int(ftop["NOD"].size / num_facets)
                assert np.isclose(ftop["NOD"].size % num_facets, 0)
                local_facet_range = compute_local_range(comm, num_facets)
                conn_data = np.zeros((local_facet_range[1]-local_facet_range[0], num_facet_dofs), dtype=np.int64)
                for i in range(num_facet_dofs):
                    conn_data[:,i] = ftop["NOD"][local_facet_range[0]+i*num_facets:local_facet_range[1]+i*num_facets]

                permutation = np.argsort(
                    dolfinx.cpp.io.perm_vtk(ft, conn_data.shape[1])
                )
                facet_conn = conn_data[:, permutation] - 1
                facet_values = ftop["FAM"][slice(*local_facet_range)]
                facet_data = (facet_conn, facet_values)

        else:
            raise RuntimeError(
                f"We do not support time dependent meshes from {filename}"
            )
    if comm.size == 1:
        assert coordinates.shape[0] == len(np.unique(np.hstack(topologies)))
        assert min(np.unique(np.hstack(topologies))) == 0
        assert (
            max(np.unique(np.hstack(topologies)))
            == len(np.unique(np.hstack(topologies))) - 1
        )

    part = create_cell_partitioner(dolfinx.mesh.GhostMode.none)
    _cpp_mesh = create_mesh(comm, topologies, coord_els, coordinates, part, 2)
    mesh = dolfinx.mesh.Mesh(_cpp_mesh, None)
    if len(topologies) == 1:
        cell_data, cell_values = dolfinx.io.distribute_entity_data(
            mesh,
            mesh.topology.dim,
            topologies[0].reshape(-1, mesh.geometry.cmap.dim),
            top_markers[0],
        )
        ct = dolfinx.mesh.meshtags_from_entities(
            mesh,
            mesh.topology.dim,
            dolfinx.graph.adjacencylist(cell_data),
            cell_values.astype(np.int32),
        )
        ct.name = "CellTags"
        if facet_data is not None:

            fdata, fvalues = dolfinx.io.distribute_entity_data(
                mesh, mesh.topology.dim - 1, *facet_data
            )
            mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
            ft = dolfinx.mesh.meshtags_from_entities(
                mesh,
                mesh.topology.dim - 1,
                dolfinx.graph.adjacencylist(fdata),
                fvalues.astype(np.int32),
            )
            ft.name = "FacetTags"
            return mesh, ct, ft
        return mesh, ct, None
    else:
        return mesh, None, None


if __name__ == "__main__":
    filename = Path("Mesh_2tetraRegions.med")
    mesh, ct, ft = read_med_mesh(filename, MPI.COMM_WORLD)

    dolfinx.io.vtkhdf.write_mesh("test_med.vtkhdf", mesh)

    with dolfinx.io.XDMFFile(mesh.comm, "test_med.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        if ct is not None:
            xdmf.write_meshtags(ct, mesh.geometry)
        if ft is not None:
            xdmf.write_meshtags(ft, mesh.geometry)

    with dolfinx.io.XDMFFile(mesh.comm, "test_med.xdmf", "r") as xdmf:
        mesh = xdmf.read_mesh(name="mesh")
        ct = xdmf.read_meshtags(mesh, name="CellTags")
        mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
        ft = xdmf.read_meshtags(mesh, name="FacetTags")