diff --git a/README.md b/README.md index d270878..0a6ad8d 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 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/_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..d6bc4fb --- /dev/null +++ b/docs/backends/opencarp.rst @@ -0,0 +1,61 @@ +.. _opencarpbackend: + +openCARP-backend +---------------- + +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. + +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"]``. + +With ``markers={"ENDO": (1, 2), "EPI": (2, 2)}``, a heart mesh gives ``heart.pts``, +``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. 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). +- 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. 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``: + +.. 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 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``. + +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 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..05458f8 --- /dev/null +++ b/src/io4dolfinx/backends/opencarp/__init__.py @@ -0,0 +1,3 @@ +from . import backend + +__all__ = ["backend"] diff --git a/src/io4dolfinx/backends/opencarp/backend.py b/src/io4dolfinx/backends/opencarp/backend.py new file mode 100644 index 0000000..4a0d97d --- /dev/null +++ b/src/io4dolfinx/backends/opencarp/backend.py @@ -0,0 +1,1198 @@ +""" +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. 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 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. + +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 +(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, + 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. 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). + +SPDX License identifier: MIT + +Copyright: Henrik Finsberg and Simula Research Laboratory +""" + +import glob +import json +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") + +# 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" + +# 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_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 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: + 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 _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: + """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 _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_]: + """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 _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[str, npt.NDArray[np.int64], npt.NDArray[np.int32]]: + """Read the elements of an ``.elem`` file. + + Args: + path: The ``.elem`` file + local_range: Range of elements to read, all if ``None`` + + Returns: + Tuple containing: + - 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, num_nodes + 2) + if other := sorted(set(rows[:, 0]) - {elem_type}): + raise NotImplementedError( + f"Meshes with several element types are not supported, {path} has " + f"{[elem_type, *other]}." + ) + 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 elem_type, 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 _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: + """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 first 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) + matched, first = np.unique(inverse[: len(reference)], return_index=True) + lookup[matched] = first + return lookup[inverse[len(reference) :]] + + +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. + + 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 a 3D mesh to point out of an element they are a face of. + + 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 + 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 facets, in the order of ``facets`` + """ + local_faces = _faces[elem_type] + faces = cells[:, local_faces].reshape(-1, local_faces.shape[1]) + 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 = 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 + + +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) + return facets[np.lexsort(facets.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 mesh to ``.pts`` and ``.elem``. + + 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 + comm: MPI communicator used in storage + 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 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 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 + 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): + 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) + _remove_marker_files(filename, list(_read_marker_index(filename))) + _write_marker_index(filename, {}) + + _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 mesh from ``.pts`` and ``.elem``. + + 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 + 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) + 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[:, np.argsort(_dolfinx_to_opencarp[elem_type])], + cell_type=_cell_types[elem_type], + 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. + - 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"]``). + 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. 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 + comm: MPI communicator used in storage + 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) + 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) + 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, 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") + 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]] + facets = _oriented_facets(x, elem_type, cells, facets) + for value, name in zip(marker_values, marker_names): + marked = values == value + nodes = np.unique(indices[marked]) + _write_rows( + _marker_file(filename, name, ".vtx"), + [str(len(nodes)), args["vtx_domain"]], + nodes, + "%d", + ) + if write_surf: + 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) + + +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 + 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 ``.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 + """ + args = get_default_backend_args(backend_args) + if name == "cell": + elem = _mesh_file(filename, ".elem") + check_file_exists(elem) + 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) + index = _read_marker_index(filename) + if name == "facet": + 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"{_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 + 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. 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) + 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"], 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) + result = _on_root(comm, read) + if result is not None: + indices, values = result + 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. A 2D mesh gets the x- and y-components. + """ + 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])}" + ) + + _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) + 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) + 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( + 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. + + 2D directions get ``z = 0``. + + 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 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 + backend_args: The backend arguments + """ + if array_data.type != "Cell": + 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.") + 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}.") + 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 = vectors / norms[:, None] + + if mode == FileMode.append: + 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) + + _on_root(comm, write) diff --git a/src/io4dolfinx/checkpointing.py b/src/io4dolfinx/checkpointing.py index 1822873..bc287c4 100644 --- a/src/io4dolfinx/checkpointing.py +++ b/src/io4dolfinx/checkpointing.py @@ -25,6 +25,7 @@ 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 @@ -724,6 +725,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. @@ -737,6 +739,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}") @@ -748,8 +752,21 @@ 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 = send_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", ) logger.debug( f"Using {backend} backend with arguments {backend_args} and mode {mode} to write point data" @@ -768,6 +785,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,20 +798,42 @@ 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}") + 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) 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 + mesh = V.mesh + 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 = values[V.dofmap.list[:num_owned_cells, 0]] + assert isinstance(comm, MPI.Intracomm) + 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, 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", ) logger.debug( f"Using {backend} backend with arguments {backend_args} and mode {mode} to write cell data" 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 new file mode 100644 index 0000000..3f94e81 --- /dev/null +++ b/tests/test_opencarp.py @@ -0,0 +1,567 @@ +from pathlib import Path + +from mpi4py import MPI + +import basix.ufl +import dolfinx +import numpy as np +import pytest +import ufl + +import io4dolfinx + +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) +orderings = pytest.mark.parametrize( + "on_input_mesh", [True, False], ids=["input_ordering", "dolfinx_ordering"] +) + + +def markers(tdim: int): + return {"left": (1, tdim - 1), "right": (2, tdim - 1)} + + +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.ones_like(x[0]), 3 + x[2], x[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, 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) + 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, (gdim,))) + fiber = dolfinx.fem.Function(V, name="f0") + fiber.interpolate(lambda x: fiber_expr(x)[:gdim]) + sheet = dolfinx.fem.Function(V, name="s0") + sheet.interpolate(lambda x: sheet_expr(x)[:gdim]) + 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, + on_input_mesh: bool = True, + backend_args: dict | None = None, +): + mesh, cell_tags, facet_tags, fiber, sheet = create_problem(comm, cell_type) + 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]]: + return [line.split() for line in path.read_text().splitlines()] + + +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] + 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 + + +@orderings +@cell_types +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, on_input_mesh) + 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 == num_nodes + x = np.array(pts[1:], dtype=np.float64) + 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 == 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 oriented as by openCARP's mesher, filling the cube or square + p = x[cells] + 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) + 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 == (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"mesh_{name}.vtx") + assert vtx[1] == ["extra"] + nodes = np.array(vtx[2:], dtype=np.int64).reshape(-1) + assert int(vtx[0][0]) == len(nodes) == num_marked_nodes + assert np.allclose(x[nodes, 0], plane) + + 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): + """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: + parallel.parent.mkdir() + serial.parent.mkdir() + MPI.COMM_WORLD.Barrier() + write_problem(parallel, MPI.COMM_WORLD, cell_type) + if MPI.COMM_WORLD.rank == 0: + 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, .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() + + +@orderings +@cell_types +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] + 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) + 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) + + 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)) == 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 + ) + 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, 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") + 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, 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 / "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(3)} + ) + 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 / "mesh" + mesh, _, facet_tags, fiber, _ = create_problem(MPI.COMM_WORLD) + 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:]} + 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 / "mesh" + mesh, _, _, fiber, _ = create_problem(comm) + + with pytest.raises(ValueError, match="scale"): + io4dolfinx.write_mesh_input_order(filename, mesh, backend="opencarp") + + 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"): + 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 + 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 + ) + + +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}, + ) + + +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) diff --git a/tests/test_write_data_input_order.py b/tests/test_write_data_input_order.py new file mode 100644 index 0000000..17d98a3 --- /dev/null +++ b/tests/test_write_data_input_order.py @@ -0,0 +1,56 @@ +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 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 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)