Add openCARP backend by finsberg · Pull Request #93 · scientificcomputing/io4dolfinx

finsberg · 2026-10-28T06:20:49Z

Summary

Adds an opencarp backend that writes and reads linear meshes in the ASCII mesh format of the cardiac simulator openCARP. It supports triangles (Tr) and quadrilaterals (Qd) in 2D, and tetrahedra (Tt) and hexahedra (Hx) in 3D.

File	Content
`<base>.pts`	Node coordinates (µm)
`<base>.elem`	Elements with region tags (from cell tags)
`<base>.lon`	Fiber, or fiber + sheet, direction per element
`<base>_<marker>.vtx`	Nodes of each facet or vertex marker
`<base>_<marker>.surf`	Facets of each marker of a 3D mesh (`Tr` or `Qd`)
`<base>.markers.json`	Name, value and dimension of each marker (io4dolfinx only, openCARP ignores it)

There is no backend-specific function. Files are written and read with the generic API and backend="opencarp":

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)

mesh = io4dolfinx.read_mesh("heart.elem", MPI.COMM_WORLD, backend="opencarp", backend_args=args)
facets = io4dolfinx.read_meshtags("heart.elem", mesh, "facet", backend="opencarp")
f0 = io4dolfinx.read_cell_data("heart.elem", "fiber", mesh, backend="opencarp")

Changes to existing code

write_cell_data and write_point_data get an on_input_mesh: bool = False option. With it, values are written in the cell or node ordering of the input mesh, matching a mesh written with write_mesh_input_order. The default leaves the other backends unchanged.
The redistribution to input order is a new helper, comm_helpers.send_to_input_order, built on exchange_to_owners. original_checkpoint.create_original_mesh_data now uses it too, which replaces about 130 lines of hand-written neighbourhood communication. Please check this part: it is the only change to code that predates the backend. The original-checkpoint tests pass.
write_cell_data rejects functions in a quadrature space with a clear error, since dolfinx cannot interpolate them into DG-0.

Design

Ordering. With write_mesh_input_order and on_input_mesh=True, the files use the input ordering and are byte-identical on 1–4 ranks. With write_mesh and the default on_input_mesh=False, they use the DOLFINx ordering. Rank 0 gathers and writes every file.
Units. openCARP uses µm, so scale (the factor from mesh coordinates to µm) is required for writing and reading. There is no default.
Orientation.
- Elements get the orientation of openCARP's mesher: positive volume in 3D, counter-clockwise in 2D.
- .surf facets point out of the mesh. A facet between two elements points out of the element whose centre comes first in (x, y, z) order, so the result doesn't depend on element ordering.
Markers. Marker files are listed in <base>.markers.json. As a result:
- writing the mesh again removes the old marker files;
- writing tags again replaces the markers of that dimension;
- read_meshtags(..., "facet") reads only this mesh's markers, with their values, so it skips the files of a mesh like heart_lv that sits next to heart;
- values that would give the same file name are rejected.
A single marker can also be read by name, including marker files written by other tools.
Fibers are interpolated to DG-0 and normalised. 2D directions get z = 0. Complex data is rejected, and .lon row counts are checked against .elem on read.
.vtx header is extra by default. openCARP 19 warns that intra is wrong for indices into the whole mesh.
Out of scope: binary formats, mixed element types, higher-order geometry and .igb simulation output.

Testing

tests/test_opencarp.py and tests/test_write_data_input_order.py don't need openCARP. They run serially and under MPI (2, 3 and 4 ranks), and cover:

file content for all four element types, in both orderings: headers, counts, positive orientation, unit fibers, .vtx nodes, outward .surf normals;
byte-identical output between COMM_WORLD and COMM_SELF writes;
write → read round trips of mesh, region tags, facet and vertex markers, and fibers and sheets;
regression tests for marker files of other meshes, stale marker files, marker name clashes, interior facet orientation, .lon length, complex data, and flat 2D meshes with 3D coordinates;
errors that must be raised on every rank.

