3D node redistribution: the MMPDE mover drives tetrahedral meshes (capstone round 2)

docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md

@@ -573,6 +573,29 @@ out of scope here, as in 2D.
   metric-variation term is 2·d analytic metric evaluations per iteration —
   measure, and if hot, evaluate the Shepard-baked metric's gradient instead.
 
+**Round-2a status: DONE (2026-07-23, worktree `mmpde-3d`).** The rim
+locks fell as predicted (`_tet_cells` / new `_signed_volumes` /
+`fact = d!`; guards flipped; `mesh_metric_mismatch` and the
+`follow_metric` polish gained their tet forms). The one REAL blocker
+was not on the list: `_pinned_mask` stopped its label closure at
+edges, and 3D gmsh labels tag *faces only*, so on any 3D mesh **no
+boundary vertex was classified at all** — nothing pinned, nothing
+slipped, the whole boundary drifted freely into the interior. Every
+tagged non-vertex point now closes down to its vertices (bit-identical
+in 2D, where an edge's closure is exactly its endpoints). With that
+fixed, first light on the dipping-fault-plane box: monotone energy
+descent, zero folds, boundary nodes held **exactly** on their faces
+under slip, and h-grading ratio 1.12 vs 1.30 for the identical 2D
+problem — same family, compressed by the dimension exponent (the
+ideal equidistribution ratio for this metric is 4× in 2D but only
+2.5× in 3D, h ∝ ρ^(−1/d)); repeated calls hold a stable equilibrium,
+matching the 2D "maintains, does not compound" behaviour.
+`follow_metric` runs end-to-end on 3D. Full mover/adapt family sweep:
+215 green serially; the negative 3D-raises contract tests flipped to
+capability tests (`test_0764`, `test_0850_mesh_smoothing`). Spherical
+slip validation deferred to round 3b with the rest of the curved-
+boundary work.
+
 ### Round 2b — 3D MMPDE, parallel (same worktree, gated separately)
 
 The mover's parallel machinery (coordinate-DM `localToGlobal(ADD_VALUES)`
@@ -583,6 +606,15 @@ bit-identical; the known ~1e-4%-level step-cap partition drift documented in
 `mmpde.py:487-495` applies unchanged). Any 3D-specific divergence is a bug to
 fix, not a new mechanism to build.
 
+**Round-2b status: DONE (2026-07-23).** Gates, not code, exactly as
+predicted — zero 3D-specific changes were needed. The dipping-plane
+box gate at np=1/2/4: runs clean at every rank count (exercising the
+post-#379 collective sync machinery in 3D), zero folds everywhere,
+boundary nodes held exactly on their faces, and the graded
+near/far h medians agree across partitions to ≤2.5% each (grading
+ratio 1.121 / 1.088 / 1.077) — the same few-percent equilibrium
+drift the rank-local step cap produces in the 2D combination gate.
+
 ### Round 3a — unified adapt + redistribute workflow
 
 Composition semantics need one design commitment: **redistribute-then-adapt**

src/underworld3/discretisation/discretisation_mesh.py

@@ -3119,11 +3119,8 @@ 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.
+        # _pinned_mask closes every tagged point (edge in 2D, face in 3D)
+        # down to its vertices, so face-only 3D labels classify correctly.
         is_bnd = _pinned_mask(dm, all_labels)
 
         slip_labels, free_labels = self._resolve_slip_spec(slip_spec)
@@ -6159,12 +6156,12 @@ def redistribute_nodes(self, metric, *, verbose=False, **kwargs):
         mesh in place.
 
         This method is how each mesh type controls whether (and how) it
-        can be modified: the base implementation supports **2D simplex
-        (triangle) meshes**, where it drives the Huang–Kamenski
-        variational MMPDE mover (non-folding by construction,
-        parallel-safe; scalar metric → isotropic equidistribution,
-        tensor metric → anisotropic clustering and alignment).
-        Quadrilateral / hexahedral meshes, 3D meshes and constrained
+        can be modified: the base implementation supports **2D
+        (triangle) and 3D (tetrahedral) simplex meshes**, where it
+        drives the Huang–Kamenski variational MMPDE mover (non-folding
+        by construction, parallel-safe; scalar metric → isotropic
+        equidistribution, tensor metric → anisotropic clustering and
+        alignment). Quadrilateral / hexahedral meshes and constrained
         manifolds raise ``NotImplementedError`` — no mover is
         implemented for them yet.
 
@@ -6202,16 +6199,16 @@ def redistribute_nodes(self, metric, *, verbose=False, **kwargs):
                 f"manifold meshes (mesh.dim={self.dim} != mesh.cdim="
                 f"{self.cdim}): every node would have to be constrained "
                 "to the surface. Implemented today: 2D simplex (triangle) "
-                "meshes, via the MMPDE mover.")
-        if not bool(self.dm.isSimplex()) or self.cdim != 2:
+                "and 3D simplex (tetrahedral) meshes, via the MMPDE mover.")
+        if not bool(self.dm.isSimplex()) or self.cdim not in (2, 3):
             kind = "simplex" if bool(self.dm.isSimplex()) else "tensor-product (quad/hex)"
             raise NotImplementedError(
                 f"node redistribution is not implemented for {self.cdim}D "
                 f"{kind} meshes. Implemented today: 2D simplex (triangle) "
-                "meshes, via the MMPDE mover (its 3D / quad discretization "
-                "does not exist yet). To add resolution instead, use "
-                "mesh.adapt(metric_field, max_levels=...) — a topology "
-                "change.")
+                "and 3D simplex (tetrahedral) meshes, via the MMPDE mover "
+                "(no quad/hex discretization exists). To add resolution "
+                "instead, use mesh.adapt(metric_field, max_levels=...) — "
+                "a topology change.")
         from underworld3.meshing.smoothing import smooth_mesh_interior
         smooth_mesh_interior(self, metric=metric, method="mmpde",
                              verbose=verbose, **kwargs)

