diff --git a/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md b/docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md index 933a847a..df444220 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,44 @@ 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. +* *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 | round | new-code size | risk | notes | diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index f0c02025..4370b9fa 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) @@ -6516,6 +6523,94 @@ 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. + # 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)] + + 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 = 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 # base hierarchy keeps the originals — without this carry, an adapt @@ -6588,6 +6683,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 +6724,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 +6778,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/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py index a1f47080..2708560d 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 8b05acf1..d3ec0679 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 36cbd6ff..496bf222 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/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 9023f610..ccd9b07c 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 " diff --git a/src/underworld3/utilities/nvb.py b/src/underworld3/utilities/nvb.py index 2860fb3b..b76c9cd8 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_0836_nvb_graded_adapt.py b/tests/test_0836_nvb_graded_adapt.py index b32a9a96..f9e3bebb 100644 --- a/tests/test_0836_nvb_graded_adapt.py +++ b/tests/test_0836_nvb_graded_adapt.py @@ -300,3 +300,85 @@ 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) + 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) + 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 + 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 diff --git a/tests/test_0840_nvb_3d_serial_adapt.py b/tests/test_0840_nvb_3d_serial_adapt.py index 78443859..ac0216a4 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 @@ -187,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}")