From 31eb3e1085148b6a4c14ef30755b5782348fce2d Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 16 Aug 2026 17:51:05 +1000
Subject: [PATCH 1/5] Clip the fault assembly against the mesh's own boundary,
 built as OCC geometry
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The assembly clip (place_thin_volume, 2-D and 3-D) now intersects with the
domain the mesh actually has: unshared support-1 facets are gathered
collectively, stitched by exact coordinate identity, and rebuilt in OCC —
collinear runs (coplanar triangles, one dimension up) compressed to their
corners, so a box constructs the same four-sided / six-faced tool the
analytic addRectangle/addBox clip did, and the boolean does not imprint the
wall's mesh spacing onto the clipped faces. The post-mesh snap becomes a
pure normal-component move onto the boundary's own facets, which on an
axis-aligned wall reproduces the wall-plane snap exactly.

The cut lands on the mesh's own facets (measured 3.5e-17 from the chords,
against the 8.6e-3 sagitta a cut on the smooth circle would leave), so no
corefinement library is needed and clip-then-mesh order is preserved.

Box regression oracle unchanged: test_0854, test_0855 serial and at
np=2/np=4. New tool unit tests in test_0859 assert exactness both ways —
OCC mass equals the mesh's summed cell measure to round-off, and differs
from the smooth-geometry measure by far more (the negative control).

Outcrop-band identification still runs against the axis-aligned box frame;
generalising the band to the boundary trace is the next stage.

Underworld development team with AI support from Claude Code
---
 src/underworld3/utilities/place_surface.py | 397 ++++++++++++++++++---
 tests/test_0859_domain_boundary_tool.py    | 151 ++++++++
 2 files changed, 502 insertions(+), 46 deletions(-)
 create mode 100644 tests/test_0859_domain_boundary_tool.py

diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py
index 6d3a1e39..ff8efc18 100644
--- a/src/underworld3/utilities/place_surface.py
+++ b/src/underworld3/utilities/place_surface.py
@@ -2677,6 +2677,305 @@ def place_sheet(dm, points, triangles, label=CUT_LABEL, label_value=1,
 # ``~/+Simulations/mesh_reconnection_study/thin_volume_spike.py`` — widths
 # h, h/2, h/4 and junction angles down to 10 degrees, all gated.
 
+# ------------------------------------------- the domain boundary as a tool
+
+def _domain_boundary_facets(dm):
+    """The DOMAIN boundary as one small global complex, identical everywhere.
+
+    A facet of the domain boundary has support 1 AND is unshared — a
+    partition-seam facet also has local support 1, and clipping against a
+    seam would eat the mesh differently at every rank count (the
+    :func:`_true_wall_vertex_mask` distinction). Each boundary facet lives
+    on exactly one rank, so the gathered set holds each facet once; the
+    facets are stitched by exact coordinate identity, which is sound
+    because a shared vertex's coordinates are copies of the same bytes on
+    every rank. COLLECTIVE.
+
+    Returns ``(verts, facets)``: coordinates ``(nv, dim)`` and vertex-index
+    facets ``(nf, dim)`` — edges in 2-D, triangles in 3-D.
+    """
+    dim = dm.getDimension()
+    vS, vE = dm.getDepthStratum(0)
+    pStart, _pEnd = dm.getChart()
+    X = _coords(dm)[: vE - vS]
+    shared = _shared_point_flags(dm).astype(bool)
+    corners = []
+    for f in range(*dm.getHeightStratum(1)):
+        if len(dm.getSupport(f)) == 1 and not shared[f - pStart]:
+            vv = [int(q) - vS for q in dm.getTransitiveClosure(f)[0]
+                  if vS <= int(q) < vE]
+            corners.append(X[vv])
+    local = (np.asarray(corners, dtype=float) if corners
+             else np.zeros((0, dim, dim), dtype=float))
+    pieces = [g for g in uw.mpi.comm.allgather(local) if len(g)]
+    if not pieces:
+        raise RuntimeError("the mesh has no domain boundary facet")
+    stack = np.concatenate(pieces)
+    verts, inverse = np.unique(stack.reshape(-1, dim), axis=0,
+                               return_inverse=True)
+    return verts, inverse.reshape(-1, dim).astype(np.int64)
+
+
+def _domain_loops_2d(dm):
+    """The 2-D domain boundary as closed coordinate loops. COLLECTIVE."""
+    verts, edges = _domain_boundary_facets(dm)
+    loops = _skin_loops([(int(a), int(b)) for a, b in edges],
+                        what="the domain boundary")
+    return [verts[loop] for loop in loops]
+
+
+def _compress_collinear_loop(loop_xy):
+    """The loop's vertices where the boundary actually TURNS.
+
+    The OCC tool needs the boundary's geometry, not its segmentation: a box
+    wall's collinear run collapses to its corners, so a box builds the same
+    four-sided (six-faced, one dimension up) tool the analytic box gave,
+    and the boolean does not imprint the wall's mesh spacing onto the
+    clipped faces. The compressed polygon carries the identical point set,
+    so the cut itself is unchanged.
+    """
+    P = np.asarray(loop_xy, dtype=float)
+    e1 = P - np.roll(P, 1, axis=0)
+    e2 = np.roll(P, -1, axis=0) - P
+    if P.shape[1] == 2:
+        cross = np.abs(e1[:, 0] * e2[:, 1] - e1[:, 1] * e2[:, 0])
+    else:
+        cross = np.linalg.norm(np.cross(e1, e2), axis=1)
+    turn = cross > (1e-12 * np.linalg.norm(e1, axis=1)
+                    * np.linalg.norm(e2, axis=1))
+    if turn.sum() < 3:
+        raise RuntimeError(
+            "a domain boundary loop compresses to fewer than three corners")
+    return P[turn]
+
+
+def _occ_domain_2d(occ, loops):
+    """The domain as ONE OCC plane surface built from its boundary loops.
+
+    ``loops`` are the (compressed) boundary polygons of the mesh's own
+    vertices; the largest-area loop is the exterior, the rest are holes —
+    an annulus arrives natively. The tool is exactly the discrete domain,
+    so a boolean against it lands on the mesh's own facets (measured:
+    3.5e-17 from the chords, against the 8.6e-3 sagitta a cut on the smooth
+    circle would have left). Returns the surface tag; the caller
+    synchronizes.
+    """
+    def area(P):
+        return 0.5 * float(P[:, 0] @ np.roll(P[:, 1], -1)
+                           - P[:, 1] @ np.roll(P[:, 0], -1))
+
+    order = sorted(range(len(loops)), key=lambda k: -abs(area(loops[k])))
+    rings = []
+    for k in order:
+        P = loops[k]
+        pts = [occ.addPoint(x, y, 0.0) for x, y in P]
+        lines = [occ.addLine(pts[i], pts[(i + 1) % len(pts)])
+                 for i in range(len(pts))]
+        rings.append(occ.addCurveLoop(lines))
+    return occ.addPlaneSurface(rings)
+
+
+def _snap_to_boundary_2d(xy, loops, tol=1e-9):
+    """Snap nodes within ``tol`` of the domain boundary ONTO it, exactly.
+
+    OCC's clipped edges sit within rounding of the tool; the band logic and
+    the sew need EXACT membership. A node near a boundary corner takes the
+    corner; a node near a segment loses only its normal component — on an
+    axis-aligned wall that is exactly the wall-plane snap the box clip
+    used, the tangential coordinate untouched.
+    """
+    for P in loops:
+        for v in P:
+            xy[np.linalg.norm(xy - v, axis=1) < tol] = v
+        n = len(P)
+        for i in range(n):
+            A, B = P[i], P[(i + 1) % n]
+            e = B - A
+            length = float(np.linalg.norm(e))
+            nrm = np.array([-e[1], e[0]]) / length
+            off = (xy - A) @ nrm
+            u = ((xy - A) @ e) / (length * length)
+            near = (np.abs(off) < tol) & (u > 0.0) & (u < 1.0)
+            xy[near] -= off[near, None] * nrm
+    return xy
+
+
+def _occ_domain_3d(occ, verts, tris):
+    """The domain as ONE OCC solid built from its boundary triangles.
+
+    Adjacent coplanar triangles merge into single planar faces first — a
+    box wall becomes one rectangle, so a box builds the same six-faced tool
+    the analytic box gave, and the boolean does not imprint the wall's mesh
+    spacing onto the clipped faces — and each merged face's rim compresses
+    to its corners. The chain shared by two merged faces lies in both
+    planes, hence on a straight line, so both faces compress it to the same
+    endpoints; a vertex where three or more faces meet is always kept, or
+    two faces whose planes cross the same line would disagree about it.
+    Faces share OCC points and lines, so the shells they close into are
+    topologically sewn. Returns ``(volume_tag, planes)`` with ``planes`` a
+    list of ``(anchor, unit_normal)`` per merged face, for the snap; the
+    caller synchronizes.
+    """
+    tris = np.asarray(tris, dtype=np.int64)
+    T = verts[tris]
+    raw_n = np.cross(T[:, 1] - T[:, 0], T[:, 2] - T[:, 0])
+    nn = np.linalg.norm(raw_n, axis=1)
+    if (nn == 0.0).any():
+        raise RuntimeError("a domain boundary triangle is degenerate")
+    unit_n = raw_n / nn[:, None]
+
+    parent = list(range(len(tris)))
+
+    def find(a):
+        while parent[a] != a:
+            parent[a] = parent[parent[a]]
+            a = parent[a]
+        return a
+
+    edge_first = {}
+    for t, tri in enumerate(tris):
+        opp = {frozenset((int(tri[0]), int(tri[1]))): int(tri[2]),
+               frozenset((int(tri[1]), int(tri[2]))): int(tri[0]),
+               frozenset((int(tri[2]), int(tri[0]))): int(tri[1])}
+        for key, far in opp.items():
+            if key not in edge_first:
+                edge_first[key] = t
+                continue
+            s = edge_first[key]
+            scale = max(nn[t], nn[s]) ** 0.5
+            coplanar = (np.linalg.norm(np.cross(unit_n[t], unit_n[s]))
+                        < 1e-12
+                        and abs((verts[far] - T[s, 0]) @ unit_n[s])
+                        < 1e-12 * scale)
+            if coplanar:
+                parent[find(t)] = find(s)
+
+    regions = {}
+    for t in range(len(tris)):
+        regions.setdefault(find(t), []).append(t)
+    regions_at = {}
+    for t, tri in enumerate(tris):
+        r = find(t)
+        for v in tri:
+            regions_at.setdefault(int(v), set()).add(r)
+
+    point_tag, line_tag = {}, {}
+
+    def point(v):
+        if v not in point_tag:
+            point_tag[v] = occ.addPoint(*verts[v])
+        return point_tag[v]
+
+    def line(a, b):
+        key = (a, b) if a < b else (b, a)
+        if key not in line_tag:
+            line_tag[key] = occ.addLine(point(key[0]), point(key[1]))
+        return line_tag[key]
+
+    def compress_rim(loop):
+        P = verts[np.asarray(loop)]
+        e1 = P - np.roll(P, 1, axis=0)
+        e2 = np.roll(P, -1, axis=0) - P
+        cross = np.linalg.norm(np.cross(e1, e2), axis=1)
+        turn = cross > (1e-12 * np.linalg.norm(e1, axis=1)
+                        * np.linalg.norm(e2, axis=1))
+        junction = np.array([len(regions_at[int(v)]) >= 3 for v in loop])
+        keep = turn | junction
+        if keep.sum() < 3:
+            raise RuntimeError(
+                "a domain face's rim compresses to fewer than three "
+                "corners")
+        return [int(v) for v, k in zip(loop, keep) if k]
+
+    from collections import Counter
+
+    face_of_region = {}
+    planes = []
+    for r, members in regions.items():
+        rim = Counter()
+        for t in members:
+            a, b, c = (int(v) for v in tris[t])
+            for e in (frozenset((a, b)), frozenset((b, c)),
+                      frozenset((c, a))):
+                rim[e] += 1
+        rim_edges = [tuple(e) for e, k in rim.items() if k == 1]
+        loops = _skin_loops(rim_edges, what="a domain face's rim")
+        m = unit_n[members[0]]
+        anchor = verts[int(tris[members[0]][0])]
+        planes.append((anchor, m))
+        # In-plane coordinates for the outer-vs-hole ranking only.
+        u = T[members[0], 1] - T[members[0], 0]
+        u = u / np.linalg.norm(u)
+        w = np.cross(m, u)
+
+        def rim_area(loop):
+            Q = verts[np.asarray(loop)] - anchor
+            x, y = Q @ u, Q @ w
+            return 0.5 * float(x @ np.roll(y, -1) - y @ np.roll(x, -1))
+
+        loops = sorted((compress_rim(lp) for lp in loops),
+                       key=lambda lp: -abs(rim_area(lp)))
+        rings = []
+        for lp in loops:
+            rings.append(occ.addCurveLoop(
+                [line(lp[i], lp[(i + 1) % len(lp)])
+                 for i in range(len(lp))]))
+        face_of_region[r] = occ.addPlaneSurface(rings)
+
+    # Shells: connected components of the triangulation; the component
+    # enclosing the greatest volume is the outer boundary, the others are
+    # interior cavities (a spherical shell's inner surface).
+    sparent = list(range(len(tris)))
+
+    def sfind(a):
+        while sparent[a] != a:
+            sparent[a] = sparent[sparent[a]]
+            a = sparent[a]
+        return a
+
+    edge_seen = {}
+    for t, tri in enumerate(tris):
+        a, b, c = (int(v) for v in tri)
+        for e in (frozenset((a, b)), frozenset((b, c)), frozenset((c, a))):
+            if e in edge_seen:
+                sparent[sfind(t)] = sfind(edge_seen[e])
+            else:
+                edge_seen[e] = t
+    shells, shell_tris = {}, {}
+    for t in range(len(tris)):
+        s = sfind(t)
+        shells.setdefault(s, set()).add(find(t))
+        shell_tris.setdefault(s, []).append(t)
+    enclosed = {}
+    for s, ts in shell_tris.items():
+        Q = T[ts]
+        enclosed[s] = abs(float(np.einsum(
+            "ij,ij->i", np.cross(Q[:, 0], Q[:, 1]), Q[:, 2]).sum()) / 6.0)
+    order = sorted(shells, key=lambda s: -enclosed[s])
+    loops3 = [occ.addSurfaceLoop([face_of_region[r] for r in shells[s]])
+              for s in order]
+    return occ.addVolume(loops3), planes
+
+
+def _snap_to_boundary_3d(xyz, verts, tris, planes, tol=1e-9):
+    """Snap nodes within ``tol`` of the domain boundary ONTO it, exactly.
+
+    Boundary corners first, then each merged face's plane, each a pure
+    normal-component move — on an axis-aligned wall that is exactly the
+    wall-plane snap the box clip used, the in-plane coordinates untouched.
+    A node near two planes (a domain edge) is snapped by both passes, as
+    the box path snapped both axes at a corner.
+    """
+    for v in np.unique(tris):
+        p = verts[int(v)]
+        xyz[np.linalg.norm(xyz - p, axis=1) < tol] = p
+    for anchor, m in planes:
+        off = (xyz - anchor) @ m
+        near = np.abs(off) < tol
+        xyz[near] -= off[near, None] * m
+    return xyz
+
+
 def _patch_frame(patch):
     """Unit normal of a planar patch, with planarity asserted."""
     P = np.asarray(patch, dtype=float)
@@ -2693,16 +2992,18 @@ def _patch_frame(patch):
     return n
 
 
-def _occ_assembly_3d(patches, width, size, box=None, assembly="fuse"):
+def _occ_assembly_3d(patches, width, size, domain=None, assembly="fuse"):
     """Thicken each planar patch by ±width/2, resolve overlaps, mesh.
 
     ``assembly`` is :func:`place_thin_volume`'s: ``"fuse"`` returns the union
     as one solid, ``"fragment"`` keeps every overlap piece.
 
-    ``box = (lo, hi)`` applies the specify-long contract: the thickened
-    solids are INTERSECTED with the domain box, so patches may protrude —
-    an assembly reaching the top surface leaves its clipped face exactly in
-    the wall plane (snapped to the plane value after meshing, defensively).
+    ``domain = (verts, tris)`` — the mesh's boundary complex
+    (:func:`_domain_boundary_facets`) — applies the specify-long contract:
+    the thickened solids are INTERSECTED with the domain built as OCC
+    geometry from those very facets, so patches may protrude — an assembly
+    reaching a boundary leaves its clipped face exactly on the boundary's
+    own facets (snapped onto their planes after meshing, defensively).
 
     Returns ``(points, tets, cad_volume)`` — the assembly mesh in its own
     numbering, and the CAD volume of the (clipped) pieces, against which the
@@ -2736,11 +3037,12 @@ def _occ_assembly_3d(patches, width, size, box=None, assembly="fuse"):
                 occ.fuse([(3, solids[0])], [(3, t) for t in solids[1:]])
             else:
                 occ.fragment([(3, solids[0])], [(3, t) for t in solids[1:]])
-        if box is not None:
-            lo, hi = (np.asarray(b, dtype=float) for b in box)
+        planes = None
+        if domain is not None:
+            dom_verts, dom_tris = domain
             occ.synchronize()
             solids = [t for _d, t in gmsh.model.getEntities(3)]
-            tool = occ.addBox(*lo, *(hi - lo))
+            tool, planes = _occ_domain_3d(occ, dom_verts, dom_tris)
             occ.intersect([(3, t) for t in solids], [(3, tool)])
         occ.synchronize()
 
@@ -2764,15 +3066,11 @@ def _occ_assembly_3d(patches, width, size, box=None, assembly="fuse"):
                                          dtype=np.int64).reshape(-1, 4))
         if not tets:
             raise RuntimeError("the assembly meshed to no tetrahedra")
-        if box is not None:
-            # OCC's clipped faces sit within rounding of the plane; the
-            # band logic and the cap's node sharing need EXACT plane
-            # values, so snap defensively.
-            lo, hi = (np.asarray(b, dtype=float) for b in box)
-            for axis in range(3):
-                for value in (lo[axis], hi[axis]):
-                    near = np.abs(xyz[:, axis] - value) < 1e-9
-                    xyz[near, axis] = value
+        if planes is not None:
+            # OCC's clipped faces sit within rounding of the boundary; the
+            # band logic and the cap's node sharing need EXACT membership,
+            # so snap defensively.
+            xyz = _snap_to_boundary_3d(xyz, dom_verts, dom_tris, planes)
         return xyz, np.vstack(tets), float(cad_volume)
     finally:
         gmsh.finalize()
@@ -3041,7 +3339,7 @@ def cap_tag_of(k):
         gmsh.finalize()
 
 
-def _occ_assembly_2d(polylines, width, size, assembly="fuse", box=None):
+def _occ_assembly_2d(polylines, width, size, assembly="fuse", domain=None):
     """Thicken each polyline into a ribbon, resolve overlaps, mesh.
 
     The 2-D thin volume: a ribbon is the mitred outline of one polyline, and
@@ -3052,12 +3350,14 @@ def _occ_assembly_2d(polylines, width, size, assembly="fuse", box=None):
     honour. Returns ``(points, triangles, cad_area)``, the area being that of
     the resolved faces — the union — under either choice.
 
-    ``box = (lo, hi)`` applies the specify-long contract one dimension down
-    from :func:`_occ_assembly_3d`: the resolved faces are INTERSECTED with
-    the domain rectangle, so polylines may protrude — a ribbon reaching the
-    top surface leaves its clipped edge exactly in the wall line (snapped to
-    the line value after meshing, defensively). ``cad_area`` is then the
-    clipped area, so the meshed-vs-CAD gate holds unchanged.
+    ``domain`` — the mesh's boundary loops (:func:`_domain_loops_2d`) —
+    applies the specify-long contract one dimension down from
+    :func:`_occ_assembly_3d`: the resolved faces are INTERSECTED with the
+    domain built as OCC geometry from those very loops, so polylines may
+    protrude — a ribbon reaching a boundary leaves its clipped edge exactly
+    on the boundary's own facets (snapped onto them after meshing,
+    defensively). ``cad_area`` is then the clipped area, so the
+    meshed-vs-CAD gate holds unchanged.
     """
     import gmsh
 
@@ -3125,12 +3425,13 @@ def outline(P):
                 occ.fuse([(2, surfs[0])], [(2, t) for t in surfs[1:]])
             else:
                 occ.fragment([(2, surfs[0])], [(2, t) for t in surfs[1:]])
-        if box is not None:
-            lo, hi = (np.asarray(b, dtype=float) for b in box)
+        comp = None
+        if domain is not None:
+            comp = [_compress_collinear_loop(np.asarray(L, dtype=float))
+                    for L in domain]
             occ.synchronize()
             faces = [t for _d, t in gmsh.model.getEntities(2)]
-            tool = occ.addRectangle(lo[0], lo[1], 0.0,
-                                    hi[0] - lo[0], hi[1] - lo[1])
+            tool = _occ_domain_2d(occ, comp)
             occ.intersect([(2, t) for t in faces], [(2, tool)])
         occ.synchronize()
 
@@ -3153,15 +3454,11 @@ def outline(P):
                                          dtype=np.int64).reshape(-1, 3))
         if not tris:
             raise RuntimeError("the ribbon assembly meshed to no triangles")
-        if box is not None:
-            # OCC's clipped edges sit within rounding of the wall line; the
-            # band logic and the wall relabel need EXACT line values, so
+        if comp is not None:
+            # OCC's clipped edges sit within rounding of the boundary; the
+            # band logic and the wall relabel need EXACT membership, so
             # snap defensively (the 3-D path does the same on its planes).
-            lo, hi = (np.asarray(b, dtype=float) for b in box)
-            for axis in range(2):
-                for value in (lo[axis], hi[axis]):
-                    near = np.abs(xy[:, axis] - value) < 1e-9
-                    xy[near, axis] = value
+            xy = _snap_to_boundary_2d(xy, comp)
         return xy, np.vstack(tris), float(cad_area)
     finally:
         gmsh.finalize()
@@ -3179,19 +3476,19 @@ def _segments_distance(X, pts, edges):
     return best
 
 
-def _skin_loops(skin_edges):
+def _skin_loops(skin_edges, what="the assembly's skin"):
     """Order a 2-D skin's edges into closed loops of vertex ids.
 
     A manifold skin gives every vertex exactly two incident edges; anything
