Use the nested hierarchy for MG transfers on adapt children (topological prolongation)

[lmoresi]

Design note: docs/developer/design/nested-vs-geometric-mg-transfers.md (commit e984a3b).

The inconsistency

mesh.adapt() maintains an exact refinement hierarchy — that is what newest-vertex bisection with conforming closure is for, and what the escape/halo cost buys. The MG transfer then discards it and re-derives an approximation geometrically:

tri  = Delaunay(coarse_coords)          # re-triangulate the coarse DOF cloud
simp = tri.find_simplex(fine_coords)    # locate each fine DOF in one simplex

We pay for exact hierarchy maintenance and then do not use it.

This is the root cause of #424, not a separate concern: the barycentric builder has only local support, so moving fine coordinates (relax(), snapping, free-surface deformation) can leave a coarse DOF with no fine image -> zero column -> singular PᵀAP.

The nested transfer is trivial

Every fine vertex is either an inherited coarse vertex (weight 1 on itself) or the midpoint of a coarse edge (weights ½, ½ on its endpoints). It is:

• exact for the un-deformed nested mesh (it IS the FE embedding);
• structurally full rank — every coarse DOF appears with weight 1 in its own inherited fine DOF, so 3D: shape relaxation can break the custom-P FMG build (zero-column transfer) — silent GAMG fallback #424 is impossible by construction;
• immune to node motion (topological, not geometric);
• 1–2 non-zeros per fine row, no Delaunay, no point location.

Both engines already record it: edge2mid in NVBMesh and TaggedBisectionMesh, and the native nvb_transform.c is a real DMPlexTransform with the same split-edge bookkeeping. Separately, petsc_dm_set_regular_refinement — our wrapper for the flag PETSc checks to take its exact nested interpolation path — appears exactly once in the repo: its own definition. Nothing calls it.

Scope: keep BOTH paths

The geometric builder stays, for the non-nested pairs it was asked for (moved base vs adapted child, unrelated meshes, and the planned non-hierarchical variant with newly added seed points). The change is to stop using it where a topological relation exists.

Selection: topological when the parent/child map is present, geometric otherwise. The RBF retry from #424 then becomes a fallback that should essentially never fire on adapt children.

The hard part — capture ordering

The map must be in DM point numbering, and the bridge from engine vertex ids to DM points is coordinate matching, which is exact only while midpoints are still exact float averages of their parents. Snapping and relaxation destroy that:

• native path: _nvbx.refine() -> DM, then snap_level_boundaries(). There is a clean window.
• cell-list path: the snap is applied to nvb.coords before to_dm(), so by the time a DM exists the window has closed.

Preferred fix: have to_dm() return its engine-id -> DM-point map directly, removing the dependence on coordinate matching altogether.

Proposed shape

• child._adapt_prolongation[k] = (parents, weights), parents shape (n_fine, 2) in DM vertex numbering; inherited vertex is (c, c) with weights (1, 0).
• custom_mg gains a "topological" builder consuming it, expanded from vertices to DOFs for the field basis.

To measure, not assume

With snap-every-generation a midpoint is moved onto a curved boundary, so it is no longer at ½(a+b) physically, and the ½,½ transfer is then not the FE-exact interpolation there — the geometric builder is. MG needs full rank and good representation of smooth functions, not FE-exactness, so ½,½ should be fine (it is the classical geometric-MG choice). But on strongly curved or relaxed meshes geometric may be more accurate. Accuracy vs robustness — measure on the annulus and spherical-shell cases before settling the default.

[lmoresi]

Measured what the geometric builder actually computes, prompted by "does it assume bisection of edges?" — no, it assumes nothing about refinement, and that turns out to matter more than expected.

barycentric_prolongation triangulates the coarse DOF point cloud from scratch (Delaunay(coarse_coords)) and never consults the coarse mesh cells. Its docstring claimed each fine DOF is interpolated by "the coarse element that contains it" — FE-exact. That only holds when the Delaunay triangulation coincides with the mesh. Measured agreement:

case	mesh cells	Delaunay simplices	agreement
2D uniform base, P1	168	168	100%
2D adapt child, P1	238	238	90.8%
3D uniform base, P1	800	812	58.8%
3D adapt child, P1	4685	5201	17.1%
3D adapt child, P2	4685	39597	—

So it is linear interpolation over an ad-hoc triangulation, not the FE embedding. Linear reproduction is why MG still converges well; but the coarse P2 space is not represented exactly, and the 58.8% on a uniform base reflects PETSc 1:8 tetrahedral splits not being Delaunay (same fixed-diagonal property as the 148 degree dihedral).

Predicted risk, not yet measured: Delaunay simplices do not respect mesh topology. On a non-convex domain or one with an internal interface (a fault), a simplex can bridge the discontinuity and smear the coarse correction across it. Worth a direct test given the fault work.

