Anisotropic mesh-adaptation movers + mesh-owned boundary tangent-slip + MPI robustness fixes

docs/advanced/mesh-adaptation.md

@@ -379,8 +379,8 @@ For the mathematically inclined, see the [Developer Design Document](../develope
 
 When you want to concentrate resolution on an evolving feature
 **every timestep** without re-meshing — keeping the topology and
-all field data intact — use the anisotropic metric mover instead
-of `mesh.adapt`:
+all field data intact — use `smooth_mesh_interior` (the node-moving
+mover) instead of `mesh.adapt`:
 
 ```python
 import underworld3 as uw
@@ -396,10 +396,41 @@ from underworld3.meshing import (
 rho = metric_density_from_gradient(mesh, T, amp=8.0)
 
 # Move the nodes to that metric (topology / DOFs / variables
-# all preserved — no transfer needed).
-smooth_mesh_interior(
-    mesh, metric=rho, method="anisotropic",
-    method_kwargs=dict(aniso_cap=2.0, relax=0.2, n_outer=12))
+# all preserved — no transfer needed). method="mmpde" is the
+# DEFAULT and may be omitted; shown here for clarity.
+smooth_mesh_interior(mesh, metric=rho, method="mmpde",
+                     boundary_slip=True)
+```
+
+```{tip}
+**`method="mmpde"` is the default mover** (since this release): the
+variational moving-mesh adaptation of Huang & Kamenski. It is
+dimension-general (2D/3D), matrix-free (no PETSc solve — small
+per-cell dense algebra plus a parallel `Vec` assembly), provably
+non-folding, and — uniquely among the movers here — genuinely
+*clusters and aligns* to an **anisotropic tensor** metric. It is
+both the most capable and the most straightforward to reason about,
+which is why it is now the default. Pass a **scalar** density (as
+above; it is promoted to the isotropic tensor `ρ·I`) or a `d×d`
+**tensor** metric (e.g. {py:func}`fault_metric_tensor`, or a
+`grad T`-aligned boundary-layer tensor) to get true thin-across /
+long-along refinement. Full design + derivation:
+{doc}`/developer/design/anisotropic-mmpde-mover`.
+
+The earlier movers remain available via `method=`:
+`"spring"` (fast volumetric equant-cell smoother), `"ma"`
+(isotropic Monge–Ampère), `"ot"` (linear OT-improvement step),
+`"anisotropic"` (decoupled-Winslow tensor smoother — reshapes but
+does not cluster). Use them only when you specifically need their
+behaviour; `"mmpde"` supersedes `"anisotropic"` for fault / front
+refinement.
+
+Key `mmpde` knobs (via `method_kwargs`): `p` (functional exponent,
+1.5–2), `theta` (Huang alignment/equidistribution balance, 1/3),
+`step_frac` (per-node move cap, 0.2), `tol` (scale-relative
+convergence exit), `metric_eval` (`"rbf"` default — fast baked
+metric interpolation). `boundary_slip=True` lets boundary nodes
+slide tangentially (needed for surface-reaching features).
 ```
 
 `metric_density_from_gradient` builds

docs/developer/design/anisotropic-mmpde-mover.md

@@ -0,0 +1,276 @@
+# Anisotropic MMPDE Mesh Mover (variational, Huang–Kamenski)
+
+**Status:** design + validated numpy prototype (2026-05-30); UW3 port
+pending as `smooth_mesh_interior(..., method="mmpde")`.
+
+## Why this mover exists
+
+The existing movers cannot produce a **thin refined strip aligned to a
+codimension-1 feature** (a fault, an interface):
+
+- **`"ma"` (Monge–Ampère, `_winslow_elliptic`)** genuinely *clusters*,
+  but only isotropically. A scalar value-metric peaked on a fault
+  refines a *disk around the densest part* — the cluster lands at the
+  metric's centre of gravity, not *on* the line (measured: only
+  ~60–80 % of refined cells within ~0.75 h of a fault, the rest pulled
+  toward the middle/root). Right tool for *tangentially-uniform*
+  features (thermal boundary layers, plumes); wrong *shape* for a fault.
+- **`"anisotropic"` (decoupled forward-Winslow, `_winslow_anisotropic`)**
+  uses a tensor metric but solves a *linear* M-weighted Laplacian for
+  each physical coordinate independently. That is a **smoother, not a
+  clusterer** — it reshapes/aligns cells but does not concentrate them
+  (measured band median-area/global ≈ 0.9 ≈ uniform). It also has no
+  non-folding guarantee: at fault-grade anisotropy (≳6:1) the map folds
+  one cell, and the *global* signed-area backtrack then throttles the
+  whole step to protect it, so the mesh freezes (`scale → 0`) while
+  reporting "converged" — the paradoxical "at the stability limit but
+  nothing moves." The inner solve is exact (MUMPS LU, SPD); the failure
+  is the *formulation/strategy*, not the linear solve.
+  See `ADAPTIVE_MESHING_DESIGN.md` and `mesh-adaptation-formulation.md`.
+
+The MMPDE mover fixes this structurally. It is the **variational moving
+mesh** method of Huang & Russell, in the direct simplex discretization
+of Huang & Kamenski. It generates the physical mesh as the image of a
+**fixed computational (reference) mesh** under the *inverse* coordinate
+map, minimizing a meshing functional that combines **equidistribution**
+and **alignment** to a metric tensor `M`. Because the functional has a
+barrier (`G → ∞` as `det 𝕁 → 0`) it is provably **non-folding** (Huang
+& Kamenski 2018), and because it is the inverse map of a convex
+computational domain it genuinely *clusters and aligns* — a thin strip
+*on* the feature line.
+
+Validated (numpy prototype): refinement **on the fault line**
+(on-fault fraction 0.95–0.99 vs 0.6–0.8 for `"ma"`), **0 crushed
+cells**, **monotone-convergent, never folds** (`scale = 1` throughout),
+generalizes to **multiple / crossing / curved** faults and to an
+**evolving (moving) metric**.
+
+## References
+
+- W. Huang & L. Kamenski, *A geometric discretization and a simple
+  implementation for variational mesh generation and adaptation*,
+  J. Comput. Phys. 301 (2015) 322–337. **doi:10.1016/j.jcp.2015.07.015**
+  (arXiv:1410.7872). — the discretization and the analytic nodal
+  velocities implemented here.
+- W. Huang & L. Kamenski, *On the mesh nonsingularity of the moving mesh
+  PDE method*, Math. Comp. 87 (2018) 1887–1911. **doi:10.1090/mcom/3271**
+  (arXiv:1512.04971). — the non-folding guarantee.
+- W. Huang, *Variational mesh adaptation: isotropy and equidistribution*,
+  J. Comput. Phys. 174 (2001) 903–924. **doi:10.1006/jcph.2001.6878** —
+  the functional and its `p`, `θ` parameters.
+- W. Huang & R. D. Russell, *Adaptive Moving Mesh Methods*, Springer AMS
+  174 (2011). **doi:10.1007/978-1-4419-7916-2**.
+
+## Formulation (general `d`; simplex meshes)
+
+Per element `K` with local physical vertices `x0, x1, x2` and the
+corresponding **fixed computational** vertices `ξ0, ξ1, ξ2`:
+
+```
+E    = [x1-x0,  x2-x0]          physical edge matrix (columns)
+Ehat = [ξ1-ξ0,  ξ2-ξ0]         computational edge matrix  (FIXED reference)
+𝕁    = Ehat · E^{-1}            Jacobian of the inverse map ξ(x)   (eq 17)
+r    = det 𝕁 = det Ehat / det E
+M    = M(x_K)                   SPD metric at the element centroid
+S    = tr(𝕁 M^{-1} 𝕁^T)
+```
+
+**Huang's functional** (eq 3; with `d = 2`, `dp/2 = p`):
+
+```
+G = θ · √det(M) · S^p  +  (1 - 2θ) · 2^p · r^p · det(M)^{(1-p)/2}
+I_h = Σ_K |K| · G_K                                              (eq 6)
+```
+
+`θ ∈ (0, ½]` balances **alignment** (1st term) vs **equidistribution**
+(2nd term); `p ≥ 1`. Coercive/polyconvex (unique minimizer) for
+`0 < θ ≤ ½`, `dp ≥ 2`, `p ≥ 1`.
+
+**Derivatives** (eq 16):
+
+```
+∂G/∂𝕁 = 2 p θ √det(M) · S^{p-1} · M^{-1} 𝕁^T
+∂G/∂r = p (1 - 2θ) 2^p · det(M)^{(1-p)/2} · r^{p-1}
+```
+
+**Physical-coordinate nodal velocity** (Appendix A, eqs 39–41). The
+descent velocity `v_i = −∂I_h/∂x_i` is assembled from per-element local
+velocities; for the local non-`0` vertices,
+
+```
+[v1; v2] = −G E^{-1} + E^{-1} (∂G/∂𝕁) Ehat E^{-1} + (∂G/∂r) r E^{-1}
+v0       = −(v1 + v2)
+∂I_h/∂x_i = − Σ_{K ∋ i} |K| v^K_{i}
+```
+
+### The metric-variation term is ESSENTIAL (key gotcha)
+
+Equations 39–41 as written treat `M = M(x_K)` as moving with the cell
+centroid, contributing a **`∂G/∂M : ∂M/∂x`** term. Dropping it
+("frozen-M") is *wrong wherever `M` varies sharply* — i.e. **on the
+fault**, which is exactly where it matters. Measured: frozen-M gradient
+is **65–330 % wrong** vs finite differences for a sharp fault metric,
+and the resulting flow does **not cluster** (band ≈ uniform, energy
+wanders). Including it restores agreement to **1e-8**.
+
+```
+∂G/∂M = θ √det(M) [ ½ S^q M^{-1} − q S^{q-1} M^{-1} 𝕁^T 𝕁 M^{-1} ]
+      + (1-2θ) d^q r^p · (1-p)/2 · det(M)^{(1-p)/2} · M^{-1}     (symmetric)
+```
+(general `d`, `q = dp/2`; for `d=2`, `q=p`, `d^q=2^p`.)
+
+assembled per vertex as `∂I_h/∂x_i += Σ_{K∋i} (|K|/(d+1)) · tr(∂G/∂M ·
+∂_c M)` (the `1/(d+1)` because `∂x_K/∂x_i = 1/(d+1)`), with `∂M/∂x` from
+the analytic metric (centroid finite difference is fine).
+
+**Lesson:** finite-difference-validate any hand-derived mesh gradient
+before trusting it. The prototype's `mmpde.py __main__` does exactly
+this (const-M and varying-M, all `p/θ`).
+
+## Time integration (the MMPDE)
+
+Gradient flow `∂ξ/∂t = −(P/τ) ∂I_h/∂ξ`, rewritten in physical
+coordinates (eq 39):
+
+```
+dx_i/dt = (P_i / τ) Σ_{K ∋ i} |K| v^K_i ,   P_i = det(M(x_i))^{(p-1)/2}
+```
+
+`P_i` (eq 24) makes the flow invariant under `M → cM` (scale-free
+node concentration). Discretized as **explicit Euler with two
+safeguards** (validated):
+
+1. **Per-node step cap**: limit each node's move to `step_frac · h_i`
+   (`h_i` = min incident edge). Prevents single-step overshoot
+   (the boundary-crush mechanism); `step_frac ≈ 0.2`.
+2. **Energy line-search backtrack**: accept a step only if it produces
+   **no fold** *and* **decreases `I_h`** (halve `scale` up to ~20×).
+   This makes the descent **monotone** — it reaches the true minimizer
+   instead of oscillating around it (an early non-monotone version
+   produced run-to-run-variable, over-stated refinement).
+
+`τ` sets the move scale; with the line-search, `τ` is non-critical
+(`τ = 1`). Convergence ≈ a few hundred explicit steps for `cellSize`
+0.04; the line-search crawls to small `scale` near the minimum (a
+candidate for acceleration in the port).
+
+## Boundary conditions
+
+- **Pinned** interior-only boundary: boundary nodes excluded from the
+  move (`free = ~is_bnd`). Simplest; the ring cannot open to admit a
+  surface-reaching feature.
+- **Tangential slip** (recommended for surface-reaching faults): include
+  boundary nodes in the move but **remove the outward-normal component**
+  of their velocity, then snap them back onto the surface (`project`,
+  e.g. fixed `|r|` for an annulus). `free = (~is_bnd) | slip`.
+  - Trade-off (measured, surface-reaching fault, `across=100`): slip
+    lets the ring **open** to admit the fault (finer band at the
+    surface, 0.30 → 0.26) **but** costs boundary-row angle quality
+    (min angle 20° → 9°, on-fault 0.88). It is a real knob, **not
+    free** — localize slip to the fault root and/or temper its strength.
+  - Use the projected PETSc/`Gamma_P1` normal in UW3 (the generic
+    boundary normal used for slip BCs), consistent with the existing
+    `_build_slip_projector`.
+
+## Behaviour and tuning (validated, numpy prototype)
+
+| Knob | Effect |
+|---|---|
+| `across`-ratio of `M` | primary strength; 9 → 100 deepens band 0.79 → 0.44, on-fault → 1.0, 0 crushed. `≳400` over-shoots (refinement drifts off-fault). Sweet spot ~100. |
+| `p` (with `θ`) | `p` 1.5 → 3 sharpens the band; pair higher `p` with smaller `θ` (e.g. `θ = 1/6`). |
+| node count (base `h`) | sets **absolute** on-fault cell size (≈ linear in `h`): `cellSize` 0.04 → 0.013 gives `h_fault` 0.017 → 0.006. Use to get real resolution; the *ratio* is capped by the fixed budget (~1.5–2× median, the standard r-adapt cap). |
+| cumulative reference-reset | re-running with `ref = current mesh` pushes past the single-run cap (band 0.62 → 0.19 over 3 rounds) but **degrades quality** (min angle 24° → 0.9°, crushed cells appear). Use sparingly. |
+
+**Discriminant:** judge with `n_crushed` (cells with area < 0.02 ·
+global median) and the metric-aware *radial/tangential* extent — **not**
+min-angle, which over-counts legitimate thin anisotropic cells (a
+resolved strip looks like "slivers" to an isotropic detector). See
+`alignment.py` in the prototype scratch.
+
+## UW3 port plan (`method="mmpde"`)
+
+**Architectural rule: PETSc-native and parallel-safe by construction.**
+The numpy prototype was a *validation* vehicle only. The element-level
+algebra (per-`K` `d×d` matrices `E`, `Ehat`, `𝕁`, and `G`, `∂G/∂𝕁`,
+`∂G/∂r`, `∂G/∂M`) is genuinely local and may stay vectorised NumPy over
+the rank-local element block — that is not a parallel hazard. Everything
+that **couples across vertices or ranks** must go through PETSc `Vec` /
+DM operations, never a rank-local `np.add.at` into a global array. The
+prototype's `np.add.at` assembly is serial-only and must NOT be ported
+as-is.
+
+1. Add `_winslow_mmpde(mesh, metric, pinned_labels, verbose, **kw)` to
+   `src/underworld3/meshing/smoothing.py`, **dimension-general (`d = 2`
+   and `3`) from the start** — the method and all formulas above are
+   general `d` (paper validates 3D), and UW3 already has the 3D
+   infrastructure (`_tet_cells`, `_signed_volumes`, and a 3D branch in
+   `_boundary_vertex_normals` / `_build_slip_projector` for tangent-plane
+   slip). Use `cdim` everywhere: `(d+1)`-vertex cells, `d×d` edge
+   matrices / metric, `det`/`inv` via batched `numpy.linalg` (not a
+   hand-coded `2×2`), signed *volume* via `_tet_cells`/`_signed_volumes`
+   in 3D. Do **not** raise `NotImplementedError` for 3D the way
+   `_winslow_anisotropic` does — that 2D-only limitation is what this
+   mover supersedes. Element terms (eqs 16, 40–41 + `∂G/∂M`, with
+   `q = d·p/2`) computed per rank-local cell; the **velocity assembly is
+   a PETSc Vec**, not a numpy array:
+   - assemble `Σ_{K∋i} |K| v^K_i` into a global `Vec` with `ADD_VALUES`
+     (DM section / `petsc_dm.localToGlobal(..., ADD_VALUES)`), so cells
+     straddling a partition boundary correctly contribute to off-rank
+     vertices. This is the same ghost-summation pattern flagged for the
+     lumped-V source in `_winslow_elliptic` (whose `np.add.at` is a known
+     serial-only TODO) — do it right here from the start.
+   - the `P_i` balancing and the final coordinate update act on the
+     assembled global Vec, then scatter back to the local (ghosted)
+     coordinate vector.
+2. **All scalar tests/norms are collective reductions** (`uw.mpi.comm`
+   `allreduce`):
+   - energy `I_h = Σ_K |K| G_K` — sum over owned cells then `allreduce`
+     (count each cell once; use the owned-cell mask, not ghosts);
+   - the line-search predicates — *"min signed area > floor"* and *"I_h
+     decreased"* — must be **global** (`MIN` / the globally-summed `I_h`),
+     so every rank takes the same accept/backtrack branch in lockstep
+     (otherwise ranks desynchronise on `scale`);
+   - the convergence norm `max|Δx|` — `allreduce(MAX)`.
+3. Reuse existing parallel-aware infrastructure: `_tri_cells`,
+   `_signed_areas`, `_min_incident_edge`, and `_ot_adapt._build_slip_projector`
+   / `_resolve_slip` for the slip normal (`Gamma_P1`). The slip
+   projection is per-vertex local (no coupling) once the velocity Vec is
+   assembled and ghost-updated. `Gamma_P1` is already projected/parallel.
+4. **Metric `M` and `∂M/∂x`** via `uw.function.evaluate` at element
+   centroids (already parallel-aware) — do **not** hand-roll a coordinate
+   loop. `M` is a `d×d` sympy matrix / `VarType.TENSOR` MeshVariable via
+   the existing `supplied_D`-style entry routed by `smooth_mesh_interior`.
+   Eulerian re-eval each step is safe (`M` anchored to fixed feature
+   geometry).
+5. The **fixed computational reference** = mesh coordinates at the first
+   call, cached as a *ghosted* coordinate Vec on the mesh (like
+   `_ot_adapt_reference_coords`), so each rank has its halo. For an
+   *evolving* feature, keep the uniform reference and re-relax each adapt
+   event (validated serially: tracks a moving fault cleanly).
+6. `method_kwargs`: `p` (1.5–2), `theta` (1/3), `tau` (1), `n_steps`,
+   `step_frac` (0.2), `slip` (bool/mask), `area_floor_frac` (0.01).
+7. Cross-link from `ADAPTIVE_MESHING_DESIGN.md` /
+   `mesh-adaptation-formulation.md`. **Regressions must cover both
+   dimensions and parallel** (decision 2026-05-30: port 2D+3D, validate
+   both directly in UW3 — no separate 3D numpy prototype):
+   - **2D serial** (Tier-A): uniform `M` ⇒ near no-op; single fault ⇒
+     on-fault band, 0 crushed.
+   - **3D serial**: uniform `M` ⇒ near no-op on a tet mesh; a planar
+     fault ⇒ on-plane refined slab, 0 crushed. 3D is *derived* here but
+     **not yet numerically validated**, and its decoupled non-folding
+     margin is tighter than 2D — treat this as a first-class acceptance
+     test, not an afterthought.
+   - **`np>=2` in each dimension**: matches the serial result to solver
+     tolerance (same final coords up to partition-independent reduction
+     order) — the assembly/ghost path is exactly where serial-only bugs
+     hide.
+
+## Open items
+
+- Acceleration: the line-search takes tiny end-steps near convergence;
+  an accelerated / semi-implicit step could cut iteration count. Any such
+  scheme must keep its global-reduction predicates collective (item 2).
+- Slip localization: temper/localize slip to the fault root to keep the
+  finer-at-surface benefit without the global boundary-angle cost.
+- Parallel correctness is a **release gate**, not an open item: the port
+  is not "done" until the `np>=2` regression (item 7) matches serial.