From 693760c7395ff211cfb045e4a475b8fdd882afd9 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 13:31:04 +1100
Subject: [PATCH 01/13] Test interpolation within one mesh on a rank that owns
 nothing (#902)

PETSc's DMPlexComputeInterpolatorNested skips its preallocation pass on a
rank with no local rows and no local columns but fills from that rank's
cells all the same. A shared point belongs to the highest rank holding it,
so a rank whose cells sit wholly inside higher ranks' regions owns nothing,
and its entries reach owners that never preallocated them: the owner raises
"New nonzero caused a malloc" inside MatAssemblyEnd, and the rank that did
not raise waits in the assembly's collective. On kaiju this stopped custom-P
multigrid at np 32 and 52 (Mesh._get_coords_for_basis on the 124-cell
coarsest level of the faulted annulus).

The test forces the case with a shell partition: rank 0 keeps one interior
cell of a 4 x 4 grid and owns none of its vertices and edges. It pins the
degree 1 and 2 node coordinates and the projection between degrees
(project_to_degree, which interpolates between a sub-DM and a clone of one
plex the same way). On development at np 2 the degree-2 coordinates raise
"New nonzero at (88,32)" on rank 1 and the projection fails after it; a
thread timeout turns a hang into a failure.

Underworld development team with AI support from Claude Code
---
 ...erpolation_rank_owning_nothing_parallel.py | 146 ++++++++++++++++++
 1 file changed, 146 insertions(+)
 create mode 100644 tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py

diff --git a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
new file mode 100644
index 000000000..608951139
--- /dev/null
+++ b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
@@ -0,0 +1,146 @@
+"""Interpolation within one mesh on a rank that holds cells but owns none of their DOFs.
+
+PETSc's ``DMPlexComputeInterpolatorNested`` builds the matrix between two DMs on
+one plex (a coordinate DM and a basis DM, a variable's sub-DM and a scratch
+clone). It skips its preallocation pass on a rank with no local rows and no
+local columns, but runs its fill pass over that rank's cells all the same. A
+shared point belongs to the highest rank that holds it, so a rank whose cells
+sit wholly inside higher ranks' regions owns nothing. Its entries then reach
+owners that never preallocated them: the owner raises "New nonzero caused a
+malloc" inside ``MatAssemblyEnd``, and the rank that did not raise waits in the
+assembly's own collective for ever. Found on kaiju at np 32 and 52 on the
+124-cell coarsest level of the faulted annulus, where custom-P multigrid reads
+each level's node coordinates (#902), and on a 3-D notch mesh at np 16.
+
+A shell partition gives rank 0 one interior cell of a 4 x 4 grid of triangles
+and the other ranks the rest, so rank 0 owns none of its cell's vertices and
+edges. The test checks that case is set up, then that the node coordinates and
+the projection between degrees come back right on every rank. Before the fix
+the degree-2 coordinates hung, so each test carries a thread timeout, which
+fails the run instead of stalling it.
+
+    mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
+"""
+import numpy as np
+import pytest
+from mpi4py import MPI
+from petsc4py import PETSc
+
+import underworld3 as uw
+from underworld3.function import project_to_degree
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, pytest.mark.tier_b,
+              pytest.mark.timeout(120, method="thread")]
+
+N = 4
+
+
+def _owned_nothing_mesh():
+    """The 4 x 4 grid of triangles, all on rank 0, partitioned so rank 0 keeps one
+    interior cell (every edge shared with a cell of another rank)."""
+    comm = PETSc.COMM_WORLD
+    rank, size = comm.getRank(), comm.getSize()
+    if rank == 0:
+        xs = np.linspace(0.0, 1.0, N + 1)
+        coords = np.array([[x, y] for y in xs for x in xs])
+        cells = []
+        for j in range(N):
+            for i in range(N):
+                a = j * (N + 1) + i
+                cells += [[a, a + 1, a + N + 2], [a, a + N + 2, a + N + 1]]
+        cells = np.asarray(cells, dtype=PETSc.IntType)
+    else:
+        coords = np.zeros((0, 2))
+        cells = np.zeros((0, 3), dtype=PETSc.IntType)
+    dm = PETSc.DMPlex().createFromCellList(2, cells, coords, interpolate=True, comm=comm)
+    part = dm.getPartitioner()
+    part.setType(PETSc.Partitioner.Type.SHELL)
+    if rank == 0:
+        cS, cE = dm.getHeightStratum(0)
+        interior = [c for c in range(cS, cE)
+                    if all(dm.getSupportSize(e) == 2 for e in dm.getCone(c))]
+        c0 = interior[len(interior) // 2]
+        rest = np.array_split(np.asarray([c for c in range(cS, cE) if c != c0]), size - 1)
+        part.setShellPartition(size, [1] + [len(r) for r in rest],
+                               [c0] + [int(c) for r in rest for c in r])
+    else:
+        part.setShellPartition(size, [0] * size, [])
+    dm.distribute()
+    return uw.discretisation.Mesh(dm, degree=1, qdegree=2)
+
+
+@pytest.fixture(scope="module")
+def mesh():
+    return _owned_nothing_mesh()
+
+
+def _all(flag):
+    return uw.mpi.comm.allreduce(bool(flag), op=MPI.LAND)
+
+
+def _cell_vertex_coords(mesh):
+    """Each local cell's vertex coordinates, read from the coordinate section."""
+    dm = mesh.dm
+    vS, vE = dm.getDepthStratum(0)
+    sec = dm.getCoordinateSection()
+    x = dm.getCoordinatesLocal().array
+    cS, cE = dm.getHeightStratum(0)
+    out = []
+    for c in range(cS, cE):
+        pts = dm.getTransitiveClosure(c)[0]
+        out.append(np.array([x[sec.getOffset(v):sec.getOffset(v) + 2]
+                             for v in pts if vS <= v < vE]))
+    return out
+
+
+def _sorted_rows(a):
+    a = np.round(np.asarray(a), 12)
+    return a[np.lexsort(a.T[::-1])]
+
+
+def test_rank_zero_holds_a_cell_and_owns_nothing(mesh):
+    """The case under test: rank 0 holds one cell and owns no vertex or edge."""
+    cS, cE = mesh.dm.getHeightStratum(0)
+    owned = 0
+    for d in (0, 1):
+        pS, pE = mesh.dm.getDepthStratum(d)
+        sec = mesh.dm.getPointSF()
+        leaves = set(np.asarray(sec.getGraph()[1]).tolist())
+        owned += sum(1 for p in range(pS, pE) if p not in leaves)
+    if uw.mpi.rank == 0:
+        assert cE - cS == 1 and owned == 0, (cE - cS, owned)
+    assert _all(True)
+
+
+@pytest.mark.parametrize("degree", [1, 2])
+def test_node_coordinates_on_every_rank(mesh, degree):
+    """Each cell's nodes are its vertices (degree 1), and its vertices and edge
+    midpoints (degree 2)."""
+    coords = mesh._get_coords_for_basis(degree, True)
+    rows = mesh._cell_node_indices(degree, True)
+    ok = True
+    for k, verts in enumerate(_cell_vertex_coords(mesh)):
+        expected = verts if degree == 1 else np.vstack(
+            [verts, 0.5 * (verts[[0, 1, 2]] + verts[[1, 2, 0]])])
+        got = coords[rows[k]]
+        ok &= np.allclose(_sorted_rows(got), _sorted_rows(expected), atol=1e-12)
+    assert _all(ok)
+
+
+def test_projection_between_degrees_on_every_rank(mesh):
+    """P1 to P2 of a linear field and P2 to P1 of a quadratic one are exact at the
+    nodes."""
+    s1 = uw.discretisation.MeshVariable("s1_1107", mesh, 1, degree=1)
+    s2 = uw.discretisation.MeshVariable("s2_1107", mesh, 1, degree=2)
+    x1 = np.asarray(s1.coords)
+    x2 = np.asarray(s2.coords)
+    s1.array[...] = (x1[:, 0] + 2.0 * x1[:, 1]).reshape(-1, 1, 1)
+    s2.array[...] = (x2[:, 0] ** 2).reshape(-1, 1, 1)
+
+    up = project_to_degree(s1, target_degree=2)
+    down = project_to_degree(s2, target_degree=1)
+    n2 = mesh._get_coords_for_basis(2, True)
+    n1 = mesh._get_coords_for_basis(1, True)
+    ok = (np.allclose(up[:, 0], n2[:, 0] + 2.0 * n2[:, 1], atol=1e-12)
+          and np.allclose(down[:, 0], n1[:, 0] ** 2, atol=1e-12))
+    assert _all(ok)

From fbcd2ef91cd8d0a3ac275fcd30d24781f0062fac Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 13:45:46 +1100
Subject: [PATCH 02/13] Test that a transfer-builder failure on one rank is
 retried on every rank (#886)

build_transfers retries a failed barycentric build with the RBF builder and
degrades to the default preconditioner after that, each rank on its own
verdict: a rank whose builder raised retried alone while the others went on
into the transfer's MPIAIJ assembly. The test makes the barycentric builder
raise on rank 1 only, on a single-field fmg Poisson at np 2. On the current
branch it times out with rank 0 inside _assemble_node_transfer; the contract
is the RBF rescue on every rank, multigrid installed, the solve converged.

Underworld development team with AI support from Claude Code
---
 ..._mg_builder_failure_collective_parallel.py | 52 +++++++++++++++++++
 1 file changed, 52 insertions(+)
 create mode 100644 tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py

diff --git a/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
new file mode 100644
index 000000000..2e87f1a87
--- /dev/null
+++ b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
@@ -0,0 +1,52 @@
+"""A transfer builder that fails on one rank is retried on every rank (#886).
+
+``build_transfers`` retries a failed barycentric build with the RBF builder, and
+degrades to the solver's default preconditioner when that fails too. The verdict
+was each rank's own: a rank whose builder raised retried alone while the others
+went on into the transfer's assembly, and the ranks met in different collectives
+(on kaiju at np 32 and 52, after #902's assembly failure, the run hung).
+
+The barycentric builder is rank-local work, so its verdict can be agreed right
+after it: a failure on any rank is raised on every rank, and the retry is taken
+together. Here the builder raises on rank 1 only. Contract: every rank records
+the RBF rescue, geometric multigrid is installed, and the solve converges. A
+thread timeout turns the old hang into a failure.
+
+    mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
+"""
+import pytest
+
+import underworld3 as uw
+from underworld3.utilities import custom_mg
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, pytest.mark.tier_b,
+              pytest.mark.timeout(300, method="thread")]
+
+
+def test_a_builder_failure_on_one_rank_is_retried_on_every_rank(monkeypatch):
+    barycentric = custom_mg._BUILDERS["barycentric"]
+
+    def fails_on_rank_one(coarse_coords, fine_coords):
+        if uw.mpi.rank == 1:
+            raise RuntimeError("synthetic barycentric failure on rank 1 (test probe)")
+        return barycentric(coarse_coords, fine_coords, repair_orphans=False)
+
+    monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", fails_on_rank_one)
+    box = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, refinement=2,
+        regular=True, qdegree=3)
+    T = uw.discretisation.MeshVariable("T1108", box, 1, degree=2)
+    p = uw.systems.Poisson(box, T)
+    p.constitutive_model = uw.constitutive_models.DiffusionModel
+    p.constitutive_model.Parameters.diffusivity = 1.0
+    p.add_dirichlet_bc(0.0, "Bottom")
+    p.add_dirichlet_bc(1.0, "Top")
+    p.preconditioner = "fmg"
+    with pytest.warns(UserWarning, match="retrying with the 'rbf' builder"):
+        p.solve()
+
+    assert p.snes.getConvergedReason() > 0
+    assert p.snes.getKSP().getPC().getType() == "mg"
+    rec = p.pc_fallbacks["custom_mg.transfer_builder"]
+    assert rec["reason"] == "build_failed" and "rank 1" in rec["detail"], rec
+    assert "custom_mg.build" not in p.pc_fallbacks

From 37ac408a99d330c381502d6f005e94ce27438b61 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 13:46:07 +1100
Subject: [PATCH 03/13] Interpolate within one mesh cell by cell, not through
 PETSc's matrix (#902)

Mesh._get_coords_for_basis and function.project_to_degree interpolated
between two DMs on one plex with createInterpolation, which takes PETSc's
DMPlexComputeInterpolatorNested. That function skips its preallocation on
a rank with no local rows and no local columns but fills from that rank's
cells, so on a rank that holds cells and owns none of their nodes the
owners raise inside MatAssemblyEnd and the other ranks wait in its
collective. custom-P multigrid reads every level's node coordinates, and
this stopped kaiju's faulted annulus at np 32 and 52.

The element matrix PETSc applies is the same for every cell (it comes from
the reference element alone). utilities.closure_interpolation builds it on
one reference cell, on one rank, applies it to each local cell's closure,
and hands each node its owner's value through a local-to-global and back.
That covers every field whose edge and face points carry at most one node
(any discontinuous field, continuous degree 1 and 2). A field with more
nodes per edge, ordered by the edge's orientation, keeps PETSc's matrix,
and on a partition where that matrix cannot be assembled it raises on
every rank instead of hanging.

Against PETSc's matrix at np 1 and 4, on triangles, quads, tets, hexes, a
degree-2 annulus and a spherical shell, continuous and discontinuous
degree 1 and 2: identical bits for continuous P1/P2 on simplices and Q1,
round-off (at most 1.6e-15) elsewhere, 2-7x faster. test_1107 passes at
np 2 and 4.

The remaining interpolations of this kind are marked TODO(BUG): the
degree >= 3 trace coordinates of boundary_flux and PETSc's projection of a
degree >= 2 coordinate field in Mesh setup.

Underworld development team with AI support from Claude Code
---
 .../discretisation/discretisation_mesh.py     |  41 ++---
 src/underworld3/function/field_projection.py  |  47 ++---
 src/underworld3/utilities/boundary_flux.py    |   5 +
 .../utilities/closure_interpolation.py        | 160 ++++++++++++++++++
 4 files changed, 197 insertions(+), 56 deletions(-)
 create mode 100644 src/underworld3/utilities/closure_interpolation.py

diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 199c6a119..9dcce4cca 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -2755,6 +2755,10 @@ def nuke_coords_and_rebuild(
             # DMPlexComputeBdIntegral to segfault/deadlock (issue #96).
             # createCoordinateSpace builds the FE internally with correct
             # subspace initialisation.
+            # TODO(BUG): #902 — for a degree >= 2 coordinate field PETSc projects
+            # the coordinates through an interpolation matrix between two DMs on
+            # this plex, which cannot be assembled on a rank that holds cells but
+            # owns none of their nodes (the job hangs). Degree 1 is a copy.
             self.dm.createCoordinateSpace(self.degree, False, True)
 
             # Issue #96 fix: Force coordinate field creation and strip
@@ -5729,36 +5733,23 @@ def _cell_node_indices(self, degree, continuous):
         return self._cell_node_array[key]
 
     def _get_coords_for_basis(self, degree, continuous):
+        """The coordinates of the nodes of a degree-``degree`` Lagrange field, an
+        ``(n, cdim)`` array in the local layout of :meth:`_basis_coordinate_dm`
+        (ghost nodes included).
+
+        The mesh's coordinate field is interpolated cell by cell, not through
+        PETSc's interpolation matrix, which cannot be assembled on a rank that
+        holds cells but owns none of their nodes (#902): see
+        ``utilities.closure_interpolation``.
         """
-        This function returns the vertex array for the
-        provided variable. If the array does not already exist,
-        it is first created and then returned.
-        """
+        from underworld3.utilities.closure_interpolation import interpolate_within_plex
 
         dmold = self.dm.getCoordinateDM()
         dmnew = self._basis_coordinate_dm(degree, continuous)
-
-        matInterp, vecScale = dmold.createInterpolation(dmnew)
-        coordsOld = self.dm.getCoordinates()
-        coordsNewL = dmnew.getLocalVec()
-        coordsNewG = dmnew.getGlobalVec()
-        matInterp.mult(coordsOld, coordsNewG)
-        dmnew.globalToLocal(coordsNewG, coordsNewL)
-
-        arr = coordsNewL.array
-        arrcopy = arr.reshape(-1, self.cdim).copy()
-
-        dmnew.restoreGlobalVec(coordsNewG)
-        dmnew.restoreLocalVec(coordsNewL)
-        # Clean up the PETSc interpolation objects built above. Without this
-        # they accumulate until Python GC runs — noticeable in long adapt
-        # loops that re-fill the coord cache per variable.
-        matInterp.destroy()
-        if vecScale is not None:
-            vecScale.destroy()
+        coords = interpolate_within_plex(
+            self.dm, dmold, self.dm.getCoordinatesLocal().array, dmnew, self.cdim)
         dmnew.destroy()
-
-        return arrcopy
+        return coords
 
     def _coord_rows_for_points(self, nav_dm, points):
         """Row indices into the navigation coordinate array (``_nav_coords``)
diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py
index d8922776e..37d9fb1b4 100644
--- a/src/underworld3/function/field_projection.py
+++ b/src/underworld3/function/field_projection.py
@@ -1,8 +1,8 @@
 """Project or prolong MeshVariable data to a different polynomial degree.
 
-Uses a scratch PETSc DM with ``createInterpolation`` — no MeshVariable
-is created, the mesh DM is never modified, and all scratch objects are
-destroyed before returning.
+Interpolates onto a scratch PETSc DM, cell by cell
+(``utilities.closure_interpolation``) — no MeshVariable is created, the mesh
+DM is never modified, and all scratch objects are destroyed before returning.
 
 Typical use cases
 -----------------
@@ -19,6 +19,8 @@
 import numpy as np
 from petsc4py import PETSc
 
+from underworld3.utilities.closure_interpolation import interpolate_within_plex
+
 if TYPE_CHECKING:
     from underworld3.discretisation import MeshVariable
 
@@ -51,12 +53,12 @@ def project_to_degree(
 
     Notes
     -----
-    This creates a transient scratch DM, builds the PETSc interpolation
-    matrix, applies it, and destroys everything.  The mesh DM and all
-    existing MeshVariables are completely untouched.
+    This creates a transient scratch DM, applies PETSc's element
+    interpolation matrix to each cell, and destroys everything.  The mesh DM
+    and all existing MeshVariables are completely untouched.
 
-    For ``target_degree == mesh_var.degree`` the interpolation matrix is
-    the identity and the result matches the source data exactly.
+    For ``target_degree == mesh_var.degree`` the element matrix is the
+    identity and the result matches the source data exactly.
     """
 
     mesh = mesh_var.mesh
@@ -80,31 +82,14 @@ def project_to_degree(
     # --- source sub-DM for this variable's field ---
     iset, subdm_src = mesh.dm.createSubDM(mesh_var.field_id)
 
-    # --- interpolation matrix ---
-    # PETSc: source_subdm.createInterpolation(target_dm) returns a matrix
-    # whose mult() maps source global vec → target global vec.
-    interp, _scale = subdm_src.createInterpolation(dm_scratch)
-
-    # --- apply ---
-    g_src = subdm_src.createGlobalVec()
-    subdm_src.localToGlobal(mesh_var._lvec, g_src)
-
-    g_dst = dm_scratch.createGlobalVec()
-    interp.mult(g_src, g_dst)
-
-    if include_ghosts:
-        l_dst = dm_scratch.createLocalVec()
-        dm_scratch.globalToLocal(g_dst, l_dst)
-        result = l_dst.array.reshape(-1, nc).copy()
-        l_dst.destroy()
-    else:
-        result = g_dst.array.reshape(-1, nc).copy()
+    # Cell by cell rather than through PETSc's interpolation matrix, which
+    # cannot be assembled on a rank that holds cells but owns none of their
+    # degrees of freedom (#902).
+    result = interpolate_within_plex(
+        mesh.dm, subdm_src, mesh_var._lvec.array, dm_scratch, nc, ghosts=include_ghosts)
 
-    # --- cleanup ---
-    for obj in (interp, g_src, g_dst, dm_scratch, fe_target, subdm_src, iset):
+    for obj in (dm_scratch, fe_target, subdm_src, iset):
         obj.destroy()
-    if _scale is not None:
-        _scale.destroy()
 
     return result
 
diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py
index 6e5d707b2..2483a8287 100644
--- a/src/underworld3/utilities/boundary_flux.py
+++ b/src/underworld3/utilities/boundary_flux.py
@@ -154,6 +154,11 @@ def _trace_interior_coords(solver, degree):
         mesh.dim, cdim, mesh.isSimplex, mesh.qdegree, prefix, PETSc.COMM_SELF)
     dmnew.setField(0, fe)
     dmnew.createDS()
+    # TODO(BUG): #902 — this matrix cannot be assembled on a rank that holds cells
+    # but owns none of their nodes (PETSc skips that rank's preallocation and the
+    # job hangs). A degree >= 3 trace has two or more nodes per edge, ordered by
+    # the edge's orientation, so utilities.closure_interpolation's cell-by-cell
+    # route does not cover it yet.
     mat_interp, vec_scale = dmold.createInterpolation(dmnew)
     coords_new_g = dmnew.getGlobalVec()
     coords_new_l = dmnew.getLocalVec()
diff --git a/src/underworld3/utilities/closure_interpolation.py b/src/underworld3/utilities/closure_interpolation.py
new file mode 100644
index 000000000..883589045
--- /dev/null
+++ b/src/underworld3/utilities/closure_interpolation.py
@@ -0,0 +1,160 @@
+"""Interpolation between two fields on one plex, applied cell by cell (#902).
+
+Two DMs that share a plex (a coordinate DM and a basis DM built from it, a
+variable's sub-DM and a scratch clone of the mesh) are interpolated by PETSc's
+``DMPlexComputeInterpolatorNested``. Its element matrix comes from the reference
+element alone, so it is the same for every cell, and the parallel matrix it
+assembles is only that element matrix scattered over the cells' closures. Its
+preallocation pass is skipped on a rank with no local rows and no local columns,
+while its fill pass still runs over that rank's cells. A shared point belongs to
+the highest rank holding it, so a rank whose cells lie wholly inside higher
+ranks' regions owns nothing: its entries reach owners that never preallocated
+them, the owner raises inside ``MatAssemblyEnd``, and the other ranks wait in
+that assembly's collective for ever. Measured on kaiju at np 32 and 52 on a
+124-cell coarse level.
+
+Here PETSc builds the element matrix on a single reference cell, on one rank,
+and it is applied to each local cell's closure values directly. The values at a
+node are those of its owner (the result is sent to the owners and back), as
+they were when the matrix was multiplied.
+"""
+import numpy as np
+from mpi4py import MPI
+from petsc4py import PETSc
+
+from underworld3.utilities._plex_topology import closure_dof_rows
+
+_REFERENCE_CELLS = {
+    # (dim, simplex): vertices of one cell, any valid order (the element
+    # matrix depends on the cell's topology, not on its geometry)
+    (2, True): [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
+    (3, True): [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
+    (2, False): [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]],
+    (3, False): [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0],
+                 [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 1.0], [0.0, 1.0, 1.0]],
+}
+
+
+def _reference_element_matrix(fe_src, fe_dst, dim, cdim, simplex, block):
+    """PETSc's element interpolation matrix from ``fe_src`` to ``fe_dst``, in
+    closure order of nodes (one row per destination node, one column per source
+    node, the same for every component), or ``None`` when a point other than the
+    cell carries more than one node: those nodes are ordered by the point's
+    orientation, which differs from cell to cell."""
+    verts = np.asarray(_REFERENCE_CELLS[(dim, simplex)])
+    coords = np.zeros((verts.shape[0], cdim))
+    coords[:, :dim] = verts
+    cells = np.arange(verts.shape[0], dtype=PETSc.IntType).reshape(1, -1)
+    ref = PETSc.DMPlex().createFromCellList(dim, cells, coords, interpolate=True,
+                                            comm=PETSc.COMM_SELF)
+    pair = []
+    for fe in (fe_src, fe_dst):
+        rdm = ref.clone()
+        rdm.setField(0, fe)
+        rdm.createDS()
+        pair.append(rdm)
+    try:
+        cell = ref.getHeightStratum(0)[0]
+        points = ref.getTransitiveClosure(cell)[0]
+        for rdm in pair:
+            sec = rdm.getLocalSection()
+            if any(sec.getDof(p) > block for p in points if p != cell):
+                return None
+        interp, scale = pair[0].createInterpolation(pair[1])
+        full = interp.convert("dense").getDenseArray().copy()
+        interp.destroy()
+        if scale is not None:
+            scale.destroy()
+        rows_src, rows_dst = (closure_dof_rows(rdm, rdm.getLocalSection(), [cell], block)[0]
+                              for rdm in pair)
+        return full[np.ix_(rows_dst * block, rows_src * block)]
+    finally:
+        for rdm in pair:
+            rdm.destroy()
+        ref.destroy()
+
+
+def interpolate_within_plex(plex, dm_src, x_src, dm_dst, block, ghosts=True):
+    """Field 0 of ``dm_dst`` interpolated from ``x_src``, the local array of
+    field 0 of ``dm_src``, where both DMs share the topology of ``plex`` (the
+    mesh's DMPlex; a sub-DM is a plain DM in petsc4py, without the plex queries).
+
+    The result is ``dm_dst``'s local array, shaped ``(n, block)`` (``ghosts``
+    True) or its owned part (``ghosts`` False), with each node's value computed
+    by the rank that owns it. Collective on the DMs' communicator.
+    """
+    fe_src = dm_src.getField(0)[0]
+    fe_dst = dm_dst.getField(0)[0]
+    element = _reference_element_matrix(fe_src, fe_dst, plex.getDimension(),
+                                        plex.getCoordinateDim(), plex.isSimplex(), block)
+
+    cS, cE = plex.getHeightStratum(0)
+    cells = np.arange(cS, cE)
+    if element is not None and cells.size:
+        rows_src = closure_dof_rows(dm_src, dm_src.getLocalSection(), cells, block)
+        rows_dst = closure_dof_rows(dm_dst, dm_dst.getLocalSection(), cells, block)
+        uniform = (rows_src.shape[1] == element.shape[1]
+                   and rows_dst.shape[1] == element.shape[0]
+                   and (rows_src >= 0).all() and (rows_dst >= 0).all())
+    else:
+        uniform = element is not None
+    # one verdict for every rank: the two routes below are collective
+    mpi_comm = dm_src.comm.tompi4py()
+    if not mpi_comm.allreduce(bool(uniform), op=MPI.LAND):
+        return _interpolate_by_matrix(plex, dm_src, x_src, dm_dst, block, ghosts)
+
+    lvec = dm_dst.getLocalVec()
+    gvec = dm_dst.getGlobalVec()
+    out = lvec.getArray().reshape(-1, block)
+    if cells.size:
+        x = np.asarray(x_src).reshape(-1, block)
+        out[rows_dst.ravel()] = np.einsum("ab,cbk->cak", element,
+                                          x[rows_src]).reshape(-1, block)
+    # INSERT copies only the points a rank owns, so each node keeps its owner's value
+    dm_dst.localToGlobal(lvec, gvec, addv=PETSc.InsertMode.INSERT_VALUES)
+    if ghosts:
+        dm_dst.globalToLocal(gvec, lvec, addv=PETSc.InsertMode.INSERT_VALUES)
+        result = lvec.getArray().reshape(-1, block).copy()
+    else:
+        result = gvec.getArray().reshape(-1, block).copy()
+    dm_dst.restoreGlobalVec(gvec)
+    dm_dst.restoreLocalVec(lvec)
+    return result
+
+
+def _interpolate_by_matrix(plex, dm_src, x_src, dm_dst, block, ghosts):
+    """PETSc's own interpolation matrix, for the fields whose edge or face nodes
+    are ordered by orientation. It fails on a rank that holds cells but owns no
+    row and no column, so that partition raises here, on every rank, instead."""
+    rows = dm_dst.getGlobalSection().getConstrainedStorageSize()
+    cols = dm_src.getGlobalSection().getConstrainedStorageSize()
+    cS, cE = plex.getHeightStratum(0)
+    mpi_comm = dm_src.comm.tompi4py()
+    if mpi_comm.allreduce(bool(cE > cS and rows == 0 and cols == 0), op=MPI.LOR):
+        raise NotImplementedError(
+            "interpolation between two fields of one mesh: a rank holds cells but "
+            "owns none of their degrees of freedom, where PETSc's interpolation "
+            "matrix cannot be assembled (#902), and this field has more than one "
+            "node on an edge or face. Run on fewer ranks or repartition the mesh.")
+    gsrc = dm_src.getGlobalVec()
+    lsrc = dm_src.getLocalVec()
+    lsrc.getArray()[:] = np.asarray(x_src).ravel()
+    dm_src.localToGlobal(lsrc, gsrc, addv=PETSc.InsertMode.INSERT_VALUES)
+    interp, scale = dm_src.createInterpolation(dm_dst)
+    gdst = dm_dst.getGlobalVec()
+    interp.mult(gsrc, gdst)
+    if ghosts:
+        ldst = dm_dst.getLocalVec()
+        dm_dst.globalToLocal(gdst, ldst, addv=PETSc.InsertMode.INSERT_VALUES)
+        result = ldst.getArray().reshape(-1, block).copy()
+        dm_dst.restoreLocalVec(ldst)
+    else:
+        result = gdst.getArray().reshape(-1, block).copy()
+    dm_dst.restoreGlobalVec(gdst)
+    dm_src.restoreLocalVec(lsrc)
+    dm_src.restoreGlobalVec(gsrc)
+    interp.destroy()
+    if scale is not None:
+        scale.destroy()
+    return result
+

From f4b9128acf45fbab084384f74bc966bc47fdf433 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 13:49:47 +1100
Subject: [PATCH 04/13] Agree a transfer-builder failure on every rank before
 retrying (#886)

In parallel build_transfers caught any exception from the hierarchy build
on each rank alone, and that rank retried with the RBF builder (or degraded
to the default preconditioner) while the others went on into the transfer's
assembly: the ranks met in different collectives and the run hung.

The node-transfer builders are rank-local work, so their verdict is agreed
right after them (_agreed_node_transfer, used by both parallel builds): if
the builder raised on any rank, TransferBuildError is raised on every rank,
with each failing rank's message. In parallel build_transfers retries or
degrades only on that error; any other exception is one rank's, the others
are inside the build's collectives, and it is raised where it happens.
Serial keeps the retry for every failure.

test_1108 (the barycentric builder raising on rank 1 only) passes at np 2
and 4: every rank takes the RBF rescue and the solve converges on
geometric multigrid.

Underworld development team with AI support from Claude Code
---
 src/underworld3/utilities/custom_mg.py | 42 ++++++++++++++++++++++----
 1 file changed, 36 insertions(+), 6 deletions(-)

diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py
index 53c77d1bf..666f03902 100644
--- a/src/underworld3/utilities/custom_mg.py
+++ b/src/underworld3/utilities/custom_mg.py
@@ -224,6 +224,11 @@ def _drop_structural_zeros(P_csr, tol=1e-12):
 _BUILDERS = {"barycentric": barycentric_prolongation, "rbf": rbf_prolongation}
 
 
+class TransferBuildError(RuntimeError):
+    """A node-transfer builder failed on at least one rank; raised on every rank,
+    so the ranks can retry or degrade together (#886)."""
+
+
 def _require_serial(where):
     """Custom-P transfers are serial-only (experimental). The reduced maps use
     rank-local DOF indices and the prolongations assemble as serial AIJ; at np>1
@@ -906,6 +911,25 @@ def _node_transfer_no_repair(builder, coarse_coords, fine_coords):
     return builder(coarse_coords, fine_coords)
 
 
+def _agreed_node_transfer(builder, coarse_coords, fine_coords, comm):
+    """The rank-local node transfer of a parallel build, with one verdict for
+    every rank: if the builder raised on any rank, :class:`TransferBuildError`
+    is raised on all of them. The builders do no communication, so every rank
+    reaches the verdict; a rank that kept its own error would go on to retry
+    while the others assemble the transfer, and the two would meet in different
+    collectives (#886)."""
+    problem = None
+    try:
+        Pn = _drop_structural_zeros(
+            _node_transfer_no_repair(builder, coarse_coords, fine_coords))
+    except Exception as exc:    # sanctioned: the verdict below raises it on every rank
+        Pn, problem = None, f"rank {comm.getRank()}: {exc}"
+    problems = [p for p in comm.tompi4py().allgather(problem) if p is not None]
+    if problems:
+        raise TransferBuildError("; ".join(problems))
+    return Pn
+
+
 def _assemble_node_transfer(Pn, coarse_cols, fine_layout, coarse_layout, ncomp,
                             comm):
     """A node-level transfer ``Pn`` (scipy CSR, fine local nodes by coarse
@@ -955,8 +979,7 @@ def _build_parallel_transfer(coarse_coords, fine_coords, coarse_layout,
     coarse ``l2g`` to global columns (off-rank columns are fine for MPIAIJ).
     Constrained coarse DOFs (``l2g == -1``) drop out -> reduced->reduced."""
     l2g_c = coarse_layout.l2g
-    Pn = _drop_structural_zeros(
-        _node_transfer_no_repair(builder, coarse_coords, fine_coords))  # (n_f_nodes, n_c_nodes), local
+    Pn = _agreed_node_transfer(builder, coarse_coords, fine_coords, comm)  # (n_f_nodes, n_c_nodes), local
     return _assemble_node_transfer(
         Pn, np.asarray(l2g_c, dtype=np.int64).reshape(-1, ncomp), fine_layout,
         coarse_layout, ncomp, comm)
@@ -998,8 +1021,7 @@ def _build_crosspart_transfer(coarse_coords, fine_coords, coarse_layout,
     coords_u, cols_u = _gather_coarse_cloud(coarse_coords, coarse_layout,
                                             ncomp, comm)
 
-    Pn = _drop_structural_zeros(
-        _node_transfer_no_repair(builder, coords_u, fine_coords))   # (n_f_local, Nu)
+    Pn = _agreed_node_transfer(builder, coords_u, fine_coords, comm)   # (n_f_local, Nu)
     return _assemble_node_transfer(Pn, np.asarray(cols_u, dtype=np.int64),
                                    fine_layout, coarse_layout, ncomp, comm)
 
@@ -1961,7 +1983,9 @@ def build_transfers(solver, field_id=None):
     is a DEMAND (build errors raise — the user registered it explicitly); the
     other two are PREFERENCES, built through the same opportunistic arm
     (barycentric, RBF retry, degrade to the solver's default preconditioner —
-    every step recorded in ``solver.pc_fallbacks``).
+    every step recorded in ``solver.pc_fallbacks``). In parallel only a builder
+    failure, which every rank sees (:class:`TransferBuildError`), is retried
+    or degraded; any other failure is raised where it happens (#886).
 
     A refinement child carries ``mesh._custom_mg_coarse_meshes`` (the static
     coarse tail), so a :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]``
@@ -2072,6 +2096,12 @@ def build_transfers(solver, field_id=None):
     # barycentric fell back to GAMG at 23 its while RBF kept pc=mg at 2.
     # Falling back to GAMG loses the hierarchy entirely, so try the cheaper
     # degradation first. See #424.
+    # In parallel the retry and the degrade are collective decisions, so they
+    # are taken only on a verdict every rank holds: a builder failure, which
+    # is raised on every rank. Any other exception is one rank's: the others
+    # are inside the build's collectives, and a rank that caught it and
+    # retried would meet them in different ones (#886). It is raised.
+    _retryable = TransferBuildError if solver.dm.comm.getSize() > 1 else Exception
     _attempts = [builder] + (["rbf"] if builder != "rbf" else [])
     h = Ps = None
     for _i, _b in enumerate(_attempts):
@@ -2084,7 +2114,7 @@ def build_transfers(solver, field_id=None):
         try:
             Ps = h.build(solver)
             break
-        except Exception as exc:            # pragma: no cover - defensive
+        except _retryable as exc:   # sanctioned: retried or degraded, and recorded
             import warnings
             if _i + 1 < len(_attempts):
                 solver._record_pc_fallback(

From 89037a1668bf53d5aab2d95337a0e674e7eba143 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 14:30:30 +1100
Subject: [PATCH 05/13] Interpolate within one mesh through UW3's copy of
 PETSc's nested interpolator (#902)

The cell-by-cell Python route (37ac408a) is replaced by a copy of PETSc's
DMPlexComputeInterpolatorNested (v3.25.0, BSD-2-Clause), trimmed to the
same-plex finite-element case, with its preallocation loop run over the same
cells as its fill loop: UW_DMCreateInterpolationWithinPlex in
cython/petsc_compat.h, petsc_dm_create_interpolation_within_plex in
petsc_discretisation. PETSc has guarded only the preallocation by the local
matrix sizes since b4ec2976050 (v3.24.0; still so in v3.25.6, v3.26.0 and
main), so a rank that holds cells but owns no rows and no columns fills
entries nobody preallocated.

Mesh._get_coords_for_basis, function.project_to_degree and
boundary_flux._trace_interior_coords call it in place of createInterpolation.
It works with any PETSc build, covers degree 3 (two nodes per edge, which the
Python route could not), and reproduces PETSc's matrix exactly: max |A - B| 0
and the same nonzeros on triangles, quads, tets, hexes, a degree-2 annulus and
a spherical shell, continuous and discontinuous degree 1 to 3, at np 1 and 4,
so the node coordinates are bit-identical to development.

Every DM section is created on every rank before anything rank-dependent.
The Python route created them only on ranks holding cells, so an empty rank
went on to the next collective while the others waited (the review found
test_0774/0781 at np 2 and test_0779/1103 at np 4 hanging); all four pass.

test_1107 now also pins degree 3 on the forced partition, and its partition
check is one verdict for every rank. PETSc's own projection of a degree >= 2
coordinate field in Mesh setup stays marked TODO(BUG).

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h         | 141 +++++++++++++++
 .../cython/petsc_discretisation.pyx           |  18 ++
 src/underworld3/cython/petsc_extras.pxi       |   1 +
 .../discretisation/discretisation_mesh.py     |  32 +++-
 src/underworld3/function/field_projection.py  |  49 ++++--
 src/underworld3/utilities/boundary_flux.py    |  14 +-
 .../utilities/closure_interpolation.py        | 160 ------------------
 ...erpolation_rank_owning_nothing_parallel.py |  55 ++++--
 8 files changed, 263 insertions(+), 207 deletions(-)
 delete mode 100644 src/underworld3/utilities/closure_interpolation.py

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 74075e043..5eb5993d6 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -366,3 +366,144 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, DM *sandbox)
     *sandbox = sdm;
     PetscFunctionReturn(PETSC_SUCCESS);
 }
+
+// Interpolation between two DMs on ONE plex (a coordinate DM and a basis DM
+// cloned from it, a variable's sub-DM and a scratch clone of the mesh): the
+// matrix DMCreateInterpolation builds through DMPlexComputeInterpolatorNested
+// (isRefined = PETSC_FALSE), copied from PETSc v3.25.0 (BSD-2-Clause) with one
+// change (#902). PETSc runs its preallocation loop only on a rank with local
+// rows or columns (`if (locRows || locCols)`, since b4ec2976050, v3.24.0) but
+// its fill loop on every rank that holds cells. A shared point belongs to the
+// highest rank holding it, so a rank whose cells lie inside higher ranks'
+// regions owns no rows and no columns; its cells' entries then reach owners
+// that never preallocated them ("New nonzero caused a malloc" in
+// MatAssemblyEnd) and the ranks that did not raise wait in the assembly's
+// collective. Here both loops run over the same cells. Finite-element fields
+// only, which are the only ones UW3 interpolates this way. Drop this for
+// DMCreateInterpolation once PETSc runs the two loops alike.
+PetscErrorCode UW_DMCreateInterpolationWithinPlex(DM dmc, DM dmf, Mat *In)
+{
+    PetscSection fsection, fglobalSection, csection, cglobalSection;
+    PetscScalar *elemMat, *vals;
+    PetscFE     *fe;
+    PetscInt     dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
+    PetscInt     cTotDim = 0, rTotDim = 0, m, n, locRows, locCols;
+    MatType      mtype;
+    PetscBool    ismatis;
+    Mat          preallocator;
+
+    PetscFunctionBeginUser;
+    PetscCall(DMGetGlobalSection(dmf, &fglobalSection));
+    PetscCall(PetscSectionGetConstrainedStorageSize(fglobalSection, &m));
+    PetscCall(DMGetGlobalSection(dmc, &cglobalSection));
+    PetscCall(PetscSectionGetConstrainedStorageSize(cglobalSection, &n));
+    PetscCall(DMGetMatType(dmc, &mtype));
+    PetscCall(PetscStrcmp(mtype, MATIS, &ismatis));
+    PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), In));
+    PetscCall(MatSetSizes(*In, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
+    PetscCall(MatSetType(*In, ismatis ? MATAIJ : mtype));
+
+    PetscCall(DMGetDimension(dmf, &dim));
+    PetscCall(DMGetLocalSection(dmf, &fsection));
+    PetscCall(DMGetLocalSection(dmc, &csection));
+    PetscCall(PetscSectionGetNumFields(fsection, &Nf));
+    PetscCall(DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd));
+    PetscCall(PetscCalloc1(Nf, &fe));
+    for (f = 0; f < Nf; ++f) {
+        PetscObject  obj, objc;
+        PetscClassId id, idc;
+        PetscInt     rNb, cNb;
+
+        PetscCall(DMGetField(dmf, f, NULL, &obj));
+        PetscCall(PetscObjectGetClassId(obj, &id));
+        PetscCall(DMGetField(dmc, f, NULL, &objc));
+        PetscCall(PetscObjectGetClassId(objc, &idc));
+        PetscCheck(id == PETSCFE_CLASSID && idc == PETSCFE_CLASSID, PetscObjectComm((PetscObject)dmc),
+                   PETSC_ERR_SUP, "UW_DMCreateInterpolationWithinPlex: field %" PetscInt_FMT " is not a PetscFE", f);
+        fe[f] = (PetscFE)obj;
+        PetscCall(PetscFEGetDimension(fe[f], &rNb));
+        PetscCall(PetscFEGetDimension((PetscFE)objc, &cNb));
+        rTotDim += rNb;
+        cTotDim += cNb;
+    }
+    PetscCall(PetscCalloc1(rTotDim * cTotDim, &elemMat));
+    for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
+        PetscDualSpace   Qref;
+        PetscQuadrature  quad;
+        const PetscReal *qpoints, *qweights;
+        PetscReal       *points;
+        PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
+
+        // the fine field's dual basis: every functional's quadrature points
+        PetscCall(PetscFEGetDualSpace(fe[fieldI], &Qref));
+        PetscCall(PetscFEGetNumComponents(fe[fieldI], &Nc));
+        PetscCall(PetscDualSpaceGetDimension(Qref, &fpdim));
+        for (i = 0; i < fpdim; ++i) {
+            PetscCall(PetscDualSpaceGetFunctional(Qref, i, &quad));
+            PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Np, NULL, NULL));
+            npoints += Np;
+        }
+        PetscCall(PetscMalloc1(npoints * dim, &points));
+        for (i = 0, k = 0; i < fpdim; ++i) {
+            PetscCall(PetscDualSpaceGetFunctional(Qref, i, &quad));
+            PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Np, &qpoints, NULL));
+            for (p = 0; p < Np; ++p, ++k)
+                for (d = 0; d < dim; ++d) points[k * dim + d] = qpoints[p * dim + d];
+        }
+        for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
+            PetscObject     obj;
+            PetscFE         fec;
+            PetscTabulation T = NULL;
+            PetscInt        NcJ, cpdim, j, qNc;
+
+            PetscCall(DMGetField(dmc, fieldJ, NULL, &obj));
+            fec = (PetscFE)obj;
+            PetscCall(PetscFEGetNumComponents(fec, &NcJ));
+            PetscCall(PetscFEGetDimension(fec, &cpdim));
+            // fields only interpolate themselves: the coarse basis evaluated at
+            // the fine functionals' points, weighted by the functionals
+            if (fieldI == fieldJ) {
+                PetscCheck(Nc == NcJ, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %" PetscInt_FMT " does not match coarse field %" PetscInt_FMT, Nc, NcJ);
+                PetscCall(PetscFECreateTabulation(fec, 1, npoints, points, 0, &T));
+                for (i = 0, k = 0; i < fpdim; ++i) {
+                    PetscCall(PetscDualSpaceGetFunctional(Qref, i, &quad));
+                    PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Np, NULL, &qweights));
+                    PetscCheck(qNc == NcJ, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %" PetscInt_FMT " does not match coarse field %" PetscInt_FMT, qNc, NcJ);
+                    for (p = 0; p < Np; ++p, ++k)
+                        for (j = 0; j < cpdim; ++j)
+                            for (c = 0; c < Nc; ++c) elemMat[(offsetI + i) * cTotDim + offsetJ + j] += T->T[0][k * cpdim * NcJ + j * Nc + c] * qweights[p * qNc + c];
+                }
+                PetscCall(PetscTabulationDestroy(&T));
+            }
+            offsetJ += cpdim;
+        }
+        offsetI += fpdim;
+        PetscCall(PetscFree(points));
+    }
+
+    // Preallocate from every cell this rank holds, as the fill below does,
+    // whether or not the rank owns any rows or columns (the one change).
+    PetscCall(MatGetLocalSize(*In, &locRows, &locCols));
+    PetscCall(MatCreate(PetscObjectComm((PetscObject)*In), &preallocator));
+    PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
+    PetscCall(MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE));
+    PetscCall(MatSetUp(preallocator));
+    PetscCall(PetscCalloc1(rTotDim * cTotDim, &vals));
+    for (c = cStart; c < cEnd; ++c)
+        PetscCall(DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, PETSC_FALSE, dmc, csection, cglobalSection, PETSC_FALSE, preallocator, c, vals, INSERT_VALUES));
+    PetscCall(PetscFree(vals));
+    PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
+    PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
+    PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *In));
+    PetscCall(MatDestroy(&preallocator));
+
+    PetscCall(MatZeroEntries(*In));
+    for (c = cStart; c < cEnd; ++c)
+        PetscCall(DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, PETSC_FALSE, dmc, csection, cglobalSection, PETSC_FALSE, *In, c, elemMat, INSERT_VALUES));
+    PetscCall(PetscFree(fe));
+    PetscCall(PetscFree(elemMat));
+    PetscCall(MatAssemblyBegin(*In, MAT_FINAL_ASSEMBLY));
+    PetscCall(MatAssemblyEnd(*In, MAT_FINAL_ASSEMBLY));
+    PetscCall(MatViewFromOptions(*In, NULL, "-interp_mat_view"));
+    PetscFunctionReturn(PETSC_SUCCESS);
+}
diff --git a/src/underworld3/cython/petsc_discretisation.pyx b/src/underworld3/cython/petsc_discretisation.pyx
index a15613057..ffa52d935 100644
--- a/src/underworld3/cython/petsc_discretisation.pyx
+++ b/src/underworld3/cython/petsc_discretisation.pyx
@@ -181,6 +181,24 @@ def petsc_dm_set_regular_refinement(dm, regular=True):
         return
 
 