This strengthens the case for stage 2. A parent-cell transfer is the genuine FE embedding, exact at any degree, and cannot bridge across mesh features because it never leaves the parent cell. That is a correctness argument, not just an efficiency one.

Docstring corrected in the linked commit.

[lmoresi]

Assessment: is a local linear-exact interpolant the right transfer for the non-hierarchical direction?

Context: PR #430 adds a sparse, linear-exact, topology-free local interpolant — polyharmonic r² log r with an affine tail, solved per target point on its nnn-nearest-neighbour stencil (utilities/rbf_stencil.py, reachable as KDTree.rbf_interpolator_local(..., order=1)). This issue's scope note anticipates a "planned non-hierarchical variant with newly added seed points", where points are placed freely and the mesh re-triangulated, so no parent/child relation exists. This is an assessment of that fit, not an implementation.

It fits the requirements, on paper

For an MG prolongation the transfer must be sparse, must reproduce linear fields (that is what makes the coarse-grid correction represent smooth error), and — for the re-triangulation case — must not need a topological relation between the two point sets. The local scheme is all three:

• sparse: exactly nnn non-zeros per row, measured, versus nnz/row == n_coarse for the current rbf_prolongation
• linear-exact: |Σw − 1| ≤ 4.4e-16, |Σ w_j x_j − x*| ≤ 5e-15, measured
• topology-free: it needs a kd-tree query and nothing else — no Delaunay, no find_simplex, no parent map

Compared with the two builders already in custom_mg.py: it is as sparse as barycentric_prolongation but needs no triangulation and no out-of-hull special case, and it is as linear-exact as rbf_prolongation without being dense.

The caveat that matters — it does NOT fix #424

This is the part worth stating clearly, because "sparse and linear-exact" reads as if it did.

#424 is a zero-column failure: a coarse DOF with no fine image, giving a singular PᵀAP. A row-wise kNN builder guarantees nnn non-zeros in every row — it guarantees nothing about columns. A coarse DOF that is not among the nnn nearest neighbours of any fine DOF still gets a zero column, and the Galerkin operator is still singular.

It makes the failure rarer than the barycentric builder does (a 2(dim+1) stencil reaches further than "the containing simplex", so isolated coarse points are more likely to be picked up by some fine row), but rarer is not impossible, and on a graded or strongly relaxed mesh the coarse DOFs most at risk of orphaning are exactly the ones in sparse regions where the fine points are far away. Any adoption here needs the same explicit zero-column check and repair that the current builders need — the geometric character of the interpolant is not what fixes it.

The topological transfer proposed in this issue is structurally different in exactly this respect: every coarse DOF appears with weight 1 in its own inherited fine DOF, so full column rank is guaranteed by construction rather than hoped for. That remains the right answer wherever a parent/child map exists, and this assessment does not compete with it.

Also worth measuring before adopting

• Cost. ~4.5x the inverse-distance builder per target point (measured on proxy refresh in Local RBF interpolator with a linear-reproduction guarantee #430: 24 ms vs 5 ms for 4.9k targets in 2D, 112 ms vs 27 ms for 12.7k in 3D). A prolongation is built once per level rather than once per solve, so this is likely irrelevant here — but it is a dense solve per row and should be confirmed, not assumed, at production coarse-grid sizes.
• Signed weights. The weights are not positive, so the prolongation is not a partition of unity in the "convex combination" sense (the row sums are still exactly 1). Classical geometric MG tolerates this — ½, ½ edge-midpoint interpolation is also not positivity-preserving in any stronger sense — but the interaction with the BC-per-level reduction should be checked rather than assumed.
• Degenerate stencils. Handled (rank pre-screen, reproduction check, inverse-distance fallback, single counted warning, never NaN), but on a coarse cloud the collinear/coplanar case is more likely than on a particle cloud, so the fallback fraction should be logged rather than left silent.

Summary

Good candidate for the non-hierarchical / re-triangulation direction specifically, where its topology-freedom is the whole point and the alternatives need a triangulation. Not a substitute for the nested topological transfer where a parent/child map exists, and not a fix for #424 — that needs a column check either way.

Not implemented; flagged here so the option is on record with its limits attached.

[lmoresi]

Stage 2 blocker identified precisely, and the route through it. Session prompt written (~/.claude/plans/stage2-any-degree-nested-mg-transfer.md).

What Stage 2 needs

For each fine DOF, its parent coarse CELL — then the weights are the coarse basis evaluated inside that cell. No point location, so no Delaunay, no spatial index, and structurally no orphaned coarse DOF.

Routes ruled out

1. Inferring the parent cell from the recorded vertex prolongation (intersecting incident-cell sets over the coarse vertices each fine vertex references). Some fine cells reference 5 coarse vertices and a tet has 4, so there is no candidate parent. Value-correctness of the vertex transfer does not imply unique cell attribution: a vertex on a shared face can be expressed from either adjacent cell, since P1 is continuous there. Declines on 3 of 6 3D passes.
2. DMPlexLocatePoints — petsc4py does not expose locatePoints on DMPlex (no *ocate* methods at all).

The route that works

nvb_transform.c already uses the required call, for tagged-state propagation (~lines 1501–1527):

DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_TETRAHEDRON,
                              DM_POLYTOPE_TETRAHEDRON, c, r, &pNew);