-    else is a defect of the assembly mesh and is refused.
+    else is a defect of the mesh the edges bound and is refused.
     """
     adj = {}
     for a, b in skin_edges:
         adj.setdefault(int(a), []).append(int(b))
         adj.setdefault(int(b), []).append(int(a))
     if any(len(v) != 2 for v in adj.values()):
-        raise RuntimeError("the assembly's skin is not a set of closed "
-                           "loops; the layer mesh is defective")
+        raise RuntimeError(f"{what} is not a set of closed loops; "
+                           "the mesh it bounds is defective")
     loops, seen = [], set()
     for start in sorted(adj):
         if start in seen:
@@ -3966,8 +4263,11 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
     """
     comm = uw.mpi.comm
 
-    # The (axis-aligned) domain box, collectively — the clip target and the
-    # wall lines the band is identified against.
+    # The domain's own boundary, collectively — the clip target. The
+    # axis-aligned box survives alongside it as the frame the outcrop BAND
+    # is identified against (wall lines); a general-boundary band is the
+    # trace-labelling stage, not this one.
+    domain_loops = _domain_loops_2d(dm)
     Xb = _coords(dm)
     lo_hi = np.array([Xb.min(axis=0) if len(Xb) else np.full(2, np.inf),
                       -(Xb.max(axis=0)) if len(Xb) else np.full(2, np.inf)])
@@ -3979,7 +4279,7 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
     if comm.rank == 0:
         try:
             asm_pts, asm_tris, cad_area = _occ_assembly_2d(
-                polylines, width, size, assembly, box=(box_lo, box_hi))
+                polylines, width, size, assembly, domain=domain_loops)
             P = asm_pts[asm_tris]
             twice = ((P[:, 1, 0] - P[:, 0, 0]) * (P[:, 2, 1] - P[:, 0, 1])
                      - (P[:, 1, 1] - P[:, 0, 1]) * (P[:, 2, 0] - P[:, 0, 0]))
@@ -4430,8 +4730,12 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1,
     comm = uw.mpi.comm
 
     # The specify-long contract: patches may extend PAST the domain; the
-    # assembly is clipped against the (axis-aligned) box, and a clipped
-    # face left in a wall plane becomes the zone's OUTCROP BAND.
+    # assembly is clipped against the mesh's own boundary, and a clipped
+    # face left in a wall plane becomes the zone's OUTCROP BAND. The
+    # axis-aligned box survives alongside as the frame the band is
+    # identified against; a general-boundary band is the trace-labelling
+    # stage, not this one.
+    dom_verts, dom_tris = _domain_boundary_facets(dm)
     Xb = _coords(dm)
     lo_hi = np.array([Xb.min(axis=0) if len(Xb) else np.full(3, np.inf),
                       -(Xb.max(axis=0)) if len(Xb) else np.full(3, np.inf)])
@@ -4444,7 +4748,8 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1,
     if comm.rank == 0:
         try:
             asm_pts, asm_tets, cad_vol = _occ_assembly_3d(
-                patches, width, size, box=(box_lo, box_hi), assembly=assembly)
+                patches, width, size, domain=(dom_verts, dom_tris),
+                assembly=assembly)
             v6 = np.einsum(
                 "ij,ij->i",
                 np.cross(asm_pts[asm_tets][:, 1] - asm_pts[asm_tets][:, 0],
diff --git a/tests/test_0859_domain_boundary_tool.py b/tests/test_0859_domain_boundary_tool.py
new file mode 100644
index 00000000..9a0c0637
--- /dev/null
+++ b/tests/test_0859_domain_boundary_tool.py
@@ -0,0 +1,151 @@
+"""The domain-boundary tool (:mod:`place_surface`): the mesh's own boundary
+rebuilt as OCC geometry, so an assembly can be clipped against the domain a
+mesh ACTUALLY has rather than against an analytic box.
+
+The mesh boundary IS the domain — there is no circle an annulus is failing
+to be — so the tool must reproduce the discrete boundary exactly: the loop
+(shell, in 3-D) of the mesh's own boundary facets, collinear (coplanar) runs
+compressed to their corners so a box rebuilds the four-sided (six-faced)
+tool the analytic clip used. Exactness is asserted two ways in each
+dimension: the tool's OCC mass equals the mesh's own summed cell measure to
+round-off, and — the negative control — it DIFFERS from the smooth-geometry
+measure by far more than that tolerance, which proves the tool is the
+polygon and not the circle it approximates.
+"""
+import gmsh
+import numpy as np
+import pytest
+
+import underworld3 as uw
+from underworld3.utilities.place_surface import (
+    _compress_collinear_loop, _domain_boundary_facets, _domain_loops_2d,
+    _occ_domain_2d, _occ_domain_3d, _snap_to_boundary_2d,
+    _snap_to_boundary_3d, cell_areas, _owned_cell_volume)
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_b,
+              pytest.mark.skipif(uw.mpi.size > 1,
+                                 reason="serial suite; the parallel form "
+                                        "joins ptest_0855")]
+
+
+def _shoelace(P):
+    return 0.5 * float(P[:, 0] @ np.roll(P[:, 1], -1)
+                       - P[:, 1] @ np.roll(P[:, 0], -1))
+
+
+def _occ_mass_2d(loops):
+    gmsh.initialize()
+    gmsh.option.setNumber("General.Terminal", 0)
+    try:
+        gmsh.model.add("uw_domain_tool_test")
+        occ = gmsh.model.occ
+        tool = _occ_domain_2d(occ, [_compress_collinear_loop(L)
+                                    for L in loops])
+        occ.synchronize()
+        return float(occ.getMass(2, tool))
+    finally:
+        gmsh.finalize()
+
+
+def _occ_mass_3d(verts, tris):
+    gmsh.initialize()
+    gmsh.option.setNumber("General.Terminal", 0)
+    try:
+        gmsh.model.add("uw_domain_tool_test")
+        occ = gmsh.model.occ
+        vol, planes = _occ_domain_3d(occ, verts, tris)
+        occ.synchronize()
+        return float(occ.getMass(3, vol)), planes
+    finally:
+        gmsh.finalize()
+
+
+# ------------------------------------------------------------------------ 2-D
+
+def test_a_box_boundary_compresses_to_its_four_corners():
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.11,
+        regular=False, qdegree=2)
+    loops = _domain_loops_2d(mesh.dm)
+    assert len(loops) == 1
+    comp = _compress_collinear_loop(loops[0])
+    assert len(comp) == 4
+    assert (sorted(map(tuple, comp))
+            == [(0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0)])
+    assert _occ_mass_2d(loops) == pytest.approx(1.0, rel=1e-12)
+
+
+def test_the_annulus_tool_is_the_polygon_not_the_circle():
+    mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5,
+                              cellSize=0.2, qdegree=2)
+    loops = _domain_loops_2d(mesh.dm)
+    assert len(loops) == 2
+
+    mesh_area = float(cell_areas(mesh.dm).sum())
+    loop_area = sum(abs(_shoelace(np.asarray(L))) for L in loops)
+    outer, inner = sorted((abs(_shoelace(np.asarray(L))) for L in loops),
+                          reverse=True)
+    assert outer - inner == pytest.approx(mesh_area, rel=1e-12)
+
+    mass = _occ_mass_2d(loops)
+    assert mass == pytest.approx(mesh_area, rel=1e-9)
+    # Negative control: the polygon is NOT the smooth annulus. If the tool
+    # ever silently becomes the circle, this fires before anything subtle.
+    smooth = np.pi * (1.0 ** 2 - 0.5 ** 2)
+    assert abs(mass - smooth) > 1e3 * abs(mass - mesh_area)
+    assert loop_area > 0.0
+
+
+def test_the_2d_snap_restores_wall_values_exactly():
+    square = [np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]
+    xy = np.array([[0.3, 1e-10],          # near the bottom wall
+                   [1.0 - 1e-10, 0.7],    # near the right wall
+                   [1e-10, 1e-10],        # near a corner
+                   [0.5, 2e-9],           # OUTSIDE the snap tolerance
+                   [0.4, 0.6]])           # interior
+    out = _snap_to_boundary_2d(xy.copy(), square)
+    assert out[0, 1] == 0.0 and out[0, 0] == 0.3
+    assert out[1, 0] == 1.0 and out[1, 1] == 0.7
+    assert out[2, 0] == 0.0 and out[2, 1] == 0.0
+    assert out[3, 1] == 2e-9              # untouched: the control
+    assert np.array_equal(out[4], xy[4])
+
+
+# ------------------------------------------------------------------------ 3-D
+
+def test_a_cube_boundary_merges_to_its_six_faces():
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0),
+        cellSize=0.35, regular=False, qdegree=2)
+    verts, tris = _domain_boundary_facets(mesh.dm)
+    mass, planes = _occ_mass_3d(verts, tris)
+    assert len(planes) == 6
+    assert mass == pytest.approx(1.0, rel=1e-12)
+
+
+def test_the_spherical_shell_tool_is_the_polyhedron_not_the_sphere():
+    mesh = uw.meshing.SphericalShell(radiusOuter=1.0, radiusInner=0.5,
+                                     cellSize=0.4, qdegree=2)
+    verts, tris = _domain_boundary_facets(mesh.dm)
+    mesh_volume = _owned_cell_volume(mesh.dm)
+    mass, _planes = _occ_mass_3d(verts, tris)
+    assert mass == pytest.approx(mesh_volume, rel=1e-9)
+    # Negative control, as in 2-D: the discrete shell, not the smooth one.
+    smooth = 4.0 / 3.0 * np.pi * (1.0 ** 3 - 0.5 ** 3)
+    assert abs(mass - smooth) > 1e3 * abs(mass - mesh_volume)
+
+
+def test_the_3d_snap_restores_wall_values_exactly():
+    verts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0],
+                      [0.0, 1.0, 0.0]])
+    tris = np.array([[0, 1, 2], [0, 2, 3]])
+    planes = [(verts[0], np.array([0.0, 0.0, 1.0]))]
+    xyz = np.array([[0.3, 0.4, 1e-10],
+                    [0.6, 0.2, -1e-10],
+                    [0.5, 0.5, 2e-9],
+                    [0.5, 0.5, 0.5]])
+    out = _snap_to_boundary_3d(xyz.copy(), verts, tris, planes)
+    assert out[0, 2] == 0.0 and out[0, 0] == 0.3 and out[0, 1] == 0.4
+    assert out[1, 2] == 0.0
+    assert out[2, 2] == 2e-9              # untouched: the control
+    assert np.array_equal(out[3], xyz[3])