+def petsc_dm_create_interpolation_within_plex(dmc, dmf):
+        """The interpolation matrix from the field of ``dmc`` to the field of
+        ``dmf``, two DMs on ONE plex (clones or sub-DMs of one mesh DM).
+
+        It is the matrix ``dmc.createInterpolation(dmf)`` returns, assembled
+        also on a rank that holds cells but owns none of their degrees of
+        freedom, where PETSc's own preallocation skips the rank and the
+        assembly fails on some ranks and hangs the others (#902; see
+        ``UW_DMCreateInterpolationWithinPlex`` in ``petsc_compat.h``).
+        Collective.
+        """
+        cdef DM c_dmc = dmc
+        cdef DM c_dmf = dmf
+        cdef Mat interp = PETSc.Mat()
+        CHKERRQ( UW_DMCreateInterpolationWithinPlex(c_dmc.dm, c_dmf.dm, &interp.mat) )
+        return interp
+
+
 
 # This is not cython, does it need to be here or in discretisation.py ?
 
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 4c827d2e2..4e39315b7 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -58,6 +58,7 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode UW_DMForceCoordinateField(PetscDM)
     PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*)
     PetscErrorCode UW_DMCreateBdIntegralSandbox(PetscDM, PetscDM*)
+    PetscErrorCode UW_DMCreateInterpolationWithinPlex(PetscDM, PetscDM, PetscMat*)
 
 cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[])
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 9dcce4cca..25fc70028 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -5737,19 +5737,35 @@ def _get_coords_for_basis(self, degree, continuous):
         ``(n, cdim)`` array in the local layout of :meth:`_basis_coordinate_dm`
         (ghost nodes included).
 