mapping coarse cell c, replica r, to its fine child, including the distinction between a bisected cell (2 children) and an identity cell (1).

The non-obvious part: UWNVBRefine is a drain loop — one call applies several transforms in sequence (while (nC > 0), DMPlexTransformApply ~line 1461). So one refine "pass" is not one bisection level, which is exactly why cascade vertices sit at quarter points of coarse edges — and the parent map must be composed across drain iterations, not read from a single transform.

Any-degree transfer without hand-coding bases

For a Lagrange element the basis is dual to its nodal points, so per parent cell: build the monomial matrix over the cell DOF coordinates, factor once, and solve for each fine DOF inside it. Exact for P1/P2/P3 alike, no per-degree formulae.

The bar

Do not assume the nested transfer wins. Measured here: Stokes at eta contrast 1e6 runs 3 iterations flat on the geometric builder (GAMG needs 38–53); and on 3D Poisson P1 the nested transfer was consistently one iteration worse than geometric (5,6,5 vs 4,5,4), likely because it is near-injection (1.0–1.4 nnz/row) where barycentric spreads over ~4. Judge Stage 2 on iterations and wall-clock — its coarse operators are ~3x sparser — and report a loss as a result rather than burying it.

[lmoresi]

Stage 3 added — the local RBF makes the geometric path viable without Delaunay

PR #430 (feature/linear-rbf, open, mergeable) lands
underworld3.utilities.rbf_stencil.linear_exact_weights — a sparse,
linear-exact, topology-free local interpolant. Pure NumPy, no PETSc.
Polyharmonic r^2 log r + affine tail, small saddle system per target.
Measured |sum(w) - 1| <= 4.4e-16, |sum(w_j x_j) - x*| <= 5e-15, nnn
non-zeros per row.

That closes a real gap: UW3 had sparse-but-constant-only (Shepard IDW) and
linear-exact-but-dense (custom_mg.rbf_prolongation, nnz/row == n_coarse),
and nothing that was both.

Stage 3 = replace barycentric_prolongation with it. That removes the
scipy.spatial.Delaunay dependency from the path Stokes actually uses (velocity
is P2, so every production Stokes solve on an adapt child goes through it), and
removes the topology-blind Delaunay simplices that can bridge an internal
interface such as a fault.

Traps recorded so they are not rediscovered

• A kNN builder does NOT fix 3D: shape relaxation can break the custom-P FMG build (zero-column transfer) — silent GAMG fallback #424. Row-wise kNN guarantees nnn non-zeros
  per row and nothing about columns; a coarse DOF can still have no fine
  image -> zero column -> singular PtAP. The orphan repair must stay and needs
  its own test.
• nnn >= dim+2 is enforced; at dim+1 the affine block goes square, the RBF
  block inert, and it collapses to bare barycentric interpolation on the
  neighbour simplex. Use 2*(dim+1).
• Degenerate stencils return zeros + a mask (not NaN) — but a zero row destroys
  partition of unity, so a fallback is mandatory for a prolongation.
• Rank-local, no MPI inside; hook the gather into _build_crosspart_transfer.
• Use .coords_nd, never .coords (Swarm proxy refresh is broken whenever reference quantities are set (dimensional coords into a non-dimensional kd-tree) #426).

The hypothesis worth testing

Stencil width appears to drive MG convergence more than exactness does. On
3D Poisson P1 the exact nested transfer (1.0-1.4 nnz/row, near-injection) was
consistently one iteration worse than barycentric (~4 nnz/row): 5,6,5 vs
4,5,4. The local RBF at nnn = 2*(dim+1) gives 6 (2D) / 8 (3D) per row —
wider than barycentric — so it may beat both. Treat nnn as a tunable and
sweep it; wider damps better but densifies the coarse operators, so there should
be an optimum. Unmeasured.

Bar: barycentric gives 3 iterations flat on Stokes at eta contrast 1e6 (GAMG
38-53). Match or beat, and report wall-clock as well as counts.

Ordering

Stages 2 and 3 are complementary, not alternatives. Where a parent/child map
exists the topological transfer is better — exact at any degree, cannot bridge
mesh features, structurally free of zero columns. Stage 3 covers everything
else and is the enabling piece for the non-hierarchical re-triangulation
direction.

Plan: ~/.claude/plans/mg-transfer-stages-2-and-3.md