From b0b8ad6aa5c246949e3ba5963d963ae0be927691 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Tue, 27 Oct 2026 20:10:27 +0000 Subject: [PATCH 1/8] Add openCARP backend MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Write and read linear tetrahedral meshes in openCARP's ASCII format: nodes (.pts), elements with region tags (.elem), fiber and sheet directions (.lon), and facet or vertex markers (.vtx, .surf). Everything is written in the input ordering and by rank 0, so the files do not depend on the number of processes. Tetrahedra are oriented with positive volume and .surf triangles point out of the mesh. Coordinates are converted to µm with a required scale factor. io4dolfinx.backends.opencarp.write_opencarp writes the mesh, tags and fibers in one call. To support this, MeshData, MeshTagsData and ArrayData get an on_input_mesh flag, and write_cell_data gets an on_input_mesh option that sends cell values to input cell order. Co-Authored-By: Claude Opus 5.5 --- README.md | 1 + _toc.yml | 1 + docs/backends/opencarp.rst | 74 ++ src/io4dolfinx/backends/__init__.py | 5 + src/io4dolfinx/backends/opencarp/__init__.py | 4 + src/io4dolfinx/backends/opencarp/backend.py | 921 +++++++++++++++++++ src/io4dolfinx/backends/opencarp/wrapper.py | 97 ++ src/io4dolfinx/checkpointing.py | 30 +- src/io4dolfinx/original_checkpoint.py | 1 + src/io4dolfinx/structures.py | 7 + tests/test_opencarp.py | 271 ++++++ 11 files changed, 1411 insertions(+), 1 deletion(-) create mode 100644 docs/backends/opencarp.rst create mode 100644 src/io4dolfinx/backends/opencarp/__init__.py create mode 100644 src/io4dolfinx/backends/opencarp/backend.py create mode 100644 src/io4dolfinx/backends/opencarp/wrapper.py create mode 100644 tests/test_opencarp.py diff --git a/README.md b/README.md index d270878..6a906de 100644 --- a/README.md +++ b/README.md @@ -71,6 +71,7 @@ Most meshing formats supports associating data with the nodes of the mesh (the m - [VTKHDF](./docs/backends/vtkhdf.rst): The new scalable format from VTK, called [VTKHDF](https://docs.vtk.org/en/latest/vtk_file_formats/vtkhdf_file_format/index.html) is supported by the `vtkhdf` backend. - [XDMF](./docs/backends/xdmf.rst) (eXtensible Model Data Format): `.xdmf`. The `xdmf` backend supports the `HDF5` encoding, to ensure performance in parallel. - [Gmsh](./docs/backends/gmsh.rst): `.msh` files (MSH 2.2 and 4.1, ASCII and binary) are read and written with the [Gmsh](https://gmsh.info) Python API by the `gmsh` backend, including physical groups as mesh tags and time-dependent point and cell data. +- [openCARP](./docs/backends/opencarp.rst): Linear tetrahedral meshes with region tags, fibers and markers are written to and read from the ASCII format of the cardiac simulator [openCARP](https://opencarp.org) (`.pts`, `.elem`, `.lon`, `.vtx`, `.surf`) by the `opencarp` backend, or in one call with `io4dolfinx.backends.opencarp.write_opencarp`. - [PyVista](./docs/backends/pyvista.rst) (IO backend is meshio): The [pyvista](https://pyvista.org/) backend uses {py:func}`pyvista.read` to read in meshes, point data and cell data. `pyvista` relies on [meshio](https://github.com/nschloe/meshio) for most reading operations (including the XDMF ascii format). diff --git a/_toc.yml b/_toc.yml index 229ba5c..b8569d4 100644 --- a/_toc.yml +++ b/_toc.yml @@ -40,6 +40,7 @@ parts: - file: "docs/backends/vtkhdf" - file: "docs/backends/med" - file: "docs/backends/gmsh" + - file: "docs/backends/opencarp" - caption: Contributing chapters: diff --git a/docs/backends/opencarp.rst b/docs/backends/opencarp.rst new file mode 100644 index 0000000..b878c5c --- /dev/null +++ b/docs/backends/opencarp.rst @@ -0,0 +1,74 @@ +.. _opencarpbackend: + +openCARP-backend +---------------- + +Writes and reads linear tetrahedral meshes in the ASCII mesh format of the cardiac simulator +`openCARP `_: node coordinates (``.pts``), elements with region tags +(``.elem``), fiber and sheet directions per element (``.lon``), and markers as vertex (``.vtx``) +and surface (``.surf``) files. The backend needs no extra dependencies. + +Nodes and elements are written in the ordering of the input mesh, gathered to and written by +rank 0, so the files are the same whatever the number of processes. openCARP expects coordinates +in µm, so the factor converting the mesh coordinates to µm must always be given. + +:func:`io4dolfinx.backends.opencarp.write_opencarp` writes everything in one call: + +.. code-block:: python + + from io4dolfinx.backends.opencarp import write_opencarp + + write_opencarp( + "heart", + mesh, + scale=1000.0, # mesh in mm + cell_tags=cell_tags, + facet_tags=facet_tags, + fiber=f0, + sheet=s0, + markers={"ENDO": (1, 2), "EPI": (2, 2)}, + ) + +This gives ``heart.pts``, ``heart.elem``, ``heart.lon``, and ``heart_ENDO.vtx``, +``heart_ENDO.surf``, ``heart_EPI.vtx`` and ``heart_EPI.surf``. + +- Elements with negative volume have two nodes swapped, so all elements are positively oriented. +- Cells without a cell tag get region ``0`` (set ``backend_args={"default_tag": ...}`` to change it). +- Fibers and sheets are interpolated into cell midpoints and normalised. Fields in a quadrature + space cannot be written; generate them in a P-1 or DG-0 space. +- Each facet marker value gets a ``.vtx`` file with its nodes and a ``.surf`` file with its facets, + oriented out of the mesh. The files are named from ``markers``, or by the value if it has no name. + Pass ``backend_args={"surf": False}`` to skip the ``.surf`` files. + +The same files can be written with the generic functions, as long as everything is in the input +ordering: + +.. code-block:: python + + args = {"scale": 1000.0} + io4dolfinx.write_mesh_input_order("heart", mesh, backend="opencarp", backend_args=args) + io4dolfinx.write_meshtags("heart", mesh, cell_tags, backend="opencarp", on_input_mesh=True) + io4dolfinx.write_cell_data( + "heart", f0, 0.0, io4dolfinx.FileMode.write, args, "opencarp", on_input_mesh=True + ) + io4dolfinx.write_cell_data( + "heart", s0, 0.0, io4dolfinx.FileMode.append, args, "opencarp", on_input_mesh=True + ) + +Reading takes the path of any of the mesh files. Coordinates are divided by ``scale``: + +.. code-block:: python + + args = {"scale": 1000.0, "markers": {"ENDO": (1, 2), "EPI": (2, 2)}} + mesh = io4dolfinx.read_mesh("heart.elem", MPI.COMM_WORLD, backend="opencarp", backend_args=args) + regions = io4dolfinx.read_meshtags("heart.elem", mesh, "cell", backend="opencarp") + facets = io4dolfinx.read_meshtags("heart.elem", mesh, "facet", backend="opencarp", backend_args=args) + endo = io4dolfinx.read_meshtags("heart.elem", mesh, "ENDO", backend="opencarp", backend_args=args) + f0 = io4dolfinx.read_cell_data("heart.elem", "fiber", mesh, backend="opencarp") + +``"facet"`` reads all ``.surf`` files of the mesh, and a marker name reads its ``.surf`` file, or its +``.vtx`` file as vertex tags if there is no ``.surf``. + +.. automodule:: io4dolfinx.backends.opencarp.backend + :members: + :exclude-members: read_attributes, read_timestamps, write_attributes, read_dofmap, read_dofs, read_cell_perms, write_function, read_legacy_mesh, snapshot_checkpoint, read_hdf5_array, read_point_data diff --git a/src/io4dolfinx/backends/__init__.py b/src/io4dolfinx/backends/__init__.py index e82f1b2..9718495 100644 --- a/src/io4dolfinx/backends/__init__.py +++ b/src/io4dolfinx/backends/__init__.py @@ -512,6 +512,10 @@ def get_backend(backend: str | None = None) -> IOBackend: from .gmsh import backend as GMSHInterface return GMSHInterface + elif backend == "opencarp": + from .opencarp import backend as OpenCARPInterface + + return OpenCARPInterface else: return import_module(backend) @@ -525,6 +529,7 @@ def get_backend(backend: str | None = None) -> IOBackend: "exodus", "med", "gmsh", + "opencarp", ) diff --git a/src/io4dolfinx/backends/opencarp/__init__.py b/src/io4dolfinx/backends/opencarp/__init__.py new file mode 100644 index 0000000..0032034 --- /dev/null +++ b/src/io4dolfinx/backends/opencarp/__init__.py @@ -0,0 +1,4 @@ +from . import backend +from .wrapper import write_opencarp + +__all__ = ["backend", "write_opencarp"] diff --git a/src/io4dolfinx/backends/opencarp/backend.py b/src/io4dolfinx/backends/opencarp/backend.py new file mode 100644 index 0000000..b1fedd3 --- /dev/null +++ b/src/io4dolfinx/backends/opencarp/backend.py @@ -0,0 +1,921 @@ +""" +openCARP interface to io4dolfinx. + +openCARP stores a mesh as a set of ASCII files sharing a base name, with 0-based +node indices and coordinates in µm. See +https://opencarp.org/documentation/user-manual (section "Mesh files"). + +- ``.pts``: number of nodes, then one ``x y z`` line per node. +- ``.elem``: number of elements, then one ``Tt n0 n1 n2 n3 tag`` line per + linear tetrahedron. The tag is the region of the element; in DOLFINx terms the + value of a cell :class:`dolfinx.mesh.MeshTags`. +- ``.lon``: number of directions (1 for fibers, 2 for fibers and sheets), + then one ``fx fy fz [sx sy sz]`` line per element, in the order of ``.elem``. +- ``_.vtx``: number of nodes, the domain (``extra`` for indices into the whole mesh) + and one node index per line, for stimuli and boundary conditions. +- ``_.surf``: number of triangles, then one ``Tr n0 n1 n2`` line + per triangle. + +Nodes and cells are written in the ordering of the input mesh, so the files do +not depend on the number of processes. Each file is written by rank 0. + +Backend arguments: + +- ``"scale"`` (required for meshes): factor converting mesh coordinates to µm, + e.g. ``1000.0`` for a mesh in mm. Coordinates are divided by it when read. +- ``"default_tag"``: region tag of cells without a cell tag (default ``0``). +- ``"markers"``: map from marker name to value, or to ``(value, dim)`` as used by + cardiac-geometriesx. Names the ``.vtx`` and ``.surf`` files of a marker value, + which otherwise use the value itself. +- ``"surf"``: write ``.surf`` files for facet tags (default ``True``). +- ``"vtx_domain"``: domain written to ``.vtx`` files (default ``"extra"``, meaning + indices into the whole mesh, as written; ``"intra"`` means the myocardium submesh). + +SPDX License identifier: MIT + +Copyright: Henrik Finsberg and Simula Research Laboratory +""" + +import glob +from pathlib import Path +from typing import Any, Callable, TypeVar + +from mpi4py import MPI + +import basix +import dolfinx +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 FileMode, ReadMode + +read_mode = ReadMode.parallel + +T = TypeVar("T") + +# Suffixes stripped from a filename to get the base name of the mesh files +_mesh_suffixes = (".pts", ".elem", ".lon") + +# Number of nodes of the entities written for each topological dimension +_nodes_per_entity = {0: 1, 2: 3, 3: 4} + +# Local faces of a tetrahedron, face ``i`` being opposite vertex ``i`` +_tet_faces = np.array([[1, 2, 3], [0, 2, 3], [0, 1, 3], [0, 1, 2]]) + +# Format of floating point values, enough digits for µm coordinates of a heart +_float_format = "%.10g" + +# Names of the directions in ``.lon``, in column order +_directions = ("fiber", "sheet") + + +def _basename(filename: Path | str) -> Path: + """Get the base name shared by the files of a mesh. + + Args: + filename: Base name, or the path of its ``.pts``, ``.elem`` or ``.lon`` file + + Returns: + The base name + """ + path = Path(filename) + return path.with_suffix("") if path.suffix in _mesh_suffixes else path + + +def _mesh_file(filename: Path | str, suffix: str) -> Path: + """Get the path of one of the mesh files, e.g. ``.elem``.""" + return Path(f"{_basename(filename)}{suffix}") + + +def _marker_file(filename: Path | str, name: str, suffix: str) -> Path: + """Get the path of the file of a marker, e.g. ``_.vtx``.""" + return Path(f"{_basename(filename)}_{name}{suffix}") + + +def _get_scale(backend_args: dict[str, Any]) -> float: + """Get the factor converting mesh coordinates to µm. + + It has no default, as a wrong unit gives a mesh that openCARP loads without + complaint but simulates at the wrong size. + """ + if backend_args.get("scale") is None: + raise ValueError( + "The openCARP backend needs backend_args['scale'], the factor converting mesh " + "coordinates to µm (e.g. 1000.0 for a mesh in mm, 1.0 for a mesh in µm)." + ) + scale = float(backend_args["scale"]) + if not scale > 0: + raise ValueError(f"backend_args['scale'] must be positive, got {scale}.") + return scale + + +def _marker_names(markers: dict[str, Any], dim: int) -> dict[int, str]: + """Get the name of each marker value of entities of a given dimension. + + Args: + markers: Map from name to value, or to ``(value, dim)`` + dim: Topological dimension of the entities + + Returns: + Map from value to name. Markers with another dimension are left out. + """ + names = {} + for name, marker in markers.items(): + if isinstance(marker, (int, np.integer)): + value, marker_dim = int(marker), None + else: + value = int(marker[0]) + marker_dim = int(marker[1]) if len(marker) > 1 else None + if marker_dim is None or marker_dim == dim: + names[value] = str(name) + return names + + +def _marker_value(markers: dict[str, Any], name: str) -> int: + """Get the value of a marker read from ``_.vtx`` or ``.surf``. + + Args: + markers: Map from name to value, or to ``(value, dim)`` + name: Name of the marker in the filename + + Returns: + The value in ``markers``, else the name itself if it is an integer (as + written without markers), else ``1``. + """ + if name in markers: + marker = markers[name] + return int(marker) if isinstance(marker, (int, np.integer)) else int(marker[0]) + try: + return int(name) + except ValueError: + return 1 + + +def _on_root(comm: MPI.Comm, fn: Callable[[], T]) -> T | None: + """Run a function on rank 0 and raise any error it raises on all ranks. + + Returns: + The result of ``fn`` on rank 0, ``None`` on the other ranks + """ + result, error = None, None + if comm.rank == 0: + try: + result = fn() + except Exception as e: + error = e + error = comm.bcast(error, root=0) + if error is not None: + raise error + return result + + +def _gather_rows( + comm: MPI.Comm, local: npt.NDArray, local_range: tuple[int, int] | list[int], num_global: int +) -> npt.NDArray | None: + """Gather rows distributed in contiguous chunks to rank 0. + + Returns: + All ``num_global`` rows on rank 0, ``None`` on the other ranks + """ + chunks = comm.gather((int(local_range[0]), np.asarray(local)), root=0) + if comm.rank != 0: + return None + assert chunks is not None + first = chunks[0][1] + rows = np.empty((num_global, *first.shape[1:]), dtype=first.dtype) + for start, chunk in chunks: + rows[start : start + len(chunk)] = chunk + return rows + + +def _gather_concatenate(comm: MPI.Comm, local: npt.NDArray) -> npt.NDArray | None: + """Concatenate an array over all ranks on rank 0, in rank order.""" + chunks = comm.gather(local, root=0) + return None if chunks is None else np.concatenate(chunks) + + +def _read_count(path: Path) -> int: + """Read the number on the first line of a file.""" + with open(path) as f: + return int(f.readline().split()[0]) + + +def _read_rows( + path: Path, local_range: tuple[int, int] | list[int], num_columns: int, header: int = 1 +) -> npt.NDArray[np.str_]: + """Read a range of rows after the header of a file, as strings. + + Args: + path: The file + local_range: Range ``[start, end)`` of rows to read + num_columns: Number of columns, only used to shape an empty range + header: Number of header lines + + Returns: + The rows, of shape ``(end - start, num_columns)`` + """ + start, end = int(local_range[0]), int(local_range[1]) + if end == start: + return np.zeros((0, num_columns), dtype=np.str_) + return np.loadtxt(path, dtype=np.str_, skiprows=header + start, max_rows=end - start, ndmin=2) + + +def _read_elem( + path: Path, local_range: tuple[int, int] | list[int] | None = None +) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.int32]]: + """Read the tetrahedra of an ``.elem`` file. + + Args: + path: The ``.elem`` file + local_range: Range of elements to read, all if ``None`` + + Returns: + Tuple containing: + - Nodes of each element, of shape ``(num_elements, 4)`` + - Region tag of each element, ``0`` if the file has none + """ + if local_range is None: + local_range = (0, _read_count(path)) + rows = _read_rows(path, local_range, 6) + if unsupported := sorted(set(rows[:, 0]) - {"Tt"}): + raise NotImplementedError( + f"Only linear tetrahedra (Tt) are supported, {path} has {unsupported}." + ) + cells = rows[:, 1:5].astype(np.int64) + if rows.shape[1] > 5: + tags = rows[:, 5].astype(np.int32) + else: + tags = np.zeros(len(cells), dtype=np.int32) + return cells, tags + + +def _write_rows(path: Path, header: list[str], rows: npt.NDArray, fmt: str): + """Write header lines followed by one formatted line per row.""" + with open(path, "w") as f: + f.write("\n".join(header) + "\n") + np.savetxt(f, rows, fmt=fmt) + + +def _write_elem(path: Path, cells: npt.NDArray[np.int64], tags: npt.NDArray[np.int32]): + """Write linear tetrahedra with their region tags to an ``.elem`` file.""" + _write_rows(path, [str(len(cells))], np.column_stack([cells, tags]), "Tt %d %d %d %d %d") + + +def _match_rows(reference: npt.NDArray[np.int64], query: npt.NDArray[np.int64]) -> npt.NDArray: + """Find the rows of ``reference`` holding the same nodes as each row of ``query``. + + Returns: + Index into ``reference`` of each query row, ``-1`` if it has no match. If + several rows match, the last one is used. + """ + if len(query) == 0: + return np.zeros(0, dtype=np.int64) + rows = np.vstack([np.sort(reference, axis=1), np.sort(query, axis=1)]) + _, inverse = np.unique(rows, axis=0, return_inverse=True) + inverse = inverse.reshape(-1) + lookup = np.full(inverse.max() + 1, -1, dtype=np.int64) + lookup[inverse[: len(reference)]] = np.arange(len(reference)) + return lookup[inverse[len(reference) :]] + + +def _signed_volumes(x: npt.NDArray[np.floating], cells: npt.NDArray[np.int64]) -> npt.NDArray: + """Get six times the signed volume of each tetrahedron.""" + p = x[cells] + normal = np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]) + return np.einsum("ij,ij->i", normal, p[:, 3] - p[:, 0]) + + +def _oriented_triangles( + x: npt.NDArray[np.floating], cells: npt.NDArray[np.int64], triangles: npt.NDArray[np.int64] +) -> npt.NDArray[np.int64]: + """Orient facets of the mesh and sort them in a canonical order. + + Each triangle gets the normal pointing out of a tetrahedron it is a face of, + starts with its smallest node and the triangles are sorted, so the result does + not depend on how the facets were distributed. + + Args: + x: Node coordinates + cells: Nodes of each tetrahedron + triangles: Nodes of each facet + + Returns: + The oriented triangles + """ + faces = cells[:, _tet_faces].reshape(-1, 3) + match = _match_rows(faces, triangles) + if np.any(match < 0): + raise ValueError(f"{np.count_nonzero(match < 0)} tagged facets are not faces of the mesh.") + opposite = cells.reshape(-1)[match] + + triangles = triangles.copy() + p = x[triangles] + normal = np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]) + inward = np.einsum("ij,ij->i", normal, x[opposite] - p[:, 0]) > 0 + triangles[inward, 1], triangles[inward, 2] = triangles[inward, 2], triangles[inward, 1] + + # Rotating the nodes keeps the orientation + shift = np.argmin(triangles, axis=1)[:, None] + triangles = np.take_along_axis(triangles, (shift + np.arange(3)) % 3, axis=1) + return triangles[np.lexsort(triangles.T[::-1])] + + +def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]: + args = dict(arguments or {}) + args.setdefault("default_tag", 0) + args.setdefault("markers", {}) + args.setdefault("surf", True) + args.setdefault("vtx_domain", "extra") + 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 openCARP 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 openCARP 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 openCARP 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 openCARP 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 openCARP 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 linear tetrahedral mesh to ``.pts`` and ``.elem``. + + Every cell gets the region tag ``backend_args["default_tag"]``, until cell tags + are written with :func:`write_meshtags`. Tetrahedra with negative volume have + two nodes swapped, so that all elements are positively oriented. A ``.lon`` + from an earlier mesh is removed. + + Args: + filename: Base name of the files, or the path of one of them + comm: MPI communicator used in storage + mesh: Mesh data in input ordering, see + :func:`io4dolfinx.write_mesh_input_order` + backend_args: Arguments to backend, ``"scale"`` is required + mode: Only :attr:`FileMode.write` is supported + time: Not used + """ + args = get_default_backend_args(backend_args) + if not mesh.on_input_mesh: + raise ValueError( + "The openCARP backend writes meshes in the input ordering, so that the files do " + "not depend on the number of processes. Use io4dolfinx.write_mesh_input_order or " + "io4dolfinx.backends.opencarp.write_opencarp." + ) + if mesh.cell_type != "tetrahedron" or mesh.degree != 1: + raise NotImplementedError( + "The openCARP backend only supports linear tetrahedral meshes, got a " + f"{mesh.cell_type} mesh of degree {mesh.degree}." + ) + if mesh.local_geometry.shape[1] != 3: + raise ValueError(f"Expected 3D coordinates, got {mesh.local_geometry.shape[1]}D.") + if mode != FileMode.write: + raise NotImplementedError("The openCARP backend can only write a mesh to new files.") + scale = _get_scale(args) + + x = _gather_rows(comm, mesh.local_geometry, mesh.local_geometry_pos, mesh.num_nodes_global) + cells = _gather_rows(comm, mesh.local_topology, mesh.local_topology_pos, mesh.num_cells_global) + + def write(): + assert x is not None and cells is not None + volumes = _signed_volumes(x, cells) + if np.any(volumes == 0): + raise ValueError(f"The mesh has {np.count_nonzero(volumes == 0)} degenerate cells.") + flip = volumes < 0 + cells[flip, 1], cells[flip, 2] = cells[flip, 2], cells[flip, 1] + + _write_rows(_mesh_file(filename, ".pts"), [str(len(x))], x * scale, _float_format) + tags = np.full(len(cells), args["default_tag"], dtype=np.int32) + _write_elem(_mesh_file(filename, ".elem"), cells, tags) + _mesh_file(filename, ".lon").unlink(missing_ok=True) + + _on_root(comm, write) + + +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 a linear tetrahedral mesh from ``.pts`` and ``.elem``. + + Node and cell data is read in contiguous chunks on each process. + + Args: + filename: Base name of the files, or the path of one of them + comm: The MPI communciator to distribute the mesh over + time: Not used + read_from_partition: Read mesh with partition from file (not supported) + backend_args: Arguments to backend, ``"scale"`` is required + Returns: + The mesh topology, geometry, UFL domain and partition function + """ + if read_from_partition: + raise NotImplementedError("The openCARP backend cannot read partitioning data.") + args = get_default_backend_args(backend_args) + scale = _get_scale(args) + pts = _mesh_file(filename, ".pts") + elem = _mesh_file(filename, ".elem") + check_file_exists(pts) + check_file_exists(elem) + + x = _read_rows(pts, compute_local_range(comm, _read_count(pts)), 3).astype(np.float64) + cells, _ = _read_elem(elem, compute_local_range(comm, _read_count(elem))) + return ReadMeshData( + cells=cells, + cell_type="tetrahedron", + x=x / scale, + lvar=int(basix.LagrangeVariant.equispaced), + degree=1, + partition_graph=None, + ) + + +def write_meshtags( + filename: str | Path, + comm: MPI.Comm, + data: MeshTagsData, + backend_args: dict[str, Any] | None = None, +): + """Write mesh tags to the openCARP files of a mesh written earlier. + + - Cell tags set the region tags in ``.elem``. Untagged cells keep their tag. + - Facet tags are written, per value, to ``_.vtx`` with the nodes of + the tagged facets, and ``_.surf`` with the facets as triangles + with outward normals (if ``backend_args["surf"]``). + - Vertex tags are written, per value, to ``_.vtx``. + + ```` is the name of the value in ``backend_args["markers"]``, or the + value itself. + + Args: + filename: Base name of the files, or the path of one of them + comm: MPI communicator used in storage + data: Mesh tags with nodes in input ordering, see the ``on_input_mesh`` + argument of :func:`io4dolfinx.write_meshtags` + backend_args: Arguments to backend + """ + args = get_default_backend_args(backend_args) + if not data.on_input_mesh: + raise ValueError( + "The openCARP backend writes meshtags in the input node ordering. Use " + "io4dolfinx.write_meshtags(..., on_input_mesh=True) or " + "io4dolfinx.backends.opencarp.write_opencarp." + ) + if data.dim not in _nodes_per_entity: + raise NotImplementedError( + f"The openCARP backend can write cell, facet and vertex tags, got dimension {data.dim}." + ) + num_nodes = _nodes_per_entity[data.dim] + indices = _gather_concatenate( + comm, np.asarray(data.indices, dtype=np.int64).reshape(-1, num_nodes) + ) + values = _gather_concatenate(comm, np.asarray(data.values, dtype=np.int32)) + elem = _mesh_file(filename, ".elem") + + def write(): + assert indices is not None and values is not None + check_file_exists(elem) + cells, tags = _read_elem(elem) + if data.dim == 3: + rows = _match_rows(cells, indices) + if np.any(rows < 0): + raise ValueError(f"{np.count_nonzero(rows < 0)} tagged cells are not in {elem}.") + tags[rows] = values + _write_elem(elem, cells, tags) + return + + write_surf = data.dim == 2 and args["surf"] + if write_surf: + pts = _mesh_file(filename, ".pts") + x = _read_rows(pts, (0, _read_count(pts)), 3).astype(np.float64) + names = _marker_names(args["markers"], data.dim) + for value in np.unique(values): + entities = indices[values == value] + name = names.get(int(value), str(value)) + nodes = np.unique(entities) + _write_rows( + _marker_file(filename, name, ".vtx"), + [str(len(nodes)), args["vtx_domain"]], + nodes, + "%d", + ) + if write_surf: + triangles = _oriented_triangles(x, cells, entities) + _write_rows( + _marker_file(filename, name, ".surf"), + [str(len(triangles))], + triangles, + "Tr %d %d %d", + ) + + _on_root(comm, write) + + +def read_meshtags_data( + filename: str | Path, + comm: MPI.Comm, + name: str, + backend_args: dict[str, Any] | None = None, +) -> MeshTagsData: + """Read mesh tags from openCARP files. + + Args: + filename: Base name of the files, or the path of one of them + comm: MPI communicator used in storage + name: ``"cell"`` for the region tags in ``.elem``, ``"facet"`` for the + triangles of all ``_*.surf`` files, otherwise the name of a marker: + the triangles of ``_.surf`` or, if there is none, the nodes of + ``_.vtx``. + backend_args: Arguments to backend. The value of a marker is taken from + ``"markers"``, else from its name if that is an integer, else it is ``1``. + + Returns: + Internal data structure for the mesh tags read from file + """ + args = get_default_backend_args(backend_args) + if name == "cell": + elem = _mesh_file(filename, ".elem") + check_file_exists(elem) + cells, tags = _read_elem(elem, compute_local_range(comm, _read_count(elem))) + return MeshTagsData(name=name, values=tags, indices=cells, dim=3) + + base = _basename(filename) + if name == "facet": + files = sorted(Path(p) for p in glob.glob(f"{glob.escape(str(base))}_*.surf")) + if len(files) == 0: + raise FileNotFoundError(f"Could not find any {base}_*.surf files.") + dim = 2 + elif (surf := _marker_file(filename, name, ".surf")).exists(): + files, dim = [surf], 2 + elif (vtx := _marker_file(filename, name, ".vtx")).exists(): + files, dim = [vtx], 0 + else: + available = sorted( + p.name[len(base.name) + 1 :].rsplit(".", 1)[0] + for p in base.parent.glob(f"{glob.escape(base.name)}_*") + if p.suffix in (".surf", ".vtx") + ) + raise KeyError(f"Could not find marker '{name}' for {base}. Available: {available}") + + # Marker files are small, so rank 0 reads them all + num_nodes = _nodes_per_entity[dim] + indices = np.zeros((0, num_nodes), dtype=np.int64) + values = np.zeros(0, dtype=np.int32) + if comm.rank == 0: + all_indices, all_values = [], [] + for path in files: + marker = path.name[len(base.name) + 1 : -len(path.suffix)] + if dim == 2: + rows = _read_rows(path, (0, _read_count(path)), 4) + entities = rows[:, 1:4].astype(np.int64) + else: + rows = _read_rows(path, (0, _read_count(path)), 1, header=2) + entities = rows.astype(np.int64).reshape(-1, 1) + all_indices.append(entities) + all_values.append(np.full(len(entities), _marker_value(args["markers"], marker))) + indices = np.vstack(all_indices) + values = np.concatenate(all_values).astype(np.int32) + return MeshTagsData(name=name, values=values, indices=indices, 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 openCARP 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 openCARP 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 openCARP 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 openCARP backend cannot write functions. Use io4dolfinx.write_cell_data for fibers." + ) + + +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 openCARP 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 openCARP 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 the fiber or sheet direction of each element from ``.lon``. + + Args: + filename: Base name of the files, or the path of one of them + name: ``"fiber"`` or ``"sheet"`` + comm: Communicator to launch IO on. + time: Not used + 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. + """ + lon = _mesh_file(filename, ".lon") + elem = _mesh_file(filename, ".elem") + check_file_exists(lon) + check_file_exists(elem) + num_directions = _read_count(lon) + if name not in _directions[:num_directions]: + raise KeyError( + f"Could not find '{name}' in {lon}. Available: {list(_directions[:num_directions])}" + ) + local_range = compute_local_range(comm, _read_count(elem)) + cells, _ = _read_elem(elem, local_range) + rows = _read_rows(lon, local_range, 3 * num_directions).astype(np.float64) + column = 3 * _directions.index(name) + return cells, np.ascontiguousarray(rows[:, column : column + 3]) + + +def read_function_names( + filename: Path | str, comm: MPI.Comm, backend_args: dict[str, Any] | None +) -> list[str]: + """Read the names of the directions in ``.lon``. + + Args: + filename: Base name of the files, or the path of one of them + comm: MPI communicator to launch IO on. + backend_args: Arguments to backend + + Returns: + ``["fiber"]`` or ``["fiber", "sheet"]`` + """ + lon = _mesh_file(filename, ".lon") + check_file_exists(lon) + return list(_directions[: _read_count(lon)]) + + +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 direction per element to ``.lon``, normalised to unit length. + + Writing with :attr:`FileMode.write` gives the fiber direction. Writing again + with :attr:`FileMode.append` adds the sheet direction. + + Args: + filename: Base name of the files, or the path of one of them + array_data: Cell data in input ordering, see the ``on_input_mesh`` argument + of :func:`io4dolfinx.write_cell_data` + comm: MPI communicator to open the file with + time: Not used + mode: Write for the fiber, append for the sheet direction + backend_args: The backend arguments + """ + if array_data.type != "Cell": + raise NotImplementedError( + "The openCARP backend can only write cell data (fibers and sheets)." + ) + if not array_data.on_input_mesh: + raise ValueError( + "The openCARP backend writes cell data in the input cell ordering. Use " + "io4dolfinx.write_cell_data(..., on_input_mesh=True) or " + "io4dolfinx.backends.opencarp.write_opencarp." + ) + num_cells, num_components = array_data.global_shape + if num_components != 3: + raise ValueError(f"Expected a 3D vector per cell, got {num_components} components.") + directions = _gather_rows(comm, array_data.values, array_data.local_range, num_cells) + lon = _mesh_file(filename, ".lon") + elem = _mesh_file(filename, ".elem") + + def write(): + assert directions is not None + check_file_exists(elem) + if (num_elements := _read_count(elem)) != num_cells: + raise ValueError(f"Got data for {num_cells} cells, but {elem} has {num_elements}.") + norms = np.linalg.norm(directions, axis=1) + if np.any(norms == 0): + raise ValueError(f"{np.count_nonzero(norms == 0)} cells have a zero direction.") + rows = directions / norms[:, None] + + if mode == FileMode.append: + check_file_exists(lon) + if _read_count(lon) != 1: + raise ValueError(f"{lon} already holds fibers and sheets.") + fibers = _read_rows(lon, (0, num_cells), 3).astype(np.float64) + rows = np.hstack([fibers, rows]) + _write_rows(lon, [str(rows.shape[1] // 3)], rows, _float_format) + + _on_root(comm, write) diff --git a/src/io4dolfinx/backends/opencarp/wrapper.py b/src/io4dolfinx/backends/opencarp/wrapper.py new file mode 100644 index 0000000..00afa78 --- /dev/null +++ b/src/io4dolfinx/backends/opencarp/wrapper.py @@ -0,0 +1,97 @@ +# Copyright (C) 2026 Henrik Finsberg +# +# This file is part of io4dolfinx +# +# SPDX-License-Identifier: MIT + +"""Write a mesh, its markers and fibers to openCARP in one call.""" + +from __future__ import annotations + +import logging +from pathlib import Path +from typing import Any + +import dolfinx + +from ...checkpointing import write_cell_data, write_meshtags +from ...original_checkpoint import write_mesh_input_order +from .. import FileMode + +__all__ = ["write_opencarp"] +logger = logging.getLogger(__name__) + + +def _check_direction(u: dolfinx.fem.Function, name: str): + """Check that a fiber or sheet field can be interpolated into cell midpoints.""" + if u.function_space.ufl_element().family_name == "quadrature": + raise ValueError( + f"Cannot write the {name} field {u.name!r} from a quadrature space to openCARP, as it " + "cannot be interpolated into DG-0. Generate the fibers in a P-1 or DG-0 space." + ) + if u.ufl_shape != (3,): + raise ValueError(f"Expected the {name} field to be a 3D vector, got shape {u.ufl_shape}.") + + +def write_opencarp( + filename: Path | str, + mesh: dolfinx.mesh.Mesh, + scale: float, + cell_tags: dolfinx.mesh.MeshTags | None = None, + facet_tags: dolfinx.mesh.MeshTags | None = None, + fiber: dolfinx.fem.Function | None = None, + sheet: dolfinx.fem.Function | None = None, + markers: dict[str, Any] | None = None, + backend_args: dict[str, Any] | None = None, +): + """Write a linear tetrahedral mesh with markers and fibers to openCARP files. + + Writes ``.pts`` and ``.elem``, with the region tag of each element from + ``cell_tags``; ``.lon`` with the fiber (and sheet) direction of each element; + and ``_.vtx`` and ``_.surf`` for each value of + ``facet_tags``. Nodes and elements are in the ordering of the input mesh, so the + files do not depend on the number of processes. See + :mod:`io4dolfinx.backends.opencarp.backend` for the format. + + Args: + filename: Base name of the files + mesh: The mesh + scale: Factor converting mesh coordinates to µm, e.g. ``1000.0`` for a mesh in mm + cell_tags: Region tag of each cell. Untagged cells get + ``backend_args["default_tag"]`` (``0`` by default). + facet_tags: Facet markers, e.g. for stimuli and boundary conditions + fiber: Fiber direction, interpolated into cell midpoints and normalised + sheet: Sheet direction, written next to the fibers + markers: Map from marker name to value, or to ``(value, dim)``, naming the + files of each facet marker + backend_args: Further arguments to the openCARP backend + """ + if sheet is not None and fiber is None: + raise ValueError("Cannot write a sheet direction without a fiber direction.") + for u, name in ((fiber, "fiber"), (sheet, "sheet")): + if u is not None: + _check_direction(u, name) + + args = dict(backend_args or {}) + args["scale"] = scale + if markers is not None: + args["markers"] = markers + + logger.debug(f"Writing openCARP mesh to {filename} with arguments {args}") + write_mesh_input_order(filename, mesh, backend="opencarp", backend_args=args) + for tags in (cell_tags, facet_tags): + if tags is not None: + write_meshtags( + filename, mesh, tags, backend="opencarp", backend_args=args, on_input_mesh=True + ) + for u, mode in ((fiber, FileMode.write), (sheet, FileMode.append)): + if u is not None: + write_cell_data( + filename, + u, + time=0.0, + mode=mode, + backend_args=args, + backend="opencarp", + on_input_mesh=True, + ) diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index 1822873..df74d97 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -22,6 +22,7 @@ from .backends import FileMode, ReadMode, get_backend from .comm_helpers import ( + exchange_to_owners, send_and_recv_cell_perm, send_dofmap_and_recv_values, send_dofs_and_recv_values, @@ -237,6 +238,7 @@ def write_meshtags( local_start=local_start, dim=meshtags.dim, cell_type=tag_ct, + on_input_mesh=on_input_mesh, ) # Get backend and default arguments @@ -768,6 +770,7 @@ def write_cell_data( backend_args: dict[str, Any] | None, backend: str = "vtkhdf", mesh_name: str | None = None, + on_input_mesh: bool = False, ): """Write function to file by interpolating into cell midpoints. @@ -780,6 +783,8 @@ def write_cell_data( backend_args: The backend arguments mesh_name: Name of the mesh in the file to write the data against. Defaults to the mesh written without an explicit name. + on_input_mesh: If True, the cell values are written in the cell ordering of the + input mesh, matching a mesh written with :func:`write_mesh_input_order`. """ backend_args = _with_mesh_name(backend_args, mesh_name) logger.debug(f"Writing cell data to {filename} for function {u.name} at time {time}") @@ -792,8 +797,31 @@ def write_cell_data( num_dofs_local = V.dofmap.index_map.size_local data = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[:num_dofs_local] + if on_input_mesh: + # Send the value of each owned cell to the process holding its input cell index, + # so that each process gets a contiguous chunk of the input cell ordering + mesh = V.mesh + tdim = mesh.topology.dim + num_owned_cells = mesh.topology.index_map(tdim).size_local + num_cells_global = mesh.topology.index_map(tdim).size_global + input_cells = mesh.topology.original_cell_index[:num_owned_cells].astype(np.int64) + cell_values = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[ + V.dofmap.list[:num_owned_cells, 0] + ] + owners = index_owner(comm, input_cells, num_cells_global) + (recv_cells, recv_values), *_ = exchange_to_owners(comm, owners, [input_cells, cell_values]) + input_range = compute_local_range(comm, num_cells_global) + data = np.empty_like(recv_values) + data[recv_cells - input_range[0]] = recv_values + local_range = (input_range[0], input_range[1]) + ad = ArrayData( - name=v_out.name, values=data, global_shape=data_shape, local_range=local_range, type="Cell" + name=v_out.name, + values=data, + global_shape=data_shape, + local_range=local_range, + type="Cell", + on_input_mesh=on_input_mesh, ) logger.debug( f"Using {backend} backend with arguments {backend_args} and mode {mode} to write cell data" diff --git a/src/io4dolfinx/original_checkpoint.py b/src/io4dolfinx/original_checkpoint.py index 39d1dc8..e9a0263 100644 --- a/src/io4dolfinx/original_checkpoint.py +++ b/src/io4dolfinx/original_checkpoint.py @@ -204,6 +204,7 @@ def create_original_mesh_data(mesh: dolfinx.mesh.Mesh) -> MeshData: ownership_offset=None, partition_range=None, partition_global=None, + on_input_mesh=True, ) diff --git a/src/io4dolfinx/structures.py b/src/io4dolfinx/structures.py index e19e97e..83349a1 100644 --- a/src/io4dolfinx/structures.py +++ b/src/io4dolfinx/structures.py @@ -45,6 +45,9 @@ class MeshData: None #: Local insert position for partitioning information ) partition_global: int | None = None #: Global size of partitioning array + #: Nodes and cells use the ordering of the input mesh rather than the DOLFINx global + #: ordering, which depends on the partitioning + on_input_mesh: bool = False @dataclass @@ -98,6 +101,8 @@ class MeshTagsData: # Optional info to help visualization cell_type: str | None = None #: The cell type + #: Entity nodes use the node ordering of the input mesh + on_input_mesh: bool = False @dataclass @@ -107,3 +112,5 @@ class ArrayData: global_shape: tuple[int, int] # Global shape of data local_range: tuple[int, int] # Start and stop of local array in global array type: str # Type of array-data + #: Rows use the cell ordering of the input mesh (only for cell data) + on_input_mesh: bool = False diff --git a/tests/test_opencarp.py b/tests/test_opencarp.py new file mode 100644 index 0000000..0c84fe1 --- /dev/null +++ b/tests/test_opencarp.py @@ -0,0 +1,271 @@ +from pathlib import Path + +from mpi4py import MPI + +import basix.ufl +import dolfinx +import numpy as np +import pytest +import ufl + +import io4dolfinx +from io4dolfinx.backends.opencarp import write_opencarp + +markers = {"left": (1, 2), "right": (2, 2)} + + +def fiber_expr(x): + return np.vstack([2 + x[1], x[0], np.zeros_like(x[0])]) + + +def sheet_expr(x): + return np.vstack([np.zeros_like(x[0]), np.zeros_like(x[0]), 3 + x[2]]) + + +def normalised(values): + return values / np.linalg.norm(values, axis=0) + + +def create_problem(comm: MPI.Intracomm): + """Unit cube with region 1 (x < 0.4) and 2, facet markers on x=0 and x=1, and + linear fiber and sheet fields (exact in P-1).""" + mesh = dolfinx.mesh.create_unit_cube(comm, 3, 3, 3) + tdim = mesh.topology.dim + mesh.topology.create_connectivity(tdim - 1, tdim) + + cells = np.arange(mesh.topology.index_map(tdim).size_local, dtype=np.int32) + midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, cells) + regions = np.where(midpoints[:, 0] < 0.4, 1, 2).astype(np.int32) + cell_tags = dolfinx.mesh.meshtags(mesh, tdim, cells, regions) + + left = dolfinx.mesh.locate_entities_boundary(mesh, tdim - 1, lambda x: np.isclose(x[0], 0)) + right = dolfinx.mesh.locate_entities_boundary(mesh, tdim - 1, lambda x: np.isclose(x[0], 1)) + facets = np.concatenate([left, right]) + values = np.concatenate([np.full_like(left, 1), np.full_like(right, 2)]).astype(np.int32) + order = np.argsort(facets) + facet_tags = dolfinx.mesh.meshtags(mesh, tdim - 1, facets[order], values[order]) + + V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (3,))) + fiber = dolfinx.fem.Function(V, name="f0") + fiber.interpolate(fiber_expr) + sheet = dolfinx.fem.Function(V, name="s0") + sheet.interpolate(sheet_expr) + return mesh, cell_tags, facet_tags, fiber, sheet + + +def write_problem(filename: Path, comm: MPI.Intracomm, backend_args=None): + mesh, cell_tags, facet_tags, fiber, sheet = create_problem(comm) + write_opencarp( + filename, + mesh, + scale=1000.0, + cell_tags=cell_tags, + facet_tags=facet_tags, + fiber=fiber, + sheet=sheet, + markers=markers, + backend_args=backend_args, + ) + + +def read_lines(path: Path) -> list[list[str]]: + return [line.split() for line in path.read_text().splitlines()] + + +def test_write_content(shared_tmp_path): + filename = shared_tmp_path / "cube" + write_problem(filename, MPI.COMM_WORLD) + if MPI.COMM_WORLD.rank != 0: + return + + pts = read_lines(filename.with_suffix(".pts")) + assert int(pts[0][0]) == len(pts) - 1 == 64 + x = np.array(pts[1:], dtype=np.float64) + # Scaled from the unit cube to µm + assert np.allclose(x.min(axis=0), 0) and np.allclose(x.max(axis=0), 1000) + + elem = read_lines(filename.with_suffix(".elem")) + assert int(elem[0][0]) == len(elem) - 1 == 6 * 27 + assert {row[0] for row in elem[1:]} == {"Tt"} + cells = np.array([row[1:5] for row in elem[1:]], dtype=np.int64) + tags = np.array([row[5] for row in elem[1:]], dtype=np.int32) + assert cells.min() == 0 and cells.max() == 63 + + # All elements positively oriented, filling the cube + p = x[cells] + volumes = np.einsum( + "ij,ij->i", np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]), p[:, 3] - p[:, 0] + ) + assert np.all(volumes > 0) + assert np.isclose(volumes.sum() / 6, 1000.0**3) + + # Region tags follow the midpoints + midpoints = p.mean(axis=1) + np.testing.assert_array_equal(tags, np.where(midpoints[:, 0] < 400, 1, 2)) + + # One unit fiber and sheet per element, in element order + lon = read_lines(filename.with_suffix(".lon")) + assert lon[0] == ["2"] + directions = np.array(lon[1:], dtype=np.float64) + assert directions.shape == (len(cells), 6) + np.testing.assert_allclose(directions[:, :3], normalised(fiber_expr(midpoints.T / 1000)).T) + np.testing.assert_allclose(directions[:, 3:], normalised(sheet_expr(midpoints.T / 1000)).T) + + for name, plane in (("left", 0.0), ("right", 1000.0)): + vtx = read_lines(filename.parent / f"cube_{name}.vtx") + assert vtx[1] == ["extra"] + nodes = np.array(vtx[2:], dtype=np.int64).reshape(-1) + assert int(vtx[0][0]) == len(nodes) == 16 + assert np.allclose(x[nodes, 0], plane) + + surf = read_lines(filename.parent / f"cube_{name}.surf") + assert int(surf[0][0]) == len(surf) - 1 == 2 * 9 + assert {row[0] for row in surf[1:]} == {"Tr"} + triangles = np.array([row[1:] for row in surf[1:]], dtype=np.int64) + np.testing.assert_array_equal(np.unique(triangles), np.sort(nodes)) + # Normals point out of the cube + q = x[triangles] + normals = np.cross(q[:, 1] - q[:, 0], q[:, 2] - q[:, 0]) + assert np.all(np.sign(normals[:, 0]) == (-1 if plane == 0 else 1)) + + +def test_files_independent_of_partitioning(shared_tmp_path): + parallel = shared_tmp_path / "parallel" / "cube" + serial = shared_tmp_path / "serial" / "cube" + if MPI.COMM_WORLD.rank == 0: + parallel.parent.mkdir() + serial.parent.mkdir() + MPI.COMM_WORLD.Barrier() + write_problem(parallel, MPI.COMM_WORLD) + if MPI.COMM_WORLD.rank == 0: + write_problem(serial, MPI.COMM_SELF) + files = sorted(p.name for p in serial.parent.iterdir()) + assert files == sorted(p.name for p in parallel.parent.iterdir()) + assert len(files) == 7 + for name in files: + assert (serial.parent / name).read_bytes() == (parallel.parent / name).read_bytes() + + +def test_read_write(shared_tmp_path): + write_problem(shared_tmp_path / "cube", MPI.COMM_WORLD) + filename = shared_tmp_path / "cube.elem" + comm = MPI.COMM_WORLD + backend_args = {"scale": 1000.0, "markers": markers} + + mesh = io4dolfinx.read_mesh(filename, comm, backend="opencarp", backend_args=backend_args) + tdim = mesh.topology.dim + assert mesh.topology.index_map(tdim).size_global == 6 * 27 + assert mesh.geometry.index_map().size_global == 64 + volume = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1 * ufl.dx(domain=mesh))) + assert np.isclose(comm.allreduce(volume, op=MPI.SUM), 1.0) + + cell_tags = io4dolfinx.read_meshtags(filename, mesh, "cell", backend="opencarp") + midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, cell_tags.indices) + np.testing.assert_array_equal(cell_tags.values, np.where(midpoints[:, 0] < 0.4, 1, 2)) + num_owned = mesh.topology.index_map(tdim).size_local + assert comm.allreduce(np.count_nonzero(cell_tags.indices < num_owned)) == 6 * 27 + + facet_tags = io4dolfinx.read_meshtags( + filename, mesh, "facet", backend="opencarp", backend_args=backend_args + ) + midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim - 1, facet_tags.indices) + np.testing.assert_allclose(midpoints[:, 0], np.where(facet_tags.values == 1, 0.0, 1.0)) + num_owned = mesh.topology.index_map(tdim - 1).size_local + assert comm.allreduce(np.count_nonzero(facet_tags.indices < num_owned)) == 2 * 18 + + right = io4dolfinx.read_meshtags( + filename, mesh, "right", backend="opencarp", backend_args=backend_args + ) + assert right.dim == tdim - 1 + assert np.all(right.values == 2) + assert np.allclose(dolfinx.mesh.compute_midpoints(mesh, tdim - 1, right.indices)[:, 0], 1) + + assert io4dolfinx.read_function_names(filename, comm, backend="opencarp") == ["fiber", "sheet"] + for name, expr in (("fiber", fiber_expr), ("sheet", sheet_expr)): + u = io4dolfinx.read_cell_data(filename, name, mesh, backend="opencarp") + num_cells = mesh.topology.index_map(tdim).size_local + cells = np.arange(num_cells, dtype=np.int32) + midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, cells) + values = u.x.array.reshape(-1, 3)[u.function_space.dofmap.list[cells, 0]] + np.testing.assert_allclose(values, normalised(expr(midpoints.T)).T) + + +def test_read_vertex_markers(shared_tmp_path): + """Without ``.surf`` files a marker is read from its ``.vtx`` file, as vertex tags.""" + write_problem(shared_tmp_path / "cube", MPI.COMM_WORLD, backend_args={"surf": False}) + filename = shared_tmp_path / "cube.elem" + assert not (shared_tmp_path / "cube_left.surf").exists() + + mesh = io4dolfinx.read_mesh( + filename, MPI.COMM_WORLD, backend="opencarp", backend_args={"scale": 1000.0} + ) + left = io4dolfinx.read_meshtags( + filename, mesh, "left", backend="opencarp", backend_args={"markers": markers} + ) + assert left.dim == 0 + assert np.all(left.values == 1) + assert np.allclose(dolfinx.mesh.compute_midpoints(mesh, 0, left.indices)[:, 0], 0) + num_owned = mesh.topology.index_map(0).size_local + assert MPI.COMM_WORLD.allreduce(np.count_nonzero(left.indices < num_owned)) == 16 + + +def test_fiber_only_and_unnamed_markers(shared_tmp_path): + """Markers without names are written and read back by value.""" + filename = shared_tmp_path / "cube" + mesh, _, facet_tags, fiber, _ = create_problem(MPI.COMM_WORLD) + write_opencarp(filename, mesh, scale=1.0, facet_tags=facet_tags, fiber=fiber) + if MPI.COMM_WORLD.rank == 0: + assert read_lines(filename.with_suffix(".lon"))[0] == ["1"] + tags = {row[5] for row in read_lines(filename.with_suffix(".elem"))[1:]} + assert tags == {"0"} + + filename = filename.with_suffix(".elem") + assert io4dolfinx.read_function_names(filename, MPI.COMM_WORLD, backend="opencarp") == ["fiber"] + in_mesh = io4dolfinx.read_mesh( + filename, MPI.COMM_WORLD, backend="opencarp", backend_args={"scale": 1.0} + ) + right = io4dolfinx.read_meshtags(filename, in_mesh, "2", backend="opencarp") + assert np.all(right.values == 2) + with pytest.raises(KeyError): + io4dolfinx.read_cell_data(filename, "sheet", in_mesh, backend="opencarp") + with pytest.raises(KeyError): + io4dolfinx.read_meshtags(filename, in_mesh, "missing", backend="opencarp") + + +def test_errors(shared_tmp_path): + comm = MPI.COMM_WORLD + filename = shared_tmp_path / "cube" + mesh, cell_tags, _, fiber, sheet = create_problem(comm) + + with pytest.raises(ValueError, match="scale"): + io4dolfinx.write_mesh_input_order(filename, mesh, backend="opencarp") + with pytest.raises(ValueError, match="input ordering"): + io4dolfinx.write_mesh(filename, mesh, backend="opencarp", backend_args={"scale": 1.0}) + with pytest.raises(ValueError, match="without a fiber"): + write_opencarp(filename, mesh, scale=1.0, sheet=sheet) + + Q = dolfinx.fem.functionspace( + mesh, basix.ufl.quadrature_element(mesh.basix_cell(), value_shape=(3,), degree=2) + ) + with pytest.raises(ValueError, match="quadrature"): + write_opencarp(filename, mesh, scale=1.0, fiber=dolfinx.fem.Function(Q)) + + hex_mesh = dolfinx.mesh.create_unit_cube( + comm, 2, 2, 2, cell_type=dolfinx.mesh.CellType.hexahedron + ) + with pytest.raises(NotImplementedError, match="tetrahedral"): + write_opencarp(filename, hex_mesh, scale=1.0) + + # Errors found when writing on rank 0 are raised on all ranks + write_opencarp(filename, mesh, scale=1.0) + with pytest.raises(ValueError, match="input node ordering"): + io4dolfinx.write_meshtags(filename, mesh, cell_tags, backend="opencarp") + zero = dolfinx.fem.Function(fiber.function_space) + with pytest.raises(ValueError, match="zero direction"): + io4dolfinx.write_cell_data( + filename, zero, 0.0, io4dolfinx.FileMode.write, None, "opencarp", on_input_mesh=True + ) + with pytest.raises(ValueError, match="input cell ordering"): + io4dolfinx.write_cell_data( + filename, fiber, 0.0, io4dolfinx.FileMode.write, None, "opencarp" + ) From f473e09450b8c08bef828efbc46418568a6ed888 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 28 Oct 2026 06:34:13 +0000 Subject: [PATCH 2/8] Fix mypy --- src/io4dolfinx/checkpointing.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index df74d97..9300f0d 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -809,6 +809,7 @@ def write_cell_data( V.dofmap.list[:num_owned_cells, 0] ] owners = index_owner(comm, input_cells, num_cells_global) + assert isinstance(comm, MPI.Intracomm) (recv_cells, recv_values), *_ = exchange_to_owners(comm, owners, [input_cells, cell_values]) input_range = compute_local_range(comm, num_cells_global) data = np.empty_like(recv_values) From 786766e654756bb89f1a562034225e3607936d3f Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 28 Oct 2026 08:36:23 +0000 Subject: [PATCH 3/8] Add on_input_mesh option to write_point_data Like write_cell_data, write_point_data can now send node values in the node ordering of the input mesh, matching write_mesh_input_order. Both use a shared helper that moves values to the process holding their input index. Co-Authored-By: Claude Opus 5.5 --- src/io4dolfinx/checkpointing.py | 57 ++++++++++++++++++++++----- tests/test_write_data_input_order.py | 58 ++++++++++++++++++++++++++++ 2 files changed, 105 insertions(+), 10 deletions(-) create mode 100644 tests/test_write_data_input_order.py diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index 9300f0d..7f8ed63 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -81,6 +81,33 @@ def _with_mesh_name( return args +def _to_input_order( + comm: MPI.Intracomm, + input_indices: npt.NDArray[np.int64], + values: npt.NDArray, + num_global: int, +) -> tuple[npt.NDArray, tuple[int, int]]: + """Send values to the process holding their input index, so that each process + gets a contiguous chunk of the input ordering. + + Args: + comm: The MPI communicator + input_indices: Input index of each owned entity on this process + values: Value of each owned entity, of shape ``(len(input_indices), ...)`` + num_global: Global number of entities + + Returns: + The values of the entities in this process' chunk, in input order, and the range + of the chunk + """ + owners = index_owner(comm, input_indices, num_global) + (recv_indices, recv_values), *_ = exchange_to_owners(comm, owners, [input_indices, values]) + local_range = compute_local_range(comm, num_global) + data = np.empty_like(recv_values) + data[recv_indices - local_range[0]] = recv_values + return data, (local_range[0], local_range[1]) + + def write_attributes( filename: Path | str, comm: MPI.Comm, @@ -726,6 +753,7 @@ def write_point_data( backend_args: dict[str, Any] | None, backend: str = "vtkhdf", mesh_name: str | None = None, + on_input_mesh: bool = False, ): """Write function to file by interpolating into geometry nodes. @@ -739,6 +767,8 @@ def write_point_data( backend: Which backend to use. mesh_name: Name of the mesh in the file to write the data against. Defaults to the mesh written without an explicit name. + on_input_mesh: If True, the node values are written in the node ordering of the + input mesh, matching a mesh written with :func:`write_mesh_input_order`. """ backend_args = _with_mesh_name(backend_args, mesh_name) logger.debug(f"Writing point data to {filename} for function {u.name} at time {time}") @@ -750,8 +780,22 @@ def write_point_data( local_range = V.dofmap.index_map.local_range num_dofs_local = V.dofmap.index_map.size_local data = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[:num_dofs_local] + + if on_input_mesh: + # The dofs of the geometry space are the geometry nodes + input_nodes = np.asarray( + V.mesh.geometry.input_global_indices[:num_dofs_local], dtype=np.int64 + ) + assert isinstance(comm, MPI.Intracomm) + data, local_range = _to_input_order(comm, input_nodes, data, data_shape[0]) + ad = ArrayData( - name=v_out.name, values=data, global_shape=data_shape, local_range=local_range, type="Point" + name=v_out.name, + values=data, + global_shape=data_shape, + local_range=local_range, + type="Point", + on_input_mesh=on_input_mesh, ) logger.debug( f"Using {backend} backend with arguments {backend_args} and mode {mode} to write point data" @@ -798,23 +842,16 @@ def write_cell_data( data = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[:num_dofs_local] if on_input_mesh: - # Send the value of each owned cell to the process holding its input cell index, - # so that each process gets a contiguous chunk of the input cell ordering + # DG-0 dofs of owned cells are owned, so send the value of each owned cell mesh = V.mesh tdim = mesh.topology.dim num_owned_cells = mesh.topology.index_map(tdim).size_local - num_cells_global = mesh.topology.index_map(tdim).size_global input_cells = mesh.topology.original_cell_index[:num_owned_cells].astype(np.int64) cell_values = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[ V.dofmap.list[:num_owned_cells, 0] ] - owners = index_owner(comm, input_cells, num_cells_global) assert isinstance(comm, MPI.Intracomm) - (recv_cells, recv_values), *_ = exchange_to_owners(comm, owners, [input_cells, cell_values]) - input_range = compute_local_range(comm, num_cells_global) - data = np.empty_like(recv_values) - data[recv_cells - input_range[0]] = recv_values - local_range = (input_range[0], input_range[1]) + data, local_range = _to_input_order(comm, input_cells, cell_values, data_shape[0]) ad = ArrayData( name=v_out.name, diff --git a/tests/test_write_data_input_order.py b/tests/test_write_data_input_order.py new file mode 100644 index 0000000..beecde3 --- /dev/null +++ b/tests/test_write_data_input_order.py @@ -0,0 +1,58 @@ +from types import SimpleNamespace + +from mpi4py import MPI + +import dolfinx +import numpy as np +import pytest + +import io4dolfinx +from io4dolfinx import checkpointing +from io4dolfinx.original_checkpoint import create_original_mesh_data + + +def f(x): + return np.vstack([x[0] + 2 * x[1], 3 * x[2] - x[0], np.ones_like(x[0])]) + + +@pytest.fixture +def captured(monkeypatch): + """Capture the data passed to the backend, keyed by type.""" + data = {} + + def write_data(filename, array_data, **kwargs): + data[array_data.type] = array_data + + monkeypatch.setattr( + checkpointing, "get_backend", lambda backend=None: SimpleNamespace(write_data=write_data) + ) + return data + + +@pytest.mark.parametrize( + "cell_type", [dolfinx.mesh.CellType.tetrahedron, dolfinx.mesh.CellType.hexahedron] +) +def test_data_on_input_mesh(captured, cell_type): + """Point and cell data written with ``on_input_mesh`` follow the nodes and cells of + :func:`io4dolfinx.write_mesh_input_order`.""" + mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 4, 3, 5, cell_type=cell_type) + V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (3,))) + u = dolfinx.fem.Function(V) + u.interpolate(f) + + for writer in (io4dolfinx.write_point_data, io4dolfinx.write_cell_data): + writer("unused", u, 0.0, io4dolfinx.FileMode.write, None, "any", on_input_mesh=True) + mesh_data = create_original_mesh_data(mesh) + + point_data = captured["Point"] + assert point_data.on_input_mesh + assert tuple(point_data.local_range) == tuple(mesh_data.local_geometry_pos) + np.testing.assert_allclose(point_data.values, f(mesh_data.local_geometry.T).T, atol=1e-12) + + # P-1 interpolated into DG-0 is the mean of the vertex values for these cells + cell_data = captured["Cell"] + assert cell_data.on_input_mesh + assert tuple(cell_data.local_range) == tuple(mesh_data.local_topology_pos) + x = np.vstack(MPI.COMM_WORLD.allgather(mesh_data.local_geometry)) + midpoints = x[mesh_data.local_topology].mean(axis=1) + np.testing.assert_allclose(cell_data.values, f(midpoints.T).T, atol=1e-12) From d073dccf2a2a802fbb118e3f312b132e902be2a0 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 28 Oct 2026 10:03:44 +0000 Subject: [PATCH 4/8] Support triangle, quadrilateral and hexahedral meshes in openCARP backend Write and read Tr and Qd elements for 2D meshes and Hx elements for 3D meshes, besides Tt. Node orderings and orientations follow openCARP's mesher: counter-clockwise in 2D, positive volume in 3D, and hexahedra listing one face followed by the nodes opposite nodes 0, 3, 2 and 1. 2D meshes are written with z = 0 and read back with 2D coordinates and 2D fibers. Their markers are written as .vtx files only. Facet markers of hexahedral meshes are written as outward Qd faces in .surf files. Errors raised while reading chunks on some ranks are now raised on all ranks, so that no rank waits in a later collective call. Co-Authored-By: Claude Opus 5.5 --- README.md | 2 +- docs/backends/opencarp.rst | 20 +- src/io4dolfinx/backends/opencarp/backend.py | 390 ++++++++++++++------ src/io4dolfinx/backends/opencarp/wrapper.py | 13 +- tests/test_opencarp.py | 260 ++++++++----- 5 files changed, 466 insertions(+), 219 deletions(-) diff --git a/README.md b/README.md index 6a906de..b439e34 100644 --- a/README.md +++ b/README.md @@ -71,7 +71,7 @@ Most meshing formats supports associating data with the nodes of the mesh (the m - [VTKHDF](./docs/backends/vtkhdf.rst): The new scalable format from VTK, called [VTKHDF](https://docs.vtk.org/en/latest/vtk_file_formats/vtkhdf_file_format/index.html) is supported by the `vtkhdf` backend. - [XDMF](./docs/backends/xdmf.rst) (eXtensible Model Data Format): `.xdmf`. The `xdmf` backend supports the `HDF5` encoding, to ensure performance in parallel. - [Gmsh](./docs/backends/gmsh.rst): `.msh` files (MSH 2.2 and 4.1, ASCII and binary) are read and written with the [Gmsh](https://gmsh.info) Python API by the `gmsh` backend, including physical groups as mesh tags and time-dependent point and cell data. -- [openCARP](./docs/backends/opencarp.rst): Linear tetrahedral meshes with region tags, fibers and markers are written to and read from the ASCII format of the cardiac simulator [openCARP](https://opencarp.org) (`.pts`, `.elem`, `.lon`, `.vtx`, `.surf`) by the `opencarp` backend, or in one call with `io4dolfinx.backends.opencarp.write_opencarp`. +- [openCARP](./docs/backends/opencarp.rst): Linear triangle, quadrilateral, tetrahedral and hexahedral meshes with region tags, fibers and markers are written to and read from the ASCII format of the cardiac simulator [openCARP](https://opencarp.org) (`.pts`, `.elem`, `.lon`, `.vtx`, `.surf`) by the `opencarp` backend, or in one call with `io4dolfinx.backends.opencarp.write_opencarp`. - [PyVista](./docs/backends/pyvista.rst) (IO backend is meshio): The [pyvista](https://pyvista.org/) backend uses {py:func}`pyvista.read` to read in meshes, point data and cell data. `pyvista` relies on [meshio](https://github.com/nschloe/meshio) for most reading operations (including the XDMF ascii format). diff --git a/docs/backends/opencarp.rst b/docs/backends/opencarp.rst index b878c5c..16bc1a8 100644 --- a/docs/backends/opencarp.rst +++ b/docs/backends/opencarp.rst @@ -3,7 +3,7 @@ openCARP-backend ---------------- -Writes and reads linear tetrahedral meshes in the ASCII mesh format of the cardiac simulator +Writes and reads linear meshes in the ASCII mesh format of the cardiac simulator `openCARP `_: node coordinates (``.pts``), elements with region tags (``.elem``), fiber and sheet directions per element (``.lon``), and markers as vertex (``.vtx``) and surface (``.surf``) files. The backend needs no extra dependencies. @@ -32,12 +32,18 @@ in µm, so the factor converting the mesh coordinates to µm must always be give This gives ``heart.pts``, ``heart.elem``, ``heart.lon``, and ``heart_ENDO.vtx``, ``heart_ENDO.surf``, ``heart_EPI.vtx`` and ``heart_EPI.surf``. -- Elements with negative volume have two nodes swapped, so all elements are positively oriented. +- Triangle (``Tr``) and quadrilateral (``Qd``) meshes in 2D, and tetrahedral (``Tt``) and hexahedral + (``Hx``) meshes in 3D are supported. 2D meshes are written with ``z = 0``, and read back with 2D + coordinates. +- Elements get the orientation used by openCARP's ``mesher``: positive volume in 3D, + counter-clockwise in 2D. - Cells without a cell tag get region ``0`` (set ``backend_args={"default_tag": ...}`` to change it). -- Fibers and sheets are interpolated into cell midpoints and normalised. Fields in a quadrature - space cannot be written; generate them in a P-1 or DG-0 space. -- Each facet marker value gets a ``.vtx`` file with its nodes and a ``.surf`` file with its facets, - oriented out of the mesh. The files are named from ``markers``, or by the value if it has no name. +- Fibers and sheets are interpolated into cell midpoints and normalised. 2D directions get + ``z = 0``, and are read back with two components. Fields in a quadrature space cannot be + written; generate them in a P-1 or DG-0 space. +- Each facet marker value gets a ``.vtx`` file with its nodes and, for 3D meshes, a ``.surf`` file + with its facets (``Tr`` or ``Qd``), oriented out of the mesh. The files are named from + ``markers``, or by the value if it has no name. Pass ``backend_args={"surf": False}`` to skip the ``.surf`` files. The same files can be written with the generic functions, as long as everything is in the input @@ -67,7 +73,7 @@ Reading takes the path of any of the mesh files. Coordinates are divided by ``sc f0 = io4dolfinx.read_cell_data("heart.elem", "fiber", mesh, backend="opencarp") ``"facet"`` reads all ``.surf`` files of the mesh, and a marker name reads its ``.surf`` file, or its -``.vtx`` file as vertex tags if there is no ``.surf``. +``.vtx`` file as vertex tags if there is no ``.surf`` (as for 2D meshes). .. automodule:: io4dolfinx.backends.opencarp.backend :members: diff --git a/src/io4dolfinx/backends/opencarp/backend.py b/src/io4dolfinx/backends/opencarp/backend.py index b1fedd3..5142176 100644 --- a/src/io4dolfinx/backends/opencarp/backend.py +++ b/src/io4dolfinx/backends/opencarp/backend.py @@ -5,20 +5,27 @@ node indices and coordinates in µm. See https://opencarp.org/documentation/user-manual (section "Mesh files"). -- ``.pts``: number of nodes, then one ``x y z`` line per node. -- ``.elem``: number of elements, then one ``Tt n0 n1 n2 n3 tag`` line per - linear tetrahedron. The tag is the region of the element; in DOLFINx terms the +- ``.pts``: number of nodes, then one ``x y z`` line per node. 2D meshes have + ``z = 0``. +- ``.elem``: number of elements, then one line per element with its type, its + nodes and its tag, e.g. ``Tt n0 n1 n2 n3 tag``. Linear triangles (``Tr``), + quadrilaterals (``Qd``), tetrahedra (``Tt``) and hexahedra (``Hx``) are supported, + one type per mesh. The tag is the region of the element; in DOLFINx terms the value of a cell :class:`dolfinx.mesh.MeshTags`. - ``.lon``: number of directions (1 for fibers, 2 for fibers and sheets), then one ``fx fy fz [sx sy sz]`` line per element, in the order of ``.elem``. - ``_.vtx``: number of nodes, the domain (``extra`` for indices into the whole mesh) and one node index per line, for stimuli and boundary conditions. -- ``_.surf``: number of triangles, then one ``Tr n0 n1 n2`` line - per triangle. +- ``_.surf``: number of faces, then one ``Tr n0 n1 n2`` or + ``Qd n0 n1 n2 n3`` line per face of a 3D mesh. Nodes and cells are written in the ordering of the input mesh, so the files do not depend on the number of processes. Each file is written by rank 0. +Elements are written with the orientation used by openCARP's ``mesher``: positive +volume in 3D, counter-clockwise in the xy-plane in 2D. Hexahedra list one face +(nodes 0-3) followed by the nodes opposite nodes 0, 3, 2 and 1. + Backend arguments: - ``"scale"`` (required for meshes): factor converting mesh coordinates to µm, @@ -27,7 +34,7 @@ - ``"markers"``: map from marker name to value, or to ``(value, dim)`` as used by cardiac-geometriesx. Names the ``.vtx`` and ``.surf`` files of a marker value, which otherwise use the value itself. -- ``"surf"``: write ``.surf`` files for facet tags (default ``True``). +- ``"surf"``: write ``.surf`` files for facet tags of 3D meshes (default ``True``). - ``"vtx_domain"``: domain written to ``.vtx`` files (default ``"extra"``, meaning indices into the whole mesh, as written; ``"intra"`` means the myocardium submesh). @@ -58,11 +65,43 @@ # Suffixes stripped from a filename to get the base name of the mesh files _mesh_suffixes = (".pts", ".elem", ".lon") -# Number of nodes of the entities written for each topological dimension -_nodes_per_entity = {0: 1, 2: 3, 3: 4} - -# Local faces of a tetrahedron, face ``i`` being opposite vertex ``i`` -_tet_faces = np.array([[1, 2, 3], [0, 2, 3], [0, 1, 3], [0, 1, 2]]) +# openCARP element type of each supported DOLFINx cell type, and the reverse +_element_types = { + "triangle": "Tr", + "quadrilateral": "Qd", + "tetrahedron": "Tt", + "hexahedron": "Hx", +} +_cell_types = {element: cell for cell, element in _element_types.items()} + +# Topological dimension of each openCARP element type +_tdim = {"Tr": 2, "Qd": 2, "Tt": 3, "Hx": 3} + +# Node permutation ``p`` with ``opencarp[i] = dolfinx[p[i]]`` for each element type +_dolfinx_to_opencarp = { + "Tr": np.array([0, 1, 2]), + "Qd": np.array([0, 1, 3, 2]), + "Tt": np.array([0, 1, 2, 3]), + "Hx": np.array([0, 1, 3, 2, 4, 6, 7, 5]), +} + +# Node permutation mirroring an element in openCARP ordering, flipping its orientation +_mirror = { + "Tr": np.array([0, 2, 1]), + "Qd": np.array([0, 3, 2, 1]), + "Tt": np.array([0, 2, 1, 3]), + "Hx": np.array([0, 3, 2, 1, 4, 7, 6, 5]), +} + +# Local faces of 3D elements in openCARP ordering, with their nodes in cyclic order, +# and the openCARP element type of the faces +_faces = { + "Tt": np.array([[1, 2, 3], [0, 2, 3], [0, 1, 3], [0, 1, 2]]), + "Hx": np.array( + [[0, 1, 2, 3], [4, 5, 6, 7], [0, 1, 7, 4], [1, 2, 6, 7], [2, 3, 5, 6], [3, 0, 4, 5]] + ), +} +_face_types = {"Tt": "Tr", "Hx": "Qd"} # Format of floating point values, enough digits for µm coordinates of a heart _float_format = "%.10g" @@ -171,6 +210,20 @@ def _on_root(comm: MPI.Comm, fn: Callable[[], T]) -> T | None: return result +def _on_all(comm: MPI.Comm, fn: Callable[[], T]) -> T: + """Run a function on all ranks and raise an error raised on any rank on all ranks, + so that no rank is left waiting in a later collective call.""" + result, error = None, None + try: + result = fn() + except Exception as e: + error = e + errors = [e for e in comm.allgather(error) if e is not None] + if errors: + raise errors[0] + return result # type: ignore[return-value] + + def _gather_rows( comm: MPI.Comm, local: npt.NDArray, local_range: tuple[int, int] | list[int], num_global: int ) -> npt.NDArray | None: @@ -222,10 +275,24 @@ def _read_rows( return np.loadtxt(path, dtype=np.str_, skiprows=header + start, max_rows=end - start, ndmin=2) +def _elem_type(path: Path) -> str: + """Read the openCARP element type of the first element of an ``.elem`` file.""" + with open(path) as f: + f.readline() + line = f.readline().split() + if len(line) == 0: + raise ValueError(f"{path} has no elements.") + if line[0] not in _tdim: + raise NotImplementedError( + f"openCARP element type {line[0]} in {path} is not supported, only {list(_tdim)}." + ) + return line[0] + + def _read_elem( path: Path, local_range: tuple[int, int] | list[int] | None = None -) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.int32]]: - """Read the tetrahedra of an ``.elem`` file. +) -> tuple[str, npt.NDArray[np.int64], npt.NDArray[np.int32]]: + """Read the elements of an ``.elem`` file. Args: path: The ``.elem`` file @@ -233,22 +300,26 @@ def _read_elem( Returns: Tuple containing: - - Nodes of each element, of shape ``(num_elements, 4)`` + - The openCARP element type + - Nodes of each element in openCARP ordering - Region tag of each element, ``0`` if the file has none """ + elem_type = _elem_type(path) + num_nodes = len(_dolfinx_to_opencarp[elem_type]) if local_range is None: local_range = (0, _read_count(path)) - rows = _read_rows(path, local_range, 6) - if unsupported := sorted(set(rows[:, 0]) - {"Tt"}): + rows = _read_rows(path, local_range, num_nodes + 2) + if other := sorted(set(rows[:, 0]) - {elem_type}): raise NotImplementedError( - f"Only linear tetrahedra (Tt) are supported, {path} has {unsupported}." + f"Meshes with several element types are not supported, {path} has " + f"{[elem_type, *other]}." ) - cells = rows[:, 1:5].astype(np.int64) - if rows.shape[1] > 5: - tags = rows[:, 5].astype(np.int32) + cells = rows[:, 1 : 1 + num_nodes].astype(np.int64) + if rows.shape[1] > 1 + num_nodes: + tags = rows[:, 1 + num_nodes].astype(np.int32) else: tags = np.zeros(len(cells), dtype=np.int32) - return cells, tags + return elem_type, cells, tags def _write_rows(path: Path, header: list[str], rows: npt.NDArray, fmt: str): @@ -258,9 +329,17 @@ def _write_rows(path: Path, header: list[str], rows: npt.NDArray, fmt: str): np.savetxt(f, rows, fmt=fmt) -def _write_elem(path: Path, cells: npt.NDArray[np.int64], tags: npt.NDArray[np.int32]): - """Write linear tetrahedra with their region tags to an ``.elem`` file.""" - _write_rows(path, [str(len(cells))], np.column_stack([cells, tags]), "Tt %d %d %d %d %d") +def _element_format(elem_type: str, num_columns: int) -> str: + """Get the format of an element line with ``num_columns`` integers.""" + return " ".join([elem_type] + ["%d"] * num_columns) + + +def _write_elem( + path: Path, elem_type: str, cells: npt.NDArray[np.int64], tags: npt.NDArray[np.int32] +): + """Write elements with their region tags to an ``.elem`` file.""" + rows = np.column_stack([cells, tags]) + _write_rows(path, [str(len(cells))], rows, _element_format(elem_type, rows.shape[1])) def _match_rows(reference: npt.NDArray[np.int64], query: npt.NDArray[np.int64]) -> npt.NDArray: @@ -280,46 +359,73 @@ def _match_rows(reference: npt.NDArray[np.int64], query: npt.NDArray[np.int64]) return lookup[inverse[len(reference) :]] -def _signed_volumes(x: npt.NDArray[np.floating], cells: npt.NDArray[np.int64]) -> npt.NDArray: - """Get six times the signed volume of each tetrahedron.""" - p = x[cells] - normal = np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]) - return np.einsum("ij,ij->i", normal, p[:, 3] - p[:, 0]) - +def _orientation( + elem_type: str, x: npt.NDArray[np.floating], cells: npt.NDArray[np.int64] +) -> npt.NDArray: + """Get the orientation of each element in openCARP ordering. -def _oriented_triangles( - x: npt.NDArray[np.floating], cells: npt.NDArray[np.int64], triangles: npt.NDArray[np.int64] + Returns: + The triple product of the edges at node 0 for 3D elements, the z-component of + the cross product of the edges at node 0 for 2D elements. Positive for the + orientation of openCARP's ``mesher``. + """ + p = x[cells] + # Nodes 1 and 2 of simplices, 1 and 3 of quadrilaterals and hexahedra, are the + # neighbours of node 0 in the first face + second = 2 if elem_type in ("Tr", "Tt") else 3 + normal = np.cross(p[:, 1] - p[:, 0], p[:, second] - p[:, 0]) + if _tdim[elem_type] == 2: + return normal[:, 2] + third = 3 if elem_type == "Tt" else 4 + return np.einsum("ij,ij->i", normal, p[:, third] - p[:, 0]) + + +def _facet_normals(x: npt.NDArray[np.floating], facets: npt.NDArray[np.int64]) -> npt.NDArray: + """Get a normal of each triangle or quadrilateral, with nodes in cyclic order.""" + p = x[facets] + if facets.shape[1] == 3: + return np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]) + return np.cross(p[:, 2] - p[:, 0], p[:, 3] - p[:, 1]) + + +def _oriented_facets( + x: npt.NDArray[np.floating], + elem_type: str, + cells: npt.NDArray[np.int64], + facets: npt.NDArray[np.int64], ) -> npt.NDArray[np.int64]: - """Orient facets of the mesh and sort them in a canonical order. + """Orient facets of a 3D mesh and sort them in a canonical order. - Each triangle gets the normal pointing out of a tetrahedron it is a face of, - starts with its smallest node and the triangles are sorted, so the result does - not depend on how the facets were distributed. + Each facet gets the normal pointing out of an element it is a face of, starts with + its smallest node and the facets are sorted, so the result does not depend on how + the facets were distributed. Args: x: Node coordinates - cells: Nodes of each tetrahedron - triangles: Nodes of each facet + elem_type: openCARP element type of the mesh + cells: Nodes of each element in openCARP ordering + facets: Nodes of each facet in cyclic order Returns: - The oriented triangles + The oriented facets """ - faces = cells[:, _tet_faces].reshape(-1, 3) - match = _match_rows(faces, triangles) + local_faces = _faces[elem_type] + faces = cells[:, local_faces].reshape(-1, local_faces.shape[1]) + match = _match_rows(faces, facets) if np.any(match < 0): raise ValueError(f"{np.count_nonzero(match < 0)} tagged facets are not faces of the mesh.") - opposite = cells.reshape(-1)[match] + cell_centres = x[cells[match // len(local_faces)]].mean(axis=1) - triangles = triangles.copy() - p = x[triangles] - normal = np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]) - inward = np.einsum("ij,ij->i", normal, x[opposite] - p[:, 0]) > 0 - triangles[inward, 1], triangles[inward, 2] = triangles[inward, 2], triangles[inward, 1] + facets = facets.copy() + outward = x[facets].mean(axis=1) - cell_centres + inward = np.einsum("ij,ij->i", _facet_normals(x, facets), outward) < 0 + facets[inward] = facets[inward][:, _mirror[_face_types[elem_type]]] # Rotating the nodes keeps the orientation - shift = np.argmin(triangles, axis=1)[:, None] - triangles = np.take_along_axis(triangles, (shift + np.arange(3)) % 3, axis=1) - return triangles[np.lexsort(triangles.T[::-1])] + num_nodes = facets.shape[1] + shift = np.argmin(facets, axis=1)[:, None] + facets = np.take_along_axis(facets, (shift + np.arange(num_nodes)) % num_nodes, axis=1) + return facets[np.lexsort(facets.T[::-1])] def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]: @@ -437,12 +543,13 @@ def write_mesh( mode: FileMode = FileMode.write, time: float = 0.0, ): - """Write a linear tetrahedral mesh to ``.pts`` and ``.elem``. + """Write a linear mesh to ``.pts`` and ``.elem``. - Every cell gets the region tag ``backend_args["default_tag"]``, until cell tags - are written with :func:`write_meshtags`. Tetrahedra with negative volume have - two nodes swapped, so that all elements are positively oriented. A ``.lon`` - from an earlier mesh is removed. + Supports triangle and quadrilateral meshes in 2D, and tetrahedral and hexahedral + meshes in 3D. Every cell gets the region tag ``backend_args["default_tag"]``, until + cell tags are written with :func:`write_meshtags`. Elements with the opposite + orientation of openCARP's ``mesher`` are mirrored. A ``.lon`` from an + earlier mesh is removed. Args: filename: Base name of the files, or the path of one of them @@ -460,13 +567,18 @@ def write_mesh( "not depend on the number of processes. Use io4dolfinx.write_mesh_input_order or " "io4dolfinx.backends.opencarp.write_opencarp." ) - if mesh.cell_type != "tetrahedron" or mesh.degree != 1: + if mesh.cell_type not in _element_types or mesh.degree != 1: + raise NotImplementedError( + "The openCARP backend supports linear triangle, quadrilateral, tetrahedral and " + f"hexahedral meshes, got a {mesh.cell_type} mesh of degree {mesh.degree}." + ) + elem_type = _element_types[mesh.cell_type] + gdim = mesh.local_geometry.shape[1] + if gdim != _tdim[elem_type]: raise NotImplementedError( - "The openCARP backend only supports linear tetrahedral meshes, got a " - f"{mesh.cell_type} mesh of degree {mesh.degree}." + f"The openCARP backend needs {_tdim[elem_type]}D coordinates for a " + f"{mesh.cell_type} mesh, got {gdim}D." ) - if mesh.local_geometry.shape[1] != 3: - raise ValueError(f"Expected 3D coordinates, got {mesh.local_geometry.shape[1]}D.") if mode != FileMode.write: raise NotImplementedError("The openCARP backend can only write a mesh to new files.") scale = _get_scale(args) @@ -476,15 +588,18 @@ def write_mesh( def write(): assert x is not None and cells is not None - volumes = _signed_volumes(x, cells) - if np.any(volumes == 0): - raise ValueError(f"The mesh has {np.count_nonzero(volumes == 0)} degenerate cells.") - flip = volumes < 0 - cells[flip, 1], cells[flip, 2] = cells[flip, 2], cells[flip, 1] - - _write_rows(_mesh_file(filename, ".pts"), [str(len(x))], x * scale, _float_format) - tags = np.full(len(cells), args["default_tag"], dtype=np.int32) - _write_elem(_mesh_file(filename, ".elem"), cells, tags) + points = np.zeros((len(x), 3), dtype=x.dtype) + points[:, :gdim] = x + elements = cells[:, _dolfinx_to_opencarp[elem_type]] + orientation = _orientation(elem_type, points, elements) + if np.any(orientation == 0): + raise ValueError(f"The mesh has {np.count_nonzero(orientation == 0)} degenerate cells.") + flip = orientation < 0 + elements[flip] = elements[flip][:, _mirror[elem_type]] + + _write_rows(_mesh_file(filename, ".pts"), [str(len(points))], points * scale, _float_format) + tags = np.full(len(elements), args["default_tag"], dtype=np.int32) + _write_elem(_mesh_file(filename, ".elem"), elem_type, elements, tags) _mesh_file(filename, ".lon").unlink(missing_ok=True) _on_root(comm, write) @@ -497,9 +612,10 @@ def read_mesh_data( read_from_partition: bool = False, backend_args: dict[str, Any] | None = None, ) -> ReadMeshData: - """Read a linear tetrahedral mesh from ``.pts`` and ``.elem``. + """Read a linear mesh from ``.pts`` and ``.elem``. - Node and cell data is read in contiguous chunks on each process. + Node and cell data is read in contiguous chunks on each process. A mesh of 2D + elements must lie in a plane of constant ``z``, and gets 2D coordinates. Args: filename: Base name of the files, or the path of one of them @@ -518,12 +634,26 @@ def read_mesh_data( elem = _mesh_file(filename, ".elem") check_file_exists(pts) check_file_exists(elem) - - x = _read_rows(pts, compute_local_range(comm, _read_count(pts)), 3).astype(np.float64) - cells, _ = _read_elem(elem, compute_local_range(comm, _read_count(elem))) + elem_type = _elem_type(elem) + + def read(): + x = _read_rows(pts, compute_local_range(comm, _read_count(pts)), 3).astype(np.float64) + _, cells, _ = _read_elem(elem, compute_local_range(comm, _read_count(elem))) + return x, cells + + x, cells = _on_all(comm, read) + if _tdim[elem_type] == 2: + z_min = comm.allreduce(x[:, 2].min(initial=np.inf), op=MPI.MIN) + z_max = comm.allreduce(x[:, 2].max(initial=-np.inf), op=MPI.MAX) + if z_max > z_min: + raise NotImplementedError( + f"{elem} has 2D elements that do not lie in a plane of constant z. Surface " + "meshes in 3D are not supported." + ) + x = x[:, :2] return ReadMeshData( - cells=cells, - cell_type="tetrahedron", + cells=cells[:, np.argsort(_dolfinx_to_opencarp[elem_type])], + cell_type=_cell_types[elem_type], x=x / scale, lvar=int(basix.LagrangeVariant.equispaced), degree=1, @@ -540,10 +670,10 @@ def write_meshtags( """Write mesh tags to the openCARP files of a mesh written earlier. - Cell tags set the region tags in ``.elem``. Untagged cells keep their tag. - - Facet tags are written, per value, to ``_.vtx`` with the nodes of - the tagged facets, and ``_.surf`` with the facets as triangles - with outward normals (if ``backend_args["surf"]``). - - Vertex tags are written, per value, to ``_.vtx``. + - Other tags are written, per value, to ``_.vtx`` with the nodes of + the tagged entities. + - Facet tags of a 3D mesh are also written to ``_.surf``, as + triangles or quadrilaterals with outward normals (if ``backend_args["surf"]``). ```` is the name of the value in ``backend_args["markers"]``, or the value itself. @@ -562,38 +692,36 @@ def write_meshtags( "io4dolfinx.write_meshtags(..., on_input_mesh=True) or " "io4dolfinx.backends.opencarp.write_opencarp." ) - if data.dim not in _nodes_per_entity: - raise NotImplementedError( - f"The openCARP backend can write cell, facet and vertex tags, got dimension {data.dim}." - ) - num_nodes = _nodes_per_entity[data.dim] - indices = _gather_concatenate( - comm, np.asarray(data.indices, dtype=np.int64).reshape(-1, num_nodes) - ) + local_indices = np.asarray(data.indices, dtype=np.int64) + num_nodes = data.num_dofs_per_entity or local_indices.shape[-1] + indices = _gather_concatenate(comm, local_indices.reshape(-1, num_nodes)) values = _gather_concatenate(comm, np.asarray(data.values, dtype=np.int32)) elem = _mesh_file(filename, ".elem") def write(): assert indices is not None and values is not None check_file_exists(elem) - cells, tags = _read_elem(elem) - if data.dim == 3: + elem_type, cells, tags = _read_elem(elem) + if data.dim == _tdim[elem_type]: rows = _match_rows(cells, indices) if np.any(rows < 0): raise ValueError(f"{np.count_nonzero(rows < 0)} tagged cells are not in {elem}.") tags[rows] = values - _write_elem(elem, cells, tags) + _write_elem(elem, elem_type, cells, tags) return - write_surf = data.dim == 2 and args["surf"] + write_surf = args["surf"] and _tdim[elem_type] == 3 and data.dim == 2 if write_surf: pts = _mesh_file(filename, ".pts") x = _read_rows(pts, (0, _read_count(pts)), 3).astype(np.float64) + face_type = _face_types[elem_type] + # DOLFINx orders the nodes of a quadrilateral facet as a tensor product + facets = indices[:, _dolfinx_to_opencarp[face_type]] names = _marker_names(args["markers"], data.dim) for value in np.unique(values): - entities = indices[values == value] + marked = values == value name = names.get(int(value), str(value)) - nodes = np.unique(entities) + nodes = np.unique(indices[marked]) _write_rows( _marker_file(filename, name, ".vtx"), [str(len(nodes)), args["vtx_domain"]], @@ -601,12 +729,12 @@ def write(): "%d", ) if write_surf: - triangles = _oriented_triangles(x, cells, entities) + faces = _oriented_facets(x, elem_type, cells, facets[marked]) _write_rows( _marker_file(filename, name, ".surf"), - [str(len(triangles))], - triangles, - "Tr %d %d %d", + [str(len(faces))], + faces, + _element_format(face_type, faces.shape[1]), ) _on_root(comm, write) @@ -624,8 +752,8 @@ def read_meshtags_data( filename: Base name of the files, or the path of one of them comm: MPI communicator used in storage name: ``"cell"`` for the region tags in ``.elem``, ``"facet"`` for the - triangles of all ``_*.surf`` files, otherwise the name of a marker: - the triangles of ``_.surf`` or, if there is none, the nodes of + faces of all ``_*.surf`` files, otherwise the name of a marker: + the faces of ``_.surf`` or, if there is none, the nodes of ``_.vtx``. backend_args: Arguments to backend. The value of a marker is taken from ``"markers"``, else from its name if that is an integer, else it is ``1``. @@ -637,8 +765,11 @@ def read_meshtags_data( if name == "cell": elem = _mesh_file(filename, ".elem") check_file_exists(elem) - cells, tags = _read_elem(elem, compute_local_range(comm, _read_count(elem))) - return MeshTagsData(name=name, values=tags, indices=cells, dim=3) + elem_type, cells, tags = _on_all( + comm, lambda: _read_elem(elem, compute_local_range(comm, _read_count(elem))) + ) + cells = cells[:, np.argsort(_dolfinx_to_opencarp[elem_type])] + return MeshTagsData(name=name, values=tags, indices=cells, dim=_tdim[elem_type]) base = _basename(filename) if name == "facet": @@ -658,24 +789,31 @@ def read_meshtags_data( ) raise KeyError(f"Could not find marker '{name}' for {base}. Available: {available}") - # Marker files are small, so rank 0 reads them all - num_nodes = _nodes_per_entity[dim] - indices = np.zeros((0, num_nodes), dtype=np.int64) - values = np.zeros(0, dtype=np.int32) - if comm.rank == 0: + # Marker files are small, so rank 0 reads them all. read_meshtags reshapes the + # empty indices of the other ranks. + def read(): all_indices, all_values = [], [] for path in files: marker = path.name[len(base.name) + 1 : -len(path.suffix)] if dim == 2: rows = _read_rows(path, (0, _read_count(path)), 4) - entities = rows[:, 1:4].astype(np.int64) + face_types = set(rows[:, 0]) + if len(face_types) > 1 or not face_types <= {"Tr", "Qd"}: + raise NotImplementedError(f"Unsupported face types {face_types} in {path}.") + face_type = face_types.pop() if face_types else "Tr" + perm = np.argsort(_dolfinx_to_opencarp[face_type]) + entities = rows[:, 1:].astype(np.int64)[:, perm] else: rows = _read_rows(path, (0, _read_count(path)), 1, header=2) entities = rows.astype(np.int64).reshape(-1, 1) all_indices.append(entities) all_values.append(np.full(len(entities), _marker_value(args["markers"], marker))) - indices = np.vstack(all_indices) - values = np.concatenate(all_values).astype(np.int32) + return np.vstack(all_indices), np.concatenate(all_values).astype(np.int32) + + indices, values = np.zeros((0, 1), dtype=np.int64), np.zeros(0, dtype=np.int32) + result = _on_root(comm, read) + if result is not None: + indices, values = result return MeshTagsData(name=name, values=values, indices=indices, dim=dim) @@ -825,7 +963,7 @@ def read_cell_data( Returns: A tuple (topology, dofs) where topology contains the vertex indices of the cells, dofs the degrees of - freedom within that cell. + freedom within that cell. A 2D mesh gets the x- and y-components. """ lon = _mesh_file(filename, ".lon") elem = _mesh_file(filename, ".elem") @@ -836,11 +974,18 @@ def read_cell_data( raise KeyError( f"Could not find '{name}' in {lon}. Available: {list(_directions[:num_directions])}" ) - local_range = compute_local_range(comm, _read_count(elem)) - cells, _ = _read_elem(elem, local_range) - rows = _read_rows(lon, local_range, 3 * num_directions).astype(np.float64) + + def read(): + local_range = compute_local_range(comm, _read_count(elem)) + elem_type, cells, _ = _read_elem(elem, local_range) + rows = _read_rows(lon, local_range, 3 * num_directions).astype(np.float64) + return elem_type, cells, rows + + elem_type, cells, rows = _on_all(comm, read) column = 3 * _directions.index(name) - return cells, np.ascontiguousarray(rows[:, column : column + 3]) + dofs = rows[:, column : column + _tdim[elem_type]] + cells = cells[:, np.argsort(_dolfinx_to_opencarp[elem_type])] + return cells, np.ascontiguousarray(dofs) def read_function_names( @@ -871,6 +1016,8 @@ def write_data( ): """Write a direction per element to ``.lon``, normalised to unit length. + 2D directions get ``z = 0``. + Writing with :attr:`FileMode.write` gives the fiber direction. Writing again with :attr:`FileMode.append` adds the sheet direction. @@ -894,8 +1041,8 @@ def write_data( "io4dolfinx.backends.opencarp.write_opencarp." ) num_cells, num_components = array_data.global_shape - if num_components != 3: - raise ValueError(f"Expected a 3D vector per cell, got {num_components} components.") + if num_components not in (2, 3): + raise ValueError(f"Expected a 2D or 3D vector per cell, got {num_components} components.") directions = _gather_rows(comm, array_data.values, array_data.local_range, num_cells) lon = _mesh_file(filename, ".lon") elem = _mesh_file(filename, ".elem") @@ -905,10 +1052,15 @@ def write(): check_file_exists(elem) if (num_elements := _read_count(elem)) != num_cells: raise ValueError(f"Got data for {num_cells} cells, but {elem} has {num_elements}.") - norms = np.linalg.norm(directions, axis=1) + if num_components < (tdim := _tdim[_elem_type(elem)]): + raise ValueError(f"Expected a {tdim}D vector per cell, got {num_components}D.") + # openCARP directions are 3D, also for 2D meshes + vectors = np.zeros((num_cells, 3), dtype=np.float64) + vectors[:, :num_components] = directions + norms = np.linalg.norm(vectors, axis=1) if np.any(norms == 0): raise ValueError(f"{np.count_nonzero(norms == 0)} cells have a zero direction.") - rows = directions / norms[:, None] + rows = vectors / norms[:, None] if mode == FileMode.append: check_file_exists(lon) diff --git a/src/io4dolfinx/backends/opencarp/wrapper.py b/src/io4dolfinx/backends/opencarp/wrapper.py index 00afa78..8f41672 100644 --- a/src/io4dolfinx/backends/opencarp/wrapper.py +++ b/src/io4dolfinx/backends/opencarp/wrapper.py @@ -29,8 +29,11 @@ def _check_direction(u: dolfinx.fem.Function, name: str): f"Cannot write the {name} field {u.name!r} from a quadrature space to openCARP, as it " "cannot be interpolated into DG-0. Generate the fibers in a P-1 or DG-0 space." ) - if u.ufl_shape != (3,): - raise ValueError(f"Expected the {name} field to be a 3D vector, got shape {u.ufl_shape}.") + gdim = u.function_space.mesh.geometry.dim + if u.ufl_shape not in ((gdim,), (3,)): + raise ValueError( + f"Expected the {name} field to be a {gdim}D or 3D vector, got shape {u.ufl_shape}." + ) def write_opencarp( @@ -44,13 +47,13 @@ def write_opencarp( markers: dict[str, Any] | None = None, backend_args: dict[str, Any] | None = None, ): - """Write a linear tetrahedral mesh with markers and fibers to openCARP files. + """Write a linear mesh with markers and fibers to openCARP files. Writes ``.pts`` and ``.elem``, with the region tag of each element from ``cell_tags``; ``.lon`` with the fiber (and sheet) direction of each element; and ``_.vtx`` and ``_.surf`` for each value of - ``facet_tags``. Nodes and elements are in the ordering of the input mesh, so the - files do not depend on the number of processes. See + ``facet_tags`` (``.surf`` only for 3D meshes). Nodes and elements are in the + ordering of the input mesh, so the files do not depend on the number of processes. See :mod:`io4dolfinx.backends.opencarp.backend` for the format. Args: diff --git a/tests/test_opencarp.py b/tests/test_opencarp.py index 0c84fe1..5090fce 100644 --- a/tests/test_opencarp.py +++ b/tests/test_opencarp.py @@ -11,7 +11,21 @@ import io4dolfinx from io4dolfinx.backends.opencarp import write_opencarp -markers = {"left": (1, 2), "right": (2, 2)} +CellType = dolfinx.mesh.CellType + +# openCARP element type, and number of cells, nodes, marked facets and marked nodes +# per marker of the meshes of create_problem +meshes = { + CellType.tetrahedron: ("Tt", 6 * 27, 64, 2 * 9, 16), + CellType.hexahedron: ("Hx", 27, 64, 9, 16), + CellType.triangle: ("Tr", 2 * 9, 16, 3, 4), + CellType.quadrilateral: ("Qd", 9, 16, 3, 4), +} +cell_types = pytest.mark.parametrize("cell_type", list(meshes), ids=lambda c: c.name) + + +def markers(tdim: int): + return {"left": (1, tdim - 1), "right": (2, tdim - 1)} def fiber_expr(x): @@ -19,18 +33,26 @@ def fiber_expr(x): def sheet_expr(x): - return np.vstack([np.zeros_like(x[0]), np.zeros_like(x[0]), 3 + x[2]]) + return np.vstack([np.ones_like(x[0]), 3 + x[2], x[0]]) -def normalised(values): - return values / np.linalg.norm(values, axis=0) +def expected_direction(expr, points, gdim): + """Normalised direction at points, with the components a mesh of dimension gdim + stores and a zero z-component in 2D, as written to ``.lon``.""" + values = np.zeros((3, len(points))) + values[:gdim] = expr(points.T)[:gdim] + return (values / np.linalg.norm(values, axis=0)).T -def create_problem(comm: MPI.Intracomm): - """Unit cube with region 1 (x < 0.4) and 2, facet markers on x=0 and x=1, and - linear fiber and sheet fields (exact in P-1).""" - mesh = dolfinx.mesh.create_unit_cube(comm, 3, 3, 3) +def create_problem(comm: MPI.Intracomm, cell_type: dolfinx.mesh.CellType = CellType.tetrahedron): + """Unit cube or square with region 1 (x < 0.4) and 2, facet markers on x=0 and x=1, + and linear fiber and sheet fields (exact in P-1).""" + if dolfinx.mesh.cell_dim(cell_type) == 3: + mesh = dolfinx.mesh.create_unit_cube(comm, 3, 3, 3, cell_type=cell_type) + else: + mesh = dolfinx.mesh.create_unit_square(comm, 3, 3, cell_type=cell_type) tdim = mesh.topology.dim + gdim = mesh.geometry.dim mesh.topology.create_connectivity(tdim - 1, tdim) cells = np.arange(mesh.topology.index_map(tdim).size_local, dtype=np.int32) @@ -45,16 +67,21 @@ def create_problem(comm: MPI.Intracomm): order = np.argsort(facets) facet_tags = dolfinx.mesh.meshtags(mesh, tdim - 1, facets[order], values[order]) - V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (3,))) + V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (gdim,))) fiber = dolfinx.fem.Function(V, name="f0") - fiber.interpolate(fiber_expr) + fiber.interpolate(lambda x: fiber_expr(x)[:gdim]) sheet = dolfinx.fem.Function(V, name="s0") - sheet.interpolate(sheet_expr) + sheet.interpolate(lambda x: sheet_expr(x)[:gdim]) return mesh, cell_tags, facet_tags, fiber, sheet -def write_problem(filename: Path, comm: MPI.Intracomm, backend_args=None): - mesh, cell_tags, facet_tags, fiber, sheet = create_problem(comm) +def write_problem( + filename: Path, + comm: MPI.Intracomm, + cell_type: dolfinx.mesh.CellType = CellType.tetrahedron, + **kwargs, +): + mesh, cell_tags, facet_tags, fiber, sheet = create_problem(comm, cell_type) write_opencarp( filename, mesh, @@ -63,8 +90,8 @@ def write_problem(filename: Path, comm: MPI.Intracomm, backend_args=None): facet_tags=facet_tags, fiber=fiber, sheet=sheet, - markers=markers, - backend_args=backend_args, + markers=markers(mesh.topology.dim), + **kwargs, ) @@ -72,32 +99,54 @@ def read_lines(path: Path) -> list[list[str]]: return [line.split() for line in path.read_text().splitlines()] -def test_write_content(shared_tmp_path): - filename = shared_tmp_path / "cube" - write_problem(filename, MPI.COMM_WORLD) +def orientation(elem_type: str, p: np.ndarray) -> np.ndarray: + """Signed measure of each element with nodes p in openCARP ordering, positive for + the orientation of openCARP's mesher, and the factor to its volume or area.""" + e1 = p[:, 1] - p[:, 0] + e2 = p[:, 2 if elem_type in ("Tr", "Tt") else 3] - p[:, 0] + if elem_type in ("Tr", "Qd"): + return np.cross(e1, e2)[:, 2], 2 if elem_type == "Tr" else 1 + e3 = p[:, 3 if elem_type == "Tt" else 4] - p[:, 0] + return np.einsum("ij,ij->i", np.cross(e1, e2), e3), 6 if elem_type == "Tt" else 1 + + +@cell_types +def test_write_content(shared_tmp_path, cell_type): + filename = shared_tmp_path / "mesh" + write_problem(filename, MPI.COMM_WORLD, cell_type) if MPI.COMM_WORLD.rank != 0: return + elem_type, num_cells, num_nodes, num_facets, num_marked_nodes = meshes[cell_type] + tdim = dolfinx.mesh.cell_dim(cell_type) pts = read_lines(filename.with_suffix(".pts")) - assert int(pts[0][0]) == len(pts) - 1 == 64 + assert int(pts[0][0]) == len(pts) - 1 == num_nodes x = np.array(pts[1:], dtype=np.float64) - # Scaled from the unit cube to µm - assert np.allclose(x.min(axis=0), 0) and np.allclose(x.max(axis=0), 1000) + assert x.shape == (num_nodes, 3) + # Scaled from the unit cube or square to µm, z = 0 in 2D + assert np.allclose(x.min(axis=0)[:tdim], 0) and np.allclose(x.max(axis=0)[:tdim], 1000) + assert np.all(x[:, tdim:] == 0) elem = read_lines(filename.with_suffix(".elem")) - assert int(elem[0][0]) == len(elem) - 1 == 6 * 27 - assert {row[0] for row in elem[1:]} == {"Tt"} - cells = np.array([row[1:5] for row in elem[1:]], dtype=np.int64) - tags = np.array([row[5] for row in elem[1:]], dtype=np.int32) - assert cells.min() == 0 and cells.max() == 63 + assert int(elem[0][0]) == len(elem) - 1 == num_cells + assert {row[0] for row in elem[1:]} == {elem_type} + cells = np.array([row[1:-1] for row in elem[1:]], dtype=np.int64) + tags = np.array([row[-1] for row in elem[1:]], dtype=np.int32) + assert cells.min() == 0 and cells.max() == num_nodes - 1 - # All elements positively oriented, filling the cube + # All elements oriented as by openCARP's mesher, filling the cube or square p = x[cells] - volumes = np.einsum( - "ij,ij->i", np.cross(p[:, 1] - p[:, 0], p[:, 2] - p[:, 0]), p[:, 3] - p[:, 0] - ) - assert np.all(volumes > 0) - assert np.isclose(volumes.sum() / 6, 1000.0**3) + measure, factor = orientation(elem_type, p) + assert np.all(measure > 0) + assert np.isclose(measure.sum() / factor, 1000.0**tdim) + if elem_type == "Qd": + np.testing.assert_allclose(p[:, 2] - p[:, 0], (p[:, 1] - p[:, 0]) + (p[:, 3] - p[:, 0])) + if elem_type == "Hx": + e1, e2, e3 = p[:, 1] - p[:, 0], p[:, 3] - p[:, 0], p[:, 4] - p[:, 0] + np.testing.assert_allclose(p[:, 2] - p[:, 0], e1 + e2) + np.testing.assert_allclose(p[:, 5] - p[:, 0], e2 + e3) + np.testing.assert_allclose(p[:, 6] - p[:, 0], e1 + e2 + e3) + np.testing.assert_allclose(p[:, 7] - p[:, 0], e1 + e3) # Region tags follow the midpoints midpoints = p.mean(axis=1) @@ -107,55 +156,72 @@ def test_write_content(shared_tmp_path): lon = read_lines(filename.with_suffix(".lon")) assert lon[0] == ["2"] directions = np.array(lon[1:], dtype=np.float64) - assert directions.shape == (len(cells), 6) - np.testing.assert_allclose(directions[:, :3], normalised(fiber_expr(midpoints.T / 1000)).T) - np.testing.assert_allclose(directions[:, 3:], normalised(sheet_expr(midpoints.T / 1000)).T) + assert directions.shape == (num_cells, 6) + for column, expr in ((0, fiber_expr), (3, sheet_expr)): + expected = expected_direction(expr, midpoints / 1000, tdim) + np.testing.assert_allclose(directions[:, column : column + 3], expected, atol=1e-9) for name, plane in (("left", 0.0), ("right", 1000.0)): - vtx = read_lines(filename.parent / f"cube_{name}.vtx") + vtx = read_lines(filename.parent / f"mesh_{name}.vtx") assert vtx[1] == ["extra"] nodes = np.array(vtx[2:], dtype=np.int64).reshape(-1) - assert int(vtx[0][0]) == len(nodes) == 16 + assert int(vtx[0][0]) == len(nodes) == num_marked_nodes assert np.allclose(x[nodes, 0], plane) - surf = read_lines(filename.parent / f"cube_{name}.surf") - assert int(surf[0][0]) == len(surf) - 1 == 2 * 9 - assert {row[0] for row in surf[1:]} == {"Tr"} - triangles = np.array([row[1:] for row in surf[1:]], dtype=np.int64) - np.testing.assert_array_equal(np.unique(triangles), np.sort(nodes)) - # Normals point out of the cube - q = x[triangles] - normals = np.cross(q[:, 1] - q[:, 0], q[:, 2] - q[:, 0]) - assert np.all(np.sign(normals[:, 0]) == (-1 if plane == 0 else 1)) - - -def test_files_independent_of_partitioning(shared_tmp_path): - parallel = shared_tmp_path / "parallel" / "cube" - serial = shared_tmp_path / "serial" / "cube" + surf_file = filename.parent / f"mesh_{name}.surf" + if tdim == 2: + assert not surf_file.exists() + continue + surf = read_lines(surf_file) + assert int(surf[0][0]) == len(surf) - 1 == num_facets + face_type = {"Tt": "Tr", "Hx": "Qd"}[elem_type] + assert {row[0] for row in surf[1:]} == {face_type} + faces = np.array([row[1:] for row in surf[1:]], dtype=np.int64) + np.testing.assert_array_equal(np.unique(faces), np.sort(nodes)) + # Normals point out of the cube, and quadrilateral nodes are in cyclic order + q = x[faces] + if face_type == "Tr": + normals = np.cross(q[:, 1] - q[:, 0], q[:, 2] - q[:, 0]) / 2 + area = (1000 / 3) ** 2 / 2 + else: + normals = np.cross(q[:, 2] - q[:, 0], q[:, 3] - q[:, 1]) / 2 + area = (1000 / 3) ** 2 + np.testing.assert_allclose(normals[:, 0], (-1 if plane == 0 else 1) * area) + + +@cell_types +def test_files_independent_of_partitioning(shared_tmp_path, cell_type): + parallel = shared_tmp_path / "parallel" / "mesh" + serial = shared_tmp_path / "serial" / "mesh" if MPI.COMM_WORLD.rank == 0: parallel.parent.mkdir() serial.parent.mkdir() MPI.COMM_WORLD.Barrier() - write_problem(parallel, MPI.COMM_WORLD) + write_problem(parallel, MPI.COMM_WORLD, cell_type) if MPI.COMM_WORLD.rank == 0: - write_problem(serial, MPI.COMM_SELF) + write_problem(serial, MPI.COMM_SELF, cell_type) files = sorted(p.name for p in serial.parent.iterdir()) assert files == sorted(p.name for p in parallel.parent.iterdir()) - assert len(files) == 7 + # .pts, .elem, .lon and a .vtx (and in 3D a .surf) per marker + assert len(files) == (7 if dolfinx.mesh.cell_dim(cell_type) == 3 else 5) for name in files: assert (serial.parent / name).read_bytes() == (parallel.parent / name).read_bytes() -def test_read_write(shared_tmp_path): - write_problem(shared_tmp_path / "cube", MPI.COMM_WORLD) - filename = shared_tmp_path / "cube.elem" +@cell_types +def test_read_write(shared_tmp_path, cell_type): + write_problem(shared_tmp_path / "mesh", MPI.COMM_WORLD, cell_type) + filename = shared_tmp_path / "mesh.elem" comm = MPI.COMM_WORLD - backend_args = {"scale": 1000.0, "markers": markers} + _, num_cells, num_nodes, num_facets, num_marked_nodes = meshes[cell_type] + tdim = dolfinx.mesh.cell_dim(cell_type) + backend_args = {"scale": 1000.0, "markers": markers(tdim)} mesh = io4dolfinx.read_mesh(filename, comm, backend="opencarp", backend_args=backend_args) - tdim = mesh.topology.dim - assert mesh.topology.index_map(tdim).size_global == 6 * 27 - assert mesh.geometry.index_map().size_global == 64 + assert mesh.topology.cell_type == cell_type + assert mesh.geometry.dim == tdim + assert mesh.topology.index_map(tdim).size_global == num_cells + assert mesh.geometry.index_map().size_global == num_nodes volume = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1 * ufl.dx(domain=mesh))) assert np.isclose(comm.allreduce(volume, op=MPI.SUM), 1.0) @@ -163,44 +229,50 @@ def test_read_write(shared_tmp_path): midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, cell_tags.indices) np.testing.assert_array_equal(cell_tags.values, np.where(midpoints[:, 0] < 0.4, 1, 2)) num_owned = mesh.topology.index_map(tdim).size_local - assert comm.allreduce(np.count_nonzero(cell_tags.indices < num_owned)) == 6 * 27 - - facet_tags = io4dolfinx.read_meshtags( - filename, mesh, "facet", backend="opencarp", backend_args=backend_args - ) - midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim - 1, facet_tags.indices) - np.testing.assert_allclose(midpoints[:, 0], np.where(facet_tags.values == 1, 0.0, 1.0)) - num_owned = mesh.topology.index_map(tdim - 1).size_local - assert comm.allreduce(np.count_nonzero(facet_tags.indices < num_owned)) == 2 * 18 + assert comm.allreduce(np.count_nonzero(cell_tags.indices < num_owned)) == num_cells + # Facets are read from .surf files in 3D, marked nodes from .vtx files in 2D right = io4dolfinx.read_meshtags( filename, mesh, "right", backend="opencarp", backend_args=backend_args ) - assert right.dim == tdim - 1 + dim = tdim - 1 if tdim == 3 else 0 + assert right.dim == dim assert np.all(right.values == 2) - assert np.allclose(dolfinx.mesh.compute_midpoints(mesh, tdim - 1, right.indices)[:, 0], 1) + assert np.allclose(dolfinx.mesh.compute_midpoints(mesh, dim, right.indices)[:, 0], 1) + num_owned = mesh.topology.index_map(dim).size_local + num_marked = num_facets if tdim == 3 else num_marked_nodes + assert comm.allreduce(np.count_nonzero(right.indices < num_owned)) == num_marked + if tdim == 3: + facet_tags = io4dolfinx.read_meshtags( + filename, mesh, "facet", backend="opencarp", backend_args=backend_args + ) + midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim - 1, facet_tags.indices) + np.testing.assert_allclose(midpoints[:, 0], np.where(facet_tags.values == 1, 0.0, 1.0)) + num_owned = mesh.topology.index_map(tdim - 1).size_local + assert comm.allreduce(np.count_nonzero(facet_tags.indices < num_owned)) == 2 * num_facets assert io4dolfinx.read_function_names(filename, comm, backend="opencarp") == ["fiber", "sheet"] for name, expr in (("fiber", fiber_expr), ("sheet", sheet_expr)): u = io4dolfinx.read_cell_data(filename, name, mesh, backend="opencarp") - num_cells = mesh.topology.index_map(tdim).size_local - cells = np.arange(num_cells, dtype=np.int32) + assert u.ufl_shape == (tdim,) + cells = np.arange(mesh.topology.index_map(tdim).size_local, dtype=np.int32) midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, cells) - values = u.x.array.reshape(-1, 3)[u.function_space.dofmap.list[cells, 0]] - np.testing.assert_allclose(values, normalised(expr(midpoints.T)).T) + values = u.x.array.reshape(-1, tdim)[u.function_space.dofmap.list[cells, 0]] + expected = expected_direction(expr, midpoints, tdim)[:, :tdim] + np.testing.assert_allclose(values, expected, atol=1e-9) def test_read_vertex_markers(shared_tmp_path): """Without ``.surf`` files a marker is read from its ``.vtx`` file, as vertex tags.""" - write_problem(shared_tmp_path / "cube", MPI.COMM_WORLD, backend_args={"surf": False}) - filename = shared_tmp_path / "cube.elem" - assert not (shared_tmp_path / "cube_left.surf").exists() + write_problem(shared_tmp_path / "mesh", MPI.COMM_WORLD, backend_args={"surf": False}) + filename = shared_tmp_path / "mesh.elem" + assert not (shared_tmp_path / "mesh_left.surf").exists() mesh = io4dolfinx.read_mesh( filename, MPI.COMM_WORLD, backend="opencarp", backend_args={"scale": 1000.0} ) left = io4dolfinx.read_meshtags( - filename, mesh, "left", backend="opencarp", backend_args={"markers": markers} + filename, mesh, "left", backend="opencarp", backend_args={"markers": markers(3)} ) assert left.dim == 0 assert np.all(left.values == 1) @@ -211,7 +283,7 @@ def test_read_vertex_markers(shared_tmp_path): def test_fiber_only_and_unnamed_markers(shared_tmp_path): """Markers without names are written and read back by value.""" - filename = shared_tmp_path / "cube" + filename = shared_tmp_path / "mesh" mesh, _, facet_tags, fiber, _ = create_problem(MPI.COMM_WORLD) write_opencarp(filename, mesh, scale=1.0, facet_tags=facet_tags, fiber=fiber) if MPI.COMM_WORLD.rank == 0: @@ -234,7 +306,7 @@ def test_fiber_only_and_unnamed_markers(shared_tmp_path): def test_errors(shared_tmp_path): comm = MPI.COMM_WORLD - filename = shared_tmp_path / "cube" + filename = shared_tmp_path / "mesh" mesh, cell_tags, _, fiber, sheet = create_problem(comm) with pytest.raises(ValueError, match="scale"): @@ -250,11 +322,9 @@ def test_errors(shared_tmp_path): with pytest.raises(ValueError, match="quadrature"): write_opencarp(filename, mesh, scale=1.0, fiber=dolfinx.fem.Function(Q)) - hex_mesh = dolfinx.mesh.create_unit_cube( - comm, 2, 2, 2, cell_type=dolfinx.mesh.CellType.hexahedron - ) - with pytest.raises(NotImplementedError, match="tetrahedral"): - write_opencarp(filename, hex_mesh, scale=1.0) + interval_mesh = dolfinx.mesh.create_unit_interval(comm, 4) + with pytest.raises(NotImplementedError, match="interval"): + write_opencarp(filename, interval_mesh, scale=1.0) # Errors found when writing on rank 0 are raised on all ranks write_opencarp(filename, mesh, scale=1.0) @@ -269,3 +339,19 @@ def test_errors(shared_tmp_path): io4dolfinx.write_cell_data( filename, fiber, 0.0, io4dolfinx.FileMode.write, None, "opencarp" ) + + +def test_2d_surface_not_supported(shared_tmp_path): + """Triangles that do not lie in a plane of constant z (surfaces in 3D) are rejected.""" + filename = shared_tmp_path / "surface" + if MPI.COMM_WORLD.rank == 0: + filename.with_suffix(".pts").write_text("3\n0 0 0\n1 0 0\n0 1 1\n") + filename.with_suffix(".elem").write_text("1\nTr 0 1 2 0\n") + MPI.COMM_WORLD.Barrier() + with pytest.raises(NotImplementedError, match="constant z"): + io4dolfinx.read_mesh( + filename.with_suffix(".elem"), + MPI.COMM_WORLD, + backend="opencarp", + backend_args={"scale": 1.0}, + ) From 682c77751c43c49eb0fd72bc81ef2f411151292a Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 28 Oct 2026 12:08:48 +0000 Subject: [PATCH 5/8] Fix return type of orientation helper in openCARP tests Co-Authored-By: Claude Opus 5.5 --- tests/test_opencarp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_opencarp.py b/tests/test_opencarp.py index 5090fce..5a1f5b8 100644 --- a/tests/test_opencarp.py +++ b/tests/test_opencarp.py @@ -99,7 +99,7 @@ def read_lines(path: Path) -> list[list[str]]: return [line.split() for line in path.read_text().splitlines()] -def orientation(elem_type: str, p: np.ndarray) -> np.ndarray: +def orientation(elem_type: str, p: np.ndarray) -> tuple[np.ndarray, int]: """Signed measure of each element with nodes p in openCARP ordering, positive for the orientation of openCARP's mesher, and the factor to its volume or area.""" e1 = p[:, 1] - p[:, 0] From f218d9bbacf28e0913b91650ccef7b5bbbc7d4c3 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 30 Oct 2026 14:15:54 +0000 Subject: [PATCH 6/8] Write openCARP files with the generic API in either ordering Remove io4dolfinx.backends.opencarp.write_opencarp. The backend module documentation now describes how to write a mesh, its tags and its fiber and sheet directions with write_mesh_input_order (or write_mesh), write_meshtags and write_cell_data. The backend accepts data in the DOLFINx ordering as well as the input ordering, as long as every step uses the same one. Only the input ordering gives files that do not depend on the number of processes. The on_input_mesh fields of MeshData, MeshTagsData and ArrayData are no longer needed and are removed. write_cell_data raises a clear error for functions in a quadrature space, which cannot be interpolated into DG-0. The openCARP tests run for both orderings. Co-Authored-By: Claude Opus 5.5 --- README.md | 2 +- docs/backends/opencarp.rst | 45 ++------ src/io4dolfinx/backends/opencarp/__init__.py | 3 +- src/io4dolfinx/backends/opencarp/backend.py | 74 ++++++++----- src/io4dolfinx/backends/opencarp/wrapper.py | 100 ------------------ src/io4dolfinx/checkpointing.py | 8 +- src/io4dolfinx/original_checkpoint.py | 1 - src/io4dolfinx/structures.py | 7 -- tests/test_opencarp.py | 103 ++++++++++++------- tests/test_write_data_input_order.py | 2 - 10 files changed, 131 insertions(+), 214 deletions(-) delete mode 100644 src/io4dolfinx/backends/opencarp/wrapper.py diff --git a/README.md b/README.md index b439e34..0a6ad8d 100644 --- a/README.md +++ b/README.md @@ -71,7 +71,7 @@ Most meshing formats supports associating data with the nodes of the mesh (the m - [VTKHDF](./docs/backends/vtkhdf.rst): The new scalable format from VTK, called [VTKHDF](https://docs.vtk.org/en/latest/vtk_file_formats/vtkhdf_file_format/index.html) is supported by the `vtkhdf` backend. - [XDMF](./docs/backends/xdmf.rst) (eXtensible Model Data Format): `.xdmf`. The `xdmf` backend supports the `HDF5` encoding, to ensure performance in parallel. - [Gmsh](./docs/backends/gmsh.rst): `.msh` files (MSH 2.2 and 4.1, ASCII and binary) are read and written with the [Gmsh](https://gmsh.info) Python API by the `gmsh` backend, including physical groups as mesh tags and time-dependent point and cell data. -- [openCARP](./docs/backends/opencarp.rst): Linear triangle, quadrilateral, tetrahedral and hexahedral meshes with region tags, fibers and markers are written to and read from the ASCII format of the cardiac simulator [openCARP](https://opencarp.org) (`.pts`, `.elem`, `.lon`, `.vtx`, `.surf`) by the `opencarp` backend, or in one call with `io4dolfinx.backends.opencarp.write_opencarp`. +- [openCARP](./docs/backends/opencarp.rst): Linear triangle, quadrilateral, tetrahedral and hexahedral meshes with region tags, fibers and markers are written to and read from the ASCII format of the cardiac simulator [openCARP](https://opencarp.org) (`.pts`, `.elem`, `.lon`, `.vtx`, `.surf`) by the `opencarp` backend. - [PyVista](./docs/backends/pyvista.rst) (IO backend is meshio): The [pyvista](https://pyvista.org/) backend uses {py:func}`pyvista.read` to read in meshes, point data and cell data. `pyvista` relies on [meshio](https://github.com/nschloe/meshio) for most reading operations (including the XDMF ascii format). diff --git a/docs/backends/opencarp.rst b/docs/backends/opencarp.rst index 16bc1a8..f4bbb4e 100644 --- a/docs/backends/opencarp.rst +++ b/docs/backends/opencarp.rst @@ -8,29 +8,15 @@ Writes and reads linear meshes in the ASCII mesh format of the cardiac simulator (``.elem``), fiber and sheet directions per element (``.lon``), and markers as vertex (``.vtx``) and surface (``.surf``) files. The backend needs no extra dependencies. -Nodes and elements are written in the ordering of the input mesh, gathered to and written by -rank 0, so the files are the same whatever the number of processes. openCARP expects coordinates -in µm, so the factor converting the mesh coordinates to µm must always be given. +The mesh, its tags and its fiber and sheet directions are written with +:func:`io4dolfinx.write_mesh_input_order` (or :func:`io4dolfinx.write_mesh`), +:func:`io4dolfinx.write_meshtags` and :func:`io4dolfinx.write_cell_data`, as shown in the module +documentation below. Rank 0 writes all files. openCARP expects coordinates in µm, so the factor +converting the mesh coordinates to µm must always be given as ``backend_args["scale"]``. -:func:`io4dolfinx.backends.opencarp.write_opencarp` writes everything in one call: - -.. code-block:: python - - from io4dolfinx.backends.opencarp import write_opencarp - - write_opencarp( - "heart", - mesh, - scale=1000.0, # mesh in mm - cell_tags=cell_tags, - facet_tags=facet_tags, - fiber=f0, - sheet=s0, - markers={"ENDO": (1, 2), "EPI": (2, 2)}, - ) - -This gives ``heart.pts``, ``heart.elem``, ``heart.lon``, and ``heart_ENDO.vtx``, -``heart_ENDO.surf``, ``heart_EPI.vtx`` and ``heart_EPI.surf``. +With ``markers={"ENDO": (1, 2), "EPI": (2, 2)}``, a heart mesh gives ``heart.pts``, +``heart.elem``, ``heart.lon``, and ``heart_ENDO.vtx``, ``heart_ENDO.surf``, ``heart_EPI.vtx`` and +``heart_EPI.surf``. - Triangle (``Tr``) and quadrilateral (``Qd``) meshes in 2D, and tetrahedral (``Tt``) and hexahedral (``Hx``) meshes in 3D are supported. 2D meshes are written with ``z = 0``, and read back with 2D @@ -46,21 +32,6 @@ This gives ``heart.pts``, ``heart.elem``, ``heart.lon``, and ``heart_ENDO.vtx``, ``markers``, or by the value if it has no name. Pass ``backend_args={"surf": False}`` to skip the ``.surf`` files. -The same files can be written with the generic functions, as long as everything is in the input -ordering: - -.. code-block:: python - - args = {"scale": 1000.0} - io4dolfinx.write_mesh_input_order("heart", mesh, backend="opencarp", backend_args=args) - io4dolfinx.write_meshtags("heart", mesh, cell_tags, backend="opencarp", on_input_mesh=True) - io4dolfinx.write_cell_data( - "heart", f0, 0.0, io4dolfinx.FileMode.write, args, "opencarp", on_input_mesh=True - ) - io4dolfinx.write_cell_data( - "heart", s0, 0.0, io4dolfinx.FileMode.append, args, "opencarp", on_input_mesh=True - ) - Reading takes the path of any of the mesh files. Coordinates are divided by ``scale``: .. code-block:: python diff --git a/src/io4dolfinx/backends/opencarp/__init__.py b/src/io4dolfinx/backends/opencarp/__init__.py index 0032034..05458f8 100644 --- a/src/io4dolfinx/backends/opencarp/__init__.py +++ b/src/io4dolfinx/backends/opencarp/__init__.py @@ -1,4 +1,3 @@ from . import backend -from .wrapper import write_opencarp -__all__ = ["backend", "write_opencarp"] +__all__ = ["backend"] diff --git a/src/io4dolfinx/backends/opencarp/backend.py b/src/io4dolfinx/backends/opencarp/backend.py index 5142176..fb1b0ca 100644 --- a/src/io4dolfinx/backends/opencarp/backend.py +++ b/src/io4dolfinx/backends/opencarp/backend.py @@ -19,8 +19,48 @@ - ``_.surf``: number of faces, then one ``Tr n0 n1 n2`` or ``Qd n0 n1 n2 n3`` line per face of a 3D mesh. -Nodes and cells are written in the ordering of the input mesh, so the files do -not depend on the number of processes. Each file is written by rank 0. +Each file is written by rank 0. + +Writing a mesh with markers and fibers takes the generic io4dolfinx functions, called +in this order with ``backend="opencarp"``: + +.. code-block:: python + + args = {"scale": 1000.0, "markers": {"ENDO": (1, 2), "EPI": (2, 2)}} # mesh in mm + io4dolfinx.write_mesh_input_order("heart", mesh, backend="opencarp", backend_args=args) + for tags in (cell_tags, facet_tags): + io4dolfinx.write_meshtags( + "heart", mesh, tags, backend="opencarp", backend_args=args, on_input_mesh=True + ) + io4dolfinx.write_cell_data( + "heart", f0, 0.0, io4dolfinx.FileMode.write, args, "opencarp", on_input_mesh=True + ) + io4dolfinx.write_cell_data( + "heart", s0, 0.0, io4dolfinx.FileMode.append, args, "opencarp", on_input_mesh=True + ) + +- The mesh comes first, as tags and directions are written against ``.elem``. +- Cell tags set the regions in ``.elem``, facet tags give a ``.vtx`` (and in + 3D a ``.surf``) file per marker value. +- :func:`io4dolfinx.write_cell_data` with :attr:`FileMode.write + ` writes the fibers, and with :attr:`FileMode.append + ` adds the sheets. Fields are interpolated into DG-0, + so fields in a quadrature space must first be generated in a P-1 or DG-0 space. + +The ordering of nodes and elements must be the same in every step: + +- With :func:`io4dolfinx.write_mesh_input_order` and ``on_input_mesh=True`` (as + above), the files use the ordering of the input mesh and are the same whatever the + number of processes. Use this for reproducible files. +- With :func:`io4dolfinx.write_mesh` and ``on_input_mesh=False`` (the default of + :func:`io4dolfinx.write_meshtags` and :func:`io4dolfinx.write_cell_data`), the files + use the DOLFINx ordering, which depends on the partitioning. + +Mixing the two gives directions and markers on the wrong elements and nodes. + +The files are read with :func:`io4dolfinx.read_mesh`, :func:`io4dolfinx.read_meshtags` +and :func:`io4dolfinx.read_cell_data`, given the path of any of the mesh files, e.g. +``"heart.elem"``. Elements are written with the orientation used by openCARP's ``mesher``: positive volume in 3D, counter-clockwise in the xy-plane in 2D. Hexahedra list one face @@ -554,19 +594,13 @@ def write_mesh( Args: filename: Base name of the files, or the path of one of them comm: MPI communicator used in storage - mesh: Mesh data in input ordering, see - :func:`io4dolfinx.write_mesh_input_order` + mesh: Mesh data, from :func:`io4dolfinx.write_mesh_input_order` or + :func:`io4dolfinx.write_mesh` backend_args: Arguments to backend, ``"scale"`` is required mode: Only :attr:`FileMode.write` is supported time: Not used """ args = get_default_backend_args(backend_args) - if not mesh.on_input_mesh: - raise ValueError( - "The openCARP backend writes meshes in the input ordering, so that the files do " - "not depend on the number of processes. Use io4dolfinx.write_mesh_input_order or " - "io4dolfinx.backends.opencarp.write_opencarp." - ) if mesh.cell_type not in _element_types or mesh.degree != 1: raise NotImplementedError( "The openCARP backend supports linear triangle, quadrilateral, tetrahedral and " @@ -681,17 +715,11 @@ def write_meshtags( Args: filename: Base name of the files, or the path of one of them comm: MPI communicator used in storage - data: Mesh tags with nodes in input ordering, see the ``on_input_mesh`` - argument of :func:`io4dolfinx.write_meshtags` + data: Mesh tags with nodes in the ordering the mesh was written in, see the + ``on_input_mesh`` argument of :func:`io4dolfinx.write_meshtags` backend_args: Arguments to backend """ args = get_default_backend_args(backend_args) - if not data.on_input_mesh: - raise ValueError( - "The openCARP backend writes meshtags in the input node ordering. Use " - "io4dolfinx.write_meshtags(..., on_input_mesh=True) or " - "io4dolfinx.backends.opencarp.write_opencarp." - ) local_indices = np.asarray(data.indices, dtype=np.int64) num_nodes = data.num_dofs_per_entity or local_indices.shape[-1] indices = _gather_concatenate(comm, local_indices.reshape(-1, num_nodes)) @@ -1023,8 +1051,8 @@ def write_data( Args: filename: Base name of the files, or the path of one of them - array_data: Cell data in input ordering, see the ``on_input_mesh`` argument - of :func:`io4dolfinx.write_cell_data` + array_data: Cell data in the ordering the mesh was written in, see the + ``on_input_mesh`` argument of :func:`io4dolfinx.write_cell_data` comm: MPI communicator to open the file with time: Not used mode: Write for the fiber, append for the sheet direction @@ -1034,12 +1062,6 @@ def write_data( raise NotImplementedError( "The openCARP backend can only write cell data (fibers and sheets)." ) - if not array_data.on_input_mesh: - raise ValueError( - "The openCARP backend writes cell data in the input cell ordering. Use " - "io4dolfinx.write_cell_data(..., on_input_mesh=True) or " - "io4dolfinx.backends.opencarp.write_opencarp." - ) num_cells, num_components = array_data.global_shape if num_components not in (2, 3): raise ValueError(f"Expected a 2D or 3D vector per cell, got {num_components} components.") diff --git a/src/io4dolfinx/backends/opencarp/wrapper.py b/src/io4dolfinx/backends/opencarp/wrapper.py deleted file mode 100644 index 8f41672..0000000 --- a/src/io4dolfinx/backends/opencarp/wrapper.py +++ /dev/null @@ -1,100 +0,0 @@ -# Copyright (C) 2026 Henrik Finsberg -# -# This file is part of io4dolfinx -# -# SPDX-License-Identifier: MIT - -"""Write a mesh, its markers and fibers to openCARP in one call.""" - -from __future__ import annotations - -import logging -from pathlib import Path -from typing import Any - -import dolfinx - -from ...checkpointing import write_cell_data, write_meshtags -from ...original_checkpoint import write_mesh_input_order -from .. import FileMode - -__all__ = ["write_opencarp"] -logger = logging.getLogger(__name__) - - -def _check_direction(u: dolfinx.fem.Function, name: str): - """Check that a fiber or sheet field can be interpolated into cell midpoints.""" - if u.function_space.ufl_element().family_name == "quadrature": - raise ValueError( - f"Cannot write the {name} field {u.name!r} from a quadrature space to openCARP, as it " - "cannot be interpolated into DG-0. Generate the fibers in a P-1 or DG-0 space." - ) - gdim = u.function_space.mesh.geometry.dim - if u.ufl_shape not in ((gdim,), (3,)): - raise ValueError( - f"Expected the {name} field to be a {gdim}D or 3D vector, got shape {u.ufl_shape}." - ) - - -def write_opencarp( - filename: Path | str, - mesh: dolfinx.mesh.Mesh, - scale: float, - cell_tags: dolfinx.mesh.MeshTags | None = None, - facet_tags: dolfinx.mesh.MeshTags | None = None, - fiber: dolfinx.fem.Function | None = None, - sheet: dolfinx.fem.Function | None = None, - markers: dict[str, Any] | None = None, - backend_args: dict[str, Any] | None = None, -): - """Write a linear mesh with markers and fibers to openCARP files. - - Writes ``.pts`` and ``.elem``, with the region tag of each element from - ``cell_tags``; ``.lon`` with the fiber (and sheet) direction of each element; - and ``_.vtx`` and ``_.surf`` for each value of - ``facet_tags`` (``.surf`` only for 3D meshes). Nodes and elements are in the - ordering of the input mesh, so the files do not depend on the number of processes. See - :mod:`io4dolfinx.backends.opencarp.backend` for the format. - - Args: - filename: Base name of the files - mesh: The mesh - scale: Factor converting mesh coordinates to µm, e.g. ``1000.0`` for a mesh in mm - cell_tags: Region tag of each cell. Untagged cells get - ``backend_args["default_tag"]`` (``0`` by default). - facet_tags: Facet markers, e.g. for stimuli and boundary conditions - fiber: Fiber direction, interpolated into cell midpoints and normalised - sheet: Sheet direction, written next to the fibers - markers: Map from marker name to value, or to ``(value, dim)``, naming the - files of each facet marker - backend_args: Further arguments to the openCARP backend - """ - if sheet is not None and fiber is None: - raise ValueError("Cannot write a sheet direction without a fiber direction.") - for u, name in ((fiber, "fiber"), (sheet, "sheet")): - if u is not None: - _check_direction(u, name) - - args = dict(backend_args or {}) - args["scale"] = scale - if markers is not None: - args["markers"] = markers - - logger.debug(f"Writing openCARP mesh to {filename} with arguments {args}") - write_mesh_input_order(filename, mesh, backend="opencarp", backend_args=args) - for tags in (cell_tags, facet_tags): - if tags is not None: - write_meshtags( - filename, mesh, tags, backend="opencarp", backend_args=args, on_input_mesh=True - ) - for u, mode in ((fiber, FileMode.write), (sheet, FileMode.append)): - if u is not None: - write_cell_data( - filename, - u, - time=0.0, - mode=mode, - backend_args=args, - backend="opencarp", - on_input_mesh=True, - ) diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index 7f8ed63..9c83558 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -265,7 +265,6 @@ def write_meshtags( local_start=local_start, dim=meshtags.dim, cell_type=tag_ct, - on_input_mesh=on_input_mesh, ) # Get backend and default arguments @@ -795,7 +794,6 @@ def write_point_data( global_shape=data_shape, local_range=local_range, type="Point", - on_input_mesh=on_input_mesh, ) logger.debug( f"Using {backend} backend with arguments {backend_args} and mode {mode} to write point data" @@ -832,6 +830,11 @@ def write_cell_data( """ backend_args = _with_mesh_name(backend_args, mesh_name) logger.debug(f"Writing cell data to {filename} for function {u.name} at time {time}") + if u.function_space.ufl_element().family_name == "quadrature": + raise ValueError( + f"Cannot write {u.name!r} as cell data, as a function in a quadrature space cannot " + "be interpolated into DG-0." + ) V = dolfinx.fem.functionspace(u.function_space.mesh, ("DG", 0, u.ufl_shape)) v_out = dolfinx.fem.Function(V, name=u.name, dtype=u.x.array.dtype) v_out.interpolate(u) @@ -859,7 +862,6 @@ def write_cell_data( global_shape=data_shape, local_range=local_range, type="Cell", - on_input_mesh=on_input_mesh, ) logger.debug( f"Using {backend} backend with arguments {backend_args} and mode {mode} to write cell data" diff --git a/src/io4dolfinx/original_checkpoint.py b/src/io4dolfinx/original_checkpoint.py index e9a0263..39d1dc8 100644 --- a/src/io4dolfinx/original_checkpoint.py +++ b/src/io4dolfinx/original_checkpoint.py @@ -204,7 +204,6 @@ def create_original_mesh_data(mesh: dolfinx.mesh.Mesh) -> MeshData: ownership_offset=None, partition_range=None, partition_global=None, - on_input_mesh=True, ) diff --git a/src/io4dolfinx/structures.py b/src/io4dolfinx/structures.py index 83349a1..e19e97e 100644 --- a/src/io4dolfinx/structures.py +++ b/src/io4dolfinx/structures.py @@ -45,9 +45,6 @@ class MeshData: None #: Local insert position for partitioning information ) partition_global: int | None = None #: Global size of partitioning array - #: Nodes and cells use the ordering of the input mesh rather than the DOLFINx global - #: ordering, which depends on the partitioning - on_input_mesh: bool = False @dataclass @@ -101,8 +98,6 @@ class MeshTagsData: # Optional info to help visualization cell_type: str | None = None #: The cell type - #: Entity nodes use the node ordering of the input mesh - on_input_mesh: bool = False @dataclass @@ -112,5 +107,3 @@ class ArrayData: global_shape: tuple[int, int] # Global shape of data local_range: tuple[int, int] # Start and stop of local array in global array type: str # Type of array-data - #: Rows use the cell ordering of the input mesh (only for cell data) - on_input_mesh: bool = False diff --git a/tests/test_opencarp.py b/tests/test_opencarp.py index 5a1f5b8..35f8abb 100644 --- a/tests/test_opencarp.py +++ b/tests/test_opencarp.py @@ -9,7 +9,6 @@ import ufl import io4dolfinx -from io4dolfinx.backends.opencarp import write_opencarp CellType = dolfinx.mesh.CellType @@ -22,6 +21,9 @@ CellType.quadrilateral: ("Qd", 9, 16, 3, 4), } cell_types = pytest.mark.parametrize("cell_type", list(meshes), ids=lambda c: c.name) +orderings = pytest.mark.parametrize( + "on_input_mesh", [True, False], ids=["input_ordering", "dolfinx_ordering"] +) def markers(tdim: int): @@ -75,24 +77,50 @@ def create_problem(comm: MPI.Intracomm, cell_type: dolfinx.mesh.CellType = CellT return mesh, cell_tags, facet_tags, fiber, sheet +def write_files( + filename: Path, + mesh: dolfinx.mesh.Mesh, + backend_args: dict, + on_input_mesh: bool, + cell_tags: dolfinx.mesh.MeshTags | None = None, + facet_tags: dolfinx.mesh.MeshTags | None = None, + fiber: dolfinx.fem.Function | None = None, + sheet: dolfinx.fem.Function | None = None, +): + """Write openCARP files as described in the backend documentation.""" + if on_input_mesh: + io4dolfinx.write_mesh_input_order( + filename, mesh, backend="opencarp", backend_args=backend_args + ) + else: + io4dolfinx.write_mesh(filename, mesh, backend="opencarp", backend_args=backend_args) + for tags in (cell_tags, facet_tags): + if tags is not None: + io4dolfinx.write_meshtags( + filename, + mesh, + tags, + backend="opencarp", + backend_args=backend_args, + on_input_mesh=on_input_mesh, + ) + for u, mode in ((fiber, io4dolfinx.FileMode.write), (sheet, io4dolfinx.FileMode.append)): + if u is not None: + io4dolfinx.write_cell_data( + filename, u, 0.0, mode, backend_args, "opencarp", on_input_mesh=on_input_mesh + ) + + def write_problem( filename: Path, comm: MPI.Intracomm, cell_type: dolfinx.mesh.CellType = CellType.tetrahedron, - **kwargs, + on_input_mesh: bool = True, + backend_args: dict | None = None, ): mesh, cell_tags, facet_tags, fiber, sheet = create_problem(comm, cell_type) - write_opencarp( - filename, - mesh, - scale=1000.0, - cell_tags=cell_tags, - facet_tags=facet_tags, - fiber=fiber, - sheet=sheet, - markers=markers(mesh.topology.dim), - **kwargs, - ) + args = {"scale": 1000.0, "markers": markers(mesh.topology.dim), **(backend_args or {})} + write_files(filename, mesh, args, on_input_mesh, cell_tags, facet_tags, fiber, sheet) def read_lines(path: Path) -> list[list[str]]: @@ -110,10 +138,11 @@ def orientation(elem_type: str, p: np.ndarray) -> tuple[np.ndarray, int]: return np.einsum("ij,ij->i", np.cross(e1, e2), e3), 6 if elem_type == "Tt" else 1 +@orderings @cell_types -def test_write_content(shared_tmp_path, cell_type): +def test_write_content(shared_tmp_path, cell_type, on_input_mesh): filename = shared_tmp_path / "mesh" - write_problem(filename, MPI.COMM_WORLD, cell_type) + write_problem(filename, MPI.COMM_WORLD, cell_type, on_input_mesh) if MPI.COMM_WORLD.rank != 0: return elem_type, num_cells, num_nodes, num_facets, num_marked_nodes = meshes[cell_type] @@ -191,6 +220,7 @@ def test_write_content(shared_tmp_path, cell_type): @cell_types def test_files_independent_of_partitioning(shared_tmp_path, cell_type): + """With the input ordering, the files do not depend on the number of processes.""" parallel = shared_tmp_path / "parallel" / "mesh" serial = shared_tmp_path / "serial" / "mesh" if MPI.COMM_WORLD.rank == 0: @@ -208,9 +238,10 @@ def test_files_independent_of_partitioning(shared_tmp_path, cell_type): assert (serial.parent / name).read_bytes() == (parallel.parent / name).read_bytes() +@orderings @cell_types -def test_read_write(shared_tmp_path, cell_type): - write_problem(shared_tmp_path / "mesh", MPI.COMM_WORLD, cell_type) +def test_read_write(shared_tmp_path, cell_type, on_input_mesh): + write_problem(shared_tmp_path / "mesh", MPI.COMM_WORLD, cell_type, on_input_mesh) filename = shared_tmp_path / "mesh.elem" comm = MPI.COMM_WORLD _, num_cells, num_nodes, num_facets, num_marked_nodes = meshes[cell_type] @@ -285,7 +316,7 @@ def test_fiber_only_and_unnamed_markers(shared_tmp_path): """Markers without names are written and read back by value.""" filename = shared_tmp_path / "mesh" mesh, _, facet_tags, fiber, _ = create_problem(MPI.COMM_WORLD) - write_opencarp(filename, mesh, scale=1.0, facet_tags=facet_tags, fiber=fiber) + write_files(filename, mesh, {"scale": 1.0}, True, facet_tags=facet_tags, fiber=fiber) if MPI.COMM_WORLD.rank == 0: assert read_lines(filename.with_suffix(".lon"))[0] == ["1"] tags = {row[5] for row in read_lines(filename.with_suffix(".elem"))[1:]} @@ -307,38 +338,40 @@ def test_fiber_only_and_unnamed_markers(shared_tmp_path): def test_errors(shared_tmp_path): comm = MPI.COMM_WORLD filename = shared_tmp_path / "mesh" - mesh, cell_tags, _, fiber, sheet = create_problem(comm) + mesh, _, _, fiber, _ = create_problem(comm) with pytest.raises(ValueError, match="scale"): io4dolfinx.write_mesh_input_order(filename, mesh, backend="opencarp") - with pytest.raises(ValueError, match="input ordering"): - io4dolfinx.write_mesh(filename, mesh, backend="opencarp", backend_args={"scale": 1.0}) - with pytest.raises(ValueError, match="without a fiber"): - write_opencarp(filename, mesh, scale=1.0, sheet=sheet) + interval_mesh = dolfinx.mesh.create_unit_interval(comm, 4) + with pytest.raises(NotImplementedError, match="interval"): + io4dolfinx.write_mesh_input_order( + filename, interval_mesh, backend="opencarp", backend_args={"scale": 1.0} + ) + + io4dolfinx.write_mesh_input_order( + filename, mesh, backend="opencarp", backend_args={"scale": 1.0} + ) Q = dolfinx.fem.functionspace( mesh, basix.ufl.quadrature_element(mesh.basix_cell(), value_shape=(3,), degree=2) ) with pytest.raises(ValueError, match="quadrature"): - write_opencarp(filename, mesh, scale=1.0, fiber=dolfinx.fem.Function(Q)) - - interval_mesh = dolfinx.mesh.create_unit_interval(comm, 4) - with pytest.raises(NotImplementedError, match="interval"): - write_opencarp(filename, interval_mesh, scale=1.0) + io4dolfinx.write_cell_data( + filename, + dolfinx.fem.Function(Q), + 0.0, + io4dolfinx.FileMode.write, + None, + "opencarp", + on_input_mesh=True, + ) # Errors found when writing on rank 0 are raised on all ranks - write_opencarp(filename, mesh, scale=1.0) - with pytest.raises(ValueError, match="input node ordering"): - io4dolfinx.write_meshtags(filename, mesh, cell_tags, backend="opencarp") zero = dolfinx.fem.Function(fiber.function_space) with pytest.raises(ValueError, match="zero direction"): io4dolfinx.write_cell_data( filename, zero, 0.0, io4dolfinx.FileMode.write, None, "opencarp", on_input_mesh=True ) - with pytest.raises(ValueError, match="input cell ordering"): - io4dolfinx.write_cell_data( - filename, fiber, 0.0, io4dolfinx.FileMode.write, None, "opencarp" - ) def test_2d_surface_not_supported(shared_tmp_path): diff --git a/tests/test_write_data_input_order.py b/tests/test_write_data_input_order.py index beecde3..17d98a3 100644 --- a/tests/test_write_data_input_order.py +++ b/tests/test_write_data_input_order.py @@ -45,13 +45,11 @@ def test_data_on_input_mesh(captured, cell_type): mesh_data = create_original_mesh_data(mesh) point_data = captured["Point"] - assert point_data.on_input_mesh assert tuple(point_data.local_range) == tuple(mesh_data.local_geometry_pos) np.testing.assert_allclose(point_data.values, f(mesh_data.local_geometry.T).T, atol=1e-12) # P-1 interpolated into DG-0 is the mean of the vertex values for these cells cell_data = captured["Cell"] - assert cell_data.on_input_mesh assert tuple(cell_data.local_range) == tuple(mesh_data.local_topology_pos) x = np.vstack(MPI.COMM_WORLD.allgather(mesh_data.local_geometry)) midpoints = x[mesh_data.local_topology].mean(axis=1) From a04351349a286874fbaf6a0f4bf8ecad329f149e Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 2 Nov 2026 09:27:45 +0000 Subject: [PATCH 7/8] Fix openCARP marker files and share input-order redistribution - Record marker files in .markers.json. write_mesh removes the marker files of an earlier mesh, "facet" reads only the markers of the mesh with their values, and values giving the same file name are rejected. - Orient facets between two elements out of the element whose centre comes first in (x, y, z) order, independent of element ordering, and build the facet lookup once for all markers. - Check that .lon has a direction per element when reading or appending. - Reject complex cell data, and write 2D meshes with 3D coordinates in a plane of constant z. - Move the redistribution to input order into send_to_input_order and use it in create_original_mesh_data. Co-Authored-By: Claude Opus 5.5 --- docs/backends/opencarp.rst | 20 ++- src/io4dolfinx/backends/opencarp/backend.py | 173 +++++++++++++++---- src/io4dolfinx/checkpointing.py | 44 +---- src/io4dolfinx/comm_helpers.py | 37 +++- src/io4dolfinx/original_checkpoint.py | 157 ++--------------- tests/test_opencarp.py | 181 +++++++++++++++++++- 6 files changed, 391 insertions(+), 221 deletions(-) diff --git a/docs/backends/opencarp.rst b/docs/backends/opencarp.rst index f4bbb4e..e8c4ae3 100644 --- a/docs/backends/opencarp.rst +++ b/docs/backends/opencarp.rst @@ -15,12 +15,14 @@ documentation below. Rank 0 writes all files. openCARP expects coordinates in µ converting the mesh coordinates to µm must always be given as ``backend_args["scale"]``. With ``markers={"ENDO": (1, 2), "EPI": (2, 2)}``, a heart mesh gives ``heart.pts``, -``heart.elem``, ``heart.lon``, and ``heart_ENDO.vtx``, ``heart_ENDO.surf``, ``heart_EPI.vtx`` and -``heart_EPI.surf``. +``heart.elem``, ``heart.lon``, ``heart_ENDO.vtx``, ``heart_ENDO.surf``, ``heart_EPI.vtx`` and +``heart_EPI.surf``, and ``heart.markers.json`` with the name, value and dimension of each marker. +Writing the mesh again removes the ``.lon`` and marker files of the old mesh, and writing facet tags +again replaces the facet markers. - Triangle (``Tr``) and quadrilateral (``Qd``) meshes in 2D, and tetrahedral (``Tt``) and hexahedral (``Hx``) meshes in 3D are supported. 2D meshes are written with ``z = 0``, and read back with 2D - coordinates. + coordinates. 2D meshes with 3D coordinates in a plane of constant ``z`` are also written. - Elements get the orientation used by openCARP's ``mesher``: positive volume in 3D, counter-clockwise in 2D. - Cells without a cell tag get region ``0`` (set ``backend_args={"default_tag": ...}`` to change it). @@ -28,8 +30,10 @@ With ``markers={"ENDO": (1, 2), "EPI": (2, 2)}``, a heart mesh gives ``heart.pts ``z = 0``, and are read back with two components. Fields in a quadrature space cannot be written; generate them in a P-1 or DG-0 space. - Each facet marker value gets a ``.vtx`` file with its nodes and, for 3D meshes, a ``.surf`` file - with its facets (``Tr`` or ``Qd``), oriented out of the mesh. The files are named from - ``markers``, or by the value if it has no name. + with its facets (``Tr`` or ``Qd``), oriented out of the mesh. A facet between two elements + points out of the element whose centre comes first in lexicographic ``(x, y, z)`` order. The files + are named from ``markers``, or by the value if it has no name; values that would give the same + file name are rejected. Pass ``backend_args={"surf": False}`` to skip the ``.surf`` files. Reading takes the path of any of the mesh files. Coordinates are divided by ``scale``: @@ -43,8 +47,10 @@ Reading takes the path of any of the mesh files. Coordinates are divided by ``sc endo = io4dolfinx.read_meshtags("heart.elem", mesh, "ENDO", backend="opencarp", backend_args=args) f0 = io4dolfinx.read_cell_data("heart.elem", "fiber", mesh, backend="opencarp") -``"facet"`` reads all ``.surf`` files of the mesh, and a marker name reads its ``.surf`` file, or its -``.vtx`` file as vertex tags if there is no ``.surf`` (as for 2D meshes). +``"facet"`` reads the ``.surf`` files of the facet markers listed in ``heart.markers.json``, with +their values. A marker name reads its ``.surf`` file, or its ``.vtx`` file as vertex tags if there is +no ``.surf`` (as for 2D meshes). This also works for marker files written by other tools; their value +is taken from ``markers``, else from the name if it is an integer, else it is ``1``. .. automodule:: io4dolfinx.backends.opencarp.backend :members: diff --git a/src/io4dolfinx/backends/opencarp/backend.py b/src/io4dolfinx/backends/opencarp/backend.py index fb1b0ca..4a0d97d 100644 --- a/src/io4dolfinx/backends/opencarp/backend.py +++ b/src/io4dolfinx/backends/opencarp/backend.py @@ -18,6 +18,9 @@ and one node index per line, for stimuli and boundary conditions. - ``_.surf``: number of faces, then one ``Tr n0 n1 n2`` or ``Qd n0 n1 n2 n3`` line per face of a 3D mesh. +- ``.markers.json``: the name, value and dimension of each marker written by + io4dolfinx, so that reading gives the original values, and writing the mesh again + removes the marker files of the old mesh. openCARP does not read it. Each file is written by rank 0. @@ -73,7 +76,8 @@ - ``"default_tag"``: region tag of cells without a cell tag (default ``0``). - ``"markers"``: map from marker name to value, or to ``(value, dim)`` as used by cardiac-geometriesx. Names the ``.vtx`` and ``.surf`` files of a marker value, - which otherwise use the value itself. + which otherwise use the value itself. When reading, it takes precedence over the + values in ``.markers.json``. - ``"surf"``: write ``.surf`` files for facet tags of 3D meshes (default ``True``). - ``"vtx_domain"``: domain written to ``.vtx`` files (default ``"extra"``, meaning indices into the whole mesh, as written; ``"intra"`` means the myocardium submesh). @@ -84,6 +88,7 @@ """ import glob +import json from pathlib import Path from typing import Any, Callable, TypeVar @@ -212,20 +217,57 @@ def _marker_names(markers: dict[str, Any], dim: int) -> dict[int, str]: return names -def _marker_value(markers: dict[str, Any], name: str) -> int: +def _marker_index_file(filename: Path | str) -> Path: + """Get the path of the index of the markers of a mesh, ``.markers.json``.""" + return _mesh_file(filename, ".markers.json") + + +def _read_marker_index(filename: Path | str) -> dict[str, dict[str, Any]]: + """Read the index of the markers written for a mesh. + + Returns: + Map from marker name to its ``"value"``, ``"dim"`` and whether it has a + ``.surf`` file (``"surf"``). Empty if there is no index. + """ + path = _marker_index_file(filename) + if not path.exists(): + return {} + return json.loads(path.read_text()) + + +def _write_marker_index(filename: Path | str, index: dict[str, dict[str, Any]]): + """Write the index of the markers of a mesh, or remove it if there are none.""" + path = _marker_index_file(filename) + if index: + path.write_text(json.dumps(index, indent=2, sort_keys=True) + "\n") + else: + path.unlink(missing_ok=True) + + +def _remove_marker_files(filename: Path | str, names: list[str]): + """Remove the ``.vtx`` and ``.surf`` files of markers.""" + for name in names: + for suffix in (".vtx", ".surf"): + _marker_file(filename, name, suffix).unlink(missing_ok=True) + + +def _marker_value(markers: dict[str, Any], index: dict[str, dict[str, Any]], name: str) -> int: """Get the value of a marker read from ``_.vtx`` or ``.surf``. Args: markers: Map from name to value, or to ``(value, dim)`` + index: The index of the markers of the mesh, see :func:`_read_marker_index` name: Name of the marker in the filename Returns: - The value in ``markers``, else the name itself if it is an integer (as - written without markers), else ``1``. + The value in ``markers``, else in ``index``, else the name itself if it is an + integer, else ``1``. """ if name in markers: marker = markers[name] return int(marker) if isinstance(marker, (int, np.integer)) else int(marker[0]) + if name in index: + return int(index[name]["value"]) try: return int(name) except ValueError: @@ -295,6 +337,20 @@ def _read_count(path: Path) -> int: return int(f.readline().split()[0]) +def _count_rows(path: Path, header: int = 1) -> int: + """Count the non-empty lines after the header of a file.""" + with open(path) as f: + return sum(1 for line in f if line.strip()) - header + + +def _check_lon(lon: Path, num_elements: int): + """Check that ``.lon`` has a direction line per element.""" + if (num_rows := _count_rows(lon)) != num_elements: + raise ValueError( + f"{lon} has {num_rows} directions, but the mesh has {num_elements} elements." + ) + + def _read_rows( path: Path, local_range: tuple[int, int] | list[int], num_columns: int, header: int = 1 ) -> npt.NDArray[np.str_]: @@ -387,7 +443,7 @@ def _match_rows(reference: npt.NDArray[np.int64], query: npt.NDArray[np.int64]) Returns: Index into ``reference`` of each query row, ``-1`` if it has no match. If - several rows match, the last one is used. + several rows match, the first one is used. """ if len(query) == 0: return np.zeros(0, dtype=np.int64) @@ -395,7 +451,8 @@ def _match_rows(reference: npt.NDArray[np.int64], query: npt.NDArray[np.int64]) _, inverse = np.unique(rows, axis=0, return_inverse=True) inverse = inverse.reshape(-1) lookup = np.full(inverse.max() + 1, -1, dtype=np.int64) - lookup[inverse[: len(reference)]] = np.arange(len(reference)) + matched, first = np.unique(inverse[: len(reference)], return_index=True) + lookup[matched] = first return lookup[inverse[len(reference) :]] @@ -434,11 +491,11 @@ def _oriented_facets( cells: npt.NDArray[np.int64], facets: npt.NDArray[np.int64], ) -> npt.NDArray[np.int64]: - """Orient facets of a 3D mesh and sort them in a canonical order. + """Orient facets of a 3D mesh to point out of an element they are a face of. - Each facet gets the normal pointing out of an element it is a face of, starts with - its smallest node and the facets are sorted, so the result does not depend on how - the facets were distributed. + A facet between two elements points out of the element whose centre comes first + in lexicographic ``(x, y, z)`` order, so the result does not depend on the + ordering of the elements. Args: x: Node coordinates @@ -447,21 +504,27 @@ def _oriented_facets( facets: Nodes of each facet in cyclic order Returns: - The oriented facets + The oriented facets, in the order of ``facets`` """ local_faces = _faces[elem_type] faces = cells[:, local_faces].reshape(-1, local_faces.shape[1]) - match = _match_rows(faces, facets) + face_cell_centres = np.repeat(x[cells].mean(axis=1), len(local_faces), axis=0) + order = np.lexsort(face_cell_centres.T[::-1]) + match = _match_rows(faces[order], facets) if np.any(match < 0): raise ValueError(f"{np.count_nonzero(match < 0)} tagged facets are not faces of the mesh.") - cell_centres = x[cells[match // len(local_faces)]].mean(axis=1) + cell_centres = face_cell_centres[order][match] facets = facets.copy() outward = x[facets].mean(axis=1) - cell_centres inward = np.einsum("ij,ij->i", _facet_normals(x, facets), outward) < 0 facets[inward] = facets[inward][:, _mirror[_face_types[elem_type]]] + return facets + - # Rotating the nodes keeps the orientation +def _sort_facets(facets: npt.NDArray[np.int64]) -> npt.NDArray[np.int64]: + """Start each facet with its smallest node, keeping its orientation, and sort the + facets, so the result does not depend on how the facets were distributed.""" num_nodes = facets.shape[1] shift = np.argmin(facets, axis=1)[:, None] facets = np.take_along_axis(facets, (shift + np.arange(num_nodes)) % num_nodes, axis=1) @@ -585,11 +648,12 @@ def write_mesh( ): """Write a linear mesh to ``.pts`` and ``.elem``. - Supports triangle and quadrilateral meshes in 2D, and tetrahedral and hexahedral - meshes in 3D. Every cell gets the region tag ``backend_args["default_tag"]``, until - cell tags are written with :func:`write_meshtags`. Elements with the opposite - orientation of openCARP's ``mesher`` are mirrored. A ``.lon`` from an - earlier mesh is removed. + Supports triangle and quadrilateral meshes in 2D (also with 3D coordinates in a + plane of constant z), and tetrahedral and hexahedral meshes in 3D. Every cell gets + the region tag ``backend_args["default_tag"]``, until cell tags are written with + :func:`write_meshtags`. Elements with the opposite orientation of openCARP's + ``mesher`` are mirrored. The ``.lon`` and marker files of an earlier mesh + are removed. Args: filename: Base name of the files, or the path of one of them @@ -608,11 +672,6 @@ def write_mesh( ) elem_type = _element_types[mesh.cell_type] gdim = mesh.local_geometry.shape[1] - if gdim != _tdim[elem_type]: - raise NotImplementedError( - f"The openCARP backend needs {_tdim[elem_type]}D coordinates for a " - f"{mesh.cell_type} mesh, got {gdim}D." - ) if mode != FileMode.write: raise NotImplementedError("The openCARP backend can only write a mesh to new files.") scale = _get_scale(args) @@ -624,6 +683,11 @@ def write(): assert x is not None and cells is not None points = np.zeros((len(x), 3), dtype=x.dtype) points[:, :gdim] = x + if _tdim[elem_type] == 2 and np.ptp(points[:, 2]) > 0: + raise NotImplementedError( + f"The {mesh.cell_type} mesh does not lie in a plane of constant z. Surface " + "meshes in 3D are not supported." + ) elements = cells[:, _dolfinx_to_opencarp[elem_type]] orientation = _orientation(elem_type, points, elements) if np.any(orientation == 0): @@ -635,6 +699,8 @@ def write(): tags = np.full(len(elements), args["default_tag"], dtype=np.int32) _write_elem(_mesh_file(filename, ".elem"), elem_type, elements, tags) _mesh_file(filename, ".lon").unlink(missing_ok=True) + _remove_marker_files(filename, list(_read_marker_index(filename))) + _write_marker_index(filename, {}) _on_root(comm, write) @@ -708,9 +774,12 @@ def write_meshtags( the tagged entities. - Facet tags of a 3D mesh are also written to ``_.surf``, as triangles or quadrilaterals with outward normals (if ``backend_args["surf"]``). + A facet between two elements points out of the element whose centre comes + first in lexicographic ``(x, y, z)`` order. ```` is the name of the value in ``backend_args["markers"]``, or the - value itself. + value itself. The markers are recorded in ``.markers.json``, and replace + the markers of the same dimension written earlier. Args: filename: Base name of the files, or the path of one of them @@ -738,6 +807,24 @@ def write(): _write_elem(elem, elem_type, cells, tags) return + names = _marker_names(args["markers"], data.dim) + marker_values = np.unique(values) + marker_names = [names.get(int(value), str(value)) for value in marker_values] + if len(set(marker_names)) < len(marker_names): + raise ValueError( + f"Marker values {marker_values.tolist()} of dimension {data.dim} get the " + f"file names {marker_names}, which are not unique." + ) + index = _read_marker_index(filename) + replaced = [name for name, marker in index.items() if marker["dim"] == data.dim] + if clash := sorted(set(marker_names) & (set(index) - set(replaced))): + raise ValueError( + f"Markers {clash} of dimension {data.dim} have the same name as markers of " + f"another dimension in {_marker_index_file(filename)}." + ) + _remove_marker_files(filename, replaced) + index = {name: marker for name, marker in index.items() if name not in replaced} + write_surf = args["surf"] and _tdim[elem_type] == 3 and data.dim == 2 if write_surf: pts = _mesh_file(filename, ".pts") @@ -745,10 +832,9 @@ def write(): face_type = _face_types[elem_type] # DOLFINx orders the nodes of a quadrilateral facet as a tensor product facets = indices[:, _dolfinx_to_opencarp[face_type]] - names = _marker_names(args["markers"], data.dim) - for value in np.unique(values): + facets = _oriented_facets(x, elem_type, cells, facets) + for value, name in zip(marker_values, marker_names): marked = values == value - name = names.get(int(value), str(value)) nodes = np.unique(indices[marked]) _write_rows( _marker_file(filename, name, ".vtx"), @@ -757,13 +843,15 @@ def write(): "%d", ) if write_surf: - faces = _oriented_facets(x, elem_type, cells, facets[marked]) + faces = _sort_facets(facets[marked]) _write_rows( _marker_file(filename, name, ".surf"), [str(len(faces))], faces, _element_format(face_type, faces.shape[1]), ) + index[name] = {"value": int(value), "dim": int(data.dim), "surf": bool(write_surf)} + _write_marker_index(filename, index) _on_root(comm, write) @@ -780,11 +868,13 @@ def read_meshtags_data( filename: Base name of the files, or the path of one of them comm: MPI communicator used in storage name: ``"cell"`` for the region tags in ``.elem``, ``"facet"`` for the - faces of all ``_*.surf`` files, otherwise the name of a marker: - the faces of ``_.surf`` or, if there is none, the nodes of + faces of the ``.surf`` files of all facet markers in + ``.markers.json``, otherwise the name of a marker: the faces of + ``_.surf`` or, if there is none, the nodes of ``_.vtx``. backend_args: Arguments to backend. The value of a marker is taken from - ``"markers"``, else from its name if that is an integer, else it is ``1``. + ``"markers"``, else from ``.markers.json``, else from its name if + that is an integer, else it is ``1``. Returns: Internal data structure for the mesh tags read from file @@ -800,10 +890,18 @@ def read_meshtags_data( return MeshTagsData(name=name, values=tags, indices=cells, dim=_tdim[elem_type]) base = _basename(filename) + index = _read_marker_index(filename) if name == "facet": - files = sorted(Path(p) for p in glob.glob(f"{glob.escape(str(base))}_*.surf")) + files = [ + _marker_file(filename, marker, ".surf") + for marker, entry in sorted(index.items()) + if entry["dim"] == 2 and entry["surf"] + ] if len(files) == 0: - raise FileNotFoundError(f"Could not find any {base}_*.surf files.") + raise FileNotFoundError( + f"{_marker_index_file(filename)} lists no facet markers with .surf files. " + "Read the markers of other files by name." + ) dim = 2 elif (surf := _marker_file(filename, name, ".surf")).exists(): files, dim = [surf], 2 @@ -835,7 +933,7 @@ def read(): rows = _read_rows(path, (0, _read_count(path)), 1, header=2) entities = rows.astype(np.int64).reshape(-1, 1) all_indices.append(entities) - all_values.append(np.full(len(entities), _marker_value(args["markers"], marker))) + all_values.append(np.full(len(entities), _marker_value(args["markers"], index, marker))) return np.vstack(all_indices), np.concatenate(all_values).astype(np.int32) indices, values = np.zeros((0, 1), dtype=np.int64), np.zeros(0, dtype=np.int32) @@ -1003,6 +1101,8 @@ def read_cell_data( f"Could not find '{name}' in {lon}. Available: {list(_directions[:num_directions])}" ) + _on_root(comm, lambda: _check_lon(lon, _read_count(elem))) + def read(): local_range = compute_local_range(comm, _read_count(elem)) elem_type, cells, _ = _read_elem(elem, local_range) @@ -1062,6 +1162,8 @@ def write_data( raise NotImplementedError( "The openCARP backend can only write cell data (fibers and sheets)." ) + if np.iscomplexobj(array_data.values): + raise NotImplementedError("openCARP directions are real, got complex cell data.") num_cells, num_components = array_data.global_shape if num_components not in (2, 3): raise ValueError(f"Expected a 2D or 3D vector per cell, got {num_components} components.") @@ -1088,6 +1190,7 @@ def write(): check_file_exists(lon) if _read_count(lon) != 1: raise ValueError(f"{lon} already holds fibers and sheets.") + _check_lon(lon, num_cells) fibers = _read_rows(lon, (0, num_cells), 3).astype(np.float64) rows = np.hstack([fibers, rows]) _write_rows(lon, [str(rows.shape[1] // 3)], rows, _float_format) diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index 9c83558..bc287c4 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -22,10 +22,10 @@ from .backends import FileMode, ReadMode, get_backend from .comm_helpers import ( - exchange_to_owners, send_and_recv_cell_perm, send_dofmap_and_recv_values, send_dofs_and_recv_values, + send_to_input_order, ) from .readers import create_geometry_function_space from .structures import ArrayData, FunctionData, MeshTagsData @@ -81,33 +81,6 @@ def _with_mesh_name( return args -def _to_input_order( - comm: MPI.Intracomm, - input_indices: npt.NDArray[np.int64], - values: npt.NDArray, - num_global: int, -) -> tuple[npt.NDArray, tuple[int, int]]: - """Send values to the process holding their input index, so that each process - gets a contiguous chunk of the input ordering. - - Args: - comm: The MPI communicator - input_indices: Input index of each owned entity on this process - values: Value of each owned entity, of shape ``(len(input_indices), ...)`` - num_global: Global number of entities - - Returns: - The values of the entities in this process' chunk, in input order, and the range - of the chunk - """ - owners = index_owner(comm, input_indices, num_global) - (recv_indices, recv_values), *_ = exchange_to_owners(comm, owners, [input_indices, values]) - local_range = compute_local_range(comm, num_global) - data = np.empty_like(recv_values) - data[recv_indices - local_range[0]] = recv_values - return data, (local_range[0], local_range[1]) - - def write_attributes( filename: Path | str, comm: MPI.Comm, @@ -786,7 +759,7 @@ def write_point_data( V.mesh.geometry.input_global_indices[:num_dofs_local], dtype=np.int64 ) assert isinstance(comm, MPI.Intracomm) - data, local_range = _to_input_order(comm, input_nodes, data, data_shape[0]) + data, local_range = send_to_input_order(comm, input_nodes, data, data_shape[0]) ad = ArrayData( name=v_out.name, @@ -840,9 +813,7 @@ def write_cell_data( v_out.interpolate(u) comm = v_out.function_space.mesh.comm data_shape = (V.dofmap.index_map.size_global, V.dofmap.index_map_bs) - local_range = V.dofmap.index_map.local_range - num_dofs_local = V.dofmap.index_map.size_local - data = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[:num_dofs_local] + values = v_out.x.array.reshape(-1, V.dofmap.index_map_bs) if on_input_mesh: # DG-0 dofs of owned cells are owned, so send the value of each owned cell @@ -850,11 +821,12 @@ def write_cell_data( tdim = mesh.topology.dim num_owned_cells = mesh.topology.index_map(tdim).size_local input_cells = mesh.topology.original_cell_index[:num_owned_cells].astype(np.int64) - cell_values = v_out.x.array.reshape(-1, V.dofmap.index_map_bs)[ - V.dofmap.list[:num_owned_cells, 0] - ] + cell_values = values[V.dofmap.list[:num_owned_cells, 0]] assert isinstance(comm, MPI.Intracomm) - data, local_range = _to_input_order(comm, input_cells, cell_values, data_shape[0]) + data, local_range = send_to_input_order(comm, input_cells, cell_values, data_shape[0]) + else: + local_range = V.dofmap.index_map.local_range + data = values[: V.dofmap.index_map.size_local] ad = ArrayData( name=v_out.name, diff --git a/src/io4dolfinx/comm_helpers.py b/src/io4dolfinx/comm_helpers.py index 1036b9f..a36ac77 100644 --- a/src/io4dolfinx/comm_helpers.py +++ b/src/io4dolfinx/comm_helpers.py @@ -9,7 +9,12 @@ import numpy as np import numpy.typing as npt -from .utils import compute_insert_position, compute_local_range, valid_function_types +from .utils import ( + compute_insert_position, + compute_local_range, + index_owner, + valid_function_types, +) __all__ = [ "send_dofmap_and_recv_values", @@ -17,6 +22,7 @@ "send_dofs_and_recv_values", "neighbourhood_ranks", "exchange_to_owners", + "send_to_input_order", "numpy_to_mpi", "all_to_all", ] @@ -383,3 +389,32 @@ def exchange_to_owners( received.append(out) forward.Free() return received, sources, send_counts, recv_counts, insert_position + + +def send_to_input_order( + comm: MPI.Intracomm, + input_indices: npt.NDArray[np.int64], + values: npt.NDArray, + num_global: int, +) -> tuple[npt.NDArray, tuple[int, int]]: + """Send values to the process holding their input index, so that each process + gets a contiguous chunk of the input ordering. + + Args: + comm: The MPI communicator + input_indices: Input index of each owned entity on this process + values: Value of each owned entity, of shape ``(len(input_indices), ...)`` + num_global: Global number of entities + + Returns: + The values of the entities in this process' chunk, in input order, and the range + of the chunk + """ + input_indices = np.asarray(input_indices, dtype=np.int64) + owners = index_owner(comm, input_indices, num_global) + (recv_indices, recv_values), *_ = exchange_to_owners(comm, owners, [input_indices, values]) + local_range = compute_local_range(comm, num_global) + assert len(recv_indices) == local_range[1] - local_range[0] + data = np.empty_like(recv_values) + data[recv_indices - local_range[0]] = recv_values + return data, (int(local_range[0]), int(local_range[1])) diff --git a/src/io4dolfinx/original_checkpoint.py b/src/io4dolfinx/original_checkpoint.py index 39d1dc8..56a4b5d 100644 --- a/src/io4dolfinx/original_checkpoint.py +++ b/src/io4dolfinx/original_checkpoint.py @@ -17,7 +17,7 @@ import numpy as np from .backends import FileMode, get_backend -from .comm_helpers import all_to_all, numpy_to_mpi +from .comm_helpers import all_to_all, send_to_input_order from .structures import FunctionData, MeshData from .utils import ( compute_insert_position, @@ -35,155 +35,32 @@ def create_original_mesh_data(mesh: dolfinx.mesh.Mesh) -> MeshData: """ Store data locally on output process """ + assert isinstance(mesh.comm, MPI.Intracomm) - # 1. Send cell indices owned by current process to the process which owned its input - - # Get the input cell index for cells owned by this process + # 1. Send the topology of each owned cell, in original node indices, to the process + # owning its original cell index num_owned_cells = mesh.topology.index_map(mesh.topology.dim).size_local - original_cell_index = mesh.topology.original_cell_index[:num_owned_cells] - - # Compute owner of cells on this process based on the original cell index num_cells_global = mesh.topology.index_map(mesh.topology.dim).size_global - output_cell_owner = index_owner(mesh.comm, original_cell_index, num_cells_global) - local_cell_range = compute_local_range(mesh.comm, num_cells_global) - - # Compute outgoing edges from current process to outputting process - # Computes the number of cells sent to each process at the same time - cell_destinations, _send_cells_per_proc = np.unique(output_cell_owner, return_counts=True) - send_cells_per_proc = _send_cells_per_proc.astype(np.int32) - del _send_cells_per_proc - assert isinstance(mesh.comm, MPI.Intracomm) - cell_to_output_comm = mesh.comm.Create_dist_graph( - [mesh.comm.rank], - [len(cell_destinations)], - cell_destinations.tolist(), - reorder=False, - ) - cell_sources, cell_dests, _ = cell_to_output_comm.Get_dist_neighbors() - assert np.allclose(cell_dests, cell_destinations) - - # Compute number of recieving cells - recv_cells_per_proc = np.zeros_like(cell_sources, dtype=np.int32) - send_cells_per_proc = send_cells_per_proc.astype(np.int32) - all_to_all(cell_to_output_comm, send_cells_per_proc, recv_cells_per_proc) - - assert recv_cells_per_proc.sum() == local_cell_range[1] - local_cell_range[0] - # Pack and send cell indices (used for mapping topology dofmap later) - cell_insert_position = compute_insert_position( - output_cell_owner, cell_destinations, send_cells_per_proc - ) - send_cells = np.empty_like(cell_insert_position, dtype=np.int64) - send_cells[cell_insert_position] = original_cell_index - recv_cells = np.empty(recv_cells_per_proc.sum(), dtype=np.int64) - send_cells_msg = [send_cells, send_cells_per_proc, MPI.INT64_T] - recv_cells_msg = [recv_cells, recv_cells_per_proc, MPI.INT64_T] - cell_to_output_comm.Neighbor_alltoallv(send_cells_msg, recv_cells_msg) - del send_cells_msg, recv_cells_msg, send_cells - - # Map received cells to the local index - local_cell_index = recv_cells - local_cell_range[0] - - # 2. Create dofmap based on original geometry indices and re-order in the same order as original - # cell indices on output process - - # Get original node index for all nodes (including ghosts) and convert dofmap to these indices + original_cell_index = mesh.topology.original_cell_index[:num_owned_cells] original_node_index = mesh.geometry.input_global_indices - _, num_nodes_per_cell = fenicsx_compat.dofmap(mesh).shape local_geometry_dofmap = fenicsx_compat.dofmap(mesh)[:num_owned_cells, :] - global_geometry_dofmap = original_node_index[local_geometry_dofmap.reshape(-1)] - - # Unroll insert position for geometry dofmap - dofmap_insert_position = unroll_insert_position(cell_insert_position, num_nodes_per_cell) - - # Create and commmnicate connecitivity in original geometry indices - send_geometry_dofmap = np.empty_like(dofmap_insert_position, dtype=np.int64) - send_geometry_dofmap[dofmap_insert_position] = global_geometry_dofmap - del global_geometry_dofmap - send_sizes_dofmap = send_cells_per_proc * num_nodes_per_cell - recv_sizes_dofmap = recv_cells_per_proc * num_nodes_per_cell - recv_geometry_dofmap = np.empty(recv_sizes_dofmap.sum(), dtype=np.int64) - send_geometry_dofmap_msg = [send_geometry_dofmap, send_sizes_dofmap, MPI.INT64_T] - recv_geometry_dofmap_msg = [recv_geometry_dofmap, recv_sizes_dofmap, MPI.INT64_T] - cell_to_output_comm.Neighbor_alltoallv(send_geometry_dofmap_msg, recv_geometry_dofmap_msg) - del send_geometry_dofmap_msg, recv_geometry_dofmap_msg - - # Reshape dofmap and sort by original cell index - recv_dofmap = recv_geometry_dofmap.reshape(-1, num_nodes_per_cell) - sorted_recv_dofmap = np.empty_like(recv_dofmap) - sorted_recv_dofmap[local_cell_index] = recv_dofmap - - # 3. Move geometry coordinates to input process - # Compute outgoing edges from current process and create neighbourhood communicator - # Also create number of outgoing cells at the same time - num_owned_nodes = mesh.geometry.index_map().size_local - num_nodes_global = mesh.geometry.index_map().size_global - output_node_owner = index_owner( - mesh.comm, original_node_index[:num_owned_nodes], num_nodes_global - ) - - node_destinations, _send_nodes_per_proc = np.unique(output_node_owner, return_counts=True) - send_nodes_per_proc = _send_nodes_per_proc.astype(np.int32) - del _send_nodes_per_proc - - geometry_to_owner_comm = mesh.comm.Create_dist_graph( - [mesh.comm.rank], - [len(node_destinations)], - node_destinations.tolist(), - reorder=False, + global_geometry_dofmap = original_node_index[local_geometry_dofmap].astype(np.int64) + sorted_recv_dofmap, local_cell_range = send_to_input_order( + mesh.comm, original_cell_index, global_geometry_dofmap, num_cells_global ) - node_sources, node_dests, _ = geometry_to_owner_comm.Get_dist_neighbors() - assert np.allclose(node_dests, node_destinations) - - # Compute send node insert positions - send_nodes_position = compute_insert_position( - output_node_owner, node_destinations, send_nodes_per_proc - ) - unrolled_nodes_positiion = unroll_insert_position(send_nodes_position, 3) - - send_coordinates = np.empty_like(unrolled_nodes_positiion, dtype=mesh.geometry.x.dtype) - send_coordinates[unrolled_nodes_positiion] = mesh.geometry.x[:num_owned_nodes, :].reshape(-1) - - # Send and recieve geometry sizes - send_coordinate_sizes = (send_nodes_per_proc * 3).astype(np.int32) - recv_coordinate_sizes = np.zeros_like(node_sources, dtype=np.int32) - all_to_all(geometry_to_owner_comm, send_coordinate_sizes, recv_coordinate_sizes) - - # Send node coordinates - recv_coordinates = np.empty(recv_coordinate_sizes.sum(), dtype=mesh.geometry.x.dtype) - mpi_type = numpy_to_mpi[recv_coordinates.dtype.type] - send_coord_msg = [send_coordinates, send_coordinate_sizes, mpi_type] - recv_coord_msg = [recv_coordinates, recv_coordinate_sizes, mpi_type] - geometry_to_owner_comm.Neighbor_alltoallv(send_coord_msg, recv_coord_msg) - del send_coord_msg, recv_coord_msg - - # Send node ordering for reordering the coordinates on output process - send_nodes = np.empty(num_owned_nodes, dtype=np.int64) - send_nodes[send_nodes_position] = original_node_index[:num_owned_nodes] - - recv_indices = np.empty(recv_coordinate_sizes.sum() // 3, dtype=np.int64) - send_nodes_msg = [send_nodes, send_nodes_per_proc, MPI.INT64_T] - recv_nodes_msg = [recv_indices, recv_coordinate_sizes // 3, MPI.INT64_T] - geometry_to_owner_comm.Neighbor_alltoallv(send_nodes_msg, recv_nodes_msg) - - # Compute local ording of received nodes - local_node_range = compute_local_range(mesh.comm, num_nodes_global) - recv_indices -= local_node_range[0] - - # Sort geometry based on input index and strip to gdim + # 2. Send the coordinates of each owned node to the process owning its original index gdim = mesh.geometry.dim - recv_nodes = recv_coordinates.reshape(-1, 3) - _geometry = np.empty(recv_nodes.shape, dtype=mesh.geometry.x.dtype) - _geometry[recv_indices, :] = recv_nodes - geometry = _geometry[:, :gdim].copy() - del _geometry, recv_nodes - - assert local_node_range[1] - local_node_range[0] == geometry.shape[0] + num_owned_nodes = mesh.geometry.index_map().size_local + num_nodes_global = mesh.geometry.index_map().size_global + geometry, local_node_range = send_to_input_order( + mesh.comm, + original_node_index[:num_owned_nodes], + mesh.geometry.x[:num_owned_nodes, :gdim].copy(), + num_nodes_global, + ) cmap = fenicsx_compat.cmap(mesh) - cell_to_output_comm.Free() - geometry_to_owner_comm.Free() - # NOTE: Could in theory store partitioning information, but would not work nicely # as one would need to read this data rather than the xdmffile. # NOTE: Local geometry type hint skip is only required on DOLFINX<0.10 where diff --git a/tests/test_opencarp.py b/tests/test_opencarp.py index 35f8abb..3f94e81 100644 --- a/tests/test_opencarp.py +++ b/tests/test_opencarp.py @@ -232,8 +232,8 @@ def test_files_independent_of_partitioning(shared_tmp_path, cell_type): write_problem(serial, MPI.COMM_SELF, cell_type) files = sorted(p.name for p in serial.parent.iterdir()) assert files == sorted(p.name for p in parallel.parent.iterdir()) - # .pts, .elem, .lon and a .vtx (and in 3D a .surf) per marker - assert len(files) == (7 if dolfinx.mesh.cell_dim(cell_type) == 3 else 5) + # .pts, .elem, .lon, .markers.json and a .vtx (and in 3D a .surf) per marker + assert len(files) == (8 if dolfinx.mesh.cell_dim(cell_type) == 3 else 6) for name in files: assert (serial.parent / name).read_bytes() == (parallel.parent / name).read_bytes() @@ -388,3 +388,180 @@ def test_2d_surface_not_supported(shared_tmp_path): backend="opencarp", backend_args={"scale": 1.0}, ) + + +def read_facets(filename: Path) -> tuple[list[str], np.ndarray, np.ndarray]: + """Read the facet types, and their nodes and normals in µm, of a ``.surf`` file.""" + surf = read_lines(filename) + x = np.array(read_lines(filename.parent / "mesh.pts")[1:], dtype=np.float64) + faces = np.array([row[1:] for row in surf[1:]], dtype=np.int64) + q = x[faces] + if faces.shape[1] == 3: + normals = np.cross(q[:, 1] - q[:, 0], q[:, 2] - q[:, 0]) + else: + normals = np.cross(q[:, 2] - q[:, 0], q[:, 3] - q[:, 1]) + return [row[0] for row in surf[1:]], faces, normals + + +def test_facet_markers_of_mesh_only(shared_tmp_path): + """Reading the facets of a mesh skips the marker files of a mesh whose base name + starts with the same name, and gets the values written without ``markers``.""" + comm = MPI.COMM_WORLD + write_problem(shared_tmp_path / "mesh", comm) + mesh, _, facet_tags, _, _ = create_problem(comm) + other = {"scale": 1000.0, "markers": {"ENDO": (5, 2)}} + tags = dolfinx.mesh.meshtags(mesh, 2, facet_tags.indices, np.full_like(facet_tags.values, 5)) + write_files(shared_tmp_path / "mesh_lv", mesh, other, True, facet_tags=tags) + assert (shared_tmp_path / "mesh_lv_ENDO.surf").exists() + + filename = shared_tmp_path / "mesh.elem" + in_mesh = io4dolfinx.read_mesh( + filename, comm, backend="opencarp", backend_args={"scale": 1000.0} + ) + facets = io4dolfinx.read_meshtags(filename, in_mesh, "facet", backend="opencarp") + midpoints = dolfinx.mesh.compute_midpoints(in_mesh, 2, facets.indices) + np.testing.assert_allclose(midpoints[:, 0], np.where(facets.values == 1, 0.0, 1.0)) + num_owned = in_mesh.topology.index_map(2).size_local + assert ( + comm.allreduce(np.count_nonzero(facets.indices < num_owned)) + == 2 * meshes[CellType.tetrahedron][3] + ) + + +def test_rewrite_removes_old_markers(shared_tmp_path): + """Writing a mesh again removes the marker files of the old mesh, and writing + facet tags again replaces the facet markers.""" + comm = MPI.COMM_WORLD + filename = shared_tmp_path / "mesh" + write_problem(filename, comm) + assert (shared_tmp_path / "mesh_right.surf").exists() + + mesh, _, facet_tags, _, _ = create_problem(comm) + left = facet_tags.find(1) + left_tags = dolfinx.mesh.meshtags(mesh, 2, left, np.ones_like(left)) + write_files(filename, mesh, {"scale": 1000.0, "markers": markers(3)}, True) + assert not any(shared_tmp_path.glob("mesh_*")) + assert not (shared_tmp_path / "mesh.markers.json").exists() + + write_problem(filename, comm) + io4dolfinx.write_meshtags( + filename, + mesh, + left_tags, + backend="opencarp", + backend_args={"markers": markers(3)}, + on_input_mesh=True, + ) + assert sorted(p.name for p in shared_tmp_path.glob("mesh_*")) == [ + "mesh_left.surf", + "mesh_left.vtx", + ] + + +def test_marker_name_clash(shared_tmp_path): + """Marker values that would write the same files are rejected.""" + comm = MPI.COMM_WORLD + filename = shared_tmp_path / "mesh" + mesh, _, facet_tags, _, _ = create_problem(comm) + io4dolfinx.write_mesh_input_order( + filename, mesh, backend="opencarp", backend_args={"scale": 1.0} + ) + # Value 1 is named "2", the same name as the unnamed value 2 + with pytest.raises(ValueError, match="not unique"): + io4dolfinx.write_meshtags( + filename, + mesh, + facet_tags, + backend="opencarp", + backend_args={"markers": {"2": (1, 2)}}, + on_input_mesh=True, + ) + + io4dolfinx.write_meshtags(filename, mesh, facet_tags, backend="opencarp", on_input_mesh=True) + mesh.topology.create_connectivity(2, 0) + vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, facet_tags.find(1), 2, 0) + vertex_tags = dolfinx.mesh.meshtags(mesh, 0, vertices, np.ones_like(vertices)) + with pytest.raises(ValueError, match="another dimension"): + io4dolfinx.write_meshtags( + filename, mesh, vertex_tags, backend="opencarp", on_input_mesh=True + ) + + +@orderings +@pytest.mark.parametrize("cell_type", [CellType.tetrahedron, CellType.hexahedron]) +def test_interior_facet_orientation(shared_tmp_path, cell_type, on_input_mesh): + """A facet between two elements points out of the element whose centre has the + smallest x, whatever the ordering of the elements.""" + comm = MPI.COMM_WORLD + filename = shared_tmp_path / "mesh" + mesh, *_ = create_problem(comm, cell_type) + facets = dolfinx.mesh.locate_entities(mesh, 2, lambda x: np.isclose(x[0], 1 / 3)) + tags = dolfinx.mesh.meshtags(mesh, 2, facets, np.full_like(facets, 3)) + write_files(filename, mesh, {"scale": 1000.0}, on_input_mesh, facet_tags=tags) + if comm.rank == 0: + _, faces, normals = read_facets(shared_tmp_path / "mesh_3.surf") + assert len(faces) == meshes[cell_type][3] + assert np.all(normals[:, 0] > 0) + np.testing.assert_allclose(normals[:, 1:], 0, atol=1e-6) + + +def test_lon_rows_checked(shared_tmp_path): + """A ``.lon`` without a direction per element is rejected when read.""" + comm = MPI.COMM_WORLD + filename = shared_tmp_path / "mesh" + write_problem(filename, comm) + if comm.rank == 0: + lines = filename.with_suffix(".lon").read_text().splitlines() + filename.with_suffix(".lon").write_text("\n".join(lines[:-1]) + "\n") + comm.Barrier() + mesh = io4dolfinx.read_mesh( + filename.with_suffix(".elem"), comm, backend="opencarp", backend_args={"scale": 1000.0} + ) + with pytest.raises(ValueError, match="directions"): + io4dolfinx.read_cell_data(filename.with_suffix(".elem"), "fiber", mesh, backend="opencarp") + + +def test_complex_directions_rejected(shared_tmp_path): + from io4dolfinx.backends.opencarp import backend + from io4dolfinx.structures import ArrayData + + comm = MPI.COMM_WORLD + num_cells = 2 * comm.size + data = ArrayData( + name="f0", + values=np.ones((2, 3), dtype=np.complex128), + global_shape=(num_cells, 3), + local_range=(2 * comm.rank, 2 * comm.rank + 2), + type="Cell", + ) + with pytest.raises(NotImplementedError, match="complex"): + backend.write_data(shared_tmp_path / "mesh", data, comm, 0.0, io4dolfinx.FileMode.write, {}) + + +def test_flat_2d_mesh_with_3d_coordinates(shared_tmp_path): + """Triangles with 3D coordinates in a plane of constant z are written, and read back + as a 2D mesh. Triangles out of such a plane are rejected.""" + comm = MPI.COMM_WORLD + filename = shared_tmp_path / "mesh" + domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))) + cells = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.int64) + for z, ok in ((np.array([2.0, 2.0, 2.0, 2.0]), True), (np.array([0, 0, 0, 1.0]), False)): + x = np.column_stack([[0, 1, 0, 1], [0, 0, 1, 1], z]).astype(np.float64) + if comm.rank != 0: + cells, x = cells[:0], x[:0] + mesh = dolfinx.mesh.create_mesh(comm, cells=cells, x=x, e=domain) + if not ok: + with pytest.raises(NotImplementedError, match="constant z"): + io4dolfinx.write_mesh_input_order( + filename, mesh, backend="opencarp", backend_args={"scale": 1.0} + ) + continue + io4dolfinx.write_mesh_input_order( + filename, mesh, backend="opencarp", backend_args={"scale": 1.0} + ) + in_mesh = io4dolfinx.read_mesh( + filename.with_suffix(".elem"), comm, backend="opencarp", backend_args={"scale": 1.0} + ) + assert in_mesh.geometry.dim == 2 + area = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1 * ufl.dx(domain=in_mesh))) + assert np.isclose(comm.allreduce(area, op=MPI.SUM), 1.0) From d54b952c206a238cb18eb8e174240dc4937cbc8f Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 2 Nov 2026 09:33:05 +0000 Subject: [PATCH 8/8] Link to the openCARP examples repository in the docs Co-Authored-By: Claude Opus 5.5 --- docs/backends/opencarp.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/backends/opencarp.rst b/docs/backends/opencarp.rst index e8c4ae3..d6bc4fb 100644 --- a/docs/backends/opencarp.rst +++ b/docs/backends/opencarp.rst @@ -52,6 +52,10 @@ their values. A marker name reads its ``.surf`` file, or its ``.vtx`` file as ve no ``.surf`` (as for 2D meshes). This also works for marker files written by other tools; their value is taken from ``markers``, else from the name if it is an integer, else it is ``1``. +Examples that write a mesh with this backend, run openCARP on it and read the results back (a planar +wave in a square, and the Niederer et al. (2011) benchmark) are in +`opencarp-io4dolfinx-example `_. + .. automodule:: io4dolfinx.backends.opencarp.backend :members: :exclude-members: read_attributes, read_timestamps, write_attributes, read_dofmap, read_dofs, read_cell_perms, write_function, read_legacy_mesh, snapshot_checkpoint, read_hdf5_array, read_point_data