From 9b87c4ea604a8374ba67f65489c5378b617156bc Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 16 Aug 2026 17:57:17 +1000
Subject: [PATCH 2/5] Identify the outcrop band as the boundary TRACE and label
 it
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The band split no longer tests coordinate equality against wall-plane
values: _split_skin_trace takes the domain's own boundary complex and
keeps a skin facet when every vertex and its centroid lie on it — exact
after the post-clip snap, and defined on any boundary. The single-wall
frame the carve/sew still runs in is derived afterwards by
_trace_wall_code, which keeps the established refusals (more than one
wall; and, until the general carve/sew exists, a non-axis-aligned
outcrop).

The band facets — the intersection itself, never a partition of the wall
into named pieces — now carry 'label_trace' alongside the skin and wall
labels, full closures included, in 2-D and 3-D, and the facet count is
reported as info['n_trace_facets']. The outcrop tests assert the trace
coincides with the band, with the interior twin as the negative control.

Underworld development team with AI support from Claude Code
---
 src/underworld3/utilities/place_surface.py | 132 ++++++++++++++-------
 tests/test_0855_place_thin_volume.py       |  32 +++--
 2 files changed, 117 insertions(+), 47 deletions(-)

diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py
index ff8efc18..8a7b416b 100644
--- a/src/underworld3/utilities/place_surface.py
+++ b/src/underworld3/utilities/place_surface.py
@@ -3099,30 +3099,56 @@ def _assembly_skin(points, cells):
     return points[node_ids], skin_local, node_ids
 
 
-def _split_skin_band(skin_xyz, skin_facets, box_lo, box_hi):
-    """Split a skin into its interior part and the wall BAND — the strip of
-    the assembly's clipped face lying exactly in a wall plane (the zone's
-    outcrop). One wall only; more refuses (box-edge bands are not built).
-    Returns ``(interior_idx, band_idx, wall_code)`` with wall_code None
-    when nothing touches a wall.
-
-    Dimension-free: facets are triangles against wall planes in 3-D, edges
-    against wall lines in 2-D.
+def _split_skin_trace(skin_xyz, skin_facets, dom_verts, dom_facets,
+                      tol=1e-9):
+    """Split a skin into its interior part and the boundary TRACE — the
+    facets of the assembly's clipped face lying ON the domain's own
+    boundary facets (the zone's outcrop). Membership is metric but
+    unambiguous: the post-clip snap puts a clipped vertex within rounding
+    of the boundary complex, while an interior skin vertex is a layer
+    mesh-size away. Every vertex AND the centroid must be on the boundary,
+    so a facet spanning between two boundary contacts through the interior
+    cannot pass.
+
+    Dimension-free: edges against boundary edges in 2-D, triangles against
+    boundary triangles in 3-D. Returns ``(interior_idx, trace_idx)``.
     """
-    codes = np.full(len(skin_facets), -1, dtype=np.int64)
+    distance = (_segments_distance if skin_xyz.shape[1] == 2
+                else _sheet_distance)
+    d_vertex = distance(skin_xyz, dom_verts, dom_facets)
+    corners = d_vertex[skin_facets].max(axis=1)
+    d_centre = distance(skin_xyz[skin_facets].mean(axis=1),
+                        dom_verts, dom_facets)
+    on = (corners < tol) & (d_centre < tol)
+    return np.flatnonzero(~on), np.flatnonzero(on)
+
+
+def _trace_wall_code(skin_xyz, skin_facets, trace_idx, box_lo, box_hi):
+    """The single axis-aligned wall the trace lies in, as the legacy
+    ``2*axis + side`` code — the frame the carve and the sew still run
+    against. One wall only; more refuses (box-edge bands are not built),
+    and a trace on a non-axis-aligned boundary refuses until the general
+    carve/sew exists.
+    """
+    codes = np.full(len(trace_idx), -1, dtype=np.int64)
     for axis in range(skin_xyz.shape[1]):
         for side, value in ((0, box_lo[axis]), (1, box_hi[axis])):
-            on = (skin_xyz[skin_facets][:, :, axis] == value).all(axis=1)
+            on = (skin_xyz[skin_facets[trace_idx]][:, :, axis]
+                  == value).all(axis=1)
             codes[on] = 2 * axis + side
-    walls = set(int(c) for c in codes[codes >= 0])
-    if not walls:
-        return np.arange(len(skin_facets)), np.empty(0, dtype=np.int64), None
+    if (codes < 0).any():
+        # TODO(DESIGN): the general carve/sew (arc-ordered splice, cap
+        # meshed on the boundary facets) lifts this; the trace split above
+        # already finds the band on any boundary.
+        raise NotImplementedError(
+            "the zone outcrops on a non-axis-aligned boundary; the "
+            "general carve/sew is not yet built.")
+    walls = set(int(c) for c in codes)
     if len(walls) > 1:
         raise NotImplementedError(
             "the zone meets more than one domain wall; box-edge outcrop "
             "bands are not built.")
-    band = np.flatnonzero(codes >= 0)
-    return np.flatnonzero(codes < 0), band, walls.pop()
+    return walls.pop()
 
 
 def _single_loop(edges, what):
@@ -4263,11 +4289,13 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
     """
     comm = uw.mpi.comm
 
-    # The domain's own boundary, collectively — the clip target. The
-    # axis-aligned box survives alongside it as the frame the outcrop BAND
-    # is identified against (wall lines); a general-boundary band is the
-    # trace-labelling stage, not this one.
-    domain_loops = _domain_loops_2d(dm)
+    # The domain's own boundary, collectively — the clip target and the
+    # complex the outcrop TRACE is identified against. The axis-aligned box
+    # survives alongside as the frame the carve/sew still runs in.
+    dom_verts, dom_edges = _domain_boundary_facets(dm)
+    domain_loops = [dom_verts[loop] for loop in
+                    _skin_loops([(int(a), int(b)) for a, b in dom_edges],
+                                what="the domain boundary")]
     Xb = _coords(dm)
     lo_hi = np.array([Xb.min(axis=0) if len(Xb) else np.full(2, np.inf),
                       -(Xb.max(axis=0)) if len(Xb) else np.full(2, np.inf)])
@@ -4303,11 +4331,15 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
                   for a, b in skin_local]
     loops_asm = _skin_loops(skin_edges)
 
-    # The outcrop band: skin edges lying exactly in ONE wall line. The
+    # The outcrop band: the skin's trace on the domain boundary. The
     # outcropping loop splits into band + interior chain; loops away from
-    # the wall stay holes of the fill.
-    _int_idx, band_idx, wall_code = _split_skin_band(
-        asm_pts, np.asarray(skin_edges, dtype=np.int64), box_lo, box_hi)
+    # the boundary stay holes of the fill.
+    skin_edge_arr = np.asarray(skin_edges, dtype=np.int64)
+    _int_idx, band_idx = _split_skin_trace(asm_pts, skin_edge_arr,
+                                           dom_verts, dom_edges)
+    wall_code = (None if not len(band_idx)
+                 else _trace_wall_code(asm_pts, skin_edge_arr, band_idx,
+                                       box_lo, box_hi))
     open_wall = None
     chain_asm = None
     hole_loops = loops_asm
@@ -4503,7 +4535,8 @@ def mixed(v):
         dm_work, drop_arr, victims_arr, made, placed)
 
     skin_label = label + "_skin"
-    for name in (label, skin_label):
+    trace_label = label + "_trace"
+    for name in (label, skin_label, trace_label):
         if not new_dm.hasLabel(name):
             new_dm.createLabel(name)
     n_zone_local = 0
@@ -4553,17 +4586,23 @@ def mixed(v):
                          int(point_map[corner_r + vS - pStart])]]
             wall_ids += [[int(placed_points[a]), int(placed_points[b])]
                          for a, b in (tuple(p) for p in band_pairs)]
-            for ids in wall_ids:
+            # The first two entries are the splice segments (the 2-D cap);
+            # the rest are the band — the TRACE, the intersection itself,
+            # labelled as such so the model can form whatever unions it
+            # needs (never a partition of the wall into named halves).
+            out_trace = new_dm.getLabel(trace_label)
+            for k, ids in enumerate(wall_ids):
                 joined = new_dm.getFullJoin(ids)
                 if len(joined) != 1:
                     failure = ("an outcrop wall edge is not an edge of the "
                                "sewn mesh; the splice or band was not "
                                "sewn.")
                     break
-                for name, val in pairs:
-                    lab = new_dm.getLabel(name)
-                    for q in new_dm.getTransitiveClosure(int(joined[0]))[0]:
-                        lab.setValue(int(q), int(val))
+                for q in new_dm.getTransitiveClosure(int(joined[0]))[0]:
+                    for name, val in pairs:
+                        new_dm.getLabel(name).setValue(int(q), int(val))
+                    if k >= 2:
+                        out_trace.setValue(int(q), int(label_value))
                 n_wall_local += 1
         failures = comm.allgather(failure)
         real = [f for f in failures if f]
