diff --git a/.claude/skills/adaptive-meshing/SKILL.md b/.claude/skills/adaptive-meshing/SKILL.md index 16802e28..36102de0 100644 --- a/.claude/skills/adaptive-meshing/SKILL.md +++ b/.claude/skills/adaptive-meshing/SKILL.md @@ -64,8 +64,9 @@ uw.meshing.smooth_mesh_interior(mesh, metric=M, method="mmpde", slip_surfaces=True, skip_threshold=None) # tensor metric: do the skip check yourself ``` -Pitfalls that make it "not work": `method="anisotropic"`/`"ot"` (shred/sliver — use -mmpde); injecting `relax`/`n_outer` (starves mmpde's CG); `strategy=` instead of +Pitfalls that make it "not work": `method="anisotropic"`/`"ot"`/`"spring"`/`"ma"` +(RETIRED 2026-07 — they now raise ValueError; mmpde is the default and only +metric mover); injecting `relax`/`n_outer` (starves mmpde's CG); `strategy=` instead of `refinement=R` (under-grades); a scalar bump for a fault (refines a fat corridor, leaves the centre coarse); signed `Surface.distance.sym` for `d` (bleeds along the line extension — use a direct unsigned distance). Full rationale + the rest of the diff --git a/docs/advanced/mesh-adaptation.md b/docs/advanced/mesh-adaptation.md index 0eb94178..86a27859 100644 --- a/docs/advanced/mesh-adaptation.md +++ b/docs/advanced/mesh-adaptation.md @@ -65,7 +65,7 @@ for each edge vector $\mathbf{e}$. Edges that are too long get subdivided; regio UW3 offers **two complementary** ways to put resolution where it is needed: -| | `mesh.remesh(...)` (this page) | `smooth_mesh_interior(method="anisotropic")` | +| | `mesh.remesh(...)` (this page) | `smooth_mesh_interior(method="mmpde")` | |---|---|---| | Mechanism | **Re-mesh** (MMG): insert/remove/retriangulate | **Redistribute** the existing nodes (move only) | | Node budget | *Changes* — targets an **absolute** edge length `h` | **Fixed** — relative redistribution to a target *density* | @@ -405,9 +405,10 @@ smooth_mesh_interior(mesh, metric=rho, method="mmpde", ```{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 +matrix-free (no PETSc solve — small per-cell dense algebra plus a +parallel `Vec` assembly), provably non-folding, currently 2D +(triangle meshes; the method itself is dimension-general but the +3D discretization is not implemented), and 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 @@ -417,13 +418,13 @@ above; it is promoted to the isotropic tensor `ρ·I`) or a `d×d` 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. +The earlier movers — `"spring"` (volumetric equant-cell smoother), +`"ma"` (isotropic Monge–Ampère), `"ot"` (linear OT-improvement +step) and `"anisotropic"` (decoupled-Winslow tensor smoother) — +were **retired in 2026-07**: `"mmpde"` with a scalar metric +reproduces their isotropic equidistribution, and with a tensor +metric it clusters and aligns where they could not. The retired +spellings now raise a `ValueError` pointing here. Key `mmpde` knobs (via `method_kwargs`): `p` (functional exponent, 1.5–2), `theta` (Huang alignment/equidistribution balance, 1/3), @@ -439,9 +440,6 @@ T| - g_{lo})/(g_{hi}-g_{lo}),0,1\big)$ with $g_{lo},g_{hi}$ the lo/hi percentiles of $|\nabla T|$ — deliberately the same shape as {py:func}`underworld3.adaptivity.metric_from_gradient`, so the *intent* you express is identical whichever family you choose. -The mover then builds a gradient-derived **anisotropic tensor** -metric internally and solves an M-weighted Laplace (Winslow) -coordinate map. ```{important} This is a **gradient** metric: it resolves where the field @@ -455,13 +453,8 @@ on general non-separable features and on cell-alignment / quality (it never produces slivers). ``` -Key knobs (via `method_kwargs`): `aniso_cap` (max cell anisotropy -— the binding stability lever; ≈2 robust, ≳6 folds), `relax` -(damping), `n_outer` (composed damped steps), `linear_solver` -(`"direct"` MUMPS, or `"gamg"` for the parallel-scalable path — -validated bit-parity). The full mathematical derivation (OT / -Monge–Ampère, the metric-tensor / Winslow mover, dynamic field -handling, Nusselt) is in +The full mathematical derivation (including the retired OT / +Monge–Ampère and Winslow movers, kept as an R&D record) is in {doc}`/developer/design/mesh-adaptation-formulation`; operational detail in {doc}`/developer/subsystems/mesh-metric-redistribution`; the dated R&D log in diff --git a/docs/advanced/multigrid-preconditioning.md b/docs/advanced/multigrid-preconditioning.md index fd836958..e85d9c9f 100644 --- a/docs/advanced/multigrid-preconditioning.md +++ b/docs/advanced/multigrid-preconditioning.md @@ -83,8 +83,8 @@ re-discretising on the coarse mesh. `pc_mg_galerkin = both` does this. This is the key reason to prefer FMG when you adapt: -- The coordinate-deforming adaptation movers (Winslow / anisotropic / - `OT_adapt` / `follow_metric`) **preserve mesh topology**. The refinement +- The coordinate-deforming adaptation movers (the MMPDE mover / + `follow_metric`) **preserve mesh topology**. The refinement hierarchy survives them, so geometric multigrid keeps working as the mesh deforms — precisely where GAMG struggles with the resulting anisotropy. - A **true remesh** (a topology change) collapses the hierarchy to a single diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index fe933594..39060345 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,30 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### Retired Interior Movers — MMPDE Is the Mover (July 2026) + +**The superseded fixed-topology interior movers were retired** (maintainer +ruling 2026-07): the spring-equilibrium, Monge–Ampère, OT-improvement-step +and anisotropic-Winslow movers were deleted, together with `mesh.OT_adapt()` +(built on the OT step; closes #346, whose latent MPI deadlock dies with the +spring mover, and #353, whose `strategy=` TypeError dies with the dispatch). + +- `smooth_mesh_interior(method=...)` now defaults to **`"mmpde"`** (was + `"spring"`) — a sanctioned behaviour change: with a scalar metric the + MMPDE mover reproduces the retired movers' isotropic equidistribution + (the isotropic-metric equivalence), and with a tensor metric it clusters + and aligns where they could not. Retired spellings raise a `ValueError` + naming the replacement. +- `follow_metric(...)` (the two-knob adapter) now drives the MMPDE mover; + `mesh.OT_adapt()` raises a `RuntimeError` tombstone pointing at + `follow_metric` / `smooth_mesh_interior` / `mesh.adapt`. +- The graph-Laplacian Jacobi smoother and the Taubin surface-field smoother + (`smooth_surface_field`) are separate, current tools and are unchanged. +- The boundary-facet / boundary-slip primitives shared with surviving code + moved from `meshing/_ot_adapt.py` into `meshing/smoothing/graph.py`; + the style-gate allowlist shrank by the deleted files' entries. + + ### July 2026 Quality Campaign — Audit, Style Charter, Remediation Waves (July 2026) **A systematic post-development-burst quality campaign**: six adversarially diff --git a/docs/developer/design/ma-newton-cofactor-exploration.md b/docs/developer/design/ma-newton-cofactor-exploration.md index 72acfc89..b226984e 100644 --- a/docs/developer/design/ma-newton-cofactor-exploration.md +++ b/docs/developer/design/ma-newton-cofactor-exploration.md @@ -1,5 +1,10 @@ # Monge–Ampère mesh redistribution: Newton/cofactor linearisation +> **RETIRED movers note (2026-07):** the spring / Monge-Ampère / OT-step / +> anisotropic-Winslow interior movers this document discusses were retired +> (superseded by `method="mmpde"`, the default). Kept as the R&D record; +> the code lives in git history before the retirement commit. + > **Status**: exploration (Phase 0), `feature/winslow-mesh-smoother`, > 2026-05-17. Companion to > `docs/developer/subsystems/mesh-metric-redistribution.md` (the diff --git a/docs/developer/design/mesh-adaptation-formulation.md b/docs/developer/design/mesh-adaptation-formulation.md index 40f03cf2..2012b035 100644 --- a/docs/developer/design/mesh-adaptation-formulation.md +++ b/docs/developer/design/mesh-adaptation-formulation.md @@ -1,5 +1,10 @@ # Mesh adaptation by metric-driven node redistribution — mathematical formulation +> **RETIRED movers note (2026-07):** the spring / Monge-Ampère / OT-step / +> anisotropic-Winslow interior movers this document discusses were retired +> (superseded by `method="mmpde"`, the default). Kept as the R&D record; +> the code lives in git history before the retirement commit. + **Status**: Reference (current, 2026-05) — the mathematical formulation for the implemented `smooth_mesh_interior` family (`meshing/smoothing.py`). > **Scope.** This is the self-contained *mathematical* reference for the diff --git a/docs/developer/design/ot-adapt-api-proposal.md b/docs/developer/design/ot-adapt-api-proposal.md index 216cd2d7..104aae35 100644 --- a/docs/developer/design/ot-adapt-api-proposal.md +++ b/docs/developer/design/ot-adapt-api-proposal.md @@ -1,9 +1,14 @@ --- title: "Mesh.OT_adapt() — public API proposal" date: 2026-05-24 -status: proposal +status: retired (2026-07 — OT_adapt removed; superseded by the MMPDE mover) --- +> **RETIRED movers note (2026-07):** the spring / Monge-Ampère / OT-step / +> anisotropic-Winslow interior movers this document discusses were retired +> (superseded by `method="mmpde"`, the default). Kept as the R&D record; +> the code lives in git history before the retirement commit. + # `mesh.OT_adapt()` — public API proposal ## Background diff --git a/docs/developer/subsystems/mesh-metric-redistribution.md b/docs/developer/subsystems/mesh-metric-redistribution.md index e4ce5a99..29390cf4 100644 --- a/docs/developer/subsystems/mesh-metric-redistribution.md +++ b/docs/developer/subsystems/mesh-metric-redistribution.md @@ -7,18 +7,25 @@ which remeshes / changes topology). > **Mathematics:** the full derivations (optimal-transport / > Monge–Ampère, the volumetric spring, the anisotropic -> metric-tensor / Winslow mover, the gradient-metric construction, -> dynamic field handling, and the Nusselt diagnostic) are in +> metric-tensor / Winslow mover — all retired 2026-07 — the +> gradient-metric construction, dynamic field handling, and the +> Nusselt diagnostic) are in > {doc}`/developer/design/mesh-adaptation-formulation`. This page is > the operational guide. +> **Retirement note (2026-07):** the `spring`, `ma`, `ot` and +> `anisotropic` movers described in the historical sections below +> were retired — superseded by `method="mmpde"` (the default), +> whose scalar-metric behaviour reproduces their isotropic +> equidistribution. The retired spellings raise a `ValueError`. +> Sections describing them are kept as the R&D record. + ```python import underworld3 as uw from underworld3.meshing import smooth_mesh_interior -smooth_mesh_interior(mesh, metric=f, method="spring") # fast -smooth_mesh_interior(mesh, metric=f, method="ma") # robust -smooth_mesh_interior(mesh, metric=f, method="anisotropic") # cleanest, aligned +smooth_mesh_interior(mesh, metric=f) # mmpde (default) +smooth_mesh_interior(mesh, metric=M_tensor) # tensor metric: aligned ``` ## When to use it diff --git a/docs/developer/subsystems/meshing.md b/docs/developer/subsystems/meshing.md index 9681fbd3..dfcc4441 100644 --- a/docs/developer/subsystems/meshing.md +++ b/docs/developer/subsystems/meshing.md @@ -64,8 +64,8 @@ This section needs: - [Metric-driven mesh redistribution](mesh-metric-redistribution.md) — topology-preserving node redistribution toward a target - size/density field (`smooth_mesh_interior`; spring & - Monge–Ampère methods). Restores the grading of a deformed + size/density field (`smooth_mesh_interior`; the variational + MMPDE mover). Restores the grading of a deformed adapted mesh or bunches nodes ~2× at a feature; contrast `mesh.adapt()` which remeshes. diff --git a/scripts/deprecated_pattern_allowlist.txt b/scripts/deprecated_pattern_allowlist.txt index 005c1efc..fa55e67f 100644 --- a/scripts/deprecated_pattern_allowlist.txt +++ b/scripts/deprecated_pattern_allowlist.txt @@ -27,9 +27,8 @@ src/underworld3/tests/test_adaptivity_metrics.py:access-context src/underworld3/discretisation/discretisation_mesh.py:mesh-data # --- hedging-name (Charter S3) --------------------------------------------- -# Four `_do_move()` local closures in the mesh-mover machinery. +# `_do_move()` local closures in the mesh-mover machinery (three sites). src/underworld3/discretisation/discretisation_mesh.py:hedging-name -src/underworld3/meshing/_ot_adapt.py:hedging-name src/underworld3/meshing/smoothing/api.py:hedging-name # --- except-pass (Charter S4) ---------------------------------------------- @@ -46,9 +45,6 @@ src/underworld3/function/_function.pyx:except-pass src/underworld3/function/expressions.py:except-pass src/underworld3/function/quantities.py:except-pass src/underworld3/function/unit_conversion.py:except-pass -src/underworld3/meshing/_ot_adapt.py:except-pass -src/underworld3/meshing/smoothing/anisotropic.py:except-pass -src/underworld3/meshing/smoothing/monge_ampere.py:except-pass src/underworld3/meshing/surfaces.py:except-pass src/underworld3/model.py:except-pass src/underworld3/scaling/_scaling.py:except-pass diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 5317095a..d33fc64e 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -495,7 +495,7 @@ def __init__( # Operator on_remesh(ctx) hooks (SemiLagrangian / Lagrangian DDt, # solver-coupled history transfers). Stored as weakrefs so a # forgotten operator does not keep the mesh holding it alive. - # The adapt op (smooth_mesh_interior / OT_adapt / follow_metric) + # The adapt op (smooth_mesh_interior / follow_metric) # fires these after the generic per-variable REMAP pass; see # discretisation/remesh.py. self._remesh_hooks = [] @@ -3067,7 +3067,7 @@ def boundary_slip(self, slip_spec=True, reference_coords=None, """ from underworld3.meshing.smoothing import ( _pinned_mask, _auto_pinned_labels, _owned_vertex_mask) - from underworld3.meshing._ot_adapt import _boundary_facets + from underworld3.meshing.smoothing import _boundary_facets from underworld3.meshing.bounding_surface import BoundingSurface dm = self.dm @@ -3267,7 +3267,7 @@ def __del__(self): def register_remesh_hook(self, op): """Register an operator's ``on_remesh(ctx)`` callback. - Called by the adapt op (``smooth_mesh_interior``, ``OT_adapt``, + Called by the adapt op (``smooth_mesh_interior``, ``follow_metric``) after the generic per-variable REMAP pass. ``op`` must expose an ``on_remesh(ctx)`` method; ``ctx`` is a :class:`~underworld3.discretisation.remesh.RemeshContext` with @@ -3338,7 +3338,7 @@ def _assert_coord_mutation_allowed(self): " Use instead:\n" " • mesh.deform(new_coords, dt=…) — impose an arbitrary " "node displacement (free surface / prescribed motion)\n" - " • mesh.remesh(metric) / mesh.OT_adapt(field)\n" + " • mesh.remesh(metric) / mesh.adapt(metric_field)\n" " • uw.meshing.smooth_mesh_interior(…) / uw.meshing.follow_metric(…)\n" " These route the field + SL/DDt-history transfer " "(remesh_with_field_transfer). For trusted scheme-internal trial " @@ -6046,160 +6046,23 @@ def _increment_mesh_version(self): self._mesh_version += 1 print(f"Mesh version manually incremented to {self._mesh_version}") - @timing.routine_timer_decorator - def OT_adapt( - self, - field, - *, - refinement=3.0, - coarsening="auto", - grad_smoothing_length="auto", - metric_choice="front-following", - fields_to_remap=None, - fields_to_zero=None, - skip_threshold=None, - reference_coords=None, - verbose=False, - ): - r"""Adapt the mesh in place so cell sizes track ``|∇field|``, using - the validated optimal-transport reset pattern. - - Each call resets the mesh to a cached reference (the initial uniform - coordinates), FE-remaps ``field`` onto that clean canvas, builds a - gradient-density metric, runs the OT mover, and FE-remaps the - requested fields onto the adapted positions. Resetting every event - (rather than composing adaptations across time steps) is what keeps - the mover sliver-free over long runs. The "reset" is internal — from - the caller's point of view this just tracks the moving feature. - - Topology is preserved (vertex count, DOF maps, rank partition - unchanged); only coordinates move. Registered solvers are marked for - rebuild via ``_deform_mesh``. - - Reference coordinates - --------------------- - The reset target is snapshotted lazily on the **first** call as - ``self._ot_adapt_reference_coords`` (a copy of the current - ``mesh.X.coords``) and reused thereafter. - - .. warning:: - If the mesh is deformed by something other than ``OT_adapt`` - between calls (e.g. a manual ``mesh._deform_mesh(...)`` or a - resume that loads a *deformed* snapshot), the cached reference no - longer matches the intended pristine state. Use - :meth:`OT_adapt_reset_reference` to re-baseline, or pass an - explicit ``reference_coords`` for a one-off override. - - Parameters - ---------- - field : MeshVariable - Scalar field whose gradient drives refinement (typically ``T``). - Always FE-remapped onto the adapted mesh. - refinement : float, default 3.0 - Cell-size envelope ``h0/refinement`` for the densest cells. - Validated range 1.5–5; 3 is the Nu sweet spot. - coarsening : float or "auto", default "auto" - ``"auto"`` = budget-conserving ``refinement**(1/d)``. - grad_smoothing_length : "auto", None, float, or Pint Quantity, default "auto" - Screened-Poisson de-noising length for ``|∇field|`` before the - metric is built — the most effective sliver lever; without it, - production refinement chases sub-cell gradient noise. - ``"auto"`` (default) ≈ the mesh's uniform cell size (mean edge - length) — the validated setting. ``None`` turns it off. A number - or Pint length sets ``L`` explicitly; **user-supplied lengths are - unit-aware** (non-dimensionalised via the projection), so pass a - Pint quantity (or a non-dimensional number) — ``≈ h0`` is mild, - ``≈ 2·h0`` stronger. - metric_choice : {"front-following", "gradient-uniform"}, default "front-following" - fields_to_remap : list of MeshVariable, optional - Extra fields to FE-remap onto the adapted positions (``field`` is - always remapped). ``None`` ⇒ just ``field``. - fields_to_zero : list of MeshVariable, optional - Fields to zero after the adapt (e.g. ``[V, P]`` for a cold - restart of the flow solve). - skip_threshold : float, optional - If the mesh is already aligned with the metric (misalignment - below this; see :func:`~underworld3.meshing.mesh_metric_mismatch`), - skip the adapt and return ``False``. ``None`` ⇒ always adapt. - reference_coords : array, optional - One-off override of the reset target (does not update the cache). - verbose : bool, default False - - Returns - ------- - bool - ``True`` if the mesh was adapted, ``False`` if the - ``skip_threshold`` check short-circuited it. - - Notes - ----- - Boundary nodes slide tangentially and stay on the boundary for - radial coordinate systems (Annulus / shell), using the projected - boundary normal ``mesh.Gamma_P1``. Cartesian boundaries are pinned - (the vertex-evaluated normal is degenerate there). + def OT_adapt(self, *args, **kwargs): + """RETIRED (2026-07). The optimal-transport reset adapt was + superseded by the variational MMPDE mover; calling this raises + RuntimeError. - Constrained-manifold meshes (``mesh.dim != mesh.cdim``, e.g. a 2D - spherical surface in 3D) are **not supported**: the OT mover would - have to constrain *every* node to the surface, not just boundary - nodes. See ``docs/developer/design/ot-adapt-api-proposal.md``. - - Examples - -------- - >>> mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5, - ... cellSize=1/16, qdegree=3) - >>> T = uw.discretisation.MeshVariable("T", mesh, 1, degree=3) - >>> # ... initialise T ... - >>> mesh.OT_adapt(T, refinement=3.0, fields_to_remap=[T]) - - See Also - -------- - OT_adapt_reset_reference : Re-baseline the reset reference coords. - underworld3.meshing.follow_metric : The single-shot anisotropic mover. - adapt : Topology-changing MMG remeshing (different mechanism). + Use ``uw.meshing.follow_metric(mesh, field, refinement=...)`` (the + two-knob gradient-following adapter) or + ``uw.meshing.smooth_mesh_interior(mesh, metric=..., method="mmpde")`` + for fixed-topology node redistribution; use :meth:`adapt` when more + resolution than a fixed node budget can provide is needed. """ - if self.dim != self.cdim: - raise NotImplementedError( - "OT_adapt is not supported on constrained-manifold meshes " - f"(mesh.dim={self.dim} != mesh.cdim={self.cdim}). The OT " - "mover would need to constrain every node to the surface, " - "not just boundary nodes — see " - "docs/developer/design/ot-adapt-api-proposal.md." - ) - if (not hasattr(self, "_ot_adapt_reference_coords") - or self._ot_adapt_reference_coords is None): - # Lazy snapshot of the reset target on first call. - self._ot_adapt_reference_coords = numpy.asarray( - self.X.coords).copy() - - from underworld3.meshing._ot_adapt import _ot_adapt_step - - return _ot_adapt_step( - self, field, - refinement=refinement, - coarsening=coarsening, - grad_smoothing_length=grad_smoothing_length, - metric_choice=metric_choice, - fields_to_remap=fields_to_remap, - fields_to_zero=fields_to_zero, - skip_threshold=skip_threshold, - reference_coords=reference_coords, - verbose=verbose, - ) - - def OT_adapt_reset_reference(self, coords=None): - r"""Re-baseline the reference coordinates used by :meth:`OT_adapt`. - - ``coords=None`` re-snapshots the current ``mesh.X.coords`` as the new - reset target; passing explicit ``coords`` (e.g. the initial uniform - mesh loaded from a checkpoint) sets those instead. Use on resume, - when the loaded mesh is in a deformed state and the cache would - otherwise lazily initialise from it. - """ - if coords is None: - self._ot_adapt_reference_coords = numpy.asarray( - self.X.coords).copy() - else: - self._ot_adapt_reference_coords = numpy.asarray(coords).copy() + raise RuntimeError( + "mesh.OT_adapt was retired (2026-07): the OT reset adapt was " + "superseded by the variational MMPDE mover. Use " + "uw.meshing.follow_metric(mesh, field, refinement=...) or " + "uw.meshing.smooth_mesh_interior(mesh, metric=..., " + "method='mmpde'); for a topology change use mesh.adapt(...).") @timing.routine_timer_decorator def _wrap_coarse_level(self, dm): diff --git a/src/underworld3/discretisation/remesh.py b/src/underworld3/discretisation/remesh.py index 9dc84c63..ecc23f68 100644 --- a/src/underworld3/discretisation/remesh.py +++ b/src/underworld3/discretisation/remesh.py @@ -18,8 +18,7 @@ displacement, dt, and a bound interpolator. Passed to operator ``on_remesh`` hooks (Phase 2 ALE uses this). * :func:`remesh_with_field_transfer` — the helper used by the adapt - entry points (``smooth_mesh_interior``, ``OT_adapt``, - ``follow_metric``). Takes a closure ``do_move`` that runs the mover + entry points (``smooth_mesh_interior``, ``follow_metric``). Takes a closure ``do_move`` that runs the mover (calling :meth:`Mesh._deform_mesh` repeatedly); the helper handles snapshot + transfer + operator-hook dispatch. @@ -133,9 +132,10 @@ class RemeshContext: ``scratch`` is a free-form dict where adapt ops or hooks publish flags / per-step state. Two keys are in use: - * ``"ale_opt_out"`` — published by ``OT_adapt`` - (:mod:`underworld3.meshing._ot_adapt`, from inside ``do_move`` via - ``mesh._remesh_pending_scratch``) on a reset adapt; consumed by + * ``"ale_opt_out"`` — published (from inside ``do_move`` via + ``mesh._remesh_pending_scratch``) by any adapt whose node motion + is a discrete jump rather than a smooth displacement (the retired + ``OT_adapt`` reset was the original publisher); consumed by the ``DuDt`` ``on_remesh`` hook (:mod:`underworld3.systems.ddt`), which falls back to REMAP for its managed history vars. * ``"v_mesh"`` — the Phase 2 ALE convention: an ALE-style operator @@ -327,15 +327,15 @@ def remesh_with_field_transfer( leave stale values at the partition seams (the documented failure mode that motivated this design). """ - # Re-entrancy guard: composite adapt ops (OT_adapt) wrap the whole - # reset+build+smooth pipeline once at the outer level; inner movers + # Re-entrancy guard: composite adapt ops (follow_metric) wrap the + # whole build+smooth pipeline once at the outer level; inner movers # (smooth_mesh_interior called from inside that pipeline) consult # this flag and skip their own wrap. if getattr(mesh, "_in_remesh_transfer", False): # Nested call: run the mover only. The OUTER wrapper owns # snapshot / transfer / hook dispatch, and it already set # mesh._remesh_pending_scratch, so an inner adapt op can still - # publish flags there (e.g. OT_adapt marking a reset adapt). + # publish flags there (e.g. ale_opt_out on a discrete-jump adapt). # The True return here is unconditional and meaningless — nested # callers cannot use it to tell whether the mesh actually moved; # only the outer call's return value carries that information. @@ -420,8 +420,8 @@ def _remesh_with_field_transfer_impl( # Caller-supplied zero list (e.g. V, P after a topology-preserving # mover, when the flow solve wants a cold start). This is NOT the - # REINIT policy — it is a user knob preserved from the old OT_adapt - # API. REINIT is for *framework-stamped* stateless vars. + # REINIT policy — it is a user knob preserved from the retired + # OT_adapt API. REINIT is for *framework-stamped* stateless vars. if extra_zero: for var in extra_zero: _write_var_data(var, 0.0) @@ -462,8 +462,8 @@ def remap_var_set(mesh, vars_, old_X, new_X, old_data, *, verbose=False): Used by the generic per-variable pass in :func:`remesh_with_field_transfer`, and public so an operator's ``on_remesh`` hook can force-REMAP its CARRY-managed vars on an - adapt that is ALE-incompatible (e.g. an OT_adapt reset, where the - linear ``Δx/dt → v_mesh`` interpretation breaks down). + adapt that is ALE-incompatible (a discrete node-position jump, + where the linear ``Δx/dt → v_mesh`` interpretation breaks down). Contract: on entry the mesh is at ``new_X`` and each var's ``.data`` may hold *either* the original snapshot value (CARRY: diff --git a/src/underworld3/meshing/_ot_adapt.py b/src/underworld3/meshing/_ot_adapt.py deleted file mode 100644 index 515b66d6..00000000 --- a/src/underworld3/meshing/_ot_adapt.py +++ /dev/null @@ -1,423 +0,0 @@ -r"""Optimal-transport mesh adaptation — the validated reset-to-uniform step. - -This module factors the production pattern that was inlined in -``scripts/stagnant_lid_adapt_loop.py`` (the ``ot-reset`` branch) into a -reusable library function. The public entry point is :meth:`Mesh.OT_adapt` -(see ``discretisation/discretisation_mesh.py``); this module holds the -algorithm and the boundary-slip helpers it shares with the OT mover -(``_ot_improvement_step`` in the ``underworld3.meshing.smoothing`` -package, ``smoothing/monge_ampere.py``). - -The algorithm, per adapt event: - -1. Reset the mesh to its reference (IC uniform) coordinates. -2. FE-remap the driving ``field`` onto the reference-mesh DOFs. -3. Build the gradient density metric ``ρ`` on that clean canvas. -4. Run the OT mover from the uniform canvas (``smooth_mesh_interior``, - ``method="ot"``). -5. FE-remap the requested fields onto the adapted positions and zero any - fields flagged for a cold restart. - -The "reset every event" discipline is load-bearing: carrying mesh state -*across* time steps is the broken incremental pattern (slivers lock in). -Composition *within* an adapt is fine. See -``docs/developer/design/ot-adapt-api-proposal.md`` and the -``project_ot_reset_validated`` memory note. - -Boundary slip uses the mesh's **projected boundary normals** -(``mesh.Gamma_P1`` / ``mesh._update_projected_normals``) — the symbolic -``mesh.Gamma`` projected to a P1 vector field and normalised. This is the -general, free-surface-ready normal source: it is re-projected on demand here -because the projected field goes stale every time the mesh deforms. No -per-mesh-class normal code is used. Nodes whose projected normal is -degenerate (box corners, or an occasional unlocatable vertex) are pinned -rather than slipped. -""" - -import numpy as np - -import underworld3 as uw - -# Validated OT-mover constants (2026-05-23/24 investigation). These are -# deliberately *not* exposed on the public OT_adapt signature — they are the -# settled production point, not user dials. -_OT_N_OUTER = 5 -_OT_RELAX = 0.1 -_OT_STEP_FRAC = 0.3 - - -def _is_radial_coords(mesh) -> bool: - """True for coordinate systems with a radial boundary (the snap-back - target is a fixed ``|r|``). Cartesian boundaries are flat — zeroing the - normal displacement keeps nodes on the face, so no snap-back is needed.""" - from underworld3.coordinates import CoordinateSystemType as CT - - return mesh.CoordinateSystem.coordinate_type in ( - CT.CYLINDRICAL2D, - CT.CYLINDRICAL3D, - CT.SPHERICAL, - CT.GEOGRAPHIC, - ) - - -def _auto_grad_smoothing_length(mesh): - """The mesh's characteristic (uniform) cell size — mean edge length, - parallel-safe — returned as a unit-aware length when the mesh carries - coordinate units, else a bare (non-dimensional) float. Used as the - default ``grad_smoothing_length`` so gradient de-noising is on by - default at a scale comparable to the grid (the validated production - setting); ``None`` turns it off.""" - from underworld3.meshing.smoothing import _edge_pairs - - ep = _edge_pairs(mesh.dm) - X = np.asarray(mesh.X.coords) - if ep.shape[0]: - h0 = float(np.linalg.norm( - X[ep[:, 1]] - X[ep[:, 0]], axis=1).mean()) - else: - h0 = 1.0 - if uw.mpi.size > 1: - h0 = uw.mpi.comm.allreduce(h0) / uw.mpi.size - units = getattr(mesh.X, "units", None) - return h0 if units is None else h0 * units - - -def _slip_normals(mesh, boundary_coords: np.ndarray): - """Unit outward normals at ``boundary_coords`` from the projected - boundary-normal field. - - Re-projects ``mesh._projected_normals`` (``mesh.Gamma_P1``) first so the - normals reflect the mesh's *current* coordinates — the projected field is - stale after any deform. Returns ``(normals, valid)`` where ``normals`` is - ``(k, cdim)`` and ``valid`` is a boolean mask; ``valid`` is ``False`` for - nodes with a degenerate (zero / non-finite) normal (e.g. box corners - where opposing face normals cancel, or an occasional unlocatable vertex). - Such nodes should be pinned, not slipped. - """ - cdim = mesh.cdim - n = np.zeros((boundary_coords.shape[0], cdim)) - try: - mesh._update_projected_normals() - n = np.asarray( - uw.function.evaluate(mesh.Gamma_P1, boundary_coords) - ).reshape(-1, cdim) - except Exception: - # Projection unavailable / degenerate on this mesh — fall back to - # all-pinned boundaries (valid stays all-False below). - n = np.zeros((boundary_coords.shape[0], cdim)) - mag = np.linalg.norm(n, axis=1) - valid = np.isfinite(mag) & (mag > 0.5) - out = np.zeros_like(n) - out[valid] = n[valid] / mag[valid, None] - return out, valid - - -def _ot_adapt_step( - mesh, - field, - *, - refinement=3.0, - coarsening="auto", - grad_smoothing_length="auto", - metric_choice="front-following", - fields_to_remap=None, - fields_to_zero=None, - skip_threshold=None, - reference_coords=None, - verbose=False, -) -> bool: - r"""Run one OT-reset adapt event. Returns ``True`` if the mesh moved, - ``False`` if the skip-on-aligned check short-circuited. - - See the module docstring for the algorithm. ``field`` is the scalar - MeshVariable whose gradient drives refinement; it is always FE-remapped - onto the adapted mesh. ``reference_coords`` overrides the reset target - for this call only (defaults to ``mesh._ot_adapt_reference_coords``). - - ``grad_smoothing_length`` de-noises ``|∇field|`` before the metric is - built: ``"auto"`` (default) ≈ the mesh's uniform cell size — the - validated setting that keeps the metric clean at production refinement; - ``None`` turns it off; a number or Pint length sets it explicitly - (user-supplied lengths are unit-aware via the projection's - non-dimensionalisation). - """ - cdim = mesh.cdim - ref_R = float(refinement) - coar = coarsening - if coar != "auto": - coar = float(coar) - # Resolve the gradient de-noising length: "auto" ≈ uniform grid size. - if isinstance(grad_smoothing_length, str): - if grad_smoothing_length.strip().lower() != "auto": - raise ValueError( - "grad_smoothing_length string must be 'auto'; got " - f"{grad_smoothing_length!r}. Pass None (off) or a " - "unit-aware length.") - grad_smoothing_length = _auto_grad_smoothing_length(mesh) - # R for the alignment clamp matches follow_metric: max(refine, coarsen). - coar_val = (ref_R ** (1.0 / cdim)) if coar == "auto" else float(coar) - R_clamp = max(ref_R, coar_val) - - if reference_coords is not None: - ref_X = np.asarray(reference_coords) - else: - ref_X = np.asarray(mesh._ot_adapt_reference_coords) - - # For radial coordinate systems (where boundary slip is used), create the - # projected-normal field up front — before the metric builder / OT mover - # set up any solver DM. Creating that MeshVariable mid-mover would stale - # those DM handles (see project_uw3_smoother_footguns). Cartesian meshes - # pin their boundary (no slip), so no normal field is needed there. - if _is_radial_coords(mesh): - try: - mesh._update_projected_normals() - except Exception: - pass - - # --- skip-on-aligned ------------------------------------------------- - if skip_threshold is not None: - rho_now = uw.meshing.metric_density_from_gradient( - mesh, field, refinement=ref_R, coarsening=coar, - metric_choice=metric_choice, - gradient_smoothing_length=grad_smoothing_length, - degree=1, name="ot_adapt_skip") - mm = uw.meshing.mesh_metric_mismatch( - mesh, rho_now, resolution_ratio=R_clamp) - if mm["misalignment"] < float(skip_threshold): - if verbose: - uw.pprint( - f" OT_adapt: skip — misalignment " - f"{mm['misalignment']:.3f} < {float(skip_threshold):.3f}") - return False - - # Phase-1 remesh redesign: the snapshot/move/transfer dance is now - # owned by the adapt op via remesh_with_field_transfer. The closure - # below performs the reset-to-reference + metric-canvas write + - # OT-mover steps. The helper snapshots `field` (and every other - # REMAP variable on the mesh — including hidden solver history) at - # entry, runs the closure (which may clobber `field` for the metric - # canvas — that write is INTENDED to be discarded by the helper's - # post-move transfer), then performs ONE deform-back / - # global_evaluate / deform-forward pair to bring every REMAP var - # onto the adapted positions. Fields the user previously listed in - # ``fields_to_remap`` are now transferred automatically; the kwarg - # is preserved for API compatibility (vars must be REMAP-policy, - # which is the default — so listing them is a no-op). - from underworld3.discretisation.remesh import ( - remesh_with_field_transfer) - - def _do_move(): - # Phase-2 ALE opt-out: the OT reset-to-reference step is a - # discrete jump in node positions, not a smooth displacement, - # so the linear ``v_mesh = Δx / dt`` interpretation that - # SemiLagrangian.on_remesh uses for ALE is meaningless here. - # Publish a flag so DDt hooks fall back to Phase-1 REMAP for - # this adapt; the mesh's _remesh_pending_scratch dict is the - # pre-fire channel into ctx.scratch. - if hasattr(mesh, "_remesh_pending_scratch"): - scratch = getattr(mesh, "_remesh_pending_scratch", None) - if scratch is not None: - scratch["ale_opt_out"] = True - - old_X_local = np.asarray(mesh.X.coords).copy() - # --- step 1: capture `field` at the reference-mesh DOF positions - mesh._deform_mesh(ref_X) - ref_field_coords = np.asarray(field.coords).copy() - mesh._deform_mesh(old_X_local) - field_at_ref = np.asarray( - uw.function.global_evaluate( - field.sym[0], ref_field_coords)).reshape(-1) - # --- step 2: load the reference (clean) mesh with the remapped field - mesh._deform_mesh(ref_X) - field.data[:, 0] = field_at_ref - # --- step 3: build the gradient metric + run the OT mover - rho = uw.meshing.metric_density_from_gradient( - mesh, field, refinement=ref_R, coarsening=coar, - metric_choice=metric_choice, - gradient_smoothing_length=grad_smoothing_length, - degree=1, name="ot_adapt") - uw.meshing.smooth_mesh_interior( - mesh, metric=rho, method="ot", boundary_slip=True, - method_kwargs=dict(n_outer=_OT_N_OUTER, relax=_OT_RELAX, - step_frac=_OT_STEP_FRAC), - verbose=verbose) - - return remesh_with_field_transfer( - mesh, _do_move, - extra_zero=fields_to_zero, - verbose=verbose, - ) - -# ===== grafted from feature/elliptic-ma: slip helpers for mmpde mover ===== -def _boundary_facets(mesh, cdim): - """Boundary facets + opposite cell-vertex, found from the cell topology. - - For each cell, every facet (edge in 2D, triangle in 3D) is a candidate - boundary facet; one that occurs in **exactly one** cell is on the - boundary. Returns ``(facets, opp)`` where ``facets`` is ``(n_bnd, k)`` - (``k=2`` for 2D edges, ``k=3`` for 3D triangles) and ``opp`` is the - cell vertex opposite each facet — used to orient the facet normal - outward. Returns ``(None, None)`` for non-simplicial meshes. - """ - from underworld3.meshing.smoothing import _tri_cells, _tet_cells - if cdim == 2: - cells = _tri_cells(mesh.dm) - if cells is None: - return None, None - rows = [] - for k in range(3): - v0 = cells[:, k]; v1 = cells[:, (k + 1) % 3] - vopp = cells[:, (k + 2) % 3] - vmin = np.minimum(v0, v1); vmax = np.maximum(v0, v1) - rows.append(np.column_stack([vmin, vmax, vopp])) - e = np.vstack(rows) - idx = np.lexsort((e[:, 1], e[:, 0])) - e = e[idx] - same_prev = np.zeros(len(e), dtype=bool) - same_prev[1:] = ((e[1:, 0] == e[:-1, 0]) - & (e[1:, 1] == e[:-1, 1])) - same_next = np.zeros(len(e), dtype=bool) - same_next[:-1] = same_prev[1:] - bnd_mask = (~same_prev) & (~same_next) - bnd = e[bnd_mask] - return bnd[:, :2], bnd[:, 2] - if cdim == 3: - cells = _tet_cells(mesh.dm) - if cells is None: - return None, None - rows = [] - for k in range(4): - others = [(k + 1) % 4, (k + 2) % 4, (k + 3) % 4] - tri = np.sort(np.column_stack( - [cells[:, others[0]], cells[:, others[1]], - cells[:, others[2]]]), axis=1) - rows.append(np.column_stack([tri, cells[:, k]])) - f = np.vstack(rows) - idx = np.lexsort((f[:, 2], f[:, 1], f[:, 0])) - f = f[idx] - same_prev = np.zeros(len(f), dtype=bool) - same_prev[1:] = ((f[1:, 0] == f[:-1, 0]) - & (f[1:, 1] == f[:-1, 1]) - & (f[1:, 2] == f[:-1, 2])) - same_next = np.zeros(len(f), dtype=bool) - same_next[:-1] = same_prev[1:] - bnd_mask = (~same_prev) & (~same_next) - bnd = f[bnd_mask] - return bnd[:, :3], bnd[:, 3] - return None, None - - -def _all_boundary_labels(mesh): - """Named codim-1 boundary labels of the mesh, skipping the synthetic / - non-geometric ones (``All_Boundaries``, ``Null_Boundary``, and the - Annulus single-point ``Centre`` pseudo-label that hard-aborts PETSc).""" - skip = {"All_Boundaries", "Null_Boundary", "Centre"} - out = [] - try: - names = [b.name for b in mesh.boundaries] - except Exception: - names = [] - for nm in names: - if nm in skip: - continue - out.append(nm) - return tuple(out) - - -def _resolve_slip(mesh, slip_spec): - """Resolve the ``slip_spec`` (the value passed as ``boundary_slip`` / - ``slip_surfaces``) into a tuple of named slip-surface labels, and - pre-touch ``mesh.Gamma_P1`` so the projected-normal field ``_n_proj`` - exists BEFORE any mover builds its solver DM (creating that MeshVariable - mid-mover would stale the DM handle — see project_uw3_smoother_footguns; - the matrix-free ``mmpde`` mover has no such DM but the elliptic / - anisotropic movers do). - - Accepted forms (back-compatible): - * ``True`` / truthy / legacy ``'ring'``,``'box'`` strings → ALL named - codim-1 boundary surfaces slip. - * ``False`` / ``None`` / ``[]`` → no slip (pin all boundaries). - * a label name, or a list of label names → only those surfaces slip. - * a ``dict`` ``{label: snap_bool}`` → those labels slip; ``snap_bool`` - is the per-surface return-to-bounds flag (``False`` = FREE surface, - slip but do not snap back). The dict keys are the slip labels. - - Returns the tuple of slip-surface label names (possibly empty). - """ - if slip_spec is None or slip_spec is False: - return () - if slip_spec is True: - labels = _all_boundary_labels(mesh) - elif isinstance(slip_spec, dict): - labels = tuple(slip_spec.keys()) - elif isinstance(slip_spec, str): - s = slip_spec.strip().lower() - if s in ("ring", "box", "axes", "axis", "true", "on", "1", "all"): - labels = _all_boundary_labels(mesh) - elif s in ("false", "off", "0", "none", ""): - return () - else: - labels = (slip_spec,) # a single explicit label name - else: - # an iterable of label names - labels = tuple(slip_spec) - if labels: - # Pre-create the projected-normal field (footgun-safe; see docstring). - try: - _ = mesh.Gamma_P1 - except Exception: - pass - return labels - - -def _nearest_on_facets_2d(pts, seg): - """Closest point on a set of 2D line segments. ``pts`` (m,2), - ``seg`` (nf,2,2). Returns (m,2) closest points (over all segments).""" - a = seg[:, 0]; b = seg[:, 1] # (nf,2) - ab = b - a - ab2 = np.einsum('fi,fi->f', ab, ab) - ab2 = np.where(ab2 > 1.0e-30, ab2, 1.0) - out = np.empty_like(pts) - for i, p in enumerate(pts): - t = np.clip(((p - a) * ab).sum(axis=1) / ab2, 0.0, 1.0) - proj = a + t[:, None] * ab # (nf,2) - d2 = ((proj - p) ** 2).sum(axis=1) - out[i] = proj[d2.argmin()] - return out - - -def _nearest_on_facets_3d(pts, tri): - """Closest point on a set of 3D triangles. ``pts`` (m,3), - ``tri`` (nf,3,3). Returns (m,3). Per-point loop, vectorised over - triangles via the standard region-based closest-point algorithm.""" - A = tri[:, 0]; B = tri[:, 1]; C = tri[:, 2] - AB = B - A; AC = C - A - out = np.empty_like(pts) - for i, p in enumerate(pts): - AP = p - A - d1 = np.einsum('fi,fi->f', AB, AP) - d2 = np.einsum('fi,fi->f', AC, AP) - BP = p - B - d3 = np.einsum('fi,fi->f', AB, BP) - d4 = np.einsum('fi,fi->f', AC, BP) - CP = p - C - d5 = np.einsum('fi,fi->f', AB, CP) - d6 = np.einsum('fi,fi->f', AC, CP) - va = d3 * d6 - d5 * d4 - vb = d5 * d2 - d1 * d6 - vc = d1 * d4 - d3 * d2 - denom = va + vb + vc - denom = np.where(np.abs(denom) > 1.0e-30, denom, 1.0) - v = vb / denom - w = vc / denom - # interior barycentric point; clamp handles edge/vertex regions well - # enough for a small return-to-bounds correction on convex surfaces. - v = np.clip(v, 0.0, 1.0); w = np.clip(w, 0.0, 1.0) - s = v + w - over = s > 1.0 - v = np.where(over, v / np.where(s > 0, s, 1.0), v) - w = np.where(over, w / np.where(s > 0, s, 1.0), w) - proj = A + v[:, None] * AB + w[:, None] * AC - dd = ((proj - p) ** 2).sum(axis=1) - out[i] = proj[dd.argmin()] - return out diff --git a/src/underworld3/meshing/bounding_surface.py b/src/underworld3/meshing/bounding_surface.py index 99923e47..8b05acf1 100644 --- a/src/underworld3/meshing/bounding_surface.py +++ b/src/underworld3/meshing/bounding_surface.py @@ -12,8 +12,8 @@ ``plane`` restore are analytic; ``facet`` (nearest reference facet) and the ``free`` live-surface follow are follow-ups (the orchestrator pins labels with no analytic surface). The module depends only on the primitive -``_ot_adapt`` helpers that exist on ``development`` (``_slip_normals``), not on -the mover feature branch's unified projector. +boundary helpers in ``meshing.smoothing.graph`` (``_slip_normals``, +``_nearest_on_facets_*``). """ import numpy as np @@ -140,7 +140,7 @@ def normals(self, coords): ``_gamma_p1_at_vertices`` did in git history) and re-solve only for free surfaces. See docs/developer/design/boundary-slip-strategy.md. """ - from underworld3.meshing._ot_adapt import _slip_normals + from underworld3.meshing.smoothing import _slip_normals return _slip_normals(self._mesh, np.ascontiguousarray(coords, dtype=float)) # -- tangent slide ------------------------------------------------------- @@ -182,7 +182,7 @@ def restore(self, coords): if self.kind == "facet": if coords.shape[0] == 0: return coords - from underworld3.meshing._ot_adapt import ( + from underworld3.meshing.smoothing import ( _nearest_on_facets_2d, _nearest_on_facets_3d) if coords.shape[1] == 2: return _nearest_on_facets_2d(coords, self.reference_facets) diff --git a/src/underworld3/meshing/smoothing/__init__.py b/src/underworld3/meshing/smoothing/__init__.py index 2739e546..1071e098 100644 --- a/src/underworld3/meshing/smoothing/__init__.py +++ b/src/underworld3/meshing/smoothing/__init__.py @@ -8,7 +8,7 @@ :func:`metric_density_from_gradient` (metric builder) and :func:`mesh_metric_mismatch` (alignment diagnostic). -The movers, one paragraph each: +The movers: * **Graph-Laplacian Jacobi** (``metric=None``, the no-metric default): each interior vertex is blended toward the mean of its edge @@ -18,47 +18,27 @@ Mat assembled with GLOBAL vertex indices, so partition-boundary rows are complete and results are bit-identical at any rank count. -* **Spring equilibrium** (``method="spring"``, the default metric - path; :func:`_spring_equilibrium_mover`): every edge is an - equal-rest-length spring (shape regulariser, equant cells) plus a - per-cell target-area constraint ``A0 ∝ 1/ρ_tgt`` (the size grading), - minimised to mechanical equilibrium by preconditioned nonlinear CG - (Polak–Ribière⁺) with a fold-rejecting Armijo line search. Fast and - robust — the workhorse for restoring a design grading. Limitation: - edge forces are accumulated rank-locally, so this path is - serial-exact (partition-boundary nodes under-count forces under - MPI; a future PR can assemble forces cross-rank like the Jacobi - adjacency Mat). - -* **Monge–Ampère** (``method="ma"``; :func:`_monge_ampere_mover`): - Benamou–Froese–Oberman convex-branch MA equidistribution with a - variationally-recovered Hessian. Highest-fidelity *isotropic* - refinement, ~60× costlier than the spring; preserved for reference, - not a production choice — see the MA section banner below for the - 2026-05-16 investigation summary (all variants cap at the same - fixed-topology grading). - -* **OT improvement step** (``method="ot"``; - :func:`_ot_improvement_step`): one linear weighted-Poisson - equidistribution-flow step with respect to the *current* mesh, - freely composable. Incomplete (boundary slip is gated to radial - geometries) and deprecated — superseded by ``"mmpde"`` with a - scalar metric. - -* **Anisotropic Winslow** (``method="anisotropic"``; - :func:`_winslow_anisotropic`): the one genuine Winslow smooth here — - an M-weighted Laplace solve of the coordinate map with an - eigen-clamped, gradient-derived metric *tensor*. Reshapes cells - (short across a feature, long along it) and removes slivers; single - primary knob ``resolution_ratio``. Improves alignment / quality, - not grading magnitude. - -* **MMPDE** (``method="mmpde"``; :func:`_mmpde_mover`): Huang–Kamenski - variational moving mesh driven by a full tensor (or scalar) metric — - non-folding by construction, and it genuinely clusters and aligns to - the metric. **Recommended for production adaptive meshing**, though - it is not the API default (``method`` defaults to ``"spring"``). - Currently 2D (triangle meshes) only; parallel-safe. +* **MMPDE** (``method="mmpde"``, the metric-path default; + :func:`_mmpde_mover`): Huang–Kamenski variational moving mesh driven + by a full tensor (or scalar) metric — non-folding by construction, + and it genuinely clusters and aligns to the metric. The production + mover for metric-driven node redistribution. Currently 2D (triangle + meshes) only; parallel-safe. + +The **Taubin surface smoother** (:func:`smooth_surface_field`, in +``graph.py``) is a separate, current tool: it smooths a FIELD on a +boundary/submesh, not the mesh coordinates. + +Retired movers (2026-07 maintainer ruling): the **spring-equilibrium**, +**Monge–Ampère**, **OT-improvement-step** and **anisotropic-Winslow** +interior movers were superseded by ``method="mmpde"`` — with a scalar +(isotropic) metric the MMPDE mover reaches the same equidistributed +grading (the isotropic-metric equivalence), and with a tensor metric it +does what the anisotropic mover could not (genuine clustering and +alignment). Their ``method=`` spellings now raise a ``ValueError`` +naming the retirement. Git history and +``docs/developer/design/anisotropic-mmpde-mover.md`` / +``ma-newton-cofactor-exploration.md`` record the retired algorithms. With a fixed node count no mover exceeds ≈1.3–1.8× deep/near grading (the exact optimal-transport ~10× needs *more nodes* — a topology @@ -66,22 +46,20 @@ ``docs/developer/subsystems/mesh-metric-redistribution.md`` and ``docs/developer/design/anisotropic-mmpde-mover.md``. -Package layout (READ-04 split of the former 4,500-line module): +Package layout (READ-04 split of the former 4,500-line module, pruned +2026-07): -* ``graph.py`` — topology / masks / parallel reductions / shared - mover primitives + :func:`smooth_surface_field` -* ``spring.py`` — the spring-equilibrium mover -* ``monge_ampere.py`` — MA machinery, the OT step, solver wiring -* ``anisotropic.py`` — the tensor (true Winslow) mover +* ``graph.py`` — topology / masks / parallel reductions / boundary + facet + slip primitives + :func:`smooth_surface_field` * ``mmpde.py`` — the variational MMPDE mover * ``metrics.py`` — strategies, metric builder, alignment diagnostic * ``api.py`` — :func:`smooth_mesh_interior` dispatch and :func:`follow_metric` -Every name that was module-level in the old ``smoothing.py`` is -re-exported here, so ``from underworld3.meshing.smoothing import X`` -and ``smoothing.X`` keep working for the private cross-module surface -(``_edge_pairs``, ``_tri_cells``, ``_pinned_mask``, ...). +Every name that is module-level in the submodules is re-exported here, +so ``from underworld3.meshing.smoothing import X`` and ``smoothing.X`` +keep working for the private cross-module surface (``_edge_pairs``, +``_tri_cells``, ``_pinned_mask``, ...). """ import warnings @@ -93,7 +71,6 @@ _build_scalar_dm, _build_adjacency_matrix, _build_local_to_owned_map, - _min_incident_edge, _min_incident_edge_nd, _owned_cell_mask, _tri_cells, @@ -101,13 +78,16 @@ _signed_areas, _edge_pairs, _mean_edge_length, - _cap_step_to_edge_fraction, - _backtracked_move, - _reweight_displacement_radial_tangential, _global_sum, _global_min, _global_max, _global_mean, + _slip_normals, + _boundary_facets, + _all_boundary_labels, + _resolve_slip, + _nearest_on_facets_2d, + _nearest_on_facets_3d, smooth_surface_field, ) from .metrics import ( @@ -116,56 +96,51 @@ mesh_metric_mismatch, metric_density_from_gradient, ) -from .spring import _spring_equilibrium_mover -from .monge_ampere import ( - _use_direct_solver, - _use_iterative_solver, - _solver_wiring, - _warm_start_krylov, - _patch_volumes, - _lumped_vertex_volumes, - _hessian_recovery_class, - _monge_ampere_mover, - _ot_improvement_step, -) -from .anisotropic import _winslow_anisotropic from .mmpde import _mmpde_mover from .api import ( smooth_mesh_interior, _smooth_mesh_interior_bare, follow_metric, + _RETIRED_MOVER_MESSAGE, ) # --------------------------------------------------------------------------- -# One-cycle deprecated aliases (renamed 2026-07, READ-06): the ``_winslow_`` -# prefix was a misnomer on four of the five prefixed movers — only -# ``_winslow_anisotropic`` actually solves a Winslow (M-weighted Laplace) -# coordinate map. The old names are kept for one release cycle because they -# appear in exploratory scripts (scripts/) and external user scripts. +# Retired-mover tombstones (2026-07 retirement ruling). The spring / MA / +# OT / anisotropic-Winslow movers were deleted outright (superseded by +# ``method="mmpde"``); these callables catch scripts that imported the old +# private names — including the READ-06 ``_winslow_*`` aliases, which were +# introduced less than one release cycle before the retirement — and raise +# the same retirement message as the ``method=`` dispatch. # --------------------------------------------------------------------------- -def _deprecated_mover_alias(new_func, old_name): - """Wrap ``new_func`` so calls through ``old_name`` still work but emit - a DeprecationWarning naming the replacement.""" - import functools - - @functools.wraps(new_func) - def _alias(*args, **kwargs): - warnings.warn( - f"{old_name} was renamed to {new_func.__name__} (READ-06: the " - f"mover is not a Winslow smooth); the old name is a one-cycle " - f"deprecated alias.", - DeprecationWarning, stacklevel=2) - return new_func(*args, **kwargs) - - return _alias - - -_winslow_spring = _deprecated_mover_alias( - _spring_equilibrium_mover, "_winslow_spring") -_winslow_elliptic = _deprecated_mover_alias( - _monge_ampere_mover, "_winslow_elliptic") -_winslow_equidistribute = _deprecated_mover_alias( - _ot_improvement_step, "_winslow_equidistribute") -_winslow_mmpde = _deprecated_mover_alias( - _mmpde_mover, "_winslow_mmpde") +def _retired_mover_stub(old_name): + """A callable that raises the retirement ValueError for ``old_name``.""" + + def _stub(*args, **kwargs): + raise ValueError( + f"{old_name} was retired (2026-07): {_RETIRED_MOVER_MESSAGE}") + + _stub.__name__ = old_name + _stub.__doc__ = (f"RETIRED (2026-07). {old_name} was superseded by the " + f"MMPDE mover; calling it raises ValueError.") + return _stub + + +_spring_equilibrium_mover = _retired_mover_stub("_spring_equilibrium_mover") +_monge_ampere_mover = _retired_mover_stub("_monge_ampere_mover") +_ot_improvement_step = _retired_mover_stub("_ot_improvement_step") +_winslow_anisotropic = _retired_mover_stub("_winslow_anisotropic") +_winslow_spring = _retired_mover_stub("_winslow_spring") +_winslow_elliptic = _retired_mover_stub("_winslow_elliptic") +_winslow_equidistribute = _retired_mover_stub("_winslow_equidistribute") + + +def _winslow_mmpde(*args, **kwargs): + """One-cycle deprecated alias for :func:`_mmpde_mover` (READ-06 + rename, 2026-07): identical behaviour plus one DeprecationWarning.""" + warnings.warn( + "_winslow_mmpde was renamed to _mmpde_mover (READ-06: the mover " + "is not a Winslow smooth); the old name is a one-cycle deprecated " + "alias.", + DeprecationWarning, stacklevel=2) + return _mmpde_mover(*args, **kwargs) diff --git a/src/underworld3/meshing/smoothing/anisotropic.py b/src/underworld3/meshing/smoothing/anisotropic.py deleted file mode 100644 index bb3a1065..00000000 --- a/src/underworld3/meshing/smoothing/anisotropic.py +++ /dev/null @@ -1,604 +0,0 @@ -"""The anisotropic tensor-metric mover — the one genuine Winslow -(M-weighted Laplace) coordinate-map smooth in the package. See the -package docstring for the module map. -""" - -import numpy as np - -import underworld3 as uw - -from .graph import (_tri_cells, _signed_areas, _mean_edge_length, - _backtracked_move, - _reweight_displacement_radial_tangential, - _global_sum, _global_min, _global_max, _global_mean) -from .monge_ampere import _solver_wiring, _warm_start_krylov - - -# Cached anisotropic-mover state keyed by (mesh-id, pinned-labels, -# topology, solver, φ-order, slip): the ∇ρ projector, the -# eigen-clamped metric-tensor field D, and the cdim displacement -# Poisson solvers (all sharing the tensor operator _c = D). Rebuilt -# on a topology change (a new key). -_ANISO_CACHE: dict = {} - - -# Per-(mesh,config) running state for the equidistribution -# normaliser's temporal damping: the EMA of ln G carried across -# adaptation events (same key as _ANISO_CACHE). Empty ⇒ first -# event seeds it. Only touched in the resolution_ratio>1 regime. -_GEMA_STATE: dict = {} - - -def _winslow_anisotropic(mesh, metric, pinned_labels, verbose, - n_outer=12, relax=0.2, beta=200.0, - resolution_ratio=1.0, - geom_mean_smoothing=0.25, - aniso_to_base=False, - aniso_cap=2.0, coarsen_cap=1.0, - boundary_slip=False, - linear_solver="direct", phi_degree=2, - move_anisotropy=None, metric_role="M", - outer_tol=1.0e-4, - rest_size_cap_max=None, - rest_size_cap_min=None, - rest_spring_K=1.0, - h0_override=None, - rest_coords_override=None, - metric_refresh_per_iter=False): - r"""Anisotropic metric-tensor mesh redistribution — approach (3). - - The settled scalar equidistribution paths (``_spring_equilibrium_mover``, - ``_monge_ampere_mover``) cannot do coherent *anisotropic* bulk - transport on a fixed topology — a scalar potential is isotropic, - so an annulus radial feature over-collapses one pinned-boundary - sliver layer while the tangential edges sit frozen (see the - project memory + the design doc's angular-OT section). This is - the **tensor** mover: it solves the M-weighted Laplace smooth of - the coordinate map with an *anisotropic* metric tensor, so cells - are reshaped (short across the feature, long along it) and the - slivers / wasted isotropic resolution are removed. - - Construction (verified — ``scripts/ma_metric_tensor_viz.py``): - from a scalar density ``ρ`` (typically Lagrangian - ``f(r0.sym)``), the *projected* gradient ``∇ρ`` (a first - derivative only — UW3-clean) builds, per node, - - .. math:: - - M \;=\; \tfrac1{h_0^2}\!\left[\,I - + \beta\,\hat g\hat g^{\mathsf T} - (|\nabla\rho|/\nabla\rho_{\mathrm{ref}})^2\right], - - eigen-clamped so the spacing ratio ``≤ aniso_cap`` (``≤8:1`` by - default). The eigenframe **auto-aligns to the feature** from the - Cartesian ``∇ρ`` alone — no ``(r,θ)`` frame is specified. - - Mover: solve, per physical coordinate component ``c``, the - displacement form of the M-weighted Laplace (Winslow) map - - .. math:: - - \nabla\!\cdot(D\,\nabla u_c) \;=\; - -\,\nabla\!\cdot(D\,e_c) - \;=\; -\textstyle\sum_j \partial_j D_{jc}, - \qquad u_c = 0 \text{ on the pinned boundary}, - - with ``D = M`` (the eigen-clamped metric). Then - ``ψ_c = x_c + u_c`` is exactly the M-harmonic coordinate map - ``∇·(D∇ψ_c)=0``, ``ψ=x`` on the boundary; the direct Winslow - smoother clusters nodes where ``D`` is large (fine spacing), so - ``D = M`` grades the mesh toward the metric. The two components - share the **same** tensor operator (``_c = D``, the - ``_CofDiff``-style ``DiffusionModel`` pattern) and the - factor-once-reuse direct solver. **Linear** — one solve per - component per outer step, no Picard (much cheaper than the BFO - ``_monge_ampere_mover``). Homogeneous Dirichlet ``u=0`` on the - pinned boundary makes the per-component operator non-singular — - no ``constant_nullspace``, side-stepping the GAMG-pure-Neumann - fragility entirely (``boundary_slip=True`` falls back to the - pure-Neumann + ring-projection treatment of - ``_monge_ampere_mover``). ``n_outer`` composes the map (re-project - ``∇ρ`` / rebuild ``D`` on the moved mesh — the standard MMPDE - outer iteration). Reuses ``_monge_ampere_mover``'s coherent global - signed-area backtrack, ``boundary_slip`` and ``move_anisotropy``. - - .. warning:: - - (3) improves cell **alignment / quality** and removes the - slivers + wasted isotropic resolution; it does **not** beat - the fixed node-count grading cap (≈1.5–1.8× for an 8–20× - density target — that needs ``mesh.adapt``, a topology - change). For a *separable* feature the explicit 1-D OT - (``scripts/ma_analytic_check.py`` / - ``ma_angular_ot_target.py``) is exact and strictly cheaper; - (3) earns its keep on the general **non-separable** case. - Validate with anisotropy-aware diagnostics - (radial/tangential edge split + minA/meanA, *not* the - anisotropy-blind d/n). - - Parameters mirror ``_monge_ampere_mover`` where shared. - - The **decoupled direct** Winslow form (each physical coordinate - M-harmonic, independently) has no Rado–Kneser–Choquet - non-folding guarantee, so its stable regime is bounded by the - metric anisotropy/contrast. Empirically (interior radial - feature, the validation arc) there is a clean Pareto frontier: - - * ``aniso_cap=2``, ``relax≈0.1–0.2`` → minA/meanA ≈ 0.5 (a - near-pristine, valid, feature-aligned mesh — cleaner than the - isotropic MA ≈0.18 / spring ≈0.25 which sliver), modest 2:1 - cell alignment. **The robust default.** - * higher ``aniso_cap`` is only stable with a *gentler* ``relax`` - + more ``n_outer`` (cap 4 needs relax ≈0.05, n_outer ≳25 → - minA ≈0.35, sharper alignment). ``aniso_cap ≳ 6`` folds the - decoupled map regardless — it would need the coupled / inverse - Winslow (the heavy MMPDE, out of this prototype's scope). - - **Single-knob model (`resolution_ratio` R).** The gradient-only - metric ``M ⪰ base·I`` is *refine-only* (keeps only ``∇ρ``, - discards ρ's magnitude ⇒ flat cells pinned at ``h0``, cannot - release nodes, the steepest feature scavenges the budget). The - fix makes the isotropic density a genuinely **equidistributed** - field ``s = base·ρ/G`` (``G`` = geometric mean of ρ on the - near-uniform undeformed D mesh ⇒ ``⟨ln s⟩=ln base``, node budget - centred). Refine (``s>base``) and coarsen (``s cStart else 0 - phi_degree = int(phi_degree) - aux_degree = max(1, phi_degree - 1) - key = (id(mesh), pinned_labels, pEnd - pStart, cEnd - cStart, - cone_size, linear_solver, phi_degree, bool(boundary_slip)) - - cache = _ANISO_CACHE.get(key) - if cache is None: - _wire = _solver_wiring(linear_solver) - - X = mesh.CoordinateSystem.X - # Projected ∇ρ — first derivative only (UW3-clean), the - # same construction verified in ma_metric_tensor_viz. ρ may - # be Lagrangian f(r0.sym): metric.diff(X) then differentiates - # through the frozen r0 field (FE ∂r0/∂x), so ∇ρ is - # re-evaluated on the moved mesh each outer step (MMPDE). - grho = uw.discretisation.MeshVariable( - f"aniso_grho_{id(mesh)}", mesh, - vtype=uw.VarType.VECTOR, degree=aux_degree, - continuous=True) - gproj = uw.systems.Vector_Projection(mesh, grho) - gproj.smoothing = 0.0 - gproj.uw_function = sympy.Matrix( - [metric.diff(X[i]) for i in range(cdim)]).T - _wire(gproj, elliptic=False) - - # Eigen-clamped metric tensor field D (filled numerically - # per outer step). Init to the identity so an unsolved D is - # a harmless isotropic operator. - Df = uw.discretisation.MeshVariable( - f"aniso_D_{id(mesh)}", mesh, - vtype=uw.VarType.TENSOR, degree=aux_degree, - continuous=True) - Df.array[:, 0, 0] = 1.0 - Df.array[:, 1, 1] = 1.0 - Df.array[:, 0, 1] = 0.0 - Df.array[:, 1, 0] = 0.0 - Dsym = Df.sym # 2×2 sympy Matrix (stable) - - class _TensorDiff(uw.constitutive_models.DiffusionModel): - def _build_c_tensor(self): - self._c = Dsym - - # boundary_slip ⇒ pure-Neumann per component (constant - # nullspace, ring-projected in the move — exactly the - # _monge_ampere_mover slip treatment). Default (pinned) ⇒ - # homogeneous Dirichlet u=0 → non-singular, no nullspace. - singular = bool(boundary_slip) - usolvers, ufields = [], [] - for c in range(cdim): - uc = uw.discretisation.MeshVariable( - f"aniso_u{c}_{id(mesh)}", mesh, - vtype=uw.VarType.SCALAR, degree=phi_degree, - continuous=True) - ps = uw.systems.Poisson(mesh, uc) - ps.constitutive_model = _TensorDiff - # f_c = div(column c of D) = Σ_j ∂D_{jc}/∂x_j. UW3 - # SNES_Poisson is F0=-f ⇒ strong form ∇·(D∇u)=-ps.f; - # we want ∇·(D∇u_c) = -div_c ⇒ ps.f = +div_c. (First - # derivative of the projected D field — UW3-legal.) - src = sympy.Integer(0) - for j in range(cdim): - src = src + Dsym[j, c].diff(X[j]) - ps.f = sympy.Matrix([[src]]) - if singular: - ps.constant_nullspace = True - else: - for lbl in pinned_labels: - try: - ps.add_dirichlet_bc(0.0, lbl) - except Exception: - pass - _wire(ps, singular=singular, elliptic=True) - usolvers.append(ps) - ufields.append(uc) - - _ANISO_CACHE[key] = (grho, gproj, Df, usolvers, ufields) - else: - grho, gproj, Df, usolvers, ufields = cache - - zero_init_guess = not _warm_start_krylov(linear_solver) - - # ---- build the eigen-clamped metric tensor field D ONCE ------ - # on the *undeformed* mesh (the design metric), then hold it - # fixed and Lagrangian (the field rides material points through - # _deform_mesh, exactly as _spring_equilibrium_mover computes its - # rest-lengths / A0 once). Re-projecting ∇ρ on the progressively - # distorted mesh inside the outer loop is a positive feedback — - # D blows up on squashed cells → catastrophic over-collapse - # (verified failure mode). With D fixed the outer loop is a - # *stable damped fixed-point iteration* of one linear operator - # toward the M-harmonic map; no feedback. - dm = mesh.dm - # `old0` is the SPRING REST reference — vertices get pulled - # toward these positions when a cell exceeds the size caps. - # If the caller passes `rest_coords_override`, use that - # (typically the truly-undeformed mesh coords captured at - # the first adapt). Falling back to the entry-state of THIS - # call makes the spring "preserve" each successive refined - # state instead of pulling back to undeformed — the third - # leg of the compounding-refinement bug (2026-05-22). - if rest_coords_override is not None: - old0 = np.asarray(rest_coords_override).copy() - else: - old0 = np.asarray(mesh.X.coords).copy() - # h0 = undeformed mean edge length; `h0_override` lets the caller - # supply the value cached at the FIRST adapt (see the - # compounding-refinement note on _mean_edge_length). - if h0_override is not None: - h0 = float(h0_override) - else: - h0 = _mean_edge_length(dm, old0) - # CRITICAL no-op guard: uniform ρ ⇒ ∇ρ ≡ 0, but the L2 - # projection of the zero function leaves ~1e-18 round-off. - # Normalising by that noisy max would make (|∇ρ|/gref)² ~ O(1) - # from pure round-off → a fabricated huge anisotropy and a - # spurious move. Any *real* feature gradient is O(AMP/WIDTH) - # ~ O(1–100); g_eps=1e-9 is ~9 orders above projection noise - # and ~10 below the weakest meaningful feature, so AMP=0 is an - # exact isotropic no-op while AMP>0 is bit-identical to the - # verified ma_metric_tensor_viz construction. - g_eps = 1.0e-9 - base = 1.0 / h0 ** 2 - # Metric-build state shared with the closure below: bound (via - # ``nonlocal``) by ``_build_M_tensor()``, whose pre-loop call is the - # single build site (the per-iteration refresh path rebinds the same - # names when ``metric_refresh_per_iter`` is on). - Dcoords = gvec = gn = gmax = gref = None - - # --- isotropic density: which redistribution model ------------ - # Three regimes, in precedence order: - # - # (1) ``resolution_ratio > 1`` → SINGLE-KNOB EQUIDISTRIBUTION - # (the primary, documented API). The isotropic density is - # ``s = base·ρ/G`` with ``G`` the geometric mean of ρ on - # the (near-uniform, *undeformed*) D mesh, so - # ``⟨ln s⟩ = ln base``: the node budget is centred and - # refine ⇄ coarsen are **complementary by the conservation - # law itself** — there is no coarsening parameter. The - # eigen-clamp ``[base/R², base·R²]`` (cells ∈ - # ``[h0/R, h0·R]``) is a pure safety rail set by the one - # knob ``R``. M-harmonic is scale-invariant, so the - # normalisation *constant* is irrelevant to the realised - # mesh — only ρ's spatial *ratio* and the clamp matter; - # the geometric-mean centring just places the band - # symmetrically so the clamp bites tails, not the bulk. - # - # (2) ``coarsen_cap > 1`` (legacy expert override, not the - # documented API) → the earlier ad-hoc - # ``s = base·cc^(q-1)`` law. Preserved **bit-for-bit** so - # every historical ``a16c*`` result still reproduces. - # - # (3) otherwise → refine-only metric (``s ≡ base``), - # **bit-identical** to the validated historical default. - # ``resolution_ratio = 1`` (the default) lands here ⇒ an - # exact no-op vs. all prior results. - def _build_M_tensor(): - """Compute the metric tensor field Df from the current - metric and mesh state. Mutates Dout-equivalent into Df. - Called once before the iteration loop, and (when - metric_refresh_per_iter=True) also at the start of each - outer iteration to re-query the metric against the - deformed mesh.""" - nonlocal Dcoords, gvec, gn, gmax, gref - Dcoords = np.asarray(Df.coords) # picks up deformed mesh - gproj.solve() - gvec = np.asarray( - uw.function.evaluate(grho.sym, Dcoords) - ).reshape(-1, cdim) - gn = np.linalg.norm(gvec, axis=1) - # Local max first, THEN the (collective) reduction — every rank - # must participate even if it owns no D-mesh points. - gmax = float(gn.max()) if gn.size else 0.0 - gmax = _global_max(gmax) - gref = gmax if gmax > g_eps else 1.0 - # Density branches (same as legacy code path) - if resolution_ratio > 1.0: - R_ = float(resolution_ratio) - rho_v_ = np.asarray( - uw.function.evaluate(metric, Dcoords) - ).reshape(-1) - s_log_ = np.log(np.clip(rho_v_, 1.0e-12, None)) - if uw.mpi.size > 1: - tot = _global_sum(s_log_.sum()) - cnt = _global_sum(s_log_.size) - ln_g_ = tot / max(cnt, 1) - else: - ln_g_ = float(s_log_.mean()) - a_ = float(geom_mean_smoothing) - if 0.0 < a_ < 1.0: - prev = _GEMA_STATE.get(key) - if prev is not None: - ln_g_ = a_ * ln_g_ + (1.0 - a_) * prev - _GEMA_STATE[key] = ln_g_ - iso_ = base * np.exp(s_log_ - ln_g_) - lam_lo_ = base / R_ ** 2 - lam_hi_ = base * R_ ** 2 - aniso_keyed_ = (np.full(Dcoords.shape[0], base) - if aniso_to_base else iso_) - elif coarsen_cap > 1.0: - rho_v_ = np.asarray( - uw.function.evaluate(metric, Dcoords) - ).reshape(-1) - r_lo_ = float(np.percentile(rho_v_, 10.0)) - r_hi_ = float(np.percentile(rho_v_, 90.0)) - r_lo_ = _global_min(r_lo_) - r_hi_ = _global_max(r_hi_) - q_ = np.clip( - (rho_v_ - r_lo_) / max(r_hi_ - r_lo_, 1e-30), - 0.0, 1.0) - iso_ = base * float(coarsen_cap) ** (q_ - 1.0) - lam_lo_ = base / float(coarsen_cap) - lam_hi_ = 1.0 / (h0 / np.sqrt(aniso_cap)) ** 2 - aniso_keyed_ = np.full(Dcoords.shape[0], base) - else: - iso_ = np.full(Dcoords.shape[0], base) - lam_lo_ = base - lam_hi_ = 1.0 / (h0 / np.sqrt(aniso_cap)) ** 2 - aniso_keyed_ = np.full(Dcoords.shape[0], base) - # Assemble M tensor and write to Df - Dout_ = np.empty((Dcoords.shape[0], 2, 2)) - eye2_ = np.eye(2) - for ii in range(Dcoords.shape[0]): - g_ = gvec[ii] - gni_ = gn[ii] - bi_ = iso_[ii] - ai_ = aniso_keyed_[ii] - if gni_ > g_eps and gmax > g_eps: - gh_ = g_ / gni_ - M_ = bi_ * eye2_ + ai_ * beta * (gni_ / gref) ** 2 \ - * np.outer(gh_, gh_) - else: - M_ = bi_ * eye2_ - w_, V_ = np.linalg.eigh(M_) - w_ = np.clip(w_, lam_lo_, lam_hi_) - if metric_role == "Minv": - w_ = 1.0 / w_ - Dout_[ii] = (V_ * w_) @ V_.T - Df.array[:, 0, 0] = Dout_[:, 0, 0] - Df.array[:, 0, 1] = Dout_[:, 0, 1] - Df.array[:, 1, 0] = Dout_[:, 1, 0] - Df.array[:, 1, 1] = Dout_[:, 1, 1] - - # Build D once here, on the undeformed mesh. (This call replaced a - # ~100-line inline duplicate of the closure body — READ-02.) - _build_M_tensor() - - # Pre-compute the undeformed-mesh median cell area, used by the - # backtrack's sliver guard. Captured ONCE before the iteration - # loop so the floor doesn't shrink as cells refine — the same - # absolute floor is enforced throughout. - _tris_for_a0 = _tri_cells(mesh.dm) - if _tris_for_a0 is not None and _tris_for_a0.size: - _a0_undeformed_med = float(np.median(np.abs( - _signed_areas(old0, _tris_for_a0)))) - else: - _a0_undeformed_med = 0.0 - - for outer in range(n_outer): - dm = mesh.dm - pStart, pEnd = dm.getDepthStratum(0) - n_verts = pEnd - pStart - tris = _tri_cells(dm) - old_coords = np.asarray(mesh.X.coords).copy() - _cdim = mesh.cdim - - # If requested, re-query the metric at the deformed - # mesh state and rebuild M tensor. Default off - # preserves the legacy behaviour (M frozen at first - # iteration). Used to isolate whether Eulerian - # re-querying of the metric changes the outcome. - if metric_refresh_per_iter and outer > 0: - _build_M_tensor() - - # Boundary tangential slip via the mesh-owned contract - # (boundary-slip-strategy.md): slip vertices slide tangentially and - # snap back onto their bounding surface (radial ring / plane / facet); - # non-slip, junction, and degenerate-normal vertices pin. Replaces the - # inline per-ring COM radial snap (one node/ring anchored the rotation - # gauge; the signed-area backtrack below still guards against tangle). - is_pinned, _project = mesh.boundary_slip( - boundary_slip, reference_coords=old_coords, - boundary_labels=pinned_labels) - - # D is fixed & Lagrangian (built once, above) — no - # re-projection feedback. The outer loop is a damped - # fixed-point iteration toward the fixed M-harmonic map. - - # --- solve the cdim displacement components ---------------- - disp = np.zeros_like(old_coords) - for c in range(cdim): - usolvers[c].solve(zero_init_guess=zero_init_guess) - disp[:, c] = np.asarray( - uw.function.evaluate(ufields[c].sym, old_coords) - ).reshape(-1) - - # Directional move-weighting (opt-in; default None ⇒ unchanged). - if move_anisotropy is not None and cdim == 2: - disp = _reweight_displacement_radial_tangential( - disp, old_coords, move_anisotropy) - - # --- per-cell Lagrangian rest-size spring ----------------- - # When `rest_size_cap_max` / `rest_size_cap_min` are set, - # add a restoring force to each vertex that pulls it - # toward its rest position (`old0`, captured before the - # mover started) whenever an incident cell's edge would - # overshoot the cap under the proposed move. - # - # We use **max-edge** for the coarsening cap (a cell - # grew in *any* direction beyond `h0·coarsening`) and - # **min-edge** for the refinement cap (a cell shrunk - # in *any* direction below `h0/refinement`). Both - # measures are sliver-aware — they catch anisotropic - # cells that mean-edge wouldn't flag. - # - # Motivation: the metric-mover is a local graph-Laplacian - # — nodes cannot transport across high-gradient ridges, - # so cells *adjacent* to a refinement zone absorb most - # of the freed area while cells topologically isolated - # from the refinement stay near rest size. Without a - # spring, the adjacent cells over-coarsen by ~2× the cap - # and the BL cells over-refine to thin slivers (aspect - # ratios > 10). The spring restores both by literally - # pulling nodes back along the original positions, - # weighted by how much the local cell exceeds the cap. - if (rest_size_cap_max is not None - or rest_size_cap_min is not None): - proposed = old_coords + float(relax) * disp - p = proposed[tris] - 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) - # Sliver-aware per-cell extremes: - max_h = np.maximum(np.maximum(e0, e1), e2) - min_h = np.minimum(np.minimum(e0, e1), e2) - # Per-cell fractional excess vs cap. Both ≥ 0. - # over = max(any edge)/cap_max - 1 (coarsening - # fault: at least one edge too long) - # under = cap_min / min(any edge) - 1 (refinement - # fault: at least one edge too short, i.e. sliver) - if rest_size_cap_max is not None: - over = np.maximum( - max_h / float(rest_size_cap_max) - 1.0, 0.0) - else: - over = np.zeros_like(max_h) - if rest_size_cap_min is not None: - under = np.maximum( - float(rest_size_cap_min) - / np.maximum(min_h, 1.0e-30) - 1.0, 0.0) - else: - under = np.zeros_like(min_h) - # Per-vertex restoring weight ← Σ over incident cells, - # CAPPED AT 1. Without the cap, a vertex incident on - # several violating cells accumulates restore_w > 1 - # and the spring overshoots its rest position - # (`new = old + restore_w · (rest - old)` lands past - # `rest`), pulling two vertices together and creating - # degenerate (near-zero-area) triangles. Capping at 1 - # makes the worst-case per-iteration motion "exactly - # back to rest", never further. - restore_w = np.zeros(old_coords.shape[0]) - cell_w = float(rest_spring_K) * (over + under) - np.add.at(restore_w, tris[:, 0], cell_w) - np.add.at(restore_w, tris[:, 1], cell_w) - np.add.at(restore_w, tris[:, 2], cell_w) - np.minimum(restore_w, 1.0, out=restore_w) - # Add the restoring contribution to disp. (Divide by - # relax so the downstream `step = relax · disp` gives - # the intended fraction restore_w · (rest - current).) - spring_disp = restore_w[:, None] * (old0 - old_coords) - disp = disp + spring_disp / max(float(relax), 1.0e-30) - - # Damped MMPDE step. The *direct* Winslow form (physical - # coords as M-harmonic functions of themselves) has no - # Rado–Kneser–Choquet non-folding guarantee — applied as a - # single elliptic jump it overshoots and the signed-area - # backtrack thrashes into a degenerate sliver. The standard - # remedy is to integrate the mesh PDE as a damped gradient - # flow: under-relax the displacement and compose over - # n_outer steps (the metric is re-projected each step). This - # is the exact analogue of _monge_ampere_mover's picard_relax - # (the BFO path needs ω≈0.4 or its Hessian grows unbounded). - step = float(relax) * disp - - # --- coherent global signed-area backtrack + slip + move -- - # Positive area floor: the flip-only test (`a1min > 0`) misses - # near-degenerate cells with three near-collinear vertices, so - # require min area > 1% of the **undeformed-mesh** median cell - # area (`_a0_undeformed_med`, captured before the iteration - # loop, so the same absolute floor is enforced throughout). A - # refinement of 3 in 2D legitimately shrinks cells by 3²=9× in - # area, so 1% rejects degenerate slivers (1000× smaller) - # without rejecting legitimate refinement. - free = ~is_pinned - new_coords, scale = _backtracked_move( - old_coords, step, free, tris, _project, - area_floor=0.01 * _a0_undeformed_med) - - mesh._deform_mesh(new_coords) - - d = float(np.linalg.norm( - new_coords - old_coords, axis=1).max()) - if uw.mpi.size > 1: - d = _global_sum(d ** 2) ** 0.5 - if verbose: - uw.pprint( - f" anisotropic mover outer {outer+1}/{n_outer}: " - f"h0={h0:.3e} scale={scale:.3f} " - f"max|Δx|={d:.3e}") - if d < outer_tol: - break diff --git a/src/underworld3/meshing/smoothing/api.py b/src/underworld3/meshing/smoothing/api.py index 8ec357ce..a9d3e51d 100644 --- a/src/underworld3/meshing/smoothing/api.py +++ b/src/underworld3/meshing/smoothing/api.py @@ -13,41 +13,50 @@ from .graph import (_auto_pinned_labels, _pinned_mask, _build_adjacency_matrix, _build_local_to_owned_map, - _tri_cells, _signed_areas, _mean_edge_length, + _tri_cells, _signed_areas, _global_sum, _global_min, _global_max) from .metrics import (ADAPT_STRATEGIES, _UNSET, mesh_metric_mismatch, metric_density_from_gradient) -from .spring import _spring_equilibrium_mover -from .monge_ampere import _monge_ampere_mover, _ot_improvement_step -from .anisotropic import _winslow_anisotropic from .mmpde import _mmpde_mover +# The retired ``method=`` spellings (2026-07 maintainer ruling: the +# spring / Monge-Ampère / OT-step / anisotropic-Winslow interior movers +# are superseded by the MMPDE mover — a scalar metric reproduces their +# isotropic equidistribution, a tensor metric does what they could not). +_RETIRED_METHODS = { + "spring": "spring-equilibrium", + "ma": "Monge-Ampère", "monge-ampere": "Monge-Ampère", + "monge_ampere": "Monge-Ampère", + "ot": "OT-improvement-step", "equidistribute": "OT-improvement-step", + "improve": "OT-improvement-step", + "anisotropic": "anisotropic-Winslow", "aniso": "anisotropic-Winslow", + "tensor": "anisotropic-Winslow", +} + +_RETIRED_MOVER_MESSAGE = ( + "the spring-equilibrium, Monge-Ampère, OT-step and anisotropic-Winslow " + "interior movers were retired 2026-07, superseded by the variational " + "MMPDE mover. Use smooth_mesh_interior(mesh, metric=..., " + "method='mmpde') — a scalar metric gives the same isotropic " + "equidistributed grading, a tensor metric adds genuine anisotropic " + "clustering — or uw.meshing.follow_metric(...) for the two-knob " + "refinement/coarsening form. For more resolution than a fixed node " + "budget can grade, change topology with mesh.adapt(...) instead.") + + # Cached adjacency keyed by (mesh-id, pinned-label-tuple, topology). # Rebuilt automatically when the mesh topology changes. _ADJ_CACHE: dict = {} -# Cache of the **original** (undeformed) state per mesh, -# captured the first time follow_metric is called on that mesh: -# h0 — mean edge length -# rest_coords — vertex positions (the spring's pull-back target) -# Subsequent calls reuse these references instead of measuring the -# (already-refined) current mesh, otherwise the spring's reference -# state shrinks at every adapt and the refinement compounds, -# crashing the CFL-bound dt by 2× per adapt step. -# Keyed by id(mesh). -_FOLLOW_METRIC_H0_CACHE: dict = {} -_FOLLOW_METRIC_REST_CACHE: dict = {} - - def smooth_mesh_interior( mesh, pinned_labels: Optional[Sequence[str]] = None, n_iters: int = 5, alpha: float = 0.5, metric=None, - method: str = "spring", + method: str = "mmpde", boundary_slip: bool = False, method_kwargs: Optional[dict] = None, verbose: bool = False, @@ -70,19 +79,17 @@ def smooth_mesh_interior( over ``n_iters`` sweeps. Equalises connectivity → equant cells. - **With a ``metric``** — an elastic-spring network relaxed to - equilibrium. Every edge is a linear spring with rest length - ``∝ ρ_tgt^{-1/d}`` (``ρ_tgt = metric``), scaled so the mean rest - length equals the current mean edge length (overall scale - preserved — pure redistribution). Damped Jacobi force iteration - relaxes interior nodes to force balance, with a coherent global - signed-area backtrack guaranteeing no cell inverts. The rest - length is an *absolute* target, so the mesh genuinely grades - toward spacing ``∝ ρ_tgt^{-1/d}`` (a regime the weighted - Laplacian / Jacobi cannot reach). ``n_iters`` and ``alpha`` are - ignored on this path (it has its own internal sweep budget). A + **With a ``metric``** — the variational MMPDE mover + (Huang–Kamenski; ``method="mmpde"``, the only metric-grading + method since the 2026-07 mover retirement). A scalar metric + ``ρ_tgt`` equidistributes cell size toward spacing + ``∝ ρ_tgt^{-1/d}``; a full ``d×d`` tensor metric additionally + clusters and *aligns* cells to the tensor (short across a + feature, long along it). The map is non-folding by construction. + ``n_iters`` and ``alpha`` are ignored on this path (the mover + has its own outer-iteration budget and convergence test). A Lagrangian density (``f(r0.sym)`` peaked at the original outer - radius) keeps the rest lengths fixed per material point, so the + radius) keeps the target fixed per material point, so the *design* boundary-layer grading is restored even after free-surface deformation. @@ -102,7 +109,7 @@ def smooth_mesh_interior( n_iters : int, default 5 Number of Jacobi sweeps. 5-10 is typical for surface- deformation cleanup. **Ignored when ``metric`` is given** - (the spring path has its own internal sweep budget). + (the metric path has its own outer-iteration budget). alpha : float, default 0.5 Under-relaxation in ``(0, 1]`` for the Jacobi path. 1.0 is pure Jacobi; smaller is more damped. **Ignored when @@ -117,118 +124,45 @@ def smooth_mesh_interior( positive and finite. ``None`` (default) ⇒ the graph-Laplacian Jacobi path, unchanged behaviour bit-for-bit. - method : {"spring", "ma", "anisotropic", "mmpde"}, default "spring" + method : str, default "mmpde" Metric-grading solver (ignored when ``metric is None``). - ``"mmpde"`` is the recommended production mover for adaptive - meshing; ``"ot"`` is accepted but deprecated (incomplete — - prefer ``"mmpde"`` with a scalar metric): - - * ``"spring"`` — *volumetric* elastic-spring equilibrium: - equal edge springs (shape regulariser, equant cells, no - slivers) + a per-cell area constraint - ``A0 ∝ 1/ρ_tgt`` (the size grading), minimised by - preconditioned nonlinear CG. **Fast** (~0.3 s on a - res-16 Annulus), robust, scales with the metric - amplitude; slightly anisotropic at sharp interior - features. - * ``"ma"`` — Benamou–Froese–Oberman convex-branch - **Monge–Ampère** equidistribution. Highest-fidelity - *isotropic* refinement and robust to the boundary - treatment, but ~60× costlier than the spring. - * ``"anisotropic"`` — **tensor** metric mover: an - M-weighted Laplace (Winslow) smooth of the coordinate - map with an eigen-clamped, gradient-derived *anisotropic* - metric tensor. Reshapes cells (short across a feature, - long along it) and removes the slivers / wasted isotropic - resolution the scalar paths leave near a boundary-peaked - feature. Linear (one solve/component/step — cheaper than - ``"ma"``). It improves cell **alignment / quality**, not - the grading magnitude (see the cap note below); for a - *separable* feature the explicit 1-D OT is exact and - cheaper — ``"anisotropic"`` earns its keep on the general - non-separable case. - * ``"mmpde"`` — variational moving-mesh (Huang–Kamenski - MMPDE) with a full tensor (or scalar) metric; the - recommended production mover for adaptive meshing. - **Currently 2D-only** (triangle meshes) — a 3D mesh - raises ``NotImplementedError``. - * ``"ot"`` (deprecated) — one linear OT-improvement step, - composable; boundary slip is gated to radial geometries. - Kept for the internal ``mesh.OT_adapt`` reset path; new - code should use ``"mmpde"``. - - With a fixed node count neither can exceed ≈1.3–1.8× + ``"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``. + + The historical spellings ``"spring"``, ``"ma"``, ``"ot"`` + and ``"anisotropic"`` name interior movers retired in + 2026-07 (superseded by ``"mmpde"``) and raise a + ``ValueError`` describing the replacement. + + With a fixed node count the mover cannot exceed ≈1.3–1.8× deep/near grading (the optimal-transport ≈10× needs *more nodes* — a topology change, not this smoother). See ``docs/developer/subsystems/mesh-metric-redistribution.md``. boundary_slip : bool, default False Let boundary nodes slide tangentially along their boundary (snapped back to the boundary each step — they cannot leave - it; serial circular/spherical boundaries only). Strongly - helps the spring (+~10 % grading, faster); near-no-op for - ``ma`` (its natural Neumann BC already handles the - boundary). Off by default — for a free surface the boundary - is the moving surface, so sliding interacts with the - free-surface coupling; enable per use-context. + it). Off by default — for a free surface the boundary is + the moving surface, so sliding interacts with the + free-surface coupling; enable per use-context. Prefer the + richer ``slip_surfaces`` spelling in new code. method_kwargs : dict, optional - Extra tuning forwarded to the chosen metric solver (ignored - when ``metric is None``). Keeps the shared signature clean - while exposing the per-method knobs. For - ``method="anisotropic"`` there is **one primary knob**: - - * ``resolution_ratio`` (``R``, default **1.0 = exact - no-op**) — *the* tuneable. Cells may refine to ``h0/R`` - and coarsen to ``h0·R``; the refine ⇄ coarsen split is - **not a parameter** — the isotropic density is - equidistribution-normalised (``s = base·ρ/G``, ``G`` the - geometric mean of ρ), so flat regions release exactly the - budget the fronts consume, *complementary by the - conservation law itself*. The eigen-clamp - ``[h0/R, h0·R]`` is just a safety rail. ``R=1`` ⇒ - bit-identical to the refine-only historical default (an - exact no-op vs. every prior result). ``R≈2`` is the - validated production point (clean mesh through a full - convection lifecycle, ``minA/meanA``≈0.2, genuine - plume-reaching de-resolution, settled physics intact). - One number; complementary coarsening is automatic. - * ``geom_mean_smoothing`` (``a``, default 0.25) — - *internal* temporal damping of the equidistribution - normaliser ``G`` (not a grading knob; only acts when - ``R>1``). ``G`` is recomputed from the instantaneous - field every adaptation event; in a violent transient - that lurches the whole ``ρ/G`` distribution across the - fixed clamp band → clamp-saturation → the mesh visibly - "wobbles". An EMA in log space - (``lnG ← a·lnG_now+(1−a)·lnG_prev``) keeps the band - centred: ``a=1`` ⇒ no damping (instantaneous, the - original wobbly behaviour); ``a≈0.25`` ⇒ strong damping - of the startup over-reaction + steady-state contrast - pulse. It smooths **only the one global intensity - scalar** — the spatial ρ(x) pattern still tracks the - current field every event, so the API stays single-knob - (``R``); ``a`` carries one internal scalar across events. - * ``relax`` (0.2) / ``n_outer`` (12) — damped-MMPDE - under-relaxation + composed steps (early-exit - ``outer_tol``). ``linear_solver`` (``"direct"`` | MUMPS | - ``"gamg"``, bit-parity, parallel-scalable). ``beta`` - (200) — anisotropic-bump saturation. ``move_anisotropy`` - — optional radial/tangential move reweight. - * **Expert overrides (not the documented API; only honoured - when ``resolution_ratio≤1``):** ``aniso_cap`` (2.0) and - ``coarsen_cap`` (1.0) are the legacy two-knob clamp - (``h_min=h0/√aniso_cap``, ``h_max=h0·√coarsen_cap``, - ad-hoc ``s=base·cc^(q-1)``). Retained **bit-for-bit** so - historical scripts reproduce; superseded by - ``resolution_ratio``. + Extra tuning forwarded to the mover (ignored when ``metric + is None``). The MMPDE tunables are documented on + :func:`~underworld3.meshing.smoothing.mmpde._mmpde_mover` + (``n_outer``, ``step_frac``, ``accel``, ``metric_eval``, + tolerances, ...); unknown keywords are warned about and + ignored. Example:: smooth_mesh_interior( - mesh, metric=rho, method="anisotropic", - method_kwargs=dict(resolution_ratio=2.0, - relax=0.05, n_outer=25)) + mesh, metric=rho, method="mmpde", + method_kwargs=dict(n_outer=50, accel="cg")) verbose : bool, default False - Print per-sweep (Jacobi) or periodic (spring/MA) progress. + Print per-sweep (Jacobi) or periodic (mover) progress. skip_threshold : float, optional If set, evaluate the *misalignment* between current mesh cell density and the metric (via @@ -260,12 +194,10 @@ def smooth_mesh_interior( Results are bit-identical (to a single ULP) between serial and parallel runs at any rank count. - **Spring path**: serial-exact. Edge forces are accumulated over - locally-visible edges only, so rank-partition-boundary nodes - under-count their incident forces in parallel (a future PR can - assemble the edge forces cross-rank like the Jacobi adjacency - Mat). The edge list and per-node degree are cached against the - topology key and rebuilt only on a topology change. + **MMPDE path**: parallel-safe — the velocity assembly reduces + cross-rank via the coordinate DM and all accept/backtrack + decisions are collective, so serial and parallel runs take the + same steps (see the mover's docstring). **Topology preservation**: vertex IDs, DOF mappings, and the rank partition are unchanged. Only coordinates move. Anything @@ -326,8 +258,8 @@ def smooth_mesh_interior( # snapshotted, the mover runs, and a single deform-back / # global_evaluate / deform-forward pair carries every variable onto # the adapted node positions. Re-entrancy guard - # ``mesh._in_remesh_transfer`` lets composite adapts (OT_adapt) wrap - # the whole reset+build+smooth dance once at the outer level and + # ``mesh._in_remesh_transfer`` lets composite adapts (follow_metric) + # wrap the whole build+smooth dance once at the outer level and # have this inner call skip its own wrap. if not getattr(mesh, "_in_remesh_transfer", False): from underworld3.discretisation.remesh import ( @@ -371,7 +303,7 @@ def _smooth_mesh_interior_bare( n_iters: int = 5, alpha: float = 0.5, metric=None, - method: str = "spring", + method: str = "mmpde", boundary_slip: bool = False, method_kwargs: Optional[dict] = None, verbose: bool = False, @@ -381,8 +313,8 @@ def _smooth_mesh_interior_bare( """Internal mover dispatch — no transfer, no helper wrap. Identical to the body of :func:`smooth_mesh_interior` minus the - Phase-1 transfer wrap. Composite adapt ops (``_ot_adapt_step``, - ``follow_metric``) own the wrap at their level and call this bare + Phase-1 transfer wrap. Composite adapt ops (``follow_metric``) + own the wrap at their level and call this bare form to avoid nesting the snapshot/restore dance. End-users should keep using :func:`smooth_mesh_interior`. """ @@ -411,14 +343,11 @@ def _smooth_mesh_interior_bare( method_kwargs = {} else: method_kwargs = dict(method_kwargs) - # TODO(BUG): this injection is unconditional, but only the - # 'anisotropic' and 'mmpde' movers accept resolution_ratio — - # strategy= combined with the default method='spring' (or - # 'ma'/'ot') raises TypeError at the mover call. Pre-existing - # at the Wave D base (verified by signature-binding probe); - # spring/MA/OT are retired in favour of 'mmpde' (see #346 - # context), so the fix is to inject only for movers that - # accept it (or route strategy= to a surviving mover). + # #353: strategy= now routes to the mmpde mover (the default + # method), which accepts resolution_ratio (advisory — MMPDE's + # clustering intensity comes from the metric itself). The + # retired movers that could not bind it raise their retirement + # ValueError before this kwarg reaches any call. method_kwargs.setdefault( "resolution_ratio", _s["resolution_ratio"]) if skip_threshold is _UNSET: @@ -432,8 +361,8 @@ def _smooth_mesh_interior_bare( # evaluate + a few NumPy reductions) — avoids a redundant # mover call when the mesh hasn't drifted from its target. # Mismatch is measured against the R-clamped achievable - # target (when the anisotropic mover's resolution_ratio is - # given), so a perfectly-adapted mesh measures ~0. + # target (when a resolution_ratio is given in method_kwargs), + # so a perfectly-adapted mesh measures ~0. if skip_threshold is not None: _R = mk.get("resolution_ratio", None) mm = mesh_metric_mismatch( @@ -465,48 +394,20 @@ def _smooth_mesh_interior_bare( f"threshold {float(skip_threshold):.3f}; " f"alignment r={mm['alignment']:.3f})", flush=True) - if method == "spring": - _spring_equilibrium_mover(mesh, metric, pinned_labels, verbose, - boundary_slip=boundary_slip, **mk) - elif method in ("ma", "monge-ampere", "monge_ampere"): - _monge_ampere_mover(mesh, metric, pinned_labels, verbose, - boundary_slip=boundary_slip, **mk) - elif method in ("ot", "equidistribute", "improve"): - # The OT / equidistribution mover is incomplete — e.g. its boundary - # slip is gated to radial geometries (box boundaries are pinned, not - # slid; see boundary-slip-strategy.md) — and is expected to be - # superseded by ``method='mmpde'`` with a scalar metric. This fires - # for every OT use, including the internal ``mesh.OT_adapt`` reset - # path. (Python shows a given DeprecationWarning once per location.) - warnings.warn( - "smooth_mesh_interior(method='ot'/'equidistribute'/'improve') " - "is an incomplete mesh mover (boundary slip is gated to radial " - "geometries) and is expected to be superseded by " - "method='mmpde' with a scalar metric. Prefer 'mmpde' for " - "production adaptive meshing.", - DeprecationWarning, stacklevel=2) - _ot_improvement_step(mesh, metric, pinned_labels, - verbose, - boundary_slip=boundary_slip, - **mk) - elif method in ("anisotropic", "aniso", "tensor"): - _winslow_anisotropic(mesh, metric, pinned_labels, - verbose, - boundary_slip=boundary_slip, **mk) - elif method in ("mmpde", "variational"): + if method in ("mmpde", "variational"): _mmpde_mover(mesh, metric, pinned_labels, verbose, boundary_slip=boundary_slip, **mk) + elif method in _RETIRED_METHODS: + raise ValueError( + f"smooth_mesh_interior: method={method!r} " + f"({_RETIRED_METHODS[method]}) was retired: " + f"{_RETIRED_MOVER_MESSAGE}") else: raise ValueError( f"smooth_mesh_interior: unknown method {method!r}; " - f"use 'spring' (default, fast volumetric), " - f"'ma' (Monge–Ampère, isotropic, ~60× costlier), " - f"'anisotropic' (tensor metric — reshapes cells / " - f"removes slivers; does not beat the node-count " - f"cap), 'mmpde' (variational moving mesh — the " - f"recommended production mover) or " - f"'ot' / 'equidistribute' (deprecated linear " - f"OT-improvement step).") + f"use 'mmpde' (the variational moving-mesh mover — " + f"the only metric-grading method; scalar or tensor " + f"metric).") return dm = mesh.dm @@ -621,40 +522,34 @@ def follow_metric( r"""Move the mesh's interior nodes so cell sizes follow a target derived from ``|∇field|``. - Two-knob, cell-size-envelope API for the anisotropic node mover. - The user specifies how *fine* the densest cells can get and - (optionally) how *coarse* the sparsest can get; the function - derives the metric density and invokes the mover. + Two-knob, cell-size-envelope API for the metric-driven node + mover (the variational MMPDE mover since the 2026-07 mover + retirement). The user specifies how *fine* the densest cells + can get and (optionally) how *coarse* the sparsest can get; the + function derives the metric density and invokes the mover. Cell-size envelope (approximate) -------------------------------- - The mover's eigenvalue → cell-size map is - :math:`h = h_0/\sqrt{\hat\rho}` (after the mover's - geometric-mean normalisation :math:`\hat\rho = \rho/G`), so - asking for the envelope + The metric density is built so that equidistribution maps it to + cell size as :math:`h = h_0/\sqrt{\rho}` (with + :math:`\mathrm{geomean}(\rho) = 1` by construction), so asking + for the envelope .. math:: h \;\in\; \bigl[\, h_0/\text{refinement},\; h_0\cdot\text{coarsening} \,\bigr] - corresponds to :math:`\hat\rho \in [1/\text{coarsening}^2, + corresponds to :math:`\rho \in [1/\text{coarsening}^2, \text{refinement}^2]` — note this is **dimension- - independent** (the eigenvalue λ has units of 1/length²). - - Validation on a sharp-tanh annulus test problem shows: - - * **Refinement side:** achieved :math:`h_\min` within ~5-10% - of :math:`h_0/\text{refinement}` for refinement ∈ [1.5, 3]. - * **Coarsening side:** achieved :math:`h_\max` typically - ~2× the requested :math:`h_0\cdot\text{coarsening}`. The - mover's anisotropic cells and iterative deformation map - together don't honour the eigenvalue clamp on a per-cell - basis as tightly as the refinement side. This is a known - feature of the underlying mover, not of the new API. + independent**. - The :func:`mesh_metric_mismatch` diagnostic is the right tool + The achieved envelope is *approximate*: the mover + equidistributes toward the metric under a fixed node budget, it + does not clamp per-cell sizes. Expect the refinement side to + track the request more tightly than the coarsening side. The + :func:`mesh_metric_mismatch` diagnostic is the right tool for measuring how close the achieved mesh is to the requested metric in practice. @@ -673,9 +568,8 @@ def follow_metric( \text{coarsening})} . This break point makes :math:`\mathrm{geomean}(\rho) = 1` - by construction, so the mover's :math:`G`-normalisation - leaves :math:`\rho` unshifted and the eigenvalue clamps land - on the desired envelope. Concretely: + by construction, so equidistributing :math:`\rho` lands the + cell sizes on the desired envelope. Concretely: * "front-following" (default) — log-:math:`\rho` is linear in percentile rank on each side of :math:`p^{\ast}`. Every @@ -692,7 +586,7 @@ def follow_metric( :math:`\rho = \sqrt{1 + (A\,|\nabla\text{field}|/g_{hi})^2}`, clipped to the envelope. Grades continuously from :math:`\rho = 1` in flat regions (no clip kink), giving - cleaner OT / Monge–Ampère meshes. + the smoothest equidistributed meshes. Auto coarsening (the budget-conserving default) ----------------------------------------------- @@ -747,7 +641,7 @@ def follow_metric( per-cell field change everywhere (best for advection- diffusion accuracy). ``"arc-length"`` is a smooth arc-length monitor — grades continuously from flat - regions with no clip kink (cleaner OT / Monge–Ampère + regions with no clip kink (the smoothest equidistributed meshes). skip_threshold : float, default 0.9 Misalignment threshold for the adapt-on-demand skip. If the @@ -764,7 +658,7 @@ def follow_metric( ``≈ 2 * h_0`` (background cell size). polish_max_iters : int, default 5 Maximum Jacobi (graph-Laplacian) polish iterations - applied AFTER the anisotropic mover. The polish runs + applied AFTER the mover. The polish runs adaptively: each iteration averages every interior vertex toward the mean of its edge neighbours (cell-quality cleanup), and the loop stops as soon as @@ -834,63 +728,27 @@ def follow_metric( gradient_smoothing_length=gradient_smoothing_length, name=name, ) - # Resolve auto coarsening + # Resolve auto coarsening (the R clamp fed to the skip check's + # achievable-target normalisation). if coarsening is None or coarsening == "auto": coar_val = float(refinement) ** (1.0 / mesh.cdim) else: coar_val = float(coarsening) - # Mover's `resolution_ratio` is a SYMMETRIC eigenvalue clamp - # (h ∈ [h0/R, h0·R]) — too loose for either side on its own. - # We pass R = max(refinement, coarsening) so the clamp doesn't - # bind tightly, then rely on the per-cell *rest-size spring* - # (below) to enforce the literal cell-size envelope. R = max(float(refinement), coar_val) - # The spring caps refer to h0 — the **undeformed** mean edge - # length, captured ONCE per mesh and reused (the dt-crash / - # compounding-refinement bug, 2026-05-22 — full story on the - # _FOLLOW_METRIC_H0_CACHE declaration and _mean_edge_length). - _key = id(mesh) - h0 = _FOLLOW_METRIC_H0_CACHE.get(_key) - rest_coords = _FOLLOW_METRIC_REST_CACHE.get(_key) - if h0 is None: - coords = np.asarray(mesh.X.coords) - h0 = _mean_edge_length(mesh.dm, coords) - _FOLLOW_METRIC_H0_CACHE[_key] = h0 - rest_coords = coords.copy() - _FOLLOW_METRIC_REST_CACHE[_key] = rest_coords - if verbose: - uw.pprint(f" follow_metric: captured h0={h0:.4e}, " - f"rest_coords (first call on this mesh)") - - mover_kwargs = dict( - relax=0.2, - n_outer=12, - # Per-cell Lagrangian rest-size spring: literal cell-size - # cap enforced by pulling vertices back toward their - # rest positions when an incident cell exceeds the cap. - # h0 is the undeformed mean edge length. - rest_size_cap_max=h0 * coar_val, - rest_size_cap_min=h0 / float(refinement), - rest_spring_K=1.0, - # Override the mover's internal h0 measurement (which - # would otherwise re-measure on the already-deformed - # mesh and shrink each adapt — the second leg of the - # dt-crash bug surfaced 2026-05-22). - h0_override=h0, - # Override the spring's rest-coords (and the area-floor - # baseline) so they refer to the **truly-undeformed** - # mesh. Otherwise each adapt's "rest" is the previous - # adapt's output, the spring "preserves" each successive - # refinement, and refinement compounds — third leg of - # the dt-crash bug. - rest_coords_override=rest_coords, - ) + # The metric density carries the whole cell-size envelope + # (geomean(rho)=1 by construction), so the MMPDE mover needs no + # per-call clamp translation. theta=0.5 weights Huang's functional + # fully toward EQUIDISTRIBUTION — the right balance for this + # scalar-density contract (the alignment term only matters for + # tensor metrics) and measurably tighter on the requested h_min + # envelope than the mover's general-purpose default (1/3). + mover_kwargs = dict(theta=0.5) if method_kwargs: mover_kwargs.update(method_kwargs) - # Phase-1 remesh redesign: wrap the whole anisotropic-move + polish - # pipeline in a single field-transfer pass at this composite level. + # Phase-1 remesh redesign: wrap the whole move + polish pipeline + # in a single field-transfer pass at this composite level. # The inner smooth_mesh_interior calls see ``mesh._in_remesh_transfer`` # set by the helper and skip their own wrap, so REMAP variables # (including hidden solver history) are transferred exactly once, @@ -904,7 +762,7 @@ def _do_move(): smooth_mesh_interior( mesh, metric=rho, - method="anisotropic", + method="mmpde", method_kwargs={**mover_kwargs, "resolution_ratio": R}, skip_threshold=skip_threshold, verbose=verbose, diff --git a/src/underworld3/meshing/smoothing/graph.py b/src/underworld3/meshing/smoothing/graph.py index 0cc8d83b..83f0f365 100644 --- a/src/underworld3/meshing/smoothing/graph.py +++ b/src/underworld3/meshing/smoothing/graph.py @@ -206,37 +206,13 @@ def gidx(p): return A, dm_scalar, gsection -def _min_incident_edge(dm, coords): - """Per-vertex minimum incident edge length (local-chart - v-pStart order). Used as an optional secondary per-node cap on - the spring step (the primary tangle guard is the coherent global - signed-area backtrack in ``_spring_equilibrium_mover``).""" - pStart, pEnd = dm.getDepthStratum(0) - eStart, eEnd = dm.getDepthStratum(1) - h = np.full(pEnd - pStart, np.inf) - for e in range(eStart, eEnd): - cone = dm.getCone(e) - if len(cone) != 2: - continue - v0, v1 = cone[0], cone[1] - if not (pStart <= v0 < pEnd and pStart <= v1 < pEnd): - continue - i0, i1 = v0 - pStart, v1 - pStart - L = float(np.linalg.norm(coords[i0] - coords[i1])) - if L < h[i0]: - h[i0] = L - if L < h[i1]: - h[i1] = L - return h - def _tri_cells(dm): """Triangle vertex-index triples (local-chart, v-pStart order). Returns an ``(n_tri, 3)`` int array, or ``None`` if the mesh is - not all-triangle (then the global signed-area backtrack is - skipped and only the optional per-node edge cap guards against - tangling). + not all-triangle (callers that need the cell array — fold checks, + boundary-facet extraction — then skip those steps). """ cStart, cEnd = dm.getHeightStratum(0) pStart, pEnd = dm.getDepthStratum(0) @@ -261,113 +237,16 @@ def _signed_areas(coords, tris): - (b[:, 1] - a[:, 1]) * (c[:, 0] - a[:, 0])) -def _cap_step_to_edge_fraction(step, dm, coords, step_frac): - """Per-vertex displacement cap: ``|step_i| <= step_frac * h_i`` with - ``h_i`` the shortest edge incident on vertex ``i``. - - Prevents a mover step from creating LOCAL cell folds near sharp - features (where the source is strongest) without killing the global - motion the way the coherent global signed-area backtrack does. - No-op when ``step_frac`` is ``None`` or non-finite. (The MMPDE mover - has its own, structurally different per-node cap — do not fold it in - here.)""" - if step_frac is None or not np.isfinite(step_frac): - return step - h = _min_incident_edge(dm, coords) - mag = np.linalg.norm(step, axis=1) - cap = float(step_frac) * h - clip = np.isfinite(cap) & (mag > cap) & (mag > 0.0) - sc = np.ones_like(mag) - sc[clip] = cap[clip] / mag[clip] - return step * sc[:, None] - - -def _backtracked_move(old_coords, step, free, tris, project, - area_floor=0.0): - """Coherent global signed-area backtrack shared by the MA / OT / - anisotropic movers. - - Apply ``step`` to the ``free`` vertices at a global scale, halving - the scale (up to 10 times) until no triangle inverts and none drops - to (or below) ``area_floor``. The min-area test is reduced globally - (MPI MIN) so every rank takes the same accept/backtrack branch — the - loop is collective. ``project`` re-imposes the boundary-slip - constraint on each trial (slip vertices snap back to their bounding - surface). - - ``area_floor=0.0`` reproduces the historical flip-only acceptance - (``a1min > 0``) bit-for-bit; the anisotropic mover passes a positive - floor (a fraction of the undeformed median cell area) so - near-degenerate slivers are rejected as well as inverted cells. - ``tris is None`` (non-triangle mesh) applies the step unguarded. - Returns - ------- - (new_coords, scale) : (ndarray, float) - ``scale == 0.0`` means no acceptable move was found and - ``new_coords`` equals ``old_coords``. - """ - scale = 1.0 - new_coords = old_coords.copy() - if tris is not None: - a0 = _signed_areas(old_coords, tris) - orient = np.sign(np.median(a0)) or 1.0 - for _bt in range(10): - trial = old_coords.copy() - trial[free] += scale * step[free] - trial = project(trial) - a1min = _global_min( - (_signed_areas(trial, tris) * orient).min()) - if a1min > area_floor: - new_coords = trial - break - scale *= 0.5 - else: - scale = 0.0 - new_coords = old_coords.copy() - else: - new_coords[free] += step[free] - new_coords = project(new_coords) - return new_coords, scale - - -def _reweight_displacement_radial_tangential(disp, coords, - move_anisotropy): - """Directional move-weighting (approach (2), opt-in; 2D only). - - The annulus node budget is anisotropic — radial is scarce and - pinned, tangential is abundant and free ("spare" angular nodes). A - scalar equidistribution is isotropic and cannot express "prefer - tangential"; rescale the realised displacement in the local - radial / tangential frame (``move_anisotropy = (w_r, w_θ)``) so - the same metric is met mostly by sliding nodes around rather than - crushing radially. Lightweight and solver-consistent — the mover's - operator algebra is untouched, only the realised move is - reweighted. Centre = the coordinate centroid (the origin for a - centred annulus). Degenerate radii (< 1e-30) keep a zero frame and - therefore a zero reweighted move.""" - w_r, w_t = (float(move_anisotropy[0]), - float(move_anisotropy[1])) - ctr = coords.mean(axis=0) - rv = coords - ctr - rn = np.linalg.norm(rv, axis=1) - ok = rn > 1.0e-30 - rhat = np.zeros_like(rv) - rhat[ok] = rv[ok] / rn[ok, None] - that = np.stack([-rhat[:, 1], rhat[:, 0]], axis=1) - d_r = (disp * rhat).sum(axis=1) - d_t = (disp * that).sum(axis=1) - return (w_r * d_r[:, None] * rhat - + w_t * d_t[:, None] * that) def _edge_pairs(dm): """``(n_edge, 2)`` int array of edge endpoint vertex indices in - local-chart (v - pStart) order — the spring network's bars. + local-chart (v - pStart) order. Skips edges whose endpoints are not both in the local vertex - stratum (rank-ghost incomplete edges); the spring path is - serial-exact (see module docstring).""" + stratum (rank-ghost incomplete edges), so the result is the + rank-local edge graph.""" pStart, pEnd = dm.getDepthStratum(0) eStart, eEnd = dm.getDepthStratum(1) pairs = [] @@ -562,8 +441,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 ``_ot_adapt`` (the MMPDE mover itself is - currently 2D-only).""" + 3D boundary-face extraction in :func:`_boundary_facets` (the MMPDE + mover itself is currently 2D-only).""" cStart, cEnd = dm.getHeightStratum(0) pStart, pEnd = dm.getDepthStratum(0) tets = [] @@ -606,9 +485,8 @@ def _owned_cell_mask(dm): def _min_incident_edge_nd(cells, coords): """Dimension-general shortest-incident-edge per vertex. ``cells`` is (n_cells, d+1); returns (n_verts,). Used by the MMPDE per-node step - cap. (The 2D-only ``_min_incident_edge`` reads the DM directly; this - works for tets too and takes an explicit cell array so the caller can - restrict the stencil.)""" + cap. Works for tets too and takes an explicit cell array so the + caller can restrict the stencil.""" n_verts = coords.shape[0] ncorner = cells.shape[1] v = np.full(n_verts, np.inf) @@ -619,3 +497,219 @@ def _min_incident_edge_nd(cells, coords): np.minimum.at(v, cells[:, a], e) np.minimum.at(v, cells[:, b], e) return v + + +# --------------------------------------------------------------------------- +# Boundary-facet and boundary-slip primitives. +# +# Relocated from ``meshing/_ot_adapt.py`` when the OT-reset adapt was +# retired (2026-07; the spring / Monge-Ampere / OT / anisotropic-Winslow +# movers are superseded by the MMPDE mover). These helpers are mover- +# agnostic and remain in use by the MMPDE mover, the mesh's +# ``boundary_slip`` orchestration, and ``BoundingSurface`` facet restore. +# --------------------------------------------------------------------------- + +def _slip_normals(mesh, boundary_coords): + """Unit outward normals at ``boundary_coords`` from the projected + boundary-normal field. + + Re-projects ``mesh._projected_normals`` (``mesh.Gamma_P1``) first so the + normals reflect the mesh's *current* coordinates — the projected field is + stale after any deform. Returns ``(normals, valid)`` where ``normals`` is + ``(k, cdim)`` and ``valid`` is a boolean mask; ``valid`` is ``False`` for + nodes with a degenerate (zero / non-finite) normal (e.g. box corners + where opposing face normals cancel, or an occasional unlocatable vertex). + Such nodes should be pinned, not slipped. + """ + cdim = mesh.cdim + n = np.zeros((boundary_coords.shape[0], cdim)) + try: + mesh._update_projected_normals() + n = np.asarray( + uw.function.evaluate(mesh.Gamma_P1, boundary_coords) + ).reshape(-1, cdim) + except Exception: + # Projection unavailable / degenerate on this mesh — fall back to + # all-pinned boundaries (valid stays all-False below). + n = np.zeros((boundary_coords.shape[0], cdim)) + mag = np.linalg.norm(n, axis=1) + valid = np.isfinite(mag) & (mag > 0.5) + out = np.zeros_like(n) + out[valid] = n[valid] / mag[valid, None] + return out, valid + + +def _boundary_facets(mesh, cdim): + """Boundary facets + opposite cell-vertex, found from the cell topology. + + For each cell, every facet (edge in 2D, triangle in 3D) is a candidate + boundary facet; one that occurs in **exactly one** cell is on the + boundary. Returns ``(facets, opp)`` where ``facets`` is ``(n_bnd, k)`` + (``k=2`` for 2D edges, ``k=3`` for 3D triangles) and ``opp`` is the + cell vertex opposite each facet — used to orient the facet normal + outward. Returns ``(None, None)`` for non-simplicial meshes. + """ + if cdim == 2: + cells = _tri_cells(mesh.dm) + if cells is None: + return None, None + rows = [] + for k in range(3): + v0 = cells[:, k]; v1 = cells[:, (k + 1) % 3] + vopp = cells[:, (k + 2) % 3] + vmin = np.minimum(v0, v1); vmax = np.maximum(v0, v1) + rows.append(np.column_stack([vmin, vmax, vopp])) + e = np.vstack(rows) + idx = np.lexsort((e[:, 1], e[:, 0])) + e = e[idx] + same_prev = np.zeros(len(e), dtype=bool) + same_prev[1:] = ((e[1:, 0] == e[:-1, 0]) + & (e[1:, 1] == e[:-1, 1])) + same_next = np.zeros(len(e), dtype=bool) + same_next[:-1] = same_prev[1:] + bnd_mask = (~same_prev) & (~same_next) + bnd = e[bnd_mask] + return bnd[:, :2], bnd[:, 2] + if cdim == 3: + cells = _tet_cells(mesh.dm) + if cells is None: + return None, None + rows = [] + for k in range(4): + others = [(k + 1) % 4, (k + 2) % 4, (k + 3) % 4] + tri = np.sort(np.column_stack( + [cells[:, others[0]], cells[:, others[1]], + cells[:, others[2]]]), axis=1) + rows.append(np.column_stack([tri, cells[:, k]])) + f = np.vstack(rows) + idx = np.lexsort((f[:, 2], f[:, 1], f[:, 0])) + f = f[idx] + same_prev = np.zeros(len(f), dtype=bool) + same_prev[1:] = ((f[1:, 0] == f[:-1, 0]) + & (f[1:, 1] == f[:-1, 1]) + & (f[1:, 2] == f[:-1, 2])) + same_next = np.zeros(len(f), dtype=bool) + same_next[:-1] = same_prev[1:] + bnd_mask = (~same_prev) & (~same_next) + bnd = f[bnd_mask] + return bnd[:, :3], bnd[:, 3] + return None, None + + +def _all_boundary_labels(mesh): + """Named codim-1 boundary labels of the mesh, skipping the synthetic / + non-geometric ones (``All_Boundaries``, ``Null_Boundary``, and the + Annulus single-point ``Centre`` pseudo-label that hard-aborts PETSc).""" + skip = {"All_Boundaries", "Null_Boundary", "Centre"} + out = [] + try: + names = [b.name for b in mesh.boundaries] + except Exception: + # Meshes without a boundaries enum (e.g. hand-built DMs) simply + # expose no named labels — treat as label-free. + names = [] + for nm in names: + if nm in skip: + continue + out.append(nm) + return tuple(out) + + +def _resolve_slip(mesh, slip_spec): + """Resolve the ``slip_spec`` (the value passed as ``boundary_slip`` / + ``slip_surfaces``) into a tuple of named slip-surface labels, and + pre-touch ``mesh.Gamma_P1`` so the projected-normal field ``_n_proj`` + exists BEFORE any mover builds its solver DM (creating that MeshVariable + mid-mover would stale the DM handle — see project_uw3_smoother_footguns). + + Accepted forms (back-compatible): + * ``True`` / truthy / legacy ``'ring'``,``'box'`` strings → ALL named + codim-1 boundary surfaces slip. + * ``False`` / ``None`` / ``[]`` → no slip (pin all boundaries). + * a label name, or a list of label names → only those surfaces slip. + * a ``dict`` ``{label: snap_bool}`` → those labels slip; ``snap_bool`` + is the per-surface return-to-bounds flag (``False`` = FREE surface, + slip but do not snap back). The dict keys are the slip labels. + + Returns the tuple of slip-surface label names (possibly empty). + """ + if slip_spec is None or slip_spec is False: + return () + if slip_spec is True: + labels = _all_boundary_labels(mesh) + elif isinstance(slip_spec, dict): + labels = tuple(slip_spec.keys()) + elif isinstance(slip_spec, str): + s = slip_spec.strip().lower() + if s in ("ring", "box", "axes", "axis", "true", "on", "1", "all"): + labels = _all_boundary_labels(mesh) + elif s in ("false", "off", "0", "none", ""): + return () + else: + labels = (slip_spec,) # a single explicit label name + else: + # an iterable of label names + labels = tuple(slip_spec) + if labels: + # Pre-create the projected-normal field (footgun-safe; see docstring). + try: + _ = mesh.Gamma_P1 + except Exception: + # No projected-normal field on this mesh (e.g. a manifold or + # label-free DM) — slip resolution still returns the labels; + # the caller's per-node normal lookup pins degenerate nodes. + pass + return labels + + +def _nearest_on_facets_2d(pts, seg): + """Closest point on a set of 2D line segments. ``pts`` (m,2), + ``seg`` (nf,2,2). Returns (m,2) closest points (over all segments).""" + a = seg[:, 0]; b = seg[:, 1] # (nf,2) + ab = b - a + ab2 = np.einsum('fi,fi->f', ab, ab) + ab2 = np.where(ab2 > 1.0e-30, ab2, 1.0) + out = np.empty_like(pts) + for i, p in enumerate(pts): + t = np.clip(((p - a) * ab).sum(axis=1) / ab2, 0.0, 1.0) + proj = a + t[:, None] * ab # (nf,2) + d2 = ((proj - p) ** 2).sum(axis=1) + out[i] = proj[d2.argmin()] + return out + + +def _nearest_on_facets_3d(pts, tri): + """Closest point on a set of 3D triangles. ``pts`` (m,3), + ``tri`` (nf,3,3). Returns (m,3). Per-point loop, vectorised over + triangles via the standard region-based closest-point algorithm.""" + A = tri[:, 0]; B = tri[:, 1]; C = tri[:, 2] + AB = B - A; AC = C - A + out = np.empty_like(pts) + for i, p in enumerate(pts): + AP = p - A + d1 = np.einsum('fi,fi->f', AB, AP) + d2 = np.einsum('fi,fi->f', AC, AP) + BP = p - B + d3 = np.einsum('fi,fi->f', AB, BP) + d4 = np.einsum('fi,fi->f', AC, BP) + CP = p - C + d5 = np.einsum('fi,fi->f', AB, CP) + d6 = np.einsum('fi,fi->f', AC, CP) + va = d3 * d6 - d5 * d4 + vb = d5 * d2 - d1 * d6 + vc = d1 * d4 - d3 * d2 + denom = va + vb + vc + denom = np.where(np.abs(denom) > 1.0e-30, denom, 1.0) + v = vb / denom + w = vc / denom + # interior barycentric point; clamp handles edge/vertex regions well + # enough for a small return-to-bounds correction on convex surfaces. + v = np.clip(v, 0.0, 1.0); w = np.clip(w, 0.0, 1.0) + s = v + w + over = s > 1.0 + v = np.where(over, v / np.where(s > 0, s, 1.0), v) + w = np.where(over, w / np.where(s > 0, s, 1.0), w) + proj = A + v[:, None] * AB + w[:, None] * AC + dd = ((proj - p) ** 2).sum(axis=1) + out[i] = proj[dd.argmin()] + return out diff --git a/src/underworld3/meshing/smoothing/metrics.py b/src/underworld3/meshing/smoothing/metrics.py index 58999a68..2bd0df99 100644 --- a/src/underworld3/meshing/smoothing/metrics.py +++ b/src/underworld3/meshing/smoothing/metrics.py @@ -274,8 +274,7 @@ def metric_density_from_gradient( :func:`smooth_mesh_interior`:: rho = metric_density_from_gradient(mesh, T, amp=8.0) - smooth_mesh_interior(mesh, metric=rho, - method="anisotropic") + smooth_mesh_interior(mesh, metric=rho, method="mmpde") The projector/fields are cached per ``(mesh, degree, name, topology)``, so calling this **every step** in an adaptive loop @@ -340,8 +339,8 @@ def metric_density_from_gradient( refinement only into the steepest fronts. degree : int, default 1 Polynomial degree of the projected-gradient / density - fields (1 matches the anisotropic mover's default - ``aux_degree``). + fields (P1 is what the mover's per-vertex metric + evaluation samples). name : str, optional Cache disambiguator. Pass distinct names if you build several independent gradient metrics on the *same* mesh @@ -518,7 +517,7 @@ def metric_density_from_gradient( # * "arc-length" — smooth arc-length monitor # ρ = √(1 + (A·|∇field|/g_hi)²), clipped to the envelope. # Grades continuously from ρ=1 in flat regions (no clip - # kink) → cleaner OT / Monge–Ampère meshes. + # kink) → the smoothest equidistributed meshes. # # ``coarsening="auto"`` uses the budget-conserving minimum # ``refinement^(1/d)`` — the smallest coarsening that @@ -593,7 +592,7 @@ def metric_density_from_gradient( # Smooth arc-length monitor rho = sqrt(1 + (A*ghat)^2), # ghat = |grad field|/g_hi, A = sqrt(ref^4 - 1) so rho = ref^2 at # the hi-percentile gradient. Grades continuously from rho=1 in - # flat regions (no clip kink) -> cleaner OT / Monge-Ampere meshes. + # flat regions (no clip kink) -> the smoothest equidistributed meshes. A = np.sqrt(max(ref_val ** 4 - 1.0, 0.0)) ghat = gmag / max(g_hi, 1.0e-30) rho_al = np.sqrt(1.0 + (A * ghat) ** 2) diff --git a/src/underworld3/meshing/smoothing/mmpde.py b/src/underworld3/meshing/smoothing/mmpde.py index bb5ace24..b20ab687 100644 --- a/src/underworld3/meshing/smoothing/mmpde.py +++ b/src/underworld3/meshing/smoothing/mmpde.py @@ -107,9 +107,9 @@ def _mmpde_mover(mesh, metric, pinned_labels, verbose, Because `G → ∞` as `det𝕁 → 0` the map is non-folding (Math. Comp. 87 (2018) 1887); because it is the inverse map of a convex computational domain it genuinely *clusters and aligns* to `M` — a thin strip on a - fault, not the isotropic centre-of-gravity blob the scalar MA mover - produces, and not the non-clustering smooth of the decoupled - `_winslow_anisotropic`. See + fault, not an isotropic centre-of-gravity blob, and not a + non-clustering decoupled Winslow smooth (the scalar Monge-Ampère and + anisotropic-Winslow movers it superseded were retired 2026-07). See ``docs/developer/design/anisotropic-mmpde-mover.md``. ``metric`` is the SPD `d×d` metric tensor: a sympy `Matrix` (function @@ -305,7 +305,7 @@ def _eval_M(pts): # (below). Pre-touch Gamma_P1 here so the projected-normal MeshVariable # exists before any DM snapshot (footgun-safe; redundant with the central # pre-touch in smooth_mesh_interior, kept as defence-in-depth). - from underworld3.meshing._ot_adapt import _resolve_slip + from .graph import _resolve_slip _slip_pretouch = _resolve_slip(mesh, boundary_slip) # pre-touch Gamma_P1 before DM build # Reference edge matrices (fixed) for the owned cells. diff --git a/src/underworld3/meshing/smoothing/monge_ampere.py b/src/underworld3/meshing/smoothing/monge_ampere.py deleted file mode 100644 index 6f6f4f9c..00000000 --- a/src/underworld3/meshing/smoothing/monge_ampere.py +++ /dev/null @@ -1,836 +0,0 @@ -"""Monge–Ampère (BFO convex-branch) equidistribution machinery, the -linear OT improvement step, and the mover sub-solver wiring that the -anisotropic mover also reuses. See the package docstring for the -module map. -""" - -import numpy as np - -import underworld3 as uw - -from .graph import (_tri_cells, _signed_areas, - _cap_step_to_edge_fraction, _backtracked_move, - _reweight_displacement_radial_tangential, - _global_sum, _global_max, _global_mean) - - -# ====================================================================== -# Monge–Ampère mesh-equidistribution machinery (PRESERVED, not the -# default metric path). Exhaustively investigated 2026-05-16: every -# FE-MA-potential variant (linear / recovered-Hessian smoothed & -# variational / BFO convex-branch + damping / outer composition) -# caps at deep/near ≈ 1.07 for an 8× target vs an exact ~10× — see -# the project memory and scripts/ma_*.py. Kept because (a) the -# "bit-identical across variants" result suggests a common missing -# ingredient worth understanding, and (b) the elastic-spring -# redistribution may work as a *preconditioner* for the MA solve -# (a graded starting mesh might let MA escape the weak branch) — -# an open investigation. Call _monge_ampere_mover() directly to use. -# ====================================================================== - -# Cached MA solver state keyed by (mesh-id, pinned-labels, topology): -# the φ Poisson, the variational Hessian-recovery solver, ∇φ -# projector, the ρ_cur proxy field. Rebuilt on a topology change. -_WINSLOW_CACHE: dict = {} - -# Sign of the BFO source vs UW3's SNES_Poisson convention -# (SNES_Poisson F0 = -f, strong form Δφ = -ps.f). With this sign the -# validated linear first iterate Δφ = (c/ρ_tgt - 1) grades the right -# way (nodes toward high target density). -_EQUIDIST_SIGN = -1.0 - -_HESSIAN_CLASS = None - - -# Cached state for the OT-improvement-step path (one weighted- -# Poisson per call). Keyed like the other movers; same lifetime. -_OT_CACHE: dict = {} - - -def _use_direct_solver(solver, singular=False, elliptic=True): - r"""Force a cached MA sub-solver onto a sparse **direct** factorisation - (MUMPS LU) instead of the UW3 default GMRES + GAMG. - - **Parallel safety (parallel-singular-corruption, 2026-05-27):** MUMPS - (the only parallel LU in this build) corrupts the heap when the same - factorisation path is exercised over *repeated* solves at np >= 3 — a - probabilistic SEGV/SIGBUS or MPI deadlock (the UW3 default GMRES+GAMG, which - never calls MUMPS, is clean). Since the movers re-solve in a Picard/outer - loop, MUMPS is unusable in parallel here. Under MPI this function therefore - falls back to the MUMPS-free iterative path (:func:`_use_iterative_solver`); - the direct MUMPS path below is kept for the validated **serial** efficiency - lever. ``elliptic`` is forwarded to the iterative fallback (φ-Poisson → - GAMG; mass systems → CG+Jacobi). - - Why this is the dominant MA-efficiency lever (profiled 2026-05-17, - res-16 Annulus, AMP=8, warm re-call): the Picard loop fixes the - mesh, so the φ-Poisson Laplacian and the Hessian-recovery SPD mass - matrix are *constant operators* re-solved ~40× with only the RHS - changing. With GAMG, every ``solve()`` pays a full multigrid - **setup** (the constant near-nullspace re-attach forces it) — the - Hessian solve alone was ~0.93 s/iter ≈ 37 s. These problems are - tiny (≲10⁴ DOF); MUMPS factorises in milliseconds and the per-iter - cost collapses to a back-substitution. A direct solve is also - *exact* (machine precision, tighter than the GMRES rtol), so the - Picard fixed point — hence the grading/quality — is unchanged. - - ``singular=True`` (the pure-Neumann φ Poisson): MUMPS null-pivot - detection (ICNTL(24)=1) handles the rank-1-deficient operator; the - ``constant_nullspace`` hook still removes the constant mode from - the RHS/solution, so the result is the same consistent solution - the iterative path produced — but it also eliminates the - GAMG-on-pure-Neumann ``DIVERGED_LINEAR_SOLVE`` re-solve pathology. - """ - # Parallel: MUMPS-repeated corrupts the heap (see docstring) — use the - # MUMPS-free iterative path instead. Serial keeps the fast direct solve. - if uw.mpi.size > 1: - _use_iterative_solver(solver, singular=singular, elliptic=elliptic) - return - - o = solver.petsc_options - # These three sub-problems are *linear* (φ Poisson with the Hessian - # source frozen; the SPD Hessian-recovery mass system; the ∇φ - # projection) → one KSP solve, no Newton line-search / 2nd iterate - # (which was doubling work and emitting spurious - # ``DIVERGED_LINEAR_SOLVE`` after 2 iters). - o["snes_type"] = "ksponly" - # ksponly does exactly ONE linear KSP solve (no Newton). Default - # snes_max_it leaves snes->iter=0, so if a converged-reason - # viewer is on (a user's global -snes_converged_reason, an outer - # debug flag, …) PETSc mislabels the *successful* linear solve - # as "DIVERGED_MAX_IT iterations 0" and floods the log with - # phantom failures. snes_max_it=1 ⇒ the single solve counts as - # one converged iteration ⇒ reason = CONVERGED, not a fake - # DIVERGED. Numerically inert (the KSP solve is identical) — - # purely stops these linear sub-solves masquerading as failures. - o["snes_max_it"] = 1 - # The Picard loop fixes the mesh, so the operator is **constant** - # across the ~40 inner solves — only the RHS changes. Lag the - # Jacobian (compute once, reuse) and the preconditioner (factorise - # once, reuse): every subsequent inner solve collapses to a MUMPS - # back-substitution. A fresh ``solver.solve()`` after - # ``_deform_mesh`` rebuilds the SNES (is_setup=False) so the lag - # counter resets and the operator is correctly re-factorised on the - # first solve of the next call — the reuse is confined to the loop - # where the mesh genuinely does not move. - o["snes_lag_jacobian"] = -2 - o["snes_lag_preconditioner"] = -2 - o["ksp_type"] = "preonly" - o["pc_type"] = "lu" - o["pc_factor_mat_solver_type"] = "mumps" - if singular: - o["mat_mumps_icntl_24"] = 1 # null-pivot detection - o["mat_mumps_icntl_25"] = 0 # one solution of the singular sys - # GAMG-only keys are inert once pc_type≠gamg; drop them so the - # effective option set is exactly what is documented. - for k in ("pc_gamg_type", "pc_gamg_repartition", "pc_mg_type", - "pc_gamg_agg_nsmooths", "mg_levels_ksp_max_it", - "mg_levels_ksp_converged_maxits"): - try: - o.delValue(k) - except Exception: - pass - - -def _use_iterative_solver(solver, singular=False, elliptic=True): - r"""Parallel-scalable alternative to ``_use_direct_solver``: keep - the *same factor/setup-once-reuse pattern* (the real efficiency - lever) but with an **iterative** PC so it scales beyond the - serial / modest-size regime where sparse direct factorisation is - viable (this PETSc build has only MUMPS + serial builtin LU — no - hypre / SuperLU_DIST). - - The Picard loop fixes the mesh ⇒ the operator is constant across - the ~25 inner solves; ``snes_lag_jacobian=-2`` / - ``snes_lag_preconditioner=-2`` build the PC **once per - ``_monge_ampere_mover`` call** and reuse it for every inner solve - (the GAMG hierarchy / Jacobi diagonal is *not* rebuilt per - iteration — that per-iter GAMG re-setup was the original ~0.9 s - Hessian cost). ``_deform_mesh`` resets ``is_setup`` so the lag - counter resets and the PC is correctly rebuilt on the next call's - first solve. Combined with a Krylov **warm start** from the - previous Picard φ (caller passes ``zero_init_guess=False``), the - inner solves are a handful of CG iterations on an already-built - hierarchy. - - ``elliptic=True`` (the φ-Poisson Laplacian): CG + GAMG with the - constant near-nullspace (already attached via - ``constant_nullspace`` — GAMG needs it for the pure-Neumann - operator). ``elliptic=False`` (the SPD Hessian-recovery / ∇φ mass - systems): a mass matrix is spectrally trivial — CG + Jacobi - converges in a few iterations with **no** hierarchy setup, fully - parallel; GAMG there would be wasted setup. - - Numerics: an iterative solve to a tight ``ksp_rtol`` reproduces - the BFO Picard fixed point — hence the grading — to well within - its 4-dp precision (validated against the direct path); it is a - *cost/parallelism* change, not a formulation change. - """ - o = solver.petsc_options - o["snes_type"] = "ksponly" - # See _use_direct_solver: snes_max_it=1 stops a converged-reason - # viewer mislabelling these linear ksponly sub-solves as - # "DIVERGED_MAX_IT iterations 0". Numerically inert. - o["snes_max_it"] = 1 - o["snes_lag_jacobian"] = -2 - o["snes_lag_preconditioner"] = -2 - # Krylov choice is per-operator (set in the branches below): - # * elliptic φ-Poisson → FGMRES. The UW3 DMPlex-FEM assembly + - # Neumann/nullspace handling does not guarantee an *exactly* - # symmetric operator, and the GAMG **SOR smoother is - # non-symmetric**, so the preconditioner is non-SPD — CG's - # assumptions are violated (it only "worked" here by - # robustness margin). FGMRES tolerates a non-symmetric - # operator *and* a varying/non-symmetric preconditioner. - # * mass systems (Hessian recovery, ∇φ projection) → CG: a - # consistent mass matrix with a Jacobi PC is provably SPD and - # symmetric, so CG is correct and the cheapest option. - # Inner solve inside an outer BFO Picard — it tolerates inexact - # inner solves (inexact-Picard); 1e-7 is far tighter than the - # Picard increment near convergence (~1e-4) so the fixed point — - # hence the grading — is unchanged, at a fraction of the iters a - # direct-path-matching 1e-10 would need. - o["ksp_rtol"] = 1.0e-7 - o["ksp_atol"] = 1.0e-12 - o["pc_factor_mat_solver_type"] = "" # not a direct solve - try: - o.delValue("pc_factor_mat_solver_type") - o.delValue("mat_mumps_icntl_24") - o.delValue("mat_mumps_icntl_25") - except Exception: - pass - if elliptic: - # P3 pure-Neumann Laplacian: plain agg-GAMG with a weak - # Jacobi/Chebyshev smoother needs ~280 iters here. A stronger - # SOR smoother with more sweeps + smoothed aggregation cuts - # that ~4×; the hierarchy is still built only once per call - # (lagged), so the extra setup is amortised over the ~25 - # reused inner solves. SOR ⇒ non-symmetric PC ⇒ FGMRES. - o["ksp_type"] = "fgmres" - o["ksp_gmres_restart"] = 100 # > the ~75-iter solve - o["pc_type"] = "gamg" - o["pc_gamg_type"] = "agg" - o["pc_gamg_agg_nsmooths"] = 1 - o["pc_gamg_threshold"] = 0.02 - o["mg_levels_ksp_type"] = "richardson" - o["mg_levels_pc_type"] = "sor" - o["mg_levels_ksp_max_it"] = 4 - # GAMG coarse solve. MUMPS (parallel LU) corrupts the heap over repeated - # parallel solves (parallel-singular-corruption) — so in parallel use a - # MUMPS-free coarse: `redundant` replicates the (tiny) coarse grid to - # every rank and solves it with a dense SVD, which is robust on the - # singular pure-Neumann coarse and never calls MUMPS (verified clean + - # convergent at np=5). Serial keeps the fast MUMPS coarse. - for k in ("mg_coarse_pc_factor_mat_solver_type", - "mg_coarse_redundant_pc_type"): - try: o.delValue(k) - except Exception: pass - if uw.mpi.size > 1: - o["mg_coarse_pc_type"] = "redundant" - o["mg_coarse_redundant_pc_type"] = "svd" - else: - o["mg_coarse_pc_type"] = "lu" - o["mg_coarse_pc_factor_mat_solver_type"] = "mumps" - else: - o["ksp_type"] = "cg" # consistent mass = SPD - o["pc_type"] = "jacobi" # mass matrix → trivial - for k in ("ksp_gmres_restart", "pc_gamg_type", - "pc_gamg_agg_nsmooths", "pc_gamg_threshold", - "mg_levels_ksp_type", "mg_levels_pc_type", - "mg_levels_ksp_max_it", "mg_coarse_pc_type", - "mg_coarse_pc_factor_mat_solver_type", - "mg_coarse_redundant_pc_type"): - try: - o.delValue(k) - except Exception: - pass - - -def _solver_wiring(linear_solver): - """Option-wiring function for a cached mover sub-solver. - - ``linear_solver="gamg"`` wires the iterative, parallel-scalable - option set (:func:`_use_iterative_solver`); anything else wires the - serial MUMPS factor-once set (:func:`_use_direct_solver`, which - itself falls back to the iterative path under MPI — the - MUMPS-heap-corruption guard).""" - if linear_solver == "gamg": - def _wire(s, singular=False, elliptic=True): - _use_iterative_solver(s, singular, elliptic) - else: - def _wire(s, singular=False, elliptic=True): - _use_direct_solver(s, singular, elliptic) - return _wire - - -def _warm_start_krylov(linear_solver): - """True when a mover's inner solves should WARM-START the Krylov - iteration from the previous solution (pass ``zero_init_guess=False``): - - * the GAMG path always — the hierarchy is built once (lagged) and - the solution changes slowly under relaxation, so the warm start - leaves only a handful of Krylov iterations per inner solve; - * the "direct" path under MPI — :func:`_use_direct_solver` silently - routes to the iterative solver there (MUMPS-heap-corruption - guard), so the warm start pays exactly as it does for GAMG. - - The serial direct path is an exact factorisation, indifferent to - the initial guess.""" - return (linear_solver == "gamg" - or (linear_solver == "direct" and uw.mpi.size > 1)) - - -def _patch_volumes(tris, coords, n_verts, vol_field=None): - """Per-vertex dual-patch area: a node's share (1/3) of every - incident triangle's |area|. ρ_cur ∝ 1/patch for the (opt-in, - n_outer>1) outer MA composition; at equidistribution - ``patch · ρ_tgt`` is uniform. - - This quantity is exactly the **lumped P1 mass diagonal** ``M_ii = ∫ N_i dV``. - The hand-rolled local sum below is serial-exact, but **under-counts shared - vertices on rank-partition boundaries in parallel** — each rank only adds its - own incident triangles and never sums the neighbouring rank's. So in parallel - we assemble it through the FE mass matrix instead (``_lumped_vertex_volumes``), - where PETSc does the cross-rank ``localToGlobal(ADD)`` for us. Requires the - P1 ``vol_field``; falls back to the local sum when it is not supplied. - """ - if vol_field is not None and uw.mpi.size > 1: - return _lumped_vertex_volumes(vol_field) - area = np.abs(_signed_areas(coords, tris)) / 3.0 - patch = np.zeros(n_verts, dtype=np.float64) - for k in range(3): - np.add.at(patch, tris[:, k], area) - patch[patch <= 0.0] = patch[patch > 0.0].mean() - return patch - - -def _lumped_vertex_volumes(vol_field): - r"""Parallel-correct per-vertex dual-patch volume = the lumped P1 mass - diagonal ``M_ii = ∫ N_i dV`` of ``vol_field``'s (P1, continuous, scalar) - space, assembled via the FE mass matrix so the cross-rank sum over shared - partition-boundary vertices is done by PETSc — unlike the hand-rolled local - sum in :func:`_patch_volumes`, which under-counts those vertices in parallel. - - Identity: by partition of unity (``Σ_j N_j ≡ 1``) the consistent mass matrix - has row sums ``Σ_j M_ij = ∫ N_i Σ_j N_j = ∫ N_i dV``, i.e. the lumped diagonal - is ``M·1``. - - TODO(petsc4py): PETSc has a purpose-built - ``DMCreateMassMatrixLumped(dm, &llm, &lm)`` that returns this lumped diagonal - directly (with the cross-rank ADD built in), but petsc4py (3.25) does not bind - it yet — only the *consistent* ``DM.createMassMatrix`` is exposed, hence the - ``M·1`` below. Replace this body with ``subdm.createMassMatrixLumped()`` once - petsc4py exposes that DM method. - - Returns a per-vertex numpy array in ``vol_field``'s local DOF ordering (the - same depth-0 vertex ordering the movers use for ``vol_field.array``). - """ - mesh = vol_field.mesh - indexset, subdm = mesh.dm.createSubDM(vol_field.field_id) - M = subdm.createMassMatrix(subdm) # consistent P1 mass (FE-assembled, parallel-correct) - ones = M.createVecRight() - ones.set(1.0) - lumped = M.createVecLeft() - M.mult(ones, lumped) # M·1 = row sums = lumped diagonal - lvec = subdm.getLocalVec() - subdm.globalToLocal(lumped, lvec, addv=False) - out = np.asarray(lvec.array).copy() - subdm.restoreLocalVec(lvec) - for obj in (M, ones, lumped, indexset, subdm): - try: - obj.destroy() - except Exception: - pass - pos = out > 0.0 - if not pos.all(): - out[~pos] = out[pos].mean() - return out - - -def _hessian_recovery_class(): - r"""Lazily build (and memoise) the variationally-consistent - Hessian-recovery solver class. - - Recovers ``H_ij ≈ ∂²φ/∂x_i∂x_j`` from an external scalar field - ``φ`` by the *weak* (integrated-by-parts) form — the plan's - :math:`R_H`: ``∫H_ij τ_ij + ∫(∂φ/∂x_i)(∂τ_ij/∂x_j) = 0`` ⇒ - ``H_ij = ∂²φ/∂x_i∂x_j``. Only **first** derivatives of ``φ`` - appear (UW3 forbids second derivatives of mesh-variable - functions); the operator is the SPD mass matrix (no nullspace). - Defined lazily to avoid an import cycle (meshing→systems/cython). - """ - global _HESSIAN_CLASS - if _HESSIAN_CLASS is not None: - return _HESSIAN_CLASS - - import sympy - from underworld3.cython.generic_solvers import SNES_MultiComponent - from underworld3.utilities._api_tools import Template - - class _HessianRecovery(SNES_MultiComponent): - def __init__(self, mesh, phi_field, degree=2, verbose=False): - self._phi = phi_field - super().__init__( - mesh, n_components=mesh.cdim * mesh.cdim, - degree=degree, verbose=verbose) - self._smoothing = sympy.sympify(0) - self._constitutive_model = ( - uw.constitutive_models.Constitutive_Model( - self.Unknowns)) - - def _hessian_source(self): - cdim = self.mesh.cdim - X = self.mesh.CoordinateSystem.X - phi = self._phi.sym[0] - rows = [] - for i in range(cdim): - for j in range(cdim): - row = [sympy.Integer(0)] * cdim - row[j] = phi.diff(X[i]) - rows.append(row) - return sympy.Matrix(rows) - - F0 = Template( - r"f_0\left(\mathbf{u}\right)", - lambda self: self.u.sym, - "Hessian-recovery mass term: f_0 = H.") - - F1 = Template( - r"\mathbf{F}_1\left(\mathbf{u}\right)", - lambda self: self._hessian_source(), - "Hessian-recovery weak source: F_1 = e_j ∂φ/∂x_i.") - - _HESSIAN_CLASS = _HessianRecovery - return _HESSIAN_CLASS - - -def _monge_ampere_mover(mesh, metric, pinned_labels, verbose, - n_outer=1, n_picard=25, relax=1.0, - step_frac=None, picard_relax=0.4, - outer_tol=1.0e-3, boundary_slip=False, - linear_solver="direct", phi_degree=2, - move_anisotropy=None, - target_side_rho=False): - r"""Metric-driven mesh equidistribution — Benamou–Froese–Oberman - convex-branch Monge–Ampère (PRESERVED; not the default path). - - Solves ``det(I+D²φ)=g``, ``g=c·ρ_cur/ρ_tgt``, by a damped Picard - on the convex-branch source - ``Δφ = √((φxx−φyy)²+4φxy²+4g) − 2`` (the +√ selects the Brenier - branch), with the variationally-consistent recovered Hessian - (``_hessian_recovery_class``) and the pure-Neumann - ``constant_nullspace`` φ Poisson. ``n_outer>1`` composes maps - (recompute ρ_cur from patch volumes each step). Moves nodes by - ∇φ with a coherent global signed-area backtrack. - - Efficiency (2026-05-17): the φ Poisson and the SPD Hessian-recovery - mass system are *constant operators* within the Picard loop (the - mesh is fixed; only the RHS changes). ``_use_direct_solver`` puts - both on MUMPS LU with a lagged (compute-once) factorisation, so the - inner iterations are back-substitutions — see that function's - docstring. ``n_picard`` defaults to 25: the deep/near grading is - flat from iter ≈20 (4-dp identical at AMP 8 & 20), so 40 was pure - overhead. Net: ~10× faster, grading/quality bit-for-bit unchanged. - - ``phi_degree`` defaults to **2** (was 3). The deep/near grading - is set by the φ *order*, not the solver: P2 ≡ P3 to ~3 dp across - AMP 0/2/8/20 (matches the recorded baseline; AMP=0 no-op exact; - no tangle) while P2 halves the cost (smaller matrices — also - helps the direct factorisation scale). P1 is **not** - grading-equivalent (≈1.40 vs 1.71 at AMP=8 — ~18 % weaker); P2 - is the floor. ``linear_solver="gamg"`` is an experimental, - documented-fragile parallel prototype (P3 was a major GAMG - confound; even at P2 GAMG re-solve is erratic — see the design - doc); ``"direct"`` (MUMPS, MPI-parallel) is the validated path. - - Grading: redistribution with a fixed node count reaches deep/near - ≈1.5–1.8× for an 8–20× density target (the exact OT ~10× needs - *more nodes* — a topology change, not this smoother). ``n_outer=1`` - is the safe default (AMP=0 exact no-op, never tangles). See the - project memory + scripts/ma_*.py / ma_cost_grading.py. - """ - import sympy - - pinned_labels = tuple(pinned_labels) - dm = mesh.dm - pStart, pEnd = dm.getDepthStratum(0) - cStart, cEnd = dm.getHeightStratum(0) - cone_size = dm.getConeSize(cStart) if cEnd > cStart else 0 - if linear_solver not in ("direct", "gamg"): - raise ValueError( - f"linear_solver must be 'direct' or 'gamg', " - f"got {linear_solver!r}") - phi_degree = int(phi_degree) - aux_degree = max(1, phi_degree - 1) # ∇φ / recovered-Hessian - key = (id(mesh), pinned_labels, - pEnd - pStart, cEnd - cStart, cone_size, - linear_solver, phi_degree) - - cdim = mesh.cdim - - cache = _WINSLOW_CACHE.get(key) - if cache is None: - _wire = _solver_wiring(linear_solver) - phi = uw.discretisation.MeshVariable( - f"winslow_phi_{id(mesh)}", mesh, - vtype=uw.VarType.SCALAR, degree=phi_degree, - continuous=True) - ps = uw.systems.Poisson(mesh, phi) - ps.constitutive_model = uw.constitutive_models.DiffusionModel - ps.constitutive_model.Parameters.diffusivity = 1.0 - ps.constant_nullspace = True - _wire(ps, singular=True, elliptic=True) - hsolver = _hessian_recovery_class()( - mesh, phi, degree=aux_degree, verbose=False) - hsolver.tolerance = 1.0e-6 - _wire(hsolver, elliptic=False) - vol_field = uw.discretisation.MeshVariable( - f"winslow_vol_{id(mesh)}", mesh, - vtype=uw.VarType.SCALAR, degree=1, continuous=True) - gradphi = uw.discretisation.MeshVariable( - f"winslow_gphi_{id(mesh)}", mesh, - vtype=uw.VarType.VECTOR, degree=aux_degree, - continuous=True) - gproj = uw.systems.Vector_Projection(mesh, gradphi) - gproj.smoothing = 0.0 - _wire(gproj, elliptic=False) - _WINSLOW_CACHE[key] = ( - phi, ps, gradphi, gproj, hsolver, vol_field) - else: - phi, ps, gradphi, gproj, hsolver, vol_field = cache - - X = mesh.CoordinateSystem.X - grad_phi = sympy.Matrix( - [phi.sym[0].diff(X[i]) for i in range(cdim)]).T - Hf = hsolver.u.sym - Hmat = sympy.Matrix(cdim, cdim, - lambda i, j: Hf[i * cdim + j]) - gproj.uw_function = grad_phi - omega = float(picard_relax) - - for outer in range(n_outer): - dm = mesh.dm - tris = _tri_cells(dm) - pStart, pEnd = dm.getDepthStratum(0) - n_verts = pEnd - pStart - old_coords = np.asarray(mesh.X.coords).copy() - _cdim = mesh.cdim - - # Boundary tangential slip via the mesh-owned contract - # (boundary-slip-strategy.md): MA's natural Neumann BC (∇φ·n̂=0) makes - # ∇φ tangential at the boundary, so slip vertices slide along their - # surface (radial ring / box face / facet) and snap back; non-slip, - # junction, and degenerate-normal vertices pin. Replaces the inline - # per-ring / box-edge snap (the 'ring'/'box' hint is now inferred from - # the registered bounding surfaces). - is_pinned, _project = mesh.boundary_slip( - boundary_slip, reference_coords=old_coords, - boundary_labels=pinned_labels) - - if tris is not None and n_outer > 1: - patch = _patch_volumes(tris, old_coords, n_verts, vol_field) - patch /= float(np.mean(patch)) - else: - patch = np.ones(n_verts, dtype=np.float64) - _va = vol_field.array - _va[...] = patch.reshape(_va.shape) - - rho_t = np.asarray( - uw.function.evaluate(metric, old_coords)).reshape(-1) - b = rho_t * patch - inv_sqrt_b_mean = _global_mean(np.mean(1.0 / np.sqrt(b))) - c = 1.0 / (inv_sqrt_b_mean ** 2) - - # Target-side ρ evaluation: substitute X[i] → X[i] + - # gradphi.sym[i] so ρ is queried at the moving target - # x + ∇φ(x), not the source x. Removes the phase error - # where refinement-by-size is transported away from the - # feature location by ∇φ. gradphi.sym values are updated - # each Picard iter (gproj.solve below) so the source self- - # consistently tracks the current map estimate. - if target_side_rho: - metric_target = metric.subs( - [(X[i], X[i] + gradphi.sym[i]) - for i in range(cdim)]) - else: - metric_target = metric - g = c / (metric_target * vol_field.sym[0]) - if cdim == 2: - Hxx = Hf[0] - Hxy = (Hf[1] + Hf[2]) / 2 - Hyy = Hf[3] - f_src = sympy.sqrt( - (Hxx - Hyy) ** 2 + 4 * Hxy ** 2 + 4 * g) - 2 - else: - f_src = (g - 1.0) - Hmat.det() - ps.f = sympy.Matrix([[_EQUIDIST_SIGN * f_src]]) - - hsolver.u.array[...] = 0.0 - - zero_init_guess = not _warm_start_krylov(linear_solver) - prev_change = None - # If target-side ρ is on, gradphi needs to be tracking the - # current φ inside the Picard loop (it's used by ps.f via - # the X→X+gradphi substitution). Initialise to zero so the - # first ps.solve sees ρ at source (= identity map estimate). - if target_side_rho: - gradphi.array[...] = 0.0 - for it in range(n_picard): - phi_prev = np.asarray(phi.array).copy() - ps.solve(zero_init_guess=zero_init_guess) - phi.array[...] = ((1.0 - omega) * phi_prev - + omega * np.asarray(phi.array)) - hsolver.solve() - if target_side_rho: - gproj.solve() # update target-side ρ for next iter - change = float(np.abs( - np.asarray(phi.array) - phi_prev).max()) - change = _global_max(change) - if prev_change is not None and change < 1.0e-6: - break - prev_change = change - - if not target_side_rho: - gproj.solve() - disp = np.asarray( - uw.function.evaluate(gradphi.sym, old_coords) - ).reshape(old_coords.shape) - - # Directional move-weighting (opt-in; default None ⇒ unchanged). - if move_anisotropy is not None and cdim == 2: - disp = _reweight_displacement_radial_tangential( - disp, old_coords, move_anisotropy) - - step = _cap_step_to_edge_fraction( - relax * disp, dm, old_coords, step_frac) - - free = ~is_pinned - # _project: slip → ring (∥ only) - new_coords, scale = _backtracked_move( - old_coords, step, free, tris, _project) - - mesh._deform_mesh(new_coords) - - d = float(np.linalg.norm( - new_coords - old_coords, axis=1).max()) - if uw.mpi.size > 1: - d = _global_sum(d ** 2) ** 0.5 - if verbose: - uw.pprint( - f" equidistribute MA outer {outer+1}/{n_outer}: " - f"c={c:.4f} scale={scale:.3f} max|Δx|={d:.3e}") - if d < outer_tol: - break - - -def _ot_improvement_step(mesh, metric, pinned_labels, verbose, - n_outer=1, relax=1.0, - step_frac=0.3, - outer_tol=1.0e-4, - boundary_slip=False, - linear_solver="direct", phi_degree=2): - r"""OT-improvement step: one (or a few) weighted-Poisson - equidistribution flow iterations. - - Solves on the *current* mesh - - .. math:: - - \nabla\!\cdot(\rho\,\nabla\phi) - \;=\;-\,\rho\,\log\!\bigl(V\rho/K\bigr), - \quad K=\exp(\langle\rho\log(V\rho)\rangle/\langle\rho\rangle), - \quad \nabla\phi\cdot\hat{n}=0, - - and moves nodes by ``relax · ∇φ``. ``V_i`` is the dual patch - area at vertex ``i``; the source vanishes identically at - equidistribution ``V_i\,\rho_i\equiv K``. - - Semantics: this is a *single OT improvement step* w.r.t. the - current mesh — the input mesh has no special status (it is - whatever you currently have). Calling it again from the - deformed mesh applies another improvement step. Compose - freely with spring / smoothing / anisotropic. - - Differences from ``_monge_ampere_mover`` (the convex-branch - BFO Picard): - - * Linear: one weighted-Poisson per outer iter, no inner - Picard, no Hessian recovery, no convex-branch radical. - * The source uses the *current* mesh's patch volumes; the - formulation is identically zero at equidistribution, so - iterations are self-stabilising (no over-correction). - * ρ at the current node positions (no source-vs-target - asymmetry; the iteration is on the current mesh, ρ is at - its physical position). - - Parameters mirror ``_monge_ampere_mover`` where they apply. - ``n_outer`` composes outer improvement steps; the source - drives toward zero so the per-iter motion naturally - diminishes. - """ - import sympy - - pinned_labels = tuple(pinned_labels) - dm = mesh.dm - pStart, pEnd = dm.getDepthStratum(0) - cStart, cEnd = dm.getHeightStratum(0) - cone_size = dm.getConeSize(cStart) if cEnd > cStart else 0 - if linear_solver not in ("direct", "gamg"): - raise ValueError( - f"linear_solver must be 'direct' or 'gamg', " - f"got {linear_solver!r}") - phi_degree = int(phi_degree) - aux_degree = max(1, phi_degree - 1) - cdim = mesh.cdim - if cdim != 2: - raise NotImplementedError( - "_ot_improvement_step: 2D meshes only for now.") - - # Boundary slip uses the projected boundary-normal field - # (mesh.Gamma_P1). This is reliable only for *radial* coordinate - # systems (cylindrical / spherical / geographic), where mesh.Gamma is - # the coordinate-derived radial field and evaluates cleanly at vertices. - # For Cartesian boundaries the vertex-evaluated facet normal is - # degenerate (0/0), so we pin the boundary instead of slipping with a - # garbage normal. 'ring'/'box'/'axes' are legacy aliases for slip-on. - from underworld3.meshing._ot_adapt import _is_radial_coords as _isr - if isinstance(boundary_slip, str): - _slip_req = boundary_slip.strip().lower() in ( - "ring", "box", "axes", "axis", "true", "on", "1") - else: - _slip_req = bool(boundary_slip) - _slip_on = _slip_req and _isr(mesh) - if _slip_on: - # Create / refresh the projected normals ONCE here, before the OT - # Poisson solver's DM is built — creating the _n_proj MeshVariable - # mid-mover would stale that DM handle (project_uw3_smoother_footguns). - try: - mesh._update_projected_normals() - except Exception: - _slip_on = False - - key = (id(mesh), pinned_labels, - pEnd - pStart, cEnd - cStart, cone_size, - linear_solver, phi_degree) - - cache = _OT_CACHE.get(key) - if cache is None: - _wire = _solver_wiring(linear_solver) - phi = uw.discretisation.MeshVariable( - f"ot_phi_{id(mesh)}", mesh, - vtype=uw.VarType.SCALAR, degree=phi_degree, - continuous=True) - ps = uw.systems.Poisson(mesh, phi) - ps.constitutive_model = uw.constitutive_models.DiffusionModel - # weighted diffusion: D(x) = ρ(x). Updated each outer iter - # via the symbolic metric expression (evaluated at the - # current mesh's quad pts). - ps.constitutive_model.Parameters.diffusivity = metric - ps.constant_nullspace = True - _wire(ps, singular=True, elliptic=True) - vol_field = uw.discretisation.MeshVariable( - f"ot_vol_{id(mesh)}", mesh, - vtype=uw.VarType.SCALAR, degree=1, continuous=True) - gradphi = uw.discretisation.MeshVariable( - f"ot_gphi_{id(mesh)}", mesh, - vtype=uw.VarType.VECTOR, degree=aux_degree, - continuous=True) - gproj = uw.systems.Vector_Projection(mesh, gradphi) - gproj.smoothing = 0.0 - _wire(gproj, elliptic=False) - X = mesh.CoordinateSystem.X - gradphi_sym = sympy.Matrix( - [phi.sym[0].diff(X[i]) for i in range(cdim)]).T - gproj.uw_function = gradphi_sym - _OT_CACHE[key] = (phi, ps, gradphi, gproj, vol_field) - else: - phi, ps, gradphi, gproj, vol_field = cache - - zero_init_guess = not _warm_start_krylov(linear_solver) - - for outer in range(n_outer): - dm = mesh.dm - tris = _tri_cells(dm) - pStart, pEnd = dm.getDepthStratum(0) - n_verts = pEnd - pStart - old_coords = np.asarray(mesh.X.coords).copy() - _cdim = mesh.cdim - - # Boundary tangential slip via the mesh-owned contract - # (boundary-slip-strategy.md). Slip stays gated to radial meshes via - # ``_slip_on`` (a Cartesian boundary pins — the vertex-evaluated facet - # normal is degenerate there, see above); on a radial mesh the - # registered radial surfaces do the tangent slide + |r| restore. - is_pinned, _project = mesh.boundary_slip( - boundary_slip if _slip_on else False, - reference_coords=old_coords, boundary_labels=pinned_labels) - - # --- compute V (patch volumes) on current mesh --------- - if tris is None: - patch = np.ones(n_verts, dtype=np.float64) - else: - patch = _patch_volumes(tris, old_coords, n_verts, vol_field) - # Normalise so the mean over the domain is the cell mean. - patch_mean = _global_mean(np.mean(patch)) - # Write current V values into the MeshVariable. - _va = vol_field.array - _va[...] = (patch / max(patch_mean, 1e-30)).reshape(_va.shape) - - # --- compute K = exp(<ρ log(Vρ)> / <ρ>) ---------------- - rho_at_y = np.asarray(uw.function.evaluate( - metric, old_coords)).reshape(-1) - Vrho = (patch / max(patch_mean, 1e-30)) * rho_at_y - # weighted geometric mean (zero-mean Neumann compat - # condition) — guard against Vrho≤0: - Vrho_pos = np.clip(Vrho, 1e-30, None) - wnum = float(np.sum(rho_at_y * np.log(Vrho_pos))) - wden = float(np.sum(rho_at_y)) - wnum = _global_sum(wnum) - wden = _global_sum(wden) - ln_K = wnum / max(wden, 1e-30) - K_val = float(np.exp(ln_K)) - - # --- source: f = -ρ · log(V·ρ / K) --------------------- - # SNES_Poisson convention: F0 = -f, strong form ∇·(D∇u) - # = -ps.f. We want ∇·(ρ∇φ) = -ρ·log(V·ρ/K) ⇒ ps.f = - # ρ·log(V·ρ/K). - f_src = metric * sympy.log( - metric * vol_field.sym[0] / sympy.Float(K_val)) - ps.f = sympy.Matrix([[f_src]]) - - # --- solve weighted Poisson ---------------------------- - ps.solve(zero_init_guess=zero_init_guess) - gproj.solve() - disp = np.asarray(uw.function.evaluate( - gradphi.sym, old_coords) - ).reshape(old_coords.shape) - - step = _cap_step_to_edge_fraction( - float(relax) * disp, dm, old_coords, step_frac) - - # --- coherent global signed-area backtrack ------------- - free = ~is_pinned - new_coords, scale = _backtracked_move( - old_coords, step, free, tris, _project) - - mesh._deform_mesh(new_coords) - - d = float(np.linalg.norm( - new_coords - old_coords, axis=1).max()) - if uw.mpi.size > 1: - d = _global_sum(d ** 2) ** 0.5 - - # Per-iter "imbalance" diagnostic — std of log(V·ρ/K). - imb = float(np.std(np.log(Vrho_pos) - ln_K)) - if uw.mpi.size > 1: - imb_sq = _global_sum(imb * imb) - cnt = int(_global_sum(Vrho_pos.size)) - imb = (imb_sq / max(cnt, 1)) ** 0.5 - - if verbose: - uw.pprint( - f" OT-improve outer {outer+1}/{n_outer}: " - f"K={K_val:.4f} imb={imb:.3e} " - f"scale={scale:.3f} max|Δx|={d:.3e}") - if d < outer_tol: - break diff --git a/src/underworld3/meshing/smoothing/spring.py b/src/underworld3/meshing/smoothing/spring.py deleted file mode 100644 index e02836af..00000000 --- a/src/underworld3/meshing/smoothing/spring.py +++ /dev/null @@ -1,315 +0,0 @@ -"""The spring-equilibrium mover (``method="spring"``, the default -metric path). See the package docstring for the module map. -""" - -import warnings - -import numpy as np - -import underworld3 as uw - -from .graph import (_edge_pairs, _tri_cells, _signed_areas, - _global_sum, _global_min, _global_max) - - -# Cached spring-smoother topology state keyed by (mesh-id, -# pinned-labels, topology): the edge vertex-index pairs and per-node -# incident-edge degree. Rebuilt automatically on a topology change -# (remesh / adapt / repartition), which produces a new cache key. -_SPRING_CACHE: dict = {} - - -def _spring_equilibrium_mover(mesh, metric, pinned_labels, verbose, - max_cg_iters=300, - boundary_slip=False, shape_w=1.0, size_w=8.0, - n_sweeps=None): - r"""Metric-driven mesh grading by elastic-spring equilibrium. - - Every mesh edge is a linear spring whose *rest length* is set - from the target density, - - .. math:: - - L^0_{ij} \;\propto\; \rho_{\mathrm{tgt}}^{-1/d}, - - scaled once so the total rest length equals the total current - edge length (overall scale preserved — pure redistribution). - The interior nodes are moved to the **mechanical equilibrium** - by *minimising the truss energy* - - .. math:: - - E(\mathbf{x}) \;=\; \tfrac12 \sum_{e} - \big(\,|\mathbf{x}_i-\mathbf{x}_j| - L^0_e\,\big)^2 - - over the free (non-pinned) nodes with **nonlinear conjugate - gradients** (Polak–Ribière⁺) and an Armijo line search whose - trial step is rejected if any cell would invert. Solving the - equilibrium — rather than creeping with damped Jacobi sweeps, - which stall against a per-sweep global tangle freeze — is what - lets the absolute rest-length target actually grade the mesh - toward spacing ``∝ ρ_tgt^{-1/d}``. - - ``ρ_tgt`` is Lagrangian (``metric = f(r0)`` with ``r0`` a frozen - mesh variable), so the rest lengths are fixed per material node - (computed once) and the *design* grading is restored even after - the mesh deformed. Uniform ``ρ_tgt`` ⇒ all rest lengths equal - the mean edge length ⇒ only a benign mild regularisation toward - uniform spacing (no grading change). - - ``max_cg_iters`` caps the CG iterations (CG converges far faster - than the old Jacobi sweep budget); ``n_sweeps`` is its deprecated - former name, accepted for one cycle with a DeprecationWarning. - The old ``relax`` / ``step_frac`` parameters were unused on the - equilibrium path (the CG line search controls the step and the - inversion guard) and have been removed. ``n_iters`` / ``alpha`` - do not apply. - """ - if n_sweeps is not None: - warnings.warn( - "the 'n_sweeps' argument is renamed; use max_cg_iters= " - "(it caps the nonlinear-CG iterations, not Jacobi sweeps)", - DeprecationWarning, stacklevel=2) - max_cg_iters = n_sweeps - pinned_labels = tuple(pinned_labels) - dm = mesh.dm - pStart, pEnd = dm.getDepthStratum(0) - cStart, cEnd = dm.getHeightStratum(0) - cone_size = dm.getConeSize(cStart) if cEnd > cStart else 0 - n_verts = pEnd - pStart - key = (id(mesh), pinned_labels, - n_verts, cEnd - cStart, cone_size) - - cache = _SPRING_CACHE.get(key) - if cache is None: - edges = _edge_pairs(dm) - if edges.shape[0] == 0: - return - deg = np.bincount( - edges.ravel(), minlength=n_verts).astype(np.float64) - deg[deg == 0.0] = 1.0 - _SPRING_CACHE[key] = (edges, deg) - else: - edges, deg = cache - - tris = _tri_cells(dm) - cdim = mesh.cdim - v0 = edges[:, 0] - v1 = edges[:, 1] - - coords = np.asarray(mesh.X.coords, dtype=np.float64).copy() - - # Boundary tangential slip via the mesh-owned contract - # (boundary-slip-strategy.md): each slip vertex slides tangentially and - # snaps back onto its bounding surface (radial ring / plane / facet); - # non-slip, junction, and degenerate-normal vertices pin. Replaces the - # per-ring COM radial snap (one node/ring anchored the rotation gauge); - # the global inversion guard below still blocks a slip node overtaking a - # neighbour, and tangential θ-drift is a harmless re-parameterisation. - is_pinned, _project = mesh.boundary_slip( - boundary_slip, reference_coords=coords, boundary_labels=pinned_labels) - - free = ~is_pinned - - # ===== Volumetric spring network (shape ⟂ size, decoupled) ==== - # EQUAL edge springs (uniform rest length L̄ = current mean - # edge) are a pure SHAPE regulariser → equant cells, resists - # the slivers/degeneracy the graded-edge form produced. The - # SIZE grading lives entirely in a per-CELL area ("volumetric") - # constraint: each triangle's area is driven to a target - # A0 ∝ 1/ρ_tgt (scaled so ΣA0 = Σ(initial area) ⇒ total area - # conserved, pure redistribution). Both energy terms are - # written as *relative* squared errors so the shape/size - # weights (shape_w, size_w) are pure dimensionless knobs. - e_vec = coords[v1] - coords[v0] - L_cur = np.linalg.norm(e_vec, axis=1) - sum_L = float(L_cur.sum()) - n_e = float(L_cur.size) - sum_L = _global_sum(sum_L) - n_e = _global_sum(n_e) - Lbar = sum_L / max(n_e, 1.0) # uniform edge rest length - L0 = np.full_like(L_cur, Lbar) - L0_mean = Lbar - - # Per-cell target area from ρ_tgt at the (initial) centroid. - # Lagrangian metric ⇒ computed ONCE (rides material points). - if tris is not None: - ca = coords[tris[:, 0]] - cb = coords[tris[:, 1]] - cc = coords[tris[:, 2]] - cent = (ca + cb + cc) / 3.0 - rho_c = np.asarray( - uw.function.evaluate(metric, cent)).reshape(-1) - rho_c = np.maximum(rho_c, 1.0e-30) - a_init = np.abs(_signed_areas(coords, tris)) - inv = 1.0 / rho_c - sA = float(a_init.sum()) - sI = float(inv.sum()) - sA = _global_sum(sA) - sI = _global_sum(sI) - A0 = (sA / max(sI, 1.0e-30)) * inv # ΣA0 = Σa_init - A0 = np.maximum(A0, 1.0e-30) - ti0, ti1, ti2 = tris[:, 0], tris[:, 1], tris[:, 2] - else: - A0 = None - - # ---- Solve the truss EQUILIBRIUM, not Jacobi creep ---------- - # Minimise the spring energy E(x) = ½ Σ_e (|x_i−x_j| − L0_e)² - # over the interior nodes (boundary pinned) by nonlinear - # conjugate gradients (Polak–Ribière⁺) with an Armijo line - # search whose trial step is REJECTED if any cell would invert - # — the tangle guard is inside the optimiser, so it converges to - # the true equilibrium instead of stalling against a per-sweep - # global freeze (the Jacobi relaxation's failure mode). - free_idx = np.nonzero(free)[0] - n_free = int(free_idx.size) - # TODO(BUG): rank-LOCAL early return before the collective CG loop. - # Under MPI, a rank whose local vertices are all pinned returns here - # while the other ranks proceed into _global_sum/_global_min - # collectives -> deadlock. Consistent with the documented - # serial-exact status of this path, but a latent parallel hazard; - # the exit decision should be reduced globally (as the OT/MA movers - # do) before this mover is promoted to parallel-exact. Tracked in - # issue #346 (mover retired in favour of method='mmpde'; not being - # hardened). NOTE: any collective-exit fix must also guard the - # rank-local reductions below for n_free == 0 — e.g. the dmax - # line-search seed takes .max() over d[free_idx], which raises on a - # zero-size array once a starved rank is made to participate. - if n_free == 0: - mesh._deform_mesh(coords) - return - - if tris is not None: - orient = np.sign(np.median(_signed_areas(coords, tris))) - orient = orient if orient != 0.0 else 1.0 - - def _feasible(X): - if tris is None: - return True - amin = _global_min((_signed_areas(X, tris) * orient).min()) - return amin > 0.0 - - have_area = (A0 is not None) and (cdim == 2) - - def _tri_signed(X): - a, b, c = X[ti0], X[ti1], X[ti2] - return 0.5 * ((b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - - (b[:, 1] - a[:, 1]) * (c[:, 0] - a[:, 0])) - - def _energy(X): - ev = X[v1] - X[v0] - L = np.sqrt((ev * ev).sum(axis=1)) - re = (L - Lbar) / Lbar # relative edge error - E = shape_w * _global_sum((re * re).sum()) - if have_area: - area = orient * _tri_signed(X) - ra = (area - A0) / A0 # relative area error - E += size_w * _global_sum((ra * ra).sum()) - return E - - def _energy_grad(X): - ev = X[v1] - X[v0] - L = np.sqrt((ev * ev).sum(axis=1)) - Ls = np.maximum(L, 1.0e-30) - re = (L - Lbar) / Lbar - E = shape_w * _global_sum((re * re).sum()) - G = np.zeros_like(X) - # equal-spring shape term: 2·shape_w·re/(Lbar·L)·ev - ce = (2.0 * shape_w * re / (Lbar * Ls))[:, None] - np.add.at(G, v1, ce * ev) - np.add.at(G, v0, -ce * ev) - if have_area: - a, b, c = X[ti0], X[ti1], X[ti2] - S = 0.5 * ((b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - - (b[:, 1] - a[:, 1]) * (c[:, 0] - a[:, 0])) - area = orient * S - ra = (area - A0) / A0 - E += size_w * _global_sum((ra * ra).sum()) - # ∂(area)/∂· = orient · ∂S/∂· (signed-area vertex grads) - fac = (2.0 * size_w * ra / A0 * orient)[:, None] - gA = np.empty_like(a) - gB = np.empty_like(a) - gC = np.empty_like(a) - gA[:, 0] = 0.5 * (b[:, 1] - c[:, 1]) - gA[:, 1] = 0.5 * (c[:, 0] - b[:, 0]) - gB[:, 0] = 0.5 * (c[:, 1] - a[:, 1]) - gB[:, 1] = 0.5 * (a[:, 0] - c[:, 0]) - gC[:, 0] = 0.5 * (a[:, 1] - b[:, 1]) - gC[:, 1] = 0.5 * (b[:, 0] - a[:, 0]) - np.add.at(G, ti0, fac * gA) - np.add.at(G, ti1, fac * gB) - np.add.at(G, ti2, fac * gC) - G[~free] = 0.0 - return E, G - - # Jacobi (diagonal) preconditioner: the truss Hessian is - # graph-Laplacian-structured (cond ~ (1/h)²), so plain CG crawls - # for fine meshes. M⁻¹ = diag(1/deg) — the Laplacian diagonal - # scale, free here since `deg` is already cached — clusters the - # spectrum and gives the order-of-magnitude convergence speed-up - # that turns "stuck at ~1.04" into the true graded minimum. - invdeg = (1.0 / deg)[:, None] - - X = _project(coords.copy()) - E, G = _energy_grad(X) - g0 = max(_global_sum((G * G).sum()) ** 0.5, 1.0e-30) - r = -G - s = r * invdeg - s[~free] = 0.0 - d = s.copy() - delta_new = _global_sum((r * s).sum()) - dmax = _global_max(max(float(np.linalg.norm( - d[free_idx], axis=1).max()), 1.0e-30)) - t0 = 0.5 * L0_mean / dmax - c_arm = 1.0e-4 - max_iter = int(max_cg_iters) - for it in range(max_iter): - gnorm = _global_sum((G * G).sum()) ** 0.5 - if gnorm <= 1.0e-8 * g0: - break - slope = _global_sum((G * d).sum()) # = −(r·d) - if slope >= 0.0: # not descent → restart - d = s.copy() - slope = _global_sum((G * d).sum()) - if slope >= 0.0: - break - t = t0 - accepted = False - for _ls in range(50): - Xt = X.copy() - Xt[free_idx] += t * d[free_idx] - Xt = _project(Xt) # slip nodes → boundary - if _feasible(Xt): - Et = _energy(Xt) - if Et <= E + c_arm * t * slope: - accepted = True - break - t *= 0.5 - if not accepted: - break # at equilibrium / stuck - Et, Gt = _energy_grad(Xt) - r_new = -Gt - s_new = r_new * invdeg - s_new[~free] = 0.0 - delta_old = delta_new - delta_mid = _global_sum((r_new * s).sum()) - delta_new = _global_sum((r_new * s_new).sum()) - beta = max(0.0, (delta_new - delta_mid) - / max(delta_old, 1.0e-30)) # preconditioned PR⁺ - X, E, G = Xt, Et, Gt - d = s_new + beta * d - s = s_new - t0 = min(2.0 * t, 100.0 * t0) # grow but stay sane - - if verbose and (it % 25 == 0 or it == max_iter - 1): - ev = X[v1] - X[v0] - L = np.sqrt((ev * ev).sum(axis=1)) - rms = (_global_sum(((L - L0) ** 2).sum()) - / max(_global_sum(L0.size), 1.0)) ** 0.5 - uw.pprint( - f" spring PCG iter {it+1}/{max_iter}: " - f"E={E:.4e} rms(L-L0)/L0=" - f"{rms / max(L0_mean, 1e-30):.3e} |g|={gnorm:.2e}") - - coords = X - mesh._deform_mesh(coords) diff --git a/src/underworld3/meshing/surfaces.py b/src/underworld3/meshing/surfaces.py index dbbf6a61..d906cd61 100644 --- a/src/underworld3/meshing/surfaces.py +++ b/src/underworld3/meshing/surfaces.py @@ -2707,8 +2707,7 @@ def fault_metric_tensor(mesh, faults, refinement=3.0, width="auto", base=1.0): M = uw.meshing.fault_metric_tensor(mesh, faults, refinement=3.0) uw.meshing.smooth_mesh_interior( - mesh, metric=M, method="anisotropic", boundary_slip=False, - method_kwargs=dict(n_outer=12, relax=0.4)) + mesh, metric=M, method="mmpde", boundary_slip=False) Construction — summed over every fault segment ``i`` (normal ``n_i``, point-to-segment distance ``d_i(x)``): @@ -2754,7 +2753,7 @@ def fault_metric_tensor(mesh, faults, refinement=3.0, width="auto", base=1.0): sympy.Matrix The ``2×2`` analytic metric tensor ``M(x)`` (a function of ``mesh.CoordinateSystem.X``), to pass as ``metric=`` with - ``method="anisotropic"``. + ``method="mmpde"``. """ cdim = mesh.cdim if cdim != 2: @@ -2814,20 +2813,16 @@ def fault_comb_metric(mesh, faults, cell_size, n_across=4, amplitude=6.0, tooth_width=None, combine="sum"): r"""Build a scalar **comb** metric ``ρ(x)`` that refines a band of a controlled number of roughly-**uniform** cells *across* one or more faults, - for the isotropic equidistribution mover (``method="ma"``). + for the scalar-metric equidistribution mover (``method="mmpde"``). Pass the result straight to the mover:: rho = uw.meshing.fault_comb_metric(mesh, faults, cell_size=0.006, n_across=4) - uw.meshing.smooth_mesh_interior( - mesh, metric=rho, method="ma", - method_kwargs=dict(n_outer=1, n_picard=25)) # single-shot + uw.meshing.smooth_mesh_interior(mesh, metric=rho, method="mmpde") - Use the **single-shot** map (``n_outer=1``): one Caffarelli-clean - Monge–Ampère solve, untangled by construction (no folding), with no - outer-iteration compounding and nothing to tune — the most robust - configuration, and the comb's teeth give the single map all the row + The MMPDE map is non-folding by construction, and the comb's teeth + give it all the row structure it needs (~``n_across``−1 even layers, centred). ``n_outer=2`` realises a touch more of the requested ``n_across`` (the single map is mildly node-budget-capped) at ~1.6× the cost; rarely needed. @@ -2855,9 +2850,9 @@ def fault_comb_metric(mesh, faults, cell_size, n_across=4, amplitude=6.0, Parameters ---------- mesh : Mesh - 2D mesh. (The isotropic equidistribution movers — ``ma``/``ot`` — are - 2D-only; 3D would need a 3D equidistribution mover, which does not yet - exist, so this builder is 2D-only.) + 2D mesh. (The MMPDE mover is 2D-only; 3D would need a 3D + discretization of the mover, which does not yet exist, so this + builder is 2D-only.) faults : Surface | array | list Fault geometry in mesh coordinate space — a :class:`Surface`, an ``(N>=2, 2|3)`` polyline array, or a list mixing those. Each fault's @@ -2885,13 +2880,13 @@ def fault_comb_metric(mesh, faults, cell_size, n_across=4, amplitude=6.0, ------- sympy.Expr The scalar comb metric ``ρ(x)``, to pass as ``metric=`` with - ``method="ma"``. + ``method="mmpde"``. """ cdim = mesh.cdim if cdim != 2: raise NotImplementedError( - "fault_comb_metric is 2D only (the isotropic equidistribution " - "movers are 2D; 3D needs a 3D equidistribution mover)") + "fault_comb_metric is 2D only (the MMPDE mover is 2D; 3D needs " + "a 3D discretization of the mover)") dx = float(cell_size) if not (dx > 0.0): raise ValueError(f"cell_size must be positive; got {cell_size}") @@ -2927,7 +2922,7 @@ def fault_comb_metric(mesh, faults, cell_size, n_across=4, amplitude=6.0, def compose_metrics(metrics, compose="max"): r"""Combine several scalar density metrics into one, for the - equidistribution mover (``method="ma"``). + scalar-metric equidistribution mover (``method="mmpde"``). Each item may be either a metric (a scalar sympy expression or MeshVariable) or a ``(metric, weight)`` tuple. The default ``"max"`` @@ -2952,8 +2947,7 @@ def compose_metrics(metrics, compose="max"): metric_choice="arc-length") rho_F = uw.meshing.fault_comb_metric(mesh, faults, cell_size=0.008) rho = uw.meshing.compose_metrics([(rho_T, 1.0), (rho_F, 3.0)]) # fault heavier - uw.meshing.smooth_mesh_interior(mesh, metric=rho, method="ma", - method_kwargs=dict(n_outer=1, n_picard=25)) + uw.meshing.smooth_mesh_interior(mesh, metric=rho, method="mmpde") Parameters ---------- @@ -3057,17 +3051,20 @@ def fault_metric(mesh, faults, method="ma", *, cell_size, ``n_across`` elements of size ``cell_size`` across a band around the fault(s)*. - The three movers consume **different metric objects with different + The consumers take **different metric objects with different semantics**, so this facade unifies the *intent* and emits the right - representation — it does not pretend they are interchangeable: + representation — it does not pretend they are interchangeable (the + ``method`` values name the metric KIND; the historical mover names + are kept as spellings even though the movers themselves were + superseded by MMPDE in 2026-07): =================== ============================ =========================== ``method`` returns pass to =================== ============================ =========================== ``"ma"`` (default) scalar comb density (sympy) ``smooth_mesh_interior( - method="ma")`` + method="mmpde")`` ``"anisotropic"`` 2×2 tensor (sympy Matrix) ``smooth_mesh_interior( - method="anisotropic")`` + method="mmpde")`` ``"adapt"``/``"mmg"`` ``h⁻²`` MeshVariable ``mesh.remesh(...)`` =================== ============================ =========================== @@ -3116,8 +3113,7 @@ def fault_metric(mesh, faults, method="ma", *, cell_size, # uniform-ish band, fixed topology (the slip-rheology recipe) rho = uw.meshing.fault_metric(mesh, faults, method="ma", cell_size=0.006, n_across=4) - uw.meshing.smooth_mesh_interior(mesh, metric=rho, method="ma", - method_kwargs=dict(n_outer=1, n_picard=25)) + uw.meshing.smooth_mesh_interior(mesh, metric=rho, method="mmpde") """ if cell_size is None or not (float(cell_size) > 0.0): raise ValueError("cell_size must be a positive number") diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 12563b16..e198aad5 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1883,7 +1883,7 @@ def on_remesh(self, ctx): mesh motion when reading the CARRY'd history. One-step pulse: the next solve consumes Δx and clears it. - * **Opt-out (e.g. OT_adapt reset).** When the adapt is a + * **Opt-out (discrete-jump adapts).** When the adapt is a discrete jump rather than a smooth displacement (``ctx.scratch.get("ale_opt_out")``), the linear ``Δx/dt → v_mesh`` interpretation breaks down. Fall back to diff --git a/tests/test_0641_wave_c_api_shims.py b/tests/test_0641_wave_c_api_shims.py index de9c9103..5ceb68ef 100644 --- a/tests/test_0641_wave_c_api_shims.py +++ b/tests/test_0641_wave_c_api_shims.py @@ -1,4 +1,5 @@ -"""Wave C API-harmonization contract tests (WC-01..WC-13). +"""Wave C API-harmonization contract tests (WC-01..WC-12; the WC-13 +spring-mover shim tests were removed with the 2026-07 mover retirement). Every deprecation shim introduced by the 2026-07 Wave C API harmonization carries the two-test contract here: @@ -415,35 +416,3 @@ def test_pack_sync_warns_once(self, swarm_var): assert _one_deprecation(rec) == 1 with _no_deprecation(): swarm_var.pack_uw_data_to_petsc(data) - - -# --------------------------------------------------------------------------- -# WC-13 - smoothing: dead params dropped, n_sweeps renamed max_cg_iters -# --------------------------------------------------------------------------- - -class TestSpringMoverSignature: - def test_dead_params_gone_new_name_present(self): - import inspect - - from underworld3.meshing.smoothing import _spring_equilibrium_mover - - params = inspect.signature(_spring_equilibrium_mover).parameters - assert "relax" not in params - assert "step_frac" not in params - assert "max_cg_iters" in params - assert "n_sweeps" in params # deprecated alias, one cycle - - def test_n_sweeps_alias_warns_once(self): - import numpy as np - - from underworld3.meshing.smoothing import _spring_equilibrium_mover - - tri_mesh = uw.meshing.UnstructuredSimplexBox( - minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5 - ) - pinned = ("Top", "Bottom", "Left", "Right") - with pytest.warns(DeprecationWarning, match="max_cg_iters") as rec: - _spring_equilibrium_mover(tri_mesh, sympy.Integer(1), pinned, False, n_sweeps=1) - assert _one_deprecation(rec) == 1 - with _no_deprecation(): - _spring_equilibrium_mover(tri_mesh, sympy.Integer(1), pinned, False, max_cg_iters=1) diff --git a/tests/test_0750_meshing_follow_metric.py b/tests/test_0750_meshing_follow_metric.py index f4008a26..104c63a4 100644 --- a/tests/test_0750_meshing_follow_metric.py +++ b/tests/test_0750_meshing_follow_metric.py @@ -17,11 +17,14 @@ The mover's achieved cell-size envelope is APPROXIMATE — anisotropic cells + iterative deformation map mean the eigenvalue clamp doesn't -literally bound the achieved h_max on a per-cell basis. Validation +literally bound the achieved h_max on a per-cell basis (since the +2026-07 mover retirement the backing mover is the variational MMPDE +mover with theta=0.5, pure equidistribution weighting). Validation on a sharp-tanh annulus shows: -* refinement side: achieved h_min within ~5-10% of h0/refinement -* coarsening side: achieved h_max within ~2× of h0·coarsening +* refinement side: achieved h_min within ~30-40% of h0/refinement + (platform-dependent: ~27% on macOS, ~40% on Linux CI at refinement=2) +* coarsening side: achieved h_max within ~2x of h0*coarsening The tests here use loose tolerances reflecting that empirical reality. """ @@ -189,17 +192,18 @@ def test_follow_metric_refinement_envelope_approximate(): m, T, refinement=ref, skip_threshold=None) h_min_cell, _, _ = _cell_h_stats(m) target_h_min = h0 / ref - # The spring keeps mean-edge min at or close to target. - # Allow up to ~25% over-spec (under-refinement, the safe - # side) and ~15% under-spec (over-refinement, the unsafe - # side). + # The MMPDE mover honours the envelope approximately (it + # equidistributes toward the metric, it does not clamp cells). + # Measured spread at refinement=2: 1.27 (macOS) - 1.40 (Linux + # CI) over-spec; the bounds catch over-refinement (unsafe) and + # gross under-refinement, not the platform jitter. assert h_min_cell / target_h_min > 0.85, ( f"refinement={ref}: mean-edge h_min/target = " f"{h_min_cell/target_h_min:.3f}, want > 0.85 " f"(should not over-refine past spec)") - assert h_min_cell / target_h_min < 1.30, ( + assert h_min_cell / target_h_min < 1.5, ( f"refinement={ref}: mean-edge h_min/target = " - f"{h_min_cell/target_h_min:.3f}, want < 1.30 " + f"{h_min_cell/target_h_min:.3f}, want < 1.5 " f"(should not under-refine far past spec)") @@ -266,17 +270,9 @@ def test_follow_metric_no_slivers_after_adaptive_polish(): @pytest.mark.tier_a @pytest.mark.level_1 -@pytest.mark.xfail( - reason="dev↔#202 (manifold/extract_surface) merge reconciliation. The " - "ADAPTATION is bit-identical to development (moved coords match to 8 " - "figures; evaluate is bit-identical interior + near-boundary). The only " - "difference is the FIELD re-interpolation during the remesh step near the " - "curved boundary (T-after differs ~1.8%), which nudges the second call's " - "mismatch just over skip_threshold=0.9 so it re-adapts instead of skipping. " - "A skip-optimisation only — adaptation correctness is preserved. Track in " - "the #202 nav/remesh field-transfer reconciliation.", - strict=False) def test_follow_metric_skip_threshold_skips_aligned_mesh(): + # (An xfail marker for a #202-era field-transfer nudge was removed + # 2026-07: the MMPDE-backed follow_metric passes this cleanly.) """A mesh that's already aligned (here: any uniform mesh with the field still in the natural shape) gets the skip-on-align short-circuit when skip_threshold is permissive.""" diff --git a/tests/test_0760_mesh_ot_adapt.py b/tests/test_0760_mesh_ot_adapt.py deleted file mode 100644 index c0b09a09..00000000 --- a/tests/test_0760_mesh_ot_adapt.py +++ /dev/null @@ -1,242 +0,0 @@ -"""Locks the mesh.OT_adapt() public API. - -OT_adapt() runs the validated optimal-transport reset adapt as a method -on the mesh: reset to a cached reference, FE-remap the driving field onto -the clean canvas, build a gradient metric, run the OT mover, then FE-remap -fields and zero a cold-restart set. These tests pin: - -* the mesh actually moves and the driving field's spatial pattern survives - the FE-remap (Annulus and Box); -* boundary behaviour — radial (Annulus) boundary nodes slide tangentially - but stay on the circle; Cartesian (Box) boundaries are pinned (their - vertex-evaluated normal is degenerate, so we do not slip with it); -* the lazy reference-coordinate cache and OT_adapt_reset_reference(); -* skip_threshold short-circuits an already-aligned mesh; -* constrained-manifold meshes (dim != cdim) raise NotImplementedError. - -Behaviour is validated, not bit-for-bit: the unified Gamma_N-based slip -is a different algorithm from the historical hardcoded box/ring slip. -""" -import numpy as np -import pytest - -import underworld3 as uw -from underworld3.meshing import smoothing as _sm - - -def _annulus_with_field(slope=15.0, r_bl=0.75): - """Annulus + a tanh boundary-layer T field at radius r_bl (set - directly on the DOFs — no projection solve needed).""" - m = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, - cellSize=0.08, qdegree=3) - T = uw.discretisation.MeshVariable( - "T", m, vtype=uw.VarType.SCALAR, degree=3, continuous=True) - feat = lambda c: np.tanh((np.linalg.norm(c, axis=1) - r_bl) * slope) - T.data[:, 0] = feat(np.asarray(T.coords)) - return m, T, feat - - -def _box_with_field(slope=15.0, y_bl=0.5): - m = uw.meshing.UnstructuredSimplexBox( - minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), - cellSize=0.08, qdegree=3) - T = uw.discretisation.MeshVariable( - "Tb", m, vtype=uw.VarType.SCALAR, degree=3, continuous=True) - feat = lambda c: np.tanh((c[:, 1] - y_bl) * slope) - T.data[:, 0] = feat(np.asarray(T.coords)) - return m, T, feat - - -def _boundary_mask(mesh): - from underworld3.meshing.smoothing import ( - _pinned_mask, _auto_pinned_labels) - return _pinned_mask(mesh.dm, tuple(_auto_pinned_labels(mesh))) - - -# --------------------------------------------------------------------------- -# Annulus: moves, preserves field, boundary stays on the circle -# --------------------------------------------------------------------------- -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_moves_mesh_annulus(): - m, T, _ = _annulus_with_field() - X0 = np.asarray(m.X.coords).copy() - moved = m.OT_adapt(T, refinement=3.0, fields_to_remap=[T]) - assert moved is True - X1 = np.asarray(m.X.coords) - # interior nodes moved appreciably - assert float(np.linalg.norm(X1 - X0, axis=1).max()) > 1.0e-3 - - -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_annulus_boundary_stays_on_circle(): - m, T, _ = _annulus_with_field() - is_bnd = _boundary_mask(m) - X0 = np.asarray(m.X.coords).copy() - r0 = np.linalg.norm(X0[is_bnd], axis=1) - m.OT_adapt(T, refinement=3.0, fields_to_remap=[T]) - X1 = np.asarray(m.X.coords) - r1 = np.linalg.norm(X1[is_bnd], axis=1) - # boundary nodes slid tangentially ... - bnd_disp = float(np.linalg.norm(X1[is_bnd] - X0[is_bnd], axis=1).max()) - assert bnd_disp > 1.0e-4 - # ... but stayed on their original radii (radial drift ~ rounding) - assert float(np.abs(r1 - r0).max()) < 1.0e-3 - - -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_preserves_field_pattern_annulus(): - m, T, feat = _annulus_with_field() - m.OT_adapt(T, refinement=3.0, fields_to_remap=[T]) - # T at the adapted DOFs should match the analytic feature there - err = np.abs(np.asarray(T.data)[:, 0] - feat(np.asarray(T.coords))).max() - assert float(err) < 5.0e-2 - - -# --------------------------------------------------------------------------- -# Box: interior moves, Cartesian boundary pinned, field preserved -# --------------------------------------------------------------------------- -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_box_moves_interior_slides_boundary_on_faces(): - # Generic topology-based slip is now active on Cartesian boxes too: - # boundary nodes on the box faces SLIDE tangentially; corners are pinned - # (incident face normals disagree); the box shape (axis-aligned bounding - # planes) is preserved exactly. - m, T, feat = _box_with_field() - is_bnd = _boundary_mask(m) - X0 = np.asarray(m.X.coords).copy() - moved = m.OT_adapt(T, refinement=3.0, fields_to_remap=[T]) - assert moved is True - X1 = np.asarray(m.X.coords) - # interior refined - assert float(np.linalg.norm(X1[~is_bnd] - X0[~is_bnd], axis=1).max()) > 1.0e-3 - # box CORNERS are pinned (each of the four unit-square corners is still - # present at its original position) - for cc in ([0., 0.], [0., 1.], [1., 0.], [1., 1.]): - assert np.any(np.all(np.abs(X1 - np.array(cc)) < 1.0e-9, axis=1)), \ - f"corner {cc} lost" - # face vertices stay on the bounding planes (x=0, x=1, y=0 or y=1): - # for each originally-boundary node, at least one coord is 0 or 1 - on_face = (np.minimum(X1[is_bnd], 1 - X1[is_bnd]).min(axis=1) < 1.0e-9) - assert on_face.all(), "boundary node left the box face" - # field pattern preserved within FE-remap tolerance — looser than the - # old pinned-box test (5e-2) because face slip reallocates some node - # budget from the BL interior to the boundary, mildly increasing remap - # error on a sharp tanh feature at this coarse base mesh. Still - # well-controlled (≪ the BL amplitude of 1). - err = np.abs(np.asarray(T.data)[:, 0] - feat(np.asarray(T.coords))).max() - assert float(err) < 1.5e-1 - - -# --------------------------------------------------------------------------- -# Reference-coordinate cache + reset -# --------------------------------------------------------------------------- -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_reference_cache_and_reset(): - m, T, _ = _annulus_with_field() - X0 = np.asarray(m.X.coords).copy() - # first call snapshots the (pristine) coords as the reset reference - m.OT_adapt(T, refinement=3.0, fields_to_remap=[T]) - assert np.allclose(m._ot_adapt_reference_coords, X0) - # the mesh itself has moved away from the reference - assert not np.allclose(np.asarray(m.X.coords), X0) - # reset_reference(None) re-baselines to the current (moved) coords - m.OT_adapt_reset_reference() - assert np.allclose(m._ot_adapt_reference_coords, - np.asarray(m.X.coords)) - # explicit coords override - m.OT_adapt_reset_reference(coords=X0) - assert np.allclose(m._ot_adapt_reference_coords, X0) - - -# --------------------------------------------------------------------------- -# skip_threshold short-circuits an already-aligned mesh -# --------------------------------------------------------------------------- -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_skip_threshold_returns_false(): - m, T, _ = _annulus_with_field() - X0 = np.asarray(m.X.coords).copy() - # misalignment is in [0, 1]; threshold > 1 always skips - moved = m.OT_adapt(T, refinement=3.0, fields_to_remap=[T], - skip_threshold=2.0) - assert moved is False - # mesh untouched - assert np.allclose(np.asarray(m.X.coords), X0) - - -# --------------------------------------------------------------------------- -# Constrained-manifold mesh (dim != cdim) is not supported -# --------------------------------------------------------------------------- -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_manifold_raises(): - try: - sm = uw.meshing.SegmentedSphericalSurface2D( - radius=1.0, cellSize=0.2, numSegments=6) - except Exception: - pytest.skip("manifold surface mesh not constructible in this build") - assert sm.dim != sm.cdim - T = uw.discretisation.MeshVariable( - "Tm", sm, vtype=uw.VarType.SCALAR, degree=1, continuous=True) - with pytest.raises(NotImplementedError): - sm.OT_adapt(T) - - -# --------------------------------------------------------------------------- -# grad_smoothing_length default ("auto" ≈ grid size) de-noises the metric so -# production refinement stays sliver-free; "off" (None) does not. -# --------------------------------------------------------------------------- -def _qmin(mesh): - tris = _sm._tri_cells(mesh.dm) - X = np.asarray(mesh.X.coords); p = X[tris] - e = (np.linalg.norm(p[:, 1]-p[:, 0], axis=1)**2 - + np.linalg.norm(p[:, 2]-p[:, 1], axis=1)**2 - + np.linalg.norm(p[:, 0]-p[:, 2], axis=1)**2) - A = _sm._signed_areas(X, tris) - return float((4*np.sqrt(3)*np.abs(A)/(e+1e-30)).min()) - - -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_auto_smoothing_default_avoids_slivers(): - # Sharp boundary layer: with smoothing OFF the metric chases sub-cell - # noise and slivers form; the "auto" default de-noises it. - m_off, T_off, _ = _annulus_with_field(slope=28.0) - m_off.OT_adapt(T_off, refinement=3.0, fields_to_remap=[T_off], - grad_smoothing_length=None) - q_off = _qmin(m_off) - - m_auto, T_auto, _ = _annulus_with_field(slope=28.0) - m_auto.OT_adapt(T_auto, refinement=3.0, fields_to_remap=[T_auto]) # default "auto" - q_auto = _qmin(m_auto) - - assert q_auto > q_off # auto is strictly better - assert q_auto > 0.2 # and sliver-free - - -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_smoothing_accepts_explicit_and_unit_aware_length(): - # explicit non-dimensional length - m, T, _ = _annulus_with_field() - assert m.OT_adapt(T, refinement=3.0, fields_to_remap=[T], - grad_smoothing_length=0.0625) is True - # unit-aware (Pint) length is accepted (routed through the projection's - # non-dimensionalisation) - m2, T2, _ = _annulus_with_field() - L = 0.0625 * uw.scaling.units.meter - assert m2.OT_adapt(T2, refinement=3.0, fields_to_remap=[T2], - grad_smoothing_length=L) is True - - -@pytest.mark.tier_a -@pytest.mark.level_1 -def test_ot_adapt_invalid_smoothing_string_raises(): - m, T, _ = _annulus_with_field() - with pytest.raises(ValueError): - m.OT_adapt(T, refinement=3.0, grad_smoothing_length="bogus") diff --git a/tests/test_0762_bounding_surfaces.py b/tests/test_0762_bounding_surfaces.py index e08eb4e8..fe16d308 100644 --- a/tests/test_0762_bounding_surfaces.py +++ b/tests/test_0762_bounding_surfaces.py @@ -144,7 +144,7 @@ def test_boundary_slip_facet_fallback_when_no_surface_registered(): # mesh.boundary_slip builds a transient `facet` surface from the reference # facets, so the vertices slip along the boundary polygon (the same path a # mesh loaded from file takes). See boundary-slip-strategy.md. - from underworld3.meshing._ot_adapt import ( + from underworld3.meshing.smoothing import ( _boundary_facets, _nearest_on_facets_2d) m = _annulus() m.bounding_surfaces.clear() # remove the analytic surfaces diff --git a/tests/test_0762_fault_metric_tensor.py b/tests/test_0762_fault_metric_tensor.py index aeff4191..605b57be 100644 --- a/tests/test_0762_fault_metric_tensor.py +++ b/tests/test_0762_fault_metric_tensor.py @@ -6,7 +6,7 @@ * tensor structure — at a fault the across-fault eigenvalue is base*R^2 and the along-fault eigenvalue is base; far away M = base*I; -* the supplied-tensor mover (smooth_mesh_interior method="anisotropic", +* the supplied-tensor mover (smooth_mesh_interior method="mmpde", metric=M) centres TWO close faults on their lines (|offset| < one refined cell) with the topology preserved (r-adapt, not h-adapt); * the builder accepts Surface objects equivalently to raw segments; @@ -157,8 +157,11 @@ def test_fault_comb_metric_two_faults_band(): n0 = len(np.asarray(m.X.coords)); nc0 = len(_tri_cells(m.dm)) dx = 0.006 rho = uw.meshing.fault_comb_metric(m, _SEG3, cell_size=dx, n_across=4) - _sm.smooth_mesh_interior(m, metric=rho, method="ma", boundary_slip=False, - method_kwargs=dict(n_outer=1, n_picard=25)) + # p=2: sharper equidistribution power — resolves the comb's + # sub-cell teeth to the banded target the retired MA mover reached. + _sm.smooth_mesh_interior(m, metric=rho, method="mmpde", + boundary_slip=False, + method_kwargs=dict(p=2.0)) Xa = np.asarray(m.X.coords); tris = _tri_cells(m.dm) a = _signed_areas(Xa, tris) assert len(Xa) == n0 and len(tris) == nc0 # topology preserved @@ -189,8 +192,9 @@ def test_fault_comb_metric_curved(): arc = np.array([0.2, 0.2]) + 0.42 * np.column_stack( [np.cos(phis), np.sin(phis)]) rho = uw.meshing.fault_comb_metric(m, [arc], cell_size=0.008, n_across=4) - _sm.smooth_mesh_interior(m, metric=rho, method="ma", boundary_slip=False, - method_kwargs=dict(n_outer=1, n_picard=25)) + _sm.smooth_mesh_interior(m, metric=rho, method="mmpde", + boundary_slip=False, + method_kwargs=dict(p=2.0)) Xa = np.asarray(m.X.coords); tris = _tri_cells(m.dm) a = _signed_areas(Xa, tris) assert len(Xa) == n0 and len(tris) == nc0 @@ -325,26 +329,3 @@ def test_compose_metrics_rejects_tensor(): M = uw.meshing.fault_metric_tensor(m, [_SEG3[0]], refinement=3.0, width=0.01) with pytest.raises(ValueError): uw.meshing.compose_metrics([M]) - - -@pytest.mark.tier_a -@pytest.mark.level_1 -@pytest.mark.xfail( - reason="list-of-(metric,weight) composition inside smooth_mesh_interior was " - "an elliptic-ma feature dropped in the development merge (dev's wrapper " - "passes the metric straight to the mover). Compose faults via " - "fault_metric_tensor / fault_comb_metric instead.", - strict=False) -def test_smooth_mesh_interior_list_of_metrics(): - # smooth_mesh_interior accepts a list and composes internally - m = _box(cs=1.0 / 50) - n0 = len(np.asarray(m.X.coords)); nc0 = len(_tri_cells(m.dm)) - r1 = uw.meshing.fault_comb_metric(m, [_SEG3[0]], cell_size=0.008) - r2 = uw.meshing.fault_comb_metric(m, [_SEG3[1]], cell_size=0.008) - _sm.smooth_mesh_interior(m, metric=[(r1, 1.0), (r2, 1.0)], method="ma", - boundary_slip=False, - method_kwargs=dict(n_outer=1, n_picard=25)) - Xa = np.asarray(m.X.coords); tris = _tri_cells(m.dm) - a = _signed_areas(Xa, tris) - assert len(Xa) == n0 and len(tris) == nc0 - assert int((np.sign(a) != np.sign(np.median(a))).sum()) == 0 diff --git a/tests/test_0763_boundary_slip_correctness.py b/tests/test_0763_boundary_slip_correctness.py index f5f3f4c4..cc9ebb9d 100644 --- a/tests/test_0763_boundary_slip_correctness.py +++ b/tests/test_0763_boundary_slip_correctness.py @@ -1,7 +1,7 @@ """Correctness of the mesh-owned ``mesh.boundary_slip`` contract. Step 2 of the boundary tangent-slip refactor swapped every metric mover from the -private ``_ot_adapt._build_slip_projector`` onto ``mesh.boundary_slip`` (see +private slip projector onto ``mesh.boundary_slip`` (see ``docs/developer/design/boundary-slip-strategy.md``) and removed the old projector. This test locks the replacement's behaviour directly: slip vertices land **exactly** on their analytic bounding surface (radius / plane), junctions @@ -82,7 +82,7 @@ def test_box_facet_fallback_stays_on_polygon(): """Unregistered slip labels build transient ``facet`` surfaces; projected vertices lie on the reference-facet polygon and the transient surfaces do not leak into the persistent collection.""" - from underworld3.meshing._ot_adapt import ( + from underworld3.meshing.smoothing import ( _boundary_facets, _nearest_on_facets_2d) m = _box() m.bounding_surfaces.clear() # force the facet fallback path diff --git a/tests/test_0850_mesh_smoothing.py b/tests/test_0850_mesh_smoothing.py index 037227b3..de6e328a 100644 --- a/tests/test_0850_mesh_smoothing.py +++ b/tests/test_0850_mesh_smoothing.py @@ -304,3 +304,111 @@ def test_3x3_mixed_batch_survives_batched_eigh_failure(self): good = np.ones(5, dtype=bool) good[1] = False assert np.array_equal(out[good], M[good]) + + +class TestRetiredMovers: + """2026-07 retirement ruling: the spring / Monge-Ampere / OT-step / + anisotropic-Winslow interior movers were superseded by the MMPDE + mover and deleted. The retired ``method=`` spellings raise a clear + ValueError; ``mesh.OT_adapt`` (built on the OT step) raises a clear + RuntimeError; the default ``method=`` is now ``"mmpde"`` and + ``strategy=`` routes to it (regression for #353, where strategy= + with the old default method='spring' raised TypeError).""" + + @pytest.mark.parametrize("method", ["spring", "ma", "monge-ampere", + "ot", "equidistribute", + "anisotropic", "aniso", "tensor"]) + def test_retired_method_values_raise(self, method): + import sympy + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5) + with pytest.raises(ValueError, match="retired"): + smooth_mesh_interior(mesh, metric=sympy.sympify(1), + method=method) + + def test_unknown_method_raises(self): + import sympy + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5) + with pytest.raises(ValueError, match="unknown method"): + smooth_mesh_interior(mesh, metric=sympy.sympify(1), + method="not-a-mover") + + def test_ot_adapt_tombstone_raises(self): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5) + with pytest.raises(RuntimeError, match="retired"): + mesh.OT_adapt(None) + + def test_retired_private_names_raise(self): + from underworld3.meshing import smoothing as sm + for name in ("_spring_equilibrium_mover", "_monge_ampere_mover", + "_ot_improvement_step", "_winslow_anisotropic", + "_winslow_spring", "_winslow_elliptic", + "_winslow_equidistribute"): + with pytest.raises(ValueError, match="retired"): + getattr(sm, name)() + + def test_default_method_is_mmpde(self): + """smooth_mesh_interior with a metric and NO method= runs the + MMPDE mover (the sanctioned default switch from 'spring').""" + import sympy + mesh = _box_mesh(resolution=6) + x, y = mesh.CoordinateSystem.X + rho = 1 + 4 * sympy.exp(-(((x - 0.5) ** 2 + (y - 0.5) ** 2) + / 0.05)) + before = np.asarray(mesh.X.coords).copy() + smooth_mesh_interior(mesh, metric=rho, + method_kwargs=dict(n_outer=3)) + after = np.asarray(mesh.X.coords) + assert np.all(np.isfinite(after)) + assert not np.allclose(before, after) + + def test_strategy_routes_to_mmpde_without_typeerror(self): + """#353 regression: strategy= injects resolution_ratio, which the + old default mover ('spring') could not bind (TypeError). It now + routes to the mmpde default, which accepts it.""" + import sympy + mesh = _box_mesh(resolution=6) + x, y = mesh.CoordinateSystem.X + rho = 1 + 4 * sympy.exp(-(((x - 0.5) ** 2 + (y - 0.5) ** 2) + / 0.05)) + smooth_mesh_interior(mesh, metric=rho, strategy="med", + skip_threshold=None, + method_kwargs=dict(n_outer=2)) + assert np.all(np.isfinite(np.asarray(mesh.X.coords))) + + +class TestWinslowMMPDEAliasShim: + """Wave C shim contract for the surviving READ-06 alias: the old + ``_winslow_mmpde`` spelling gives the identical result plus exactly + one DeprecationWarning; the new ``_mmpde_mover`` spelling is + warning-free.""" + + def _fresh_mesh_and_metric(self): + import sympy + mesh = _box_mesh(resolution=6) + x, y = mesh.CoordinateSystem.X + rho = 1 + 4 * sympy.exp(-(((x - 0.5) ** 2 + (y - 0.5) ** 2) + / 0.05)) + return mesh, rho + + def test_old_name_identical_result_one_warning(self): + import warnings as _w + from underworld3.meshing.smoothing import (_mmpde_mover, + _winslow_mmpde) + pinned = ("Top", "Bottom", "Left", "Right") + + mesh_new, rho_new = self._fresh_mesh_and_metric() + with _w.catch_warnings(): + _w.simplefilter("error", DeprecationWarning) + _mmpde_mover(mesh_new, rho_new, pinned, False, n_outer=3) + + mesh_old, rho_old = self._fresh_mesh_and_metric() + with pytest.warns(DeprecationWarning, match="_mmpde_mover") as rec: + _winslow_mmpde(mesh_old, rho_old, pinned, False, n_outer=3) + n_dep = sum(1 for w in rec + if issubclass(w.category, DeprecationWarning)) + assert n_dep == 1 + assert np.array_equal(np.asarray(mesh_new.X.coords), + np.asarray(mesh_old.X.coords)) diff --git a/tests/test_0855_slip_surfaces.py b/tests/test_0855_slip_surfaces.py index caf79e96..979d96b5 100644 --- a/tests/test_0855_slip_surfaces.py +++ b/tests/test_0855_slip_surfaces.py @@ -1,7 +1,7 @@ """Named-surface tangent slip for the metric movers. Locks the ``slip_surfaces`` API on ``mesh.boundary_slip`` (the mesh-owned -contract the movers now consume; the private ``_ot_adapt._build_slip_projector`` +contract the movers now consume; the private slip-projector machinery it replaced has been removed — see boundary-slip-strategy.md): * slip-vs-pin is **label-driven** — a boundary vertex slips iff it lies on @@ -18,7 +18,7 @@ import pytest import underworld3 as uw -from underworld3.meshing import _ot_adapt as ota +from underworld3.meshing import smoothing as ota from underworld3.meshing.smoothing import _pinned_mask