From df95bdcb3dd472132f61f96905ccb1f58fb7b54a Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 7 Jun 2026 10:57:24 +0100
Subject: [PATCH 1/6] docs: boundary tangent-slip strategy design

Design for a mesh-orchestrated, boundary-surface-owned tangent-slip + restore
contract. Each boundary gets a BoundingSurface object (kind radial/plane/facet/
free) that owns its geometry, state flags, and tangent_project/restore/release
methods; the mesh orchestrates cross-surface vertex classification. Replaces the
slip logic duplicated across the spring/ma/mmpde/ot movers and the runtime
_is_radial_coords heuristic. mesh.boundaries (the persisted gmsh/DMPlex labels)
is left untouched; surfaces live in a new mesh.bounding_surfaces collection.

This is the step-1 (interface) doc; implementation follows in this branch toward
a development PR. Step 2 (movers consume the API) stays on the mover feature
branch.

Underworld development team with AI support from Claude Code
---
 .../design/boundary-slip-strategy.md          | 329 ++++++++++++++++++
 1 file changed, 329 insertions(+)
 create mode 100644 docs/developer/design/boundary-slip-strategy.md

diff --git a/docs/developer/design/boundary-slip-strategy.md b/docs/developer/design/boundary-slip-strategy.md
new file mode 100644
index 00000000..bbb60036
--- /dev/null
+++ b/docs/developer/design/boundary-slip-strategy.md
@@ -0,0 +1,329 @@
+# Boundary-slip strategy: a mesh-owned tangent-slip + surface-restore contract
+
+```{note}
+**Status:** design proposal (2026-06-06). This document specifies a refactor;
+no behaviour change is intended in the first (interface) step. It gates making
+the MMPDE mover the production default — see
+`docs/developer/design/anisotropic-mmpde-mover.md` and the mover-API
+simplification work.
+```
+
+## Motivation
+
+Every metric mover that moves *boundary* vertices needs **tangential boundary
+slip**: a boundary node may slide *along* the domain surface (so the mesh can
+concentrate cells where a feature meets the wall) but must not drift *off* it
+(which would change the domain and, on a free-slip Stokes problem, leak
+`v·n ≠ 0`). Slip is therefore two operations applied to the moved boundary
+coordinates:
+
+1. **tangent-project** — remove the boundary-normal component of a node's
+   displacement, so it slides in the tangent plane;
+2. **restore-to-surface** — snap the node back onto the true surface, because
+   tangent-projecting a *finite* step off a *curved* surface leaves it a
+   sagitta inside the chord; without restoration nodes creep inward over many
+   iterations.
+
+Today this logic is **duplicated and geometry-coupled**:
+
+- The most evolved version lives in `meshing/_ot_adapt.py` as private helpers
+  (`_resolve_slip`, `_build_slip_projector`, `_slip_normals`,
+  `_is_radial_coords`, `_boundary_centre`, `_nearest_on_facets_2d/3d`). The
+  `mmpde` and `ot` movers consume it through `_build_slip_projector`.
+- The `spring` and `ma` movers carry their **own inline** per-ring radial snap
+  (`meshing/smoothing.py`, the `boundary_slip and is_bnd.any()` blocks) — a
+  special case of the same idea, written before the `_ot_adapt` unification.
+- Whether a boundary is "radial" (snap to `|r|`) vs "flat" (no snap needed) is
+  decided at **run time** by a `CoordinateSystemType` heuristic
+  (`_is_radial_coords`), and the snap **centre** is recovered every call by a
+  parallel `allreduce` over boundary coordinates — even though the constructor
+  *knew* the exact centre and radii.
+- Free / deforming surfaces (where the surface position is itself unknown) are
+  handled by a `dict {label: snap_bool}` opt-out that disables restoration —
+  a placeholder for a real kinematic restore that does not yet exist.
+
+This is fragile (four code paths to keep parallel-consistent), guesses geometry
+it could be told, and offers no clean seam for the free-surface case. Making
+MMPDE the default mover means *every* analytic-boundary mesh must slip
+correctly with no per-script `boundary_slip='ring'`/`'box'` fiddling. That
+needs a **single, mesh-owned slip contract**.
+
+## Proposal
+
+Introduce a **first-class bounding-surface object** that owns a boundary's
+geometry, state flags, and the methods to operate on it (tangent-project,
+restore, normals). Slip is then a *bounding-surface-level* capability — the
+surface knows whether it is radial, planar, free, or generic, and how to restore
+a point to itself — and the **mesh is only the orchestrator** of the
+cross-surface concerns the movers need (which vertices slip vs pin, including
+junction vertices shared by two surfaces, and composing the per-surface
+operations into one pass).
+
+This is the key correction over a mesh-level provider dict: behaviour and state
+live on the surface, not in a side table keyed by label. If one surface is
+radial and another is free, each carries its own flags and the orchestrator
+calls the correct per-surface restore for each.
+
+```{important}
+**`mesh.boundaries` is NOT repurposed.** That `Enum` carries the gmsh-style
+boundary *labels* that persist all the way through: gmsh `.msh` → PETSc DMPlex
+`DMLabel`s → the mesh hdf5 checkpoint. It is a persistence contract, not a
+scratch object — silently turning it into rich objects would risk the
+label round-trip. The bounding-surface objects live in a **separate**
+collection (`mesh.bounding_surfaces`); each *references* a boundary label by
+name but does not replace it. Reusing the `mesh.boundaries` name for the rich
+objects would require a deliberate deprecation plan and is out of scope here.
+```
+
+### Bounding-surface objects
+
+A `BoundingSurface` binds a boundary *label* (a name in `mesh.boundaries`) to
+that surface's geometry, state, and slip methods:
+
+```python
+class BoundingSurface:
+    label: str                # the gmsh/DMPlex boundary label ("Upper", "Lower", ...)
+    kind: str                 # "radial" | "plane" | "facet" | "free"
+    is_free: bool             # True ⇒ restore follows the live surface, not a fixed target
+    # geometry it was told at construction (kind-dependent):
+    #   radial: centre, radius     plane: point, normal     facet: reference_facets
+
+    def normals(self, coords): ...          # outward unit normals on THIS surface (Gamma_P1 restricted)
+    def tangent_project(self, coords): ...  # remove this surface's normal component
+    def restore(self, coords): ...          # snap back onto THIS surface (kind-specific)
+    def release(self): ...                  # flip rigid → free (free-surface activation)
+```
+
+The surface object *is* the restoration capability (the earlier
+`TangentSlipProvider` folds into it). `restore` dispatches on `kind`:
+
+| `kind` | meshes / origin | `restore` |
+|--------|-----------------|-----------|
+| `radial` | `Annulus`, `SphericalShell`, `CubedSphere`, cylinders | exact `|r|` snap about the **known** centre: `centre + r̂·radius`. Concave-safe (no chord sag). Inner/outer surfaces are separate objects with their own radius. |
+| `plane` | `UnstructuredSimplexBox`, `StructuredQuadBox` faces | zero the off-plane coordinate (axis-aligned ⇒ tangent-project already keeps it on the face; supports non-axis-aligned planes too). |
+| `facet` | loaded-from-file / internal boundaries with no analytic form | nearest point on the surface's reference facets (the current `_nearest_on_facets_*` fallback). Convex-safe; concave bias is a documented TODO. |
+| `free` | free-surface module, OR any surface that has been `release()`-d | follow the **current deformed discrete surface** (generic surface tangent), not a fixed target. `is_free=True`. |
+
+The object **encapsulates the geometry the constructor already knows**, so the
+runtime `_is_radial_coords` guess and the per-call centre `allreduce` disappear:
+`Annulus(radiusOuter, radiusInner, centre)` builds a `radial` surface on
+`"Upper"` (radius `radiusOuter`) and on `"Lower"` (radius `radiusInner`)
+directly.
+
+### Public API
+
+The mesh orchestrates the per-surface objects. The low-level contract is the
+`(is_pinned, project)` pair the movers already consume, plus an in-place
+convenience and a public registration hook (decisions confirmed in review,
+2026-06-06):
+
+```python
+# Low-level: the projector tuple the movers expect. is_pinned drives the
+# solve's pinned-DOF set; project() does tangent-project + restore per surface.
+is_pinned, project = mesh.boundary_slip(slip_spec, reference_coords=X0)
+Y = project(Y_moved)
+
+# In-place convenience (built on the above) for callers that just want coords
+# snapped back — e.g. a checkpoint reload, a diagnostic, the free-surface module.
+mesh.project_to_slip_surface(coords, slip_spec, reference_coords=X0)
+
+# Public registration so users can install a custom analytic surface object
+# (e.g. an ellipsoid) that the constructors don't know about.
+mesh.register_tangent_slip_provider(label, surface)
+```
+
+`slip_spec` keeps the back-compatible forms already accepted by
+`_resolve_slip`: `True`/`"ring"`/`"box"`/`"all"` (all named codim-1
+boundaries), a label name, a list of labels, or a `dict {label: snap_bool}`
+(`False` = free surface, slip but do not restore).
+
+The two primitives are surface-level, exposed at mesh level for convenience
+over a label (and reusable by checkpoints, diagnostics, the free-surface
+module):
+
+```python
+mesh.tangent_project(coords, labels)    # = surface.tangent_project per label
+mesh.restore_to_surface(coords, label)  # = surface.restore for that label
+```
+
+`tangent_project` is **geometry-agnostic** — it uses the projected P1
+boundary-normal field (`mesh.Gamma_P1`), which is already smooth and
+consistently oriented on curved boundaries where raw face normals are noisy.
+`restore` is **geometry-specific** and dispatches on the surface's `kind`.
+
+### Restoration is a surface-*state* question, not just initial geometry
+
+A subtlety the surface model must get right (raised in review, 2026-06-06):
+**the correct restoration target depends on the surface's current state, not
+only on how the mesh was built.** The motivating case is a free surface on a
+sphere/annulus. At construction the outer surface is rigid at `|r| = radius`, so
+the `radial` restore is exact. But once that surface is `release()`-d to deform
+(free-surface dynamics), the analytic radius is *no longer the surface* —
+snapping returned nodes to the original `|r|` would actively fight the surface
+motion. The restore must fall back to a **generic surface tangent**: keep
+returned nodes on the *current* discrete surface (the deformed boundary facets),
+not on the frozen analytic shape.
+
+So each surface object is **state-aware** along two axes:
+
+1. **Mode (the `kind`/`is_free` flag on the object).** A surface is *rigid*
+   (snap to an analytic target — `radial`/`plane`) or *free* (follow the current
+   deformed surface — the generic surface-tangent / `facet`-style restore
+   against the live boundary, `is_free=True`). `release()` flips a `radial`/
+   `plane` surface to `free` **in place on the object** — no side-table to keep
+   in sync. This is the same generic-surface mechanism as the loaded-from-file
+   `facet` fallback, reached by a *state transition* rather than by mesh type.
+2. **Reference.** A rigid surface restores relative to a captured reference
+   (`reference_coords` per adapt) for idempotence; a free surface's "reference"
+   is the *current* surface, re-read each call from live mesh state.
+
+Implication for the contract: a surface's `restore` must be able to read
+**current mesh state** (its live facets / the coordinate field), not just
+constants frozen at construction. The analytic kinds are the constant-target
+special case; `free`/`facet` are the live-target general case behind the same
+method. The free-surface case is therefore a *mode of the same surface object*,
+not a parallel code path — and "one surface radial, one free" is just two
+objects with different flags, which the orchestrator handles per surface.
+
+### Mesh-side data model
+
+`mesh.bounding_surfaces` is a **new** collection of `BoundingSurface` objects,
+keyed by boundary label — *separate from* and *additional to* `mesh.boundaries`
+(the persistent gmsh/DMPlex label `Enum`, left untouched). The objects ARE the
+slip registry — behaviour and state live on them, not in a side table:
+
+```python
+mesh.boundaries                          # unchanged: gmsh/DMPlex labels (persisted)
+mesh.bounding_surfaces["Upper"].kind     # "radial"
+mesh.bounding_surfaces["Upper"].release()  # → "free" (free-surface activation)
+```
+
+Constructors that know their geometry build `radial`/`plane` surfaces for their
+labels. A label with no analytic object (mesh loaded from file, an internal
+boundary, third-party mesh) defaults to a `facet` surface built from the current
+boundary facets — exactly today's geometry-general path, so nothing regresses.
+The surfaces are reconstructed on a checkpoint reload from the stored
+construction parameters (a loaded mesh otherwise gets only the `facet` default);
+because the labels themselves persist in `mesh.boundaries`/DMPlex, the
+surface-to-label binding is always recoverable.
+`register_tangent_slip_provider(label, surface)` lets a user install a custom
+surface object (e.g. an ellipsoid) the constructors don't know about.
+
+### Slip-vs-pin classification (unchanged, made canonical)
+
+A boundary vertex **slips iff it belongs to exactly one slip surface**.
+Vertices on a non-slip boundary (count 0), at a **junction** of two slip
+surfaces (count ≥2 — a box corner, where the normal is ambiguous), or with a
+degenerate/non-finite projected normal are **pinned**. This is the
+label-driven rule already in `_build_slip_projector` (it fixed an older
+topology classifier that spuriously pinned coarse-but-smooth curved rings); the
+refactor promotes it to the one canonical implementation.
+
+## Why this is the clean-land gate
+
+- **Mmpde-as-default needs zero per-script slip config.** With surface objects
+  built at construction, `boundary_slip=True` on any analytic-boundary mesh
+  "just works" — no `'ring'` vs `'box'` choice, no coordinate-type guess.
+- **One parallel-safe code path.** Today four movers must each keep the
+  owned-vertex projection + halo sync + collective centre consistent; the
+  refactor leaves exactly one.
+- **A real seam for free surfaces.** A surface's `free` mode is where the
+  kinematic restore lands later, behind the *same* `mesh.boundary_slip`
+  interface the movers already call — no second integration.
+
+## Migration plan (interface first, per the branching discipline)
+
+Following `docs/developer/guides/branching-strategy.md` (extract the interface,
+land it on `development`, keep the feature branch to implementation):
+
+1. **Bounding-surface objects + mesh API, no behaviour change** (lands on
+   `development`):
+   - Add `mesh.bounding_surfaces`: a **new** collection of `BoundingSurface`
+     objects keyed by boundary label, with `kind`/`is_free` + the
+     `normals`/`tangent_project`/`restore`/`release` methods. **Leave
+     `mesh.boundaries` (the persisted gmsh/DMPlex label `Enum`) untouched.**
+   - Add `mesh.boundary_slip` (orchestrator), `mesh.tangent_project`,
+     `mesh.restore_to_surface`, `mesh.project_to_slip_surface`, and
+     `register_tangent_slip_provider`.
+   - Build `radial`/`plane` surfaces in the `Annulus`, `SphericalShell`,
+     `CubedSphere`, and box constructors.
+   - Internally, `mesh.boundary_slip` *delegates to the existing `_ot_adapt`
+     helpers*, with the per-surface object supplying the geometry instead of
+     `_is_radial_coords` + the centre `allreduce`. Verify bit-identical mesh
+     trajectories before/after (the `radial` restore is already what
+     `_build_slip_projector` does for radial meshes; the only change is *where
+     the centre/radius come from*).
+2. **Movers consume the public API**: replace the `_build_slip_projector` call
+   in `mmpde`/`ot` and the **inline** radial-snap blocks in `spring`/`ma` with
+   `mesh.boundary_slip(...)`. Delete the four private duplicates. Re-validate
+   the convection harness (serial + np=5) for bit-identical behaviour.
+3. **Follow-up (separate work)**: the `facet` concave-bias cure (mean-preserving
+   / smoothness constraint — the documented TODO) and the `free`-surface restore
+   mode (`release()` + live-surface follow) for the free-surface + adaptive-mesh
+   case (cf. `project_freesurface_ale_design`).
+
+## Invariants the refactor must preserve
+
+- **Parallel safety.** Projection touches only owned vertices; the caller (the
+  mover) halo-syncs. Any collective (e.g. a free-surface global reduction) must
+  run unconditionally on every rank — the analytic surfaces remove the only
+  current collective (the centre `allreduce`) because the centre is a known
+  constant.
+- **DM-stale safety.** `mesh.Gamma_P1` (the `_projected_normals` MeshVariable)
+  must exist *before* a mover builds its solver DM (creating that MeshVariable
+  mid-mover stales the DM handle — see `project_uw3_smoother_footguns`).
+  Building surface objects / pre-touching `Gamma_P1` at construction makes this
+  automatic instead of relying on `_resolve_slip`'s pre-touch.
+- **Free surfaces** (`surface.is_free`, or a `dict` `False` value) slip without
+  restoration.
+- **Reference surface.** A rigid surface's restore is relative to a fixed
+  reference (`reference_coords`, captured once per adapt) so repeated projection
+  is idempotent and does not accumulate drift.
+- **Units.** Surfaces store radii/points in the mesh's coordinate units.
+
+## Decisions (review, 2026-06-06)
+
+- **Locus.** Slip behaviour + state live on **bounding-surface objects** in a
+  new `mesh.bounding_surfaces` collection, not a mesh-level provider table and
+  **not** by repurposing `mesh.boundaries` (the persisted gmsh/DMPlex labels,
+  left untouched). The mesh orchestrates cross-surface concerns (vertex
+  slip-vs-pin classification, junctions, composition).
+- **API surface.** Low-level `(is_pinned, project)` (the movers' contract) **plus**
+  an in-place `mesh.project_to_slip_surface(coords, spec)` convenience built on
+  it.
+- **User-registerable surfaces.** Yes —
+  `mesh.register_tangent_slip_provider(label, surface)` is public, so users can
+  install a custom surface object (e.g. an ellipsoid).
+- **Land plan.** Step 1 (bounding-surface objects + mesh API, internally
+  delegating to the existing helpers, bit-identical) lands on `development` as
+  its own PR; the mover-side swap (step 2) stays on the feature branch — clean
+  API/impl split.
+
+## Still open (for review)
+
+1. **Naming.** Method `boundary_slip` (matches the existing kwarg) vs
+   `slip_project` vs `tangent_slip`; the surface class `BoundingSurface` and its
+   collection `mesh.bounding_surfaces` (chosen to keep `mesh.boundaries` safe);
+   the registration method `register_tangent_slip_provider` (now installs a
+   surface object — reconcile vs `register_bounding_surface`?); the `kind`
+   values `radial`/`plane`/`facet`/`free`.
+2. **Concave non-analytic restore.** Defer the `facet`-kind concave-bias cure to
+   the follow-up, or block on it? (Radial surfaces take the exact analytic
+   branch and are immune; the bias only bites a concave *non-analytic* surface,
+   which no current production case hits.)
+
+## Deferred cases (handle after the simple analytic geometries)
+
+- **Geographic meshes are an odd case** (flagged in review, 2026-06-06). The
+  `GEOGRAPHIC` coordinate system mixes a radial (depth) direction with
+  lon/lat surface coordinates and a non-trivial metric; "tangent to the
+  surface" and "the projected normal" are not the plain Cartesian operations the
+  `radial`/`Gamma_P1` path assumes. **Do not** try to cover it in the first cut
+  — get `Annulus`/`SphericalShell`/box working under the new contract first, then
+  design a `geographic` surface `kind` against the settled interface. Until then
+  geographic meshes keep the current `_is_radial_coords` → radial-snap behaviour
+  via the fallback path (they classify as radial today).
+- **Free-surface restoration** (the surface `free` mode and the `release()`
+  rigid→free transition described above) is the primary follow-up
+  (cf. `project_freesurface_ale_design`).
+```

From 85e3152ab56003af280bbbd69b201dc462e2a217 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 7 Jun 2026 12:41:53 +0100
Subject: [PATCH 2/6] feat(meshing): BoundingSurface objects +
 mesh.boundary_slip API (step 1)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Additive, self-contained boundary tangent-slip API per
docs/developer/design/boundary-slip-strategy.md. No behaviour change: the
development movers do not call the new API, so nothing existing is affected.

- New meshing/bounding_surface.py: BoundingSurface (kind radial/plane/facet/
  free) owns a boundary's geometry, slip state, and normals/tangent_project/
  restore/release methods. radial restore = exact |r| snap about a known
  centre; plane restore = orthogonal projection onto a face; release() flips a
  rigid surface to free (restore becomes a no-op — the free-surface follow). +
  register_radial_surfaces / register_box_face_surfaces constructor helpers.
- Mesh gains a SEPARATE mesh.bounding_surfaces collection (mesh.boundaries — the
  persisted gmsh/DMPlex labels — is left untouched), plus the orchestrator
  mesh.boundary_slip (returns (is_pinned, project)), the per-surface delegates
  mesh.tangent_project / restore_to_surface, the in-place
  project_to_slip_surface convenience, and register_tangent_slip_provider.
  Slip-vs-pin is label-driven: a vertex slips iff it lies on exactly one
  registered analytic surface; junctions / unregistered / degenerate-normal
  vertices are pinned (the step-1 safe default; facet restore is a follow-up).
- Annulus / SphericalShell / CubedSphere register radial surfaces (Upper/Lower
  at the known radii about the origin); UnstructuredSimplexBox / StructuredQuadBox
  register plane surfaces for their axis-aligned faces.
- tests/test_0762_bounding_surfaces.py (11 tests): constructor registration,
  mesh.boundaries untouched, radial/plane restore, release, register + type
  check, slip orchestrator keeps slipped vertices on the boundary while leaving
  interior alone, pin fallback when no surface is registered.

The mover-side swap to this API (step 2) lands on the mover feature branch,
which has the unified projector; the development movers are untouched here.

Underworld development team with AI support from Claude Code
---
 .../design/boundary-slip-strategy.md          |  52 +++--
 .../discretisation/discretisation_mesh.py     | 169 ++++++++++++++
 src/underworld3/meshing/annulus.py            |   7 +
 src/underworld3/meshing/bounding_surface.py   | 210 ++++++++++++++++++
 src/underworld3/meshing/cartesian.py          |  10 +
 src/underworld3/meshing/spherical.py          |  14 ++
 tests/test_0762_bounding_surfaces.py          | 162 ++++++++++++++
 7 files changed, 607 insertions(+), 17 deletions(-)
 create mode 100644 src/underworld3/meshing/bounding_surface.py
 create mode 100644 tests/test_0762_bounding_surfaces.py

diff --git a/docs/developer/design/boundary-slip-strategy.md b/docs/developer/design/boundary-slip-strategy.md
index bbb60036..85f04a41 100644
--- a/docs/developer/design/boundary-slip-strategy.md
+++ b/docs/developer/design/boundary-slip-strategy.md
@@ -236,8 +236,21 @@ refactor promotes it to the one canonical implementation.
 Following `docs/developer/guides/branching-strategy.md` (extract the interface,
 land it on `development`, keep the feature branch to implementation):
 
-1. **Bounding-surface objects + mesh API, no behaviour change** (lands on
-   `development`):
+```{note}
+**Implementation reality (2026-06-07).** The unified slip projector
+(`_build_slip_projector`, `_resolve_slip`, `_gamma_p1_at_vertices`,
+`_nearest_on_facets_*`) lives on the mover **feature branch**, not on
+`development`. On `development` `_ot_adapt.py` has only the *primitive* helpers
+(`_slip_normals`, `_boundary_centre`, `_is_radial_coords`) and the movers use
+their older inline slip. So step 1 on `development` is implemented as a
+**self-contained, additive** API (it does not depend on the feature-branch
+projector), and step 2's bit-identical claim is interpreted as
+*machine-precision identical* — the analytic centre differs from the feature
+branch's boundary-COM `allreduce` only at round-off.
+```
+
+1. **Bounding-surface objects + mesh API, additive — no behaviour change**
+   (lands on `development`):
    - Add `mesh.bounding_surfaces`: a **new** collection of `BoundingSurface`
      objects keyed by boundary label, with `kind`/`is_free` + the
      `normals`/`tangent_project`/`restore`/`release` methods. **Leave
@@ -246,21 +259,26 @@ land it on `development`, keep the feature branch to implementation):
      `mesh.restore_to_surface`, `mesh.project_to_slip_surface`, and
      `register_tangent_slip_provider`.
    - Build `radial`/`plane` surfaces in the `Annulus`, `SphericalShell`,
-     `CubedSphere`, and box constructors.
-   - Internally, `mesh.boundary_slip` *delegates to the existing `_ot_adapt`
-     helpers*, with the per-surface object supplying the geometry instead of
-     `_is_radial_coords` + the centre `allreduce`. Verify bit-identical mesh
-     trajectories before/after (the `radial` restore is already what
-     `_build_slip_projector` does for radial meshes; the only change is *where
-     the centre/radius come from*).
-2. **Movers consume the public API**: replace the `_build_slip_projector` call
-   in `mmpde`/`ot` and the **inline** radial-snap blocks in `spring`/`ma` with
-   `mesh.boundary_slip(...)`. Delete the four private duplicates. Re-validate
-   the convection harness (serial + np=5) for bit-identical behaviour.
-3. **Follow-up (separate work)**: the `facet` concave-bias cure (mean-preserving
-   / smoothness constraint — the documented TODO) and the `free`-surface restore
-   mode (`release()` + live-surface follow) for the free-surface + adaptive-mesh
-   case (cf. `project_freesurface_ale_design`).
+     `CubedSphere`, and box constructors. **Self-contained**: `radial`/`plane`
+     `restore` are direct (analytic centre/radius, plane projection); `normals`
+     reuses the primitive `_slip_normals` (Gamma_P1); a label with **no**
+     analytic surface is **pinned** (safe default) — `facet` restore is a
+     follow-up. No dependence on the feature-branch projector.
+   - Because the `development` movers do not call the new API, step 1 cannot
+     change any existing trajectory (additive). Validate with **new tests**:
+     constructors register the right `kind`/geometry; `radial.restore` lands
+     points on `|r|`, `plane.restore` on the face; `release()` flips to `free`;
+     `mesh.boundaries` is unchanged.
+2. **Movers consume the public API** (on the mover feature branch, which has the
+   unified projector): replace the `_build_slip_projector` call in `mmpde`/`ot`
+   and the **inline** radial-snap blocks in `spring`/`ma` with
+   `mesh.boundary_slip(...)`. Delete the private duplicates. Re-validate the
+   convection harness (serial + np=5) for *machine-precision-identical*
+   behaviour (the centre source changes COM → analytic).
+3. **Follow-up (separate work)**: the `facet` restore + its concave-bias cure
+   (mean-preserving / smoothness constraint — the documented TODO) and the
+   `free`-surface restore mode (`release()` + live-surface follow) for the
+   free-surface + adaptive-mesh case (cf. `project_freesurface_ale_design`).
 
 ## Invariants the refactor must preserve
 
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index da6020f1..23c5dd3f 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -442,6 +442,11 @@ class replacement_boundaries(Enum):
 
         self.filename = filename
         self.boundaries = boundaries
+        # Bounding-surface objects (tangent-slip + restore), keyed by boundary
+        # label. SEPARATE from self.boundaries (the persisted gmsh/DMPlex
+        # labelling, untouched). Populated by analytic-geometry constructors;
+        # see docs/developer/design/boundary-slip-strategy.md.
+        self._bounding_surfaces = {}
         self.boundary_normals = boundary_normals
         self.regions = regions
         self.parent = None       # Set by extract_region() for submeshes
@@ -2007,6 +2012,170 @@ def Gamma_P1(self):
             self._update_projected_normals()
         return self._projected_normals.sym
 
+    # ===================================================================
+    #  Bounding surfaces — per-surface tangent-slip + restore.
+    #  See docs/developer/design/boundary-slip-strategy.md. SEPARATE from
+    #  self.boundaries (the persisted gmsh/DMPlex labelling, untouched).
+    # ===================================================================
+    @property
+    def bounding_surfaces(self):
+        """Mapping ``{label: BoundingSurface}`` of this mesh's registered
+        bounding-surface objects (tangent-slip + restore).
+
+        This is a NEW collection, separate from and additional to
+        :attr:`boundaries` (the persisted gmsh/DMPlex label ``Enum``, left
+        untouched). Populated by analytic-geometry constructors (Annulus,
+        SphericalShell, CubedSphere, box meshes); user-extendable via
+        :meth:`register_tangent_slip_provider`.
+        """
+        if not hasattr(self, "_bounding_surfaces") or self._bounding_surfaces is None:
+            self._bounding_surfaces = {}
+        return self._bounding_surfaces
+
+    def register_tangent_slip_provider(self, label, surface):
+        """Install a :class:`BoundingSurface` object for a boundary ``label``
+        (separate from the persisted ``mesh.boundaries`` labelling).
+
+        Lets a user declare a custom analytic surface (e.g. an ellipsoid) the
+        constructors don't know about, or replace one (free-surface release).
+        """
+        from underworld3.meshing.bounding_surface import BoundingSurface
+        if not isinstance(surface, BoundingSurface):
+            raise TypeError(
+                "surface must be a BoundingSurface; got "
+                f"{type(surface).__name__}")
+        self.bounding_surfaces[str(label)] = surface
+        return surface
+
+    def _resolve_slip_spec(self, slip_spec):
+        """Resolve a ``slip_spec`` to ``(slip_labels tuple, free_labels set)``.
+
+        Back-compatible forms: ``True``/``"all"``/``"ring"``/``"box"`` → all
+        geometric boundary labels; ``False``/``None`` → none; a label name; a
+        list of labels; a ``dict {label: snap_bool}`` (``False`` = free
+        surface, slip but do not restore).
+        """
+        from underworld3.meshing.smoothing import _auto_pinned_labels
+        geo = _auto_pinned_labels(self)
+        if slip_spec is None or slip_spec is False:
+            return (), set()
+        if slip_spec is True:
+            return tuple(geo), set()
+        if isinstance(slip_spec, dict):
+            return tuple(slip_spec.keys()), {k for k, v in slip_spec.items() if not v}
+        if isinstance(slip_spec, str):
+            s = slip_spec.strip().lower()
+            if s in ("true", "on", "1", "all", "ring", "box", "axes", "axis"):
+                return tuple(geo), set()
+            if s in ("false", "off", "0", "none", ""):
+                return (), set()
+            return (slip_spec,), set()
+        return tuple(slip_spec), set()
+
+    def restore_to_surface(self, coords, label):
+        """Snap ``coords`` onto the named bounding surface (delegates to the
+        surface object's ``restore``)."""
+        return self.bounding_surfaces[str(label)].restore(
+            numpy.asarray(coords, dtype=float))
+
+    def tangent_project(self, coords, label, reference):
+        """Tangent-slide ``coords`` (displacement measured from ``reference``)
+        on the named bounding surface (delegates to the surface object)."""
+        return self.bounding_surfaces[str(label)].tangent_project(
+            numpy.asarray(coords, dtype=float),
+            numpy.asarray(reference, dtype=float))
+
+    def boundary_slip(self, slip_spec=True, reference_coords=None,
+                      boundary_labels=None):
+        """Build ``(is_pinned, project)`` for tangent slip on this mesh's
+        bounding surfaces — the orchestrator the metric movers call.
+
+        See ``docs/developer/design/boundary-slip-strategy.md``. The mesh
+        classifies which vertices slip vs pin (the cross-surface concern); each
+        surface object owns its tangent-project + restore.
+
+        A vertex **slips** iff it lies on **exactly one** slip surface that has a
+        registered analytic :class:`BoundingSurface`; non-boundary, junction
+        (≥2 surfaces), unregistered-surface, or degenerate-normal vertices are
+        **pinned** (the step-1 safe default — ``facet`` restore is a follow-up).
+
+        Parameters
+        ----------
+        slip_spec :
+            See :meth:`_resolve_slip_spec`. Default ``True`` (all surfaces).
+        reference_coords : ndarray, optional
+            Fixed reference vertex positions (local-chart vertex order) that
+            displacements are measured from. Defaults to ``mesh.X.coords``.
+        boundary_labels : iterable of str, optional
+            Boundary labels defining the boundary (``is_bnd``). Defaults to all
+            geometric boundary labels; pass a mover's pinned set for parity.
+
+        Returns
+        -------
+        (is_pinned, project) : (ndarray[bool], callable)
+            ``is_pinned`` is the per-vertex pinned mask (local-chart order);
+            ``project(Y)`` slides+restores the slip vertices of ``Y`` in place
+            and returns it.
+        """
+        from underworld3.meshing.smoothing import _pinned_mask, _auto_pinned_labels
+
+        dm = self.dm
+        pStart, pEnd = dm.getDepthStratum(0)
+        n_verts = pEnd - pStart
+        if reference_coords is None:
+            reference_coords = numpy.asarray(self.X.coords, dtype=float)
+        ref = numpy.asarray(reference_coords, dtype=float)
+
+        all_labels = (tuple(boundary_labels) if boundary_labels is not None
+                      else _auto_pinned_labels(self))
+        is_bnd = _pinned_mask(dm, all_labels)
+
+        slip_labels, free_labels = self._resolve_slip_spec(slip_spec)
+        surf = self.bounding_surfaces
+        # Only labels with a registered analytic surface can slip (step 1).
+        usable = [lab for lab in slip_labels if lab in surf]
+        masks = {lab: _pinned_mask(dm, (lab,)) for lab in usable}
+        count = numpy.zeros(n_verts, dtype=int)
+        for m in masks.values():
+            count += m.astype(int)
+        slip_mask = is_bnd & (count == 1)
+        is_pinned = is_bnd & ~slip_mask
+        vert_label = numpy.empty(n_verts, dtype=object)
+        for lab, m in masks.items():
+            vert_label[m & slip_mask] = lab
+
+        slip_b = numpy.nonzero(slip_mask)[0]
+        if slip_b.size == 0:
+            return is_pinned, (lambda Y: Y)
+        old_slip = ref[slip_b]
+        labels_b = vert_label[slip_b]
+
+        def project(Y):
+            Y = numpy.asarray(Y, dtype=float)
+            for lab in usable:
+                sel = labels_b == lab
+                if not sel.any():
+                    continue
+                bs = surf[lab]
+                idx = slip_b[sel]
+                slid = bs.tangent_project(Y[idx], old_slip[sel])
+                # FREE surfaces (dict spec False) slide but do not restore.
+                Y[idx] = slid if lab in free_labels else bs.restore(slid)
+            return Y
+
+        return is_pinned, project
+
+    def project_to_slip_surface(self, coords, slip_spec=True,
+                                reference_coords=None, boundary_labels=None):
+        """In-place convenience over :meth:`boundary_slip`: slide+restore the
+        slip vertices of ``coords`` (a full local-chart vertex array) and return
+        it. For callers that just want coordinates snapped back (a checkpoint
+        reload, a diagnostic, the free-surface module)."""
+        _is_pinned, project = self.boundary_slip(
+            slip_spec, reference_coords=reference_coords,
+            boundary_labels=boundary_labels)
+        return project(numpy.asarray(coords, dtype=float))
+
     @timing.routine_timer_decorator
     def update_lvec(self):
         """
diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py
index b495b72f..a1f47080 100644
--- a/src/underworld3/meshing/annulus.py
+++ b/src/underworld3/meshing/annulus.py
@@ -542,6 +542,13 @@ class boundary_normals(Enum):
     x, y = new_mesh.X
     new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])]
 
+    # Bounding-surface objects for tangent slip: Upper/Lower are radial about the
+    # origin at the known radii (see docs/developer/design/boundary-slip-strategy.md).
+    from underworld3.meshing.bounding_surface import register_radial_surfaces
+    register_radial_surfaces(
+        new_mesh, centre=(0.0, 0.0),
+        label_radius={"Upper": radiusOuter, "Lower": radiusInner})
+
     return new_mesh
 
 
diff --git a/src/underworld3/meshing/bounding_surface.py b/src/underworld3/meshing/bounding_surface.py
new file mode 100644
index 00000000..927a0f79
--- /dev/null
+++ b/src/underworld3/meshing/bounding_surface.py
@@ -0,0 +1,210 @@
+"""Boundary-surface objects: per-surface tangent-slip + restore.
+
+See ``docs/developer/design/boundary-slip-strategy.md``.
+
+A :class:`BoundingSurface` binds a boundary *label* (a name in
+``mesh.boundaries`` — the persisted gmsh/DMPlex labelling, which is **not**
+replaced here) to that surface's geometry, slip state, and the tangent-project /
+restore operations for it. Surfaces live in the separate
+``mesh.bounding_surfaces`` collection.
+
+This is the step-1 (additive, self-contained) implementation: ``radial`` and
+``plane`` restore are analytic; ``facet`` (nearest reference facet) and the
+``free`` live-surface follow are follow-ups (the orchestrator pins labels with
+no analytic surface). The module depends only on the primitive
+``_ot_adapt`` helpers that exist on ``development`` (``_slip_normals``), not on
+the mover feature branch's unified projector.
+"""
+import numpy as np
+
+import underworld3 as uw
+
+_VALID_KINDS = ("radial", "plane", "facet", "free")
+
+
+def _unit(v):
+    v = np.asarray(v, dtype=float).ravel()
+    n = np.linalg.norm(v)
+    return v / n if n > 0 else v
+
+
+def _as_float(x):
+    """Coerce a radius/length that may be a UWQuantity / pint Quantity to a
+    bare non-dimensional float in the mesh's coordinate units."""
+    try:
+        return float(x)
+    except (TypeError, ValueError):
+        nd = uw.non_dimensionalise(x)
+        return float(getattr(nd, "magnitude", getattr(nd, "value", nd)))
+
+
+class BoundingSurface:
+    """One bounding surface of a mesh — its geometry, slip state, and methods.
+
+    Parameters
+    ----------
+    mesh : underworld3 mesh
+        The mesh this surface belongs to.
+    label : str
+        The boundary label (a name in ``mesh.boundaries``) this surface is for.
+    kind : {"radial", "plane", "facet", "free"}
+        ``radial`` — snap to ``|r| = radius`` about ``centre`` (annulus / sphere
+        / cylinder); ``plane`` — project onto the plane through ``point`` with
+        unit ``normal`` (box face); ``facet`` — nearest reference facet
+        (follow-up); ``free`` — follow the live deformed surface, analytic
+        restore is a no-op (follow-up).
+    centre, radius : for ``radial``.
+    point, normal : for ``plane``.
+    is_free : bool
+        If True (or ``kind == "free"``), :meth:`restore` is a no-op — the
+        surface follows the live discrete boundary rather than a fixed target.
+    """
+
+    def __init__(self, mesh, label, kind, *, centre=None, radius=None,
+                 point=None, normal=None, is_free=False):
+        if kind not in _VALID_KINDS:
+            raise ValueError(
+                f"BoundingSurface kind must be one of {_VALID_KINDS}; got {kind!r}")
+        self._mesh = mesh
+        self.label = str(label)
+        self.kind = kind
+        self.is_free = bool(is_free) or kind == "free"
+        self.centre = None if centre is None else np.asarray(centre, dtype=float).ravel()
+        self.radius = None if radius is None else _as_float(radius)
+        self.point = None if point is None else np.asarray(point, dtype=float).ravel()
+        self.normal = None if normal is None else _unit(normal)
+        if kind == "radial" and (self.centre is None or self.radius is None):
+            raise ValueError("radial BoundingSurface requires centre and radius")
+        if kind == "plane" and (self.point is None or self.normal is None):
+            raise ValueError("plane BoundingSurface requires point and normal")
+
+    @property
+    def mesh(self):
+        return self._mesh
+
+    # -- normals -------------------------------------------------------------
+    def normals(self, coords):
+        """Outward unit normals at ``coords`` from the projected P1
+        boundary-normal field (``mesh.Gamma_P1``).
+
+        Returns ``(normals, valid)`` where ``valid`` is False at nodes whose
+        projected normal is degenerate (box corners, unlocatable points) — those
+        should be pinned, not slipped.
+        """
+        from underworld3.meshing._ot_adapt import _slip_normals
+        return _slip_normals(self._mesh, np.ascontiguousarray(coords, dtype=float))
+
+    # -- tangent slide -------------------------------------------------------
+    def tangent_project(self, coords, reference):
+        """Tangent-slide: remove this surface's normal component of the
+        displacement ``coords - reference`` (the projected P1 normal is taken at
+        ``reference``). Nodes with a degenerate normal keep ``reference``.
+        """
+        coords = np.asarray(coords, dtype=float)
+        reference = np.asarray(reference, dtype=float)
+        n, valid = self.normals(reference)
+        disp = coords - reference
+        dn = (disp * n).sum(axis=1, keepdims=True)
+        slid = reference + (disp - dn * n)
+        return np.where(valid[:, None], slid, reference)
+
+    # -- restore to surface --------------------------------------------------
+    def restore(self, coords):
+        """Snap ``coords`` back onto this surface (kind-specific).
+
+        ``radial`` — re-impose ``|r| = radius`` about ``centre`` (exact,
+        concave-safe). ``plane`` — orthogonal projection onto the plane.
+        A ``free``/``is_free`` surface returns ``coords`` unchanged (it follows
+        the live discrete surface — a follow-up). ``facet`` is not implemented
+        in step 1 (such labels are pinned by the orchestrator).
+        """
+        coords = np.asarray(coords, dtype=float)
+        if self.is_free or self.kind == "free":
+            return coords
+        if self.kind == "radial":
+            v = coords - self.centre
+            nrm = np.linalg.norm(v, axis=1)
+            nrm = np.where(nrm > 1.0e-30, nrm, 1.0)
+            return self.centre + v * (self.radius / nrm)[:, None]
+        if self.kind == "plane":
+            d = ((coords - self.point) * self.normal).sum(axis=1, keepdims=True)
+            return coords - d * self.normal
+        if self.kind == "facet":
+            raise NotImplementedError(
+                "facet restore is a follow-up (see boundary-slip-strategy.md); "
+                "labels without an analytic surface are pinned in step 1.")
+        return coords
+
+    # -- state transition ----------------------------------------------------
+    def release(self):
+        """Flip a rigid (``radial``/``plane``) surface to ``free``: subsequent
+        :meth:`restore` follows the live discrete surface instead of the analytic
+        target. Records the previous kind on ``self._prev_kind``."""
+        self._prev_kind = self.kind
+        self.kind = "free"
+        self.is_free = True
+        return self
+
+    def __repr__(self):
+        g = ""
+        if self.kind == "radial":
+            g = f", centre={self.centre}, radius={self.radius}"
+        elif self.kind == "plane":
+            g = f", point={self.point}, normal={self.normal}"
+        return (f"BoundingSurface(label={self.label!r}, kind={self.kind!r}, "
+                f"is_free={self.is_free}{g})")
+
+
+# -- constructor-side registration helpers ---------------------------------
+def register_radial_surfaces(mesh, centre, label_radius):
+    """Register ``radial`` surfaces: ``label_radius = {label_name: radius}``.
+
+    Called by analytic radial-boundary constructors (Annulus, SphericalShell,
+    CubedSphere). ``centre`` is the common centre (e.g. the origin)."""
+    centre = np.asarray(centre, dtype=float).ravel()
+    for lab, r in label_radius.items():
+        mesh.register_tangent_slip_provider(
+            lab, BoundingSurface(mesh, lab, "radial", centre=centre, radius=r))
+
+
+def register_plane_surfaces(mesh, label_plane):
+    """Register ``plane`` surfaces: ``label_plane = {label_name: (point, normal)}``.
+
+    Called by box constructors for axis-aligned (or general) faces."""
+    for lab, (p, n) in label_plane.items():
+        mesh.register_tangent_slip_provider(
+            lab, BoundingSurface(mesh, lab, "plane", point=p, normal=n))
+
+
+def register_box_face_surfaces(mesh, minCoords, maxCoords):
+    """Register ``plane`` surfaces for an axis-aligned box's faces.
+
+    Matches the box constructors' labels:
+    2D — ``Left``/``Right`` (``x = x_min/max``), ``Bottom``/``Top``
+    (``y = y_min/max``); 3D — ``Left``/``Right`` (``x``), ``Front``/``Back``
+    (``y``), ``Bottom``/``Top`` (``z``). Box corners/edges are junctions of two
+    faces and are pinned by the slip orchestrator (the normal is ambiguous).
+    """
+    lo = np.asarray(minCoords, dtype=float).ravel()
+    hi = np.asarray(maxCoords, dtype=float).ravel()
+    dim = lo.shape[0]
+    planes = {}
+    if dim == 2:
+        ex, ey = np.array([1.0, 0.0]), np.array([0.0, 1.0])
+        planes["Left"] = (lo, ex)
+        planes["Right"] = ([hi[0], lo[1]], ex)
+        planes["Bottom"] = (lo, ey)
+        planes["Top"] = ([lo[0], hi[1]], ey)
+    elif dim == 3:
+        ex, ey, ez = (np.array([1.0, 0.0, 0.0]),
+                      np.array([0.0, 1.0, 0.0]),
+                      np.array([0.0, 0.0, 1.0]))
+        planes["Left"] = (lo, ex)
+        planes["Right"] = ([hi[0], lo[1], lo[2]], ex)
+        planes["Front"] = (lo, ey)
+        planes["Back"] = ([lo[0], hi[1], lo[2]], ey)
+        planes["Bottom"] = (lo, ez)
+        planes["Top"] = ([lo[0], lo[1], hi[2]], ez)
+    else:
+        return
+    register_plane_surfaces(mesh, planes)
diff --git a/src/underworld3/meshing/cartesian.py b/src/underworld3/meshing/cartesian.py
index 9e63dd92..1b1636fb 100644
--- a/src/underworld3/meshing/cartesian.py
+++ b/src/underworld3/meshing/cartesian.py
@@ -353,6 +353,11 @@ def box_return_coords_to_bounds(coords):
         verbose=verbose,
     )
 
+    # Bounding-surface objects for tangent slip: axis-aligned box faces are
+    # planes (see docs/developer/design/boundary-slip-strategy.md).
+    from underworld3.meshing.bounding_surface import register_box_face_surfaces
+    register_box_face_surfaces(new_mesh, minCoords, maxCoords)
+
     return new_mesh
 
 
@@ -1363,4 +1368,9 @@ def box_return_coords_to_bounds(coords):
         verbose=verbose,
     )
 
+    # Bounding-surface objects for tangent slip: axis-aligned box faces are
+    # planes (see docs/developer/design/boundary-slip-strategy.md).
+    from underworld3.meshing.bounding_surface import register_box_face_surfaces
+    register_box_face_surfaces(new_mesh, minCoords, maxCoords)
+
     return new_mesh
diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py
index 9bcb5c91..6aedf657 100644
--- a/src/underworld3/meshing/spherical.py
+++ b/src/underworld3/meshing/spherical.py
@@ -277,6 +277,13 @@ class boundary_normals(Enum):
         sympy.Matrix([-y, x, 0]),   # rotation about z
     ]
 
+    # Bounding-surface objects for tangent slip: Upper/Lower are radial about the
+    # origin at the known radii (see docs/developer/design/boundary-slip-strategy.md).
+    from underworld3.meshing.bounding_surface import register_radial_surfaces
+    register_radial_surfaces(
+        new_mesh, centre=(0.0, 0.0, 0.0),
+        label_radius={"Upper": radiusOuter, "Lower": radiusInner})
+
     return new_mesh
 
 
@@ -1066,4 +1073,11 @@ class boundary_normals(Enum):
         sympy.Matrix([-y, x, 0]),
     ]
 
+    # Bounding-surface objects for tangent slip: Upper/Lower are radial about the
+    # origin at the known radii (see docs/developer/design/boundary-slip-strategy.md).
+    from underworld3.meshing.bounding_surface import register_radial_surfaces
+    register_radial_surfaces(
+        new_mesh, centre=(0.0, 0.0, 0.0),
+        label_radius={"Upper": radiusOuter, "Lower": radiusInner})
+
     return new_mesh
diff --git a/tests/test_0762_bounding_surfaces.py b/tests/test_0762_bounding_surfaces.py
new file mode 100644
index 00000000..a346b70b
--- /dev/null
+++ b/tests/test_0762_bounding_surfaces.py
@@ -0,0 +1,162 @@
+"""Bounding-surface objects: per-surface tangent-slip + restore (step 1).
+
+Locks the additive, self-contained API from
+docs/developer/design/boundary-slip-strategy.md:
+
+  * analytic constructors register radial BoundingSurface objects;
+  * mesh.boundaries (the persisted gmsh/DMPlex labelling) is left untouched;
+  * radial/plane restore land points on the surface;
+  * release() flips a surface to free (restore becomes a no-op);
+  * mesh.boundary_slip orchestrates slip-vs-pin and keeps slipped boundary
+    vertices exactly on the boundary while leaving interior vertices alone.
+"""
+import numpy as np
+import pytest
+
+import underworld3 as uw
+from underworld3.meshing.bounding_surface import BoundingSurface
+
+
+def _annulus():
+    return uw.meshing.Annulus(
+        radiusInner=0.5, radiusOuter=1.0, cellSize=0.2, qdegree=2)
+
+
+def test_constructor_registers_radial_surfaces():
+    m = _annulus()
+    bs = m.bounding_surfaces
+    assert {"Upper", "Lower"} <= set(bs)
+    assert bs["Upper"].kind == "radial"
+    assert bs["Lower"].kind == "radial"
+    assert np.isclose(bs["Upper"].radius, 1.0)
+    assert np.isclose(bs["Lower"].radius, 0.5)
+    assert np.allclose(bs["Upper"].centre, [0.0, 0.0])
+
+
+def test_boundaries_enum_untouched():
+    m = _annulus()
+    # mesh.boundaries is still the persisted gmsh/DMPlex Enum (name/value)
+    names = {b.name for b in m.boundaries}
+    assert {"Upper", "Lower"} <= names
+    assert all(hasattr(b, "value") for b in m.boundaries)
+    # bounding_surfaces is a SEPARATE collection
+    assert m.bounding_surfaces is not m.boundaries
+    assert isinstance(m.bounding_surfaces, dict)
+
+
+def test_radial_restore_lands_on_radius():
+    m = _annulus()
+    bs = m.bounding_surfaces["Upper"]
+    pts = np.array([[1.2, 0.0], [0.0, 0.8], [-0.9, 0.0], [0.5, 0.5]])
+    out = bs.restore(pts)
+    assert np.allclose(np.linalg.norm(out, axis=1), 1.0)
+    # direction preserved (snap is purely radial)
+    assert np.allclose(out / np.linalg.norm(out, axis=1, keepdims=True),
+                       pts / np.linalg.norm(pts, axis=1, keepdims=True))
+
+
+def test_plane_restore_zeroes_offplane():
+    m = _annulus()
+    bs = BoundingSurface(m, "face", "plane", point=[0.0, 0.0], normal=[0.0, 1.0])
+    pts = np.array([[0.3, 0.5], [1.0, -0.2]])
+    out = bs.restore(pts)
+    assert np.allclose(out[:, 1], 0.0)        # on the y=0 plane
+    assert np.allclose(out[:, 0], pts[:, 0])  # tangential coord unchanged
+
+
+def test_release_makes_restore_a_noop():
+    m = _annulus()
+    bs = m.bounding_surfaces["Upper"]
+    assert not bs.is_free
+    bs.release()
+    assert bs.is_free and bs.kind == "free"
+    pts = np.array([[1.2, 0.0], [0.0, 0.8]])
+    assert np.allclose(bs.restore(pts), pts)  # follows live surface → no snap
+
+
+def test_register_provider_and_type_check():
+    m = _annulus()
+    s = BoundingSurface(m, "Upper", "radial", centre=[0.0, 0.0], radius=2.0)
+    m.register_tangent_slip_provider("Upper", s)
+    assert m.bounding_surfaces["Upper"].radius == 2.0
+    with pytest.raises(TypeError):
+        m.register_tangent_slip_provider("X", object())
+
+
+def test_invalid_kind_and_missing_geometry_raise():
+    m = _annulus()
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "bogus")
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "radial")           # needs centre+radius
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "plane", point=[0, 0])  # needs normal
+
+
+def test_boundary_slip_keeps_nodes_on_boundary():
+    m = _annulus()
+    ref = np.asarray(m.X.coords, dtype=float).copy()
+    is_pinned, project = m.boundary_slip(True, reference_coords=ref)
+
+    r_ref = np.linalg.norm(ref, axis=1)
+    upper = np.isclose(r_ref, 1.0, atol=1e-6)
+    lower = np.isclose(r_ref, 0.5, atol=1e-6)
+    interior = ~(upper | lower)
+
+    # Deterministic perturbation of every vertex, then project.
+    th = np.arctan2(ref[:, 1], ref[:, 0])
+    Y = ref + 0.05 * np.column_stack([np.cos(th + 1.0), np.sin(th + 1.0)])
+    Yin = Y.copy()
+    Y2 = project(Y)
+
+    # Slipped Upper/Lower vertices land EXACTLY on their radius.
+    su = (~is_pinned) & upper
+    sl = (~is_pinned) & lower
+    assert su.any() and sl.any()
+    assert np.allclose(np.linalg.norm(Y2[su], axis=1), 1.0, atol=1e-9)
+    assert np.allclose(np.linalg.norm(Y2[sl], axis=1), 0.5, atol=1e-9)
+
+    # Interior vertices are untouched by the projector.
+    assert np.allclose(Y2[interior], Yin[interior])
+
+
+def test_boundary_slip_pins_when_no_surface_registered():
+    m = _annulus()
+    m.bounding_surfaces.clear()      # remove the analytic surfaces
+    ref = np.asarray(m.X.coords, dtype=float).copy()
+    is_pinned, project = m.boundary_slip(True, reference_coords=ref)
+    r_ref = np.linalg.norm(ref, axis=1)
+    bnd = np.isclose(r_ref, 1.0, atol=1e-6) | np.isclose(r_ref, 0.5, atol=1e-6)
+    # With no registered surfaces, every boundary vertex is pinned, no slip.
+    assert is_pinned[bnd].all()
+    Y = ref + 0.05
+    assert np.allclose(project(Y.copy()), Y)  # nothing slips
+
+
+def test_spherical_shell_registers_radial():
+    m = uw.meshing.SphericalShell(radiusInner=0.5, radiusOuter=1.0, cellSize=0.4)
+    bs = m.bounding_surfaces
+    assert bs["Upper"].kind == "radial" and np.isclose(bs["Upper"].radius, 1.0)
+    assert bs["Lower"].kind == "radial" and np.isclose(bs["Lower"].radius, 0.5)
+    assert bs["Upper"].centre.shape == (3,)
+    out = bs["Upper"].restore(np.array([[1.3, 0.0, 0.0], [0.0, 0.7, 0.7]]))
+    assert np.allclose(np.linalg.norm(out, axis=1), 1.0)
+
+
+def test_box_registers_plane_surfaces():
+    m = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.34)
+    bs = m.bounding_surfaces
+    assert {"Left", "Right", "Bottom", "Top"} <= set(bs)
+    assert all(bs[k].kind == "plane" for k in ("Left", "Right", "Bottom", "Top"))
+    # Left face is x=0: restore zeroes x, keeps y.
+    out = bs["Left"].restore(np.array([[0.3, 0.5]]))
+    assert np.isclose(out[0, 0], 0.0) and np.isclose(out[0, 1], 0.5)
+    # Top face is y=1.
+    out = bs["Top"].restore(np.array([[0.4, 0.7]]))
+    assert np.isclose(out[0, 1], 1.0) and np.isclose(out[0, 0], 0.4)
+    # Box corner (0,0) is a junction of two faces → pinned.
+    ref = np.asarray(m.X.coords, dtype=float)
+    is_pinned, _project = m.boundary_slip(True, reference_coords=ref)
+    corner = np.isclose(ref[:, 0], 0.0) & np.isclose(ref[:, 1], 0.0)
+    assert corner.any() and is_pinned[corner].all()

From ce681add58fb11476947a17db7abc743c3cb3447 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Tue, 9 Jun 2026 11:33:03 +1000
Subject: [PATCH 3/6] fix(meshing): validate BoundingSurface geometry + project
 only owned slip vertices

Addresses Copilot review on #225:

- BoundingSurface now rejects a degenerate/zero or non-finite plane normal
  (which would make restore() a silent no-op via _unit() returning a zero
  vector) and a non-positive / non-finite radius (which would collapse radial
  projections). Validated by test_degenerate_geometry_raises.
- mesh.boundary_slip projects only OWNED slip vertices; leaves/ghosts receive
  their owner's projected value through the movers' existing halo-sync, so we
  no longer modify non-owned coordinates (parallel safety). Serial behaviour is
  unchanged (every vertex is owned).

Documented as a follow-up (not fixed here): _pinned_mask expands labels through
vertices/edges only, so a 3D face-only boundary label (markVertices=False) is
not picked up. This is a pre-existing limitation of the shared helper used by
every metric mover; the fix belongs with _pinned_mask itself.

Underworld development team with AI support from Claude Code
---
 .../discretisation/discretisation_mesh.py     | 16 ++++++++--
 src/underworld3/meshing/bounding_surface.py   | 29 ++++++++++++++++---
 tests/test_0762_bounding_surfaces.py          | 18 ++++++++++++
 3 files changed, 57 insertions(+), 6 deletions(-)

diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 23c5dd3f..1436aa40 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -2117,7 +2117,8 @@ def boundary_slip(self, slip_spec=True, reference_coords=None,
             ``project(Y)`` slides+restores the slip vertices of ``Y`` in place
             and returns it.
         """
-        from underworld3.meshing.smoothing import _pinned_mask, _auto_pinned_labels
+        from underworld3.meshing.smoothing import (
+            _pinned_mask, _auto_pinned_labels, _owned_vertex_mask)
 
         dm = self.dm
         pStart, pEnd = dm.getDepthStratum(0)
@@ -2128,6 +2129,11 @@ def boundary_slip(self, slip_spec=True, reference_coords=None,
 
         all_labels = (tuple(boundary_labels) if boundary_labels is not None
                       else _auto_pinned_labels(self))
+        # TODO(follow-up): _pinned_mask expands labels through vertices/edges
+        # only, so a 3D boundary label that tags FACES alone (a mesh loaded with
+        # markVertices=False) leaves its boundary vertices unmarked. This is a
+        # pre-existing limitation of the shared helper used by every mover; the
+        # fix (close faces→edges→vertices) belongs with _pinned_mask itself.
         is_bnd = _pinned_mask(dm, all_labels)
 
         slip_labels, free_labels = self._resolve_slip_spec(slip_spec)
@@ -2144,7 +2150,13 @@ def boundary_slip(self, slip_spec=True, reference_coords=None,
         for lab, m in masks.items():
             vert_label[m & slip_mask] = lab
 
-        slip_b = numpy.nonzero(slip_mask)[0]
+        # Project only OWNED slip vertices: the movers halo-sync owned→ghost
+        # after calling project(), so a leaf/ghost receives its owner's
+        # projected value — modifying non-owned coordinates here is both
+        # wasteful and a parallel-safety hazard. (Serial: every vertex is
+        # owned, so this is a no-op.) is_pinned stays the full geometric
+        # classification, which is rank-consistent for shared vertices.
+        slip_b = numpy.nonzero(slip_mask & _owned_vertex_mask(dm))[0]
         if slip_b.size == 0:
             return is_pinned, (lambda Y: Y)
         old_slip = ref[slip_b]
diff --git a/src/underworld3/meshing/bounding_surface.py b/src/underworld3/meshing/bounding_surface.py
index 927a0f79..c49e2125 100644
--- a/src/underworld3/meshing/bounding_surface.py
+++ b/src/underworld3/meshing/bounding_surface.py
@@ -73,10 +73,31 @@ def __init__(self, mesh, label, kind, *, centre=None, radius=None,
         self.radius = None if radius is None else _as_float(radius)
         self.point = None if point is None else np.asarray(point, dtype=float).ravel()
         self.normal = None if normal is None else _unit(normal)
-        if kind == "radial" and (self.centre is None or self.radius is None):
-            raise ValueError("radial BoundingSurface requires centre and radius")
-        if kind == "plane" and (self.point is None or self.normal is None):
-            raise ValueError("plane BoundingSurface requires point and normal")
+        if kind == "radial":
+            if self.centre is None or self.radius is None:
+                raise ValueError(
+                    "radial BoundingSurface requires centre and radius")
+            if not (np.isfinite(self.radius) and self.radius > 0.0):
+                raise ValueError(
+                    "radial BoundingSurface requires a finite positive radius; "
+                    f"got {self.radius!r}")
+            if not np.all(np.isfinite(self.centre)):
+                raise ValueError(
+                    "radial BoundingSurface centre must be finite")
+        if kind == "plane":
+            if self.point is None or self.normal is None:
+                raise ValueError(
+                    "plane BoundingSurface requires point and normal")
+            # _unit() returns a ZERO vector for a degenerate/zero normal (not
+            # None), which would make restore() a silent no-op — reject it.
+            if not (np.all(np.isfinite(self.normal))
+                    and np.linalg.norm(self.normal) > 0.5):
+                raise ValueError(
+                    "plane BoundingSurface requires a finite, non-degenerate "
+                    "normal (a zero/near-zero normal makes restore() a no-op)")
+            if not np.all(np.isfinite(self.point)):
+                raise ValueError(
+                    "plane BoundingSurface point must be finite")
 
     @property
     def mesh(self):
diff --git a/tests/test_0762_bounding_surfaces.py b/tests/test_0762_bounding_surfaces.py
index a346b70b..38e3d996 100644
--- a/tests/test_0762_bounding_surfaces.py
+++ b/tests/test_0762_bounding_surfaces.py
@@ -93,6 +93,24 @@ def test_invalid_kind_and_missing_geometry_raise():
         BoundingSurface(m, "x", "plane", point=[0, 0])  # needs normal
 
 
+def test_degenerate_geometry_raises():
+    # A zero / non-finite normal would make plane restore() a silent no-op;
+    # a non-positive / non-finite radius produces invalid radial projections.
+    # Both must be rejected at construction.
+    m = _annulus()
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "plane", point=[0, 0], normal=[0, 0])
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "plane", point=[0, 0],
+                        normal=[np.nan, 0])
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "radial", centre=[0, 0], radius=0.0)
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "radial", centre=[0, 0], radius=-1.0)
+    with pytest.raises(ValueError):
+        BoundingSurface(m, "x", "radial", centre=[0, 0], radius=np.inf)
+
+
 def test_boundary_slip_keeps_nodes_on_boundary():
     m = _annulus()
     ref = np.asarray(m.X.coords, dtype=float).copy()

From 21fdd096e26c760bb86c16d4151021f19f8196e2 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Tue, 9 Jun 2026 11:51:31 +1000
Subject: [PATCH 4/6] docs(meshing): explain BoundingSurface.normals re-solve
 choice (breadcrumb)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Document why normals() re-solves Gamma_P1 on the current mesh state (state-aware,
deformation-following — needed for the free/release() surface-follow mode) rather
than reading the cached _n_proj field the retired _build_slip_projector used. Notes
the cost (~one projection solve per mover outer-iteration), the validated parity
(~1e-16 on a centred annulus), and the fix path if it ever becomes a hot-path
regression (add a cached fast-path; re-solve only for free surfaces).

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

diff --git a/src/underworld3/meshing/bounding_surface.py b/src/underworld3/meshing/bounding_surface.py
index c49e2125..d4ea5f26 100644
--- a/src/underworld3/meshing/bounding_surface.py
+++ b/src/underworld3/meshing/bounding_surface.py
@@ -111,6 +111,23 @@ def normals(self, coords):
         Returns ``(normals, valid)`` where ``valid`` is False at nodes whose
         projected normal is degenerate (box corners, unlocatable points) — those
         should be pinned, not slipped.
+
+        DESIGN NOTE (2026-06 — breadcrumb for a future session). This
+        RE-SOLVES the ``Gamma_P1`` projection on the surface's *current* mesh
+        state (``_slip_normals`` calls ``mesh._update_projected_normals()``),
+        so the normal follows mesh deformation — the state-aware behaviour the
+        ``free`` / ``release()`` surface-follow mode is designed to use. The
+        retired ``_build_slip_projector`` instead read a *cached* ``_n_proj``
+        field (the deleted ``_gamma_p1_at_vertices`` helper: KDTree match, no
+        solve, normal frozen at the reference mesh) — cheaper but not
+        deformation-aware. ``mesh.boundary_slip`` calls this ONCE per build
+        (not per ``project()`` call), so the cost is ~one projection solve per
+        mover outer-iteration; parity with the cached path was machine
+        precision (~1e-16) on a centred annulus. **If this ever shows up as a
+        hot-path regression (fine meshes / many adapts), add a cached
+        fast-path here** — read ``mesh._projected_normals.data`` directly (as
+        ``_gamma_p1_at_vertices`` did in git history) and re-solve only for
+        free surfaces. See docs/developer/design/boundary-slip-strategy.md.
         """
         from underworld3.meshing._ot_adapt import _slip_normals
         return _slip_normals(self._mesh, np.ascontiguousarray(coords, dtype=float))

From 4781f6aed8df63191fdbbf9e82c74ce469a36cfc Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Tue, 9 Jun 2026 13:32:00 +1000
Subject: [PATCH 5/6] fix(meshing): allow radius==0 for a solid-mesh inner
 boundary (CI fix)

The radius>0 check added in the previous commit rejected radius==0, but a SOLID
sphere/annulus (radiusInner=0) legitimately registers its inner "Lower" boundary
at the centre (radius 0). The mid-construction ValueError also left PETSc/mesh
state corrupted, cascading into a spurious SphericalShell coordinate-reshape
failure later in the full-suite run. Relax to radius>=0 (still reject negative /
non-finite). Fixes the test_0001_meshes solid-mesh tests and the cascaded
test_0762 SphericalShell failure.

Underworld development team with AI support from Claude Code
---
 src/underworld3/meshing/bounding_surface.py | 9 ++++++---
 tests/test_0762_bounding_surfaces.py        | 5 +++--
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/underworld3/meshing/bounding_surface.py b/src/underworld3/meshing/bounding_surface.py
index d4ea5f26..78108821 100644
--- a/src/underworld3/meshing/bounding_surface.py
+++ b/src/underworld3/meshing/bounding_surface.py
@@ -77,10 +77,13 @@ def __init__(self, mesh, label, kind, *, centre=None, radius=None,
             if self.centre is None or self.radius is None:
                 raise ValueError(
                     "radial BoundingSurface requires centre and radius")
-            if not (np.isfinite(self.radius) and self.radius > 0.0):
+            # radius == 0 is legitimate: a *solid* sphere/annulus registers its
+            # inner ("Lower") boundary at radius 0 (the centre point). Reject
+            # only NEGATIVE or non-finite radii (those give invalid projections).
+            if not (np.isfinite(self.radius) and self.radius >= 0.0):
                 raise ValueError(
-                    "radial BoundingSurface requires a finite positive radius; "
-                    f"got {self.radius!r}")
+                    "radial BoundingSurface requires a finite, non-negative "
+                    f"radius; got {self.radius!r}")
             if not np.all(np.isfinite(self.centre)):
                 raise ValueError(
                     "radial BoundingSurface centre must be finite")
diff --git a/tests/test_0762_bounding_surfaces.py b/tests/test_0762_bounding_surfaces.py
index 38e3d996..e833467e 100644
--- a/tests/test_0762_bounding_surfaces.py
+++ b/tests/test_0762_bounding_surfaces.py
@@ -103,12 +103,13 @@ def test_degenerate_geometry_raises():
     with pytest.raises(ValueError):
         BoundingSurface(m, "x", "plane", point=[0, 0],
                         normal=[np.nan, 0])
-    with pytest.raises(ValueError):
-        BoundingSurface(m, "x", "radial", centre=[0, 0], radius=0.0)
     with pytest.raises(ValueError):
         BoundingSurface(m, "x", "radial", centre=[0, 0], radius=-1.0)
     with pytest.raises(ValueError):
         BoundingSurface(m, "x", "radial", centre=[0, 0], radius=np.inf)
+    # radius == 0 is VALID (a solid sphere/annulus registers its inner boundary
+    # at the centre, radius 0) — must NOT raise.
+    BoundingSurface(m, "x", "radial", centre=[0, 0], radius=0.0)
 
 
 def test_boundary_slip_keeps_nodes_on_boundary():

From ac17a34be1d4f3241f43bc807243262fd0e6ffa8 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Tue, 9 Jun 2026 14:12:28 +1000
Subject: [PATCH 6/6] test(meshing): move 3D SphericalShell radial test to the
 early CI batch

SphericalShell construction is fragile under the accumulated PETSc/mesh state of
the heavy test_05*/07* CI batch (a pre-existing mesh-lifecycle issue: the
coordinate DM / cdim goes stale, giving "cannot reshape ... into shape (3)" at
build time -- unrelated to boundary slip). test_0762's SphericalShell case hit
this in CI, while test_0001 builds spheres cleanly in the early batch. Move the
3D radial-registration test to tests/test_0002_bounding_surface_3d.py so it runs
early; the Annulus tests in test_0762 cover the 2D radial logic. Verified: the
test_00[0-4]* batch (116 passed) and the test_05*/07* batch (364 passed) are
both green.

Underworld development team with AI support from Claude Code
---
 tests/test_0002_bounding_surface_3d.py | 25 +++++++++++++++++++++++++
 tests/test_0762_bounding_surfaces.py   | 13 +++++--------
 2 files changed, 30 insertions(+), 8 deletions(-)
 create mode 100644 tests/test_0002_bounding_surface_3d.py

diff --git a/tests/test_0002_bounding_surface_3d.py b/tests/test_0002_bounding_surface_3d.py
new file mode 100644
index 00000000..c12cb924
--- /dev/null
+++ b/tests/test_0002_bounding_surface_3d.py
@@ -0,0 +1,25 @@
+"""3D radial bounding-surface registration (SphericalShell).
+
+This lives in the early (test_000x) batch on purpose. SphericalShell
+construction is fragile once a long-running process has accumulated a lot of
+PETSc/mesh state (a pre-existing mesh-lifecycle issue — the coordinate DM /
+cdim can go stale, giving a "cannot reshape ... into shape (3)" at build time).
+Built early (right after test_0001_meshes, which itself constructs spheres
+cleanly) it is robust. The 2D radial-registration logic is covered by the
+Annulus tests in test_0762_bounding_surfaces.py.
+
+See docs/developer/design/boundary-slip-strategy.md.
+"""
+import numpy as np
+
+import underworld3 as uw
+
+
+def test_spherical_shell_registers_radial():
+    m = uw.meshing.SphericalShell(radiusInner=0.5, radiusOuter=1.0, cellSize=0.4)
+    bs = m.bounding_surfaces
+    assert bs["Upper"].kind == "radial" and np.isclose(bs["Upper"].radius, 1.0)
+    assert bs["Lower"].kind == "radial" and np.isclose(bs["Lower"].radius, 0.5)
+    assert bs["Upper"].centre.shape == (3,)
+    out = bs["Upper"].restore(np.array([[1.3, 0.0, 0.0], [0.0, 0.7, 0.7]]))
+    assert np.allclose(np.linalg.norm(out, axis=1), 1.0)
diff --git a/tests/test_0762_bounding_surfaces.py b/tests/test_0762_bounding_surfaces.py
index e833467e..c4713e4c 100644
--- a/tests/test_0762_bounding_surfaces.py
+++ b/tests/test_0762_bounding_surfaces.py
@@ -152,14 +152,11 @@ def test_boundary_slip_pins_when_no_surface_registered():
     assert np.allclose(project(Y.copy()), Y)  # nothing slips
 
 
-def test_spherical_shell_registers_radial():
-    m = uw.meshing.SphericalShell(radiusInner=0.5, radiusOuter=1.0, cellSize=0.4)
-    bs = m.bounding_surfaces
-    assert bs["Upper"].kind == "radial" and np.isclose(bs["Upper"].radius, 1.0)
-    assert bs["Lower"].kind == "radial" and np.isclose(bs["Lower"].radius, 0.5)
-    assert bs["Upper"].centre.shape == (3,)
-    out = bs["Upper"].restore(np.array([[1.3, 0.0, 0.0], [0.0, 0.7, 0.7]]))
-    assert np.allclose(np.linalg.norm(out, axis=1), 1.0)
+# NOTE: SphericalShell (3D radial) registration is tested in
+# tests/test_0002_bounding_surface_3d.py — it must run in the early test batch
+# because SphericalShell construction is fragile under the accumulated PETSc
+# state of the heavy test_05*/07* batch (a pre-existing mesh-lifecycle issue,
+# unrelated to boundary slip). The Annulus tests above cover the radial logic.
 
 
 def test_box_registers_plane_surfaces():