src/underworld3/meshing/smoothing/api.py

@@ -135,9 +135,9 @@ def smooth_mesh_interior(
         Metric-grading solver (ignored when ``metric is None``).
         ``"mmpde"`` (alias ``"variational"``) — variational
         moving-mesh (Huang–Kamenski MMPDE) with a full tensor (or
-        scalar) metric; non-folding by construction. **Currently
-        2D-only** (triangle meshes) — a 3D mesh raises
-        ``NotImplementedError``.
+        scalar) metric; non-folding by construction. Triangle (2D)
+        and tetrahedral (3D) simplex meshes; quad/hex meshes and
+        constrained manifolds raise ``NotImplementedError``.
 
         The historical spellings ``"spring"``, ``"ma"``, ``"ot"``
         and ``"anisotropic"`` name interior movers retired in
@@ -796,17 +796,33 @@ def _polish(moved):
         # stays intact while sliver cells get rounded out.
         # `polish_max_iters=0` disables entirely.
         if moved and polish_max_iters > 0:
-            tris_polish = _tri_cells(mesh.dm)
+            if mesh.cdim == 2:
+                cells_polish = _tri_cells(mesh.dm)
+            else:
+                from .graph import _tet_cells, _signed_volumes
+                cells_polish = _tet_cells(mesh.dm)
             for _polish_iter in range(int(polish_max_iters)):