@@ -4620,6 +4659,7 @@ def mixed(v):
     comm.Allreduce(MPI.IN_PLACE, mins, op=MPI.MIN)
     info = {"n_zone_cells": n_zone, "n_skin_faces": n_skin,
             "n_placed": n_placed, "n_removed": n_removed,
+            "n_trace_facets": int(len(band_idx)),
             "min_area": float(mins[0]), "min_angle": float(mins[1])}
     if verbose:
         uw.pprint(f"[place_thin_volume {label!r}] {n_zone} zone cells, "
@@ -4775,9 +4815,12 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1,
     asm_pts, asm_tets = comm.bcast(payload, root=0)
     skin_xyz, skin_tris, skin_node_ids = _assembly_skin(asm_pts, asm_tets)
 
-    # The outcrop band: skin triangles lying exactly in ONE wall plane.
-    interior_idx, band_idx, wall_code = _split_skin_band(
-        skin_xyz, skin_tris, box_lo, box_hi)
+    # The outcrop band: the skin's trace on the domain boundary.
+    interior_idx, band_idx = _split_skin_trace(skin_xyz, skin_tris,
+                                               dom_verts, dom_tris)
+    wall_code = (None if not len(band_idx)
+                 else _trace_wall_code(skin_xyz, skin_tris, band_idx,
+                                       box_lo, box_hi))
     open_wall = None
     band_outline = None
     if wall_code is not None:
@@ -4982,7 +5025,8 @@ def gap_code(v):
     # is invisible to the interface machinery, so the skin faces — which a
     # later placement must hold clear of — go under their own name.
     skin_label = label + "_skin"
-    for name in (label, skin_label):
+    trace_label = label + "_trace"
+    for name in (label, skin_label, trace_label):
         if not new.hasLabel(name):
             new.createLabel(name)
     n_cells_local = 0
@@ -5040,18 +5084,25 @@ def new_id(v):
 
         wall_tris = ([[new_id(int(v)) for v in t] for t in cap_out]
                      if cap_out is not None else [])
+        n_cap_tris = len(wall_tris)
         wall_tris += [[int(placed_new[int(skin_row[int(v)])]) for v in t]
                       for t in skin_tris[band_idx]]
-        for ids in wall_tris:
+        # The cap's annulus first, then the band — the TRACE, the
+        # intersection itself, labelled as such so the model can form
+        # whatever unions it needs (never a partition of the wall into
+        # named pieces).
+        out_trace = new.getLabel(trace_label)
+        for k, ids in enumerate(wall_tris):
             joined = new.getFullJoin(ids)
             if len(joined) != 1:
                 failure = ("an outcrop wall triangle is not a face of the "
                            "sewn mesh; the cap or band was not sewn.")
                 break
-            for name, val in pairs:
-                lab = new.getLabel(name)
-                for q in new.getTransitiveClosure(int(joined[0]))[0]:
-                    lab.setValue(int(q), int(val))
+            for q in new.getTransitiveClosure(int(joined[0]))[0]:
+                for name, val in pairs:
+                    new.getLabel(name).setValue(int(q), int(val))
+                if k >= n_cap_tris:
+                    out_trace.setValue(int(q), int(label_value))
             n_wall_local += 1
     failures = comm.allgather(failure)
     real = [f for f in failures if f]
@@ -5110,6 +5161,7 @@ def new_id(v):
 
     info = {"n_zone_cells": n_zone, "n_skin_faces": n_skin_faces,
             "n_placed": n_placed, "n_removed": n_removed,
+            "n_trace_facets": int(len(band_idx)),
             "min_volume": float(min_vol[0])}
     if verbose:
         uw.pprint(f"[place_thin_volume {label!r}] {info['n_zone_cells']} "
diff --git a/tests/test_0855_place_thin_volume.py b/tests/test_0855_place_thin_volume.py
index 595f9ac0..22bff3b7 100644
--- a/tests/test_0855_place_thin_volume.py
+++ b/tests/test_0855_place_thin_volume.py
@@ -307,8 +307,9 @@ def test_an_outcropping_zone_leaves_a_band_on_the_surface():
                     ).reshape(-1, 3)[: vE - vS]
     top_label = new.getLabel("Top")
     skin_label = new.getLabel("Zone_skin")
+    trace_label = new.getLabel("Zone_trace")
     tv = bounds["Top"].value
-    n_band = 0
+    n_band = n_trace = 0
     for f in range(fS, fE):
         if len(new.getSupport(f)) != 1:
             continue
@@ -318,7 +319,12 @@ def test_an_outcropping_zone_leaves_a_band_on_the_surface():
             assert top_label.getValue(f) == tv
             if skin_label.getValue(f) == zv:
                 n_band += 1
+            if trace_label.getValue(f) == zv:
+                n_trace += 1
     assert n_band > 0, "the zone left no band on the surface"
+    assert n_trace == n_band, (
+        "the trace label does not coincide with the band")
+    assert info["n_trace_facets"] == n_trace
 
     mesh = uw.discretisation.Mesh(
         new, simplex=True, qdegree=3, boundaries=bounds,
@@ -343,17 +349,21 @@ def test_an_outcropping_zone_leaves_a_band_on_the_surface():
 # ------------------------------------------------------------ 2-D outcrop
 
 def _top_edge_census_2d(new, skin_value):
-    """Boundary edges in the top wall line: (total, band, missing Top).
+    """Boundary edges in the top wall line: (total, band, missing Top,
+    trace).
 
     The band is an edge carrying BOTH the zone's skin label and the wall's;
     an edge missing the Top label is a hole the relabel left in the wall.
+    The trace label marks the intersection itself and must coincide with
+    the band.
     """
     fS, fE = new.getHeightStratum(1)
     vS, vE = new.getDepthStratum(0)
     Xn = np.asarray(new.getCoordinatesLocal().array).reshape(-1, 2)[: vE - vS]
     top = new.getLabel("Top")
     skin = new.getLabel("Zone_skin")
-    n_top = n_band = n_bare = 0
+    trace = new.getLabel("Zone_trace")
+    n_top = n_band = n_bare = n_trace = 0
     for f in range(fS, fE):
         if len(new.getSupport(f)) != 1:
             continue
@@ -365,7 +375,9 @@ def _top_edge_census_2d(new, skin_value):
                 n_bare += 1
             if skin.getValue(f) == skin_value:
                 n_band += 1
-    return n_top, n_band, n_bare
+            if trace is not None and trace.getValue(f) == skin_value:
+                n_trace += 1
+    return n_top, n_band, n_bare, n_trace
 
 
 def test_an_outcropping_ribbon_leaves_a_band_on_the_surface():
@@ -388,22 +400,28 @@ def test_an_outcropping_ribbon_leaves_a_band_on_the_surface():
     bounds = base._boundaries_with("Zone")
     zv = bounds["Zone"].value
 
-    interior, _ = place_thin_volume(
+    interior, info0 = place_thin_volume(
         base.dm, [np.array([[0.35, 0.40], [0.60, 0.80]])], width=0.03,
         label="Zone", label_value=zv)
-    n_top0, n_band0, _bare0 = _top_edge_census_2d(interior, zv)
+    n_top0, n_band0, _bare0, n_trace0 = _top_edge_census_2d(interior, zv)
     assert n_top0 > 0 and n_band0 == 0, (
         "the census counted a band on an interior ribbon; it cannot "
         "validate the outcrop")
+    assert n_trace0 == 0 and info0["n_trace_facets"] == 0, (
+        "an interior ribbon carries a trace; the trace label cannot "
+        "validate the outcrop")
 
     before = float(cell_areas(base.dm).sum())
     line = np.array([[0.35, 0.40], [0.60, 1.10]])   # past the top wall
     new, info = place_thin_volume(base.dm, [line], width=0.03,
                                   label="Zone", label_value=zv)
     assert info["n_zone_cells"] > 0
-    n_top, n_band, n_bare = _top_edge_census_2d(new, zv)
+    n_top, n_band, n_bare, n_trace = _top_edge_census_2d(new, zv)
     assert n_band > 0, "the ribbon left no band on the surface"
     assert n_bare == 0, "the relabel left top-wall edges without Top"
+    assert n_trace == n_band, (
+        "the trace label does not coincide with the band")
+    assert info["n_trace_facets"] == n_trace
     assert float(cell_areas(new).sum()) == pytest.approx(before, rel=1e-12)
 
     mesh = uw.discretisation.Mesh(

From 19129589c38bcf7f88e283940c64f6039c1901fd Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 16 Aug 2026 18:10:06 +1000
Subject: [PATCH 3/5] Generalise the 2-D outcrop carve/sew to the mesh's own
 boundary; annulus outcrop lands
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The 2-D outcrop machinery no longer assumes a straight axis-aligned wall:

- The ring splice detects the wall span topologically (the mesh's own
  boundary-edge pairs), and coverage and chain orientation are arc
  projections onto the removed span's polyline, so a curved boundary
  orders exactly as a straight wall did.
- Which boundary vertices the carve may delete is decided by
  _outcrop_frame_2d: a vertex goes when the band re-provides its geometry
  (it lies ON the band) or when it lies on a band-touched straight
  segment, where the splice's cap stays collinear; a compressed corner
  off the band is the domain's shape itself and is protected. On a box
  this reproduces the wall-plane rule exactly.
- One contiguous band only (_refuse_multiple_bands); two ribbons out two
  walls and an arch out the same wall twice keep their refusals.
- The boolean imprints tool corners into the clipped face; a clipped
  corner within 0.1*size of an imprinted domain vertex is merged into it
  (_collapse_boundary_imprints) — measured 0.18-degree sliver restored to
  29.4 degrees against the 29.6 interior twin, with the domain's shape
  and area untouched since the merge moves the band end along the
  boundary segment.

New test: a radial ribbon out of the annulus's outer boundary embeds,
with the interior twin as the negative control; area conserved to 1e-12,
every boundary edge relabelled, trace coincides with the band, minimum
angle floored at 15 degrees, and the P2 oracle exact through the zone.
Box tests unchanged, serial and np=4. The 3-D cap keeps the box frame;
generalising it is the next stage.

Underworld development team with AI support from Claude Code
---
 src/underworld3/utilities/place_surface.py | 259 +++++++++++++++------
 tests/test_0855_place_thin_volume.py       |  88 +++++++
 2 files changed, 281 insertions(+), 66 deletions(-)

diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py
index 8a7b416b..411b01b5 100644
--- a/src/underworld3/utilities/place_surface.py
+++ b/src/underworld3/utilities/place_surface.py
@@ -3123,6 +3123,121 @@ def _split_skin_trace(skin_xyz, skin_facets, dom_verts, dom_facets,
     return np.flatnonzero(~on), np.flatnonzero(on)
 
 
+def _collapse_boundary_imprints(asm_pts, asm_tris, loops, delta):
+    """Merge a clipped-corner node into the domain vertex it almost hits.
+
+    The boolean imprints the tool's corners into the clipped face. Where
+    the assembly's own clipped corner lands within ``delta`` of an
+    imprinted domain vertex, two forced nodes sit a hair apart on the same
+    boundary segment and the layer meshes a sliver spanning them
+    (measured: 0.18 degrees on an annulus radial outcrop, against 29.6 for
+    its interior twin). Taking the domain vertex moves the band's end
+    ALONG the boundary, so the domain's shape and area are untouched; the
+    assembly's side edge tilts by under ``delta`` across its last cell —
+    the end-snap contract, one level up. Runs after the meshed-vs-CAD
+    gate: the collapse is our own exact bookkeeping, not the mesher's.
+    """
+    comp = [_compress_collinear_loop(L) for L in loops]
+    seg_pts, seg_edges = [], []
+    for L in comp:
+        base = len(seg_pts)
+        n = len(L)
+        seg_edges += [(base + i, base + (i + 1) % n) for i in range(n)]
+        seg_pts += list(L)
+    on_bound = (_segments_distance(asm_pts, np.asarray(seg_pts), seg_edges)
+                < 1e-9)
+    remap = np.arange(len(asm_pts))
+    for c in np.vstack(loops):
+        hit = np.flatnonzero((asm_pts == c).all(axis=1))
+        if not len(hit):
+            continue
+        j = int(hit[0])
+        d = np.linalg.norm(asm_pts - c, axis=1)
+        for m in np.flatnonzero((d > 0.0) & (d < delta) & on_bound):
+            remap[int(m)] = j
+    if (remap == np.arange(len(asm_pts))).all():
+        return asm_pts, asm_tris
+    tris = remap[asm_tris]
+    degenerate = ((tris[:, 0] == tris[:, 1]) | (tris[:, 1] == tris[:, 2])
+                  | (tris[:, 2] == tris[:, 0]))
+    tris = tris[~degenerate]
+    used = np.unique(tris)
+    compact = np.full(len(asm_pts), -1, dtype=np.int64)
+    compact[used] = np.arange(len(used))
+    return asm_pts[used], compact[tris]
+
+
+def _refuse_multiple_bands(band_facets):
+    """One contiguous band only: a single carve-and-splice is built.
+
+    Two ribbons out of two walls, or one arch out of the same wall twice,
+    both arrive here as more than one connected component of band facets.
+    """
+    parent = {}
+
+    def find(v):
+        while parent[v] != v:
+            parent[v] = parent[parent[v]]
+            v = parent[v]
+        return v
+
+    for facet in band_facets:
+        for v in facet:
+            parent.setdefault(int(v), int(v))
+        roots = {find(int(v)) for v in facet}
+        anchor = roots.pop()
+        for r in roots:
+            parent[r] = anchor
+    if len({find(v) for v in parent}) > 1:
+        raise NotImplementedError(
+            "the zone meets the domain boundary in more than one band; a "
+            "multiply-outcropping zone is not built.")
+
+
+def _outcrop_frame_2d(comp_loops, asm_pts, band_pairs, X, on_wall):
+    """Which mesh boundary vertices an outcrop may DELETE, and which ring
+    vertices count as boundary contact AWAY from the outcrop.
+
+    A boundary vertex may go when the band re-provides its geometry (the
+    vertex lies ON the band — the assembly holds a node at its exact
+    position), or when it lies on a band-touched straight segment of the
+    compressed boundary, where the splice's cap segments stay collinear
+    with the segment and the domain shape is preserved. A compressed
+    CORNER off the band is the domain's shape itself and is never deleted.
+    On a box this reproduces the wall-plane rule exactly: the touched
+    segment is the whole wall, and the box corners are the protected ones.
+
+    Returns ``(deletable, near_touched)`` — masks over the mesh vertices.
+    """
+    band_nodes = sorted({int(v) for pair in band_pairs for v in pair})
+    B = asm_pts[band_nodes]
+    seg_pts, seg_edges, corners = [], [], []
+    for L in comp_loops:
+        corners.append(L)
+        n = len(L)
+        for i in range(n):
+            A, C = L[i], L[(i + 1) % n]
+            e = C - A
+            u = np.clip(((B - A) @ e) / float(e @ e), 0.0, 1.0)
+            d = np.linalg.norm(B - (A + u[:, None] * e), axis=1)
+            if (d < 1e-9).any():
+                seg_edges.append((len(seg_pts), len(seg_pts) + 1))
+                seg_pts += [A, C]
+    if not seg_edges:
+        raise RuntimeError(
+            "the outcrop band lies on no boundary segment; the clip and "
+            "the boundary disagree")
+    d_seg = _segments_distance(X, np.asarray(seg_pts), seg_edges)
+    d_band = _segments_distance(X, asm_pts,
+                                [tuple(p) for p in band_pairs])
+    is_corner = np.zeros(len(X), dtype=bool)
+    for c in np.vstack(corners):
+        is_corner |= np.linalg.norm(X - c, axis=1) < 1e-9
+    near_touched = d_seg < 1e-9
+    deletable = on_wall & near_touched & (~is_corner | (d_band < 1e-9))
+    return deletable, near_touched
+
+
 def _trace_wall_code(skin_xyz, skin_facets, trace_idx, box_lo, box_hi):
     """The single axis-aligned wall the trace lies in, as the legacy
     ``2*axis + side`` code — the frame the carve and the sew still run
@@ -3572,27 +3687,47 @@ def _outcrop_chain_2d(loops, band_pairs):
     return chain, holes
 
 
-def _outcrop_ring_splice(ring, X, open_wall, chain_ids, chain_t):
-    """Replace the cavity ring's wall span with the ribbon's interior chain.
-
-    The raw ring of an outcropping carve runs ALONG the open wall through
-    vertices about to be deleted. The fill's boundary instead descends
-    around the ribbon: the ring's single contiguous run of wall edges is
-    removed and the chain is spliced between the run's surviving end
-    vertices, oriented to meet them. The two splice segments are the 2-D
-    cap — what remains of the 3-D wall annulus one dimension down.
-
-    ``chain_ids`` are the chain's rows in the fill's combined numbering,
-    ``chain_t`` the along-wall coordinates of its two ends. Returns
-    ``(spliced_ring, removed_wall_pairs)``, the removed pairs as old vertex
-    rows for the wall-label discovery. A second wall run refuses — the
-    carve spilled onto the wall away from the outcrop.
+def _arc_project(polyline, p):
+    """``(distance, arc length)`` of the closest point on an open polyline."""
+    best_d, best_s = np.inf, 0.0
+    s0 = 0.0
+    for i in range(len(polyline) - 1):
+        A, B = polyline[i], polyline[i + 1]
+        e = B - A
+        length = float(np.linalg.norm(e))
+        u = float(np.clip((p - A) @ e / (length * length), 0.0, 1.0))
+        d = float(np.linalg.norm(p - (A + u * e)))
+        if d < best_d:
+            best_d, best_s = d, s0 + u * length
+        s0 += length
+    return best_d, best_s
+
+
+def _outcrop_ring_splice(ring, X, boundary_pairs, chain_ids, chain_ends):
+    """Replace the cavity ring's boundary span with the ribbon's interior
+    chain.
+
+    The raw ring of an outcropping carve runs ALONG the domain boundary
+    through vertices about to be deleted. The fill's boundary instead
+    descends around the ribbon: the ring's single contiguous run of
+    boundary edges is removed and the chain is spliced between the run's
+    surviving end vertices, oriented to meet them. The two splice segments
+    are the 2-D cap — what remains of the 3-D wall annulus one dimension
+    down.
+
+    ``boundary_pairs`` is the mesh's own boundary-edge vertex pairs, so
+    the span is a topological run, not a coordinate test.  ``chain_ids``
+    are the chain's rows in the fill's combined numbering, ``chain_ends``
+    the coordinates of its two end nodes; coverage and orientation are
+    arc projections onto the removed span's own polyline, so a curved
+    boundary orders exactly as a straight wall did. Returns
+    ``(spliced_ring, removed_wall_pairs)``, the removed pairs as old
+    vertex rows for the wall-label discovery. A second boundary run
+    refuses — the carve spilled onto the boundary away from the outcrop.
     """
-    axis, value = open_wall
-    t = 1 - axis
     n = len(ring)
-    onw = [X[v][axis] == value for v in ring]
-    wall_edge = [onw[i] and onw[(i + 1) % n] for i in range(n)]
+    wall_edge = [frozenset((ring[i], ring[(i + 1) % n])) in boundary_pairs
+                 for i in range(n)]
     if not any(wall_edge):
         raise RuntimeError(
             "the outcrop band does not meet the cavity's wall span; raise "
@@ -3604,20 +3739,24 @@ def _outcrop_ring_splice(ring, X, open_wall, chain_ids, chain_t):
             "reduce `clearance` or move the zone off the wall.")
     i0 = next(i for i in range(n) if wall_edge[i] and not wall_edge[i - 1])
     k = sum(wall_edge)
-    corner_l, corner_r = ring[i0], ring[(i0 + k) % n]
-    lo_t, hi_t = sorted((float(X[corner_l][t]), float(X[corner_r][t])))
-    if not (lo_t < min(chain_t) and max(chain_t) < hi_t):
+    removed = [(ring[(i0 + j) % n], ring[(i0 + j + 1) % n])
+               for j in range(k)]
+    span = X[np.asarray([ring[(i0 + j) % n] for j in range(k + 1)])]
+    ends = [_arc_project(span, np.asarray(c, dtype=float))
+            for c in chain_ends]
+    strictly_inside = all(
+        d < 1e-9
+        and np.linalg.norm(c - span[0]) > 1e-9
+        and np.linalg.norm(c - span[-1]) > 1e-9
+        for (d, _s), c in zip(ends, chain_ends))
+    if not strictly_inside:
         raise RuntimeError(
             "the cavity's wall span does not cover the outcrop band; raise "
             "`clearance`.")
-    removed = [(ring[(i0 + j) % n], ring[(i0 + j + 1) % n])
-               for j in range(k)]
     # The surviving arc, corner_r around to corner_l, then the chain with
     # its corner_l-side end first — the loop's orientation is preserved.
     tail = [ring[(i0 + k + j) % n] for j in range(n - k + 1)]
-    tl = float(X[corner_l][t])
-    seq = (chain_ids if abs(tl - chain_t[0]) <= abs(tl - chain_t[1])
-           else chain_ids[::-1])
+    seq = chain_ids if ends[0][1] <= ends[1][1] else chain_ids[::-1]
     return tail + list(seq), removed
 
 
@@ -4290,17 +4429,11 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
     comm = uw.mpi.comm
 
     # The domain's own boundary, collectively — the clip target and the
-    # complex the outcrop TRACE is identified against. The axis-aligned box
-    # survives alongside as the frame the carve/sew still runs in.
+    # complex the outcrop TRACE is identified against.
     dom_verts, dom_edges = _domain_boundary_facets(dm)
     domain_loops = [dom_verts[loop] for loop in
                     _skin_loops([(int(a), int(b)) for a, b in dom_edges],
                                 what="the domain boundary")]
-    Xb = _coords(dm)
-    lo_hi = np.array([Xb.min(axis=0) if len(Xb) else np.full(2, np.inf),
-                      -(Xb.max(axis=0)) if len(Xb) else np.full(2, np.inf)])
-    comm.Allreduce(MPI.IN_PLACE, lo_hi, op=MPI.MIN)
-    box_lo, box_hi = lo_hi[0], -lo_hi[1]
 
     failure = None
     payload = None
@@ -4317,6 +4450,8 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
                     f"the ribbon assembly meshed to area {mesh_area:.12e} "
                     f"against CAD {cad_area:.12e}; the layer mesh does not "
                     "fill its own outlines.")
+            asm_pts, asm_tris = _collapse_boundary_imprints(
+                asm_pts, asm_tris, domain_loops, 0.1 * size)
             payload = (asm_pts, asm_tris)
         except Exception as exc:
             failure = f"{type(exc).__name__}: {exc}"
@@ -4337,19 +4472,14 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
     skin_edge_arr = np.asarray(skin_edges, dtype=np.int64)
     _int_idx, band_idx = _split_skin_trace(asm_pts, skin_edge_arr,
                                            dom_verts, dom_edges)
-    wall_code = (None if not len(band_idx)
-                 else _trace_wall_code(asm_pts, skin_edge_arr, band_idx,
-                                       box_lo, box_hi))
-    open_wall = None
+    outcropping = bool(len(band_idx))
     chain_asm = None
     hole_loops = loops_asm
     band_pairs = set()
-    if wall_code is not None:
-        w_axis, w_side = divmod(int(wall_code), 2)
-        open_wall = (w_axis,
-                     float(box_hi[w_axis] if w_side else box_lo[w_axis]))
+    if outcropping:
+        _refuse_multiple_bands(skin_edge_arr[band_idx])
         band_pairs = {frozenset((int(a), int(b)))
-                      for a, b in np.asarray(skin_edges)[band_idx]}
+                      for a, b in skin_edge_arr[band_idx]}
         chain_asm, hole_loops = _outcrop_chain_2d(loops_asm, band_pairs)
 
     vS, vE = dm.getDepthStratum(0)
@@ -4396,16 +4526,12 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
             beside_held = np.zeros(len(X), dtype=bool)
             beside_held[cells[held_c].ravel()] = True
             deletable_wall = np.zeros(len(X), dtype=bool)
-            if open_wall is not None:
-                # An outcrop deletes wall vertices — in the open wall's
-                # line only. A vertex in a second wall line is a domain
-                # corner and stays protected: deleting it for one wall
-                # would breach the other.
-                deletable_wall = on_wall & (X[:, open_wall[0]]
-                                            == open_wall[1])
-                other = 1 - open_wall[0]
-                for v2 in (box_lo[other], box_hi[other]):
-                    deletable_wall &= ~(X[:, other] == v2)
+            near_touched = None
+            if outcropping:
+                comp_loops = [_compress_collinear_loop(L)
+                              for L in domain_loops]
+                deletable_wall, near_touched = _outcrop_frame_2d(
+                    comp_loops, asm_pts, band_pairs, X, on_wall)
             protected = (on_wall & ~deletable_wall) | held_v | beside_held
             victim = (d_skin < reach_v) & ~protected
 
@@ -4427,25 +4553,26 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
 
             ring, drop = _ring_growing(cells, drop, held_c)
             removed_wall = []
-            if open_wall is None:
+            if not outcropping:
                 if on_wall[np.asarray(ring)].any():
                     raise RuntimeError(
                         "the ribbon's cavity reached the domain wall; the "
                         "volume must be interior, with clearance to spare")
             else:
+                boundary_pairs = {frozenset((a, b))
+                                  for _e, a, b in _boundary_edges(dm_work)}
                 chain_ids = [len(X) + int(v) for v in chain_asm]
-                t_ax = 1 - open_wall[0]
-                chain_t = (float(asm_pts[chain_asm[0]][t_ax]),
-                           float(asm_pts[chain_asm[-1]][t_ax]))
+                chain_ends = (asm_pts[chain_asm[0]],
+                              asm_pts[chain_asm[-1]])
                 ring, removed_wall = _outcrop_ring_splice(
-                    ring, X, open_wall, chain_ids, chain_t)
-                off_open = [v for v in ring if v < len(X)
-                            and on_wall[v]
-                            and X[v][open_wall[0]] != open_wall[1]]
-                if off_open:
+                    ring, X, boundary_pairs, chain_ids, chain_ends)
+                off_band = [v for v in ring if v < len(X)
+                            and on_wall[v] and not near_touched[v]]
+                if off_band:
                     raise RuntimeError(
-                        "the ribbon's cavity reached a second domain "
-                        "wall; only a one-wall outcrop is built.")
+                        "the ribbon's cavity reached the domain boundary "
+                        "away from the outcrop band; only a one-band "
+                        "outcrop is built.")
             old_ring = np.asarray([v for v in ring if v < len(X)])
             if victim[old_ring].any():
                 raise RuntimeError(
@@ -4503,7 +4630,7 @@ def mixed(v):
                     "place_thin_volume internal: the gathered region "
                     "touches a shared point; the gather mask under-reached.")
             outcrop = None
-            if open_wall is not None:
+            if outcropping:
                 # The splice ends: the surviving corner vertices and the
                 # chain ends they meet — the chain is the ring's tail, in
                 # the orientation the splice chose.
@@ -4572,7 +4699,7 @@ def mixed(v):
     # and the band itself — is relabelled EXPLICITLY with what the deleted
     # wall span carried, full closures included.
     n_wall_local = 0
-    if open_wall is not None:
+    if outcropping:
         pairs = comm.bcast(outcrop[0] if comm.rank == target else None,
                            root=target)
         for name, val in pairs:
diff --git a/tests/test_0855_place_thin_volume.py b/tests/test_0855_place_thin_volume.py
index 22bff3b7..ec0745ff 100644
--- a/tests/test_0855_place_thin_volume.py
+++ b/tests/test_0855_place_thin_volume.py
@@ -445,6 +445,94 @@ def test_an_outcropping_ribbon_leaves_a_band_on_the_surface():
         f"= {float(err.max()):.3e}")
 
 
+def _annulus_boundary_census(new, zone_value):
+    """Boundary edges of the annulus: (total, trace, missing a wall label).
+
+    An edge missing both Lower and Upper is a hole the relabel left in the
+    boundary; the trace marks the outcrop band.
+    """
+    fS, fE = new.getHeightStratum(1)
+    lower = new.getLabel("Lower")
+    upper = new.getLabel("Upper")
+    trace = new.getLabel("Zone_trace")
+    n_bound = n_trace = n_bare = 0
+    for f in range(fS, fE):
+        if len(new.getSupport(f)) != 1:
+            continue
+        n_bound += 1
+        if lower.getValue(f) < 0 and upper.getValue(f) < 0:
+            n_bare += 1
+        if trace.getValue(f) == zone_value:
+            n_trace += 1
+    return n_bound, n_trace, n_bare
+
+
+def test_an_outcropping_ribbon_on_the_annulus_leaves_a_trace():
+    """The general boundary: a radial ribbon out of the annulus's outer
+    boundary. There is no circle the annulus is failing to be — the clip
+    must land on the mesh's own chords, so the domain area is conserved to
+    round-off, the trace edges carry the Upper label as well as the trace,
+    and the P2 oracle stays exact through the zone.
+
+    The negative control comes first: an interior twin of the same ribbon
+    carries no trace, so the trace counted on the outcrop is produced by
+    the outcrop.
+    """
+    import sympy
+    from underworld3.utilities.line_cut import cell_areas
+
+    base = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5,
+                              cellSize=0.06, qdegree=3)
+    bounds = base._boundaries_with("Zone")
+    zv = bounds["Zone"].value
+    theta = 0.4
+    ray = np.array([np.cos(theta), np.sin(theta)])
+
+    interior, info0 = place_thin_volume(
+        base.dm, [np.array([0.62 * ray, 0.88 * ray])], width=0.03,
+        label="Zone", label_value=zv)
+    _nb0, n_trace0, _bare0 = _annulus_boundary_census(interior, zv)
+    assert info0["n_trace_facets"] == 0 and n_trace0 == 0, (
+        "an interior ribbon carries a trace; the census cannot validate "
+        "the outcrop")
+
+    before = float(cell_areas(base.dm).sum())
+    line = np.array([0.62 * ray, 1.30 * ray])       # past the outer boundary
+    new, info = place_thin_volume(base.dm, [line], width=0.03,
+                                  label="Zone", label_value=zv)
+    assert info["n_zone_cells"] > 0
+    assert info["n_trace_facets"] > 0, "the ribbon left no trace"
+    n_bound, n_trace, n_bare = _annulus_boundary_census(new, zv)
+    assert n_trace == info["n_trace_facets"]
+    assert n_bare == 0, "the relabel left boundary edges without a label"
+    assert float(cell_areas(new).sum()) == pytest.approx(before, rel=1e-12)
+    # A clipped corner landing near an imprinted domain vertex once meshed
+    # a 0.18-degree sliver; the imprint collapse restores the interior
+    # twin's quality (measured 29.4 against 29.6 degrees). The floor is
+    # far above any sliver and far below the healthy fill.
+    assert info["min_angle"] > 15.0
+
+    mesh = uw.discretisation.Mesh(
+        new, simplex=True, qdegree=3, boundaries=bounds,
+        coordinate_system_type=base.CoordinateSystem.coordinate_type)
+    x, y = mesh.X
+    exact = x**2 + y**2
+    t = uw.discretisation.MeshVariable("T_ann", mesh, 1, degree=2)
+    poisson = uw.systems.Poisson(mesh, u_Field=t)
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = 1.0
+    poisson.f = -4.0
+    for wall in ("Lower", "Upper"):
+        poisson.add_dirichlet_bc(sympy.Matrix([exact]), wall)
+    poisson.tolerance = 1e-11
+    poisson.solve()
+    X = np.asarray(t.coords)
+    err = np.abs(np.asarray(t.data[:, 0]) - (X[:, 0]**2 + X[:, 1]**2))
+    assert float(err.max()) < 1e-8, (
+        f"the outcropped annulus assembles a wrong operator: "
+        f"max |u - exact| = {float(err.max()):.3e}")
+
+
 def test_a_ribbon_stopping_short_of_the_wall_still_refuses():
     """No band, no outcrop: the interior contract keeps its refusal.
 

From ef5dbd21292d9729fa51657f7794f2943aab62f6 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 16 Aug 2026 18:15:45 +1000
Subject: [PATCH 4/5] Parallel coverage for the general-boundary outcrop; trace
 count is partition-independent
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The annulus outcrop runs through the same gather-first mechanism at
np=2, 3 and 4, and the info dict — trace count included — is identical
at np=1, 2 and 4 (measured: trace=1, zone=44, skin=24, placed=60,
removed=17 at every rank count), because the boundary complex is
gathered from unshared support-1 facets and sorted into a canonical
order before the clip.

The new ptest's boundary census skips SHARED facets, not just SF
leaves: a partition-seam facet also has local support 1 on its owner
side, and a leaves-only probe reports phantom bare boundary edges.

Underworld development team with AI support from Claude Code
---
 .../ptest_0855_place_thin_volume_parallel.py  | 52 +++++++++++++++++++
 1 file changed, 52 insertions(+)

diff --git a/tests/parallel/ptest_0855_place_thin_volume_parallel.py b/tests/parallel/ptest_0855_place_thin_volume_parallel.py
index bac9f6ec..8dfd7f71 100644
--- a/tests/parallel/ptest_0855_place_thin_volume_parallel.py
+++ b/tests/parallel/ptest_0855_place_thin_volume_parallel.py
@@ -135,6 +135,58 @@ def test_a_2d_outcropping_ribbon_embeds_in_parallel():
     assert n_bare == 0, "the relabel left top-wall edges without Top"
 
 
+def test_a_2d_outcrop_on_the_annulus_embeds_in_parallel():
+    """The general-boundary outcrop at np>=2, on the annulus.
+
+    The domain loops are gathered from UNSHARED support-1 edges — a
+    partition seam also has local support 1, and taking it would eat the
+    mesh at the seam differently at every rank count — so the clip, the
+    trace and the splice are functions of the mesh, not of the partition.
+    The gathered complex is sorted into a canonical order, so the info
+    dict (trace count included) is identical at any rank count; the
+    np=1 value is asserted equal by running the serial suite.
+    """
+    comm = uw.mpi.comm
+    mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5,
+                              cellSize=0.06, qdegree=2)
+    theta = 0.4
+    ray = np.array([np.cos(theta), np.sin(theta)])
+    line = np.array([0.62 * ray, 1.30 * ray])
+    new, info = place_thin_volume(mesh.dm, [line], width=0.03,
+                                  label="Zone", label_value=5)
+    assert info["n_zone_cells"] > 0
+    assert info["n_trace_facets"] > 0, "the ribbon left no trace"
+    assert _owned_label_count(new, "Zone", 5) == info["n_zone_cells"]
+    assert _owned_label_count(new, "Zone_skin", 5) == info["n_skin_faces"]
+    gathered = comm.allgather(info)
+    assert all(g == gathered[0] for g in gathered)
+
+    # Every true boundary edge still carries a wall label, and the trace
+    # edges among them number exactly what the info advertises. A domain
+    # boundary edge is support 1 AND unshared — a partition-seam edge also
+    # has local support 1, and counting it would report phantom bare edges
+    # (the _true_wall_vertex_mask distinction, met here by the probe).
+    from underworld3.utilities.place_surface import _shared_point_flags
+    pStart, _pEnd = new.getChart()
+    shared = _shared_point_flags(new).astype(bool)
+    fS, fE = new.getHeightStratum(1)
+    lower = new.getLabel("Lower")
+    upper = new.getLabel("Upper")
+    trace = new.getLabel("Zone_trace")
+    n_trace = n_bare = 0
+    for f in range(fS, fE):
+        if shared[f - pStart] or len(new.getSupport(f)) != 1:
+            continue
+        if lower.getValue(f) < 0 and upper.getValue(f) < 0:
+            n_bare += 1
+        if trace.getValue(f) == 5:
+            n_trace += 1
+    n_trace = int(comm.allreduce(n_trace, op=MPI.SUM))
+    n_bare = int(comm.allreduce(n_bare, op=MPI.SUM))
+    assert n_trace == info["n_trace_facets"]
+    assert n_bare == 0, "the relabel left boundary edges without a label"
+
+
 def test_2d_refusals_are_collective():
     """A ribbon stopping short of the wall refuses IDENTICALLY everywhere."""
     comm = uw.mpi.comm

From 04cca8ed3d636a3ac867b3c46027f21238d27d1b Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 16 Aug 2026 19:32:18 +1000
Subject: [PATCH 5/5] A ribbon out through a box corner embeds; wall labels
 restored per removed segment
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The general carve/sew subsumes the box-corner outcrop the wall-code
frame used to refuse: the corner vertex lies ON the band, so it is
deletable and the band re-provides its exact position. What the corner
exposed was the label harvest, not the geometry — labels common to ALL
removed wall edges are empty across a corner (Top one side, Right the
other), leaving the new wall edges bare. The harvest is now per removed
edge, and each splice/band edge takes the labels of the removed segment
it lies on, so the corner node itself recovers both walls' labels
through the two closures. Uniform spans (one wall, the annulus loop)
are unchanged.

Pinned by a new test: the diagonal corner ribbon embeds with area
conserved to 1e-12, the trace present on both walls, no bare wall edge,
and minimum angle above 15 degrees. Two separate bands stay refused.

Underworld development team with AI support from Claude Code
---
 src/underworld3/utilities/place_surface.py | 74 ++++++++++++++--------
 tests/test_0855_place_thin_volume.py       | 50 +++++++++++++++
 2 files changed, 99 insertions(+), 25 deletions(-)

diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py
index 411b01b5..e285dddb 100644
--- a/src/underworld3/utilities/place_surface.py
+++ b/src/underworld3/utilities/place_surface.py
@@ -4588,26 +4588,29 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value,
                     "volume must be interior, with clearance to spare")
             victim |= orphan
 