The full suite passes, except 4 tests in test_legacy_readers.py that fail the same way on the unchanged code. They come from the local legacy data, not this PR. ruff and mypy pass.

openCARP 19.0.0 (conda-forge) was run on files written by the backend, using the examples in finsberg/opencarp-io4dolfinx-example:

Unit square: a planar wave from a .vtx stimulus. The right-hand region has a quarter of the conductivity and conducts at about half the speed (0.24 vs 0.57 m/s), which shows openCARP read the region tags.
Niederer et al. (2011) benchmark at 0.2 mm: all nodes activate, and P8 activates at 54 ms. Conduction is 0.56 m/s along the fibers and 0.14 m/s across, which shows the .lon was read.

Results are the same serially and on 3 ranks.

Docs

New page docs/backends/opencarp.rst (linking to the examples repository), with entries in _toc.yml and the README. The module docstring documents the write recipe and the backend arguments.

🤖 Generated with Claude Code

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 <noreply@anthropic.com>

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 <noreply@anthropic.com>

…kend 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 <noreply@anthropic.com>

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

jorgensd · 2026-10-30T09:51:13Z

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)},
)

I'm not that keen on this new API entry point.
As we discussed, as long as the flag write_on_input_mesh is passed consistently in write_mesh and write_meshtags the outputted mesh should work with opencarp.

Of course having the new write_on_input_mesh flag in write_meshtags is fine, as it just add functionality.

Furthermore, I think this function belongs as a helper function somewhere else, rather than being part of the actual backend, as it is very specific..:P

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 <noreply@anthropic.com>

- Record marker files in <base>.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 <noreply@anthropic.com>

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

jorgensd · 2026-11-02T11:08:11Z

+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`),


Why do we need both?:P

jorgensd · 2026-11-02T11:09:16Z

+    geometry, local_node_range = send_to_input_order(
+        mesh.comm,
+        original_node_index[:num_owned_nodes],
+        mesh.geometry.x[:num_owned_nodes, :gdim].copy(),


Nice cleanup

jorgensd · 2026-11-02T11:10:32Z

+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)


These look like general MPI utils and should be moved to such a file?

jorgensd · 2026-11-02T11:10:53Z

+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


Should these be put in an ascii_helper somewhere central?

jorgensd · 2026-11-02T11:11:05Z

+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)


Ascii helper?

jorgensd · 2026-11-02T11:11:18Z

+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)
+


Ascii helper!

jorgensd · 2026-11-02T11:11:38Z

+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) :]]


General numpy helper?

jorgensd · 2026-11-02T11:12:50Z

+    return np.cross(p[:, 2] - p[:, 0], p[:, 3] - p[:, 1])
+
+
+def _oriented_facets(


Maybe in general, we should start putting all helper functions like this in a helper file, so the backend file stays very clean and uniform across io4dolfinx? As any backend is a folder/module anyhow.

finsberg and others added 5 commits October 27, 2026 20:10

Fix mypy

f473e09

Fix return type of orientation helper in openCARP tests

682c777

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

finsberg requested a review from jorgensd October 30, 2026 09:18

finsberg and others added 3 commits October 30, 2026 14:15

Link to the openCARP examples repository in the docs

d54b952

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

jorgensd reviewed Nov 2, 2026

View reviewed changes

Navigation Menu

Repository navigation

Add openCARP backend - #93

Add openCARP backend#93
finsberg wants to merge 8 commits into
mainfrom
finsberg/opencarp-backend

finsberg commented Oct 28, 2026 •

edited

Loading

Uh oh!

jorgensd commented Oct 30, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

jorgensd Nov 2, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

		return np.cross(p[:, 2] - p[:, 0], p[:, 3] - p[:, 1])


		def _oriented_facets(

Conversation

finsberg commented Oct 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Changes to existing code

Design

Testing

Docs

Uh oh!

jorgensd commented Oct 30, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

finsberg commented Oct 28, 2026 •

edited

Loading