-                # Check current shape quality
-                p = np.asarray(mesh.X.coords)[tris_polish]
-                e0 = np.linalg.norm(p[:, 1] - p[:, 0], axis=1)
-                e1 = np.linalg.norm(p[:, 2] - p[:, 1], axis=1)
-                e2 = np.linalg.norm(p[:, 0] - p[:, 2], axis=1)
-                A = np.abs(_signed_areas(np.asarray(mesh.X.coords),
-                                           tris_polish))
-                q = (4.0 * np.sqrt(3.0) * A
-                     / (e0 * e0 + e1 * e1 + e2 * e2 + 1.0e-30))
+                # Current worst cell-shape quality, normalised so the
+                # regular simplex scores q=1 and a degenerate sliver
+                # q→0. 2D: q = 4√3·A/Σe²  (Σ over the 3 edges).
+                # 3D: q = 6√2·V/ē³ with ē = rms of the 6 edge lengths
+                # (the volume of a regular tet of edge e is e³/(6√2)).
+                X = np.asarray(mesh.X.coords)
+                p = X[cells_polish]
+                nv = p.shape[1]
+                e2sum = np.zeros(p.shape[0])
+                n_edges = 0
+                for _a in range(nv):
+                    for _b in range(_a + 1, nv):
+                        e2sum += np.sum((p[:, _b] - p[:, _a]) ** 2, axis=1)
+                        n_edges += 1
+                if mesh.cdim == 2:
+                    A = np.abs(_signed_areas(X, cells_polish))
+                    q = 4.0 * np.sqrt(3.0) * A / (e2sum + 1.0e-30)
+                else:
+                    V = np.abs(_signed_volumes(X, cells_polish))
+                    ebar = np.sqrt(e2sum / n_edges)
+                    q = 6.0 * np.sqrt(2.0) * V / (ebar ** 3 + 1.0e-30)
                 q_min = _global_min(q.min())
                 if verbose:
                     uw.pprint(
@@ -831,9 +847,10 @@ def node_redistribution(mesh, metric, *, verbose=False, **kwargs):
     partition preserved. This is a thin dispatch to
     ``mesh.redistribute_nodes(metric, ...)`` — the mesh type controls
     whether (and how) it can be modified. The base implementation
-    supports 2D simplex (triangle) meshes via the Huang–Kamenski
-    variational MMPDE mover; unsupported mesh types (quad/hex, 3D,
-    manifolds) raise ``NotImplementedError`` stating what exists.
+    supports 2D (triangle) and 3D (tetrahedral) simplex meshes via the
+    Huang–Kamenski variational MMPDE mover; unsupported mesh types
+    (quad/hex, manifolds) raise ``NotImplementedError`` stating what
+    exists.
 
     Parameters
     ----------

src/underworld3/meshing/smoothing/graph.py

@@ -101,7 +101,6 @@ def _pinned_mask(dm, pinned_labels):
     ``Centre`` pressure-pin marker on an Annulus, whose underlying
     ``DMLabel`` has no strata and hard-crashes any query)."""
     pStart, pEnd = dm.getDepthStratum(0)
-    eStart, eEnd = dm.getDepthStratum(1)
     n_verts = pEnd - pStart
     is_pinned = np.zeros(n_verts, dtype=bool)
     for lname in pinned_labels:
@@ -127,10 +126,15 @@ def _pinned_mask(dm, pinned_labels):
                 if pStart <= idx < pEnd:
                     # Tagged vertex — pin directly.
                     is_pinned[idx - pStart] = True
-                elif eStart <= idx < eEnd:
-                    # Tagged edge — pin both endpoint vertices.
-                    cone = dm.getCone(idx)
-                    for c in cone:
+                else:
+                    # Tagged edge (2D labels) or face (3D labels tag
+                    # faces ONLY, so stopping at edges left every 3D
+                    # boundary vertex unmarked — no pin, no slip, the
+                    # boundary drifted freely): close the point down to
+                    # its vertices. For an edge the closure is exactly
+                    # its endpoints, so the 2D mask is bit-identical.
+                    closure = dm.getTransitiveClosure(idx)[0]
+                    for c in closure:
                         if pStart <= c < pEnd:
                             is_pinned[c - pStart] = True
     return is_pinned
@@ -237,6 +241,20 @@ def _signed_areas(coords, tris):
                   - (b[:, 1] - a[:, 1]) * (c[:, 0] - a[:, 0]))
 
 
+def _signed_volumes(coords, tets):
+    """Signed volume of each tetrahedron (sign = orientation) — the 3D
+    analogue of :func:`_signed_areas`, indexed the same way. Note the
+    DMPlex reference tet is det-NEGATIVE, so a healthy all-DMPlex-ordered
+    tet mesh yields uniformly negative values; callers (the MMPDE mover's
+    fold guard) already orientation-normalise with the global median
+    sign, exactly as they do for 2D."""
+    a = coords[tets[:, 0]]
+    e = np.stack([coords[tets[:, 1]] - a,
+                  coords[tets[:, 2]] - a,
+                  coords[tets[:, 3]] - a], axis=1)
+    return np.linalg.det(e) / 6.0
+
+
 
 
 
@@ -441,8 +459,8 @@ def _blend(f):
 def _tet_cells(dm):
     """Tetrahedron vertex-index quadruples (local-chart), or ``None`` if the
     mesh is not all-tet. The 3D analogue of :func:`_tri_cells` — used by the
-    3D boundary-face extraction in :func:`_boundary_facets` (the MMPDE
-    mover itself is currently 2D-only)."""
+    3D boundary-face extraction in :func:`_boundary_facets` and by the
+    MMPDE mover's 3D discretization."""
     cStart, cEnd = dm.getHeightStratum(0)
     pStart, pEnd = dm.getDepthStratum(0)
     tets = []

src/underworld3/meshing/smoothing/metrics.py

@@ -98,7 +98,8 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None):
     Parameters
     ----------
     mesh : underworld3.discretisation.Mesh
-        Triangle mesh (only 2-D for now).
+        Triangle (2D) or tetrahedral (3D) mesh; "area" throughout
+        this docstring reads as the cell measure (volume in 3D).
     metric : sympy / UW expression
         The target *density* ρ (larger ⇒ finer cells) — same
         object you would pass to ``smooth_mesh_interior``.
@@ -134,11 +135,19 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None):
     import underworld3 as _uw
 
     coords = np.asarray(mesh.X.coords)
-    tris = _tri_cells(mesh.dm)
+    if mesh.cdim == 2:
+        tris = _tri_cells(mesh.dm)
+        A_actual = None if tris is None else np.abs(
+            _signed_areas(coords, tris))
+    else:
+        from .graph import _tet_cells, _signed_volumes
+        tris = _tet_cells(mesh.dm)
+        A_actual = None if tris is None else np.abs(
+            _signed_volumes(coords, tris))
     if tris is None:
         raise NotImplementedError(
-            "mesh_metric_mismatch: triangle mesh required")
-    A_actual = np.abs(_signed_areas(coords, tris))
+            "mesh_metric_mismatch: simplex (triangle/tetrahedral) mesh "
+            "required")
     centroids = coords[tris].mean(axis=1)
     rho = np.asarray(_uw.function.evaluate(
         metric, centroids)).reshape(-1)

