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 00000000..4f4a97a1 --- /dev/null +++ b/docs/developer/design/mesh-reconnection-and-delaunay-adapt.md @@ -0,0 +1,549 @@ +# 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. + +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. + +## 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) + 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 a032ca6b..99a9eabf 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), @@ -6830,7 +6957,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. @@ -6915,12 +7042,43 @@ 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`. + 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 @@ -6990,18 +7148,35 @@ 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}") + 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 @@ -7406,6 +7581,99 @@ 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) + 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(), + 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 0d7c627e..3681d1f3 100644 --- a/src/underworld3/utilities/__init__.py +++ b/src/underworld3/utilities/__init__.py @@ -94,3 +94,5 @@ def _append_petsc_path(): from . import boundary_flux from . import custom_mg from .custom_mg import set_custom_fmg +from . import edge_split +from . import reconnect diff --git a/src/underworld3/utilities/edge_split.py b/src/underworld3/utilities/edge_split.py new file mode 100644 index 00000000..877ddb91 --- /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/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py new file mode 100644 index 00000000..a3ceb7a7 --- /dev/null +++ b/src/underworld3/utilities/reconnect.py @@ -0,0 +1,570 @@ +"""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. + +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 +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 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 + + +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 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: + 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 = 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 + + +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_0843_edge_split_parallel.py b/tests/parallel/ptest_0843_edge_split_parallel.py new file mode 100644 index 00000000..9f6be15b --- /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/parallel/ptest_0844_reconnect_parallel.py b/tests/parallel/ptest_0844_reconnect_parallel.py new file mode 100644 index 00000000..b27780b3 --- /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/parallel/ptest_0845_relax_pinned_band_parallel.py b/tests/parallel/ptest_0845_relax_pinned_band_parallel.py new file mode 100644 index 00000000..7461c4db --- /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_0843_edge_split_adapt.py b/tests/test_0843_edge_split_adapt.py new file mode 100644 index 00000000..37cc831c --- /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) diff --git a/tests/test_0844_reconnect_repair.py b/tests/test_0844_reconnect_repair.py new file mode 100644 index 00000000..f9d15e53 --- /dev/null +++ b/tests/test_0844_reconnect_repair.py @@ -0,0 +1,246 @@ +"""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_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. + + 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_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( + 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) diff --git a/tests/test_0845_relax_pinned_band.py b/tests/test_0845_relax_pinned_band.py new file mode 100644 index 00000000..9ae0a1dd --- /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