Add MED mesh format backend

.github/workflows/test_workflow.yml

@@ -34,6 +34,13 @@ env:
   # the image currently ships 5.0.10. MPICH ignores the variable, so it is
   # safe to set unconditionally. Drop once the image has OpenMPI >= 5.0.11.
   OMPI_MCA_topo: basic
+  # MEDCoupling (LGPL) writes the reference MED files in tests/test_med.py
+  MEDCOUPLING_VERSION: V9_15_0
+  MEDCOUPLING_PREFIX: /opt/medcoupling
+  # MED-fichier (LGPL) from the copy conda-forge uses, as Salome's server blocks
+  # automated downloads
+  MEDFILE_VERSION: "4.2"
+  MEDFILE_SHA256: a0193de4f0fef916ea9b83fb3e69e0717068a6bb41474ea8b6b0dece10b56d13
 jobs:
   get_image_tag:
     runs-on: ubuntu-latest
@@ -44,7 +51,8 @@ jobs:
         run: echo "image=${{ inputs.tag || env.DEFAULT_TAG }}" >> $GITHUB_OUTPUT
 
   test-code:
-    timeout-minutes: 30
+    # Leaves room for building MEDCoupling on a cache miss
+    timeout-minutes: 60
     runs-on: ubuntu-24.04
     # The container tag is now dynamic based on the input passed by the caller
     container: ghcr.io/fenics/test-env:current
@@ -109,6 +117,94 @@ jobs:
           python3 -m pip install --group dev --group test
           python3 -m pip install . --no-build-isolation
 