src/underworld3/meshing/smoothing/mmpde.py

@@ -1,15 +1,17 @@
 """The Huang–Kamenski variational MMPDE mover (recommended for
-production adaptive meshing; 2D-only). See the package docstring for
-the module map.
+production adaptive meshing; triangle and tetrahedral meshes). See the
+package docstring for the module map.
 """
 
+import math
 import warnings
 
 import numpy as np
 
 import underworld3 as uw
 
-from .graph import (_tri_cells, _signed_areas, _owned_cell_mask,
+from .graph import (_tri_cells, _tet_cells, _signed_areas, _signed_volumes,
+                    _owned_cell_mask,
                     _owned_vertex_mask, _min_incident_edge_nd,
                     _global_sum, _global_min, _global_max, _global_mean)
 
@@ -94,10 +96,10 @@ def _mmpde_mover(mesh, metric, pinned_labels, verbose,
                    **_unknown_kwargs):
     r"""Anisotropic variational moving-mesh adaptation (Huang–Kamenski
     MMPDE; the direct simplex discretization of JCP 301 (2015) 322,
-    arXiv:1410.7872). **2D (triangle meshes) only** and parallel-safe.
-    The underlying method is dimension-general, but the 3D (tetrahedral)
-    discretization has not been implemented — a 3D mesh raises
-    ``NotImplementedError`` immediately.
+    arXiv:1410.7872). Triangle (2D) and tetrahedral (3D) meshes,
+    parallel-safe. The per-element algebra is written over ``cdim``
+    throughout — the two dimensions differ only in the cell list, the
+    signed-measure primitive and the ``d!`` volume factor.
 
     Generates the physical mesh as the image of a **fixed computational
     (reference) mesh** under the inverse coordinate map, minimizing
@@ -160,14 +162,12 @@ def _mmpde_mover(mesh, metric, pinned_labels, verbose,
             stacklevel=2)
     pinned_labels = tuple(pinned_labels)
     cdim = mesh.cdim
-    if cdim != 2:
-        # Guard here, before any metric parsing or DM work, so a 3D caller
-        # gets an honest message rather than a NameError from the (never
-        # implemented) 3D discretization deeper in the mover (READ-01).
+    if cdim not in (2, 3):
+        # Guard here, before any metric parsing or DM work, so the caller
+        # gets an honest message rather than a failure deeper in the mover.
         raise NotImplementedError(
-            "MMPDE mesh movement is currently 2D-only (triangle meshes): "
-            "the 3D tetrahedral discretization of the mover has not been "
-            f"implemented. Got a mesh with cdim={cdim}.")
+            "MMPDE mesh movement supports 2D (triangle) and 3D "
+            f"(tetrahedral) simplex meshes. Got a mesh with cdim={cdim}.")
     p = float(p); theta = float(theta); tau = float(tau)
     q = cdim * p / 2.0
     dq = float(cdim) ** q
@@ -230,14 +230,19 @@ def _dM_dx(cen):
     dm = mesh.dm
     pStart, pEnd = dm.getDepthStratum(0)
     n_verts = pEnd - pStart
-    # cdim == 2 is guaranteed by the guard at the top of this function.
-    # (The former 3D branch here referenced `_signed_volumes`, which was
-    # never implemented — READ-01.)
-    cells_all = _tri_cells(dm)
-    signed_vol = _signed_areas
+    # cdim in (2, 3) is guaranteed by the guard at the top; the rest of
+    # the mover is written over cdim and needs only the right cell list,
+    # signed-measure primitive and d! volume factor.
+    cells_all = _tri_cells(dm) if cdim == 2 else _tet_cells(dm)
+    signed_vol = _signed_areas if cdim == 2 else _signed_volumes
     if cells_all is None:
+        # TODO(BUG): rank-LOCAL early return — under MPI a rank with zero
+        # local cells (or a non-simplex local patch) returns here while the
+        # other ranks enter the collective energy/line-search reductions
+        # below, hanging the job. Pre-existing 2D posture, inherited by 3D
+        # (2026-07 review note); the fix is a collective emptiness check.
         return
-    fact = 2.0                                 # d! → |K| = |detE|/d!
+    fact = float(math.factorial(cdim))         # d! → |K| = |detE|/d!
     owned_cell = _owned_cell_mask(dm)
     cells_own = cells_all[owned_cell]
     is_owned_v = _owned_vertex_mask(dm)