-            # Labels the deleted wall span carried, read before the
-            # rebuild forgets them — the splice segments and the band are
-            # relabelled with exactly these after sewing.
-            wall_pairs = []
+            # Labels the deleted wall span carried, read PER EDGE before
+            # the rebuild forgets them — a span across a domain corner
+            # changes label mid-way (Top on one side, Right on the other),
+            # so each new wall edge must take the labels of the old
+            # segment it lies on, not a set common to the whole span.
+            span_labels = []
+            span_segments = []
             if removed_wall:
-                edge_pts = [int(dm_work.getFullJoin(
-                    [int(a) + vS, int(b) + vS])[0])
-                    for a, b in removed_wall]
-                for i in range(dm_work.getNumLabels()):
-                    name = dm_work.getLabelName(i)
-                    if name in reconnect._TOPOLOGY_LABELS:
-                        continue
-                    lab = dm_work.getLabel(name)
-                    values = lab.getValueIS()
-                    if values is None:
-                        continue
-                    for val in values.getIndices():
-                        if all(lab.getValue(p) == int(val)
-                               for p in edge_pts):
-                            wall_pairs.append((name, int(val)))
+                names = [dm_work.getLabelName(i)
+                         for i in range(dm_work.getNumLabels())
+                         if dm_work.getLabelName(i)
+                         not in reconnect._TOPOLOGY_LABELS]
+                for a, b in removed_wall:
+                    p = int(dm_work.getFullJoin([int(a) + vS,
+                                                 int(b) + vS])[0])
+                    pairs_p = []
+                    for name in names:
+                        val = dm_work.getLabel(name).getValue(p)
+                        if val >= 0:
+                            pairs_p.append((name, int(val)))
+                    span_labels.append(pairs_p)
+                    span_segments.append((X[int(a)].copy(),
+                                          X[int(b)].copy()))
 
             Xall = np.vstack([X, asm_pts])
             holes = [[len(X) + int(v) for v in loop] for loop in hole_loops]