+      - name: Get Python version
+        id: python
+        run: echo "version=$(python3 -c 'import sys; print(f"{sys.version_info.major}.{sys.version_info.minor}")')" >> $GITHUB_OUTPUT
+
+      # MEDCoupling has no wheels for recent Python versions, so it is built from
+      # source, with MED-fichier, against the image's MPI HDF5. Ubuntu's
+      # MED-fichier 4.1.0 lacks API that MEDCoupling 9.15 uses.
+      - name: Install MEDCoupling dependencies
+        run: |
+          apt-get update
+          apt-get install -y --no-install-recommends swig
+
+      - name: Cache MEDCoupling
+        id: cache-medcoupling
+        uses: actions/cache@v6
+        with:
+          path: ${{ env.MEDCOUPLING_PREFIX }}
+          key: medcoupling-${{ env.MEDCOUPLING_VERSION }}-med${{ env.MEDFILE_VERSION }}-py${{ steps.python.outputs.version }}-${{ runner.arch }}
+
+      - name: Build MED-fichier
+        if: steps.cache-medcoupling.outputs.cache-hit != 'true'
+        run: |
+          wget -q -O /tmp/med.tar.gz https://github.com/Krande/med-copy/archive/refs/tags/${MEDFILE_VERSION}.tar.gz
+          echo "${MEDFILE_SHA256}  /tmp/med.tar.gz" | sha256sum -c -
+          mkdir /tmp/med && tar -xzf /tmp/med.tar.gz -C /tmp/med --strip-components=1
+          # MED-fichier 4.2 only accepts HDF5 1.10.x. Accept newer versions as Debian
+          # does (debian/patches/hdf5-1.14.patch); files are still written in the
+          # HDF5 1.8 format.
+          cd /tmp/med
+          sed -i 's/^\( *\)IF (NOT HDF_VERSION_MAJOR_REF EQUAL 1 OR NOT HDF_VERSION_MINOR_REF EQUAL 10 OR NOT HDF_VERSION_RELEASE_REF GREATER 1)/\1IF (HDF_VERSION_REF VERSION_LESS 1.10.2)/' config/cmake_files/medMacros.cmake
+          sed -i '/#if H5_VERS_MINOR > 10/,/#endif/s/#error/#warning/' src/ci/MEDfileCompatibility.c src/hdfi/*.c
+          sed -i 's/(_hmineur == HDF_VERSION_MINOR_REF  )/(_hmineur >= HDF_VERSION_MINOR_REF  )/' src/ci/MEDfileCompatibility.c
+          grep -q "VERSION_LESS 1.10.2" config/cmake_files/medMacros.cmake
+          cd -
+          # The 1.10 API keeps H5Oget_info and friends at their 1.10 signatures, as
+          # in Debian's build. Fortran is off: MEDCoupling only uses the C API, and
+          # without it med_int defaults to int.
+          cmake -G Ninja -S /tmp/med -B /tmp/med-build \
+            -DCMAKE_BUILD_TYPE=Release \
+            -DCMAKE_C_FLAGS="-DH5_USE_110_API" \
+            -DCMAKE_Fortran_COMPILER=NOTFOUND \
+            -DCMAKE_INSTALL_PREFIX=${MEDCOUPLING_PREFIX} \
+            -DHDF5_ROOT_DIR=${HDF5_ROOT} \
+            -DMEDFILE_BUILD_TESTS=OFF \
+            -DMEDFILE_INSTALL_DOC=OFF \
+            -DMEDFILE_BUILD_PYTHON=OFF
+          cmake --build /tmp/med-build
+          cmake --install /tmp/med-build
+
+      - name: Build MEDCoupling
+        if: steps.cache-medcoupling.outputs.cache-hit != 'true'
+        run: |
+          for repo in configuration medcoupling; do
+            git clone --depth 1 -b ${MEDCOUPLING_VERSION} https://github.com/SalomePlatform/${repo}.git /tmp/${repo}
+          done
+          # MEDCoupling requires CMake 2.8.11, which CMake 4 rejects; policies at 3.5
+          # also keep the legacy FindPythonInterp/FindPythonLibs that Salome uses.
+          # SALOME_USE_MPI is required for a parallel HDF5; ParaMEDMEM stays off.
+          # 32-bit ids match MED-fichier's default 32-bit med_int.
+          # Makefiles, not Ninja: the SWIG-generated .py files are compiled to .pyc
+          # by a rule that Ninja rejects, as they are not declared build outputs.
+          cmake -G "Unix Makefiles" -S /tmp/medcoupling -B /tmp/medcoupling-build \
+            -DCMAKE_POLICY_VERSION_MINIMUM=3.5 \
+            -DCMAKE_BUILD_TYPE=Release \
+            -DCMAKE_INSTALL_PREFIX=${MEDCOUPLING_PREFIX} \
+            -DCONFIGURATION_ROOT_DIR=/tmp/configuration \
+            -DMEDFILE_ROOT_DIR=${MEDCOUPLING_PREFIX} \
+            -DSALOME_USE_MPI=ON \
+            -DMEDCOUPLING_USE_MPI=OFF \
+            -DMEDCOUPLING_USE_64BIT_IDS=OFF \
+            -DMEDCOUPLING_ENABLE_PARTITIONER=OFF \
+            -DMEDCOUPLING_ENABLE_RENUMBER=OFF \
+            -DMEDCOUPLING_MEDLOADER_USE_XDR=OFF \
+            -DMEDCOUPLING_BUILD_TESTS=OFF \
+            -DMEDCOUPLING_BUILD_PY_TESTS=OFF \
+            -DMEDCOUPLING_BUILD_DOC=OFF \
+            -DMEDCOUPLING_WITH_FILE_EXAMPLES=OFF
+          cmake --build /tmp/medcoupling-build --parallel $(nproc)
+          cmake --install /tmp/medcoupling-build
+
+      - name: Make MEDCoupling importable
+        run: |
+          echo "PYTHONPATH=${MEDCOUPLING_PREFIX}/lib/python${{ steps.python.outputs.version }}/site-packages${PYTHONPATH:+:${PYTHONPATH}}" >> $GITHUB_ENV
+          echo "LD_LIBRARY_PATH=${MEDCOUPLING_PREFIX}/lib${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}" >> $GITHUB_ENV
+
+      - name: Show MEDCoupling version
+        run: python3 -c "import medcoupling; print(medcoupling.__version__)"
+
       - name: Show adios2 version
         run: python3 -c "import adios2; print(adios2.__version__)"
 

_toc.yml

@@ -38,6 +38,7 @@ parts:
     - file: "docs/backends/pyvista"
     - file: "docs/backends/xdmf"
     - file: "docs/backends/vtkhdf"
+    - file: "docs/backends/med"
     - file: "docs/backends/gmsh"
 
   - caption: Contributing

docs/backends/med.rst

@@ -0,0 +1,46 @@
+.. _medbackend:
+
+MED-backend
+-----------
+
+Reads meshes and mesh tags from `MED <https://docs.salome-platform.org/latest/dev/MEDCoupling/developer/med-file.html>`_
+files, as exported by Salome and Code_Aster. The backend is read-only and relies on :code:`h5py`,
+which must be MPI compatible to read in parallel (see the :ref:`H5PY backend<h5pybackend>`).
+
+In a MED file each entity stores a *family* id, and each family lists the named *groups* it belongs to.
+An entity in several groups gets a family of its own for that combination of groups.
+Mesh tags can be read in two ways:
+
+- ``meshtag_name="cell"`` or ``meshtag_name="facet"`` returns the family ids of the cells or facets
+  (entities without a family are skipped). MED element family ids are negative, while integral
+  subdomain ids must be non-negative, so the ids are negated: family ``-6`` gives tag ``6``.
+- Any other ``meshtag_name`` is a group name, and returns the entities of that group tagged with ``1``.
+
+.. code-block:: python
+
+    mesh = io4dolfinx.read_mesh("mesh.med", MPI.COMM_WORLD, backend="med")
+    facet_families = io4dolfinx.read_meshtags("mesh.med", mesh, meshtag_name="facet", backend="med")
+    inlet = io4dolfinx.read_meshtags("mesh.med", mesh, meshtag_name="inlet", backend="med")
+
+Use :func:`io4dolfinx.backends.med.backend.read_groups` to list the groups and their family ids.
+If a file holds several meshes, select one with ``mesh_name``.
+
+A time-dependent mesh stores one step per time. Select the step with ``time`` when reading the mesh,
+and with ``backend_args={"time": t}`` when reading its tags (the first step is used if no time is
+given). A mesh with a single step is read whatever the time.
+
+.. code-block:: python
+
+    mesh = io4dolfinx.read_mesh("mesh.med", MPI.COMM_WORLD, backend="med", time=0.5)
+    cells = io4dolfinx.read_meshtags(
+        "mesh.med", mesh, meshtag_name="cell", backend="med", backend_args={"time": 0.5}
+    )
+
+The cells must be of a single type. Linear and quadratic Lagrange cells (``SE3``, ``TR6``, ``QU9``,
+``T10``, ``P18``, ``H27``) and the quadratic serendipity cells ``QU8`` and ``H20`` are supported.
+The serendipity cells require the main branch of DOLFINx.
+MED node orderings are converted with tables derived from MEDCoupling's MED cell model.
+
+.. automodule:: io4dolfinx.backends.med.backend
+    :members:
+    :exclude-members: read_attributes, read_timestamps, write_attributes, write_mesh, write_meshtags, read_dofmap, read_dofs, read_cell_perms, write_function, read_legacy_mesh, snapshot_checkpoint, read_hdf5_array, read_point_data, read_cell_data, read_function_names, write_data

pyproject.toml

@@ -27,6 +27,7 @@ h5py = ["h5py"]
 pyvista = ["pyvista"]
 xdmf = ["h5py"]
 exodus = ["netcdf4"]
+med = ["h5py"]
 gmsh = ["gmsh"]
 
 [dependency-groups]
@@ -46,6 +47,9 @@ test = [
     "io4dolfinx[h5py]",
     "io4dolfinx[pyvista]",
     "io4dolfinx[xdmf]",
+    # Reference MED writer for tests. Wheels only exist up to CPython 3.13 on x86_64;
+    # CI builds it from source otherwise. Bump the bound when newer wheels appear.
+    "medcoupling>=9.15; python_version < '3.14' and (platform_machine == 'x86_64' or platform_machine == 'AMD64')",
     "io4dolfinx[gmsh]",
 ]
 dev = [

src/io4dolfinx/backends/__init__.py

@@ -50,7 +50,7 @@ def get_mesh_name(backend_args: dict[str, Any] | None) -> str:
 def set_default_backend(backend: str):
     """Set the global default backend for io4dolfinx."""
     global _DEFAULT_BACKEND
-    if backend not in BUILTIN_BAKENDS:
+    if backend not in BUILTIN_BACKENDS:
         try:
             get_backend(backend)
         except ImportError as e:
@@ -504,6 +504,10 @@ def get_backend(backend: str | None = None) -> IOBackend:
         from .exodus import backend as EXODUSInterface
 
         return EXODUSInterface
+    elif backend == "med":
+        from .med import backend as MEDInterface
+
+        return MEDInterface
     elif backend == "gmsh":
         from .gmsh import backend as GMSHInterface
 
@@ -512,7 +516,16 @@ def get_backend(backend: str | None = None) -> IOBackend:
         return import_module(backend)
 
 
-BUILTIN_BAKENDS = ("h5py", "adios2", "pyvista", "xdmf", "vtkhdf", "exodus", "gmsh")
+BUILTIN_BACKENDS = (
+    "h5py",
+    "adios2",
+    "pyvista",
+    "xdmf",
+    "vtkhdf",
+    "exodus",
+    "med",
+    "gmsh",
+)
 
 
 def list_builtin_backends() -> list[str]:
@@ -522,7 +535,7 @@ def list_builtin_backends() -> list[str]:
         List of available backends
     """
     lst = []
-    for backend in BUILTIN_BAKENDS:
+    for backend in BUILTIN_BACKENDS:
         try:
             get_backend(backend)
         except ImportError:

src/io4dolfinx/backends/h5py/backend.py

@@ -60,6 +60,41 @@ def h5pyfile(h5name, filemode="r", force_serial: bool = False, comm=None):
             h5file.close()
 
 
+def read_component_major_rows(
+    dataset: h5py.Dataset, num_entities_global: int, local_range: list[int]
+) -> np.ndarray:
+    """Read rows ``local_range`` of a 2D array stored component by component.
+
+    The array is a one-dimensional dataset holding all first components, then all
+    second components and so on, i.e. ``x0 x1 ... xn y0 y1 ... yn z0 z1 ... zn``
+    rather than ``x0 y0 z0 x1 y1 z1 ...`` (``MED_NO_INTERLACE`` in the MED format).
+    The rows are one block per component, read as a single hyperslab.
+
+    Args:
+        dataset: One-dimensional dataset
+        num_entities_global: Number of entities (rows) in the whole dataset
+        local_range: Range ``[start, end)`` of rows to read on this process
+
+    Returns:
+        Array of shape ``(end - start, num_components)``, where
+        ``num_components = dataset.size // num_entities_global``
+    """
+    if num_entities_global <= 0 or dataset.size == 0 or dataset.size % num_entities_global != 0:
+        raise ValueError(
+            f"Size {dataset.size} of {dataset.name} is not a positive multiple of "
+            f"{num_entities_global} entities."
+        )
+    num_components = dataset.size // num_entities_global
+    num_rows = local_range[1] - local_range[0]
+    if num_rows == 0:
+        # h5py.MultiBlockSlice rejects empty blocks
+        return np.empty((0, num_components), dtype=dataset.dtype)
+    blocks = h5py.MultiBlockSlice(
+        start=local_range[0], stride=num_entities_global, count=num_components, block=num_rows
+    )
+    return np.ascontiguousarray(dataset[blocks].reshape(num_components, num_rows).T)
+
+
 def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]:
     args = dict(arguments) if arguments else {}
     args.setdefault("legacy", False)  # If meshtags is read from legacy

src/io4dolfinx/backends/med/__init__.py

@@ -0,0 +1,3 @@
+from . import backend
+
+__all__ = ["backend"]