-        The mesh's coordinate field is interpolated cell by cell, not through
-        PETSc's interpolation matrix, which cannot be assembled on a rank that
-        holds cells but owns none of their nodes (#902): see
-        ``utilities.closure_interpolation``.
+        The interpolation matrix comes from UW3's copy of PETSc's nested
+        interpolator, which assembles also on a rank that holds cells but owns
+        none of their nodes (#902; ``petsc_dm_create_interpolation_within_plex``).
         """
-        from underworld3.utilities.closure_interpolation import interpolate_within_plex
+        from underworld3.cython.petsc_discretisation import (
+            petsc_dm_create_interpolation_within_plex,
+        )
 
         dmold = self.dm.getCoordinateDM()
         dmnew = self._basis_coordinate_dm(degree, continuous)
-        coords = interpolate_within_plex(
-            self.dm, dmold, self.dm.getCoordinatesLocal().array, dmnew, self.cdim)
+
+        matInterp = petsc_dm_create_interpolation_within_plex(dmold, dmnew)
+        coordsOld = self.dm.getCoordinates()
+        coordsNewL = dmnew.getLocalVec()
+        coordsNewG = dmnew.getGlobalVec()
+        matInterp.mult(coordsOld, coordsNewG)
+        dmnew.globalToLocal(coordsNewG, coordsNewL)
+
+        arrcopy = coordsNewL.array.reshape(-1, self.cdim).copy()
+
+        dmnew.restoreGlobalVec(coordsNewG)
+        dmnew.restoreLocalVec(coordsNewL)
+        # Clean up the PETSc interpolation objects built above. Without this
+        # they accumulate until Python GC runs — noticeable in long adapt
+        # loops that re-fill the coord cache per variable.
+        matInterp.destroy()
         dmnew.destroy()
-        return coords
+
+        return arrcopy
 
     def _coord_rows_for_points(self, nav_dm, points):
         """Row indices into the navigation coordinate array (``_nav_coords``)
diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py
index 37d9fb1b4..0c6e95ded 100644
--- a/src/underworld3/function/field_projection.py
+++ b/src/underworld3/function/field_projection.py
@@ -1,8 +1,8 @@
 """Project or prolong MeshVariable data to a different polynomial degree.
 
-Interpolates onto a scratch PETSc DM, cell by cell
-(``utilities.closure_interpolation``) — no MeshVariable is created, the mesh
-DM is never modified, and all scratch objects are destroyed before returning.
+Uses a scratch PETSc DM with ``createInterpolation`` — no MeshVariable
+is created, the mesh DM is never modified, and all scratch objects are
+destroyed before returning.
 
 Typical use cases
 -----------------
@@ -19,8 +19,6 @@
 import numpy as np
 from petsc4py import PETSc
 
-from underworld3.utilities.closure_interpolation import interpolate_within_plex
-
 if TYPE_CHECKING:
     from underworld3.discretisation import MeshVariable
 
@@ -53,12 +51,12 @@ def project_to_degree(
 
     Notes
     -----
-    This creates a transient scratch DM, applies PETSc's element
-    interpolation matrix to each cell, and destroys everything.  The mesh DM
-    and all existing MeshVariables are completely untouched.
+    This creates a transient scratch DM, builds the PETSc interpolation
+    matrix, applies it, and destroys everything.  The mesh DM and all
+    existing MeshVariables are completely untouched.
 
-    For ``target_degree == mesh_var.degree`` the element matrix is the
-    identity and the result matches the source data exactly.
+    For ``target_degree == mesh_var.degree`` the interpolation matrix is
+    the identity and the result matches the source data exactly.
     """
 
     mesh = mesh_var.mesh
@@ -82,13 +80,32 @@ def project_to_degree(
     # --- source sub-DM for this variable's field ---
     iset, subdm_src = mesh.dm.createSubDM(mesh_var.field_id)
 
-    # Cell by cell rather than through PETSc's interpolation matrix, which
-    # cannot be assembled on a rank that holds cells but owns none of their
-    # degrees of freedom (#902).
-    result = interpolate_within_plex(
-        mesh.dm, subdm_src, mesh_var._lvec.array, dm_scratch, nc, ghosts=include_ghosts)
+    # --- interpolation matrix ---
+    # Its mult() maps source global vec → target global vec. UW3's copy of
+    # PETSc's nested interpolator: PETSc's own (createInterpolation) cannot be
+    # assembled on a rank that holds cells but owns none of their DOFs (#902).
+    from underworld3.cython.petsc_discretisation import (
+        petsc_dm_create_interpolation_within_plex,
+    )
+    interp = petsc_dm_create_interpolation_within_plex(subdm_src, dm_scratch)
+
+    # --- apply ---
+    g_src = subdm_src.createGlobalVec()
+    subdm_src.localToGlobal(mesh_var._lvec, g_src)
+
+    g_dst = dm_scratch.createGlobalVec()
+    interp.mult(g_src, g_dst)
+
+    if include_ghosts:
+        l_dst = dm_scratch.createLocalVec()
+        dm_scratch.globalToLocal(g_dst, l_dst)
+        result = l_dst.array.reshape(-1, nc).copy()
+        l_dst.destroy()
+    else:
+        result = g_dst.array.reshape(-1, nc).copy()
 
-    for obj in (dm_scratch, fe_target, subdm_src, iset):
+    # --- cleanup ---
+    for obj in (interp, g_src, g_dst, dm_scratch, fe_target, subdm_src, iset):
         obj.destroy()
 
     return result
diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py
index 2483a8287..ec8654ac1 100644
--- a/src/underworld3/utilities/boundary_flux.py
+++ b/src/underworld3/utilities/boundary_flux.py
@@ -154,12 +154,12 @@ def _trace_interior_coords(solver, degree):
         mesh.dim, cdim, mesh.isSimplex, mesh.qdegree, prefix, PETSc.COMM_SELF)
     dmnew.setField(0, fe)
     dmnew.createDS()
-    # TODO(BUG): #902 — this matrix cannot be assembled on a rank that holds cells
-    # but owns none of their nodes (PETSc skips that rank's preallocation and the
-    # job hangs). A degree >= 3 trace has two or more nodes per edge, ordered by
-    # the edge's orientation, so utilities.closure_interpolation's cell-by-cell
-    # route does not cover it yet.
-    mat_interp, vec_scale = dmold.createInterpolation(dmnew)
+    # UW3's copy of PETSc's nested interpolator: PETSc's own cannot be
+    # assembled on a rank that holds cells but owns none of their nodes (#902)
+    from underworld3.cython.petsc_discretisation import (
+        petsc_dm_create_interpolation_within_plex,
+    )
+    mat_interp = petsc_dm_create_interpolation_within_plex(dmold, dmnew)
     coords_new_g = dmnew.getGlobalVec()
     coords_new_l = dmnew.getLocalVec()
     mat_interp.mult(solver.dm.getCoordinates(), coords_new_g)
@@ -176,8 +176,6 @@ def _trace_interior_coords(solver, degree):
     dmnew.restoreGlobalVec(coords_new_g)
     dmnew.restoreLocalVec(coords_new_l)
     mat_interp.destroy()
-    if vec_scale is not None:
-        vec_scale.destroy()
     fe.destroy()
     dmnew.destroy()
     return out
diff --git a/src/underworld3/utilities/closure_interpolation.py b/src/underworld3/utilities/closure_interpolation.py
deleted file mode 100644
index 883589045..000000000
--- a/src/underworld3/utilities/closure_interpolation.py
+++ /dev/null
@@ -1,160 +0,0 @@
-"""Interpolation between two fields on one plex, applied cell by cell (#902).
-
-Two DMs that share a plex (a coordinate DM and a basis DM built from it, a
-variable's sub-DM and a scratch clone of the mesh) are interpolated by PETSc's
-``DMPlexComputeInterpolatorNested``. Its element matrix comes from the reference
-element alone, so it is the same for every cell, and the parallel matrix it
-assembles is only that element matrix scattered over the cells' closures. Its
-preallocation pass is skipped on a rank with no local rows and no local columns,
-while its fill pass still runs over that rank's cells. A shared point belongs to
-the highest rank holding it, so a rank whose cells lie wholly inside higher
-ranks' regions owns nothing: its entries reach owners that never preallocated
-them, the owner raises inside ``MatAssemblyEnd``, and the other ranks wait in
-that assembly's collective for ever. Measured on kaiju at np 32 and 52 on a
-124-cell coarse level.
-
-Here PETSc builds the element matrix on a single reference cell, on one rank,
-and it is applied to each local cell's closure values directly. The values at a
-node are those of its owner (the result is sent to the owners and back), as
-they were when the matrix was multiplied.
-"""
-import numpy as np
-from mpi4py import MPI
-from petsc4py import PETSc
-
-from underworld3.utilities._plex_topology import closure_dof_rows
-
-_REFERENCE_CELLS = {
-    # (dim, simplex): vertices of one cell, any valid order (the element
-    # matrix depends on the cell's topology, not on its geometry)
-    (2, True): [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
-    (3, True): [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
-    (2, False): [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]],
-    (3, False): [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0],
-                 [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 1.0], [0.0, 1.0, 1.0]],
-}
-
-
-def _reference_element_matrix(fe_src, fe_dst, dim, cdim, simplex, block):
-    """PETSc's element interpolation matrix from ``fe_src`` to ``fe_dst``, in
-    closure order of nodes (one row per destination node, one column per source
-    node, the same for every component), or ``None`` when a point other than the
-    cell carries more than one node: those nodes are ordered by the point's
-    orientation, which differs from cell to cell."""
-    verts = np.asarray(_REFERENCE_CELLS[(dim, simplex)])
-    coords = np.zeros((verts.shape[0], cdim))
-    coords[:, :dim] = verts
-    cells = np.arange(verts.shape[0], dtype=PETSc.IntType).reshape(1, -1)
-    ref = PETSc.DMPlex().createFromCellList(dim, cells, coords, interpolate=True,
-                                            comm=PETSc.COMM_SELF)
-    pair = []
-    for fe in (fe_src, fe_dst):
-        rdm = ref.clone()
-        rdm.setField(0, fe)
-        rdm.createDS()
-        pair.append(rdm)
-    try:
-        cell = ref.getHeightStratum(0)[0]
-        points = ref.getTransitiveClosure(cell)[0]
-        for rdm in pair:
-            sec = rdm.getLocalSection()
-            if any(sec.getDof(p) > block for p in points if p != cell):
-                return None
-        interp, scale = pair[0].createInterpolation(pair[1])
-        full = interp.convert("dense").getDenseArray().copy()
-        interp.destroy()
-        if scale is not None:
-            scale.destroy()
-        rows_src, rows_dst = (closure_dof_rows(rdm, rdm.getLocalSection(), [cell], block)[0]
-                              for rdm in pair)
-        return full[np.ix_(rows_dst * block, rows_src * block)]
-    finally:
-        for rdm in pair:
-            rdm.destroy()
-        ref.destroy()
-
-
-def interpolate_within_plex(plex, dm_src, x_src, dm_dst, block, ghosts=True):
-    """Field 0 of ``dm_dst`` interpolated from ``x_src``, the local array of
-    field 0 of ``dm_src``, where both DMs share the topology of ``plex`` (the
-    mesh's DMPlex; a sub-DM is a plain DM in petsc4py, without the plex queries).
-
-    The result is ``dm_dst``'s local array, shaped ``(n, block)`` (``ghosts``
-    True) or its owned part (``ghosts`` False), with each node's value computed
-    by the rank that owns it. Collective on the DMs' communicator.
-    """
-    fe_src = dm_src.getField(0)[0]
-    fe_dst = dm_dst.getField(0)[0]
-    element = _reference_element_matrix(fe_src, fe_dst, plex.getDimension(),
-                                        plex.getCoordinateDim(), plex.isSimplex(), block)
-
-    cS, cE = plex.getHeightStratum(0)
-    cells = np.arange(cS, cE)
-    if element is not None and cells.size:
-        rows_src = closure_dof_rows(dm_src, dm_src.getLocalSection(), cells, block)
-        rows_dst = closure_dof_rows(dm_dst, dm_dst.getLocalSection(), cells, block)
-        uniform = (rows_src.shape[1] == element.shape[1]
-                   and rows_dst.shape[1] == element.shape[0]
-                   and (rows_src >= 0).all() and (rows_dst >= 0).all())
-    else:
-        uniform = element is not None
-    # one verdict for every rank: the two routes below are collective
-    mpi_comm = dm_src.comm.tompi4py()
-    if not mpi_comm.allreduce(bool(uniform), op=MPI.LAND):
-        return _interpolate_by_matrix(plex, dm_src, x_src, dm_dst, block, ghosts)
-
-    lvec = dm_dst.getLocalVec()
-    gvec = dm_dst.getGlobalVec()
-    out = lvec.getArray().reshape(-1, block)
-    if cells.size:
-        x = np.asarray(x_src).reshape(-1, block)
-        out[rows_dst.ravel()] = np.einsum("ab,cbk->cak", element,
-                                          x[rows_src]).reshape(-1, block)
-    # INSERT copies only the points a rank owns, so each node keeps its owner's value
-    dm_dst.localToGlobal(lvec, gvec, addv=PETSc.InsertMode.INSERT_VALUES)
-    if ghosts:
-        dm_dst.globalToLocal(gvec, lvec, addv=PETSc.InsertMode.INSERT_VALUES)
-        result = lvec.getArray().reshape(-1, block).copy()
-    else:
-        result = gvec.getArray().reshape(-1, block).copy()
-    dm_dst.restoreGlobalVec(gvec)
-    dm_dst.restoreLocalVec(lvec)
-    return result
-
-
-def _interpolate_by_matrix(plex, dm_src, x_src, dm_dst, block, ghosts):
-    """PETSc's own interpolation matrix, for the fields whose edge or face nodes
-    are ordered by orientation. It fails on a rank that holds cells but owns no
-    row and no column, so that partition raises here, on every rank, instead."""
-    rows = dm_dst.getGlobalSection().getConstrainedStorageSize()
-    cols = dm_src.getGlobalSection().getConstrainedStorageSize()
-    cS, cE = plex.getHeightStratum(0)
-    mpi_comm = dm_src.comm.tompi4py()
-    if mpi_comm.allreduce(bool(cE > cS and rows == 0 and cols == 0), op=MPI.LOR):
-        raise NotImplementedError(
-            "interpolation between two fields of one mesh: a rank holds cells but "
-            "owns none of their degrees of freedom, where PETSc's interpolation "
-            "matrix cannot be assembled (#902), and this field has more than one "
-            "node on an edge or face. Run on fewer ranks or repartition the mesh.")
-    gsrc = dm_src.getGlobalVec()
-    lsrc = dm_src.getLocalVec()
-    lsrc.getArray()[:] = np.asarray(x_src).ravel()
-    dm_src.localToGlobal(lsrc, gsrc, addv=PETSc.InsertMode.INSERT_VALUES)
-    interp, scale = dm_src.createInterpolation(dm_dst)
-    gdst = dm_dst.getGlobalVec()
-    interp.mult(gsrc, gdst)
-    if ghosts:
-        ldst = dm_dst.getLocalVec()
-        dm_dst.globalToLocal(gdst, ldst, addv=PETSc.InsertMode.INSERT_VALUES)
-        result = ldst.getArray().reshape(-1, block).copy()
-        dm_dst.restoreLocalVec(ldst)
-    else:
-        result = gdst.getArray().reshape(-1, block).copy()
-    dm_dst.restoreGlobalVec(gdst)
-    dm_src.restoreLocalVec(lsrc)
-    dm_src.restoreGlobalVec(gsrc)
-    interp.destroy()
-    if scale is not None:
-        scale.destroy()
-    return result
-
diff --git a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
index 608951139..ff1718fec 100644
--- a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
+++ b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
@@ -2,8 +2,9 @@
 
 PETSc's ``DMPlexComputeInterpolatorNested`` builds the matrix between two DMs on
 one plex (a coordinate DM and a basis DM, a variable's sub-DM and a scratch
-clone). It skips its preallocation pass on a rank with no local rows and no
-local columns, but runs its fill pass over that rank's cells all the same. A
+clone). Since v3.24.0 it skips its preallocation pass on a rank with no local
+rows and no local columns, but runs its fill pass over that rank's cells all
+the same. A
 shared point belongs to the highest rank that holds it, so a rank whose cells
 sit wholly inside higher ranks' regions owns nothing. Its entries then reach
 owners that never preallocated them: the owner raises "New nonzero caused a
@@ -14,8 +15,10 @@
 
 A shell partition gives rank 0 one interior cell of a 4 x 4 grid of triangles
 and the other ranks the rest, so rank 0 owns none of its cell's vertices and
-edges. The test checks that case is set up, then that the node coordinates and
-the projection between degrees come back right on every rank. Before the fix
+edges. The test checks that case is set up, then that the node coordinates
+(degree 1 to 3, the third with two nodes on each edge) and the projection
+between degrees come back right on every rank, through UW3's copy of the
+interpolator (``petsc_dm_create_interpolation_within_plex``). Before the fix
 the degree-2 coordinates hung, so each test carries a thread timeout, which
 fails the run instead of stalling it.
 
@@ -107,23 +110,45 @@ def test_rank_zero_holds_a_cell_and_owns_nothing(mesh):
         sec = mesh.dm.getPointSF()
         leaves = set(np.asarray(sec.getGraph()[1]).tolist())
         owned += sum(1 for p in range(pS, pE) if p not in leaves)
-    if uw.mpi.rank == 0:
-        assert cE - cS == 1 and owned == 0, (cE - cS, owned)
-    assert _all(True)
+    # one verdict for every rank: a rank-0-only assert would leave the others
+    # waiting in the next collective until the timeout
+    assert _all(uw.mpi.rank != 0 or (cE - cS == 1 and owned == 0))
 
 
-@pytest.mark.parametrize("degree", [1, 2])
+def _expected_nodes(verts, degree, t=None):
+    """A triangle's Lagrange nodes: its vertices, then (degree 2) its edge
+    midpoints, or (degree 3) two points on each edge, at ``t`` and ``1 - t``
+    from either end, and the centroid."""
+    a, b = verts[[0, 1, 2]], verts[[1, 2, 0]]
+    if degree == 1:
+        return verts
+    if degree == 2:
+        return np.vstack([verts, 0.5 * (a + b)])
+    return np.vstack([verts, (1 - t) * a + t * b, t * a + (1 - t) * b,
+                      verts.mean(axis=0, keepdims=True)])
+
+
+@pytest.mark.parametrize("degree", [1, 2, 3])
 def test_node_coordinates_on_every_rank(mesh, degree):
-    """Each cell's nodes are its vertices (degree 1), and its vertices and edge
-    midpoints (degree 2)."""
+    """Each cell's nodes are its vertices and, by degree, its edge midpoints,
+    or two points per edge and its centroid. Degree 3 places its edge points by
+    PETSc's node rule (``node_endpoints`` False), so ``t`` is read from the
+    first cell and must then hold on every edge of every cell, on every rank."""
     coords = mesh._get_coords_for_basis(degree, True)
     rows = mesh._cell_node_indices(degree, True)
-    ok = True
-    for k, verts in enumerate(_cell_vertex_coords(mesh)):
-        expected = verts if degree == 1 else np.vstack(
-            [verts, 0.5 * (verts[[0, 1, 2]] + verts[[1, 2, 0]])])
+    cells = _cell_vertex_coords(mesh)
+    t = None
+    if degree == 3 and cells:
+        v, nodes = cells[0], coords[rows[0]]
+        on_edge = [np.dot(x - v[0], v[1] - v[0]) / np.dot(v[1] - v[0], v[1] - v[0])
+                   for x in nodes if abs(np.cross(v[1] - v[0], x - v[0])) < 1e-12]
+        t = min(s for s in on_edge if s > 1e-12)
+    t_all = uw.mpi.comm.allgather(t)
+    ok = all(abs(s - t_all[0]) < 1e-12 for s in t_all if s is not None) if degree == 3 else True
+    for k, verts in enumerate(cells):
         got = coords[rows[k]]
-        ok &= np.allclose(_sorted_rows(got), _sorted_rows(expected), atol=1e-12)
+        ok &= np.allclose(_sorted_rows(got), _sorted_rows(_expected_nodes(verts, degree, t)),
+                          atol=1e-12)
     assert _all(ok)
 
 

From e25585878622e5f66be39dee44f61207ee742c0a Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 14:30:30 +1100
Subject: [PATCH 06/13] Agreed build checks degrade every rank together (#886)

f4b9128a let only a builder failure be retried or degraded in parallel, so the
checks in CustomMGHierarchy.build whose verdict every rank already holds (the
finest map against the operator, DOFs divisible by nodes, one component count,
all on allreduced sizes) now crashed a parallel solve that development had
degraded on every rank together, as serial still does (found in review).

They raise CollectiveBuildError, which TransferBuildError becomes: an error
raised on every rank together, so build_transfers retries or degrades on it in
parallel; a rank-local exception is still raised where it happens. test_1108
gains the case: the finest-map check made to fail on every rank degrades every
rank to the default preconditioner and the solve converges (np 2 and 4).

Underworld development team with AI support from Claude Code
---
 src/underworld3/utilities/custom_mg.py        | 39 +++++++-------
 ..._mg_builder_failure_collective_parallel.py | 52 +++++++++++++++----
 2 files changed, 62 insertions(+), 29 deletions(-)

diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py
index 666f03902..b2cbc5705 100644
--- a/src/underworld3/utilities/custom_mg.py
+++ b/src/underworld3/utilities/custom_mg.py
@@ -224,9 +224,10 @@ def _drop_structural_zeros(P_csr, tol=1e-12):
 _BUILDERS = {"barycentric": barycentric_prolongation, "rbf": rbf_prolongation}
 
 
-class TransferBuildError(RuntimeError):
-    """A node-transfer builder failed on at least one rank; raised on every rank,
-    so the ranks can retry or degrade together (#886)."""
+class CollectiveBuildError(RuntimeError):
+    """A hierarchy-build failure raised on every rank together, from a verdict
+    every rank holds (a builder that failed on any rank, a check on allreduced
+    sizes), so the ranks can retry or degrade together (#886)."""
 
 
 def _require_serial(where):
@@ -913,7 +914,7 @@ def _node_transfer_no_repair(builder, coarse_coords, fine_coords):
 
 def _agreed_node_transfer(builder, coarse_coords, fine_coords, comm):
     """The rank-local node transfer of a parallel build, with one verdict for
-    every rank: if the builder raised on any rank, :class:`TransferBuildError`
+    every rank: if the builder raised on any rank, :class:`CollectiveBuildError`
     is raised on all of them. The builders do no communication, so every rank
     reaches the verdict; a rank that kept its own error would go on to retry
     while the others assemble the transfer, and the two would meet in different
@@ -926,7 +927,7 @@ def _agreed_node_transfer(builder, coarse_coords, fine_coords, comm):
         Pn, problem = None, f"rank {comm.getRank()}: {exc}"
     problems = [p for p in comm.tompi4py().allgather(problem) if p is not None]
     if problems:
-        raise TransferBuildError("; ".join(problems))
+        raise CollectiveBuildError("; ".join(problems))
     return Pn
 
 
@@ -1517,14 +1518,15 @@ def build(self, solver):
                 level_comm = solver.dm.comm.tompi4py()
                 nodes, dofs = level_comm.allreduce(nodes), level_comm.allreduce(dofs)
             nc = dofs // nodes
-            if dofs % nodes != 0:
-                raise RuntimeError(
+            if dofs % nodes != 0:      # allreduced: every rank raises
+                raise CollectiveBuildError(
                     f"level {k}: full DOFs {dofs} not divisible by nodes {nodes}")
             coords.append(c)
             ncomp.append(nc)
 
-        if len(set(ncomp)) != 1:
-            raise RuntimeError(f"inconsistent component counts across levels: {ncomp}")
+        if len(set(ncomp)) != 1:       # from allreduced counts: every rank raises
+            raise CollectiveBuildError(
+                f"inconsistent component counts across levels: {ncomp}")
         nc = ncomp[0]
 
         # Operator-faithful check: the finest reduced map must span exactly the
@@ -1832,8 +1834,8 @@ def _assert_finest_matches_operator(solver, finest_map, parallel):
                 int(finest_map.rend - finest_map.rstart)))
         else:
             red_n = int(len(finest_map))             # r2f length = reduced global size
-        if red_n != op_n:
-            raise RuntimeError(
+        if red_n != op_n:               # allreduced: every rank raises
+            raise CollectiveBuildError(
                 f"custom_mg: finest reduced-map size {red_n} != assembled operator "
                 f"size {op_n}. The DM global section disagrees with the operator "
                 f"(an adapt-child section inconsistency); the finest transfer would "
@@ -1984,8 +1986,8 @@ def build_transfers(solver, field_id=None):
     other two are PREFERENCES, built through the same opportunistic arm
     (barycentric, RBF retry, degrade to the solver's default preconditioner —
     every step recorded in ``solver.pc_fallbacks``). In parallel only a builder
-    failure, which every rank sees (:class:`TransferBuildError`), is retried
-    or degraded; any other failure is raised where it happens (#886).
+    failure every rank raises together (:class:`CollectiveBuildError`) is
+    retried or degraded; any other failure is raised where it happens (#886).
 
     A refinement child carries ``mesh._custom_mg_coarse_meshes`` (the static
     coarse tail), so a :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]``
@@ -2097,11 +2099,12 @@ def build_transfers(solver, field_id=None):
     # Falling back to GAMG loses the hierarchy entirely, so try the cheaper
     # degradation first. See #424.
     # In parallel the retry and the degrade are collective decisions, so they
-    # are taken only on a verdict every rank holds: a builder failure, which
-    # is raised on every rank. Any other exception is one rank's: the others
-    # are inside the build's collectives, and a rank that caught it and
-    # retried would meet them in different ones (#886). It is raised.
-    _retryable = TransferBuildError if solver.dm.comm.getSize() > 1 else Exception
+    # are taken only on a verdict every rank holds (CollectiveBuildError: a
+    # builder failure, a check on allreduced sizes). Any other exception is
+    # one rank's: the others are inside the build's collectives, and a rank
+    # that caught it and retried would meet them in different ones (#886).
+    # It is raised.
+    _retryable = CollectiveBuildError if solver.dm.comm.getSize() > 1 else Exception
     _attempts = [builder] + (["rbf"] if builder != "rbf" else [])
     h = Ps = None
     for _i, _b in enumerate(_attempts):
diff --git a/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
index 2e87f1a87..4143bbcb6 100644
--- a/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
+++ b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
@@ -10,7 +10,9 @@
 after it: a failure on any rank is raised on every rank, and the retry is taken
 together. Here the builder raises on rank 1 only. Contract: every rank records
 the RBF rescue, geometric multigrid is installed, and the solve converges. A
-thread timeout turns the old hang into a failure.
+check on allreduced sizes fails on every rank at once; it degrades every rank
+to the default preconditioner together, as in serial. A thread timeout turns
+the old hang into a failure.
 
     mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
 """
@@ -23,25 +25,32 @@
               pytest.mark.timeout(300, method="thread")]
 
 
-def test_a_builder_failure_on_one_rank_is_retried_on_every_rank(monkeypatch):
-    barycentric = custom_mg._BUILDERS["barycentric"]
-
-    def fails_on_rank_one(coarse_coords, fine_coords):
-        if uw.mpi.rank == 1:
-            raise RuntimeError("synthetic barycentric failure on rank 1 (test probe)")
-        return barycentric(coarse_coords, fine_coords, repair_orphans=False)
-
-    monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", fails_on_rank_one)
+def _fmg_poisson(name):
+    """A single-field Poisson asking for geometric multigrid over the mesh's own
+    refinement hierarchy (the opportunistic, degradable build)."""
     box = uw.meshing.UnstructuredSimplexBox(
         minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, refinement=2,
         regular=True, qdegree=3)
-    T = uw.discretisation.MeshVariable("T1108", box, 1, degree=2)
+    T = uw.discretisation.MeshVariable(name, box, 1, degree=2)
     p = uw.systems.Poisson(box, T)
     p.constitutive_model = uw.constitutive_models.DiffusionModel
     p.constitutive_model.Parameters.diffusivity = 1.0
     p.add_dirichlet_bc(0.0, "Bottom")
     p.add_dirichlet_bc(1.0, "Top")
     p.preconditioner = "fmg"
+    return p
+
+
+def test_a_builder_failure_on_one_rank_is_retried_on_every_rank(monkeypatch):
+    barycentric = custom_mg._BUILDERS["barycentric"]
+
+    def fails_on_rank_one(coarse_coords, fine_coords):
+        if uw.mpi.rank == 1:
+            raise RuntimeError("synthetic barycentric failure on rank 1 (test probe)")
+        return barycentric(coarse_coords, fine_coords, repair_orphans=False)
+
+    monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", fails_on_rank_one)
+    p = _fmg_poisson("T1108")
     with pytest.warns(UserWarning, match="retrying with the 'rbf' builder"):
         p.solve()
 
@@ -50,3 +59,24 @@ def fails_on_rank_one(coarse_coords, fine_coords):
     rec = p.pc_fallbacks["custom_mg.transfer_builder"]
     assert rec["reason"] == "build_failed" and "rank 1" in rec["detail"], rec
     assert "custom_mg.build" not in p.pc_fallbacks
+
+
+def test_a_failure_every_rank_raises_degrades_every_rank(monkeypatch):
+    """The finest map's size is checked against the operator on allreduced sizes,
+    so every rank raises together. Shifting each rank's owned range by one makes
+    the check fail everywhere: the build is retried, then degraded to the
+    default preconditioner on every rank, and the solve converges on it."""
+    check = custom_mg.CustomMGHierarchy._assert_finest_matches_operator
+
+    def off_by_one(solver, finest_map, parallel):
+        return check(solver, finest_map._replace(rend=finest_map.rend + 1), parallel)
+
+    monkeypatch.setattr(custom_mg.CustomMGHierarchy, "_assert_finest_matches_operator",
+                        staticmethod(off_by_one))
+    p = _fmg_poisson("T1108b")
+    with pytest.warns(UserWarning, match="custom-P FMG build failed"):
+        p.solve()
+
+    assert p.snes.getConvergedReason() > 0
+    assert p.snes.getKSP().getPC().getType() == "gamg"
+    assert p.pc_fallbacks["custom_mg.build"]["reason"] == "build_failed"

From 96d39f7b32d4ee4995221486add4c138fe767503 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 16:57:23 +1100
Subject: [PATCH 07/13] Mark the summed discontinuous-to-continuous projection
 in project_to_degree (#905)

A discontinuous source projected to a continuous target sums the adjacent
cells at each shared node. It was already on development and is outside #902,
so this branch only marks the place.

Underworld development team with AI support from Claude Code
---
 src/underworld3/function/field_projection.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py
index 0c6e95ded..61a3456d1 100644
--- a/src/underworld3/function/field_projection.py
+++ b/src/underworld3/function/field_projection.py
@@ -89,6 +89,12 @@ def project_to_degree(
     )
     interp = petsc_dm_create_interpolation_within_plex(subdm_src, dm_scratch)
 
+    # TODO(BUG): a discontinuous source projected to a continuous target sums
+    # the adjacent cells' values at each shared node (2-7x at P1 vertices): the
+    # matrix inserts one entry per cell into different source columns. The
+    # XDMF writer is not affected (it projects a discontinuous field to DG0).
+    # See issue #905.
+
     # --- apply ---
     g_src = subdm_src.createGlobalVec()
     subdm_src.localToGlobal(mesh_var._lvec, g_src)

From fcef03fec5eb0c6ca42b53e03481628db396dc88 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 17:07:41 +1100
Subject: [PATCH 08/13] test_1107: custom-P FMG over a coarse level where a
 rank owns nothing (#902)

The coarse level is the shell-partitioned grid on which rank 0 holds one
cell and owns none of its DOFs, under the grid's refinement. The build reads
each level's node coordinates, which is where kaiju's np 32/52 runs failed.
With PETSc's createInterpolation routed back in, the test fails inside
CustomMGHierarchy.build with "New nonzero at (88,32) caused a malloc" at
np 2. On the branch it installs multigrid and agrees with GAMG to 1e-6, at
np 2 and 4.

Underworld development team with AI support from Claude Code
---
 ...erpolation_rank_owning_nothing_parallel.py | 81 +++++++++++++++----
 1 file changed, 67 insertions(+), 14 deletions(-)

diff --git a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
index ff1718fec..165a82be1 100644
--- a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
+++ b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
@@ -18,12 +18,16 @@
 edges. The test checks that case is set up, then that the node coordinates
 (degree 1 to 3, the third with two nodes on each edge) and the projection
 between degrees come back right on every rank, through UW3's copy of the
-interpolator (``petsc_dm_create_interpolation_within_plex``). Before the fix
-the degree-2 coordinates hung, so each test carries a thread timeout, which
-fails the run instead of stalling it.
+interpolator (``petsc_dm_create_interpolation_within_plex``). A custom-P FMG
+solve whose coarse level is that grid, under its refinement, pins the same
+case inside a real multigrid build. Before the fix the degree-2 coordinates
+hung, so each test carries a thread timeout, which fails the run instead of
+stalling it.
 
     mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
 """
+from enum import Enum
+
 import numpy as np
 import pytest
 from mpi4py import MPI
@@ -31,16 +35,19 @@
 
 import underworld3 as uw
 from underworld3.function import project_to_degree
+from underworld3.utilities import custom_mg
 
 pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, pytest.mark.tier_b,
               pytest.mark.timeout(120, method="thread")]
 
 N = 4
+_WALLS = Enum("Boundaries", {"Walls": 1})
 
 
-def _owned_nothing_mesh():
+def _owned_nothing_dm():
     """The 4 x 4 grid of triangles, all on rank 0, partitioned so rank 0 keeps one
-    interior cell (every edge shared with a cell of another rank)."""
+    interior cell (every edge shared with a cell of another rank). Its outer
+    edges carry the label ``Walls``."""
     comm = PETSc.COMM_WORLD
     rank, size = comm.getRank(), comm.getSize()
     if rank == 0:
@@ -56,6 +63,9 @@ def _owned_nothing_mesh():
         coords = np.zeros((0, 2))
         cells = np.zeros((0, 3), dtype=PETSc.IntType)
     dm = PETSc.DMPlex().createFromCellList(2, cells, coords, interpolate=True, comm=comm)
+    dm.createLabel("Walls")
+    dm.markBoundaryFaces("Walls", 1)
+    dm.labelComplete(dm.getLabel("Walls"))
     part = dm.getPartitioner()
     part.setType(PETSc.Partitioner.Type.SHELL)
     if rank == 0:
@@ -69,12 +79,12 @@ def _owned_nothing_mesh():
     else:
         part.setShellPartition(size, [0] * size, [])
     dm.distribute()
-    return uw.discretisation.Mesh(dm, degree=1, qdegree=2)
+    return dm
 
 
 @pytest.fixture(scope="module")
 def mesh():
-    return _owned_nothing_mesh()
+    return uw.discretisation.Mesh(_owned_nothing_dm(), degree=1, qdegree=2)
 
 
 def _all(flag):
@@ -101,18 +111,22 @@ def _sorted_rows(a):
     return a[np.lexsort(a.T[::-1])]
 
 
-def test_rank_zero_holds_a_cell_and_owns_nothing(mesh):
-    """The case under test: rank 0 holds one cell and owns no vertex or edge."""
+def _rank_zero_owns_nothing(mesh):
+    """Rank 0 holds one cell and owns no vertex or edge, as one verdict for every
+    rank: a rank-0-only assert would leave the others waiting in the next
+    collective until the timeout."""
     cS, cE = mesh.dm.getHeightStratum(0)
+    leaves = set(np.asarray(mesh.dm.getPointSF().getGraph()[1]).tolist())
     owned = 0
     for d in (0, 1):
         pS, pE = mesh.dm.getDepthStratum(d)
-        sec = mesh.dm.getPointSF()
-        leaves = set(np.asarray(sec.getGraph()[1]).tolist())
         owned += sum(1 for p in range(pS, pE) if p not in leaves)
-    # one verdict for every rank: a rank-0-only assert would leave the others
-    # waiting in the next collective until the timeout
-    assert _all(uw.mpi.rank != 0 or (cE - cS == 1 and owned == 0))
+    return _all(uw.mpi.rank != 0 or (cE - cS == 1 and owned == 0))
+
+
+def test_rank_zero_holds_a_cell_and_owns_nothing(mesh):
+    """The case under test."""
+    assert _rank_zero_owns_nothing(mesh)
 
 
 def _expected_nodes(verts, degree, t=None):
@@ -169,3 +183,42 @@ def test_projection_between_degrees_on_every_rank(mesh):
     ok = (np.allclose(up[:, 0], n2[:, 0] + 2.0 * n2[:, 1], atol=1e-12)
           and np.allclose(down[:, 0], n1[:, 0] ** 2, atol=1e-12))
     assert _all(ok)
+
+
+def _poisson_on(dm, name):
+    """A degree-2 Poisson problem, zero on the walls with a unit source."""
+    mesh = uw.discretisation.Mesh(dm.clone(), degree=1, qdegree=3, boundaries=_WALLS)
+    T = uw.discretisation.MeshVariable(name, mesh, 1, degree=2)
+    p = uw.systems.Poisson(mesh, T)
+    p.constitutive_model = uw.constitutive_models.DiffusionModel
+    p.constitutive_model.Parameters.diffusivity = 1.0
+    p.f = 1.0
+    p.add_dirichlet_bc(0.0, "Walls")
+    p.petsc_options["ksp_rtol"] = 1.0e-10
+    return p
+
+
+def test_custom_fmg_over_a_level_where_a_rank_owns_nothing():
+    """Custom-P FMG reads every level's node coordinates. Its coarse level is the
+    grid on which rank 0 owns nothing, under the grid's refinement. The build
+    installs multigrid, and the solve agrees with GAMG's."""
+    base = _owned_nothing_dm()
+    fine = base.refine()
+    coarse = uw.discretisation.Mesh(base.clone(), degree=1, qdegree=3, boundaries=_WALLS)
+    assert _rank_zero_owns_nothing(coarse)                      # the case under test
+
+    reference = _poisson_on(fine, "Tg_1107")
+    reference.preconditioner = "gamg"
+    reference.solve()
+
+    solver = _poisson_on(fine, "Tm_1107")
+    custom_mg.set_custom_fmg(solver, [coarse], builder="barycentric")
+    solver.solve()
+
+    assert solver.snes.getKSP().getPC().getType() == "mg"
+    assert solver.snes.getConvergedReason() > 0
+    assert "custom_mg.build" not in solver.pc_fallbacks
+    u, r = np.asarray(solver.Unknowns.u.data), np.asarray(reference.Unknowns.u.data)
+    num = uw.mpi.comm.allreduce(float(np.sum((u - r) ** 2)))
+    den = uw.mpi.comm.allreduce(float(np.sum(r ** 2)))
+    assert (num / den) ** 0.5 < 1.0e-6

From 7b39d2c61900eafe765b9fe7c40f7222226b3d5f Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 18:19:30 +1100
Subject: [PATCH 09/13] #886: fault_mg's agreed checks degrade a parallel
 build; the split replay agrees on any exception

Review round 2 found two gaps in the narrowed parallel catch of
build_transfers.

- provenance_transfer_parallel's two checks on allreduced counts ("a fine
  node belongs to no cell", "fine node(s) lie in no cell") raised a plain
  RuntimeError on every rank. In parallel only CollectiveBuildError is
  retried or degraded, so the solve crashed where development degraded it.
  They raise CollectiveBuildError now.
- split_tail_levels caught only ValueError, RuntimeError and
  NotImplementedError before its collective verdict, so an IndexError on
  one rank after the split's last collective escaped the verdict. It
  catches any exception now. test_1089 is parametrised over IndexError,
  which failed on the previous build ("DID NOT WARN", then a hang) and
  passes at np 2 and 3.

rotated-freeslip.md says which failures degrade in parallel and which
are raised.

Underworld development team with AI support from Claude Code
---
 docs/developer/subsystems/rotated-freeslip.md   |  7 ++++++-
 src/underworld3/utilities/fault_mg.py           | 17 +++++++++++------
 .../test_1089_fault_split_multigrid_parallel.py | 11 +++++++----
 3 files changed, 24 insertions(+), 11 deletions(-)

diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md
index 4354917c4..129bc7a74 100644
--- a/docs/developer/subsystems/rotated-freeslip.md
+++ b/docs/developer/subsystems/rotated-freeslip.md
@@ -181,7 +181,12 @@ path uses**:
 2. a **mesh-owned** coarse tail (`mesh._custom_mg_coarse_meshes`, what
    `mesh.adapt()` leaves on a refinement child) is picked up opportunistically,
    with the barycentric→RBF builder fallback, so a failed build degrades to GAMG
-   rather than crashing the solve.
+   rather than crashing the solve. In parallel the retry and the degrade happen
+   on a verdict every rank holds: a builder failure, which is agreed across
+   ranks as soon as the builder returns, or a check on allreduced sizes
+   (`custom_mg.CollectiveBuildError`). Any other exception is one rank's while
+   the others are inside the build's collectives, so it is raised rather than
+   retried alone (#886).
 
 Step 2 used to be unreachable from here. The standard path's injection hook runs
 *after* the rotated dispatch has already returned, and the rotated builder
diff --git a/src/underworld3/utilities/fault_mg.py b/src/underworld3/utilities/fault_mg.py
index a2f1b7008..6cf36f61c 100644
--- a/src/underworld3/utilities/fault_mg.py
+++ b/src/underworld3/utilities/fault_mg.py
@@ -461,9 +461,12 @@ def provenance_transfer_parallel(coarse, fine, pre, degree, continuous,
     nv = cv_c.shape[1]
     dim = Xc.shape[1]
 
+    # both checks below are on allreduced counts, so every rank raises
+    # together: build_transfers may retry or degrade on them (#886)
+    from underworld3.utilities.custom_mg import CollectiveBuildError
     owner = _fine_cell_of_node(fn, n_f)
     if m4.allreduce(int((owner < 0).sum())):
-        raise RuntimeError("fault_mg: a fine node belongs to no cell")
+        raise CollectiveBuildError("fault_mg: a fine node belongs to no cell")
     mine = np.flatnonzero(((l2g_f >= fstart) & (l2g_f < fend)).any(axis=1))
     fcent = vxy_f[cv_f].mean(axis=1)
     if pre is fine:
@@ -478,7 +481,7 @@ def provenance_transfer_parallel(coarse, fine, pre, degree, continuous,
     n_out = int(m4.allreduce(int((plam.min(axis=1) < -_INSIDE_TOL).sum())
                              if mine.size else 0))
     if n_out:
-        raise RuntimeError(
+        raise CollectiveBuildError(
             f"fault_mg: {n_out} fine node(s) lie in no cell of the level as the "
             f"add-on partition left it; the surgery was not rank-local.")
     parent = np.asarray(pre._addon_parent, dtype=np.int32)[pcell]   # (nm, 2)
@@ -612,12 +615,14 @@ def split_tail_levels(levels, fault_specs, verbose=False):
             # partition ran once, for every level): no second partition
             split = add_fault(level, entries, verbose=verbose, partition=False,
                               **kwargs)
-        except (ValueError, RuntimeError, NotImplementedError) as exc:
-            problem = str(exc)
+        except Exception as exc:   # sanctioned: agreed below, the level stays welded
+            problem = f"{type(exc).__name__}: {exc}"
         # the split's refusals are collective; a failure on one rank after its
         # last collective is not, and every rank must stop at the same level or
-        # the ranks hold different hierarchies (a rank-local raise inside the
-        # split, before its own collectives, still hangs there: issue #874)
+        # the ranks hold different hierarchies. Any exception is caught, so an
+        # IndexError on one rank meets the same verdict (#886); a rank-local
+        # raise inside the split, before its own collectives, still hangs
+        # there (issue #874)
         try:
             _collective_raise(problem)
         except RuntimeError as exc:
diff --git a/tests/parallel/test_1089_fault_split_multigrid_parallel.py b/tests/parallel/test_1089_fault_split_multigrid_parallel.py
index 2fd5489c7..4ed9792be 100644
--- a/tests/parallel/test_1089_fault_split_multigrid_parallel.py
+++ b/tests/parallel/test_1089_fault_split_multigrid_parallel.py
@@ -94,14 +94,17 @@ def centroid(q):
     assert np.abs(np.einsum("ij,ij->i", jumps, normals)).max() < 1.0e-8
 
 
-def test_a_failure_on_one_rank_stops_the_replay_on_every_rank(monkeypatch):
+@pytest.mark.parametrize("error", [RuntimeError, IndexError])
+def test_a_failure_on_one_rank_stops_the_replay_on_every_rank(monkeypatch, error):
     """A failure of the split on one rank after its last collective (a final
     consistency check, rank-local) stopped the replay of the fault on that rank
     alone: the ranks kept different split levels and the next collective on the
     hierarchy met mismatched partners. The verdict is exchanged per level, so
     every rank keeps the same split levels and the levels from the failure down
-    stay welded. (A rank-local raise BEFORE the split's own collectives still
-    hangs inside the split; issue #874.)"""
+    stay welded. Any exception meets the verdict, not only the split's own
+    refusal types: an IndexError escaped it on one rank (#886). (A rank-local
+    raise BEFORE the split's own collectives still hangs inside the split;
+    issue #874.)"""
     import underworld3.utilities.fault_split as fault_split
     original, calls = fault_split.add_fault, [0]
 
@@ -109,7 +112,7 @@ def fails_on_one_rank_after_the_second_add_on_split(mesh, *args, **kwargs):
         calls[0] += 1                 # 1: the finest; 2, 3: the add-on levels
         split = original(mesh, *args, **kwargs)
         if calls[0] == 3 and uw.mpi.rank == 1:
-            raise RuntimeError("an internal failure on one rank")
+            raise error("an internal failure on one rank")
         return split
 
     monkeypatch.setattr(fault_split, "add_fault",

From 57c0de90ed7dbd4f8aa62b837c64be5f4b2c2f80 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 18:19:47 +1100
Subject: [PATCH 10/13] #902: the mesh's coordinate space is built through
 UW3's interpolator; PETSc's notice, a one-plex check, sharper tests

Review round 2.

- Mesh(degree >= 2) on a rank that holds cells but owns none of their
  nodes failed during mesh setup. createCoordinateSpace projects the
  coordinates through PETSc's DMCreateInterpolation between the old and
  the new coordinate DM. UW_DMPlexCreateCoordinateSpaceWithinPlex installs
  the space without projecting, then copies the coordinates (same dual
  space, as PETSc does) or interpolates them with
  UW_DMCreateInterpolationWithinPlex. Coordinates injected into a higher
  coordinate dimension are not projected, as in PETSc. Periodic
  coordinates keep PETSc's own projection, which is local for isoperiodic
  ones. The TODO(BUG) at that call is gone.
- The copy of PETSc's interpolator carries PETSc's full BSD-2 notice:
  the copyright line, the conditions and the disclaimer.
- The copy refuses, on every rank, two DMs that are not on one plex.
  PETSc only dispatches there for that case. A pair on two plexes would
  otherwise give a wrong matrix without an error.
- test_1107 runs every case on two partitions.
  - Rank 0 holds an interior cell, owning no vertex or edge: PETSc's P2
    matrices fail there.
  - Rank 0 holds a boundary cell, owning its wall edge and no vertex:
    PETSc's P1 matrices fail there too, since the wall edge's vertex pair
    is in no other cell.
  - A degree-2 mesh on each partition checks that the vertices stay on
    the grid and each edge coordinate is its midpoint.
  - The FMG case asserts its 4 iterations, the same at np 2, 3 and 4
    (GAMG takes 12). Coarse nodes displaced by 0.3 of a cell take 5.
  - It reads .array, not .data.
- On PETSc's route, run one test at a time at np 2, these fail: the
  interior P2 coordinates, the projection, the degree-2 mesh and the
  FMG case; and the boundary P1 coordinates and projection. Degree 3
  passes there, since rank 0 owns its cell's interior node, so those
  cases check the copy's node placement. On the branch all 13 pass at
  np 2, 3 and 4.
- The test_1107 and test_1108 docstrings no longer claim the thread
  timeout stops a hang inside PETSc. petsc4py holds the GIL there, and
  the MPI supervisor's silence limit is what ends it.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h         | 108 +++++++++++++++-
 .../cython/petsc_discretisation.pyx           |  16 +++
 src/underworld3/cython/petsc_extras.pxi       |   1 +
 .../discretisation/discretisation_mesh.py     |  14 ++-
 ...erpolation_rank_owning_nothing_parallel.py | 119 +++++++++++++-----
 ..._mg_builder_failure_collective_parallel.py |   6 +-
 6 files changed, 223 insertions(+), 41 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 5eb5993d6..65106f44a 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -370,8 +370,35 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, DM *sandbox)
 // Interpolation between two DMs on ONE plex (a coordinate DM and a basis DM
 // cloned from it, a variable's sub-DM and a scratch clone of the mesh): the
 // matrix DMCreateInterpolation builds through DMPlexComputeInterpolatorNested
-// (isRefined = PETSC_FALSE), copied from PETSc v3.25.0 (BSD-2-Clause) with one
-// change (#902). PETSc runs its preallocation loop only on a rank with local
+// (isRefined = PETSC_FALSE), copied from PETSc v3.25.0 with one change (#902),
+// under PETSc's licence:
+//
+//   Copyright (c) 1991-2025, UChicago Argonne, LLC and the PETSc Developers and
+//   Contributors
+//   All rights reserved.
+//
+//   Redistribution and use in source and binary forms, with or without
+//   modification, are permitted provided that the following conditions are met:
+//
+//   * Redistributions of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//   * Redistributions in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+//   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+//   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+//   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+//   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+//   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+//   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+//   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+//   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+//   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+//   POSSIBILITY OF SUCH DAMAGE.
+//
+// The change: PETSc runs its preallocation loop only on a rank with local
 // rows or columns (`if (locRows || locCols)`, since b4ec2976050, v3.24.0) but
 // its fill loop on every rank that holds cells. A shared point belongs to the
 // highest rank holding it, so a rank whose cells lie inside higher ranks'
@@ -379,8 +406,10 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, DM *sandbox)
 // that never preallocated them ("New nonzero caused a malloc" in
 // MatAssemblyEnd) and the ranks that did not raise wait in the assembly's
 // collective. Here both loops run over the same cells. Finite-element fields
-// only, which are the only ones UW3 interpolates this way. Drop this for
-// DMCreateInterpolation once PETSc runs the two loops alike.
+// only, which are the only ones UW3 interpolates this way, and one plex only:
+// PETSc dispatches here on the two DMs sharing their plex, and a pair on two
+// plexes is refused on every rank. Drop this for DMCreateInterpolation once
+// PETSc runs the two loops alike.
 PetscErrorCode UW_DMCreateInterpolationWithinPlex(DM dmc, DM dmf, Mat *In)
 {
     PetscSection fsection, fglobalSection, csection, cglobalSection;
@@ -393,6 +422,18 @@ PetscErrorCode UW_DMCreateInterpolationWithinPlex(DM dmc, DM dmf, Mat *In)
     Mat          preallocator;
 
     PetscFunctionBeginUser;
+    {
+        // clones share their plex, cone array included; the verdict is
+        // agreed, since a rank holding no points cannot tell two plexes apart
+        PetscInt *conesC, *conesF, other;
+
+        PetscCall(DMPlexGetCones(dmc, &conesC));
+        PetscCall(DMPlexGetCones(dmf, &conesF));
+        other = conesC != conesF;
+        PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &other, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dmc)));
+        PetscCheck(!other, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_INCOMP,
+                   "UW_DMCreateInterpolationWithinPlex: the two DMs are not on one plex");
+    }
     PetscCall(DMGetGlobalSection(dmf, &fglobalSection));
     PetscCall(PetscSectionGetConstrainedStorageSize(fglobalSection, &m));
     PetscCall(DMGetGlobalSection(dmc, &cglobalSection));
@@ -507,3 +548,62 @@ PetscErrorCode UW_DMCreateInterpolationWithinPlex(DM dmc, DM dmf, Mat *In)
     PetscCall(MatViewFromOptions(*In, NULL, "-interp_mat_view"));
     PetscFunctionReturn(PETSC_SUCCESS);
 }
+
+// The coordinate space of a plex at `degree`, as
+// DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_TRUE) builds it,
+// with the coordinates carried across by UW_DMCreateInterpolationWithinPlex.
+// PETSc's projection (DMSetCoordinateDisc) interpolates through
+// DMCreateInterpolation between the old and the new coordinate DM, two DMs on
+// this plex, so a degree >= 2 space could not be built on a rank that holds
+// cells but owns none of their nodes (#902). As in PETSc: the same space is a
+// copy; coordinates injected into a higher coordinate dimension are not
+// projected; periodic coordinates keep PETSc's own projection, which projects
+// isoperiodic ones locally, without the matrix.
+PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(DM dm, PetscInt degree)
+{
+    const PetscReal *maxCell, *Lstart, *L;
+    PetscInt         nfaceSF = 0;
+    DM               cdmOld, cdmNew;
+    Vec              xOld, xNew, x;
+    PetscFE          discOld, disc;
+    PetscDualSpace   dsOld, ds;
+    PetscBool        same;
+
+    PetscFunctionBeginUser;
+    PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
+    PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &nfaceSF, NULL));
+    if (maxCell || L || nfaceSF) {
+        PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_TRUE));
+        PetscFunctionReturn(PETSC_SUCCESS);
+    }
+    PetscCall(DMGetCoordinateDM(dm, &cdmOld));
+    PetscCall(PetscObjectReference((PetscObject)cdmOld));
+    PetscCall(DMGetCoordinates(dm, &xOld));
+    PetscCall(PetscObjectReference((PetscObject)xOld));
+    // installs the new coordinate DM and leaves the coordinates to us
+    PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_FALSE));
+    PetscCall(DMGetCoordinates(dm, &x));
+    if (x == xOld) {
+        PetscCall(DMGetCoordinateDM(dm, &cdmNew));
+        PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
+        PetscCall(DMGetField(cdmNew, 0, NULL, (PetscObject *)&disc));
+        PetscCall(PetscFEGetDualSpace(discOld, &dsOld));
+        PetscCall(PetscFEGetDualSpace(disc, &ds));
+        PetscCall(PetscDualSpaceEqual(dsOld, ds, &same));
+        PetscCall(DMCreateGlobalVector(cdmNew, &xNew));
+        if (same) {
+            PetscCall(VecCopy(xOld, xNew));
+        } else {
+            Mat In;
+
+            PetscCall(UW_DMCreateInterpolationWithinPlex(cdmOld, cdmNew, &In));
+            PetscCall(MatMult(In, xOld, xNew));
+            PetscCall(MatDestroy(&In));
+        }
+        PetscCall(DMSetCoordinates(dm, xNew));
+        PetscCall(VecDestroy(&xNew));
+    }
+    PetscCall(VecDestroy(&xOld));
+    PetscCall(DMDestroy(&cdmOld));
+    PetscFunctionReturn(PETSC_SUCCESS);
+}
diff --git a/src/underworld3/cython/petsc_discretisation.pyx b/src/underworld3/cython/petsc_discretisation.pyx
index ffa52d935..fb32a5c40 100644
--- a/src/underworld3/cython/petsc_discretisation.pyx
+++ b/src/underworld3/cython/petsc_discretisation.pyx
@@ -199,6 +199,22 @@ def petsc_dm_create_interpolation_within_plex(dmc, dmf):
         return interp
 
 