@@ -4634,7 +4637,7 @@ def mixed(v):
                 # The splice ends: the surviving corner vertices and the
                 # chain ends they meet — the chain is the ring's tail, in
                 # the orientation the splice chose.
-                outcrop = (wall_pairs,
+                outcrop = (span_labels, span_segments,
                            int(removed_wall[0][0]),
                            int(removed_wall[-1][1]),
                            int(ring[-len(chain_asm)]) - len(X),
@@ -4700,23 +4703,44 @@ def mixed(v):
     # wall span carried, full closures included.
     n_wall_local = 0
     if outcropping:
-        pairs = comm.bcast(outcrop[0] if comm.rank == target else None,
-                           root=target)
-        for name, val in pairs:
+        label_names = comm.bcast(
+            sorted({name for pairs_p in outcrop[0] for name, _v in pairs_p})
+            if comm.rank == target else None, root=target)
+        for name in label_names:
             if not new_dm.hasLabel(name):
                 new_dm.createLabel(name)
         if comm.rank == target:
-            _wp, corner_l, corner_r, a_first, a_last = outcrop
+            (span_lab, span_seg, corner_l, corner_r,
+             a_first, a_last) = outcrop
             wall_ids = [[int(point_map[corner_l + vS - pStart]),
                          int(placed_points[a_first])],
                         [int(placed_points[a_last]),
                          int(point_map[corner_r + vS - pStart])]]
             wall_ids += [[int(placed_points[a]), int(placed_points[b])]
                          for a, b in (tuple(p) for p in band_pairs)]
+            mids = [0.5 * (X[corner_l] + asm_pts[a_first]),
+                    0.5 * (asm_pts[a_last] + X[corner_r])]
+            mids += [0.5 * (asm_pts[a] + asm_pts[b])
+                     for a, b in (tuple(p) for p in band_pairs)]
+
+            def span_labels_at(m):
+                best, at = np.inf, 0
+                for k2, (A, B) in enumerate(span_seg):
+                    e = B - A
+                    u = np.clip(float((m - A) @ e) / float(e @ e),
+                                0.0, 1.0)
+                    d = float(np.linalg.norm(m - (A + u * e)))
+                    if d < best:
+                        best, at = d, k2
+                return span_lab[at]
+
             # The first two entries are the splice segments (the 2-D cap);
             # the rest are the band — the TRACE, the intersection itself,
             # labelled as such so the model can form whatever unions it
             # needs (never a partition of the wall into named halves).
+            # Each edge takes the labels of the removed segment it lies
+            # on, so a band across a domain corner restores Top on one
+            # side and Right on the other.
             out_trace = new_dm.getLabel(trace_label)
             for k, ids in enumerate(wall_ids):
                 joined = new_dm.getFullJoin(ids)
@@ -4726,7 +4750,7 @@ def mixed(v):
                                "sewn.")
                     break
                 for q in new_dm.getTransitiveClosure(int(joined[0]))[0]:
-                    for name, val in pairs:
+                    for name, val in span_labels_at(mids[k]):
                         new_dm.getLabel(name).setValue(int(q), int(val))
                     if k >= 2:
                         out_trace.setValue(int(q), int(label_value))
diff --git a/tests/test_0855_place_thin_volume.py b/tests/test_0855_place_thin_volume.py
index ec0745ff..c076e70d 100644
--- a/tests/test_0855_place_thin_volume.py
+++ b/tests/test_0855_place_thin_volume.py
@@ -533,6 +533,56 @@ def test_an_outcropping_ribbon_on_the_annulus_leaves_a_trace():
         f"max |u - exact| = {float(err.max()):.3e}")
 
 
