From f7efbb1f6bbf0a8dea6393be2e3198e193bde242 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 09:22:07 +1000 Subject: [PATCH 1/7] fix(nvb): re-seeding the tagged state must start from a fresh label MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit DMLabelSetValue does not remove a point from its previous stratum, so overwriting an existing seed left every cell in TWO strata and readers got the OLD value back. The redistribute-then-adapt composition on a 3D mesh re-seeds the refinement state on the moved geometry — the driver then silently ran the moved mesh with the stale unmoved-order seed and the conforming drain deadlocked ('UWNVB drain stalled', found by the round-3a combination experiment; a 72-state decode sweep and a 3456-case two-tet pair sweep against the Python engine cleared the C driver itself first). write_tagged_state_label now destroys and recreates the label. Regression test locks the re-seed against a single fresh seed of the same geometry. With this, the 3D mover+adapt combination runs end to end: 33,570 cells vs 36,040 adapt-only on the dipping-fault box (finer on-fault, over-refinement fraction 46 -> 39 percent), MG 3 its, exact solution 9e-10. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/nvb.py | 11 ++++++-- tests/test_0840_nvb_3d_serial_adapt.py | 38 ++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/src/underworld3/utilities/nvb.py b/src/underworld3/utilities/nvb.py index 2860fb3bc..b76c9cd86 100644 --- a/src/underworld3/utilities/nvb.py +++ b/src/underworld3/utilities/nvb.py @@ -70,8 +70,15 @@ def write_tagged_state_label(dm): perms = {p: i for i, p in enumerate(itertools.permutations(range(4)))} cS, cE = dm.getHeightStratum(0) vS, vE = dm.getDepthStratum(0) - if not dm.hasLabel(TAGGED_STATE_LABEL): - dm.createLabel(TAGGED_STATE_LABEL) + # RE-seeding must start from a FRESH label: DMLabelSetValue does not + # remove a point from its previous stratum, so overwriting an existing + # seed leaves every cell in TWO strata and readers get the OLD value + # back — the driver then runs the new geometry with the stale seed and + # the conforming drain deadlocks (found by the moved-base composition, + # 2026-07 round 3a). + if dm.hasLabel(TAGGED_STATE_LABEL): + dm.removeLabel(TAGGED_STATE_LABEL) + dm.createLabel(TAGGED_STATE_LABEL) label = dm.getLabel(TAGGED_STATE_LABEL) comm = dm.comm.tompi4py() diff --git a/tests/test_0840_nvb_3d_serial_adapt.py b/tests/test_0840_nvb_3d_serial_adapt.py index 784438596..4e6cb8024 100644 --- a/tests/test_0840_nvb_3d_serial_adapt.py +++ b/tests/test_0840_nvb_3d_serial_adapt.py @@ -95,6 +95,44 @@ def test_bounded_closure_single_cell_is_local(): assert added < 200, f"closure drained: +{added} cells for one mark" +def test_reseed_after_deform_is_fresh(): + """DMLabelSetValue does NOT remove a point from its previous stratum, + so RE-seeding the tagged state over an existing label left every cell + in two strata and readers got the OLD (pre-deform) state back — the + driver then ran the moved geometry with the stale seed and the + conforming drain deadlocked (the redistribute-then-adapt composition + stall, round 3a). write_tagged_state_label must produce a label whose + per-cell values match a single fresh seed of the same geometry.""" + from underworld3.utilities.nvb import (write_tagged_state_label, + TAGGED_STATE_LABEL) + + base = _base3(cellSize=0.6) + dm = base.dm_hierarchy[-1].clone() + write_tagged_state_label(dm) # first seed (original coords) + + # deform enough to change the coordinate-lex coloring order + v = dm.getCoordinatesLocal().duplicate() + arr = dm.getCoordinatesLocal().array.copy().reshape(-1, 3) + interior = np.all((arr > 0.05) & (arr < 0.95), axis=1) + rng = np.random.default_rng(7) + arr[interior] += 0.08 * (rng.random((int(interior.sum()), 3)) - 0.5) + v.array[:] = arr.ravel() + dm.setCoordinatesLocal(v) + write_tagged_state_label(dm) # RE-seed on moved geometry + + fresh = dm.clone() + fresh.removeLabel(TAGGED_STATE_LABEL) + write_tagged_state_label(fresh) # single seed, same geometry + + lab, flab = dm.getLabel(TAGGED_STATE_LABEL), fresh.getLabel( + TAGGED_STATE_LABEL) + cS, cE = dm.getHeightStratum(0) + stale = [c for c in range(cS, cE) + if lab.getValue(c) != flab.getValue(c)] + assert not stale, (f"{len(stale)} cells read a stale state after " + "re-seed (double-stratum label)") + + def test_shape_regularity_bounded_similarity_classes(): # The full-depth version of this bound (2.87M cells, 9 passes) was # proven by the stage-1a oracle; the CI gate only needs the plateau From 06ca8364d889fd347ecb4771ecbc107684b26c6b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 09:33:07 +1000 Subject: [PATCH 2/7] =?UTF-8?q?fix(custom-mg):=20a=20square=20finest=20tra?= =?UTF-8?q?nsfer=20is=20legitimate=20=E2=80=94=20reject=20only=20true=20si?= =?UTF-8?q?ze=20mismatches?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The mesh-owned FMG auto-pickup rejected any finest transfer with pc >= pr as 'no coarsening'. But a refinement generation whose new vertices all land on a Dirichlet boundary adds only CONSTRAINED dofs, so the free-dof counts of the last two levels coincide and the transfer is square — redundant but correct. Boundary-focused metrics (exactly what curved domains invite) hit this every time: every annulus boundary-layer adapt was silently downgraded to the default preconditioner. With the guard relaxed to genuine row/operator mismatches, the annulus child solves with custom MG in 2 iterations (identical solution to the GAMG reference), snapped or unsnapped. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/custom_mg.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 9023f610f..ccd9b07cd 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -888,7 +888,15 @@ def auto_inject_custom_mg(solver, field_id=None): solver.snes.setUp() op_n = int(solver.snes.getJacobian()[0].getSize()[0]) pr, pc = (int(v) for v in Ps[-1].getSize()) - if op_n > 0 and (pr != op_n or pc >= pr): # rows!=op or no coarsening + # Only a genuine size mismatch disqualifies the transfer. A + # SQUARE finest transfer (pc == pr) is legitimate and common on + # boundary-focused refinement: a generation whose new vertices + # all land on a Dirichlet boundary adds only CONSTRAINED dofs, + # so the free-dof counts of the last two levels coincide — the + # level is redundant but correct, and rejecting it silently + # downgraded every curved-domain boundary-layer adapt to the + # default preconditioner (round-3b annulus finding, 2026-07). + if op_n > 0 and pr != op_n: import warnings warnings.warn( "custom_mg: mesh-owned adapt-mesh FMG transfer is incompatible " From f79ea750a897d073cf6b750ced0d8d21a2b82e94 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 09:37:26 +1000 Subject: [PATCH 3/7] =?UTF-8?q?docs(design):=20round=203a/3b=20status=20?= =?UTF-8?q?=E2=80=94=203D=20combination=20works;=20shell=20+=20chord/snap?= =?UTF-8?q?=20probes=20measured?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Underworld development team with AI support from Claude Code --- .../design/ADAPTIVITY_3D_SPHERICAL_2026-07.md | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md b/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md index 933a847a3..db8793792 100644 --- a/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md +++ b/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md @@ -626,6 +626,26 @@ worked example (3D box or annulus fault: metric-driven redistribution of the base + NVB band on the fault + Stokes FMG + advection-diffusion with field transfer across re-adaptation), plus a short how-to in `docs/advanced/`. +**Round-3a status: 3D combination WORKS (2026-07-24, worktree +`adaptivity-round3`).** First run stalled the refinement drain — and the +cause was a PETSc label-semantics trap, not the driver: `DMLabelSetValue` +does not remove a point from its previous stratum, so RE-seeding the +tagged refinement state on the moved base left every cell in two strata +and readers got the *old* (unmoved-order) state back — the driver was +silently running the moved geometry with a stale seed. (The driver +itself was first cleared by exhaustive oracle sweeps: all 72 single-tet +states × 3 passes, and all 3,456 engine-valid two-tet state pairs with +single-cell marking — native vs Python engine, all equal.) +`write_tagged_state_label` now destroys and recreates the label; +regression test in `test_0840`. Result on the dipping-fault box (ml=2): +mover+adapt gives **33,570 cells vs 36,040 adapt-only**, slightly finer +on-fault (med h 0.0282 vs 0.0292) with the over-refinement fraction +down from 46% to 39% — the same budget-spends-better behaviour as 2D. +MG on the combined child: 3 iterations, exact solution to 9e-10. +Renders: `~/+Simulations/nvb_3d_adapt_evaluation/combo3d_*.png`. The +worked example + `docs/advanced/` how-to remain the landing +deliverable. + ### Round 3b — spherical geometry + curved/internal-boundary clean-ups - **MMPDE on shells**: `SphericalShell` (solid, `dim==cdim==3`) is a valid @@ -651,6 +671,34 @@ transfer across re-adaptation), plus a short how-to in `docs/advanced/`. path is boundary-label-driven). Validation: shell response benchmark (test_1064 pattern) on an NVB child. +**Round-3b probe results (2026-07-24):** + +* *MMPDE on `SphericalShell`: validated, zero new code.* With the + registered radial `BoundingSurface`s, a boundary-layer metric toward + the inner sphere gives zero folds and holds BOTH spheres to machine + precision (radius drift ≤ 2e-16) while the radial cell-size profile + grades to ratio ~1.32 (better than the box's 1.12 — the boundary-layer + metric shape suits the mover well). +* *The chord error is real and measured*: adapted annulus children carry + an outer-boundary radius error of 1.8e-3 at base h≈0.125 + (the h²/8R chord sag), frozen at base resolution regardless of depth. +* *Snap is safe on every gate*: projecting the child's boundary vertices + onto the analytic circles gives zero folds, radius error → 2e-16, and + the custom-P MG on the snapped child is untouched (2 iterations, + solution identical to the GAMG reference) — the coordinate-built + transfers tolerate the slightly non-nested snapped boundary exactly as + designed. **The chord-vs-snap ruling is now a pure decision** — no + technical blocker either way. +* *A real MG gap found and fixed on the way*: the mesh-owned FMG + auto-pickup rejected any square finest transfer as "no coarsening" — + but a generation whose new vertices all land on a Dirichlet boundary + adds only constrained dofs, making the free-dof counts of the last two + levels coincide. Boundary-focused metrics (what curved domains invite) + hit this every time, so every annulus boundary-layer adapt was + silently falling back to the default preconditioner. Guard relaxed to + genuine row/operator mismatches; annulus children now solve with + custom MG in 2 iterations. + ## Effort and risk, honestly | round | new-code size | risk | notes | From 671a452f9bdc9fa3e3925914dd59c1b1d9fe11dc Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 10:37:05 +1000 Subject: [PATCH 4/7] feat(adapt): snap curved-boundary vertices to analytic surfaces at EVERY generation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bisection puts new boundary vertices at chord midpoints, freezing a curved domain's boundary geometry at base resolution no matter how deep the refinement (measured: radius error h_base^2/8R ~ 2e-3 on the probe annulus, unchanged by depth). Per the 2026-07-24 ruling, every refinement generation now snaps its new boundary vertices onto the mesh's registered analytic bounding surfaces (radial/plane BoundingSurface restore) — not just the returned child — so each intermediate level is a valid mesh in its own right, extractable for solvers, and the metric marks on true geometry. The serial cell-list engine snaps its OWN coordinates (they feed the next generation's marking and midpoints); the native and SBR paths snap each level DM. On plane surfaces the chord midpoint already lies in the plane, so flat-domain results are bit-unchanged (confluence gates green). Adapted annulus and SphericalShell children now hold analytic radii to <1e-12 at every level (new gates in test_0836/test_0840); the snapped children solve with custom-P MG in 2 iterations, identical to the GAMG reference. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 46 +++++++++++++++++++ tests/test_0836_nvb_graded_adapt.py | 37 +++++++++++++++ tests/test_0840_nvb_3d_serial_adapt.py | 24 ++++++++++ 3 files changed, 107 insertions(+) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index f0c02025c..78f404655 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -6516,6 +6516,32 @@ def cell_geometry(dm): vol_scaled = numpy.sqrt(numpy.abs(numpy.linalg.det(G))) return X.mean(axis=1), vol_scaled ** (1.0 / dim), cs + # Analytic boundary surfaces to snap each generation onto. New + # boundary vertices are CHORD midpoints, so without snapping the + # boundary geometry of a curved domain stays frozen at base + # resolution no matter how deep the refinement. Per the 2026-07 + # round-3b ruling, EVERY generation snaps (not just the returned + # child): each intermediate level is a valid mesh in its own + # right — extractable for solvers — and the metric marks on true + # geometry. On plane surfaces (boxes) the chord midpoint already + # lies in the plane, so the snap is exactly a no-op and the flat + # confluence gates are untouched. + _snap_surfs = [s for s in dict(self.bounding_surfaces).values() + if getattr(s, "kind", None) in ("radial", "plane") + and not getattr(s, "is_free", False)] + + def snap_level_boundaries(dm): + if not _snap_surfs: + return + from underworld3.meshing.smoothing import _pinned_mask + vec = dm.getCoordinatesLocal() + arr = vec.array.reshape(-1, self.cdim) + for s in _snap_surfs: + mask = _pinned_mask(dm, (s.label,)) + if mask.any(): + arr[mask] = s.restore(arr[mask]) + dm.setCoordinatesLocal(vec) + # Refine from the mesh's CURRENT geometry. Node redistribution # (redistribute_nodes) moves mesh.dm's coordinates while the static # base hierarchy keeps the originals — without this carry, an adapt @@ -6588,6 +6614,7 @@ def cell_geometry(dm): for cidx in marked: lab.setValue(cidx, DM_ADAPT_REFINE) current_dm = _nvbx.refine(d, "adapt") + snap_level_boundaries(current_dm) level_dms.append(current_dm) if verbose: fs, fe = current_dm.getHeightStratum(0) @@ -6628,6 +6655,24 @@ def cell_geometry(dm): marked = [int(cids[j]) for j in sel] markers_per_level.append(marked) nvb.refine(set(marked)) + # Snap INSIDE the engine: its own coordinates feed the next + # generation's marking AND the next midpoints, so snapping + # only the exported DM would leave the engine's geometry on + # the chords. + if _snap_surfs: + _val2surf = {v: s for s in _snap_surfs + for (nm, v) in carry if nm == s.label} + _sv = {} + for _fkey, _val in nvb.facet_label.items(): + _s = _val2surf.get(_val) + if _s is not None: + _sv.setdefault(id(_s), (_s, set()))[1].update(_fkey) + for _s, _vs in _sv.values(): + _idx = numpy.fromiter(_vs, dtype=numpy.int64) + _snapped = _s.restore( + numpy.array([nvb.coords[i] for i in _idx])) + for _k, _i in enumerate(_idx): + nvb.coords[_i] = _snapped[_k] level_dms.append(nvb.to_dm(boundaries=carry, regions=rcarry, comm=self.dm.comm)) if verbose: @@ -6664,6 +6709,7 @@ def cell_geometry(dm): uw.pprint(0, f"[adapt] level {level}: refining {len(cell_ids)} " f"of {ncells} cells") current_dm = custom_mg.sbr_refine(current_dm, cell_ids) + snap_level_boundaries(current_dm) level_dms.append(current_dm) if verbose: diff --git a/tests/test_0836_nvb_graded_adapt.py b/tests/test_0836_nvb_graded_adapt.py index b32a9a963..18b6fbd97 100644 --- a/tests/test_0836_nvb_graded_adapt.py +++ b/tests/test_0836_nvb_graded_adapt.py @@ -300,3 +300,40 @@ def test_nvb_3d_serial_returns_child(): child = base3.adapt(H, max_levels=1, engine="nvb") assert child.parent is base3 assert _ncell(child) > _ncell(base3) + + +def test_curved_boundary_snaps_every_generation(): + """Round-3b ruling (2026-07-24): new boundary vertices on curved + domains snap onto the registered analytic surfaces at EVERY + generation, so each intermediate level is a valid mesh in its own + right and boundary geometry converges with refinement (chords froze + it at base resolution: radius error ~h_base^2/8R).""" + import underworld3 as uw + + mesh = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, + cellSize=0.25, refinement=1, qdegree=2) + + def metric(centroids): + r = np.linalg.norm(np.asarray(centroids)[:, :2], axis=1) + return 1.0 / np.minimum(0.03 + 0.5 * np.abs(1.0 - r), 0.3) ** 2 + + child = mesh.adapt(metric, max_levels=2) + from underworld3.meshing.smoothing import _pinned_mask + + X = np.asarray(child.X.coords) + r = np.linalg.norm(X, axis=1) + for label, R in (("Lower", 0.5), ("Upper", 1.0)): + mask = _pinned_mask(child.dm, (label,)) + assert mask.any() + assert np.abs(r[mask] - R).max() < 1.0e-12, ( + f"{label} boundary not snapped: " + f"max radius error {np.abs(r[mask]-R).max():.2e}") + # ... and the intermediate MG levels are snapped too (the ruling's + # point: every level is a valid mesh) + for i, lvl in enumerate(child._custom_mg_coarse_meshes[-2:]): + dm = lvl.dm if hasattr(lvl, "dm") else lvl + Xl = dm.getCoordinatesLocal().array.reshape(-1, 2) + rl = np.linalg.norm(Xl, axis=1) + on_out = _pinned_mask(dm, ("Upper",)) + if on_out.any(): + assert np.abs(rl[on_out] - 1.0).max() < 1.0e-12 diff --git a/tests/test_0840_nvb_3d_serial_adapt.py b/tests/test_0840_nvb_3d_serial_adapt.py index 4e6cb8024..ac0216a45 100644 --- a/tests/test_0840_nvb_3d_serial_adapt.py +++ b/tests/test_0840_nvb_3d_serial_adapt.py @@ -225,3 +225,27 @@ def test_poisson_fmg_on_3d_nvb_child_matches_gamg(): err = np.linalg.norm(s.Unknowns.u.data[:, 0] - s.Unknowns.u.coords[:, 2]) / ( np.linalg.norm(s.Unknowns.u.coords[:, 2]) + 1e-30) assert err < 1e-8, f"Dirichlet labels wrong on 3D NVB child: err {err}" + + +def test_shell_boundary_snaps_every_generation(): + """3D version of the round-3b snap ruling: an adapted SphericalShell + child holds BOTH spheres to analytic radius at every level.""" + import underworld3 as uw + from underworld3.meshing.smoothing import _pinned_mask + + mesh = uw.meshing.SphericalShell(radiusInner=0.55, radiusOuter=1.0, + cellSize=0.3, refinement=1, qdegree=2) + + def metric(centroids): + r = np.linalg.norm(np.asarray(centroids), axis=1) + return 1.0 / np.minimum(0.08 + 0.6 * np.abs(r - 0.55), 0.35) ** 2 + + child = mesh.adapt(metric, max_levels=1) + X = np.asarray(child.X.coords) + r = np.linalg.norm(X, axis=1) + for label, R in (("Lower", 0.55), ("Upper", 1.0)): + mask = _pinned_mask(child.dm, (label,)) + assert mask.any() + assert np.abs(r[mask] - R).max() < 1.0e-12, ( + f"{label} sphere not snapped: " + f"max radius error {np.abs(r[mask]-R).max():.2e}") From c7735bea8d7c76e11c73e8b5705995dfb0e362bf Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 11:38:33 +1000 Subject: [PATCH 5/7] =?UTF-8?q?fix(adapt):=20review=20follow-ups=20?= =?UTF-8?q?=E2=80=94=20collective=20snap=20mask,=20fold=20guard,=20hardene?= =?UTF-8?q?d=20gates?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The adversarial review of this branch found one real hole on the capstone's target path: in 3D parallel, a boundary-edge midpoint can live on a rank whose only cells containing it are INTERIOR members of the edge's star — that rank sees no labelled face, skips the snap the face-owning rank applies, and one global vertex ends with two coordinates (rank-divergent marking, broken confluence). The snap mask is now reduced over the point SF (ADD as logical-or) before snapping: pre-snap coordinates are rank-consistent and the surface restore is a pure function of them, so a consistent mask gives a consistent snap by construction. Gate: SphericalShell adapt is exactly confluent at np=1/2/4 (12,784 cells) with both spheres at 2e-16 on every rank's chart, ghosts included. Also from review: a loud fold guard after snapping (a cell incident to a snapped vertex may not flip orientation — a base with h ~ R would otherwise tangle silently), the disjoint-surfaces limitation documented at the snap site (sequential restore does not converge to curved-curved junctions), and the intermediate-level radius gate hardened against vacuous passes. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 64 ++++++++++++++++++- tests/test_0836_nvb_graded_adapt.py | 4 ++ 2 files changed, 67 insertions(+), 1 deletion(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 78f404655..5d1d2346e 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -6526,6 +6526,12 @@ def cell_geometry(dm): # geometry. On plane surfaces (boxes) the chord midpoint already # lies in the plane, so the snap is exactly a no-op and the flat # confluence gates are untouched. + # NOTE: surfaces must be disjoint or intersect only where each + # projection is an exact no-op (concentric radial pairs; orthogonal + # planes, whose corner vertices are fixed points of both). Sequential + # restore does not converge to the intersection of two CURVED + # surfaces — junction handling as in boundary_slip is a follow-up + # for the day such a mesh registers surfaces. _snap_surfs = [s for s in dict(self.bounding_surfaces).values() if getattr(s, "kind", None) in ("radial", "plane") and not getattr(s, "is_free", False)] @@ -6534,13 +6540,69 @@ def snap_level_boundaries(dm): if not _snap_surfs: return from underworld3.meshing.smoothing import _pinned_mask + + def sync_mask(mask): + # In 3D a boundary-edge midpoint can live on a rank whose + # only cells containing it are INTERIOR members of the + # edge's star — that rank sees no labelled face and would + # skip the snap the face-owning rank applies, leaving one + # global vertex with two coordinates. Reduce the mask over + # the point SF (ADD ≡ logical-or here) so every rank + # holding the vertex agrees; pre-snap coordinates are + # rank-consistent and restore is a pure function of them, + # so a consistent mask gives a consistent snap. + if uw.mpi.size == 1: + return mask + from underworld3.meshing.smoothing.graph import ( + _build_scalar_dm) + dm_s = _build_scalar_dm(dm) + lvec = dm_s.createLocalVector() + gvec = dm_s.createGlobalVector() + lvec.array[:] = mask.astype(float) + dm_s.localToGlobal(lvec, gvec, addv=True) + dm_s.globalToLocal(gvec, lvec) + out = numpy.asarray(lvec.array) > 0.5 + lvec.destroy(); gvec.destroy(); dm_s.destroy() + return out + vec = dm.getCoordinatesLocal() arr = vec.array.reshape(-1, self.cdim) + pre = arr.copy() + snapped_any = numpy.zeros(arr.shape[0], dtype=bool) for s in _snap_surfs: - mask = _pinned_mask(dm, (s.label,)) + mask = sync_mask(_pinned_mask(dm, (s.label,))) if mask.any(): arr[mask] = s.restore(arr[mask]) + snapped_any |= mask dm.setCoordinatesLocal(vec) + if snapped_any.any(): + # A snap moves boundary vertices by the chord sagitta + # (~h²/8R); on a base coarse enough that h ≈ R this could + # invert a sliver silently. Fail loudly instead: no cell + # incident to a snapped vertex may flip orientation. + cs_, ce_ = dm.getHeightStratum(0) + vS_, vE_ = dm.getDepthStratum(0) + flipped = 0 + for c in range(cs_, ce_): + vs = [p - vS_ for p in dm.getTransitiveClosure(c)[0] + if vS_ <= p < vE_] + if not any(snapped_any[v] for v in vs): + continue + e0 = pre[vs[1:]] - pre[vs[0]] + e1 = arr[vs[1:]] - arr[vs[0]] + if (numpy.sign(numpy.linalg.det(e0)) + != numpy.sign(numpy.linalg.det(e1))): + flipped += 1 + if uw.mpi.size > 1: + from mpi4py import MPI as _MPI + flipped = uw.mpi.comm.allreduce(flipped, op=_MPI.SUM) + if flipped: + raise RuntimeError( + f"adapt: snapping boundary vertices to the analytic " + f"surfaces inverted {flipped} cell(s) — the base mesh " + "is too coarse for the boundary curvature (chord " + "sagitta ~ cell size). Refine the base mesh or adapt " + "without registered bounding surfaces.") # Refine from the mesh's CURRENT geometry. Node redistribution # (redistribute_nodes) moves mesh.dm's coordinates while the static diff --git a/tests/test_0836_nvb_graded_adapt.py b/tests/test_0836_nvb_graded_adapt.py index 18b6fbd97..05f316cd4 100644 --- a/tests/test_0836_nvb_graded_adapt.py +++ b/tests/test_0836_nvb_graded_adapt.py @@ -330,6 +330,7 @@ def metric(centroids): f"max radius error {np.abs(r[mask]-R).max():.2e}") # ... and the intermediate MG levels are snapped too (the ruling's # point: every level is a valid mesh) + checked = 0 for i, lvl in enumerate(child._custom_mg_coarse_meshes[-2:]): dm = lvl.dm if hasattr(lvl, "dm") else lvl Xl = dm.getCoordinatesLocal().array.reshape(-1, 2) @@ -337,3 +338,6 @@ def metric(centroids): on_out = _pinned_mask(dm, ("Upper",)) if on_out.any(): assert np.abs(rl[on_out] - 1.0).max() < 1.0e-12 + checked += 1 + # a label-carry regression must not make this gate silently vacuous + assert checked >= 1, "no intermediate level exposed the Upper label" From 2508596bf948445e2a90394c30de0aa7d2412d7c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 13:16:21 +1000 Subject: [PATCH 6/7] =?UTF-8?q?feat(meshing):=20interior=20interfaces=20?= =?UTF-8?q?=E2=80=94=20adapt()=20snaps=20onto=20them,=20movers=20keep=20th?= =?UTF-8?q?em=20pinned?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Embedded internal surfaces (the Internal circle/sphere of the *InternalBoundary meshes) were topologically preserved by refinement (conforming bisection + label carry) but geometrically chord-frozen, and invisible to the snap because no bounding surface was registered for them. They now register as radial surfaces flagged interior=True: - adapt() includes them in the per-generation snap, so refinement follows the true interface radius (gate: AnnulusInternalBoundary child holds the internal circle to <1e-12); - the slip machinery explicitly skips interior surfaces, so interface nodes stay FULLY pinned under the movers even with slip_surfaces=True — no normal or tangential motion. Interface motion is physics-owned (maintainer ruling, 2026-07-24; tangential slide would be a one-flag change if ever wanted). AnnulusInternalBoundary gains the refinement knob the other constructors already expose (adapt needs the hierarchy); the SphericalShellInternalBoundary registration is in place for when the gmsh blocker clears. Gate locks both behaviours: snap under adapt, bitwise-unmoved interface nodes under the mover. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 9 +++- src/underworld3/meshing/annulus.py | 14 +++++++ src/underworld3/meshing/bounding_surface.py | 20 +++++++-- src/underworld3/meshing/spherical.py | 11 +++++ tests/test_0836_nvb_graded_adapt.py | 41 +++++++++++++++++++ 5 files changed, 90 insertions(+), 5 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 5d1d2346e..4370b9fae 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3136,7 +3136,14 @@ def boundary_slip(self, slip_spec=True, reference_coords=None, # labels (dict ``False``) still slide-without-restore regardless of # kind (handled in ``project`` below). A label with no boundary facets # at all stays unusable → its vertices pin (the safe default). - surf = dict(self.bounding_surfaces) + # INTERIOR surfaces (embedded interfaces, e.g. the Internal circle/ + # sphere of the *InternalBoundary meshes) are never slip-eligible: + # interface motion is physics-owned (deform / free-surface + # machinery), so their nodes fall through to the pinned default — + # exactly as when no surface was registered. Their registration + # exists for adapt()'s refinement snapping. + surf = {lab: s for lab, s in dict(self.bounding_surfaces).items() + if not getattr(s, "interior", False)} unreg = [lab for lab in slip_labels if lab not in surf] if unreg: facets, _opp = _boundary_facets(self, cdim) diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py index a1f47080a..2708560d8 100644 --- a/src/underworld3/meshing/annulus.py +++ b/src/underworld3/meshing/annulus.py @@ -1141,6 +1141,7 @@ def AnnulusInternalBoundary( degree: int = 1, qdegree: int = 2, filename=None, + refinement=None, gmsh_verbosity=0, verbose=False, ): @@ -1403,6 +1404,7 @@ def annulus_internal_mesh_refinement_callback(dm): boundaries=boundaries, boundary_normals=None, coordinate_system_type=CoordinateSystemType.CYLINDRICAL2D, + refinement=refinement, refinement_callback=annulus_internal_mesh_refinement_callback, return_coords_to_bounds=annulus_internal_return_coords_to_bounds, verbose=verbose, @@ -1421,6 +1423,17 @@ class boundary_normals(Enum): x, y = new_mesh.X new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])] + # Radial bounding surfaces: Upper/Lower slip+snap as usual; the embedded + # Internal circle is INTERIOR — adapt() snaps refinement onto its true + # radius, but the movers keep its nodes pinned (interface motion is + # physics-owned, 2026-07 round-3b ruling). + 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, + "Internal": radiusInternal}, + interior=("Internal",)) + return new_mesh @@ -1687,6 +1700,7 @@ def annulus_internal_mesh_refinement_callback(dm): boundaries=boundaries, boundary_normals=None, coordinate_system_type=CoordinateSystemType.CYLINDRICAL2D, + refinement=refinement, refinement_callback=annulus_internal_mesh_refinement_callback, return_coords_to_bounds=annulus_internal_return_coords_to_bounds, verbose=verbose, diff --git a/src/underworld3/meshing/bounding_surface.py b/src/underworld3/meshing/bounding_surface.py index 8b05acf11..d3ec0679f 100644 --- a/src/underworld3/meshing/bounding_surface.py +++ b/src/underworld3/meshing/bounding_surface.py @@ -61,7 +61,8 @@ class BoundingSurface: """ def __init__(self, mesh, label, kind, *, centre=None, radius=None, - point=None, normal=None, reference_facets=None, is_free=False): + point=None, normal=None, reference_facets=None, is_free=False, + interior=False): if kind not in _VALID_KINDS: raise ValueError( f"BoundingSurface kind must be one of {_VALID_KINDS}; got {kind!r}") @@ -69,6 +70,12 @@ def __init__(self, mesh, label, kind, *, centre=None, radius=None, self.label = str(label) self.kind = kind self.is_free = bool(is_free) or kind == "free" + # An INTERIOR surface (an embedded interface such as the Internal + # circle/sphere of the *InternalBoundary meshes) is snapped-to by + # adapt() so refinement follows its true geometry, but is NOT + # slip-eligible: interface motion is physics-owned (deform / + # free-surface machinery), so the movers keep its nodes pinned. + self.interior = bool(interior) 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() @@ -210,15 +217,20 @@ def __repr__(self): # -- constructor-side registration helpers --------------------------------- -def register_radial_surfaces(mesh, centre, label_radius): +def register_radial_surfaces(mesh, centre, label_radius, interior=()): """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).""" + CubedSphere). ``centre`` is the common centre (e.g. the origin). Labels + named in ``interior`` are embedded interfaces: adapt() snaps refinement + onto them, but they stay pinned under the movers (interface motion is + physics-owned).""" centre = np.asarray(centre, dtype=float).ravel() + interior = set(interior) for lab, r in label_radius.items(): mesh.register_tangent_slip_provider( - lab, BoundingSurface(mesh, lab, "radial", centre=centre, radius=r)) + lab, BoundingSurface(mesh, lab, "radial", centre=centre, radius=r, + interior=lab in interior)) def register_plane_surfaces(mesh, label_plane): diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index 36cbd6ffc..496bf222b 100644 --- a/src/underworld3/meshing/spherical.py +++ b/src/underworld3/meshing/spherical.py @@ -780,6 +780,17 @@ def spherical_mesh_refinement_callback(dm): sympy.Matrix([-y, x, 0]), ] + # Radial bounding surfaces: Upper/Lower slip+snap as usual; the embedded + # Internal sphere is INTERIOR — adapt() snaps refinement onto its true + # radius, but the movers keep its nodes pinned (interface motion is + # physics-owned, 2026-07 round-3b ruling). + 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, + "Internal": radiusInternal}, + interior=("Internal",)) + return new_mesh diff --git a/tests/test_0836_nvb_graded_adapt.py b/tests/test_0836_nvb_graded_adapt.py index 05f316cd4..f9e3bebbb 100644 --- a/tests/test_0836_nvb_graded_adapt.py +++ b/tests/test_0836_nvb_graded_adapt.py @@ -341,3 +341,44 @@ def metric(centroids): checked += 1 # a label-carry regression must not make this gate silently vacuous assert checked >= 1, "no intermediate level exposed the Upper label" + + +def test_internal_interface_snaps_under_adapt_stays_pinned_under_mover(): + """Embedded interfaces (the Internal circle of AnnulusInternalBoundary): + adapt() snaps refinement onto the true interface radius (interior + surfaces are snap-eligible), while the mover keeps interface nodes + fully pinned even with slip_surfaces=True — interface motion is + physics-owned (2026-07 round-3b ruling).""" + import sympy + import underworld3 as uw + from underworld3.meshing.smoothing import _pinned_mask + + R_INT = 0.75 + mesh = uw.meshing.AnnulusInternalBoundary( + radiusOuter=1.0, radiusInternal=R_INT, radiusInner=0.5, + cellSize=0.2, refinement=1, qdegree=2) + assert mesh.bounding_surfaces["Internal"].interior + + # adapt: refine toward the interface; new interface vertices snap + def metric(centroids): + r = np.linalg.norm(np.asarray(centroids)[:, :2], axis=1) + return 1.0 / np.minimum(0.04 + 0.6 * np.abs(r - R_INT), 0.25) ** 2 + + child = mesh.adapt(metric, max_levels=1) + Xc = np.asarray(child.X.coords) + mi = _pinned_mask(child.dm, ("Internal",)) + assert mi.any() + ri = np.linalg.norm(Xc[mi], axis=1) + assert np.abs(ri - R_INT).max() < 1.0e-12, ( + f"internal interface not snapped: {np.abs(ri-R_INT).max():.2e}") + + # mover: interface nodes must not move AT ALL (pinned, not sliding) + before = np.asarray(mesh.X.coords).copy() + m0 = _pinned_mask(mesh.dm, ("Internal",)) + x, y = mesh.X + rho = 1 + 6 * sympy.exp(-(((x - 0.2) ** 2 + y**2) / 0.05)) + mesh.redistribute_nodes(rho, slip_surfaces=True, + method_kwargs=dict(n_outer=5)) + after = np.asarray(mesh.X.coords) + assert np.array_equal(after[m0], before[m0]), "interface nodes moved" + assert not np.allclose(after, before) # the rest of the mesh did From 2ca8d96eb97e0935c0ff26c7a6b55e91d96a4d15 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 13:16:47 +1000 Subject: [PATCH 7/7] docs(design): internal-interface resolution recorded (snap + pin ruling) Underworld development team with AI support from Claude Code --- .../design/ADAPTIVITY_3D_SPHERICAL_2026-07.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md b/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md index db8793792..df4442203 100644 --- a/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md +++ b/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md @@ -698,6 +698,16 @@ deliverable. silently falling back to the default preconditioner. Guard relaxed to genuine row/operator mismatches; annulus children now solve with custom MG in 2 iterations. +* *Internal interfaces resolved (maintainer discussion, 2026-07-24):* + refinement always preserved embedded interfaces topologically + (conforming bisection + label carry), but they were chord-frozen and + invisible to the snap. The Internal circle/sphere of the + *InternalBoundary meshes now registers as a radial surface flagged + `interior=True`: adapt() snaps refinement onto the true interface + radius, while the movers keep interface nodes FULLY pinned (no normal + or tangential motion, even with slip_surfaces=True) — interface + motion is physics-owned. Tangential slide along an interface is a + one-flag change if ever wanted. ## Effort and risk, honestly