+def petsc_dm_create_coordinate_space_within_plex(dm, degree):
+        """Give ``dm`` a Lagrange coordinate space of ``degree``, its coordinates
+        carried across from the current space.
+
+        It is ``dm.createCoordinateSpace(degree, False, True)``, with the
+        coordinates interpolated by UW3's copy of PETSc's interpolator, so a
+        degree >= 2 space is also built on a rank that holds cells but owns
+        none of their nodes (#902; see
+        ``UW_DMPlexCreateCoordinateSpaceWithinPlex`` in ``petsc_compat.h``).
+        Collective.
+        """
+        cdef DM c_dm = dm
+        cdef PetscInt c_degree = degree
+        CHKERRQ( UW_DMPlexCreateCoordinateSpaceWithinPlex(c_dm.dm, c_degree) )
+
+
 
 # This is not cython, does it need to be here or in discretisation.py ?
 
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 4e39315b7..337f157e3 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -59,6 +59,7 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*)
     PetscErrorCode UW_DMCreateBdIntegralSandbox(PetscDM, PetscDM*)
     PetscErrorCode UW_DMCreateInterpolationWithinPlex(PetscDM, PetscDM, PetscMat*)
+    PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(PetscDM, PetscInt)
 
 cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[])
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 25fc70028..c316dd579 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -2754,12 +2754,14 @@ def nuke_coords_and_rebuild(
             # dual space without proper point subspaces, causing
             # DMPlexComputeBdIntegral to segfault/deadlock (issue #96).
             # createCoordinateSpace builds the FE internally with correct
-            # subspace initialisation.
-            # TODO(BUG): #902 — for a degree >= 2 coordinate field PETSc projects
-            # the coordinates through an interpolation matrix between two DMs on
-            # this plex, which cannot be assembled on a rank that holds cells but
-            # owns none of their nodes (the job hangs). Degree 1 is a copy.
-            self.dm.createCoordinateSpace(self.degree, False, True)
+            # subspace initialisation. The coordinates are carried across by
+            # UW3's copy of PETSc's interpolator: PETSc's own projection of a
+            # degree >= 2 space cannot be assembled on a rank that holds cells
+            # but owns none of their nodes (#902).
+            from underworld3.cython.petsc_discretisation import (
+                petsc_dm_create_coordinate_space_within_plex,
+            )
+            petsc_dm_create_coordinate_space_within_plex(self.dm, self.degree)
 
             # Issue #96 fix: Force coordinate field creation and strip
             # boundary labels from the coordinate DM. createCoordinateSpace
diff --git a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
index 165a82be1..9ffe7d5a0 100644
--- a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
+++ b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
@@ -13,16 +13,26 @@
 124-cell coarsest level of the faulted annulus, where custom-P multigrid reads
 each level's node coordinates (#902), and on a 3-D notch mesh at np 16.
 
-A shell partition gives rank 0 one interior cell of a 4 x 4 grid of triangles
-and the other ranks the rest, so rank 0 owns none of its cell's vertices and
-edges. The test checks that case is set up, then that the node coordinates
-(degree 1 to 3, the third with two nodes on each edge) and the projection
-between degrees come back right on every rank, through UW3's copy of the
-interpolator (``petsc_dm_create_interpolation_within_plex``). A custom-P FMG
-solve whose coarse level is that grid, under its refinement, pins the same
-case inside a real multigrid build. Before the fix the degree-2 coordinates
-hung, so each test carries a thread timeout, which fails the run instead of
-stalling it.
+A shell partition gives rank 0 one cell of a 4 x 4 grid of triangles and the
+other ranks the rest. An interior cell leaves rank 0 owning none of its
+vertices and edges: PETSc's P2 matrices fail there, while its P1 matrices
+escape, because every vertex pair of the cell is also an edge of another
+rank's cell. A boundary cell leaves rank 0 owning its wall edge but none of
+its vertices: there PETSc's P1 matrices fail, since the wall edge's vertex
+pair occurs in no other cell, and its P2 matrices pass. At degree 3 rank 0
+owns its cell's interior node, so PETSc preallocates; those cases check the
+copy's node placement. The tests check each case is set up, then
+that the node coordinates (degree 1 to 3, the third with two nodes on each
+edge) and the projection between degrees come back right on every rank,
+through UW3's copy of the interpolator
+(``petsc_dm_create_interpolation_within_plex``). A degree-2 mesh on the
+partition pins the coordinate space PETSc projects in ``Mesh`` setup, and a
+custom-P FMG solve whose coarse level is the grid, under its refinement, pins
+the case inside a real multigrid build.
+
+Before the fix the ranks failed or hung inside PETSc. The thread timeout below
+stops a rank waiting in Python; a rank held inside a PETSc collective keeps the
+GIL, so a hang there is ended by the MPI supervisor's silence limit instead.
 
     mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py
 """
@@ -41,13 +51,21 @@
               pytest.mark.timeout(120, method="thread")]
 
 N = 4
+# Custom-P FMG iterations of the test below, at np 2, 3 and 4 (GAMG takes 12). A
+# transfer built from coarse nodes displaced by 0.3 of a cell takes 5.
+FMG_ITERATIONS = 4
 _WALLS = Enum("Boundaries", {"Walls": 1})
 
 
-def _owned_nothing_dm():
+def _cells_of(dm, v):
+    return {c for e in dm.getSupport(v) for c in dm.getSupport(e)}
+
+
+def _owned_nothing_dm(where="interior"):
     """The 4 x 4 grid of triangles, all on rank 0, partitioned so rank 0 keeps one
-    interior cell (every edge shared with a cell of another rank). Its outer
-    edges carry the label ``Walls``."""
+    cell: an ``interior`` one (every edge shared with a cell of another rank),
+    or a ``boundary`` one with a single wall edge and every vertex also in
+    another cell. Its outer edges carry the label ``Walls``."""
     comm = PETSc.COMM_WORLD
     rank, size = comm.getRank(), comm.getSize()
     if rank == 0:
@@ -70,9 +88,20 @@ def _owned_nothing_dm():
     part.setType(PETSc.Partitioner.Type.SHELL)
     if rank == 0:
         cS, cE = dm.getHeightStratum(0)
-        interior = [c for c in range(cS, cE)
-                    if all(dm.getSupportSize(e) == 2 for e in dm.getCone(c))]
-        c0 = interior[len(interior) // 2]
+        vS, vE = dm.getDepthStratum(0)
+
+        def walls(c):
+            return sum(dm.getSupportSize(e) == 1 for e in dm.getCone(c))
+
+        def shared_vertices(c):
+            return all(len(_cells_of(dm, v)) > 1
+                       for v in dm.getTransitiveClosure(c)[0] if vS <= v < vE)
+
+        if where == "interior":
+            chosen = [c for c in range(cS, cE) if walls(c) == 0]
+        else:
+            chosen = [c for c in range(cS, cE) if walls(c) == 1 and shared_vertices(c)]
+        c0 = chosen[len(chosen) // 2]
         rest = np.array_split(np.asarray([c for c in range(cS, cE) if c != c0]), size - 1)
         part.setShellPartition(size, [1] + [len(r) for r in rest],
                                [c0] + [int(c) for r in rest for c in r])
@@ -82,9 +111,11 @@ def _owned_nothing_dm():
     return dm
 
 
-@pytest.fixture(scope="module")
-def mesh():
-    return uw.discretisation.Mesh(_owned_nothing_dm(), degree=1, qdegree=2)
+@pytest.fixture(scope="module", params=["interior", "boundary"])
+def mesh(request):
+    mesh = uw.discretisation.Mesh(_owned_nothing_dm(request.param), degree=1, qdegree=2)
+    mesh._case = request.param
+    return mesh
 
 
 def _all(flag):
@@ -111,22 +142,23 @@ def _sorted_rows(a):
     return a[np.lexsort(a.T[::-1])]
 
 
-def _rank_zero_owns_nothing(mesh):
-    """Rank 0 holds one cell and owns no vertex or edge, as one verdict for every
-    rank: a rank-0-only assert would leave the others waiting in the next
-    collective until the timeout."""
+def _rank_zero_owns_nothing(mesh, where="interior"):
+    """Rank 0 holds one cell and owns no vertex, nor (interior cell) any edge, as
+    one verdict for every rank: a rank-0-only assert would leave the others
+    waiting in the next collective until the timeout."""
     cS, cE = mesh.dm.getHeightStratum(0)
     leaves = set(np.asarray(mesh.dm.getPointSF().getGraph()[1]).tolist())
-    owned = 0
+    owned = []
     for d in (0, 1):
         pS, pE = mesh.dm.getDepthStratum(d)
-        owned += sum(1 for p in range(pS, pE) if p not in leaves)
-    return _all(uw.mpi.rank != 0 or (cE - cS == 1 and owned == 0))
+        owned.append(sum(1 for p in range(pS, pE) if p not in leaves))
+    expected = [0, 0] if where == "interior" else [0, 1]
+    return _all(uw.mpi.rank != 0 or (cE - cS == 1 and owned == expected))
 
 
 def test_rank_zero_holds_a_cell_and_owns_nothing(mesh):
     """The case under test."""
-    assert _rank_zero_owns_nothing(mesh)
+    assert _rank_zero_owns_nothing(mesh, mesh._case)
 
 
 def _expected_nodes(verts, degree, t=None):
@@ -185,6 +217,31 @@ def test_projection_between_degrees_on_every_rank(mesh):
     assert _all(ok)
 
 
+@pytest.mark.parametrize("where", ["interior", "boundary"])
+def test_a_degree_two_mesh_builds_on_every_rank(where):
+    """``Mesh(degree=2)`` gives the plex a degree-2 coordinate space, carried
+    across from the vertex coordinates: the vertices stay on the grid and every
+    edge's coordinate is its midpoint."""
+    mesh = uw.discretisation.Mesh(_owned_nothing_dm(where), degree=2, qdegree=3)
+    assert _rank_zero_owns_nothing(mesh, where)                 # the case under test
+    dm = mesh.dm
+    sec = dm.getCoordinateSection()
+    x = dm.getCoordinatesLocal().array
+
+    def at(p):
+        return x[sec.getOffset(p):sec.getOffset(p) + 2]
+
+    ok = True
+    vS, vE = dm.getDepthStratum(0)
+    for v in range(vS, vE):
+        ok &= sec.getDof(v) == 2 and np.allclose(at(v) * N, np.round(at(v) * N), atol=1e-12)
+    eS, eE = dm.getDepthStratum(1)
+    for e in range(eS, eE):
+        a, b = dm.getCone(e)
+        ok &= sec.getDof(e) == 2 and np.allclose(at(e), 0.5 * (at(a) + at(b)), atol=1e-12)
+    assert _all(ok)
+
+
 def _poisson_on(dm, name):
     """A degree-2 Poisson problem, zero on the walls with a unit source."""
     mesh = uw.discretisation.Mesh(dm.clone(), degree=1, qdegree=3, boundaries=_WALLS)
@@ -201,7 +258,10 @@ def _poisson_on(dm, name):
 def test_custom_fmg_over_a_level_where_a_rank_owns_nothing():
     """Custom-P FMG reads every level's node coordinates. Its coarse level is the
     grid on which rank 0 owns nothing, under the grid's refinement. The build
-    installs multigrid, and the solve agrees with GAMG's."""
+    installs multigrid, the solve agrees with GAMG's, and it takes the
+    iterations it takes on a partition where every rank owns DOFs. The
+    transfer's input, each level's node coordinates, is pinned exactly by the
+    tests above."""
     base = _owned_nothing_dm()
     fine = base.refine()
     coarse = uw.discretisation.Mesh(base.clone(), degree=1, qdegree=3, boundaries=_WALLS)
@@ -218,7 +278,8 @@ def test_custom_fmg_over_a_level_where_a_rank_owns_nothing():
     assert solver.snes.getKSP().getPC().getType() == "mg"
     assert solver.snes.getConvergedReason() > 0
     assert "custom_mg.build" not in solver.pc_fallbacks
-    u, r = np.asarray(solver.Unknowns.u.data), np.asarray(reference.Unknowns.u.data)
+    assert solver.snes.getKSP().getIterationNumber() == FMG_ITERATIONS
+    u, r = np.asarray(solver.Unknowns.u.array), np.asarray(reference.Unknowns.u.array)
     num = uw.mpi.comm.allreduce(float(np.sum((u - r) ** 2)))
     den = uw.mpi.comm.allreduce(float(np.sum(r ** 2)))
     assert (num / den) ** 0.5 < 1.0e-6
diff --git a/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
index 4143bbcb6..7930a39e4 100644
--- a/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
+++ b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
@@ -11,8 +11,10 @@
 together. Here the builder raises on rank 1 only. Contract: every rank records
 the RBF rescue, geometric multigrid is installed, and the solve converges. A
 check on allreduced sizes fails on every rank at once; it degrades every rank
-to the default preconditioner together, as in serial. A thread timeout turns
-the old hang into a failure.
+to the default preconditioner together, as in serial. Before the fix the
+ranks hung in the transfer's assembly. The thread timeout stops a rank waiting
+in Python; a rank held inside a PETSc collective keeps the GIL, so a hang
+there is ended by the MPI supervisor's silence limit instead.
 
     mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py
 """

From b1481152867bd1aaceb36e1453234139ada5972b Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 18:19:49 +1100
Subject: [PATCH 11/13] Text that still named PETSc's createInterpolation where
 UW3's copy now runs (#902)

field_projection's module docstring, the XDMF guide's account of vertex
projection, and the base-level coordinate repair in _DMLevelView, which
said a non-FE coordinate field gives a zero matrix. UW3's copy refuses
such a field instead.

Underworld development team with AI support from Claude Code
---
 docs/beginner/create_xdmf.md                 |  5 +++--
 src/underworld3/function/field_projection.py |  7 ++++---
 src/underworld3/utilities/custom_mg.py       | 12 ++++++------
 3 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/docs/beginner/create_xdmf.md b/docs/beginner/create_xdmf.md
index b01f10993..18112ec74 100644
--- a/docs/beginner/create_xdmf.md
+++ b/docs/beginner/create_xdmf.md
@@ -122,8 +122,9 @@ Then:
 
 ### How vertex data is obtained
 
-The compatibility group writer uses PETSc's `createInterpolation` to project
-data to mesh vertices:
+The compatibility group writer projects data to mesh vertices with the
+interpolation matrix PETSc's `createInterpolation` builds (UW3 assembles its own
+copy of it, which also works on a rank that owns none of its cells' nodes):
 
 - **P1 continuous variables**: DOFs match mesh vertices exactly, so the owned
   global-vector data is written directly to `/vertex_fields/` — no interpolation needed.
diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py
index 61a3456d1..b0551e484 100644
--- a/src/underworld3/function/field_projection.py
+++ b/src/underworld3/function/field_projection.py
@@ -1,8 +1,9 @@
 """Project or prolong MeshVariable data to a different polynomial degree.
 
-Uses a scratch PETSc DM with ``createInterpolation`` — no MeshVariable
-is created, the mesh DM is never modified, and all scratch objects are
-destroyed before returning.
+Uses a scratch PETSc DM and the interpolation matrix between it and the
+variable's sub-DM (UW3's copy of PETSc's ``createInterpolation``, #902) — no
+MeshVariable is created, the mesh DM is never modified, and all scratch
+objects are destroyed before returning.
 
 Typical use cases
 -----------------
diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py
index b2cbc5705..b8d8b4b34 100644
--- a/src/underworld3/utilities/custom_mg.py
+++ b/src/underworld3/utilities/custom_mg.py
@@ -1895,12 +1895,12 @@ def __init__(self, dm, fine_mesh):
         self.qdegree = fine_mesh.qdegree
 
         # The gmsh-imported BASE level carries section-only coordinates (its
-        # coordinate field is a PetscContainer, no PetscFE), and
-        # DMCreateInterpolation from such a source silently returns a ZERO
-        # matrix (measured) — _get_coords_for_basis would then hand every
-        # level-0 node the coordinate (0,0) and the transfer build collapses.
-        # refine() gives the child levels an FE coordinate space; give the
-        # clone's base the same: P1 Lagrange on the identical vertex layout.
+        # coordinate field is a PetscContainer, no PetscFE). PETSc's
+        # DMCreateInterpolation from such a source silently returned a ZERO
+        # matrix (measured), and UW3's copy, which _get_coords_for_basis now
+        # uses (#902), refuses a field that is not a PetscFE. refine() gives
+        # the child levels an FE coordinate space; give the clone's base the
+        # same: P1 Lagrange on the identical vertex layout.
         cdm = self.dm.getCoordinateDM()
         field = cdm.getField(0)
         fobj = field[0] if isinstance(field, tuple) else field

From 2a9fa3ea50dc169f1fd84b1a3f33fe6836cefd78 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 22:17:02 +1100
Subject: [PATCH 12/13] #886: verdicts on their own communicator; an agreed
 non-builder failure is not retried with RBF

Review round 3.

- The split replay catches any exception before its verdict. A rank that
  raised before one of the split's own allgathers then entered the verdict's
  allgather on COMM_WORLD, and paired with the others' allgather inside the
  split. It took their data as the failure reason ("(None, 76)" for an
  IndexError) and left with a different hierarchy. Verdicts now go through
  uw.mpi.verdict_comm, a Dup of COMM_WORLD used for nothing else, so they
  meet only verdicts. The raise before a collective still hangs, at a known
  place: the raising rank in the verdict, the others in the split. That part
  is #874. The two comments that said it hangs "inside the split" now say
  where.
- A CollectiveBuildError that is not a builder failure (the provenance
  checks, the operator sizes) is no longer retried with the RBF builder.
  The provenance transfer uses no builder, so the retry repeated the same
  failure after a full rebuild, and recorded a builder fallback that never
  happened. CollectiveBuildError carries builder_failure, set by the agreed
  builder verdict only. The degrade record names only the builders tried.
- test_1089 gains the provenance degrade case. One fine node on rank 1 is
  made to belong to no cell, and every rank degrades to GAMG with one
  provenance call and no builder fallback, to the serial slip. On the
  previous build it failed with two calls ([2, 2] != [1, 1]). test_1089,
  test_1107 and test_1108 pass at np 2 and 3.

Underworld development team with AI support from Claude Code
---
 src/underworld3/mpi.py                        |  8 ++++
 src/underworld3/utilities/custom_mg.py        | 20 +++++++--
 src/underworld3/utilities/fault_contact.py    |  4 +-
 src/underworld3/utilities/fault_mg.py         |  8 ++--
 ...est_1089_fault_split_multigrid_parallel.py | 44 ++++++++++++++++++-
 5 files changed, 74 insertions(+), 10 deletions(-)

diff --git a/src/underworld3/mpi.py b/src/underworld3/mpi.py
index 97761d3d2..e804375b0 100644
--- a/src/underworld3/mpi.py
+++ b/src/underworld3/mpi.py
@@ -55,6 +55,14 @@
 size = comm.size
 rank = comm.rank
 
+# The communicator of collective error verdicts, and of nothing else. A verdict
+# exchanged on ``comm`` can pair with a collective of the same kind that
+# another rank is still inside (a rank that raised before it, the others in
+# it): the two then trade unrelated data, and the ranks leave with different
+# results. On its own communicator a verdict meets only verdicts, and a rank
+# that never reaches it leaves the others waiting at a known place (#886).
+verdict_comm = comm.Dup()
+
 # State tracking for selective execution
 _in_selective_ranks = False
 _selective_executing_ranks = None
diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py
index b8d8b4b34..8c7e7bcda 100644
--- a/src/underworld3/utilities/custom_mg.py
+++ b/src/underworld3/utilities/custom_mg.py
@@ -227,7 +227,14 @@ def _drop_structural_zeros(P_csr, tol=1e-12):
 class CollectiveBuildError(RuntimeError):
     """A hierarchy-build failure raised on every rank together, from a verdict
     every rank holds (a builder that failed on any rank, a check on allreduced
-    sizes), so the ranks can retry or degrade together (#886)."""
+    sizes), so the ranks can retry or degrade together (#886).
+    ``builder_failure`` says the transfer builder itself failed, the one case
+    a different builder can change; any other check failing again under
+    another builder, the build degrades at once."""
+
+    def __init__(self, message, builder_failure=False):
+        super().__init__(message)
+        self.builder_failure = builder_failure
 
 
 def _require_serial(where):
@@ -927,7 +934,7 @@ def _agreed_node_transfer(builder, coarse_coords, fine_coords, comm):
         Pn, problem = None, f"rank {comm.getRank()}: {exc}"
     problems = [p for p in comm.tompi4py().allgather(problem) if p is not None]
     if problems:
-        raise CollectiveBuildError("; ".join(problems))
+        raise CollectiveBuildError("; ".join(problems), builder_failure=True)
     return Pn
 
 
@@ -2119,7 +2126,12 @@ def build_transfers(solver, field_id=None):
             break
         except _retryable as exc:   # sanctioned: retried or degraded, and recorded
             import warnings
-            if _i + 1 < len(_attempts):
+            # another builder can only rescue a failure of the builder: an
+            # agreed check elsewhere (the provenance transfer, the operator
+            # sizes) would fail again after a full rebuild
+            retry = (exc.builder_failure if isinstance(exc, CollectiveBuildError)
+                     else True)
+            if retry and _i + 1 < len(_attempts):
                 solver._record_pc_fallback(
                     "custom_mg.transfer_builder",
                     requested=_b,
@@ -2139,7 +2151,7 @@ def build_transfers(solver, field_id=None):
                 continue
             solver._record_pc_fallback(
                 "custom_mg.build",
-                requested=f"custom-P geometric MG ({' -> '.join(_attempts)})",
+                requested=f"custom-P geometric MG ({' -> '.join(_attempts[:_i + 1])})",
                 installed="default preconditioner",
                 reason="build_failed",
                 detail=str(exc))
diff --git a/src/underworld3/utilities/fault_contact.py b/src/underworld3/utilities/fault_contact.py
index 014893c34..92109165d 100644
--- a/src/underworld3/utilities/fault_contact.py
+++ b/src/underworld3/utilities/fault_contact.py
@@ -70,7 +70,9 @@ def _collective_raise(problem):
     where an error was meant. COLLECTIVE — every rank must reach it.
     """
     if mpi.size > 1:
-        gathered = mpi.comm.allgather(problem)
+        # on the verdict communicator: a rank still inside a collective of
+        # the path cannot be paired with this exchange
+        gathered = mpi.verdict_comm.allgather(problem)
         problems = [p for p in gathered if p is not None]
         problem = problems[0] if problems else None
     if problem is not None:
diff --git a/src/underworld3/utilities/fault_mg.py b/src/underworld3/utilities/fault_mg.py
index 6cf36f61c..65af717cc 100644
--- a/src/underworld3/utilities/fault_mg.py
+++ b/src/underworld3/utilities/fault_mg.py
@@ -620,9 +620,11 @@ def split_tail_levels(levels, fault_specs, verbose=False):
         # the split's refusals are collective; a failure on one rank after its
         # last collective is not, and every rank must stop at the same level or
         # the ranks hold different hierarchies. Any exception is caught, so an
-        # IndexError on one rank meets the same verdict (#886); a rank-local
-        # raise inside the split, before its own collectives, still hangs
-        # there (issue #874)
+        # IndexError on one rank meets the same verdict (#886). A rank-local
+        # raise inside the split, before its own collectives, still hangs: the
+        # raising rank waits in the verdict, which has its own communicator so
+        # it cannot pair with the split's collectives, and the others wait in
+        # the split (issue #874)
         try:
             _collective_raise(problem)
         except RuntimeError as exc:
diff --git a/tests/parallel/test_1089_fault_split_multigrid_parallel.py b/tests/parallel/test_1089_fault_split_multigrid_parallel.py
index 4ed9792be..187648001 100644
--- a/tests/parallel/test_1089_fault_split_multigrid_parallel.py
+++ b/tests/parallel/test_1089_fault_split_multigrid_parallel.py
@@ -45,6 +45,46 @@ def test_split_levels_converge_like_serial():
     assert abs(slip - SERIAL_PEAK_SLIP) < 1.0e-4 * SERIAL_PEAK_SLIP, slip
 
 
+def test_an_agreed_provenance_failure_degrades_every_rank(monkeypatch):
+    """The provenance transfer checks allreduced counts ("a fine node belongs
+    to no cell"), so every rank raises CollectiveBuildError together and the
+    multigrid build degrades on every rank, not on one (#886). A different
+    transfer builder cannot help, since the provenance transfer uses none, so
+    there is no retry with the RBF builder (which repeated the same failure
+    after a full rebuild). Here one fine node on rank 1 is made to belong to no
+    cell; the solve then runs on the default preconditioner, to the same slip."""
+    from underworld3.utilities import fault_mg
+    original_owner = fault_mg._fine_cell_of_node
+    original_provenance = fault_mg.provenance_transfer_parallel
+    inside, calls = [False], [0]
+
+    def counted(*args, **kwargs):
+        calls[0] += 1
+        inside[0] = True
+        try:
+            return original_provenance(*args, **kwargs)
+        finally:
+            inside[0] = False
+
+    def one_orphan_on_rank_one(fn, n_f):
+        owner = original_owner(fn, n_f)
+        if inside[0] and uw.mpi.rank == 1 and owner.size:
+            owner = owner.copy()
+            owner[0] = -1
+        return owner
+
+    monkeypatch.setattr(fault_mg, "_fine_cell_of_node", one_orphan_on_rank_one)
+    monkeypatch.setattr(fault_mg, "provenance_transfer_parallel", counted)
+    mesh = _faulted_box(2)
+    with pytest.warns(UserWarning, match="custom-P FMG build failed"):
+        stokes, info, slip = shear_fault_solve(mesh)
+    assert info["converged"] and info["velocity_pc"] == "GAMG", info
+    assert uw.mpi.comm.allgather(calls[0]) == [1] * uw.mpi.size
+    assert "custom_mg.build" in stokes.pc_fallbacks
+    assert "custom_mg.transfer_builder" not in stokes.pc_fallbacks
+    assert abs(slip - SERIAL_PEAK_SLIP) < 1.0e-4 * SERIAL_PEAK_SLIP, slip
+
+
 def test_the_welded_tail_runs_in_parallel_too():
     info, slip = _solve(_faulted_box(2, split_levels=False))
     assert info["velocity_pc"] == "custom-FMG" and info["converged"], info
@@ -103,8 +143,8 @@ def test_a_failure_on_one_rank_stops_the_replay_on_every_rank(monkeypatch, error
     every rank keeps the same split levels and the levels from the failure down
     stay welded. Any exception meets the verdict, not only the split's own
     refusal types: an IndexError escaped it on one rank (#886). (A rank-local
-    raise BEFORE the split's own collectives still hangs inside the split;
-    issue #874.)"""
+    raise BEFORE the split's own collectives still hangs, the raising rank in
+    the verdict and the others in the split; issue #874.)"""
     import underworld3.utilities.fault_split as fault_split
     original, calls = fault_split.add_fault, [0]
 

From 57d87783879cb2ccbc34dd85c3bc7c9f7a4fe1e2 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 31 Oct 2026 22:17:05 +1100
Subject: [PATCH 13/13] #902: the coordinate space keeps PETSc's route only for
 isoperiodic meshes; no space, no change

Review round 3.

- In PETSc v3.25 only isoperiodic coordinates (DMPlexGetIsoperiodicFaceSF)
  are projected differently: locally, with no matrix. maxCell/L periodicity
  projects its non-localized coordinates through DMCreateInterpolation and
  MatMult, the step UW's copy replaces. Such meshes were sent to PETSc's
  route, which still fails on a rank owning no nodes. Now only the
  isoperiodic case keeps it.
- For a cell type PETSc builds no Lagrange space for (tensor prisms), the
  coordinate DM is not replaced, and the copy called PetscFEGetDualSpace on
  the old PetscContainer field. Now an unchanged coordinate DM leaves the
  coordinates alone, as PETSc does. No UW constructor builds such cells.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h           | 15 ++++++++-------
 src/underworld3/cython/petsc_discretisation.pyx |  3 ++-
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 65106f44a..fbaf93892 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -557,11 +557,12 @@ PetscErrorCode UW_DMCreateInterpolationWithinPlex(DM dmc, DM dmf, Mat *In)
 // this plex, so a degree >= 2 space could not be built on a rank that holds
 // cells but owns none of their nodes (#902). As in PETSc: the same space is a
 // copy; coordinates injected into a higher coordinate dimension are not
-// projected; periodic coordinates keep PETSc's own projection, which projects
-// isoperiodic ones locally, without the matrix.
+// projected; a cell type PETSc builds no Lagrange space for (tensor prisms)
+// keeps its coordinates. Isoperiodic coordinates keep PETSc's own projection,
+// which is local there, without the matrix; maxCell/L periodicity projects
+// the non-localized coordinates through the matrix, as here.
 PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(DM dm, PetscInt degree)
 {
-    const PetscReal *maxCell, *Lstart, *L;
     PetscInt         nfaceSF = 0;
     DM               cdmOld, cdmNew;
     Vec              xOld, xNew, x;
@@ -570,9 +571,8 @@ PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(DM dm, PetscInt degree)
     PetscBool        same;
 
     PetscFunctionBeginUser;
-    PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
     PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &nfaceSF, NULL));
-    if (maxCell || L || nfaceSF) {
+    if (nfaceSF) {
         PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_TRUE));
         PetscFunctionReturn(PETSC_SUCCESS);
     }
@@ -583,8 +583,9 @@ PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(DM dm, PetscInt degree)
     // installs the new coordinate DM and leaves the coordinates to us
     PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_FALSE));
     PetscCall(DMGetCoordinates(dm, &x));
-    if (x == xOld) {
-        PetscCall(DMGetCoordinateDM(dm, &cdmNew));
+    PetscCall(DMGetCoordinateDM(dm, &cdmNew));
+    // no new coordinate DM: PETSc installed no space for this cell type
+    if (x == xOld && cdmNew != cdmOld) {
         PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
         PetscCall(DMGetField(cdmNew, 0, NULL, (PetscObject *)&disc));
         PetscCall(PetscFEGetDualSpace(discOld, &dsOld));
diff --git a/src/underworld3/cython/petsc_discretisation.pyx b/src/underworld3/cython/petsc_discretisation.pyx
index fb32a5c40..364be1ecb 100644
--- a/src/underworld3/cython/petsc_discretisation.pyx
+++ b/src/underworld3/cython/petsc_discretisation.pyx
@@ -208,7 +208,8 @@ def petsc_dm_create_coordinate_space_within_plex(dm, degree):
         degree >= 2 space is also built on a rank that holds cells but owns
         none of their nodes (#902; see
         ``UW_DMPlexCreateCoordinateSpaceWithinPlex`` in ``petsc_compat.h``).
-        Collective.
+        Isoperiodic coordinates are the exception: PETSc projects them
+        locally, with no matrix, and they keep its route. Collective.
         """
         cdef DM c_dm = dm
         cdef PetscInt c_degree = degree