+def test_a_ribbon_out_through_a_box_corner_embeds():
+    """A single ribbon exiting diagonally through the (1,1) corner.
+
+    The wall-code frame used to refuse this as \"more than one domain
+    wall\"; the general carve/sew subsumes it. The corner vertex lies ON
+    the band, so it is deletable and the band re-provides its exact
+    position — one contiguous band across the corner, area conserved,
+    both walls' labels restored. Two SEPARATE bands stay refused
+    (the two-ribbon and arch tests below).
+    """
+    from underworld3.utilities.line_cut import cell_areas
+
+    base = _box2(0.05)
+    bounds = base._boundaries_with("Zone")
+    zv = bounds["Zone"].value
+    before = float(cell_areas(base.dm).sum())
+    corner = np.array([[0.6, 0.6], [1.1, 1.1]])
+    new, info = place_thin_volume(base.dm, [corner], width=0.03,
+                                  label="Zone", label_value=zv)
+    assert info["n_zone_cells"] > 0
+    assert info["n_trace_facets"] > 0, "the ribbon left no trace"
+    assert info["min_angle"] > 15.0
+    assert float(cell_areas(new).sum()) == pytest.approx(before, rel=1e-12)
+
+    # The trace edges carry a wall label each — Top or Right — and no
+    # boundary edge in either wall line is left bare.
+    fS, fE = new.getHeightStratum(1)
+    vS, vE = new.getDepthStratum(0)
+    Xn = np.asarray(new.getCoordinatesLocal().array).reshape(-1, 2)[: vE - vS]
+    top = new.getLabel("Top")
+    right = new.getLabel("Right")
+    trace = new.getLabel("Zone_trace")
+    n_trace = n_bare = 0
+    for f in range(fS, fE):
+        if len(new.getSupport(f)) != 1:
+            continue
+        verts = [int(q) - vS for q in new.getTransitiveClosure(f)[0]
+                 if vS <= int(q) < vE]
+        on_top = all(Xn[v][1] == 1.0 for v in verts)
+        on_right = all(Xn[v][0] == 1.0 for v in verts)
+        if not (on_top or on_right):
+            continue
+        if top.getValue(f) < 0 and right.getValue(f) < 0:
+            n_bare += 1
+        if trace.getValue(f) == zv:
+            n_trace += 1
+    assert n_trace == info["n_trace_facets"]
+    assert n_bare == 0, "the relabel left a wall edge without its label"
+
+
 def test_a_ribbon_stopping_short_of_the_wall_still_refuses():
     """No band, no outcrop: the interior contract keeps its refusal.