diff --git a/.github/workflows/test_workflow.yml b/.github/workflows/test_workflow.yml index 8da9490..2d67d50 100644 --- a/.github/workflows/test_workflow.yml +++ b/.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__)" diff --git a/_toc.yml b/_toc.yml index a84833c..229ba5c 100644 --- a/_toc.yml +++ b/_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 diff --git a/docs/backends/med.rst b/docs/backends/med.rst new file mode 100644 index 0000000..b3a5b90 --- /dev/null +++ b/docs/backends/med.rst @@ -0,0 +1,46 @@ +.. _medbackend: + +MED-backend +----------- + +Reads meshes and mesh tags from `MED `_ +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`). + +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 diff --git a/pyproject.toml b/pyproject.toml index 0b337e3..44aae36 100644 --- a/pyproject.toml +++ b/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 = [ diff --git a/src/io4dolfinx/backends/__init__.py b/src/io4dolfinx/backends/__init__.py index 54fa994..e82f1b2 100644 --- a/src/io4dolfinx/backends/__init__.py +++ b/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: diff --git a/src/io4dolfinx/backends/h5py/backend.py b/src/io4dolfinx/backends/h5py/backend.py index f9cc83b..be99ecb 100644 --- a/src/io4dolfinx/backends/h5py/backend.py +++ b/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 diff --git a/src/io4dolfinx/backends/med/__init__.py b/src/io4dolfinx/backends/med/__init__.py new file mode 100644 index 0000000..05458f8 --- /dev/null +++ b/src/io4dolfinx/backends/med/__init__.py @@ -0,0 +1,3 @@ +from . import backend + +__all__ = ["backend"] diff --git a/src/io4dolfinx/backends/med/backend.py b/src/io4dolfinx/backends/med/backend.py new file mode 100644 index 0000000..e74b425 --- /dev/null +++ b/src/io4dolfinx/backends/med/backend.py @@ -0,0 +1,777 @@ +""" +MED interface to io4dolfinx. + +MED is the HDF5 based mesh format of Salome and Code_Aster. See +https://docs.salome-platform.org/latest/dev/MEDCoupling/developer/med-file.html + +Families and groups, in DOLFINx terms: + +- Every entity stores one *family* id (``MAI//FAM`` for elements, negative; + ``NOE/FAM`` for nodes, positive; ``0`` for none). Per element type, these ids are + the values of a :class:`dolfinx.mesh.MeshTags`, with family 0 meaning untagged. +- ``FAS/`` maps each family id to the named *groups* it belongs to, which + ``MeshTags`` has no equivalent of. +- A family is one combination of groups, so overlapping groups (an entity in both + ``inlet`` and ``wall``) get a family of their own. A group is the union of its + families, read as a separate ``MeshTags`` with value 1. + +Element family ids are negated when read, as integral subdomain ids in DOLFINx/FFCx +must be non-negative. + +SPDX License identifier: MIT + +Copyright: Jørgen S. Dokken and Simula Research Laboratory +""" + +from pathlib import Path +from typing import Any, Literal + +from mpi4py import MPI + +import basix +import dolfinx +import h5py +import numpy as np +import numpy.typing as npt + +from ...structures import ArrayData, FunctionData, MeshData, MeshTagsData, ReadMeshData +from ...utils import check_file_exists, compute_local_range +from .. import DEFAULT_MESH_NAME, FileMode, ReadMode, get_mesh_name +from ..h5py.backend import h5pyfile, read_component_major_rows + +read_mode = ReadMode.parallel + +# MED element type -> (DOLFINx cell type, degree) +_med_cell_types = { + "PO1": (dolfinx.mesh.CellType.point, 1), + "SE2": (dolfinx.mesh.CellType.interval, 1), + "SE3": (dolfinx.mesh.CellType.interval, 2), + "TR3": (dolfinx.mesh.CellType.triangle, 1), + "TR6": (dolfinx.mesh.CellType.triangle, 2), + "TR7": (dolfinx.mesh.CellType.triangle, 2), + "QU4": (dolfinx.mesh.CellType.quadrilateral, 1), + "QU8": (dolfinx.mesh.CellType.quadrilateral, 2), + "QU9": (dolfinx.mesh.CellType.quadrilateral, 2), + "TE4": (dolfinx.mesh.CellType.tetrahedron, 1), + "T10": (dolfinx.mesh.CellType.tetrahedron, 2), + "PY5": (dolfinx.mesh.CellType.pyramid, 1), + "P13": (dolfinx.mesh.CellType.pyramid, 2), + "PE6": (dolfinx.mesh.CellType.prism, 1), + "P15": (dolfinx.mesh.CellType.prism, 2), + "P18": (dolfinx.mesh.CellType.prism, 2), + "HE8": (dolfinx.mesh.CellType.hexahedron, 1), + "H20": (dolfinx.mesh.CellType.hexahedron, 2), + "H27": (dolfinx.mesh.CellType.hexahedron, 2), +} + +# Node permutation ``p`` with ``dolfinx[i] = med[p[i]]`` for each supported MED element +# type. Derived from MEDCoupling's (LGPL) MED cell model: the nodes it adds in +# ``MEDCouplingUMesh.convertLinearCellsToQuadratic`` and the face centres of +# ``buildDescendingConnectivity`` (for ``P18``), matched to the Basix reference points. +# ``test_med_permutations_match_medcoupling`` checks the table when medcoupling is +# installed. +_med_to_dolfinx = { + "SE2": [0, 1], + "SE3": [0, 1, 2], + "TR3": [0, 1, 2], + "TR6": [0, 1, 2, 4, 5, 3], + "QU4": [0, 1, 3, 2], + "QU8": [0, 1, 3, 2, 4, 7, 5, 6], + "QU9": [0, 1, 3, 2, 4, 7, 5, 6, 8], + "TE4": [0, 1, 2, 3], + "T10": [0, 1, 2, 3, 9, 8, 5, 7, 6, 4], + "PY5": [0, 1, 3, 2, 4], + "PE6": [0, 1, 2, 3, 4, 5], + "P18": [0, 1, 2, 3, 4, 5, 6, 8, 12, 7, 13, 14, 9, 11, 10, 15, 17, 16], + "HE8": [0, 1, 3, 2, 4, 5, 7, 6], + "H20": [0, 1, 3, 2, 4, 5, 7, 6, 8, 11, 16, 9, 17, 10, 19, 18, 12, 15, 13, 14], + "H27": [0, 1, 3, 2, 4, 5, 7, 6, 8, 11, 16, 9, 17, 10, 19, 18, 12, 15, 13, 14] + + [20, 21, 24, 22, 23, 25, 26], +} + +# Layouts with the nodal serendipity coordinate element +_med_serendipity = ("QU8", "H20") + + +def geometry_layout(med_type: str) -> tuple[dolfinx.mesh.CellType, int, bool]: + """Get the coordinate element layout of a MED element type. + + Args: + med_type: MED element type, e.g. ``"T10"`` + + Returns: + Tuple containing: + - The DOLFINx cell type + - The degree of the coordinate element + - True if the coordinate element is the nodal serendipity element + (``QU8``, ``H20``), False if it is Lagrange + """ + cell_type, degree = _med_cell_types[med_type] + return cell_type, degree, med_type in _med_serendipity + + +def get_mesh_step_and_families( + h5file: h5py.File, + filename: Path | str, + backend_args: dict[str, Any] | None, + time: str | float | None, +) -> tuple[h5py.Group, h5py.Group]: + """Get the groups holding the nodes, connectivity and families of a mesh. + + In DOLFINx terms, ``NOE`` (*nœuds*) holds the geometry nodes, i.e. the coordinates + in ``mesh.geometry.x``, and ``MAI`` (*mailles*) the connectivity: one block per + element type, listing the nodes of each cell, facet, edge or point element. + + A time-dependent mesh has one group per computation step, named by its step and + iteration numbers so that the names sort in time order. A step only stores data + that changed, and links to the earlier step for the rest. + + Args: + h5file: The MED file + filename: Path to the file, used in error messages + backend_args: Arguments to backend. ``"name"`` selects the mesh; if it is not + given and the file holds a single mesh, that mesh is used. + time: Time of the step to read. Ignored for a mesh with a single step, and + ``None`` selects the first step. + + Returns: + Tuple containing: + - The step group, holding ``NOE`` and ``MAI`` + - The families of the mesh (``FAS/``) + """ + if "ENS_MAA" not in h5file.keys(): + raise ValueError(f"Could not find a set of meshes (ENS_MAA) in {filename}.") + meshes = sorted(h5file["ENS_MAA"].keys()) + name = get_mesh_name(backend_args) + if name not in meshes: + if name == DEFAULT_MESH_NAME and len(meshes) == 1: + name = meshes[0] + else: + raise KeyError(f"Could not find mesh '{name}' in {filename}. Available: {meshes}") + mesh = h5file["ENS_MAA"][name] + # Mesh type 0 is unstructured, 1 is a structured grid without connectivity + if mesh.attrs.get("TYP", 0) != 0: + raise NotImplementedError(f"Structured MED meshes are not supported ('{name}').") + if "FAS" not in h5file.keys() or name not in h5file["FAS"].keys(): + raise ValueError(f"Could not find the families (FAS) of mesh '{name}' in {filename}.") + + if "NOE" in mesh.keys(): + step = mesh + else: + steps = [mesh[key] for key in sorted(mesh.keys())] + if len(steps) == 1 or time is None: + step = steps[0] + else: + times = np.array([float(s.attrs["PDT"]) for s in steps]) + match = np.flatnonzero(np.isclose(times, float(time))) + if len(match) != 1: + raise KeyError( + f"Expected one step of mesh '{name}' at time {time} in {filename}, found " + f"{len(match)}. Available times: {times.tolist()}" + ) + step = steps[match[0]] + if unsupported := sorted(set(step["MAI"].keys()) - _med_cell_types.keys()): + raise NotImplementedError(f"MED element types {unsupported} are not supported.") + return step, h5file["FAS"][name] + + +def read_entities( + element_blocks: h5py.Group, dim: int, comm: MPI.Comm +) -> tuple[str, list[int], npt.NDArray[np.int64]] | None: + """Read the entities of one topological dimension, distributed over all processes. + + Args: + element_blocks: The group holding one connectivity block per element type + (``MAI`` in the MED file) + dim: Topological dimension of the entities + comm: MPI communicator to distribute the entities over + + Returns: + ``None`` if the mesh has no entities of dimension ``dim``, otherwise a tuple + containing: + - The MED element type of the entities, e.g. ``"TR3"`` + - The range ``[start, end)`` of entities read on this process + - Their connectivity as global input node indices in DOLFINx node order, + of shape ``(end - start, num_nodes_per_entity)`` + """ + blocks = sorted( + key + for key in element_blocks.keys() + if dolfinx.mesh.cell_dim(_med_cell_types[key][0]) == dim + ) + if len(blocks) == 0: + return None + if len(blocks) > 1: + raise NotImplementedError( + f"Mixed element types {blocks} of dimension {dim} are not supported." + ) + med_type = blocks[0] + block = element_blocks[med_type] + + # Read a contiguous chunk of the entities on each process + num_entities_global = num_entities(block) + local_range = compute_local_range(comm, num_entities_global) + nodes = read_component_major_rows(block["NOD"], num_entities_global, local_range) + + # MED node indices are 1-based + connectivity = nodes.astype(np.int64) - 1 + + # Permute the nodes of each entity from MED to DOLFINx ordering + if med_type not in _med_to_dolfinx: + raise NotImplementedError( + f"MED element type {med_type} has no Lagrange or nodal serendipity layout." + ) + perm = _med_to_dolfinx[med_type] + if connectivity.shape[1] != len(perm): + raise ValueError(f"Expected {len(perm)} nodes per {med_type}, got {connectivity.shape[1]}.") + return med_type, local_range, connectivity[:, perm] + + +def read_entity_families(block: h5py.Group, local_range: list[int]) -> npt.NDArray[np.int32]: + """Read the family id of each entity of one MED element type on this process. + + Each entity has exactly one family, so this is a tag value per entity. Which + groups a family belongs to is read with :func:`read_family_groups`. + + Args: + block: Connectivity block of the element type (``MAI/``) + local_range: Range ``[start, end)`` of entities to read on this process + + Returns: + Family ids of shape ``(end - start,)``. A block without ``FAM`` has every + entity in family 0, i.e. in no group. + """ + if "FAM" in block.keys(): + return block["FAM"][slice(*local_range)].astype(np.int32) + return np.zeros(local_range[1] - local_range[0], dtype=np.int32) + + +def num_entities(block: h5py.Group) -> int: + """Get the global number of entities of one MED element type. + + Args: + block: Connectivity block of the element type (``MAI/``) + + Returns: + The number of entities, stored as the ``NBR`` attribute of the connectivity + """ + return int(block["NOD"].attrs["NBR"]) + + +def read_family_groups( + families: h5py.Group, kind: Literal["ELEME", "NOEUD"] +) -> dict[int, list[str]]: + """Read which named groups each family of a mesh belongs to. + + A family is one combination of groups, so one family can belong to several + groups. The family of each entity is read with :func:`read_entity_families`. + + Args: + families: The families of the mesh (``FAS/``) + kind: ``"ELEME"`` (*éléments*) for element families, ``"NOEUD"`` (*nœuds*) + for node families + + Returns: + Map from family id to the names of its groups. Empty if the mesh has no + groups of this kind, in which case MED does not write ``kind``. + """ + family_groups: dict[int, list[str]] = {} + if kind not in families.keys(): + return family_groups + for family in families[kind].values(): + groups = [] + if "GRO" in family.keys(): + for row in family["GRO"]["NOM"][()]: + raw = np.asarray(row).astype(np.uint8).tobytes() + groups.append(raw.split(b"\x00")[0].decode().rstrip()) + family_groups[int(family.attrs["NUM"])] = groups + return family_groups + + +def read_groups( + filename: Path | str, comm: MPI.Comm, backend_args: dict[str, Any] | None = None +) -> dict[str, list[int]]: + """Read the element groups of a MED mesh and the families that make them up. + + Use this to find the group names accepted by :func:`read_meshtags_data`, or to + decode the (negated) family ids returned for ``"cell"`` and ``"facet"``. + + Args: + filename: Path to file + comm: MPI communicator used in storage + backend_args: Arguments to backend, ``"name"`` selects the mesh + + Returns: + Map from group name to the negated family ids that belong to the group + """ + check_file_exists(filename) + with h5pyfile(filename, "r", comm=comm) as h5file: + _, families = get_mesh_step_and_families(h5file, filename, backend_args, None) + family_groups = read_family_groups(families, "ELEME") + groups: dict[str, list[int]] = {} + for family, names in sorted(family_groups.items(), reverse=True): + for group in names: + groups.setdefault(group, []).append(-family) + return groups + + +def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]: + args = arguments or {} + return args + + +def write_attributes( + filename: Path | str, + comm: MPI.Comm, + name: str, + attributes: dict[str, np.ndarray], + backend_args: dict[str, Any] | None = None, +): + """Write attributes to file. + + Args: + filename: Path to file to write to + comm: MPI communicator used in storage + name: Name of the attributes + attributes: Dictionary of attributes to write to file + backend_args: Arguments to backend + """ + raise NotImplementedError("The MED backend cannot write attributes.") + + +def read_attributes( + filename: Path | str, + comm: MPI.Comm, + name: str, + backend_args: dict[str, Any] | None = None, +) -> dict[str, Any]: + """Read attributes from file. + + Args: + filename: Path to file to read from + comm: MPI communicator used in storage + name: Name of the attributes + backend_args: Arguments to backend + Returns: + The attributes + """ + raise NotImplementedError("The MED backend cannot read attributes.") + + +def snapshot_checkpoint( + filename: Path | str, + mode: FileMode, + u: dolfinx.fem.Function, + backend_args: dict[str, Any] | None, +): + """Create a snapshot checkpoint of a dolfinx function. + + Args: + filename: Path to file to read from + mode: File-mode to store the function + u: dolfinx function to create a snapshot checkpoint for + backend_args: Arguments to backend + """ + raise NotImplementedError("The MED backend cannot make checkpoints.") + + +def read_hdf5_array( + comm: MPI.Comm, + filename: Path | str, + group: str, + backend_args: dict[str, Any] | None, +) -> tuple[np.ndarray, int]: + """Read an array from an HDF5 file. + + Args: + comm: MPI communicator used in storage + filename: Path to file to read from + group: Group in HDF5 file where array is stored + backend_args: Arguments to backend + + Returns: + Tuple containing: + - Numpy array read from file + - Global starting point on the process. + Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc. + """ + raise NotImplementedError("The MED backend cannot read HDF5 arrays.") + + +def read_timestamps( + filename: Path | str, + comm: MPI.Comm, + function_name: str, + backend_args: dict[str, Any] | None = None, +) -> npt.NDArray[np.float64 | str]: # type: ignore[type-var] + """Read time-stamps from a checkpoint file. + + Args: + filename: Path to file + comm: MPI communicator + function_name: Name of the function to read time-stamps for + backend_args: Arguments for backend, for instance file type. + + Returns: + The time-stamps + """ + raise NotImplementedError("The MED backend cannot read timestamps.") + + +def write_mesh( + filename: Path | str, + comm: MPI.Comm, + mesh: MeshData, + backend_args: dict[str, Any] | None = None, + mode: FileMode = FileMode.write, + time: float = 0.0, +): + """Write a mesh to file. + + Args: + filename: Path to file to write to. + mesh: Internal data structure for the mesh data to save to file + comm: MPI communicator used in storage + backend_args: Arguments to backend + mode: Mode to use (write or append) + time: Time stamp + """ + raise NotImplementedError("The MED backend cannot write meshes.") + + +def read_mesh_data( + filename: Path | str, + comm: MPI.Comm, + time: str | float | None = 0.0, + read_from_partition: bool = False, + backend_args: dict[str, Any] | None = None, +) -> ReadMeshData: + """Read mesh data from a MED file. + + Only cells of the highest topological dimension are read; these must all be of + the same type. Lagrange layouts up to degree 2 (e.g. ``TR6``, ``T10``, ``H27``) and + the serendipity ``QU8`` and ``H20`` are supported. Node and cell data is read in + contiguous chunks on each process. + + Args: + filename: Path to input file + comm: The MPI communciator to distribute the mesh over + time: Time of the step to read from a time-dependent mesh. Ignored for a mesh + with a single step, and ``None`` selects the first step. + read_from_partition: Read mesh with partition from file (not supported) + backend_args: Arguments to backend, ``"name"`` selects the mesh + Returns: + The mesh topology, geometry, UFL domain and partition function + """ + if read_from_partition: + raise NotImplementedError("The MED backend cannot read partitioning data.") + check_file_exists(filename) + with h5pyfile(filename, "r", comm=comm) as h5file: + step, _ = get_mesh_step_and_families(h5file, filename, backend_args, time) + coo = step["NOE"]["COO"] + num_nodes = int(coo.attrs["NBR"]) + x = read_component_major_rows(coo, num_nodes, compute_local_range(comm, num_nodes)) + + element_blocks = step["MAI"] + tdim = max(dolfinx.mesh.cell_dim(_med_cell_types[key][0]) for key in element_blocks.keys()) + entities = read_entities(element_blocks, tdim, comm) + assert entities is not None + med_type, _, cells = entities + cell_type, degree, serendipity = geometry_layout(med_type) + return ReadMeshData( + cells=cells, + cell_type=dolfinx.mesh.to_string(cell_type), + x=x, + lvar=int(basix.LagrangeVariant.equispaced), + degree=degree, + partition_graph=None, + serendipity=serendipity, + ) + + +def write_meshtags( + filename: str | Path, + comm: MPI.Comm, + data: MeshTagsData, + backend_args: dict[str, Any] | None = None, +): + """Write mesh tags to file. + + Args: + filename: Path to file to write to + comm: MPI communicator used in storage + data: Internal data structure for the mesh tags to save to file + backend_args: Arguments to backend + """ + raise NotImplementedError("The MED backend cannot write meshtags.") + + +def read_meshtags_data( + filename: str | Path, + comm: MPI.Comm, + name: str, + backend_args: dict[str, Any] | None = None, +) -> MeshTagsData: + """Read mesh tags from a MED file. + + Args: + filename: Path to file to read from + comm: MPI communicator used in storage + name: ``"cell"`` or ``"facet"`` for the entities of the highest topological + dimension or one below, otherwise the name of a MED group (see + :func:`read_groups`) + backend_args: Arguments to backend, ``"name"`` selects the mesh and ``"time"`` the + step of a time-dependent mesh (the first step if not given) + + Returns: + Internal data structure for the mesh tags read from file. For ``"cell"`` and + ``"facet"``, each entity is tagged with its negated MED family id, and entities + without a family (id 0) are skipped. For a group, its entities are tagged with + ``1``. + """ + check_file_exists(filename) + with h5pyfile(filename, "r", comm=comm) as h5file: + time = None if backend_args is None else backend_args.get("time") + step, families = get_mesh_step_and_families(h5file, filename, backend_args, time) + element_blocks = step["MAI"] + keys = sorted(element_blocks.keys()) + tdim = max(dolfinx.mesh.cell_dim(_med_cell_types[key][0]) for key in keys) + group_families: list[int] | None = None + if name in ("cell", "facet"): + dim = tdim if name == "cell" else tdim - 1 + else: + family_groups = read_family_groups(families, "ELEME") + group_families = [f for f, groups in family_groups.items() if name in groups] + if len(group_families) == 0: + nodes = read_family_groups(families, "NOEUD") + if any(name in groups for groups in nodes.values()): + raise NotImplementedError(f"Group '{name}' only holds nodes.") + available = sorted({g for groups in family_groups.values() for g in groups}) + raise KeyError( + f"Could not find group '{name}' in {filename}. Available: {available}" + ) + + # Find the blocks holding entities of the group. FAS does not record which + # element type a family is used on, and each process only reads its chunk + # of every block, so the reduction is needed for all processes to agree on + # the blocks and dimension before the collective meshtag creation. + in_block = np.zeros(len(keys), dtype=np.int32) + for i, key in enumerate(keys): + r = compute_local_range(comm, num_entities(element_blocks[key])) + in_block[i] = np.isin( + read_entity_families(element_blocks[key], r), group_families + ).any() + comm.Allreduce(MPI.IN_PLACE, in_block, op=MPI.MAX) + blocks = [key for key, found in zip(keys, in_block) if found] + dims = {dolfinx.mesh.cell_dim(_med_cell_types[key][0]) for key in blocks} + if len(dims) != 1: + raise RuntimeError( + f"Group '{name}' must hold entities of a single dimension, got {dims}." + ) + dim = dims.pop() + + entities = read_entities(element_blocks, dim, comm) + if entities is None: + # No entities of this dimension; read_meshtags sets the number of nodes + return MeshTagsData( + name=name, + values=np.zeros(0, dtype=np.int32), + indices=np.zeros((0, 0), dtype=np.int64), + dim=dim, + ) + med_type, local_range, connectivity = entities + family = read_entity_families(element_blocks[med_type], local_range) + + if group_families is None: + keep = family != 0 + values = -family[keep] + else: + keep = np.isin(family, group_families) + values = np.ones(np.count_nonzero(keep), dtype=np.int32) + return MeshTagsData(name=name, values=values, indices=connectivity[keep], dim=dim) + + +def read_dofmap( + filename: str | Path, + comm: MPI.Comm, + name: str, + backend_args: dict[str, Any] | None, +) -> dolfinx.graph.AdjacencyList: + """Read the dofmap of a function with a given name. + + Args: + filename: Path to file to read from + comm: MPI communicator used in storage + name: Name of the function to read the dofmap for + backend_args: Arguments to backend + + Returns: + Dofmap as an :class:`dolfinx.graph.AdjacencyList` + """ + raise NotImplementedError("The MED backend cannot read dofmap.") + + +def read_dofs( + filename: str | Path, + comm: MPI.Comm, + name: str, + time: float, + backend_args: dict[str, Any] | None, +) -> tuple[npt.NDArray[np.float32 | np.float64 | np.complex64 | np.complex128], int]: + """Read the dofs (values) of a function with a given name from a given timestep. + + Args: + filename: Path to file to read from + comm: MPI communicator used in storage + name: Name of the function to read the dofs for + time: Time stamp associated with the function to read + backend_args: Arguments to backend + + Returns: + Contiguous sequence of degrees of freedom (with respect to input data) + and the global starting point on the process. + Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc. + """ + raise NotImplementedError("The MED backend cannot read dofs.") + + +def read_cell_perms( + comm: MPI.Comm, filename: Path | str, backend_args: dict[str, Any] | None +) -> npt.NDArray[np.uint32]: + """ + Read cell permutation from file with given communicator, + Split in continuous chunks based on number of cells in the input data. + + Args: + comm: MPI communicator used in storage + filename: Path to file to read from + backend_args: Arguments to backend + + Returns: + Contiguous sequence of permutations (with respect to input data) + Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc. + """ + raise NotImplementedError("The MED backend cannot read cell perms.") + + +def write_function( + filename: str | Path, + comm: MPI.Comm, + u: FunctionData, + time: float, + mode: FileMode, + backend_args: dict[str, Any] | None = None, +): + """Write a function to file. + + Args: + filename: Path to file to write to + comm: MPI communicator used in storage + u: Internal data structure for the function data to save to file + time: Time stamp associated with function + mode: File-mode to store the function + backend_args: Arguments to backend + """ + raise NotImplementedError("The MED backend cannot write functions.") + + +def read_legacy_mesh( + filename: Path | str, comm: MPI.Comm, group: str +) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.floating], str | None]: + """Read in the mesh topology, geometry and (optionally) cell type from a + legacy DOLFIN HDF5-file. + + Args: + filename: Path to file to read from + comm: MPI communicator used in storage + group: Group in HDF5 file where mesh is stored + + Returns: + Tuple containing: + - Topology as a (num_cells, num_vertices_per_cell) array of global vertex indices + - Geometry as a (num_vertices, geometric_dimension) array of vertex coordinates + - Cell type as a string (e.g. "tetrahedron") or None if not found + """ + raise NotImplementedError("The MED backend cannot read legacy meshes.") + + +def read_point_data( + filename: Path | str, + name: str, + comm: MPI.Comm, + time: float | str | None, + backend_args: dict[str, Any] | None, +) -> tuple[np.ndarray, int]: + """Read data from the nodes of a mesh. + + Args: + filename: Path to file + name: Name of point data + comm: Communicator to launch IO on. + time: The time stamp + backend_args: The backend arguments + + Returns: + Data local to process (contiguous, no mpi comm) and local start range + """ + raise NotImplementedError("The MED backend cannot read point data.") + + +def read_cell_data( + filename: Path | str, + name: str, + comm: MPI.Comm, + time: str | float | None, + backend_args: dict[str, Any] | None, +) -> tuple[npt.NDArray[np.int64], np.ndarray]: + """Read data from the cells of a mesh. + + Args: + filename: Path to file + name: Name of point data + comm: Communicator to launch IO on. + time: The time stamp + backend_args: The backend arguments + Returns: + A tuple (topology, dofs) where topology contains the + vertex indices of the cells, dofs the degrees of + freedom within that cell. + """ + raise NotImplementedError("The MED backend cannot read cell data.") + + +def read_function_names( + filename: Path | str, comm: MPI.Comm, backend_args: dict[str, Any] | None +) -> list[str]: + """Read all function names from a file. + + Args: + filename: Path to file + comm: MPI communicator to launch IO on. + backend_args: Arguments to backend + + Returns: + A list of function names. + """ + raise NotImplementedError("The MED backend cannot read function names.") + + +def write_data( + filename: Path | str, + array_data: ArrayData, + comm: MPI.Comm, + time: str | float | None, + mode: FileMode, + backend_args: dict[str, Any] | None, +): + """Write a 2D-array to file (distributed across proceses with MPI). + + Args: + filename: Path to file + array_data: Data to write to file + comm: MPI communicator to open the file with + time: Time-stamp for data. + mode: Append or write + backend_args: The backend arguments + """ + raise NotImplementedError("The MED backend cannot write data.") diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index 8f8744f..1822873 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -275,8 +275,14 @@ def read_meshtags( backend_args = backend_cls.get_default_backend_args(backend_args) data = backend_cls.read_meshtags_data(filename, mesh.comm, meshtag_name, backend_args) + indices = data.indices + if indices.shape[0] == 0: + # A backend without entities may not know the number of geometry nodes per entity + layout = fenicsx_compat.cmap(mesh).create_dof_layout() + num_nodes = len(layout.entity_closure_dofs(int(data.dim), 0)) + indices = np.zeros((0, num_nodes), dtype=np.int64) local_entities, local_values = dolfinx.io.distribute_entity_data( - mesh, int(data.dim), data.indices, data.values + mesh, int(data.dim), indices, data.values ) mesh.topology.create_connectivity(data.dim, 0) mesh.topology.create_connectivity(data.dim, mesh.topology.dim) @@ -492,14 +498,30 @@ def read_mesh( ) # Create DOLFINx mesh - element = basix.ufl.element( - basix.ElementFamily.P, - dist_in_data.cell_type, - dist_in_data.degree, - basix.LagrangeVariant(int(dist_in_data.lvar)), - shape=(dist_in_data.x.shape[1],), - dtype=dist_in_data.x.dtype, - ) + element: basix.ufl._ElementBase + if dist_in_data.serendipity: + try: + from dolfinx.io.utils import _serendipity_coordinate_element + except ImportError as e: + raise NotImplementedError( + "Meshes with serendipity geometry (8-node quadrilaterals, 20-node " + "hexahedra) require the main branch of DOLFINx." + ) from e + element = basix.ufl.blocked_element( + _serendipity_coordinate_element( + dolfinx.mesh.to_type(dist_in_data.cell_type), dist_in_data.x.dtype + ), + shape=(dist_in_data.x.shape[1],), + ) + else: + element = basix.ufl.element( + basix.ElementFamily.P, + dist_in_data.cell_type, + dist_in_data.degree, + basix.LagrangeVariant(int(dist_in_data.lvar)), + shape=(dist_in_data.x.shape[1],), + dtype=dist_in_data.x.dtype, + ) domain = ufl.Mesh(element) # `dolfinx.mesh.PartitioningFunc` is not available in all supported # versions of DOLFINx, and the accepted callback signature varies (see diff --git a/src/io4dolfinx/structures.py b/src/io4dolfinx/structures.py index 5d33f72..e19e97e 100644 --- a/src/io4dolfinx/structures.py +++ b/src/io4dolfinx/structures.py @@ -78,6 +78,9 @@ class ReadMeshData: degree: int #: The degree of the underlying Lagrange element #: Partitioning information per cell on the process partition_graph: AdjacencyList | None = None + #: Geometry uses the degree 2 serendipity element with nodes at vertices and edge + #: midpoints (8-node quadrilateral, 20-node hexahedron) instead of Lagrange + serendipity: bool = False @dataclass diff --git a/tests/test_backends.py b/tests/test_backends.py index f62a629..9055b5c 100644 --- a/tests/test_backends.py +++ b/tests/test_backends.py @@ -5,7 +5,7 @@ import pytest import io4dolfinx -from io4dolfinx.backends import BUILTIN_BAKENDS +from io4dolfinx.backends import BUILTIN_BACKENDS @pytest.fixture(autouse=True) @@ -93,7 +93,7 @@ def test_list_builtin_backends(): # Verify that all returned backends are recognized as built-in backends for backend in available_backends: assert isinstance(backend, str) - assert backend in BUILTIN_BAKENDS + assert backend in BUILTIN_BACKENDS # We can be reasonably certain that 'h5py' or 'adios2' should be # present if the io4dolfinx test suite is running successfully diff --git a/tests/test_med.py b/tests/test_med.py new file mode 100644 index 0000000..bd6387f --- /dev/null +++ b/tests/test_med.py @@ -0,0 +1,784 @@ +import urllib.error +import urllib.request +from enum import Enum + +from mpi4py import MPI + +import basix +import dolfinx +import dolfinx.io.utils +import fenicsx_compat +import numpy as np +import pytest +import ufl + +import io4dolfinx + +h5py = pytest.importorskip("h5py") + +from io4dolfinx.backends.med import backend as med_backend # noqa: E402 +from io4dolfinx.backends.med.backend import read_groups # noqa: E402 + +requires_serendipity = pytest.mark.skipif( + not hasattr(dolfinx.io.utils, "_serendipity_coordinate_element"), + reason="Serendipity geometry requires the main branch of DOLFINx", +) + +# (DOLFINx cell type, geometry degree, serendipity) -> MED element type +_dolfinx_to_med = { + ("interval", 1, False): "SE2", + ("interval", 2, False): "SE3", + ("interval", 2, True): "SE3", + ("triangle", 1, False): "TR3", + ("triangle", 2, False): "TR6", + ("quadrilateral", 1, False): "QU4", + ("quadrilateral", 2, False): "QU9", + ("quadrilateral", 2, True): "QU8", + ("tetrahedron", 1, False): "TE4", + ("tetrahedron", 2, False): "T10", + ("prism", 1, False): "PE6", + ("prism", 2, False): "P18", + ("hexahedron", 1, False): "HE8", + ("hexahedron", 2, False): "H27", + ("hexahedron", 2, True): "H20", +} + +# MED node layouts, read off the Salome connectivity diagrams +# (https://docs.salome-platform.org/latest/gui/SMESH/connectivity.html) independently +# of the permutations in the backend: the DOLFINx vertices that are the MED vertices, +# and for each further node the MED vertices it is the centroid of. +_hex_edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4)] +_hex_edges += [(0, 4), (1, 5), (2, 6), (3, 7)] +_hex_faces: list[tuple[int, ...]] = [ + (0, 1, 2, 3), + (0, 1, 5, 4), + (1, 2, 6, 5), + (2, 3, 7, 6), + (3, 0, 4, 7), +] +_hex_faces += [(4, 5, 6, 7), tuple(range(8))] +_prism_edges = [(0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (0, 3), (1, 4), (2, 5)] +_quad_edges = [(0, 1), (1, 2), (2, 3), (3, 0)] +_med_layouts = { + "SE2": ([0, 1], []), + "SE3": ([0, 1], [(0, 1)]), + "TR3": ([0, 1, 2], []), + "TR6": ([0, 1, 2], [(0, 1), (1, 2), (2, 0)]), + "QU4": ([0, 1, 3, 2], []), + "QU8": ([0, 1, 3, 2], _quad_edges), + "QU9": ([0, 1, 3, 2], _quad_edges + [(0, 1, 2, 3)]), + "TE4": ([0, 1, 2, 3], []), + "T10": ([0, 1, 2, 3], [(0, 1), (1, 2), (0, 2), (0, 3), (1, 3), (2, 3)]), + "PE6": ([0, 1, 2, 3, 4, 5], []), + # The diagram stops at PE15; the face nodes follow the side faces in order + "P18": ([0, 1, 2, 3, 4, 5], _prism_edges + [(0, 1, 4, 3), (1, 2, 5, 4), (2, 0, 3, 5)]), + "HE8": ([0, 1, 3, 2, 4, 5, 7, 6], []), + "H20": ([0, 1, 3, 2, 4, 5, 7, 6], _hex_edges), + "H27": ([0, 1, 3, 2, 4, 5, 7, 6], _hex_edges + _hex_faces), +} + +_festim_url = ( + "https://raw.githubusercontent.com/festim-dev/FESTIM-workshop/" + "e2f171484578af33ea1c4030016a5a600768a1ec/book/content/meshes/salome/Mesh_1.med" +) + + +class DownloadStatus(Enum): + success = 1 + failed = -1 + no_connection = -2 + + +def download_file_if_not_exists( + url, filename, comm: MPI.Intracomm = MPI.COMM_WORLD, rank: int = 0 +) -> DownloadStatus: + status = DownloadStatus.failed + if comm.rank == rank: + if not filename.exists(): + try: + urllib.request.urlretrieve(url, filename) + status = DownloadStatus.success + except urllib.error.URLError as e: + if str(e) == "": + status = DownloadStatus.no_connection + else: + status = DownloadStatus.failed + else: + status = DownloadStatus.success + status = comm.bcast(status, root=rank) + comm.Barrier() + return status + + +def _fixture_mesh(cell_type: dolfinx.mesh.CellType, degree: int = 1) -> dolfinx.mesh.Mesh: + """Serial mesh of the unit square/cube, with geometry of the given degree.""" + if dolfinx.mesh.cell_dim(cell_type) == 2: + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_SELF, 4, 4, cell_type) + else: + mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_SELF, 2, 2, 2, cell_type) + if degree > 1: + if hasattr(dolfinx.fem, "interpolate_geometry"): + cmap = dolfinx.fem.coordinate_element(cell_type, degree, dtype=mesh.geometry.x.dtype) + mesh = dolfinx.fem.interpolate_geometry(mesh, cmap) + else: + mesh = fenicsx_compat.reconstruct_mesh(mesh, degree) + return mesh + + +def _to_med(nodes, x, med_type): + """MED connectivity of entities given by their DOLFINx geometry nodes. + + Nodes are found by position, so ``x`` must be undeformed.""" + vertices, centroids = _med_layouts[med_type] + med_vertices = nodes[:, vertices] + columns = [med_vertices] + rows = np.arange(nodes.shape[0]) + for centroid in centroids: + target = x[med_vertices[:, centroid]].mean(axis=1) + distance = np.linalg.norm(x[nodes] - target[:, None, :], axis=2) + closest = distance.argmin(axis=1) + assert np.allclose(distance[rows, closest], 0.0) + columns.append(nodes[rows, closest][:, None]) + return np.hstack(columns) + + +def _deform(x): + """Quadratic map, represented exactly by quadratic Lagrange and serendipity geometry.""" + y = x.copy() + for i in range(x.shape[1]): + y[:, i] += 0.1 * x[:, (i + 1) % x.shape[1]] ** 2 + return y + + +def _med_blocks(mesh, cell_family, facet_markers=(), serendipity=False, deform=None): + """Coordinates and MED blocks ``{type: (0-based nodes, family)}`` of a serial mesh. + + ``cell_family`` maps cell midpoints to families and ``facet_markers`` is a sequence of + ``(boundary marker, family)``. Serendipity drops the face and interior nodes of a + quadratic mesh. ``deform`` maps the coordinates after the connectivity is built.""" + tdim = mesh.topology.dim + x = mesh.geometry.x[:, : mesh.geometry.dim].copy() + degree = fenicsx_compat.cmap(mesh).degree + cell_type = dolfinx.mesh.to_string(mesh.topology.cell_type) + dofmap = fenicsx_compat.dofmap(mesh) + med_type = _dolfinx_to_med[(cell_type, degree, serendipity)] + blocks = {med_type: (_to_med(dofmap, x, med_type), cell_family(x[dofmap].mean(axis=1)))} + if len(facet_markers) > 0: + mesh.topology.create_connectivity(tdim - 1, tdim) + facets, families = [], [] + for marker, family in facet_markers: + entities = dolfinx.mesh.locate_entities_boundary(mesh, tdim - 1, marker) + facets.append(entities) + families.append(np.full(len(entities), family, dtype=np.int32)) + nodes = dolfinx.mesh.entities_to_geometry(mesh, tdim - 1, np.hstack(facets)) + facet_type = dolfinx.cpp.mesh.cell_entity_type(mesh.topology.cell_type, tdim - 1, 0) + med_type = _dolfinx_to_med[(dolfinx.mesh.to_string(facet_type), degree, serendipity)] + blocks[med_type] = (_to_med(nodes, x, med_type), np.hstack(families)) + + # Drop nodes no entity uses + used = np.unique(np.hstack([nodes.ravel() for nodes, _ in blocks.values()])) + new_index = np.full(x.shape[0], -1, dtype=np.int64) + new_index[used] = np.arange(len(used)) + blocks = {key: (new_index[nodes], family) for key, (nodes, family) in blocks.items()} + x = x[used] + return (x if deform is None else deform(x)), blocks + + +def _write_med( + filename, + x, + blocks, + families=None, + node_families=None, + mesh_name="mesh", + mode="w", + node_num=None, + mesh_type=0, + time=0.0, + later_steps=(), +): + """Write a MED file (call on a single rank). + + ``mesh_type`` is the MED mesh type, 0 for unstructured and 1 for structured. The + first step is at ``time``. ``later_steps`` is a sequence of ``(step number, time, + coordinates)``; as MED does, a later step stores its coordinates and links its + connectivity and families to the previous step.""" + with h5py.File(filename, mode) as f: + info = f.require_group("INFOS_GENERALES") + for key, value in zip(("MAJ", "MIN", "REL"), (4, 1, 1)): + info.attrs[key] = np.int32(value) + mesh = f.require_group("ENS_MAA").create_group(mesh_name) + mesh.attrs["ESP"] = np.int32(x.shape[1]) + mesh.attrs["TYP"] = np.int32(mesh_type) + previous = None + for number, step_time, step_x in [(-1, time, x), *later_steps]: + # Steps are named by step and iteration number, the first being (-1, -1) + step = mesh.create_group(f"{number:020d}{-1:020d}") + step.attrs["NDT"] = np.int32(number) + step.attrs["NOR"] = np.int32(-1) + step.attrs["PDT"] = np.float64(step_time) + noe = step.create_group("NOE") + coo = noe.create_dataset("COO", data=step_x.flatten(order="F")) + coo.attrs["NBR"] = np.int32(step_x.shape[0]) + noe.create_dataset("FAM", data=np.zeros(step_x.shape[0], dtype=np.int32)) + if node_num is not None: + noe.create_dataset("NUM", data=node_num.astype(np.int32)) + mai = step.create_group("MAI") + for med_type, (nodes, family) in blocks.items(): + block = mai.create_group(med_type) + if previous is not None: + for name in ("NOD", "FAM"): + block[name] = h5py.SoftLink(f"{previous}/MAI/{med_type}/{name}") + continue + nod = block.create_dataset( + "NOD", data=(nodes + 1).flatten(order="F").astype(np.int32) + ) + nod.attrs["NBR"] = np.int32(nodes.shape[0]) + block.create_dataset("FAM", data=family.astype(np.int32)) + previous = step.name + + fas = f.require_group("FAS").create_group(mesh_name) + fas.create_group("FAMILLE_ZERO").attrs["NUM"] = np.int32(0) + for kind, kind_families in (("ELEME", families), ("NOEUD", node_families)): + if kind_families is None: + continue + group = fas.create_group(kind) + for family_id, names in kind_families.items(): + family = group.create_group(f"FAM_{family_id}_{'_'.join(names)}") + family.attrs["NUM"] = np.int32(family_id) + gro = family.create_group("GRO") + gro.attrs["NBR"] = np.int32(len(names)) + nom = np.zeros((len(names), 80), dtype=np.int8) + for i, name in enumerate(names): + encoded = name.encode() + nom[i, : len(encoded)] = np.frombuffer(encoded, dtype=np.int8) + gro.create_dataset("NOM", data=nom) + + +def _on_root(comm, function, *args, **kwargs): + if comm.rank == 0: + function(*args, **kwargs) + comm.Barrier() + + +def _integrate(mesh, measure, tag=None, meshtags=None): + m = ufl.Measure(measure, domain=mesh, subdomain_data=meshtags) + m = m if tag is None else m(tag) + value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1 * m)) + return mesh.comm.allreduce(value, op=MPI.SUM) + + +def _left_right(midpoints): + return np.where(midpoints[:, 0] < 0.5, -2, -3).astype(np.int32) + + +_families = {-2: ["left"], -3: ["right"], -4: ["inlet"], -5: ["outlet"]} +_facet_markers = ( + (lambda x: np.isclose(x[0], 0.0), -4), + (lambda x: np.isclose(x[0], 1.0), -5), + (lambda x: np.isclose(x[1], 0.0), 0), +) + + +@pytest.fixture +def tmp_dir(tmp_path): + return MPI.COMM_WORLD.bcast(tmp_path, root=0) + + +@pytest.mark.parametrize( + "cell_type", + [ + dolfinx.mesh.CellType.triangle, + dolfinx.mesh.CellType.quadrilateral, + dolfinx.mesh.CellType.tetrahedron, + dolfinx.mesh.CellType.hexahedron, + ], +) +def test_read_mesh_and_meshtags(tmp_dir, cell_type): + comm = MPI.COMM_WORLD + filename = tmp_dir / "mesh.med" + ref = _fixture_mesh(cell_type) + x, blocks = _med_blocks(ref, _left_right, _facet_markers) + _on_root(comm, _write_med, filename, x, blocks, _families) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + tdim = mesh.topology.dim + assert mesh.topology.cell_type == cell_type + assert mesh.topology.index_map(tdim).size_global == ref.topology.index_map(tdim).size_global + assert mesh.geometry.index_map().size_global == x.shape[0] + assert np.isclose(_integrate(mesh, "dx"), 1.0) + + ct = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="cell", backend="med") + assert np.isclose(_integrate(mesh, "dx", 2, ct), 0.5) + assert np.isclose(_integrate(mesh, "dx", 3, ct), 0.5) + + ft = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="facet", backend="med") + assert ft.dim == tdim - 1 + assert np.isclose(_integrate(mesh, "ds", 4, ft), 1.0) + assert np.isclose(_integrate(mesh, "ds", 5, ft), 1.0) + num_untagged = comm.allreduce(np.count_nonzero(ft.values == 0), op=MPI.SUM) + assert num_untagged == 0 + + left = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="left", backend="med") + assert left.dim == tdim + assert np.isclose(_integrate(mesh, "dx", 1, left), 0.5) + outlet = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="outlet", backend="med") + assert outlet.dim == tdim - 1 + assert np.isclose(_integrate(mesh, "ds", 1, outlet), 1.0) + + groups = read_groups(filename, comm) + assert groups == {"left": [2], "right": [3], "inlet": [4], "outlet": [5]} + + +def test_prism(tmp_dir): + comm = MPI.COMM_WORLD + filename = tmp_dir / "prism.med" + x, blocks = _med_blocks(_fixture_mesh(dolfinx.mesh.CellType.prism), _left_right) + _on_root(comm, _write_med, filename, x, blocks, _families) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + assert mesh.topology.cell_type == dolfinx.mesh.CellType.prism + assert np.isclose(_integrate(mesh, "dx"), 1.0) + ct = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="right", backend="med") + assert np.isclose(_integrate(mesh, "dx", 1, ct), 0.5) + + +@pytest.mark.parametrize("numbering", ["reversed", "none"]) +def test_node_numbering_ignored(tmp_dir, numbering): + """``NOE/NUM`` is optional user numbering, connectivity refers to node positions.""" + comm = MPI.COMM_WORLD + filename = tmp_dir / "numbering.med" + x, blocks = _med_blocks(_fixture_mesh(dolfinx.mesh.CellType.quadrilateral), _left_right) + node_num = np.arange(x.shape[0], 0, -1) if numbering == "reversed" else None + _on_root(comm, _write_med, filename, x, blocks, _families, node_num=node_num) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + assert np.isclose(_integrate(mesh, "dx"), 1.0) + x_mesh = ufl.SpatialCoordinate(mesh) + first_moment = dolfinx.fem.assemble_scalar(dolfinx.fem.form(x_mesh[0] * ufl.dx)) + assert np.isclose(comm.allreduce(first_moment, op=MPI.SUM), 0.5) + + +def test_overlapping_groups(tmp_dir): + comm = MPI.COMM_WORLD + filename = tmp_dir / "overlap.med" + + def family(midpoints): + return np.select([midpoints[:, 0] < 0.25, midpoints[:, 0] < 0.75], [-2, -3], -4) + + x, blocks = _med_blocks(_fixture_mesh(dolfinx.mesh.CellType.triangle), family) + families = {-2: ["A"], -3: ["A", "B"], -4: ["B"]} + _on_root(comm, _write_med, filename, x, blocks, families) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + for name in ("A", "B"): + tags = io4dolfinx.read_meshtags(filename, mesh, meshtag_name=name, backend="med") + assert np.isclose(_integrate(mesh, "dx", 1, tags), 0.75) + assert read_groups(filename, comm) == {"A": [2, 3], "B": [3, 4]} + + +def test_multiple_meshes(tmp_dir): + comm = MPI.COMM_WORLD + filename = tmp_dir / "two_meshes.med" + square = _med_blocks(_fixture_mesh(dolfinx.mesh.CellType.triangle), _left_right) + cube = _med_blocks(_fixture_mesh(dolfinx.mesh.CellType.hexahedron), _left_right) + + def write(): + _write_med(filename, square[0], square[1], _families, mesh_name="square") + _write_med(filename, cube[0], cube[1], _families, mesh_name="cube", mode="a") + + _on_root(comm, write) + + with pytest.raises(KeyError, match="Available"): + io4dolfinx.read_mesh(filename, comm, backend="med") + mesh = io4dolfinx.read_mesh(filename, comm, backend="med", mesh_name="cube") + assert mesh.topology.cell_type == dolfinx.mesh.CellType.hexahedron + ct = io4dolfinx.read_meshtags( + filename, mesh, meshtag_name="left", backend="med", mesh_name="cube" + ) + assert np.isclose(_integrate(mesh, "dx", 1, ct), 0.5) + + +def test_unsupported(tmp_dir): + comm = MPI.COMM_WORLD + + # Mixed cell types of the highest dimension + mixed = tmp_dir / "mixed.med" + x = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 1, 1]], dtype=float) + blocks = { + "TE4": (np.array([[0, 1, 3, 4]]), np.zeros(1)), + "PY5": (np.array([[0, 1, 2, 3, 5]]), np.zeros(1)), + } + _on_root(comm, _write_med, mixed, x, blocks) + with pytest.raises(NotImplementedError, match="Mixed element types"): + io4dolfinx.read_mesh(mixed, comm, backend="med") + + # Unknown group, and a group of nodes only + groups = tmp_dir / "groups.med" + x, blocks = _med_blocks(_fixture_mesh(dolfinx.mesh.CellType.triangle), _left_right) + _on_root(comm, _write_med, groups, x, blocks, _families, node_families={2: ["corner"]}) + mesh = io4dolfinx.read_mesh(groups, comm, backend="med") + with pytest.raises(KeyError, match="inlet"): + io4dolfinx.read_meshtags(groups, mesh, meshtag_name="wall", backend="med") + with pytest.raises(NotImplementedError, match="only holds nodes"): + io4dolfinx.read_meshtags(groups, mesh, meshtag_name="corner", backend="med") + + # Cell block without entries + empty = tmp_dir / "empty.med" + _on_root(comm, _write_med, empty, x, {"TR3": (np.zeros((0, 3), dtype=np.int64), np.zeros(0))}) + with pytest.raises(ValueError, match="not a positive multiple"): + io4dolfinx.read_mesh(empty, comm, backend="med") + + # Element type without a DOLFINx cell type (polygon) + polygon = tmp_dir / "polygon.med" + blocks["POG"] = (np.array([[0, 1, 2]]), np.zeros(1)) + _on_root(comm, _write_med, polygon, x, blocks) + with pytest.raises(NotImplementedError, match=r"\['POG'\] are not supported"): + io4dolfinx.read_mesh(polygon, comm, backend="med") + + +def _moments(mesh, facet_tags=None): + """Volume, a second moment and the length/area of facets tagged 4 and 5.""" + x = ufl.SpatialCoordinate(mesh) + second = dolfinx.fem.assemble_scalar(dolfinx.fem.form(x[0] * x[1] * ufl.dx)) + values = [_integrate(mesh, "dx"), mesh.comm.allreduce(second, op=MPI.SUM)] + if facet_tags is not None: + values += [_integrate(mesh, "ds", tag, facet_tags) for tag in (4, 5)] + return np.array(values) + + +@pytest.mark.parametrize( + "cell_type, serendipity", + [ + (dolfinx.mesh.CellType.triangle, False), + (dolfinx.mesh.CellType.quadrilateral, False), + pytest.param(dolfinx.mesh.CellType.quadrilateral, True, marks=requires_serendipity), + (dolfinx.mesh.CellType.tetrahedron, False), + (dolfinx.mesh.CellType.prism, False), + (dolfinx.mesh.CellType.hexahedron, False), + pytest.param(dolfinx.mesh.CellType.hexahedron, True, marks=requires_serendipity), + ], +) +def test_higher_order(tmp_dir, cell_type, serendipity): + """Curved quadratic geometry is read back exactly, so a wrong node order shows.""" + comm = MPI.COMM_WORLD + filename = tmp_dir / "curved.med" + ref = _fixture_mesh(cell_type, degree=2) + markers = _facet_markers if cell_type != dolfinx.mesh.CellType.prism else () + x, blocks = _med_blocks(ref, _left_right, markers, serendipity=serendipity, deform=_deform) + _on_root(comm, _write_med, filename, x, blocks, _families) + + # Reference: the quadratic Lagrange mesh, deformed in place. Both layouts represent + # the quadratic deformation exactly. + ref_tags = None + if len(markers) > 0: + tdim = ref.topology.dim + facets = [dolfinx.mesh.locate_entities_boundary(ref, tdim - 1, m) for m, _ in markers[:2]] + values = [np.full(len(f), tag, dtype=np.int32) for f, tag in zip(facets, (4, 5))] + order = np.argsort(np.hstack(facets)) + ref_tags = dolfinx.mesh.meshtags( + ref, tdim - 1, np.hstack(facets)[order], np.hstack(values)[order] + ) + gdim = ref.geometry.dim + ref.geometry.x[:, :gdim] = _deform(ref.geometry.x[:, :gdim]) + expected = _moments(ref, ref_tags) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + assert mesh.topology.cell_type == cell_type + assert fenicsx_compat.cmap(mesh).dim == blocks[next(iter(blocks))][0].shape[1] + ft = None + if len(markers) > 0: + ft = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="facet", backend="med") + assert np.allclose(_moments(mesh, ft), expected) + + +@pytest.mark.parametrize( + "cell_type, degree, serendipity", + [ + (dolfinx.mesh.CellType.triangle, 1, False), + (dolfinx.mesh.CellType.tetrahedron, 2, False), + pytest.param(dolfinx.mesh.CellType.hexahedron, 2, True, marks=requires_serendipity), + ], +) +def test_no_facet_block(tmp_dir, cell_type, degree, serendipity): + """Facet tags of a file without facets are empty.""" + comm = MPI.COMM_WORLD + filename = tmp_dir / "no_facets.med" + ref = _fixture_mesh(cell_type, degree=degree) + x, blocks = _med_blocks(ref, _left_right, serendipity=serendipity) + _on_root(comm, _write_med, filename, x, blocks, _families) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + ft = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="facet", backend="med") + assert ft.dim == mesh.topology.dim - 1 + assert comm.allreduce(len(ft.values), op=MPI.SUM) == 0 + + +def test_structured_mesh(tmp_dir): + """A structured grid (``TYP=1``) has no connectivity and is rejected.""" + comm = MPI.COMM_WORLD + filename = tmp_dir / "structured.med" + x = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) + _on_root(comm, _write_med, filename, x, {}, mesh_type=1) + with pytest.raises(NotImplementedError, match="Structured MED meshes"): + io4dolfinx.read_mesh(filename, comm, backend="med") + with pytest.raises(NotImplementedError, match="Structured MED meshes"): + read_groups(filename, comm) + + +@pytest.mark.parametrize("med_type", ["TR7", "P13", "P15"]) +def test_unsupported_layouts(tmp_dir, med_type): + comm = MPI.COMM_WORLD + filename = tmp_dir / "unsupported.med" + num_nodes = int(med_type[1:]) if med_type[0] == "P" else 7 + x = np.random.default_rng(0).random((num_nodes, 3 if med_type[0] == "P" else 2)) + blocks = {med_type: (np.arange(num_nodes).reshape(1, -1), np.zeros(1))} + _on_root(comm, _write_med, filename, x, blocks) + with pytest.raises(NotImplementedError, match="no Lagrange or nodal serendipity"): + io4dolfinx.read_mesh(filename, comm, backend="med") + + +def test_time_dependent(tmp_dir): + """Later steps store new coordinates and link the connectivity to the first step.""" + comm = MPI.COMM_WORLD + filename = tmp_dir / "time_dependent.med" + x, blocks = _med_blocks( + _fixture_mesh(dolfinx.mesh.CellType.triangle), _left_right, _facet_markers + ) + times = (0.0, 0.5, 1.0) + later_steps = [] + for number, t in enumerate(times[1:], start=1): + x_t = x.copy() + x_t[:, 0] *= 1 + t + later_steps.append((number, t, x_t)) + _on_root(comm, _write_med, filename, x, blocks, _families, later_steps=later_steps) + + for t in times: + mesh = io4dolfinx.read_mesh(filename, comm, backend="med", time=t) + args = {"time": t} + ct = io4dolfinx.read_meshtags(filename, mesh, "cell", backend="med", backend_args=args) + ft = io4dolfinx.read_meshtags(filename, mesh, "facet", backend="med", backend_args=args) + left = io4dolfinx.read_meshtags(filename, mesh, "left", backend="med", backend_args=args) + assert np.isclose(_integrate(mesh, "dx"), 1 + t) + assert np.isclose(_integrate(mesh, "dx", 2, ct), (1 + t) / 2) + assert np.isclose(_integrate(mesh, "dx", 1, left), (1 + t) / 2) + assert np.isclose(_integrate(mesh, "ds", 5, ft), 1.0) + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med", time=None) + assert np.isclose(_integrate(mesh, "dx"), 1.0) + with pytest.raises(KeyError, match=r"Available times: \[0.0, 0.5, 1.0\]"): + io4dolfinx.read_mesh(filename, comm, backend="med", time=0.3) + + +def _medcoupling_mesh(mc, xscale=1.0, it=-1, order=-1, t=0.0, stamp_levels=False): + """Triangulated unit square (3x3 nodes, 8 TRI3) with cell and edge groups. + + x is scaled by ``xscale``. ``it``/``order``/``t`` set the step. MEDCoupling writes + connectivity and families into a later step only if the level meshes carry its time + (``stamp_levels``); otherwise they are links to the first step.""" + xs = mc.DataArrayDouble([0.0, 0.5 * xscale, 1.0 * xscale]) + ys = mc.DataArrayDouble([0.0, 0.5, 1.0]) + cm = mc.MEDCouplingCMesh("mesh") + cm.setCoords(xs, ys) + m0 = cm.buildUnstructured() + m0.simplexize(0) + m0.setName("mesh") + m1 = m0.computeSkin() + m1.setName("mesh") + if stamp_levels: + m0.setTime(t, it, order) + m1.setTime(t, it, order) + + mm = mc.MEDFileUMesh() + mm.setName("mesh") + mm.setMeshAtLevel(0, m0) + mm.setMeshAtLevel(-1, m1) + mm.setTime(it, order, t) + + # Classify on the unscaled geometry + xc = m0.computeCellCenterOfMass()[:, 0] + xc /= xscale + left = xc.findIdsLowerThan(0.5) + left.setName("left") + right = xc.findIdsGreaterOrEqualTo(0.5) + right.setName("right") + mm.setGroupsAtLevel(0, [left, right]) + + xe = m1.computeCellCenterOfMass()[:, 0] + xe /= xscale + inlet = xe.findIdsInRange(-1e-12, 1e-12) + inlet.setName("inlet") + outlet = xe.findIdsInRange(1.0 - 1e-12, 1.0 + 1e-12) + outlet.setName("outlet") + mm.setGroupsAtLevel(-1, [inlet, outlet]) + return mm + + +# MED requires the first step of a mesh to be (iteration, order) = (-1, -1) +_medcoupling_steps = [(-1, -1, 0.0), (1, -1, 0.5), (2, -1, 1.0)] + + +def _write_medcoupling(mc, filename, variant): + if variant == "static": + _medcoupling_mesh(mc).write(str(filename), 2) + elif variant == "coords_only": + mm = _medcoupling_mesh(mc) + mm.write(str(filename), 2) + for it, order, t in _medcoupling_steps[1:]: + coords = mm.getCoords().deepCopy() + coords[:, 0] *= 1.0 + t + step = mc.MEDFileUMesh() + step.setName("mesh") + step.setCoords(coords) + step.setTime(it, order, t) + step.write(str(filename), 0) + else: + for k, (it, order, t) in enumerate(_medcoupling_steps): + mm = _medcoupling_mesh(mc, 1.0 + t, it, order, t, stamp_levels=variant == "stamped") + # Mode 2 creates the file, mode 0 appends the step + mm.write(str(filename), 2 if k == 0 else 0) + + +@pytest.mark.parametrize("variant", ["static", "time_dependent", "stamped", "coords_only"]) +def test_medcoupling_reference(tmp_dir, variant): + """Files written by MEDCoupling, Salome's MED library.""" + mc = pytest.importorskip("medcoupling") + comm = MPI.COMM_WORLD + filename = tmp_dir / f"medcoupling_{variant}.med" + _on_root(comm, _write_medcoupling, mc, filename, variant) + + # MEDCoupling numbers families -2 (left), -3 (right), -4 (boundary edges in no + # group), -5 (inlet) and -6 (outlet) + assert read_groups(filename, comm) == { + "left": [2], + "right": [3], + "inlet": [5], + "outlet": [6], + } + times = [0.0] if variant == "static" else [t for _, _, t in _medcoupling_steps] + for t in times: + mesh = io4dolfinx.read_mesh(filename, comm, backend="med", time=t) + args = {"time": t} + ct = io4dolfinx.read_meshtags(filename, mesh, "cell", backend="med", backend_args=args) + ft = io4dolfinx.read_meshtags(filename, mesh, "facet", backend="med", backend_args=args) + left = io4dolfinx.read_meshtags(filename, mesh, "left", backend="med", backend_args=args) + inlet = io4dolfinx.read_meshtags(filename, mesh, "inlet", backend="med", backend_args=args) + assert np.isclose(_integrate(mesh, "dx"), 1 + t) + assert np.isclose(_integrate(mesh, "dx", 2, ct), (1 + t) / 2) + assert np.isclose(_integrate(mesh, "dx", 1, left), (1 + t) / 2) + assert np.isclose(_integrate(mesh, "ds", 6, ft), 1.0) + assert np.isclose(_integrate(mesh, "ds", 1, inlet), 1.0) + + +# MED type -> (linear MEDCoupling type, conversion type in convertLinearCellsToQuadratic) +_medcoupling_conversions = { + "SE3": ("NORM_SEG2", 0), + "TR6": ("NORM_TRI3", 0), + "QU8": ("NORM_QUAD4", 0), + "QU9": ("NORM_QUAD4", 1), + "T10": ("NORM_TETRA4", 0), + "H20": ("NORM_HEXA8", 0), + "H27": ("NORM_HEXA8", 1), +} + + +def _check_med_node_positions(mc, med_type, x): + """Check with MEDCoupling that MED nodes at ``x`` (MED order) form a valid cell. + + MEDCoupling builds the quadratic nodes from the vertices; MEDCoupling cannot build + ``P18`` from ``PE6``, so its face midside and centre nodes are checked instead.""" + if med_type in _medcoupling_conversions: + linear, conversion = _medcoupling_conversions[med_type] + num_vertices = mc.MEDCouplingUMesh.GetNumberOfNodesOfGeometricType(getattr(mc, linear)) + dim = mc.MEDCouplingUMesh.GetDimensionOfGeometricType(getattr(mc, linear)) + m = mc.MEDCouplingUMesh("m", dim) + m.allocateCells(1) + m.insertNextCell(getattr(mc, linear), list(range(num_vertices))) + m.setCoords(mc.DataArrayDouble(x[:num_vertices].ravel().tolist(), num_vertices, 3)) + m.convertLinearCellsToQuadratic(conversion) + nodes = m.getNodeIdsOfCell(0) + # getValues, as toNumPyArray needs MEDCoupling built with NumPy + coordinates = np.reshape(m.getCoords().getValues(), (-1, 3)) + np.testing.assert_allclose(coordinates[nodes], x, atol=1e-14) + else: + assert med_type == "P18" + m = mc.MEDCouplingUMesh("m", 3) + m.allocateCells(1) + m.insertNextCell(mc.NORM_PENTA18, list(range(len(x)))) + m.setCoords(mc.DataArrayDouble(x.ravel().tolist(), len(x), 3)) + faces = m.buildDescendingConnectivity()[0] + for f in range(faces.getNumberOfCells()): + nodes = faces.getNodeIdsOfCell(f) + nc = 4 if len(nodes) >= 8 else 3 + for k in range(nc): # Midside nodes + midpoint = x[[nodes[k], nodes[(k + 1) % nc]]].mean(axis=0) + np.testing.assert_allclose(x[nodes[nc + k]], midpoint, atol=1e-14) + if len(nodes) == 9: # Face centre + np.testing.assert_allclose(x[nodes[8]], x[nodes[:4]].mean(axis=0), atol=1e-14) + + +@pytest.mark.parametrize( + "med_type", + [ + pytest.param(t, marks=requires_serendipity) if t in ("QU8", "H20") else t + for t in [*_medcoupling_conversions, "P18"] + ], +) +def test_med_permutations_match_medcoupling(med_type): + """The MED to DOLFINx node permutations agree with MEDCoupling's MED cell model.""" + mc = pytest.importorskip("medcoupling") + cell_type, degree, serendipity = med_backend.geometry_layout(med_type) + name = dolfinx.mesh.to_string(cell_type) + if serendipity: + element = dolfinx.io.utils._serendipity_coordinate_element(cell_type, np.float64) + points = element.basix_element.points + else: + points = basix.create_element( + basix.ElementFamily.P, basix.CellType[name], degree, basix.LagrangeVariant.equispaced + ).points + # MED node perm[i] is DOLFINx node i + perm = med_backend._med_to_dolfinx[med_type] + x = np.zeros((len(perm), 3)) + x[perm, : points.shape[1]] = points + _check_med_node_positions(mc, med_type, x) + + +def test_salome_file(tmp_dir): + """Salome exported mesh from the FESTIM workshop: two 100x100 squares.""" + comm = MPI.COMM_WORLD + filename = tmp_dir / "Mesh_1.med" + status = download_file_if_not_exists(_festim_url, filename) + if status == DownloadStatus.no_connection: + pytest.skip("No internet connection") + assert status == DownloadStatus.success + + mesh = io4dolfinx.read_mesh(filename, comm, backend="med") + assert mesh.topology.index_map(2).size_global == 603 + assert mesh.geometry.index_map().size_global == 347 + assert np.isclose(_integrate(mesh, "dx"), 20000.0) + + groups = read_groups(filename, comm) + assert groups == { + "left_volume": [6], + "right_volume": [7], + "left_boundary": [8], + "right_boundary": [9], + } + + ct = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="cell", backend="med") + ft = io4dolfinx.read_meshtags(filename, mesh, meshtag_name="facet", backend="med") + assert np.isclose(_integrate(mesh, "dx", 6, ct), 10000.0) + assert np.isclose(_integrate(mesh, "dx", 7, ct), 10000.0) + assert np.isclose(_integrate(mesh, "ds", 8, ft), 100.0) + assert np.isclose(_integrate(mesh, "ds", 9, ft), 100.0) + num_owned = mesh.topology.index_map(1).size_local + assert comm.allreduce(np.count_nonzero(ft.indices < num_owned), op=MPI.SUM) == 30 + + for name, measure, expected in [ + ("left_volume", "dx", 10000.0), + ("right_volume", "dx", 10000.0), + ("left_boundary", "ds", 100.0), + ("right_boundary", "ds", 100.0), + ]: + tags = io4dolfinx.read_meshtags(filename, mesh, meshtag_name=name, backend="med") + assert np.isclose(_integrate(mesh, measure, 1, tags), expected)