From b3439e1f8ae93f219b204637ffd7ef20d3f3e860 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 22:24:16 +1000 Subject: [PATCH 1/7] Add mesh.adapt(engine="edge_split"): longest-edge refinement, no closure Newest-vertex bisection picks the edge to split from a combinatorial tagging rule and then pays a conforming closure to repair the hanging nodes that choice creates. This engine splits the edge the geometry asks for -- the longest edge of every cell still coarser than the metric wants -- and needs no closure at all, because splitting an edge divides *every* incident cell at the same new vertex. There is no hanging node to repair and no longest-edge-propagation chain, so refinement cannot escape the marked region: the refined band hugs the feature instead of a halo around it. No new topology code. `uwnvb_bisect` in nvb_transform.c is already a registered DMPlexTransform driven by a per-edge label, works on triangles and tets, and is the primitive NVB uses for each sub-pass. This engine drives it from Python and therefore inherits star-forest propagation, co-partitioning, labels and coordinates for free. Marking is on the cell DIAMETER, not (d! V)^(1/d). For bisection the two shrink together and either will do; for any engine that reduces volume without shortening the longest edge they diverge badly -- a measured factor of 3.2 on a centroid-refined mesh, where the volume proxy reports the target met while the mesh is nowhere near resolved. The test asserts the diameter. Measured: 2-D 104 -> 412 cells in 7 passes and 3-D 1472 -> 3933 tets, identical at np=1/2/3/4 with no over-shared facets; through mesh.adapt, a 10-level graded MG tail with all 8 exact half-half prolongations captured (every inserted vertex is an exact float edge midpoint) and TI Stokes converging in 10 V-cycles. 3-D reaches the pass cap -- an edge is shared by more tets so fewer are independent per pass -- so the budget scales with dimension and warns rather than silently truncating. Selection is a deterministic function of geometry, not of iteration order. A greedy sweep produced a partition-dependent mesh (412/412/463/925 cells at np=1/2/3/4, every one of them conforming and individually plausible), so an edge wins only if it beats every competing candidate sharing a cell, with a midpoint-coordinate tie-break. The parallel test asserts the serial cell count because that class of defect is invisible in a serial run. Also fixes _cells_on_edge, which applied the 3-D edge -> face -> cell walk in both dimensions. In 2-D an edge *is* a face, so it asked for the support of a cell, got nothing, and reported that the edge touches no cells at all. It is not yet called from the engine, so nothing was broken, but it fails silently and the shape-repair work needs it. Underworld development team with AI support from Claude Code --- .../mesh-reconnection-and-delaunay-adapt.md | 416 ++++++++++++++++++ .../discretisation/discretisation_mesh.py | 94 +++- src/underworld3/utilities/__init__.py | 1 + src/underworld3/utilities/edge_split.py | 290 ++++++++++++ .../ptest_0843_edge_split_parallel.py | 122 +++++ tests/test_0843_edge_split_adapt.py | 188 ++++++++ 6 files changed, 1108 insertions(+), 3 deletions(-) create mode 100644 docs/developer/design/mesh-reconnection-and-delaunay-adapt.md create mode 100644 src/underworld3/utilities/edge_split.py create mode 100644 tests/parallel/ptest_0843_edge_split_parallel.py create mode 100644 tests/test_0843_edge_split_adapt.py diff --git a/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md new file mode 100644 index 000000000..a7b06dbe1 --- /dev/null +++ b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md @@ -0,0 +1,416 @@ +# Reconnection and Delaunay adapt-on-top + +Status: **investigation, 2-D prototype measured (2026-07-29).** Prototype and +raw numbers in `~/+Simulations/mesh_reconnection_study/`. No `src/` change yet. + +## Why look at this + +`mesh.adapt()` refines by **subdivision** — newest-vertex bisection +(`engine="nvb"`) or PETSc longest-edge (`engine="sbr"`). A subdivision engine +decides *where the new point goes* and *how the cells reconnect* with a single +rule and may never re-wire an existing simplex. Three limits follow: the +conforming closure refines cells the metric never asked for (measured halo 45.8 % +on the 3-D fault band); element shape is inherited from the base and can never be +improved; and there is no coarsening and no anisotropy. + +The goal of adding reconnection is **usability, not mesh quality**. MMG/ParMmg +(`mesh.remesh()`) already produces better-shaped elements than any local operator +will, and we are not trying to beat it. What it costs is control: it repartitions +the whole mesh, so every call destroys the decomposition, the point-SF, the MG +hierarchy and the parent/child lineage, at a cost that scales with the whole mesh +rather than the adapted region. A local, rank-respecting operator that carries +some load imbalance is the better trade inside a running model. + +Two constraints that used to forbid non-nested adaptation are already lifted: +custom-P MG transfers are built from **coordinates**, not nesting +(`utilities/custom_mg.py:54,156,515`), and any conforming simplex mesh keeps the +`exact` point-location capability. And `(coords, cells) → DMPlex → Mesh` is +already a production path (`utilities/nvb.py:673`, driven from +`discretisation_mesh.py:7296`). + +## Finding 1 — a flip cannot be a `DMPlexTransform` (confirmed) + +This was the assumption the whole parallel design rested on, so it was checked +first. It holds. + +`DMPlexTransformGetCone_Internal` +(`petsc/src/dm/impls/plex/transform/interface/plextransform.c:1443-1497`) finds +every cone point of every produced point by **descending from the single source +point `p`** (line 1463 `pp = p`; the loop at 1474-1490 walks `pcone[pcp]`), and +identifies the new point as `(parent point, replica)` via +`DMPlexTransformGetTargetPoint` at line 1495. A child's cone can therefore only +reference points in the **source point's own transitive closure**. A 2↔3 flip +produces tets using the apex of the *other* parent tet, outside that closure. +The `celltransform` op signature (`dmplextransformimpl.h:22`) and the +`offset[ct/rt][ctNew]` numbering scheme say the same thing: every new point +belongs to exactly one old point. `DMPlexTransform` is structurally a +**subdivision** framework. + +**Consequence.** The NVB Route-B precedent (`nvb_transform.c`, which inherits SF +propagation and the parallel closure from PETSc) does not carry over. There is no +C-transform fallback. The parallel route for reconnection must be +**freeze the seam** — forbid any cavity containing a cell incident on a shared +plex point, so every shared point is untouched, the point-SF is inherited +verbatim, and each rank rebuilds its local DM alone. No cross-rank closure, no SF +reconciliation, no collective fixpoint. + +## Finding 2 — reconnection repairs shape but not a bad point set + +The study set out to test a specific proposal: centroid (Alfeld) placement is +cheap, local and closure-free, and was ruled a shallow tool only because of +shape (3-D: max dihedral 179.6°, 70.7 % of cells below `q = 0.1`, manufactured +Poisson error stalling at 0.119→0.122→0.124→0.124). The premise was that the +1→3 star split is a bad *connectivity* choice, not a property of centroid +*placement*, and that deciding connectivity afterwards by a Delaunay criterion +would free the placement. + +**The premise is wrong.** Reconnection helps a great deal and is still not +enough. 2-D, same size field, same P1 solver, only the engine differs; in-band +interpolation error (no solve involved): + +| engine | h=0.08 | 0.04 | 0.02 | 0.01 | 0.005 | +|---|---|---|---|---|---| +| nvb (bisection) | 0.2819 | 0.1750 | 0.1029 | 0.0651 | **0.0524** | +| centroid raw | 0.2819 | 0.2285 | 0.1843 | 0.1549 | **0.1469** | +| centroid + flip | 0.2819 | 0.1827 | 0.1358 | 0.1171 | **0.1166** | + +Flips take centroid refinement from 8.06 % to **0.00 %** of cells below `q=0.1`, +and the 99th-percentile max angle from 175.5° to 132.4°. The in-band FE error +goes from *diverging* (0.898 → 1.374 with increasing DOFs) to flat (~0.78). But +the plateau is still **2.2×** bisection's. + +### The structural reason, measured + +A centroid star split leaves the parent's three edges untouched, so an original +edge inside a refined region can never be shortened by further centroid +refinement. Survival of the 178 base edges inside the refinement band: + +| engine | surviving | +|---|---| +| centroid raw | **178/178 (100 %)** | +| centroid + flip | 81/178 (46 %) | +| nvb (bisection) | 42/178 (24 %) | + +A flip *can* remove such an edge — that is exactly why flips help — but only +about half are removable, because a flip needs a convex quad and only exchanges +one diagonal for another. Delaunay optimises connectivity **given the points**; +centroid points are simply the wrong points. + +### A live measurement trap this exposed + +The production marking criterion is `h = sqrt(2A)` (`NVBMesh.centroids_h`). It is +the right proxy for bisection, which shrinks area and diameter together, and a +**misleading** one for any area-reducing split. On the centroid mesh the area +proxy reads `h = 0.0102` while the median in-band **diameter is 0.0331** — a 3.2× +overstatement. The size field is satisfied and the mesh is not resolved. Any +future engine that is not pure bisection must mark on the diameter. + +## Finding 3 — the corrected design, and it is competitive + +Keep the part of the proposal that was right (closure-free placement) and fix the +part that was wrong (place points that shorten diameters): + +> **`edge-split + flip`** — mark on the cell diameter; insert the midpoint of the +> longest edge, splitting that edge in **both** incident cells; then run the +> Delaunay flip pass. + +Splitting the shared edge in both incident cells is conforming by construction — +no hanging node, so **no conforming closure and no LEPP chain**. Unlike bisection +we never demand that the neighbour split its *own* preferred edge. The green +cells this leaves are badly shaped, and repairing them is precisely the flip +pass's job. In-band interpolation error: + +| engine | h=0.04 | 0.02 | 0.01 | 0.005 | +|---|---|---|---|---| +| nvb (bisection) | 0.1750 | 0.1029 | 0.0651 | 0.0524 | +| **edge-split + flip** | **0.1136** | **0.0677** | 0.0626 | 0.0610 | + +and in-band FE error 0.089–0.18 against bisection's 0.23–0.58. Base-edge survival +in the band is 22 %, matching bisection's 24 %. + +Marginal value of the flip pass on this engine, at the same size field: + +| h_near | cells no-flip → flip | in-band error no-flip → flip | +|---|---|---| +| 0.01 | 1935 → 1592 (**−18 %**) | 0.0670 → 0.0626 (**−7 %**) | +| 0.005 | 4083 → 3376 (**−17 %**) | 0.0658 → 0.0610 (**−7 %**) | + +Fewer cells *and* lower error. + +## Correction — the earlier 2-D comparison was confounded + +The tables above compare engines at the same nominal `h_near`. **That is not a +valid comparison** and the first version of this note drew a conclusion from it +that does not survive. + +The engines mark on different quantities: NVB on `h = sqrt(2A)` / `(6V)^(1/3)` +(its own `centroids_h`), edge-split on the **diameter**, which is always the +larger number. The same nominal target therefore asks edge-split for a finer +mesh. In 3-D the effect is severe — 36 569 tets against NVB's 9 780 for the same +`h_near`. Error must be compared **at matched DOF**, not at matched nominal +target. + +Re-run in 2-D with each engine's `h_near` bisected to land on ~1700 cells: + +| engine | h_near | cells | q_med | q<0.1 | ang p99 | in-band interp | in-band FE | +|---|---|---|---|---|---|---|---| +| nvb bisection | 0.0101 | 1718 | 0.862 | 0.00 % | 118.8 | 0.0666 | 0.2674 | +| centroid raw | 0.0089 | 1704 | 0.592 | 11.09 % | 176.2 | 0.1525 | 1.2607 | +| centroid + flip | 0.0089 | 1704 | 0.874 | 0.00 % | 132.9 | 0.1170 | 0.7858 | +| edge-split, no flip | 0.0166 | 1716 | 0.855 | 0.00 % | 136.2 | 0.0725 | 0.1318 | +| **edge-split + flip** | 0.0152 | 1676 | **0.890** | 0.00 % | **115.1** | **0.0620** | **0.0826** | + +The conclusions hold at matched DOF — edge-split + flip is 7 % better than NVB +on interpolation error and **3.2× better** on in-band FE error — but they now +rest on a fair comparison. Note also that in 2-D the flip pass is worth having: +it improves interpolation error 14 % and FE error 37 % at equal cell count. + +## Finding 5 — 3-D: placement is the whole story, reconnection is not + +The 3-D case (unit box, `cellSize=0.4`, `refinement=1`, dipping fault at 60°, +matching `~/+Simulations/nvb_3d_adapt_evaluation/`). Work-precision sweep, +in-band FE error against cell count: + +| engine | observed rate | (P1 ideal in 3-D: `N^-0.67`) | +|---|---|---| +| nvb (bisection) | `N^-0.55` | suboptimal | +| centroid | `N^-0.31` | badly suboptimal — this is the stall | +| edge-split | **`N^-0.67`** | **optimal** | +| edge-split + flip | **`N^-0.69`** | **optimal** | + +At matched DOF (~7 000 cells) edge-split gives 0.193 against NVB's ~0.35 — +**1.8× better**, and it is the only engine achieving the optimal P1 rate. + +**But reconnection buys almost nothing in 3-D.** On the good engine, at +`h_near = 0.07`: + +| | edge-split | + flip | +|---|---|---| +| cells | 26 632 | 26 557 | +| in-band FE error | 0.0862 | 0.0851 (−1.3 %) | +| cells with q<0.1 | 0.96 % | 0.94 % | +| runtime | 3 s | **116 s (39×)** | + +And on the centroid engine, flips move `q<0.1` 41.2 % → 27.8 % but leave the +error unchanged (0.6154 → 0.6192) and the max dihedral identical at 178.7°. + +**Why the 2-D result does not carry over.** 2-D Lawson flips reach the *unique* +Delaunay triangulation — a global optimum for the given points. The 3-D 2↔3 / +3↔2 flip set is a weak local search that cannot escape slivers, and a Delaunay +tetrahedralisation contains them anyway. The kernel test shows this directly: a +Delaunay tet mesh of a random cloud has `q_min = 1.9e-3` with 10 % of cells +below `q = 0.1`, and 125 quality-gated flips left `q_min` **exactly unchanged**. + +This is the outcome flagged as the real risk before the port, and it is +confirmed. The literature agrees: clearing slivers needs flips *plus* smoothing +*plus* insertion/deletion (Klingner & Shewchuk), not flips alone. + +## Finding 4 — the price of preserving the decomposition + +Freezing the seam (the synthetic partition at `x = 0.5` deliberately **crosses** +the refinement band — the worst case): + +- 13.2 % of base cells frozen; +- on `edge-split + flip` the cell-count benefit survives (3395 vs 3376) but + accuracy costs ~10 % (0.0673 vs 0.0610) — freezing eats roughly the whole + *accuracy* benefit of flipping while keeping the *cell-count* benefit; +- on the centroid engine the freeze is far more damaging (in-band FE 1.115 vs + 0.802), because that engine depends on flips for basic shape repair. + +**Design rule that falls out: reconnection must be a polish, never load-bearing.** +An engine that needs flips to be correct will suffer at a partition seam; an +engine that uses flips to be *better* degrades gracefully. This is a stronger +argument for `edge-split + flip` than its error numbers. + +## What survives reconnection downstream + +| consumer | pure flips (vertex set fixed) | point insertion | +|---|---|---| +| `child._adapt_prolongation` (exact ½,½, `nvb.py:249`) | survives — matched by *vertex coordinate identity*, which flips do not touch; must be **re-captured** after the DM rebuild, not re-indexed | invalid for new vertices → geometric builder | +| `child._adapt_parent_cells` (any-degree, #425) | invalid — a flipped cell can straddle two coarse cells | invalid | +| `barycentric` / `rbf` custom-P | fine (coordinate-based) | fine | +| `_location_capability` | stays `exact` | stays `exact` | +| boundary / region labels | preserved iff labelled facets are locked | same | +| co-partitioning invariant | preserved iff the seam layer is frozen | same | + +Secondary payoff worth measuring later: `custom_mg.py:60-89` records +Delaunay-vs-mesh cell agreement at 58.8 % (3-D uniform) and 17.1 % (3-D adapt +child). A Delaunay mesh would push that toward 100 % on a convex domain, making +the geometric P1 builder near-exact and attacking the root cause of the #424 +zero-column failure rather than the symptom. + +## Non-negotiables for any implementation + +- **Exact predicates.** Orientation and in-circle/in-sphere must be exact in + sign. Naive float determinants give inconsistent flip decisions and a + non-conforming mesh. The prototype uses a float filter with a `Fraction` + fallback; production wants Shewchuk's adaptive expansions. +- **Locked facets.** Any facet carrying a boundary label, a region interface or + a registered `Surface` must never be flipped and no cavity may swallow one. + This is the constrained-Delaunay part and it is what protects faults and + material interfaces. +- **Orientation guard.** Reject any modification producing non-positive + area/volume (precedent: the snap guard at `discretisation_mesh.py:7217`). + `to_dm` also requires CCW winding (`nvb.py:687`). +- **Mark on the diameter, not `sqrt(2A)`** — see Finding 2. +- **Never mutate a live mesh's topology in place.** Go `arrays → to_dm → new + Mesh`, as `adapt()` does; the in-place route hits the known `_nav_coords` / + face-control-point staleness traps (issues #286, #135). + +## Finding 6 — end-to-end: TI weak-plane Stokes on an edge-split child + +The engine carries a real solve. Same shear box, fault, and constitutive model as +the validated reference (`~/+Simulations/shear_box_fault_study/shear_box_fault_ti.py`); +only the refinement engine differs. Irregular base, 1056 cells, `max_levels=3`. + +| | cells | q_med | q<0.1 | in-band slip (TI − uniform) | +|---|---|---|---|---| +| nvb child | 2081 | 0.944 | 0.00 % | 0.193 | +| edge-split child | 2868 | 0.959 | 0.00 % | 0.215 | + +Both converge; the weak plane localises slip as it should +(`eta_1` verified to dip to exactly 1e-3 at the fault and recover by `|d| = 0.05`, +director = (0.866, 0.5) = the exact fault normal). + +**Multigrid via the geometric route works, with one wiring step.** Measured +velocity-block preconditioner: + +| child | velocity-block PC | +|---|---| +| nvb (from `base.adapt`) | `mg` — automatic, 8-mesh custom-P tail | +| edge-split, as built | `gamg` — falls back | +| edge-split + `set_custom_fmg(s, base._coarse_level_meshes(), field_id=0)` | **`mg`** | + +The gap is only that `_custom_mg_coarse_meshes` is attached by `_adapt_nested`, +not by `Mesh` construction, so a child built from arrays has no tail until one is +attached. One line when this is wired in as a real engine — the transfers +themselves need no work, since every inserted vertex is an exact edge midpoint. + +**What is NOT shown here.** Iteration counts. `getLinearSolveIterations()` returns +1–2 for every configuration, which cannot be right for a saddle-point solve and +means the counter is not capturing the nested KSP work. No MG-vs-GAMG performance +claim should be read off this run; the run establishes that the path *works*, not +that it is fast. + +**Measurement errors made and fixed along the way** (both would have produced a +confident wrong answer): +- Slip was first sampled at ±3·`w_mech`, *outside* the weak zone, where the two + sample points differ in `y` and the imposed simple shear dominates. The fix is + to sample at `w_mech` and subtract a uniform-viscosity solve on the same mesh, + so what is reported is the fault's contribution alone. +- The first run used `regular=True`, whose right-isoceles cells are a single + similarity class that *both* engines preserve — every mesh scored q = 0.866 and + the comparison could not discriminate. An irregular base is required. + +## Finding 7 — `relax()` is the cheapest win, and it composes with flips + +Everything above was measured **unrelaxed**. `mesh.relax()` (MMPDE in the ideal +reference frame) moves nodes without changing topology or the size distribution, +and its own docstring names the cause this study reached independently: +*"refinement chooses where new nodes go from combinatorics … never from geometry, +so a refined mesh carries needles and slivers that reflect the base mesh's +arbitrary choices"*. Flips answer that from the connectivity side; relax answers +it from the position side. + +Relax-at-end, same base and size field, 2-D: + +| engine / state | cells | q_med | q<0.1 | ang p99 | diam band | in-band interp | +|---|---|---|---|---|---|---| +| nvb | 1732 | 0.862 | 0.00 % | 118.8 | 0.0168 | 0.0651 | +| + relax | 1732 | 0.908 | 0.00 % | 108.1 | 0.0163 | **0.0535** | +| centroid | 1564 | 0.658 | 8.06 % | 175.5 | 0.0331 | 0.1549 | +| + relax | 1564 | 0.696 | 2.37 % | 170.6 | 0.0332 | **0.1791 (worse)** | +| centroid + flip | 1546 | 0.874 | 0.00 % | 132.4 | 0.0137 | 0.1171 | +| + relax | 1546 | 0.929 | 0.00 % | 113.7 | 0.0145 | **0.0869** | +| edge-split | 2580 | 0.835 | 0.00 % | 142.0 | 0.0101 | 0.0670 | +| + relax | 2580 | 0.881 | 0.00 % | 115.8 | 0.0104 | **0.0533** | +| edge-split + flip | 2238 | 0.894 | 0.00 % | 114.8 | 0.0109 | 0.0626 | +| + relax | 2238 | **0.941** | 0.00 % | **100.9** | 0.0108 | **0.0521** | + +Three things follow. + +**It is free and it is the largest single accuracy gain measured here.** Cell +counts are identical and band resolution is preserved to ~1 %, yet in-band +interpolation error drops **17–20 %** on every healthy configuration — more than +the flip pass buys (7–14 %). + +**Flips and relax compose rather than overlap.** On the centroid child, flips take +0.155 → 0.117 and relax then takes it to 0.087; each gains where the other could +not. `edge-split + flip + relax` is the best mesh measured anywhere in this study +(q_med 0.941, ang p99 100.9°, interp 0.0521) and beats relaxed NVB (0.0535). + +**Relax makes centroid refinement WORSE** — the only regression in the table +(0.1549 → 0.1791). Shape improves (q<0.1 8.06 % → 2.37 %) while accuracy +degrades: relax equalises shape and in doing so pulls nodes off where the feature +needs them, because the centroid point *set* was wrong to begin with. That is the +same conclusion as Findings 2 and 5, reached from a third independent direction — +**bad placement cannot be rescued, by connectivity or by position**. + +## Recommended next steps + +The investigation set out to add **reconnection**. What it actually found is a +better **placement** rule, twice over — and that reconnection matters in 2-D and +essentially not at all in 3-D. The recommendation follows that, not the original +premise. + +1. **Do not pursue centroid placement.** Finding 2 closes it: the limitation is + structural, and reconnection recovers less than half of it. +2. **`edge-split` is the candidate engine, and the flip pass is optional.** Mark + on the diameter, split the longest edge in both incident cells. It is + closure-free (the property that motivated centroid refinement), dimension- + general with no pattern tables, and the only engine measured at the optimal P1 + rate in 3-D. Ship the flip pass as an **opt-in polish**: worth it in 2-D + (−14 % interpolation error at equal cells), not worth it in 3-D (−1.3 % for + 39× the runtime). +3. **The 3-D sliver question is open and is not answered by flips.** If element + quality in 3-D needs to improve further, the lever is smoothing between + sweeps (`mesh.relax` exists) or insertion/deletion — not a better flip set. + Worth knowing that `edge-split` already reaches 0.96 % of cells below q=0.1 + against NVB's 1.92 %, so this may not need solving at all. +4. **Then** wire `engine="delaunay"` (better named `engine="edge-split"`) into + `_adapt_nested` alongside `"nvb"`/`"sbr"`. The MG hierarchy needs no new work: + every inserted vertex is an exact edge midpoint, so the recorded ½,½ + prolongation applies unchanged, and the geometric route already handles + everything else. + +The **repair pass** is a separate job with its own handoff plan: +`~/.claude/plans/parallel-mesh-reconnection-flips.md`. Scoped as +*bisection-artefact repair* rather than general mesh improvement: the only badly +shaped cells are those split at an edge they did not nominate, so they all lie in +the star of a newly inserted vertex and are known without search. The plan tiers +the response — strengthen the edge **selection** first (no new operator, and the +existing parallel machinery already covers it), then 2-D Lawson restricted to +new-vertex stars, then 3-D edge removal on the same stars only if a deficit +remains. + +Two corrections it carries, which matter for anyone reading Findings 5 and 4 +above: the 3-D flip verdict is **provisional**, because only 2↔3/3↔2 were tested +— the weakest operators in the family; and the seam-cost measurement froze a +fraction of *all* cells rather than of repair sites, so it is pessimistic for the +wrong reason. + +Deferred, explicitly not blocking: anisotropic (metric) predicate; edge collapse +for coarsening. + +## Open questions / caveats + +- **Depth.** The 2-D figures sit at ~3–3.7 levels (log2 of base/finest diameter) + and the sweeps reach ~4.7. Nothing here tests 6–8 levels, where a real fault + model would sit. `edge-split` has no similarity-class bound — the guarantee + newest-vertex bisection gives up front — so its quality at depth is unmeasured. +- **Halo.** `edge-split` shows a much larger "refined finer than asked" fraction + than NVB (90 % vs 42 % in 3-D). Part of that is the diameter-vs-volume marking + mismatch rather than genuine waste, and the metric was not re-tuned for the + diameter criterion. Worth separating before quoting it as leakage. +- The flip pass is O(cells) per round in pure Python and is the runtime cost in + 3-D. A production version would be incremental. + +## Related + +- `NVB_GRADED_ADAPT.md` — the current engine and why SBR cannot grade. +- `nested-vs-geometric-mg-transfers.md` — the coordinate-vs-topological transfer + trade and issue #424. +- `mesh-shape-relaxation.md:179` — the leakage convention used here (per-cell + `log2(h/h_asked)`, explicitly not a single scalar). +- memory `project_centroid_vs_bisection_refinement` — the 3-D centroid ruling this + study was testing. diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index a032ca6ba..7d99fc688 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -6915,12 +6915,22 @@ def adapt(self, metric_field, max_levels=None, node_budget=None, ``"sbr"`` (default) is the nested adapt-on-top path (the refinement engine is then chosen by ``engine``). ``"mmg"`` is a **deprecated shim** that forwards to :meth:`remesh` (in-place, returns ``self``). - engine : {"nvb", "sbr"}, optional + engine : {"nvb", "sbr", "edge_split"}, optional Advanced selector for the nested refinement engine (ignored when ``adapter="mmg"``). Default ``"nvb"`` — graded newest-vertex bisection; ``"sbr"`` is longest-edge bisection (uniform patch, still the right choice when a uniform-finest MG patch is wanted). See above. + + ``"edge_split"`` splits the **longest edge** of every cell coarser + than the metric asks for, and needs no conforming closure at all + because splitting an edge divides every cell incident on it at the + same new vertex. Refinement therefore stays inside the marked region + instead of a bounded halo around it, at the cost of giving up the + similarity-class bound that makes bisection shape-safe at arbitrary + depth. It marks on the cell **diameter** rather than + ``(dim!·vol)^(1/dim)``; see + :mod:`underworld3.utilities.edge_split`. verbose : bool Returns @@ -6990,8 +7000,9 @@ def adapt(self, metric_field, max_levels=None, node_budget=None, return self if adapter != "sbr": raise ValueError(f"adapter must be 'sbr' or 'mmg', got {adapter!r}") - if engine not in ("sbr", "nvb"): - raise ValueError(f"engine must be 'sbr' or 'nvb', got {engine!r}") + if engine not in ("sbr", "nvb", "edge_split"): + raise ValueError( + f"engine must be 'sbr', 'nvb' or 'edge_split', got {engine!r}") return self._adapt_nested( metric_field, max_levels=max_levels, node_budget=node_budget, @@ -7406,6 +7417,83 @@ def _relax_generation(engine_obj, carry, rcarry): f"-> {fe - fs} cells (rank-local)") if not level_dms: current_dm = base_finest.clone() + elif engine == "edge_split": + # Longest-edge refinement with NO conforming closure: splitting an + # edge divides every cell incident on it at the same new vertex, so + # there is no hanging node to repair and refinement cannot escape the + # marked region. Marking is on the cell DIAMETER, not (dim!·vol)^(1/dim) + # — for bisection the two shrink together, but this engine shortens + # the longest edge directly and the volume proxy would report the + # target met while the mesh is still coarse across the feature. + from underworld3.utilities import edge_split + # Independence caps a pass (no cell may carry two split edges), so a + # generation satisfies only some marked cells and the loop re-marks. + # 3D needs more passes than 2D: an edge is shared by more cells there, + # so fewer edges are independent per pass. + n_pass = 8 * dim * max_levels + current_dm = base_finest + for level in range(n_pass): + centroids, _proxy_h, cs = cell_geometry(current_dm) + if centroids.shape[0]: + M = numpy.clip(eval_metric(centroids), 1e-30, None) + h_target = 1.0 / numpy.sqrt(M) + diameter = edge_split.cell_diameters(current_dm) + sel = numpy.where(diameter > h_target)[0] + if node_budget is not None and sel.size > node_budget: + order = numpy.argsort(M[sel])[::-1] + sel = sel[order[:node_budget]] + else: + sel = numpy.empty(0, dtype=int) # rank owns no cells + + marked = [int(cs + j) for j in sel] + _coarse_for_P = current_dm + current_dm, n_split = edge_split.bisect_longest_edges( + current_dm, marked) + # n_split is global, so this stop is collective without a further + # reduction — a rank with nothing marked still enters the split. + if n_split == 0: + if verbose: + uw.pprint(0, f"[adapt] edge_split pass {level}: " + f"nothing to refine") + break + markers_per_level.append(marked) + # Every inserted vertex is the exact float midpoint of a parent + # edge, so the exact parent/child prolongation applies unchanged. + # Capture it BEFORE the snap and any relaxation move it out of + # reach of coordinate matching (#425). + from underworld3.utilities.nvb import ( + nested_prolongation_from_dms as _nested_from_dms, + nested_cell_parents as _nested_parents) + _vP = _nested_from_dms(_coarse_for_P, current_dm) + _nested_Ps.append(_vP) + _nested_parent_cells.append( + None if _vP is None + else _nested_parents(_coarse_for_P, current_dm, _vP)) + snap_level_boundaries(current_dm) + if _relax_mode == "per-generation": + _mg = Mesh(current_dm.clone(), + simplex=self.dm.isSimplex(), + coordinate_system_type=( + self.CoordinateSystem.coordinate_type), + qdegree=self.qdegree, + boundaries=self.boundaries, verbose=False) + _mg.relax(_relax_metric, **(relax_kwargs or {})) + current_dm.setCoordinatesLocal( + _mg.dm.getCoordinatesLocal()) + level_dms.append(current_dm) + if verbose: + fs, fe = current_dm.getHeightStratum(0) + uw.pprint(0, f"[adapt] edge_split pass {level}: split " + f"{n_split} edge(s) -> {fe - fs} cells " + f"(rank-local)") + else: + # Ran out of passes with cells still coarser than the metric. + # Silence here would look like a satisfied size field. + uw.pprint(0, f"[adapt] edge_split: stopped at the {n_pass}-pass " + f"cap with the metric not yet satisfied; raise " + f"max_levels if the feature needs to be finer.") + if not level_dms: + current_dm = base_finest.clone() elif engine == "nvb": # Serial cell-list engines: the slot-based NVBMesh in 2D (until # the native transform adopts the tagged rule — capstone stage diff --git a/src/underworld3/utilities/__init__.py b/src/underworld3/utilities/__init__.py index 0d7c627ec..b3e5e4ecb 100644 --- a/src/underworld3/utilities/__init__.py +++ b/src/underworld3/utilities/__init__.py @@ -94,3 +94,4 @@ def _append_petsc_path(): from . import boundary_flux from . import custom_mg from .custom_mg import set_custom_fmg +from . import edge_split diff --git a/src/underworld3/utilities/edge_split.py b/src/underworld3/utilities/edge_split.py new file mode 100644 index 000000000..877ddb91d --- /dev/null +++ b/src/underworld3/utilities/edge_split.py @@ -0,0 +1,290 @@ +"""Longest-edge refinement without a conforming closure. + +An alternative refinement engine for :meth:`Mesh.adapt`. Where newest-vertex +bisection chooses which edge to split from a combinatorial tagging rule and then +pays a *conforming closure* to repair the hanging nodes that choice creates, this +engine splits the edge the geometry asks for — the longest edge of every cell +that is still coarser than the metric wants — and needs no closure at all, +because splitting an edge divides **every** cell incident on it at the same new +vertex. There is therefore no hanging node to repair, and no +longest-edge-propagation chain: we never require a neighbour to split its *own* +preferred edge. + +Two consequences that matter for adaptation: + +* refinement does not spread beyond the cells the metric marked, so the refined + region hugs the feature rather than a bounded halo around it; +* the marking criterion is the cell **diameter**, not :math:`(d!\\,V)^{1/d}`. + For bisection the two shrink together and either will do. For any engine that + reduces volume without shortening the longest edge they diverge badly — a + measured factor of 3.2 on a centroid-refined mesh, where the volume proxy + reports the target as met while the mesh is nowhere near resolved. + +The topology, coordinates, labels and parallel star-forest are all handled by the +``uwnvb_bisect`` :c:type:`DMPlexTransform` (see +``docs/developer/design/NVB_GRADED_ADAPT.md``), which is the same primitive the +newest-vertex engine uses for each of its sub-passes. That transform bisects a +set of edges named in a per-edge label and requires them to be **pairwise +independent** — no cell may carry two marked edges in one pass — so a pass here +splits an independent subset and the caller iterates. + +Notes +----- +One pass does not necessarily satisfy every marked cell: independence caps how +many edges can be split at once. Drive it in a loop that re-marks from the +current mesh, as :meth:`Mesh.adapt` does. + +Status +------ +Wired into :meth:`Mesh.adapt` as ``engine="edge_split"``. Validated serial and +parallel in 2-D and 3-D: conforming, refinement confined to the marked region, +and the refined mesh identical at np=1/2/3/4. Tests in +``tests/test_0843_edge_split_adapt.py`` and +``tests/parallel/ptest_0843_edge_split_parallel.py``. + +Not yet done: the reconnection (flip) pass that repairs element shape. It is not +expressible as a ``DMPlexTransform`` — a flip's output cells span two parents' +closures, while a transform's children may only reference their own parent's — +so it needs a separate parallel design and is tracked outside this module. +""" + +import numpy as np +from mpi4py import MPI +from petsc4py import PETSc + +import underworld3 as uw + +_BISECT_LABEL = "uwnvb_bisect_edges" + + +def _register_transform(): + """Import the compiled extension that registers ``uwnvb_bisect`` in PETSc.""" + from underworld3.utilities import _nvb_transform # noqa: F401 (registers on import) + + +def _edge_lengths(dm): + """Length of every edge, indexed by ``edge_point - edge_start``.""" + cdim = dm.getCoordinateDim() + vS, _vE = dm.getDepthStratum(0) + eS, eE = dm.getDepthStratum(1) + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, cdim) + ends = np.array([dm.getCone(e) for e in range(eS, eE)], dtype=np.int64) - vS + d = X[ends[:, 0]] - X[ends[:, 1]] + return np.sqrt(np.einsum("ij,ij->i", d, d)) + + +def _cell_edges(dm): + """Edge points of each cell, as a list indexed by ``cell - cell_start``. + + In 2-D a cell's cone is already its edges; in 3-D the cone holds faces, so + the edges come from the transitive closure filtered to the edge stratum. + """ + eS, eE = dm.getDepthStratum(1) + cS, cE = dm.getHeightStratum(0) + if dm.getDimension() == 2: + return [np.asarray(dm.getCone(c), dtype=np.int64) for c in range(cS, cE)] + out = [] + for c in range(cS, cE): + closure = dm.getTransitiveClosure(c)[0] + out.append(np.array([p for p in closure if eS <= p < eE], dtype=np.int64)) + return out + + +def cell_diameters(dm): + """Longest edge length of every cell, in plex cell order. + + This is the quantity the interpolation error of a linear element depends on, + and the one this engine marks against. + """ + L = _edge_lengths(dm) + eS, _eE = dm.getDepthStratum(1) + return np.array([L[edges - eS].max() for edges in _cell_edges(dm)]) + + +def _sf_logical_or(dm, flag): + """Logical-OR a point-indexed flag array over the point star-forest, in place. + + Every rank holding a copy of a shared point ends up with the same value, so a + shared edge chosen for bisection anywhere is split everywhere — the condition + ``uwnvb_bisect`` needs to keep the child point star-forest conforming. + + For a plex point star-forest the leaf and root spaces are BOTH the local point + chart, so the SAME array is passed as leaf data and root data. This mirrors + ``uwnvb_sf_lor`` in ``nvb_transform.c``, which is the proven form. Gathering + the leaves into a separately-indexed buffer first — the obvious reading of the + PetscSF signature — mis-sizes the reduce and corrupts the heap. + """ + # COLLECTIVE, so every rank must reach it: a rank owning no shared point + # still has to participate or its peers block forever. Only a genuinely + # serial run may skip, and that is a communicator-size test — never a test + # of what this rank happens to own. + if uw.mpi.size == 1: + return flag + sf = dm.getPointSF() + try: + nroots, _ilocal, _iremote = sf.getGraph() + except (ValueError, TypeError): + # An unpopulated star-forest reports a negative root count that petsc4py + # cannot shape an array from; nothing is shared, so nothing to reconcile. + return flag + if nroots < 0: + return flag + + sf.reduceBegin(MPI.INT32_T, flag, flag, MPI.LOR) + sf.reduceEnd(MPI.INT32_T, flag, flag, MPI.LOR) + sf.bcastBegin(MPI.INT32_T, flag, flag, MPI.REPLACE) + sf.bcastEnd(MPI.INT32_T, flag, flag, MPI.REPLACE) + return flag + + +def _owned_count(dm, points): + """How many of ``points`` this rank owns, i.e. holds as a root not a leaf.""" + if uw.mpi.size == 1: + return len(points) + try: + _nroots, ilocal, _iremote = dm.getPointSF().getGraph() + except (ValueError, TypeError): + # Unpopulated star-forest: nothing is shared, so every point is owned. + return len(points) + if ilocal is None or len(ilocal) == 0: + return len(points) + leaves = set(int(p) for p in ilocal) + return sum(1 for p in points if int(p) not in leaves) + + +def _edge_strength(dm): + """Per-edge sort key making "the strongest candidate in a cell" well defined. + + Length decides; the midpoint coordinate breaks ties. Both are computed from + the coordinates alone, so the key is identical on every rank holding the edge + and the selection below is independent of the partition — the property that + makes the refined mesh the same at any communicator size. + """ + cdim = dm.getCoordinateDim() + vS, _vE = dm.getDepthStratum(0) + eS, eE = dm.getDepthStratum(1) + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, cdim) + ends = np.array([dm.getCone(e) for e in range(eS, eE)], dtype=np.int64) - vS + d = X[ends[:, 0]] - X[ends[:, 1]] + length = np.sqrt(np.einsum("ij,ij->i", d, d)) + mid = 0.5 * (X[ends[:, 0]] + X[ends[:, 1]]) + return length, mid + + +def _independent_edges(dm, candidates): + """The candidates that beat every competing candidate sharing a cell. + + This replaces a greedy sweep, which would depend on iteration order and + therefore on the partition (measured: 414 cells at np=1/2 but 463 at np=3 and + 925 at np=4). A candidate is *vetoed* when a stronger candidate shares one of + its cells; vetoes are OR-ed across ranks so a shared edge is judged against + the cells on both sides. What survives is independent by construction — two + edges in the same cell cannot both beat the other — and is a function of the + geometry only. + """ + eS, eE = dm.getDepthStratum(1) + cS, cE = dm.getHeightStratum(0) + pStart, pEnd = dm.getChart() + + is_candidate = np.zeros(pEnd - pStart, dtype=np.int32) + if len(candidates): + is_candidate[np.asarray(candidates, dtype=np.int64) - pStart] = 1 + _sf_logical_or(dm, is_candidate) + + length, mid = _edge_strength(dm) + veto = np.zeros(pEnd - pStart, dtype=np.int32) + edges_of = _cell_edges(dm) + for c in range(cS, cE): + edges = edges_of[c - cS] + rival = edges[is_candidate[edges - pStart] == 1] + if len(rival) < 2: + continue + keys = [(length[e - eS], *mid[e - eS]) for e in rival] + winner = rival[int(np.lexsort(np.array(keys).T[::-1])[-1])] + veto[rival[rival != winner] - pStart] = 1 + _sf_logical_or(dm, veto) + + chosen = np.flatnonzero((is_candidate == 1) & (veto == 0)) + pStart + return chosen[(chosen >= eS) & (chosen < eE)] + + +def _cells_on_edge(dm, edge): + """Cells incident on an edge — the star of the vertex a split would insert. + + The walk up from an edge is dimension-dependent and getting it wrong fails + silently rather than loudly. In 2-D an edge *is* a face, so its support is + already the cells; in 3-D the support holds faces and the cells are one level + further up. Applying the 3-D walk in 2-D asks for the support of a cell, which + is empty, so the function returns no cells at all and every caller reads "this + edge touches nothing". + """ + cS, cE = dm.getHeightStratum(0) + if dm.getDimension() == 2: + return sorted(int(c) for c in dm.getSupport(edge) if cS <= c < cE) + seen = set() + for f in dm.getSupport(edge): + for c in dm.getSupport(f): + if cS <= c < cE: + seen.add(int(c)) + return sorted(seen) + + +def bisect_longest_edges(dm, cells): + """Split the longest edge of as many of ``cells`` as one pass allows. + + Parameters + ---------- + dm : PETSc.DMPlex + Simplex mesh to refine. Not modified. + cells : array of int + Plex cell points to refine. + + Returns + ------- + refined : PETSc.DMPlex + A fresh DM, co-partitioned with ``dm`` and carrying its labels forward. + n_split : int + Number of edges bisected globally. Zero means the pass was empty and the + caller should stop. + + Notes + ----- + Independence caps one pass, so a cell marked here may still exceed the metric + afterwards. Re-mark from the returned mesh and call again. + """ + _register_transform() + + cS, _cE = dm.getHeightStratum(0) + eS, _eE = dm.getDepthStratum(1) + L = _edge_lengths(dm) + edges_of = _cell_edges(dm) + + wanted = {int(edges_of[int(c) - cS][np.argmax(L[edges_of[int(c) - cS] - eS])]) + for c in cells} + chosen = _independent_edges(dm, np.array(sorted(wanted), dtype=np.int64)) + + # Count OWNED edges only: a shared edge is held by every rank on the seam, so + # summing local counts would report it once per sharer and overstate the pass. + n_split = uw.mpi.comm.allreduce(int(_owned_count(dm, chosen)), op=MPI.SUM) + if n_split == 0: + return dm, 0 + + work = dm.clone() + work.createLabel(_BISECT_LABEL) + label = work.getLabel(_BISECT_LABEL) + label.setDefaultValue(0) + for e in chosen: + label.setValue(int(e), 1) + + transform = PETSc.DMPlexTransform().create(comm=work.comm) + transform.setType("uwnvb_bisect") + transform.setDM(work) + transform.setUp() + refined = transform.apply(work) + transform.destroy() + + # The transform copies its driving label onto the output, where it would be + # read as a stale request by the next pass. + if refined.hasLabel(_BISECT_LABEL): + refined.removeLabel(_BISECT_LABEL) + return refined, n_split diff --git a/tests/parallel/ptest_0843_edge_split_parallel.py b/tests/parallel/ptest_0843_edge_split_parallel.py new file mode 100644 index 000000000..9f6be15b3 --- /dev/null +++ b/tests/parallel/ptest_0843_edge_split_parallel.py @@ -0,0 +1,122 @@ +"""Parallel confluence of ``engine="edge_split"``. + +The refined mesh must be the SAME at any communicator size. This is the +load-bearing test for the engine: three separate defects during development +showed up here and nowhere else — + +- a collective (the star-forest reconcile) reached inside a rank-local branch, + which deadlocked as soon as one rank owned no shared point; +- an edge selection by greedy sweep, whose result depends on iteration order and + therefore on the partition (412 cells at np=1/2 but 463 at np=3 and 925 at + np=4, all conforming, all plausible-looking in isolation); +- a mis-sized ``PetscSF`` reduce buffer, which corrupted the heap only after the + second pass. + +None of them is visible in a serial run, and the first two are invisible in a +single-pass run. The reference numbers are asserted in the serial file +(``tests/test_0843_edge_split_adapt.py::test_serial_reference_for_parallel_confluence``) +so a change to the contract is visible there rather than as a mysterious +parallel failure here. + +Run with: + mpirun -n 2 python -m pytest --with-mpi tests/parallel/ptest_0843_edge_split_parallel.py + mpirun -n 3 python -m pytest --with-mpi tests/parallel/ptest_0843_edge_split_parallel.py +""" +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.utilities import edge_split + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, + pytest.mark.tier_b, pytest.mark.timeout(300)] + +# The serial reference: see the serial test file. +SERIAL_BASE_CELLS = 104 +SERIAL_REFINED_CELLS = 412 +SERIAL_PASSES = 7 + +CENTRE = np.array([0.35, 0.6]) + + +def _owned_cells(dm): + cS, cE = dm.getHeightStratum(0) + try: + _nroots, ilocal, _iremote = dm.getPointSF().getGraph() + except (ValueError, TypeError): + ilocal = None + leaves = set() if ilocal is None else {int(p) for p in ilocal} + return uw.mpi.comm.allreduce( + sum(1 for c in range(cS, cE) if c not in leaves)) + + +def _over_shared_facets(dm): + fS, fE = dm.getHeightStratum(1) + return uw.mpi.comm.allreduce( + sum(1 for f in range(fS, fE) if len(dm.getSupport(f)) > 2)) + + +def _centroids(dm): + cS, cE = dm.getHeightStratum(0) + if cE == cS: + return np.zeros((0, dm.getCoordinateDim())) + return np.array([dm.computeCellGeometryFVM(c)[1] for c in range(cS, cE)]) + + +def _h_target(cen): + d = np.linalg.norm(cen - CENTRE, axis=1) + return np.where(d < 0.2, 0.05, 0.3) + + +def test_refined_mesh_is_independent_of_the_partition(): + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.35, + refinement=1, qdegree=2) + dm = base.dm_hierarchy[-1] + assert _owned_cells(dm) == SERIAL_BASE_CELLS + + passes = 0 + while passes < 40: + cS, _cE = dm.getHeightStratum(0) + cen = _centroids(dm) + if cen.shape[0]: + sel = np.flatnonzero( + edge_split.cell_diameters(dm) > _h_target(cen)) + cS + else: + sel = np.empty(0, dtype=int) # this rank owns no cells + dm, n_split = edge_split.bisect_longest_edges(dm, sel) + if n_split == 0: + break + assert _over_shared_facets(dm) == 0, f"pass {passes} broke conformity" + passes += 1 + + assert _owned_cells(dm) == SERIAL_REFINED_CELLS, ( + f"np={uw.mpi.size} produced {_owned_cells(dm)} cells; serial gives " + f"{SERIAL_REFINED_CELLS}. The refined mesh must not depend on the " + f"partition.") + assert passes == SERIAL_PASSES + + +def test_adapt_child_is_confluent_and_carries_the_tail(): + """The full ``mesh.adapt`` path, not just the engine.""" + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2, + regular=False, refinement=2, qdegree=3) + + def metric(cen): + d = np.linalg.norm(np.asarray(cen) - np.array([0.4, 0.55]), axis=1) + return 1.0 / np.where(d < 0.2, 0.03, 0.12) ** 2 + + child = base.adapt(metric, max_levels=2, engine="edge_split") + + assert _over_shared_facets(child.dm) == 0 + tail = child._custom_mg_coarse_meshes + assert tail is not None and len(tail) >= 3 + recorded = child._adapt_prolongation + assert recorded and all(P is not None for P in recorded), ( + "the exact prolongation must survive at np>1: the inserted vertices are " + "exact float edge midpoints on every rank") + # Reported so a partition-dependent regression is legible in the log even if + # the cell count assertion below is later relaxed. + uw.pprint(0, f"[ptest_0843] np={uw.mpi.size}: child " + f"{_owned_cells(child.dm)} cells, tail {len(tail)} levels") diff --git a/tests/test_0843_edge_split_adapt.py b/tests/test_0843_edge_split_adapt.py new file mode 100644 index 000000000..37cc831cd --- /dev/null +++ b/tests/test_0843_edge_split_adapt.py @@ -0,0 +1,188 @@ +"""Longest-edge refinement without a conforming closure (``engine="edge_split"``). + +The engine (:mod:`underworld3.utilities.edge_split`) splits the longest edge of +every cell coarser than the metric asks for. Because splitting an edge divides +*every* cell incident on it at the same new vertex there is no hanging node and +no closure, so — unlike bisection — refinement cannot escape the marked region. +It drives the compiled ``uwnvb_bisect`` :c:type:`DMPlexTransform`, the same +primitive the newest-vertex engine uses for each of its sub-passes, so topology, +coordinates, labels and the parallel star-forest are PETSc's. + +What is asserted, and why each test would have caught a real defect found while +building this: + +- **conformity** — no over-shared facet, at every generation; +- **the diameter is what converges** — the size field is expressed as a + diameter, and the volume proxy ``(dim!·vol)^(1/dim)`` is NOT a substitute: it + reported the target met while the mesh was 3.2x coarser across the feature on + a non-bisection engine. A regression to the proxy passes a naive cell-count + check and fails this one; +- **no halo** — cells far from the feature are untouched. This is the property + the engine exists for, and the one a conforming closure gives up; +- **the exact prolongation survives** — every inserted vertex is the exact float + midpoint of a parent edge, so the recorded 1/2,1/2 transfer applies and the + child carries one MG level per generation; +- **partition independence** — the refined mesh is the same at any communicator + size. Three separate defects during development (a collective inside a + rank-local branch, an order-dependent greedy edge selection, and a mis-sized + star-forest reduce) all showed up here and nowhere else, so this is the + load-bearing test. The np>1 half lives in + ``tests/parallel/ptest_0843_edge_split_parallel.py``; this file records the + serial reference the parallel run must reproduce. +""" +import numpy as np +import pytest +import underworld3 as uw +from underworld3.utilities import edge_split + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _box(dim, cell_size, refinement=1): + lo = tuple([0.0] * dim) + hi = tuple([1.0] * dim) + return uw.meshing.UnstructuredSimplexBox( + minCoords=lo, maxCoords=hi, cellSize=cell_size, + refinement=refinement, qdegree=2) + + +def _centroids(dm): + cS, cE = dm.getHeightStratum(0) + if cE == cS: + return np.zeros((0, dm.getCoordinateDim())) + return np.array([dm.computeCellGeometryFVM(c)[1] for c in range(cS, cE)]) + + +def _over_shared_facets(dm): + fS, fE = dm.getHeightStratum(1) + return sum(1 for f in range(fS, fE) if len(dm.getSupport(f)) > 2) + + +def _refine_to(dm, h_of_centroid, max_passes=40): + """Drive the engine until the diameter target is met everywhere.""" + passes = 0 + while passes < max_passes: + cS, _cE = dm.getHeightStratum(0) + cen = _centroids(dm) + if cen.shape[0] == 0: + break + sel = np.flatnonzero(edge_split.cell_diameters(dm) > h_of_centroid(cen)) + cS + dm, n_split = edge_split.bisect_longest_edges(dm, sel) + if n_split == 0: + break + assert _over_shared_facets(dm) == 0, ( + f"pass {passes} left a facet shared by more than two cells") + passes += 1 + return dm, passes + + +def _disc_target(centre, radius, h_near, h_far): + def h(cen): + d = np.linalg.norm(cen - np.asarray(centre), axis=1) + return np.where(d < radius, h_near, h_far) + return h + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_conforming_and_diameter_target_met(dim): + """The mesh stays conforming and the DIAMETER reaches the target.""" + centre = np.array([0.35, 0.5] if dim == 2 else [0.35, 0.5, 0.6]) + h_near, h_far = (0.06, 0.4) if dim == 2 else (0.15, 0.5) + target = _disc_target(centre, 0.25, h_near, h_far) + + dm = _box(dim, 0.35 if dim == 2 else 0.5).dm_hierarchy[-1] + n0 = dm.getHeightStratum(0)[1] + dm, passes = _refine_to(dm, target) + + assert dm.getHeightStratum(0)[1] > n0, "no refinement happened" + assert _over_shared_facets(dm) == 0 + + cen = _centroids(dm) + inside = np.linalg.norm(cen - centre, axis=1) < 0.25 + diameter = edge_split.cell_diameters(dm) + # The engine converges the DIAMETER. The volume proxy is systematically + # smaller, so a regression to marking on it would leave these cells long. + assert diameter[inside].max() <= h_near * 1.001, ( + f"largest diameter in the target region is {diameter[inside].max():.4f}, " + f"target {h_near}") + assert passes >= 1 + + +def test_refinement_does_not_escape_the_marked_region(): + """No halo: cells far from the feature keep their original size. + + A conforming closure necessarily refines beyond the marked set; this engine + must not. Measured against the coarsest cell size of the unrefined base, so + the test states a property rather than a magic number. + """ + centre = np.array([0.3, 0.3]) + target = _disc_target(centre, 0.15, 0.04, 1.0) + + base = _box(2, 0.3) + dm0 = base.dm_hierarchy[-1] + far0 = _far_field_diameter(dm0, centre, 0.45) + dm, _passes = _refine_to(dm0, target) + far1 = _far_field_diameter(dm, centre, 0.45) + + assert far1 == pytest.approx(far0, rel=1e-12), ( + f"cells beyond r=0.45 changed size ({far0:.5f} -> {far1:.5f}); " + f"refinement escaped the marked region") + + +def _far_field_diameter(dm, centre, radius): + cen = _centroids(dm) + far = np.linalg.norm(cen - np.asarray(centre), axis=1) > radius + return float(edge_split.cell_diameters(dm)[far].max()) if far.any() else 0.0 + + +def test_adapt_returns_child_with_graded_mg_tail(): + """``mesh.adapt(engine="edge_split")`` carries the hierarchy and the exact + prolongation for every generation.""" + base = _box(2, 0.2, refinement=2) + centre = np.array([0.4, 0.55]) + + def metric(cen): + d = np.linalg.norm(np.asarray(cen) - centre, axis=1) + h = np.where(d < 0.2, 0.03, 0.12) + return 1.0 / h**2 + + child = base.adapt(metric, max_levels=2, engine="edge_split") + + assert child.parent is base + n_child = child.dm.getHeightStratum(0)[1] + assert n_child > base.dm_hierarchy[-1].getHeightStratum(0)[1] + + tail = child._custom_mg_coarse_meshes + assert tail is not None and len(tail) >= 3, ( + "the child must carry one MG level per refinement generation on top of " + "the base tail; without it the V-cycle count triples") + + # Every inserted vertex is an exact float edge midpoint, so the recorded + # 1/2,1/2 transfer must be available for EVERY generation — a None here means + # coordinate identity was lost and the geometric builder would be used. + recorded = child._adapt_prolongation + assert recorded and all(P is not None for P in recorded), ( + "exact prolongation missing for at least one generation") + + +def test_unknown_engine_is_refused(): + """The engine name is validated, so a typo cannot silently fall back.""" + base = _box(2, 0.4, refinement=1) + with pytest.raises(ValueError, match="edge_split"): + base.adapt(lambda cen: np.ones(len(cen)), max_levels=1, + engine="edgesplit") + + +def test_serial_reference_for_parallel_confluence(): + """Record the serial result the parallel test must reproduce exactly. + + Kept in the serial file deliberately: the parallel counterpart asserts + equality against these numbers, and a change here is then visible as a + change to the contract rather than as a mysterious parallel failure. + """ + centre = np.array([0.35, 0.6]) + target = _disc_target(centre, 0.2, 0.05, 0.3) + dm = _box(2, 0.35).dm_hierarchy[-1] + assert dm.getHeightStratum(0)[1] == 104 + dm, passes = _refine_to(dm, target) + assert (dm.getHeightStratum(0)[1], passes) == (412, 7) From 711581079c590bafb8a156403e62ee6561b36124 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 22:54:36 +1000 Subject: [PATCH 2/7] Add reconnection repair: mesh.adapt(engine="edge_split", repair=True) Reconnection is the missing third operation of the refine / swap / smooth triple. UW3 had refine (mesh.adapt) and smooth (mesh.relax); this is swap. Refinement chooses where a vertex goes but not how the surrounding cells reconnect, so a cell dragged into a split at an edge it did not nominate gains a thin child. This repairs that by Lawson flips. The acceptance criterion is NOT Delaunay, although these are Lawson flips. Delaunay maximises the minimum angle and says nothing about the maximum, while the P1 interpolation bound depends on the maximum angle and not the minimum (Babuska-Aziz). The two disagree in practice and not marginally: flipping a gmsh-generated mesh towards Delaunay was measured to RAISE the 99th-percentile maximum angle from 126.8 to 129.3 degrees, because gmsh optimises element shape rather than the empty-circle property and its triangulation is locally non-Delaunay exactly where it chose a better-shaped configuration. Since every UW3 mesh starts from gmsh, a repair pass that can degrade one is unusable. Gating on the angle directly makes the pass monotone: it can decline, but it cannot make a mesh worse. Measured on the production path (numbers in the study directory, see the module docstring). As a post-pass it fixes shape and only shape -- decisively on a poor base (99th-percentile maximum angle 156.0 -> 115.1 degrees on an aspect-ratio-4 base; slivers below q=0.1 3.84% -> 0.00% on a non-Delaunay one) and hardly at all on a gmsh base, with interpolation error barely moving either way. Run between refinement passes it also changes where later vertices land, because a flip changes which edge of a cell is longest, and that is worth 20-30% lower error per degree of freedom on a degraded base. The accuracy gain is therefore a placement gain that reconnection unlocks, not a connectivity gain. Parallel by the frozen seam: no cavity may contain a cell incident on a shared plex point. Measured cost 0.9-3.5% of repair sites at 56k cells and np=2..8, halving with every halving of the target size, because repair sites scale with the refined band while the sites a seam crosses stay O(1). The DM is rebuilt on the SAME point chart. A 2-D flip adds and removes no points -- the quad keeps its four vertices, five edges and two cells, and only the diagonal edge's cone and the two cell cones change -- so preserving the numbering lets the point star-forest transfer verbatim, labels transfer by point id and coordinates transfer unchanged. That removes the whole reconstruct-the-star-forest-by-matching-seam-coordinates stage, and with it the class of defect nvb._exact_vertex_map exists to refuse. Surgery on the source DM is not an option: DMPlexSymmetrize refuses to run on a plex that already has supports and nothing outside DMDestroy frees them. The cone orientation convention is derived from the edge cone every time rather than assumed, because getting it wrong does not raise -- it silently yields wrong geometry. repair is OFF by default, for one specific reason: edge_split alone produces a partition-independent mesh, identical at any communicator size, and repair gives that up, because which cavities may be flipped depends on where the partitioner drew the seam. Conformity, orientation, volume, labels and the star-forest stay exact at every rank count. Also note the 99th-percentile maximum angle recovers fully under a frozen seam but the absolute maximum does not -- a few of the worst cells sit on the seam and are exactly the untouchable ones. Orientation and in-circle sign errors produce non-conforming meshes, so the orientation predicate carries Shewchuk's static filter and DECLINES when it cannot resolve a sign. Declining is always safe here because a flip is an optimisation, never a requirement, which is what lets a filter stand in for adaptive-precision arithmetic inside a refinement loop. Repair invalidates the cell-parent map used by the any-degree nested MG transfer (a flipped cell can straddle two coarse cells), so it is set to None and a degree-2 space falls back to the geometric builder. The exact vertex prolongation survives untouched: flips move no vertex, and a P1 section numbers its DOFs from the point numbering, which is preserved. Tests: 6 serial, 4 parallel at np=2/3/4. The maximum-angle assertion is the one that caught the Delaunay criterion; the idempotence check cannot -- an inverted criterion is idempotent too, which is exactly how Delaunay passed while degrading the mesh. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 65 ++- src/underworld3/utilities/__init__.py | 1 + src/underworld3/utilities/reconnect.py | 535 ++++++++++++++++++ .../parallel/ptest_0844_reconnect_parallel.py | 157 +++++ tests/test_0844_reconnect_repair.py | 199 +++++++ 5 files changed, 951 insertions(+), 6 deletions(-) create mode 100644 src/underworld3/utilities/reconnect.py create mode 100644 tests/parallel/ptest_0844_reconnect_parallel.py create mode 100644 tests/test_0844_reconnect_repair.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 7d99fc688..799e820f0 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -6830,7 +6830,7 @@ def relax(self, metric=None, *, verbose=False, **kwargs): def adapt(self, metric_field, max_levels=None, node_budget=None, builder=None, adapter=None, engine=None, verbose=False, - relax=False, relax_kwargs=None): + relax=False, relax_kwargs=None, repair=False): r""" Nested **adapt-on-top**: return a refined **child** mesh. @@ -6931,6 +6931,27 @@ def adapt(self, metric_field, max_levels=None, node_budget=None, depth. It marks on the cell **diameter** rather than ``(dim!·vol)^(1/dim)``; see :mod:`underworld3.utilities.edge_split`. + repair : bool, default False + Run a reconnection (Lawson flip) pass after each ``edge_split`` + generation, repairing the element shapes the split leaves behind. 2-D + and ``engine="edge_split"`` only. Off by default for one specific + reason: ``edge_split`` alone produces a **partition-independent** mesh, + identical at any communicator size, and repair gives that up, because + the flips it may perform depend on where the partitioner drew the seam + (a cavity spanning two ranks cannot be flipped). Conformity, + orientation, volume, labels and the star-forest stay exact at every + rank count. + + Worth turning on when the base is poor — anisotropic, graded, relaxed + or read from a file. Measured there: 41 degrees off the 99th-percentile + maximum angle, slivers below q=0.1 from 3.84 % to 0.00 %, and 20-30 % + lower interpolation error per degree of freedom. On a well-shaped gmsh + base it costs a little time and changes little else. It also + invalidates the cell-parent map used by the any-degree nested MG + transfer (a flipped cell can straddle two coarse cells), so a degree-2 + or higher space falls back to the geometric prolongation builder; the + exact vertex prolongation is unaffected because flips move no vertex. + See :mod:`underworld3.utilities.reconnect`. verbose : bool Returns @@ -7003,16 +7024,32 @@ def adapt(self, metric_field, max_levels=None, node_budget=None, if engine not in ("sbr", "nvb", "edge_split"): raise ValueError( f"engine must be 'sbr', 'nvb' or 'edge_split', got {engine!r}") + if repair: + # Refuse rather than silently ignore: a caller asking for repair has a + # badly shaped mesh, and quietly returning an unrepaired one sends them + # looking for the problem somewhere else. + if engine != "edge_split": + raise ValueError( + f"repair=True needs engine='edge_split', got {engine!r}. The " + f"bisection engines carry a similarity-class bound that keeps " + f"child quality tied to the base, so there is nothing for a " + f"flip pass to repair.") + if self.dim != 2: + raise NotImplementedError( + "repair=True is 2-D only. In 3-D no single flip is enough — " + "the operator set has to become quality-gated edge removal. " + "See docs/developer/design/" + "mesh-reconnection-and-delaunay-adapt.md") return self._adapt_nested( metric_field, max_levels=max_levels, node_budget=node_budget, builder=builder, engine=engine, verbose=verbose, - relax=relax, relax_kwargs=relax_kwargs, + relax=relax, relax_kwargs=relax_kwargs, repair=repair, ) def _adapt_nested(self, metric_field, max_levels=2, node_budget=None, builder="barycentric", engine="nvb", verbose=False, - relax=False, relax_kwargs=None): + relax=False, relax_kwargs=None, repair=False): """Core nested adapt-on-top (SBR or NVB engine). See :meth:`adapt`.""" import math from underworld3.utilities import custom_mg @@ -7466,9 +7503,25 @@ def _relax_generation(engine_obj, carry, rcarry): nested_cell_parents as _nested_parents) _vP = _nested_from_dms(_coarse_for_P, current_dm) _nested_Ps.append(_vP) - _nested_parent_cells.append( - None if _vP is None - else _nested_parents(_coarse_for_P, current_dm, _vP)) + if repair: + # Reconnection repairs the cells the split left thin. It + # rebuilds the DM on the SAME point chart, so the vertex + # prolongation just captured stays valid (a P1 section numbers + # DOFs from the point numbering, which is preserved) — but the + # cell-parent map does not: a flipped cell can straddle two + # coarse cells, so the any-degree transfer has to fall back to + # the geometric builder. + from underworld3.utilities import reconnect + current_dm, n_flips = reconnect.flip_to_reduce_max_angle( + current_dm) + _nested_parent_cells.append(None) + if verbose: + uw.pprint(0, f"[adapt] edge_split pass {level}: repaired " + f"with {n_flips} flip(s)") + else: + _nested_parent_cells.append( + None if _vP is None + else _nested_parents(_coarse_for_P, current_dm, _vP)) snap_level_boundaries(current_dm) if _relax_mode == "per-generation": _mg = Mesh(current_dm.clone(), diff --git a/src/underworld3/utilities/__init__.py b/src/underworld3/utilities/__init__.py index b3e5e4ecb..3681d1f38 100644 --- a/src/underworld3/utilities/__init__.py +++ b/src/underworld3/utilities/__init__.py @@ -95,3 +95,4 @@ def _append_petsc_path(): from . import custom_mg from .custom_mg import set_custom_fmg from . import edge_split +from . import reconnect diff --git a/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py new file mode 100644 index 000000000..af6f9ed15 --- /dev/null +++ b/src/underworld3/utilities/reconnect.py @@ -0,0 +1,535 @@ +"""Reconnection: repair the element shapes a refinement pass leaves behind (2-D). + +Refinement engines choose *where* to put a new vertex; they do not get to choose +how the surrounding cells reconnect. :mod:`underworld3.utilities.edge_split` +splits an edge in every cell incident on it, so a cell that nominated that edge +gains two well-shaped children while a cell dragged along — split at an edge it +did not nominate — gains a thin one. Reconnection is the missing third operation +of the classical refine / swap / smooth triple: UW3 has refine (:meth:`Mesh.adapt`) +and smooth (:meth:`Mesh.relax`), and this is swap. + +What it is worth, and where the benefit actually comes from +---------------------------------------------------------- +Measured on the production path — ``edge_split`` refinement of a real DM, flat-core +size field, error over the refined core, raw numbers in +``~/+Simulations/mesh_reconnection_study/results_production_repair.txt``. + +A flip replaces two cells by two cells and inserts no vertex, so a repair pass run +**after** refinement is cell-count neutral and changes connectivity alone: + +========================== =================== ===================== +base mesh 99th-pct max angle core error, same DOFs +========================== =================== ===================== +gmsh box (the normal case) 124.7 -> 120.5 deg -0.4 % +gmsh box, regular 116.6 -> 116.6 deg 0 % +grid, aspect ratio 4 156.0 -> 115.1 deg +0.9 % +non-Delaunay (scrambled) 175.5 -> 118.0 deg -3.1 % +========================== =================== ===================== + +So as a post-pass this fixes **shape and only shape** — decisively on a poor base +(slivers below q=0.1 go 3.84 % to 0.00 %, and the aspect-ratio-4 row loses 41 +degrees of maximum angle) and hardly at all on a gmsh base. Interpolation error +barely moves either way. + +Run **between** refinement passes it does more, because a flip changes which edge +of a cell is longest and therefore where the *next* pass inserts a vertex. That +buys roughly 20-30 % lower core error per degree of freedom on a degraded base +(and 20-30 % fewer cells for the same size field), and nothing on a gmsh base. The +gain is therefore a *placement* gain that reconnection unlocks, not a connectivity +gain — the same conclusion this study reached about centroid refinement, in the +opposite direction. + +The aspect-ratio-4 row is why the pass exists at all. That base has a maximum +angle of 90 degrees — ideal for P1, since the interpolation bound depends on the +maximum angle (Babuska-Aziz) and not the minimum — and longest-edge refinement +*degrades* it to 156, because repeatedly bisecting the longest edge of a +high-aspect-ratio right triangle manufactures obtuse cells. Refinement creates the +problem; only reconnection removes it. + +Scope +----- +The pass considers every edge it is allowed to touch, not only those around +freshly inserted vertices. That is deliberate — the gains on a degraded base come +precisely from repairing connectivity the refinement did not create — but it does +mean a deliberately hand-built triangulation may be re-connected away from the +refined region, which is one reason the pass is opt-in. + +Parallel: the frozen seam +------------------------- +A flip cannot be a :c:type:`DMPlexTransform` — a child's cone may only reference +its own parent's closure, while a flip's output cells use the *other* parent's +apex — so there is no inherited star-forest propagation and the DM must be +rebuilt. It is rebuilt **on the same point chart**: a 2-D flip adds and removes no +points, since the quad keeps its four vertices, five edges and two cells and only +the diagonal edge's cone and the two cell cones change. Preserving the numbering +means the point star-forest transfers verbatim, labels transfer by point id and +coordinates transfer unchanged, with no coordinate matching anywhere. + +That holds because **no cavity may contain a cell incident on a shared plex +point**. A flip across a partition seam would need one rank's cell to reference an +edge living on another rank, which means enlarging its local chart and rebuilding +the star-forest — a much larger job. Freezing the seam instead costs a measured +0.9-3.5 % of repair sites at 56k cells and np=2..8, and that cost *halves* with +every halving of the target cell size, because repair sites scale with the refined +band while the sites a seam crosses stay O(1). + +.. warning:: + + **The repaired mesh is not partition-independent.** ``edge_split`` alone is + bit-confluent — identical at any communicator size — and repair gives that up + by construction, because which cells are frozen depends on where the + partitioner drew the seam. Conformity, orientation, volume, labels and the + star-forest remain exact at every rank count; it is the *choice* of flips near + a seam that differs. This is the same trade adapt-on-top already makes in + preferring local adaptation to global remeshing, but it is a change of contract + relative to the engine, so repair is opt-in rather than automatic. + +A related cost: the 99th-percentile maximum angle recovers fully under a frozen +seam, but the absolute maximum does not — a few of the worst cells sit on the seam +and are exactly the ones that may not be touched. + +Status +------ +2-D only. In 3-D no single flip suffices: the operator set has to become +quality-gated edge removal, and the empty-sphere property is no help either since +a Delaunay tetrahedralisation still contains slivers — measured directly, a +Delaunay tet mesh of a random cloud has 10 % of its cells below q=0.1. See +``docs/developer/design/mesh-reconnection-and-delaunay-adapt.md``. +""" + +import numpy as np +from mpi4py import MPI +from petsc4py import PETSc + +import underworld3 as uw + +# Labels PETSc maintains itself. They are rebuilt by ``stratify`` so they must not +# be copied onto a fresh plex, and an edge carrying one is not an interface. +_TOPOLOGY_LABELS = ("depth", "celltype") + +# Shewchuk's static filters (Robust Predicates, 1997) with eps = 2^-53. A +# determinant whose magnitude clears the bound has a certain sign; one that does +# not is reported as _UNCERTAIN and the caller declines to act. +# +# Declining is always safe: a flip is an optimisation, never a requirement. That +# is what lets a filter stand in for adaptive-precision arithmetic here. Exact +# arithmetic would resolve more cases and could not resolve any of them wrongly, +# but it is not needed to keep the mesh valid and it is far too slow to sit inside +# a refinement loop. What would NOT be safe is trusting a bare float determinant: +# an inconsistently signed predicate produces a non-conforming mesh. +_EPS = 1.1102230246251565e-16 +_ORIENT_BOUND = (3.0 + 16.0 * _EPS) * _EPS + +_UNCERTAIN = 0 + +# A flip must improve the pair's largest angle by at least this much, measured in +# the cosine. Without a margin, two configurations of equal quality could each +# look like an improvement on the other and the sweeps would cycle. +_MIN_GAIN = 1.0e-9 + + +def _orient2d(pa, pb, pc): + """Sign of the area of triangle ``(pa, pb, pc)``, positive for anticlockwise. + + Returns ``_UNCERTAIN`` when the filter cannot resolve the sign. + """ + acx, acy = pa[0] - pc[0], pa[1] - pc[1] + bcx, bcy = pb[0] - pc[0], pb[1] - pc[1] + left, right = acx * bcy, acy * bcx + det = left - right + if abs(det) >= _ORIENT_BOUND * (abs(left) + abs(right)): + return 1 if det > 0.0 else -1 + return _UNCERTAIN + + +def _smallest_cosine(triangles): + """The most negative cosine of any interior angle across ``triangles``. + + A monotone stand-in for "the largest angle": angle is largest exactly where + its cosine is smallest, and comparing cosines avoids an ``arccos`` per angle + and the tiny non-monotonicity its rounding would introduce near 180 degrees. + """ + worst = 1.0 + for P in triangles: + for i in range(3): + u = P[(i + 1) % 3] - P[i] + v = P[(i + 2) % 3] - P[i] + denom = np.hypot(u[0], u[1]) * np.hypot(v[0], v[1]) + if denom == 0.0: + return -1.0 + worst = min(worst, float((u[0] * v[0] + u[1] * v[1]) / denom)) + return worst + + +# --------------------------------------------------------------- topology reads + +def _coords(dm): + return np.asarray(dm.getCoordinatesLocal().array).reshape( + -1, dm.getCoordinateDim()) + + +def _shared_points(dm): + """Chart-indexed 0/1 flags marking points held by more than one rank. + + Marking the local leaves and OR-ing over the star-forest also flags the roots + on the owning side, so every rank agrees on the seam. Reuses + ``edge_split._sf_logical_or``, whose leaf/root convention — the same array + passed as both leaf and root data — is the one proven correct against + ``uwnvb_sf_lor`` in the C. + """ + from underworld3.utilities import edge_split + + pStart, pEnd = dm.getChart() + flag = np.zeros(pEnd - pStart, dtype=np.int32) + if uw.mpi.size == 1: + return flag + try: + _nroots, ilocal, _iremote = dm.getPointSF().getGraph() + except (ValueError, TypeError): + # An unpopulated star-forest reports a root count petsc4py cannot shape an + # array from. Nothing is shared, so there is no seam. + return flag + if ilocal is not None and len(ilocal): + flag[np.asarray(ilocal, dtype=np.int64) - pStart] = 1 + # COLLECTIVE, and reached on every rank: one that shares nothing still has to + # participate or its peers block. Gate on communicator size, never on what + # this rank happens to own. + edge_split._sf_logical_or(dm, flag) + return flag + + +def _labelled_points(dm): + """Chart-indexed flags for points carrying any non-topological label value. + + A labelled interior edge is an interface — a boundary, a region join, or a + registered ``Surface`` — and must never be flipped, since that is what + protects faults and material boundaries. Over-locking is safe: it declines + repair, it cannot corrupt anything. + """ + pStart, pEnd = dm.getChart() + flag = np.zeros(pEnd - pStart, dtype=bool) + for i in range(dm.getNumLabels()): + if dm.getLabelName(i) in _TOPOLOGY_LABELS: + continue + label = dm.getLabel(dm.getLabelName(i)) + values = label.getValueIS() + if values is None: + continue + for val in values.getIndices(): + points = label.getStratumIS(int(val)) + if points is None: + continue + idx = points.getIndices() + if len(idx): + flag[np.asarray(idx, dtype=np.int64) - pStart] = True + return flag + + +def _cell_regions(dm): + """Per-cell label signature, or ``None`` when every cell carries the same one. + + An edge between two cells with different region values is a material + interface even when the edge itself is unlabelled, so the signature is what + lets those edges be locked. Built from label strata rather than a per-cell + query, which would be one PETSc call per cell per label. + """ + cS, cE = dm.getHeightStratum(0) + names = [dm.getLabelName(i) for i in range(dm.getNumLabels()) + if dm.getLabelName(i) not in _TOPOLOGY_LABELS] + sig = np.zeros((cE - cS, len(names)), dtype=np.int64) + for j, name in enumerate(names): + label = dm.getLabel(name) + values = label.getValueIS() + if values is None: + continue + for val in values.getIndices(): + points = label.getStratumIS(int(val)) + if points is None: + continue + idx = np.asarray(points.getIndices(), dtype=np.int64) + cells = idx[(idx >= cS) & (idx < cE)] + if len(cells): + sig[cells - cS, j] = int(val) + if sig.shape[1] == 0 or np.all(sig == sig[0]): + return None + return sig + + +def _cell_vertices_and_seam(dm, X, shared): + """One closure pass: anticlockwise vertices of every cell, and the seam mask. + + Both need the transitive closure of every cell, so they are computed together + rather than in two passes. + """ + cS, cE = dm.getHeightStratum(0) + vS, vE = dm.getDepthStratum(0) + pStart, _pEnd = dm.getChart() + any_shared = bool(shared.any()) + + verts = np.empty((cE - cS, 3), dtype=np.int64) + frozen = np.zeros(cE - cS, dtype=bool) + for c in range(cS, cE): + closure = np.asarray(dm.getTransitiveClosure(c)[0], dtype=np.int64) + v = [int(p) for p in closure if vS <= p < vE] + if _orient2d(X[v[0] - vS], X[v[1] - vS], X[v[2] - vS]) < 0: + v = [v[0], v[2], v[1]] + verts[c - cS] = v + if any_shared: + frozen[c - cS] = bool(shared[closure - pStart].any()) + return verts, frozen + + +# ------------------------------------------------------------------- the rebuild + +def rebuild_with_cones(dm, new_cells, new_edges): + """Build a fresh plex on the **same point chart** with the given cones replaced. + + Parameters + ---------- + dm : PETSc.DMPlex + Source mesh. Not modified. + new_cells : dict + ``{cell point: (v0, v1, v2)}`` with the vertices anticlockwise. + new_edges : dict + ``{edge point: (va, vb)}``. + + Returns + ------- + PETSc.DMPlex + A new mesh whose chart, coordinates, labels and point star-forest match + the source, differing only in the replaced cones. + + Notes + ----- + Surgery on the source is not possible: ``DMPlexSymmetrize`` refuses to run on + a plex that already has supports, and nothing outside ``DMDestroy`` frees + them. Hence a fresh plex with every untouched cone copied across. + + The cone-orientation convention is derived, not assumed. For a triangle the + closure vertex order is anticlockwise, cone entry ``i`` is the edge joining + closure vertices ``i`` and ``i+1`` (mod 3), and its orientation is ``0`` when + the edge's own cone runs that way and ``-1`` when reversed. A wrong + orientation does not raise — it silently yields wrong geometry — so it is + computed from the edge cone every time. + """ + pStart, pEnd = dm.getChart() + vS, vE = dm.getDepthStratum(0) + eS, eE = dm.getDepthStratum(1) + cdim = dm.getCoordinateDim() + + new = PETSc.DMPlex().create(comm=dm.comm) + new.setDimension(dm.getDimension()) + new.setChart(pStart, pEnd) + for p in range(pStart, pEnd): + new.setConeSize(p, dm.getConeSize(p)) + new.setUp() + + # Edges first: the cell wiring below reads edge cones back to derive + # orientations, so they have to be the new ones already. + for p in range(pStart, pEnd): + if p in new_cells: + continue + if p in new_edges: + new.setCone(p, [int(v) for v in new_edges[p]]) + continue + new.setCone(p, [int(x) for x in dm.getCone(p)]) + orientation = [int(o) for o in dm.getConeOrientation(p)] + if orientation: + new.setConeOrientation(p, orientation) + + edge_of = {} + for e in range(eS, eE): + a, b = (int(v) for v in new.getCone(e)) + edge_of[(a, b) if a < b else (b, a)] = e + + for c, (v0, v1, v2) in new_cells.items(): + cone, orientation = [], [] + for x, y in ((v0, v1), (v1, v2), (v2, v0)): + e = edge_of[(x, y) if x < y else (y, x)] + cone.append(e) + orientation.append(0 if int(new.getCone(e)[0]) == x else -1) + new.setCone(c, cone) + new.setConeOrientation(c, orientation) + + new.symmetrize() + new.stratify() + + # Coordinates verbatim: the vertex points are unchanged, so this is the same + # section over the same chart holding the same values. + new.setCoordinateDim(cdim) + section = new.getCoordinateSection() + section.setNumFields(1) + section.setFieldComponents(0, cdim) + section.setChart(vS, vE) + for v in range(vS, vE): + section.setDof(v, cdim) + section.setFieldDof(v, 0, cdim) + section.setUp() + coords = PETSc.Vec().createSeq(section.getStorageSize(), + comm=PETSc.COMM_SELF) + coords.array[:] = np.asarray(dm.getCoordinatesLocal().array) + new.setCoordinatesLocal(coords) + + # Labels by point id. No coordinate matching is involved, which is the whole + # reason for preserving the numbering. + for i in range(dm.getNumLabels()): + name = dm.getLabelName(i) + if name in _TOPOLOGY_LABELS: + continue + new.createLabel(name) + source, target = dm.getLabel(name), new.getLabel(name) + values = source.getValueIS() + if values is None: + continue + for val in values.getIndices(): + points = source.getStratumIS(int(val)) + if points is None: + continue + for p in points.getIndices(): + target.setValue(int(p), int(val)) + + # The star-forest transfers verbatim: every rank preserves its numbering, so + # the remote point numbers it carries are still the right ones. + if uw.mpi.size > 1: + new.setPointSF(dm.getPointSF()) + return new + + +# ---------------------------------------------------------------- the flip pass + +def _flippable(dm, X, verts, frozen, locked, regions): + """Edges worth flipping, as ``(edge, cell_t, cell_u, p, a, q, b, gain)``. + + ``(p, a, q, b)`` is the quad anticlockwise with ``(a, b)`` the current + diagonal, so the flip replaces cells ``(p, a, b)`` and ``(a, q, b)`` by + ``(p, a, q)`` and ``(p, q, b)``. An edge qualifies on two counts: + + * the quad is **strictly convex**, so the flip cannot invert a cell. Declined + whenever the filtered orientation predicate cannot resolve a sign; + * the flip **strictly reduces the largest of the pair's six angles**. + + The second test is deliberately not the Delaunay (in-circle) criterion, even + though this is a Lawson flip. Delaunay maximises the *minimum* angle and says + nothing about the maximum, while the P1 interpolation bound depends on the + maximum angle and not the minimum (Babuska-Aziz). The two disagree in practice + and not marginally: flipping a gmsh-generated mesh towards Delaunay was + measured to *raise* the 99th-percentile maximum angle from 126.8 to 129.3 + degrees, because gmsh optimises element shape rather than the empty-circle + property and its triangulation is therefore locally non-Delaunay exactly where + it has chosen a better-shaped configuration. Since every UW3 mesh starts from + gmsh, a repair pass that can degrade such a mesh is unusable. Gating on the + angle instead makes the pass monotone by construction: it can decline, but it + cannot make a mesh worse. + """ + eS, eE = dm.getDepthStratum(1) + cS, _cE = dm.getHeightStratum(0) + vS, _vE = dm.getDepthStratum(0) + pStart, _pEnd = dm.getChart() + + out = [] + for e in range(eS, eE): + if locked[e - pStart]: + continue + support = [int(c) for c in dm.getSupport(e)] + if len(support) != 2: + continue # boundary edge: nothing to flip into + t, u = support + if frozen[t - cS] or frozen[u - cS]: + continue + if regions is not None and not np.array_equal(regions[t - cS], + regions[u - cS]): + continue # region interface + a, b = (int(v) for v in dm.getCone(e)) + p = next(v for v in verts[t - cS] if v not in (a, b)) + q = next(v for v in verts[u - cS] if v not in (a, b)) + + # Order the diagonal so the quad p-a-q-b runs anticlockwise. + if _orient2d(X[p - vS], X[a - vS], X[q - vS]) < 0: + a, b = b, a + if _orient2d(X[p - vS], X[a - vS], X[q - vS]) <= 0: + continue # not strictly convex, or unresolved + if _orient2d(X[p - vS], X[q - vS], X[b - vS]) <= 0: + continue + + Xp, Xa, Xq, Xb = X[p - vS], X[a - vS], X[q - vS], X[b - vS] + before = _smallest_cosine(((Xp, Xa, Xb), (Xa, Xq, Xb))) + after = _smallest_cosine(((Xp, Xa, Xq), (Xp, Xq, Xb))) + if after <= before + _MIN_GAIN: + continue # no shape gain worth the flip + out.append((e, t, u, int(p), a, int(q), b, after - before)) + + # Best gain first, so that when two candidate flips share a cell and only one + # can run this sweep, the sweep keeps the better of the two rather than + # whichever the edge loop happened to reach first. + out.sort(key=lambda row: -row[7]) + return out + + +def flip_to_reduce_max_angle(dm, max_sweeps=12): + """Flip edges to reduce the largest element angles, leaving the seam alone. + + Parameters + ---------- + dm : PETSc.DMPlex + A 2-D simplex mesh. Not modified. + max_sweeps : int + Cap on sweeps. Reaching it warns rather than failing silently. + + Returns + ------- + repaired : PETSc.DMPlex + A new mesh on the same point chart, or ``dm`` itself if nothing flipped. + n_flips : int + Flips performed across all ranks. + + Notes + ----- + Each sweep applies an **independent** set of flips — no two sharing a cell — + and rebuilds once, instead of rebuilding per flip. Two flips sharing a cell + would each rewire it from a stale reading of the other's result. Deferring the + loser to the next sweep costs a sweep; not deferring it costs correctness. + + Every accepted flip strictly reduces its own pair's largest angle, so the pass + cannot degrade a mesh. It is not guaranteed to reach a global optimum: reducing + one pair's largest angle can raise a neighbouring pair's, so the sequence is a + local improvement and ``max_sweeps`` is the guard against a pathological case + cycling between configurations. In practice it converges in a few sweeps. + """ + if dm.getDimension() != 2: + raise NotImplementedError( + "reconnect.flip_to_reduce_max_angle is 2-D only. In 3-D no single " + "flip is enough — the operator set has to change to quality-gated " + "edge removal, and a Delaunay tetrahedralisation still contains " + "slivers so the empty-sphere test is no help either. See " + "docs/developer/design/mesh-reconnection-and-delaunay-adapt.md") + + total = 0 + for _sweep in range(max_sweeps): + X = _coords(dm) + verts, frozen = _cell_vertices_and_seam(dm, X, _shared_points(dm)) + candidates = _flippable(dm, X, verts, frozen, _labelled_points(dm), + _cell_regions(dm)) + + claimed = set() + new_cells, new_edges = {}, {} + for e, t, u, p, a, q, b, _gain in candidates: + if t in claimed or u in claimed: + continue + claimed.update((t, u)) + new_edges[e] = (p, q) + new_cells[t] = (p, a, q) + new_cells[u] = (p, q, b) + + # COLLECTIVE, and reached on every rank: one with nothing to flip still + # has to vote or its peers block waiting for it. + n = uw.mpi.comm.allreduce(len(new_edges), op=MPI.SUM) + if n == 0: + break + dm = rebuild_with_cones(dm, new_cells, new_edges) + total += n + else: + uw.pprint(0, f"[reconnect] reached the {max_sweeps}-sweep cap with flips " + f"still pending. The mesh is valid and conforming but not " + f"fully repaired; raise max_sweeps if this matters.") + + return dm, total diff --git a/tests/parallel/ptest_0844_reconnect_parallel.py b/tests/parallel/ptest_0844_reconnect_parallel.py new file mode 100644 index 000000000..b27780b3d --- /dev/null +++ b/tests/parallel/ptest_0844_reconnect_parallel.py @@ -0,0 +1,157 @@ +"""Reconnection repair under a frozen partition seam. + +The parallel contract here is deliberately *not* partition independence. Repair +gives that up by construction: which cavities may be flipped depends on where the +partitioner drew the seam, so the flip set — and therefore the mesh — differs with +rank count. What must hold at every rank count is everything else, and that is +what this file asserts: + +* **no shared point's cone changes.** This is the freeze rule stated as a + postcondition, and it is the load-bearing one. It is what allows the rebuilt DM + to reuse the point star-forest verbatim instead of reconstructing it by matching + seam coordinates — a spatial query standing in for an identity lookup, which is + the failure mode ``nvb._exact_vertex_map`` exists to refuse; +* the **chart, cell count and total area** are invariant, globally; +* the mesh still carries a solve, which is the only real proof the labels and the + star-forest came through usable rather than merely present. + +Run with: + mpirun -n 2 python -m pytest --with-mpi tests/parallel/ptest_0844_reconnect_parallel.py + mpirun -n 3 python -m pytest --with-mpi tests/parallel/ptest_0844_reconnect_parallel.py +""" +import numpy as np +import pytest +from mpi4py import MPI + +import underworld3 as uw +from underworld3.utilities import edge_split, reconnect + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, + pytest.mark.tier_b, pytest.mark.timeout(600)] + +CENTRE = np.array([0.4, 0.55]) + + +def _refined_dm(): + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2, + regular=False, qdegree=2) + dm = base.dm + for _ in range(20): + cS, cE = dm.getHeightStratum(0) + if cE > cS: + cen = np.array([dm.computeCellGeometryFVM(c)[1] + for c in range(cS, cE)]) + d = np.linalg.norm(cen - CENTRE, axis=1) + target = np.where(d < 0.25, 0.05, 0.4) + sel = np.flatnonzero(edge_split.cell_diameters(dm) > target) + cS + else: + sel = np.empty(0, dtype=np.int64) + dm, n = edge_split.bisect_longest_edges(dm, sel) + if n == 0: + break + return dm + + +def _owned(dm, points): + try: + _nroots, ilocal, _iremote = dm.getPointSF().getGraph() + except (ValueError, TypeError): + ilocal = None + leaves = set() if ilocal is None else {int(p) for p in ilocal} + return [p for p in points if p not in leaves] + + +def _global(x, op=MPI.SUM): + return uw.mpi.comm.allreduce(x, op=op) + + +def _owned_cells_and_area(dm): + cS, cE = dm.getHeightStratum(0) + owned = _owned(dm, range(cS, cE)) + area = sum(abs(dm.computeCellGeometryFVM(c)[0]) for c in owned) + return len(owned), area + + +def test_shared_points_are_untouched(): + """The freeze rule as a postcondition — the invariant the design rests on.""" + dm = _refined_dm() + shared = reconnect._shared_points(dm) + pStart, _pEnd = dm.getChart() + idx = np.flatnonzero(shared) + pStart + assert _global(len(idx)) > 0, ( + "no point is shared, so this run cannot exercise the freeze rule") + before = {int(p): tuple(int(x) for x in dm.getCone(p)) for p in idx} + + out, nflips = reconnect.flip_to_reduce_max_angle(dm) + + assert _global(nflips, op=MPI.MAX) > 0, "nothing flipped anywhere" + for p, cone in before.items(): + assert tuple(int(x) for x in out.getCone(p)) == cone, ( + f"rank {uw.mpi.rank}: shared point {p} was rewired; the point " + f"star-forest can no longer be reused verbatim") + + +def test_geometry_and_conformity_survive(): + dm = _refined_dm() + chart = dm.getChart() + ncells, area = _owned_cells_and_area(dm) + + out, _ = reconnect.flip_to_reduce_max_angle(dm) + + assert out.getChart() == chart + ncells_after, area_after = _owned_cells_and_area(out) + assert _global(ncells_after) == _global(ncells) + assert _global(area_after) == pytest.approx(_global(area), rel=1e-12) + + fS, fE = out.getHeightStratum(1) + assert _global(sum(1 for f in range(fS, fE) + if len(out.getSupport(f)) > 2)) == 0 + + +def test_repaired_mesh_still_solves(): + """Labels and the star-forest present is not the same as usable.""" + dm = _refined_dm() + out, _ = reconnect.flip_to_reduce_max_angle(dm) + + mesh = uw.discretisation.Mesh(out, qdegree=2) + u = uw.discretisation.MeshVariable("u_par", mesh, 1, degree=1) + poisson = uw.systems.Poisson(mesh, u_Field=u) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 1.0 + poisson.add_dirichlet_bc(0.0, "All_Boundaries") + poisson.solve() + assert poisson.snes.getConvergedReason() > 0 + + # Integrating 1 exercises the assembled section over the rebuilt topology on + # every rank at once, which a rank-local area sum does not. + one = uw.discretisation.MeshVariable("one_par", mesh, 1, degree=1) + one.array[:, 0, 0] = 1.0 + assert uw.maths.Integral(mesh, one.sym[0]).evaluate() == pytest.approx( + 1.0, rel=1e-10) + + +def test_adapt_with_repair_runs_in_parallel(): + """The full ``mesh.adapt(..., repair=True)`` path.""" + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2, + regular=False, refinement=1, qdegree=3) + + def metric(centroids): + d = np.linalg.norm(np.asarray(centroids) - CENTRE, axis=1) + return 1.0 / np.where(d < 0.2, 0.04, 0.15) ** 2 + + child = base.adapt(metric, max_levels=2, engine="edge_split", repair=True) + + fS, fE = child.dm.getHeightStratum(1) + assert _global(sum(1 for f in range(fS, fE) + if len(child.dm.getSupport(f)) > 2)) == 0 + # Flips move no vertex, so the exact vertex prolongation must survive; the + # cell-parent map must NOT, because a flipped cell can straddle two coarse + # cells and using it would transfer from the wrong parent. + assert child._adapt_prolongation and all( + P is not None for P in child._adapt_prolongation) + assert all(pc is None for pc in child._adapt_parent_cells) + uw.pprint(0, f"[ptest_0844] np={uw.mpi.size}: repaired child " + f"{_global(_owned_cells_and_area(child.dm)[0])} cells") diff --git a/tests/test_0844_reconnect_repair.py b/tests/test_0844_reconnect_repair.py new file mode 100644 index 000000000..2a1e7c5e7 --- /dev/null +++ b/tests/test_0844_reconnect_repair.py @@ -0,0 +1,199 @@ +"""Reconnection (Lawson flip) repair of a refined 2-D mesh. + +The load-bearing checks here are the ones that would pass for the wrong reason if +they were written loosely: + +* the maximum angle must **improve**, and this is the check that earned its place. + The pass originally accepted a flip on the Delaunay criterion, and this + assertion is what caught that Delaunay maximises the *minimum* angle while P1 + interpolation depends on the *maximum* — flipping a gmsh mesh towards Delaunay + raised the 99th-percentile maximum angle. Assert the quantity the method claims + to improve, not a proxy for it; +* **volume conservation**, not orientation. Checking that the new cells are + positively oriented is worthless when they were built anticlockwise by + construction: the check can never fail. Equal total area is the real test; +* the **point chart is unchanged**, which is the invariant that lets the parallel + path reuse the star-forest verbatim rather than reconstructing it. + +Note what the idempotence check below does *not* prove. A second pass flipping +nothing shows the acceptance test is self-consistent, but an **inverted** criterion +is equally idempotent — it would flip every good edge once and then find nothing +more to do. That is exactly how the Delaunay criterion passed here while degrading +the mesh. Idempotence catches oscillation, not a wrong objective. +""" +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.utilities import edge_split, reconnect + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _max_angles(dm): + """Largest interior angle of every cell, in degrees. + + The maximum angle is the quantity a P1 interpolation bound depends on + (Babuska-Aziz); the minimum angle is not, so it is the wrong thing to assert. + """ + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + out = [] + for c in range(cS, cE): + v = [int(p) for p in dm.getTransitiveClosure(c)[0] if vS <= p < vE] + P = X[np.array(v) - vS] + angles = [] + for i in range(3): + u1 = P[(i + 1) % 3] - P[i] + u2 = P[(i + 2) % 3] - P[i] + cos = np.dot(u1, u2) / (np.linalg.norm(u1) * np.linalg.norm(u2)) + angles.append(np.degrees(np.arccos(np.clip(cos, -1.0, 1.0)))) + out.append(max(angles)) + return np.array(out) + + +def _signed_areas(dm): + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + out = [] + for c in range(cS, cE): + v = [int(p) for p in dm.getTransitiveClosure(c)[0] if vS <= p < vE] + a, b, d = X[np.array(v) - vS] + out.append(0.5 * ((b[0] - a[0]) * (d[1] - a[1]) + - (d[0] - a[0]) * (b[1] - a[1]))) + return np.array(out) + + +def _over_shared_facets(dm): + fS, fE = dm.getHeightStratum(1) + return sum(1 for f in range(fS, fE) if len(dm.getSupport(f)) > 2) + + +def _refined_dm(cell_size=0.3, h_near=0.05, centre=(0.35, 0.6), radius=0.2): + """A box mesh refined by ``edge_split`` — the mesh repair is meant to fix.""" + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cell_size, + regular=False, qdegree=2) + dm = base.dm + for _ in range(20): + cS, cE = dm.getHeightStratum(0) + cen = np.array([dm.computeCellGeometryFVM(c)[1] for c in range(cS, cE)]) + d = np.linalg.norm(cen - np.array(centre), axis=1) + target = np.where(d < radius, h_near, 0.4) + sel = np.flatnonzero(edge_split.cell_diameters(dm) > target) + cS + dm, n = edge_split.bisect_longest_edges(dm, sel) + if n == 0: + break + return dm + + +def test_repair_conserves_area_and_conformity(): + dm = _refined_dm() + ncells = dm.getHeightStratum(0)[1] - dm.getHeightStratum(0)[0] + chart = dm.getChart() + area = _signed_areas(dm).sum() + + out, nflips = reconnect.flip_to_reduce_max_angle(dm) + + assert nflips > 0, "nothing to repair — the fixture is not exercising the pass" + # A flip replaces two cells by two cells and adds no points, so both the cell + # count and the whole chart are invariant. The chart being invariant is what + # makes the parallel star-forest reusable. + assert out.getChart() == chart + assert out.getHeightStratum(0)[1] - out.getHeightStratum(0)[0] == ncells + assert _over_shared_facets(out) == 0 + new_areas = _signed_areas(out) + assert (new_areas > 0).all(), "repair inverted a cell" + assert new_areas.sum() == pytest.approx(area, rel=1e-13) + + +def test_repair_improves_the_maximum_angle(): + dm = _refined_dm() + before = _max_angles(dm) + out, _ = reconnect.flip_to_reduce_max_angle(dm) + after = _max_angles(out) + + assert np.percentile(after, 99) < np.percentile(before, 99) + assert after.max() <= before.max() + + +def test_second_pass_flips_nothing(): + """Idempotence: the control that catches an inconsistently signed predicate. + + A kernel that keeps finding improvements in an already-repaired mesh is + reporting a predicate bug, and that bug would be invisible in every other + check here. + """ + dm = _refined_dm() + once, n1 = reconnect.flip_to_reduce_max_angle(dm) + assert n1 > 0 + twice, n2 = reconnect.flip_to_reduce_max_angle(once) + assert n2 == 0 + assert twice is once, "a no-op pass must return the mesh it was given" + + +def test_labels_survive_and_remain_usable(): + dm = _refined_dm() + names = sorted(dm.getLabelName(i) for i in range(dm.getNumLabels())) + sizes = {} + for name in names: + if name in ("depth", "celltype"): + continue + sizes[name] = dm.getLabel(name).getStratumSize( + int(dm.getLabel(name).getValueIS().getIndices()[0])) + + out, _ = reconnect.flip_to_reduce_max_angle(dm) + + assert sorted(out.getLabelName(i) for i in range(out.getNumLabels())) == names + for name, size in sizes.items(): + label = out.getLabel(name) + assert label.getStratumSize( + int(label.getValueIS().getIndices()[0])) == size + + # Labels surviving as point sets is not the same as being usable: the real + # test is that a Dirichlet condition can still be imposed on one. + mesh = uw.discretisation.Mesh(out, qdegree=2) + u = uw.discretisation.MeshVariable("u_rec", mesh, 1, degree=1) + poisson = uw.systems.Poisson(mesh, u_Field=u) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 1.0 + poisson.add_dirichlet_bc(0.0, "All_Boundaries") + poisson.solve() + assert poisson.snes.getConvergedReason() > 0 + assert u.array[:, 0, 0].max() > 0.0 + + +def test_labelled_interior_edges_are_never_flipped(): + """A labelled interior edge is an interface and must survive untouched. + + This is what protects a fault or a material boundary from being reconnected + across. + """ + dm = _refined_dm() + eS, eE = dm.getDepthStratum(1) + interior = [e for e in range(eS, eE) if len(dm.getSupport(e)) == 2] + # Lock a slice of interior edges, including ones the pass would otherwise flip. + dm.createLabel("test_interface") + label = dm.getLabel("test_interface") + locked = interior[::7] + for e in locked: + label.setValue(e, 1) + cones = {e: tuple(int(v) for v in dm.getCone(e)) for e in locked} + + out, _ = reconnect.flip_to_reduce_max_angle(dm) + + for e, cone in cones.items(): + assert tuple(int(v) for v in out.getCone(e)) == cone, ( + f"locked interface edge {e} was flipped") + + +def test_three_dimensions_is_refused(): + """3-D must fail loudly: Delaunay is the wrong criterion, not merely untested.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), cellSize=0.5, + regular=False, qdegree=2) + with pytest.raises(NotImplementedError, match="2-D only"): + reconnect.flip_to_reduce_max_angle(mesh.dm) From 79e0a5d047aec6a619097449f313ad707f3f292c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 22:57:10 +1000 Subject: [PATCH 3/7] docs: record the repair pass and three corrections to the reconnection findings Finding 8 in the reconnection design note. Three of the earlier findings needed correcting rather than extending: - Delaunay is the wrong acceptance criterion in 2-D as well as 3-D. Finding 3 treated it as settled because Lawson flips reach the unique Delaunay triangulation; that settles the operator, not the criterion. Delaunay maximises the minimum angle while P1 interpolation depends on the maximum, and flipping a gmsh mesh towards Delaunay was measured to raise the 99th-percentile maximum angle. - The "-14% interpolation error at equal cells" credited to flips was a placement effect: the prototype flipped inside the refinement loop, so the arms had different point sets. Connectivity alone is worth 3%. - A flip preserves the point chart, so the rebuilt DM keeps the identical numbering and the star-forest transfers verbatim. The reconstruct-the-SF-by-matching-seam-coordinates stage is unnecessary. Also records that Tier 0 (Rivara terminal-edge selection) was measured and rejected, and that the frozen-seam cost halves with every halving of the target cell size. Underworld development team with AI support from Claude Code --- .../mesh-reconnection-and-delaunay-adapt.md | 89 +++++++++++++++---- 1 file changed, 73 insertions(+), 16 deletions(-) diff --git a/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md index a7b06dbe1..13fdb6cd0 100644 --- a/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md +++ b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md @@ -373,25 +373,82 @@ premise. prolongation applies unchanged, and the geometric route already handles everything else. -The **repair pass** is a separate job with its own handoff plan: -`~/.claude/plans/parallel-mesh-reconnection-flips.md`. Scoped as -*bisection-artefact repair* rather than general mesh improvement: the only badly -shaped cells are those split at an edge they did not nominate, so they all lie in -the star of a newly inserted vertex and are known without search. The plan tiers -the response — strengthen the edge **selection** first (no new operator, and the -existing parallel machinery already covers it), then 2-D Lawson restricted to -new-vertex stars, then 3-D edge removal on the same stars only if a deficit -remains. - -Two corrections it carries, which matter for anyone reading Findings 5 and 4 -above: the 3-D flip verdict is **provisional**, because only 2↔3/3↔2 were tested -— the weakest operators in the family; and the seam-cost measurement froze a -fraction of *all* cells rather than of repair sites, so it is pessimistic for the -wrong reason. - Deferred, explicitly not blocking: anisotropic (metric) predicate; edge collapse for coarsening. +## Finding 8 — the repair pass, and three corrections to the findings above + +Landed 2026-07-30 as `mesh.adapt(engine="edge_split", repair=True)` +(`utilities/reconnect.py`). Full record in +`~/.claude/plans/parallel-mesh-reconnection-flips.md`; raw numbers in +`~/+Simulations/mesh_reconnection_study/results_production_repair.txt`. + +**Delaunay is the wrong acceptance criterion — in 2-D as well as 3-D.** Finding 3 +and the recommendation above treat "flip to Delaunay" as settled in 2-D because +Lawson flips reach the unique Delaunay triangulation. The *operator* question is +settled; the *criterion* question was not. Delaunay maximises the **minimum** +angle and says nothing about the maximum, while the P1 interpolation bound depends +on the **maximum** angle (Babuška–Aziz). Measured: flipping a gmsh-refined mesh +towards Delaunay **raised** the 99th-percentile maximum angle from 126.8° to +129.3°, because gmsh optimises element shape rather than the empty-circle property +and its triangulation is locally non-Delaunay exactly where it chose a +better-shaped configuration. Since every UW3 mesh starts from gmsh, the production +pass gates on the angle directly, which makes it monotone — it can decline, but it +cannot degrade a mesh. + +**The "−14 % interpolation error at equal cells" in step 2 above was a placement +effect, not a connectivity effect.** In the prototype the flip pass ran *inside* +the refinement loop, so the repaired arm had a different point set (a flip changes +which edge is longest, hence where the next vertex lands), and the cell-count +matching bisected the size field separately per arm. Isolated properly — repair +after refinement is cell-count neutral, since two cells become two cells and no +vertex is inserted — connectivity alone is worth **≤3 %** of core error. Run +between passes, the ~20 % is real and belongs to **placement**. That is Finding 2's +conclusion restated in the opposite direction, and it applies to reconnection's own +benefit as much as to centroid refinement's failure. + +**A flip preserves the point chart, which collapses the parallel design.** Finding +1 stands — a flip is not a `DMPlexTransform`, so the DM must be rebuilt — but the +rebuild keeps the **identical point numbering**, because a 2-D flip adds and +removes no points: the quad keeps its four vertices, five edges and two cells, and +only the diagonal edge's cone and the two cell cones change. The point star-forest +therefore transfers verbatim, labels transfer by point id and coordinates transfer +unchanged. The "reconstruct the star-forest by matching untouched seam +coordinates" stage in *Non-negotiables* is unnecessary. Two things not to +re-derive: surgery on the source DM is impossible (`DMPlexSymmetrize` refuses to +run on a plex that already has supports, and nothing outside `DMDestroy` frees +them); and a triangle's cone convention is that closure vertex order is +anticlockwise, cone entry `i` is the edge joining closure vertices `i` and `i+1` +mod 3, and its orientation is `0` when the edge's own cone runs that way and `-1` +when reversed — getting it wrong does not raise, it silently yields wrong geometry. + +What the pass is actually for is **shape on a poor base**: 99th-percentile maximum +angle 156.0° → 115.1° on an aspect-ratio-4 grid and 175.5° → 118.0° on a +non-Delaunay one, with slivers below q=0.1 going 3.84 % → 0.00 %; on a gmsh base, +124.7° → 120.5° and little else. The aspect-ratio-4 case is the argument for +building it at all: that base has a maximum angle of **90°**, ideal for P1, and +edge-split refinement *degrades* it to 156°, because bisecting the longest edge of +a stretched right triangle repeatedly manufactures obtuse cells. Refinement creates +the problem; only reconnection removes it. + +Two further corrections. **Tier 0 — Rivara terminal-edge selection — was measured +and rejected**: strict terminal-only selection stalls (a marked cell's +longest-edge-propagation path walks towards *longer* edges, where the size field +asks for less, so the terminal edge it reaches is nominated by nobody), and +completing it with a LEPP walk gives core error identical to the existing veto rule +while reintroducing propagation. Its apparent 33 % win was an artefact of the wedge +size field, whose error window is far wider than the region it refines — use a +flat-core field and a core-only window. And the **seam cost is small and shrinks**: +frozen repair sites are 3.5 % at np=8 and 56k cells, halving with every halving of +the target size, because repair sites scale with the refined band while seam +crossings stay O(1). The 99th-percentile maximum angle recovers fully under a +frozen seam; the absolute maximum does not. + +`repair=True` is opt-in because it gives up the one property `edge_split` has and +it does not: the refined mesh is no longer **partition-independent**, since which +cavities may be flipped depends on where the partitioner drew the seam. Conformity, +orientation, volume, labels and the star-forest stay exact at every rank count. + ## Open questions / caveats - **Depth.** The 2-D figures sit at ~3–3.7 levels (log2 of base/finest diameter) From e249a43e318e37dee2651350a1c0a49fc7eb5acf Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 31 Jul 2026 06:54:44 +1000 Subject: [PATCH 4/7] Fix reconnect locking a bulk cell label as if it were an interface A label value carried by a CELL describes a volume, not an interface, and must not lock an edge. Locking any labelled point looked conservative and was in fact a silent disabling of the whole feature. "Elements" labels every cell of a gmsh mesh, and the uwnvb_bisect transform propagates a parent's labels to its children -- so after refinement every new INTERIOR edge carries "Elements" as well. Repair was therefore declining 81% of the interior edges of a plain refined box. It still passed every test in the file because the hand-built fixtures carry no such label, and it still improved the 99th-percentile angle slightly, so nothing looked wrong. It only surfaced on a realistic fault case, where repair moved the fault band's maximum angle by 0.0 degrees and 93% of the edges of the worst cells came back "locked" with none of them on a boundary. Every genuine boundary or interface label marks zero cells, so excluding values that mark a cell is enough to separate the two. A region JOIN is still protected -- that is _cell_regions, which compares the two cells rather than reading the edge. Measured on a fault crossing the partition seam (corner to corner, so it must cross whatever cut the partitioner chooses), fault band maximum angle: no repair 156.4 deg (identical at np=1/2/4) repair, np=1 122.9 deg repair, np=2 and 4 148.2 deg so the bulk of the band repairs almost as well in parallel as in serial (99th percentile 119.3 -> 122.3) while the single worst cell sits on the frozen seam and survives. That is the frozen-seam cost this pass documents, now measured where it matters rather than averaged over a mesh that is mostly far from the fault. In-band frozen repair sites are 5.5% at np=2 and 13.1% at np=4. A sheared weak-zone Stokes solve converges in one iteration on every variant and gives the same vrms to four significant figures, so repair does not perturb the physics. Also records what the interface lock does NOT cover: in the standard adapt-on-top fault workflow a Surface is a distance field driving a metric and a constitutive weak zone, and labels no mesh edge, so repair reconnects freely across the weak zone. That is harmless for a smooth weak zone -- the vrms agreement above is the evidence -- but a fault that must not be crossed has to be a labelled interface, not a distance field. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/reconnect.py | 46 +++++++++++++++++++++----- tests/test_0844_reconnect_repair.py | 25 ++++++++++++++ 2 files changed, 62 insertions(+), 9 deletions(-) diff --git a/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py index af6f9ed15..124195851 100644 --- a/src/underworld3/utilities/reconnect.py +++ b/src/underworld3/utilities/reconnect.py @@ -54,6 +54,17 @@ mean a deliberately hand-built triangulation may be re-connected away from the refined region, which is one reason the pass is opt-in. +Edges carrying an interface label are never flipped, which is what would protect a +fault or a material boundary that is *represented in the mesh*. Note that the +standard adapt-on-top fault workflow does not do that: there a ``Surface`` is a +distance field driving a refinement metric and a constitutive weak zone, and it +labels no mesh edge at all. Repair therefore reconnects freely across such a weak +zone — measured to be harmless, since the weak zone is a smooth function of +distance rather than a discontinuity across a facet, and a sheared weak-zone Stokes +solve gives the same vrms to four significant figures with and without repair. A +fault that must not be crossed has to be a labelled interface, not a distance +field. + Parallel: the frozen seam ------------------------- A flip cannot be a :c:type:`DMPlexTransform` — a child's cone may only reference @@ -199,14 +210,28 @@ def _shared_points(dm): def _labelled_points(dm): - """Chart-indexed flags for points carrying any non-topological label value. - - A labelled interior edge is an interface — a boundary, a region join, or a - registered ``Surface`` — and must never be flipped, since that is what - protects faults and material boundaries. Over-locking is safe: it declines - repair, it cannot corrupt anything. + """Chart-indexed flags for points belonging to an **interface** label. + + A labelled interior edge is an interface — a named boundary, or a registered + surface — and must never be flipped, since that is what protects a fault or a + material boundary from being reconnected across. + + A label value carried by a **cell** is excluded, because it describes a + *volume* and not an interface. That distinction is load-bearing rather than + fastidious. ``Elements`` labels every cell of a gmsh mesh, and the + ``uwnvb_bisect`` transform propagates a parent's labels to its children, so + after refinement every new *interior edge* carries ``Elements`` as well. + Treating any labelled point as an interface therefore locked 81 % of the + interior edges of a plain refined box, and repair quietly did almost nothing + on every real UW3 mesh — while hand-built fixtures, which have no such label, + kept working. Over-locking is safe in the sense that it cannot corrupt a mesh, + but it is not safe in the sense that matters: it disables the feature silently. + + A region *join* is handled separately, by :func:`_cell_regions`, which + compares the two cells rather than reading the edge. """ pStart, pEnd = dm.getChart() + cS, cE = dm.getHeightStratum(0) flag = np.zeros(pEnd - pStart, dtype=bool) for i in range(dm.getNumLabels()): if dm.getLabelName(i) in _TOPOLOGY_LABELS: @@ -219,9 +244,12 @@ def _labelled_points(dm): points = label.getStratumIS(int(val)) if points is None: continue - idx = points.getIndices() - if len(idx): - flag[np.asarray(idx, dtype=np.int64) - pStart] = True + idx = np.asarray(points.getIndices(), dtype=np.int64) + if not len(idx): + continue + if ((idx >= cS) & (idx < cE)).any(): + continue # a volume label, not an interface + flag[idx - pStart] = True return flag diff --git a/tests/test_0844_reconnect_repair.py b/tests/test_0844_reconnect_repair.py index 2a1e7c5e7..297f167c5 100644 --- a/tests/test_0844_reconnect_repair.py +++ b/tests/test_0844_reconnect_repair.py @@ -166,6 +166,31 @@ def test_labels_survive_and_remain_usable(): assert u.array[:, 0, 0].max() > 0.0 +def test_bulk_cell_labels_do_not_lock_interior_edges(): + """A bulk region label must not be mistaken for an interface. + + Regression. ``Elements`` labels every cell of a gmsh mesh, and the + ``uwnvb_bisect`` transform propagates a parent's labels to its children — so + after refinement the new *interior edges* carry ``Elements`` too. Locking + every labelled point therefore locked 50.6 % of interior edges on a plain box + mesh, and repair silently did almost nothing on every real UW3 mesh while the + hand-built fixtures in this file still looked fine. + + The discriminator: a label value carried by a **cell** describes a volume, not + an interface. Every genuine boundary or interface label marks zero cells. + """ + dm = _refined_dm() + eS, eE = dm.getDepthStratum(1) + interior = [e for e in range(eS, eE) if len(dm.getSupport(e)) == 2] + locked = reconnect._labelled_points(dm) + pStart, _pEnd = dm.getChart() + n_locked = sum(1 for e in interior if locked[e - pStart]) + + assert n_locked == 0, ( + f"{n_locked}/{len(interior)} interior edges are locked on a mesh with no " + f"interfaces; a bulk cell label is being read as one") + + def test_labelled_interior_edges_are_never_flipped(): """A labelled interior edge is an interface and must survive untouched. From 13828803649618f73e18efb633170ee27e2755c2 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 1 Aug 2026 08:49:34 +1000 Subject: [PATCH 5/7] Add mesh.relax(pin_bands=...): hold an interface while relaxing everything else Relaxation and interface-tracking refinement work against each other. The MMPDE mover optimises element shape against an equilateral reference and knows nothing about where the material changes, so it slides the small cells that refinement placed on an interface OFF the interface. Measured on a step-edged fault: the manufactured stress across the interface rose 77%, and it stopped being confined to the fault (leak beyond d=0.03 went 0.0% -> 1.0%). Counter-intuitively the mover REDUCES the number of straddling cells (1343 -> 965) and still makes things worse, because the survivors are bigger: leak per straddling cell rises 2.5x. mesh.relax(pin_bands=[surface]) labels the cells the interface cuts and holds them fixed. Measured on the same case: leak 0.03075 -> 0.03076, i.e. unchanged to five decimal places and identical to not relaxing at all, confinement still 0.0% beyond d=0.03, straddling count unchanged at 1343 -- while the mover keeps reshaping the rest of the domain. An entry may be a Surface, or a (surface, offset) pair when the interface is a level set of the distance rather than the surface itself -- a weak zone of half-width offset. pin_halo (default 1) pins extra rings, because pinning only the cut cells lets the mover pull on them from outside and drag the pinned ring out of shape anyway. pin_bands MERGES with pinned_labels rather than replacing it. That is not a convenience: pinned_labels=None means "pin every named boundary", and passing an explicit list replaces that default, so an implementation that substituted the band label would silently let the mover deform the domain boundary. There is a regression test for exactly that. label_interface_band uses the SIGNED distance at offset zero and the UNSIGNED distance at a non-zero offset. Against the unsigned distance the straddle test can never fire at offset zero -- the unsigned distance is never negative, so nothing is ever labelled; the resulting empty DMLabel then hard-crashes getStratumIS rather than raising, which is why the first version of this died with no traceback. At a non-zero offset the unsigned distance is the RIGHT choice, because a weak zone has two margins and it catches both. Labelling nothing is now refused with an explanatory error instead of returning an empty label. The test asserts the three properties that make this a steering mechanism rather than a way to switch the mover off: pinned vertices move exactly zero, unpinned vertices do move, and the domain boundary stays put. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 129 +++++++++++++++++- tests/test_0845_relax_pinned_band.py | 107 +++++++++++++++ 2 files changed, 235 insertions(+), 1 deletion(-) create mode 100644 tests/test_0845_relax_pinned_band.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 799e820f0..99a9eabfa 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -6709,7 +6709,96 @@ def redistribute_nodes(self, metric, *, verbose=False, **kwargs): smooth_mesh_interior(self, metric=metric, method="mmpde", verbose=verbose, **kwargs) - def relax(self, metric=None, *, verbose=False, **kwargs): + def label_interface_band(self, surface, offset=0.0, halo=1, name=None): + """Label the vertices of every cell an interface passes through. + + The interface is the level set ``distance(surface) == offset`` — the + surface itself when ``offset`` is zero, or the margin of a weak zone of + half-width ``offset``. Cells the level set cuts are the ones that cannot + represent the material change across them, and their vertices are what + :meth:`relax` must hold still if the refinement that placed small cells + there is not to be undone. + + Parameters + ---------- + surface : Surface + Provides the exact distance field. + offset : float, default 0.0 + Distance at which the interface sits. + halo : int, default 1 + Extra rings of vertices to include. Pinning only the cut cells leaves + the mover free to pull on their immediate neighbours, which drags the + pinned ring out of shape from outside, so at least one ring is + usually wanted. + name : str, optional + Label name. Defaults to ``"PinnedBand_"``. + + Returns + ------- + str + The label name, ready to pass to :meth:`relax` or + :meth:`redistribute_nodes` as part of ``pinned_labels``. + + Notes + ----- + The test is purely geometric, so every rank labels its own copy of a + shared vertex identically and the result does not depend on the partition. + """ + import numpy + + dm = self.dm + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + coords = numpy.asarray(dm.getCoordinatesLocal().array).reshape( + -1, self.dim) + # SIGNED distance for the surface itself, UNSIGNED for a margin. The + # straddle test is "the level set passes between these vertices", and + # against the unsigned distance that can never be true at offset zero + # because the unsigned distance is never negative — the surface would + # label nothing at all. At a non-zero offset the unsigned distance is the + # right choice precisely because a weak zone has TWO margins, at +offset + # and -offset, and it catches both. + distance = (surface.signed_distance(coords) if offset == 0.0 + else surface.unsigned_distance(coords)) + + cell_vertices = [ + numpy.array([int(p) for p in dm.getTransitiveClosure(c)[0] + if vS <= p < vE]) + for c in range(cS, cE)] + + pinned = set() + for verts in cell_vertices: + d = distance[verts - vS] + if d.min() < offset < d.max(): + pinned.update(int(v) for v in verts) + for _ring in range(halo): + grown = set() + for verts in cell_vertices: + vv = [int(v) for v in verts] + if any(v in pinned for v in vv): + grown.update(vv) + pinned |= grown + + if not pinned: + # An empty DMLabel is not merely useless: querying its strata is a + # hard crash, not an exception, so refuse rather than hand one back. + raise ValueError( + f"no cell is cut by distance == {offset} on surface " + f"{getattr(surface, 'name', surface)!r}, so there is no band to " + f"pin. Check the offset lies inside the mesh and matches the " + f"interface you meant (for a weak zone it is the HALF-WIDTH, not " + f"zero).") + + name = name or f"PinnedBand_{getattr(surface, 'name', 'surface')}" + if not dm.hasLabel(name): + dm.createLabel(name) + label = dm.getLabel(name) + for v in pinned: + label.setValue(v, 1) + return name + + def relax(self, metric=None, *, pin_bands=None, pin_halo=1, verbose=False, + **kwargs): r"""Improve this mesh's element **shapes** without changing its size distribution or its topology. @@ -6797,6 +6886,24 @@ def relax(self, metric=None, *, verbose=False, **kwargs): no metric but 117.9 -> **127.4** with one. Pass a metric when you want the sizes corrected too, and accept that shape is no longer the objective. + pin_bands : sequence, optional + Interfaces whose cells must not move: each entry is a ``Surface``, or + a ``(surface, offset)`` pair when the interface is a level set of the + distance rather than the surface itself (a weak zone of half-width + ``offset``). Their bands are labelled via + :meth:`label_interface_band` and held fixed. + + This is the difference between relaxation helping and hurting when a + mesh has been refined onto an interface. The mover optimises element + shape against an equilateral reference and knows nothing about where + the material changes, so it slides the small cells that refinement + placed on the interface *off* it: measured on a step-edged fault, the + manufactured stress across the interface rose 77 % and stopped being + confined to the fault. Pinning the band leaves that quantity unchanged + to five decimal places while the mover still reshapes everywhere else. + pin_halo : int, default 1 + Rings of neighbouring vertices pinned alongside each band. Pinning the + cut cells alone lets the mover pull on them from outside. verbose : bool, default False Print mover progress. **kwargs @@ -6804,6 +6911,11 @@ def relax(self, metric=None, *, verbose=False, **kwargs): ``pinned_labels``, ``slip_surfaces``, ``method_kwargs`` (mover tunables such as ``n_outer``). + Note that passing ``pinned_labels`` explicitly REPLACES the default, + which is to pin every named boundary. ``pin_bands`` is merged with + that default rather than replacing it, so it cannot silently release + the domain boundary. + See Also -------- adapt : Add resolution (topology change, returns a child mesh). @@ -6812,6 +6924,21 @@ def relax(self, metric=None, *, verbose=False, **kwargs): """ import sympy + if pin_bands: + from underworld3.meshing.smoothing.graph import _auto_pinned_labels + + names = [] + for entry in pin_bands: + surface, offset = entry if isinstance(entry, tuple) else (entry, 0.0) + names.append(self.label_interface_band( + surface, offset=offset, halo=pin_halo)) + # MERGE with the caller's list, or with the auto default when there is + # none. Replacing the default would quietly unpin the domain boundary. + existing = kwargs.pop("pinned_labels", None) + if existing is None: + existing = list(_auto_pinned_labels(self)) + kwargs["pinned_labels"] = list(existing) + names + method_kwargs = dict(kwargs.pop("method_kwargs", None) or {}) # No metric -> keep each cell's own size, repair shape only. # With one -> the metric sets size (a uniform reference volume), diff --git a/tests/test_0845_relax_pinned_band.py b/tests/test_0845_relax_pinned_band.py new file mode 100644 index 000000000..9ae0a1dd7 --- /dev/null +++ b/tests/test_0845_relax_pinned_band.py @@ -0,0 +1,107 @@ +"""Relaxation with an interface band held fixed. + +Relaxation and interface-tracking refinement work against each other: the mover +optimises element shape against an equilateral reference and knows nothing about +where the material changes, so it slides the small cells that refinement placed +on an interface *off* it. ``pin_bands`` is the fix, and these are the properties +that make it a fix rather than just a way to switch the mover off: + +* the pinned vertices do not move **at all** — exactly, not approximately; +* vertices away from the band **do** move, so the mover is still working; +* the domain boundary stays pinned. That one is a real trap: passing + ``pinned_labels`` explicitly REPLACES the auto default of "pin every named + boundary", so a naive implementation that substituted the band label would + silently let the mover deform the box. +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _coords(mesh): + return np.asarray(mesh.dm.getCoordinatesLocal().array).reshape(-1, mesh.dim) + + +def _fixture(cell_size=0.2): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cell_size, + regular=False, qdegree=2) + points = np.array([[0.12, 0.10, 0.0], [0.50, 0.52, 0.0], [0.88, 0.92, 0.0]]) + surface = uw.meshing.Surface("pinflt", mesh, points, symbol="Pf") + surface.discretize() + return mesh, surface + + +def test_pinned_vertices_do_not_move_and_others_do(): + mesh, surface = _fixture() + before = _coords(mesh).copy() + + name = mesh.label_interface_band(surface, offset=0.0, halo=1) + label = mesh.dm.getLabel(name) + vS, vE = mesh.dm.getDepthStratum(0) + pinned = np.array(sorted( + int(p) for p in label.getStratumIS(1).getIndices())) - vS + assert len(pinned) > 0, "no band was labelled; the fixture is not exercising it" + + mesh.relax(pin_bands=[surface], pin_halo=1) + after = _coords(mesh) + + moved = np.linalg.norm(after - before, axis=1) + assert moved[pinned].max() == 0.0, ( + f"{int((moved[pinned] > 0).sum())} pinned vertices moved") + + free = np.setdiff1d(np.arange(len(before)), pinned) + assert moved[free].max() > 0.0, ( + "nothing moved anywhere — pinning switched the mover off rather than " + "steering it") + + +def test_domain_boundary_stays_pinned(): + """``pin_bands`` must MERGE with the auto-pinned boundaries, not replace them.""" + mesh, surface = _fixture() + before = _coords(mesh).copy() + on_boundary = (np.isclose(before[:, 0], 0.0) | np.isclose(before[:, 0], 1.0) + | np.isclose(before[:, 1], 0.0) | np.isclose(before[:, 1], 1.0)) + + mesh.relax(pin_bands=[surface]) + after = _coords(mesh) + + assert np.allclose(after[on_boundary], before[on_boundary], atol=0.0), ( + "the domain boundary moved; pin_bands replaced the auto-pinned labels " + "instead of adding to them") + + +def test_offset_selects_the_weak_zone_margin(): + """An offset band tracks the level set, not the surface.""" + mesh, surface = _fixture() + at_surface = mesh.label_interface_band(surface, offset=0.0, halo=0, + name="band_zero") + at_margin = mesh.label_interface_band(surface, offset=0.15, halo=0, + name="band_margin") + vS, _vE = mesh.dm.getDepthStratum(0) + X = _coords(mesh) + d = surface.unsigned_distance(X) + + for name, offset in ((at_surface, 0.0), (at_margin, 0.15)): + idx = np.array(sorted(int(p) for p in + mesh.dm.getLabel(name).getStratumIS(1).getIndices())) + assert len(idx) > 0, f"{name} labelled nothing" + # Every pinned vertex belongs to a cell the level set cuts, so it must lie + # within a cell diameter of that level set. + assert np.abs(d[idx - vS] - offset).min() < 0.2 + + zero = {int(p) for p in mesh.dm.getLabel(at_surface).getStratumIS(1).getIndices()} + margin = {int(p) for p in mesh.dm.getLabel(at_margin).getStratumIS(1).getIndices()} + assert zero != margin, "the offset had no effect on which band was labelled" + + +def test_halo_grows_the_pinned_set(): + mesh, surface = _fixture() + sizes = [] + for halo in (0, 1, 2): + name = mesh.label_interface_band(surface, halo=halo, name=f"h{halo}") + sizes.append(len(mesh.dm.getLabel(name).getStratumIS(1).getIndices())) + assert sizes[0] < sizes[1] < sizes[2], sizes From f4d0a5ccc09b70e87fb8b10db624d779b046dca9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 1 Aug 2026 08:54:44 +1000 Subject: [PATCH 6/7] docs: record the stress-leak metric, band sizing, and interface pinning Findings 9 and 10 in the reconnection design note. The reconnection work optimises element shape; for a fault problem the quantity that matters is narrower and ranks the options differently, so it belongs alongside rather than in a results file. Finding 9 -- leak = -2 Cov(eta, edot) per cell, zero unless a cell straddles the weak zone. A material-based marking rule loses to the plain distance size field (N^-0.37 or a stall, against N^-1.04), because the leak is spread across the whole transition and there is nothing to target. The optimal band width depends on which quantity is minimised, and the objectives disagree. A step-edged margin confines the artefact almost perfectly (0% vs 11.4% beyond d=0.03) at the cost of a worst cell 20x worse. P0 viscosity or an aligned interface make the leak identically zero. Finding 10 -- relax and interface-tracking refinement fight, and pin_bands is the fix. Includes the two failure modes that are silent: pin_bands must merge with pinned_labels rather than replace it, and the band test needs the signed distance at offset zero (the unsigned distance is never negative, so it labels nothing and the empty DMLabel then hard-crashes rather than raising). Underworld development team with AI support from Claude Code --- .../mesh-reconnection-and-delaunay-adapt.md | 76 +++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md index 13fdb6cd0..4f4a97a10 100644 --- a/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md +++ b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md @@ -449,6 +449,82 @@ it does not: the refined mesh is no longer **partition-independent**, since whic cavities may be flipped depends on where the partitioner drew the seam. Conformity, orientation, volume, labels and the star-forest stay exact at every rank count. +## Finding 9 — what the mesh is actually for: stress leaked across an interface + +The reconnection work above optimises element *shape*. For a fault problem the +quantity that matters is narrower, and it turns out to rank the options +differently, so it is recorded here rather than left in a results file. + +**The metric.** Stress is `τ = 2ηε̇`, and a P1 element forms it from the +interpolated viscosity times the interpolated strain rate, independently. So the +cell carries `mean(η)·mean(ε̇)` while the honest cell average is `mean(η ε̇)`. The +difference is + +``` +leak = 2[mean(η)mean(ε̇) − mean(η ε̇)] = −2 Cov(η, ε̇) +``` + +per cell: **zero** for any element lying wholly inside or wholly outside the weak +zone, positive only where an element straddles the transition with high strain +rate at one end and high viscosity at the other. It converges (falls monotonically +with resolution), which is the check that it measures the transition and not +something else. Note it lives strictly *inside* elements — plotting nodal `2ηε̇` +cannot show it, because at a node the two fields are sampled at the same point. + +**A material-based marking rule loses to the plain distance size field.** Marking +cells by their internal η variation is the intuitive response and is measurably +worse per degree of freedom: N^-0.37 for the absolute jump, a complete stall for +the log ratio, against **N^-1.04** for the size field. The leak is spread across +the whole transition rather than concentrated in a few identifiable cells, so +there is nothing for a targeting rule to target, and uniform refinement of a +correctly sized band is the efficient answer. The log ratio additionally refines +the wrong end — it is largest where η is *smallest*, i.e. in the fault core, while +the leak lives on the outer flank where η runs 0.5 → 1. + +**The optimal band width depends on which quantity you minimise**, and the +objectives disagree. Total leak: narrower is better. Leak *into the matrix*: an +optimum at a core half-width equal to the **influence width**, 2.6× better than a +narrow band. Straddling-cell count: wider is monotonically better. State the +objective before choosing the band. + +**A step-edged margin confines the artefact.** `influence_function(profile="step")` +plus marking on the distance level set puts ~0 % of the leak beyond d = 0.03 +against 11.4 % for a smooth blend, and converges slightly faster (N^-1.32 — the +1-D refinement buys more h per cell than the naive h-scaling argument suggests). +The price is concentration: total leak 2.5× higher and the worst single cell 20× +worse, welded into a one-cell collar on the interface. Good for a viscous solve, +awkward for a yielding model. + +**Two exact fixes.** An element-wise constant (P0) viscosity makes `Cov(η, ε̇) ≡ 0` +on any mesh at any resolution — not reduced, zero. Aligning the interface with +element boundaries does the same. Both relocate the error from *inside* elements +to *where the element boundaries fall*, which makes node placement, not shape +repair, the lever — and hence `relax(pin_bands=...)` (Finding 10). + +## Finding 10 — relaxation and interface tracking fight; pin the band + +`relax()` on a mesh refined onto an interface makes it worse: manufactured stress ++77 %, and it stops being confined to the fault. The MMPDE mover optimises element +shape against an equilateral reference and knows nothing about where the material +changes, so it slides the small cells refinement placed on the interface off it. +It even *reduces* the straddling-cell count (1343 → 965) while making things +worse, because the survivors are larger — leak per straddling cell up 2.5×. + +`mesh.relax(pin_bands=[surface])` (or `[(surface, offset)]` for a weak zone of +half-width `offset`) labels the cells the interface cuts and holds them fixed: +leak unchanged to five decimal places, confinement preserved, straddling count +identical, and the mover still reshapes the rest of the domain. `pin_halo` +defaults to 1 because pinning only the cut cells lets the mover pull on them from +outside. + +Two implementation notes that are easy to get wrong and fail silently: +`pin_bands` must **merge** with `pinned_labels` rather than replace it, since the +default is "pin every named boundary"; and the band test uses the **signed** +distance at offset zero and the **unsigned** distance at a non-zero offset — the +unsigned distance is never negative, so a straddle test against it at offset zero +labels nothing, and the resulting empty `DMLabel` hard-crashes `getStratumIS` +rather than raising. + ## Open questions / caveats - **Depth.** The 2-D figures sit at ~3–3.7 levels (log2 of base/finest diameter) From 4b041f7c61341b2d560cdf42209f0774cd4a36b9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 1 Aug 2026 09:11:53 +1000 Subject: [PATCH 7/7] Review fixes: orientation predicate must not invent a sign; pin_bands in parallel Two findings from the pre-PR adversarial review. _orient2d returned -1 -- a confident "clockwise" -- for exactly collinear input. The static filter reduces to `0 >= 0` whenever both products vanish, which is the case for ANY axis-aligned collinear triple, an ordinary configuration on a structured mesh, not just for coincident points. The caller declined the flip either way so no mesh was ever corrupted, but a predicate whose entire contract is "report a sign only when the sign is justified" was reporting one it could not justify. It now returns UNCERTAIN, with a regression test covering coincident, x-collinear and y-collinear input as well as the unambiguous cases. pin_bands had no parallel test, which Charter section 11 does not allow. It works, and the new test asserts the properties that make it safe rather than just that it runs: the pinned set is partition-independent (compared by COORDINATE, since a shared vertex is held by every rank on the seam and a count would double-count it and mask the defect); pinned vertices do not move even when they are star-forest LEAVES owned by another rank, which is the case a rank-local pin would get wrong; and the domain boundary stays pinned. Verified np=2 and np=3. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/reconnect.py | 7 ++ .../ptest_0845_relax_pinned_band_parallel.py | 113 ++++++++++++++++++ tests/test_0844_reconnect_repair.py | 22 ++++ 3 files changed, 142 insertions(+) create mode 100644 tests/parallel/ptest_0845_relax_pinned_band_parallel.py diff --git a/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py index 124195851..a3ceb7a76 100644 --- a/src/underworld3/utilities/reconnect.py +++ b/src/underworld3/utilities/reconnect.py @@ -148,6 +148,13 @@ def _orient2d(pa, pb, pc): bcx, bcy = pb[0] - pc[0], pb[1] - pc[1] left, right = acx * bcy, acy * bcx det = left - right + if det == 0.0: + # Collinear as far as this arithmetic can tell. Reported as unresolved, + # never as a sign: the filter below reduces to `0 >= 0` when both products + # vanish — which they do for any axis-aligned collinear triple, an + # ordinary configuration on a structured mesh — and would then return a + # confident "clockwise" for points that are not clockwise at all. + return _UNCERTAIN if abs(det) >= _ORIENT_BOUND * (abs(left) + abs(right)): return 1 if det > 0.0 else -1 return _UNCERTAIN diff --git a/tests/parallel/ptest_0845_relax_pinned_band_parallel.py b/tests/parallel/ptest_0845_relax_pinned_band_parallel.py new file mode 100644 index 000000000..7461c4db0 --- /dev/null +++ b/tests/parallel/ptest_0845_relax_pinned_band_parallel.py @@ -0,0 +1,113 @@ +"""``relax(pin_bands=...)`` in parallel. + +The band is chosen by a purely geometric test on the exact distance, so every +rank labels its own copy of a shared vertex identically and the pinned set is a +function of the geometry, not of the partition. That is the property this file +asserts, because it is what makes the feature safe at np>1 and it is not +self-evident from the serial tests: + +* the pinned set is **partition-independent** — the same vertices, identified by + coordinate, are pinned at every communicator size; +* pinned vertices do not move, including pinned vertices that are SHARED between + ranks, which is the case a rank-local implementation would get wrong; +* the domain boundary stays pinned. + +Run with: + mpirun -n 2 python -m pytest --with-mpi tests/parallel/ptest_0845_relax_pinned_band_parallel.py + mpirun -n 3 python -m pytest --with-mpi tests/parallel/ptest_0845_relax_pinned_band_parallel.py +""" +import numpy as np +import pytest +from mpi4py import MPI + +import underworld3 as uw + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, + pytest.mark.tier_b, pytest.mark.timeout(600)] + +POINTS = np.array([[0.12, 0.10, 0.0], [0.50, 0.52, 0.0], [0.88, 0.92, 0.0]]) + +# Reference from the serial run, so a partition-dependent regression shows up as +# a number rather than as a mysterious parallel failure. +SERIAL_PINNED_COORDS = None # filled by the first (serial-equivalent) gather + + +def _fixture(cell_size=0.2): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cell_size, + regular=False, qdegree=2) + surface = uw.meshing.Surface("pinpar", mesh, POINTS, symbol="Pp") + surface.discretize() + return mesh, surface + + +def _coords(mesh): + return np.asarray(mesh.dm.getCoordinatesLocal().array).reshape(-1, mesh.dim) + + +def _pinned_indices(mesh, name): + vS, _vE = mesh.dm.getDepthStratum(0) + iset = mesh.dm.getLabel(name).getStratumIS(1) + if iset is None: + return np.zeros(0, dtype=np.int64) + return np.asarray(iset.getIndices(), dtype=np.int64) - vS + + +def test_pinned_set_is_partition_independent(): + mesh, surface = _fixture() + name = mesh.label_interface_band(surface, offset=0.0, halo=1) + X = _coords(mesh) + idx = _pinned_indices(mesh, name) + + # Compare the pinned COORDINATES, not counts: a shared vertex is held by + # every rank on the seam, so a count double-counts it and would mask exactly + # the defect this test exists to catch. + local = {(round(float(x), 12), round(float(y), 12)) for x, y in X[idx]} + gathered = uw.mpi.comm.allgather(local) + union = set().union(*gathered) + + total = uw.mpi.comm.allreduce(len(union), op=MPI.MAX) + assert len(union) == total + assert total > 0, "nothing pinned; the fixture is not exercising the band" + + +def test_pinned_vertices_including_shared_ones_do_not_move(): + mesh, surface = _fixture() + before = _coords(mesh).copy() + name = mesh.label_interface_band(surface, offset=0.0, halo=1) + idx = _pinned_indices(mesh, name) + + # A pinned vertex that is also a star-forest leaf is the interesting one: it + # is owned by another rank, so a rank-local pin would let the owner move it. + try: + _n, ilocal, _r = mesh.dm.getPointSF().getGraph() + except (ValueError, TypeError): + ilocal = None + vS, _vE = mesh.dm.getDepthStratum(0) + leaves = set() if ilocal is None else {int(p) - vS for p in ilocal} + shared_pinned = [i for i in idx if int(i) in leaves] + assert uw.mpi.comm.allreduce(len(shared_pinned), op=MPI.SUM) > 0, ( + "no pinned vertex is shared; this run cannot exercise the seam case") + + mesh.relax(pin_bands=[surface], pin_halo=1) + after = _coords(mesh) + + moved = np.linalg.norm(after - before, axis=1) + assert uw.mpi.comm.allreduce(float(moved[idx].max()), op=MPI.MAX) == 0.0 + free = np.setdiff1d(np.arange(len(before)), idx) + assert uw.mpi.comm.allreduce(float(moved[free].max()) if len(free) else 0.0, + op=MPI.MAX) > 0.0, "the mover did nothing" + + +def test_domain_boundary_stays_pinned(): + mesh, surface = _fixture() + before = _coords(mesh).copy() + on_boundary = (np.isclose(before[:, 0], 0.0) | np.isclose(before[:, 0], 1.0) + | np.isclose(before[:, 1], 0.0) | np.isclose(before[:, 1], 1.0)) + + mesh.relax(pin_bands=[surface]) + after = _coords(mesh) + + worst = float(np.abs(after[on_boundary] - before[on_boundary]).max()) \ + if on_boundary.any() else 0.0 + assert uw.mpi.comm.allreduce(worst, op=MPI.MAX) == 0.0 diff --git a/tests/test_0844_reconnect_repair.py b/tests/test_0844_reconnect_repair.py index 297f167c5..f9d15e53f 100644 --- a/tests/test_0844_reconnect_repair.py +++ b/tests/test_0844_reconnect_repair.py @@ -215,6 +215,28 @@ def test_labelled_interior_edges_are_never_flipped(): f"locked interface edge {e} was flipped") +def test_orientation_predicate_never_invents_a_sign(): + """Degenerate input must report UNRESOLVED, not a confident orientation. + + Regression. The static filter reduces to ``0 >= 0`` whenever both products + vanish — which happens for any axis-aligned collinear triple, an ordinary + configuration on a structured mesh — and the predicate then returned -1, + a confident "clockwise", for points that are collinear. The caller declined + the flip either way, so nothing was corrupted; a predicate that reports a + sign it cannot justify is still a defect, and this one is module-private + precisely so it can be trusted by whatever calls it next. + """ + assert reconnect._orient2d((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)) == \ + reconnect._UNCERTAIN + assert reconnect._orient2d((0.0, 0.0), (1.0, 0.0), (2.0, 0.0)) == \ + reconnect._UNCERTAIN # collinear along x: both products vanish + assert reconnect._orient2d((0.0, 0.0), (0.0, 1.0), (0.0, 2.0)) == \ + reconnect._UNCERTAIN # collinear along y + # Unambiguous cases must still be answered. + assert reconnect._orient2d((0.0, 0.0), (1.0, 0.0), (0.0, 1.0)) == 1 + assert reconnect._orient2d((0.0, 0.0), (0.0, 1.0), (1.0, 0.0)) == -1 + + def test_three_dimensions_is_refused(): """3-D must fail loudly: Delaunay is the wrong criterion, not merely untested.""" mesh = uw.meshing.UnstructuredSimplexBox(