Adapt-on-top: closure-free edge_split engine, reconnection repair, and interface-pinned relaxation - #488
Adapt-on-top: closure-free edge_split engine, reconnection repair, and interface-pinned relaxation#488lmoresi wants to merge 7 commits into
Conversation
Newest-vertex bisection picks the edge to split from a combinatorial tagging rule and then pays a conforming closure to repair the hanging nodes that choice creates. This engine splits the edge the geometry asks for -- the longest edge of every cell still coarser than the metric wants -- and needs no closure at all, because splitting an edge divides *every* incident cell at the same new vertex. There is no hanging node to repair and no longest-edge-propagation chain, so refinement cannot escape the marked region: the refined band hugs the feature instead of a halo around it. No new topology code. `uwnvb_bisect` in nvb_transform.c is already a registered DMPlexTransform driven by a per-edge label, works on triangles and tets, and is the primitive NVB uses for each sub-pass. This engine drives it from Python and therefore inherits star-forest propagation, co-partitioning, labels and coordinates for free. Marking is on the cell DIAMETER, not (d! V)^(1/d). For bisection the two shrink together and either will do; for any engine that reduces volume without shortening the longest edge they diverge badly -- a measured factor of 3.2 on a centroid-refined mesh, where the volume proxy reports the target met while the mesh is nowhere near resolved. The test asserts the diameter. Measured: 2-D 104 -> 412 cells in 7 passes and 3-D 1472 -> 3933 tets, identical at np=1/2/3/4 with no over-shared facets; through mesh.adapt, a 10-level graded MG tail with all 8 exact half-half prolongations captured (every inserted vertex is an exact float edge midpoint) and TI Stokes converging in 10 V-cycles. 3-D reaches the pass cap -- an edge is shared by more tets so fewer are independent per pass -- so the budget scales with dimension and warns rather than silently truncating. Selection is a deterministic function of geometry, not of iteration order. A greedy sweep produced a partition-dependent mesh (412/412/463/925 cells at np=1/2/3/4, every one of them conforming and individually plausible), so an edge wins only if it beats every competing candidate sharing a cell, with a midpoint-coordinate tie-break. The parallel test asserts the serial cell count because that class of defect is invisible in a serial run. Also fixes _cells_on_edge, which applied the 3-D edge -> face -> cell walk in both dimensions. In 2-D an edge *is* a face, so it asked for the support of a cell, got nothing, and reported that the edge touches no cells at all. It is not yet called from the engine, so nothing was broken, but it fails silently and the shape-repair work needs it. Underworld development team with AI support from Claude Code
Reconnection is the missing third operation of the refine / swap / smooth triple. UW3 had refine (mesh.adapt) and smooth (mesh.relax); this is swap. Refinement chooses where a vertex goes but not how the surrounding cells reconnect, so a cell dragged into a split at an edge it did not nominate gains a thin child. This repairs that by Lawson flips. The acceptance criterion is NOT Delaunay, although these are Lawson flips. Delaunay maximises the minimum angle and says nothing about the maximum, while the P1 interpolation bound depends on the maximum angle and not the minimum (Babuska-Aziz). The two disagree in practice and not marginally: flipping a gmsh-generated mesh towards Delaunay was measured to RAISE the 99th-percentile maximum angle from 126.8 to 129.3 degrees, because gmsh optimises element shape rather than the empty-circle property and its triangulation is locally non-Delaunay exactly where it chose a better-shaped configuration. Since every UW3 mesh starts from gmsh, a repair pass that can degrade one is unusable. Gating on the angle directly makes the pass monotone: it can decline, but it cannot make a mesh worse. Measured on the production path (numbers in the study directory, see the module docstring). As a post-pass it fixes shape and only shape -- decisively on a poor base (99th-percentile maximum angle 156.0 -> 115.1 degrees on an aspect-ratio-4 base; slivers below q=0.1 3.84% -> 0.00% on a non-Delaunay one) and hardly at all on a gmsh base, with interpolation error barely moving either way. Run between refinement passes it also changes where later vertices land, because a flip changes which edge of a cell is longest, and that is worth 20-30% lower error per degree of freedom on a degraded base. The accuracy gain is therefore a placement gain that reconnection unlocks, not a connectivity gain. Parallel by the frozen seam: no cavity may contain a cell incident on a shared plex point. Measured cost 0.9-3.5% of repair sites at 56k cells and np=2..8, halving with every halving of the target size, because repair sites scale with the refined band while the sites a seam crosses stay O(1). The DM is rebuilt on the SAME point chart. A 2-D flip adds and removes no points -- the quad keeps its four vertices, five edges and two cells, and only the diagonal edge's cone and the two cell cones change -- so preserving the numbering lets the point star-forest transfer verbatim, labels transfer by point id and coordinates transfer unchanged. That removes the whole reconstruct-the-star-forest-by-matching-seam-coordinates stage, and with it the class of defect nvb._exact_vertex_map exists to refuse. Surgery on the source DM is not an option: DMPlexSymmetrize refuses to run on a plex that already has supports and nothing outside DMDestroy frees them. The cone orientation convention is derived from the edge cone every time rather than assumed, because getting it wrong does not raise -- it silently yields wrong geometry. repair is OFF by default, for one specific reason: edge_split alone produces a partition-independent mesh, identical at any communicator size, and repair gives that up, because which cavities may be flipped depends on where the partitioner drew the seam. Conformity, orientation, volume, labels and the star-forest stay exact at every rank count. Also note the 99th-percentile maximum angle recovers fully under a frozen seam but the absolute maximum does not -- a few of the worst cells sit on the seam and are exactly the untouchable ones. Orientation and in-circle sign errors produce non-conforming meshes, so the orientation predicate carries Shewchuk's static filter and DECLINES when it cannot resolve a sign. Declining is always safe here because a flip is an optimisation, never a requirement, which is what lets a filter stand in for adaptive-precision arithmetic inside a refinement loop. Repair invalidates the cell-parent map used by the any-degree nested MG transfer (a flipped cell can straddle two coarse cells), so it is set to None and a degree-2 space falls back to the geometric builder. The exact vertex prolongation survives untouched: flips move no vertex, and a P1 section numbers its DOFs from the point numbering, which is preserved. Tests: 6 serial, 4 parallel at np=2/3/4. The maximum-angle assertion is the one that caught the Delaunay criterion; the idempotence check cannot -- an inverted criterion is idempotent too, which is exactly how Delaunay passed while degrading the mesh. Underworld development team with AI support from Claude Code
…n findings Finding 8 in the reconnection design note. Three of the earlier findings needed correcting rather than extending: - Delaunay is the wrong acceptance criterion in 2-D as well as 3-D. Finding 3 treated it as settled because Lawson flips reach the unique Delaunay triangulation; that settles the operator, not the criterion. Delaunay maximises the minimum angle while P1 interpolation depends on the maximum, and flipping a gmsh mesh towards Delaunay was measured to raise the 99th-percentile maximum angle. - The "-14% interpolation error at equal cells" credited to flips was a placement effect: the prototype flipped inside the refinement loop, so the arms had different point sets. Connectivity alone is worth 3%. - A flip preserves the point chart, so the rebuilt DM keeps the identical numbering and the star-forest transfers verbatim. The reconstruct-the-SF-by-matching-seam-coordinates stage is unnecessary. Also records that Tier 0 (Rivara terminal-edge selection) was measured and rejected, and that the frozen-seam cost halves with every halving of the target cell size. Underworld development team with AI support from Claude Code
A label value carried by a CELL describes a volume, not an interface, and must
not lock an edge. Locking any labelled point looked conservative and was in fact
a silent disabling of the whole feature.
"Elements" labels every cell of a gmsh mesh, and the uwnvb_bisect transform
propagates a parent's labels to its children -- so after refinement every new
INTERIOR edge carries "Elements" as well. Repair was therefore declining 81% of
the interior edges of a plain refined box. It still passed every test in the file
because the hand-built fixtures carry no such label, and it still improved the
99th-percentile angle slightly, so nothing looked wrong. It only surfaced on a
realistic fault case, where repair moved the fault band's maximum angle by 0.0
degrees and 93% of the edges of the worst cells came back "locked" with none of
them on a boundary.
Every genuine boundary or interface label marks zero cells, so excluding
values that mark a cell is enough to separate the two. A region JOIN is still
protected -- that is _cell_regions, which compares the two cells rather than
reading the edge.
Measured on a fault crossing the partition seam (corner to corner, so it must
cross whatever cut the partitioner chooses), fault band maximum angle:
no repair 156.4 deg (identical at np=1/2/4)
repair, np=1 122.9 deg
repair, np=2 and 4 148.2 deg
so the bulk of the band repairs almost as well in parallel as in serial
(99th percentile 119.3 -> 122.3) while the single worst cell sits on the frozen
seam and survives. That is the frozen-seam cost this pass documents, now measured
where it matters rather than averaged over a mesh that is mostly far from the
fault. In-band frozen repair sites are 5.5% at np=2 and 13.1% at np=4. A sheared
weak-zone Stokes solve converges in one iteration on every variant and gives the
same vrms to four significant figures, so repair does not perturb the physics.
Also records what the interface lock does NOT cover: in the standard
adapt-on-top fault workflow a Surface is a distance field driving a metric and a
constitutive weak zone, and labels no mesh edge, so repair reconnects freely
across the weak zone. That is harmless for a smooth weak zone -- the vrms
agreement above is the evidence -- but a fault that must not be crossed has to be
a labelled interface, not a distance field.
Underworld development team with AI support from Claude Code
…thing else Relaxation and interface-tracking refinement work against each other. The MMPDE mover optimises element shape against an equilateral reference and knows nothing about where the material changes, so it slides the small cells that refinement placed on an interface OFF the interface. Measured on a step-edged fault: the manufactured stress across the interface rose 77%, and it stopped being confined to the fault (leak beyond d=0.03 went 0.0% -> 1.0%). Counter-intuitively the mover REDUCES the number of straddling cells (1343 -> 965) and still makes things worse, because the survivors are bigger: leak per straddling cell rises 2.5x. mesh.relax(pin_bands=[surface]) labels the cells the interface cuts and holds them fixed. Measured on the same case: leak 0.03075 -> 0.03076, i.e. unchanged to five decimal places and identical to not relaxing at all, confinement still 0.0% beyond d=0.03, straddling count unchanged at 1343 -- while the mover keeps reshaping the rest of the domain. An entry may be a Surface, or a (surface, offset) pair when the interface is a level set of the distance rather than the surface itself -- a weak zone of half-width offset. pin_halo (default 1) pins extra rings, because pinning only the cut cells lets the mover pull on them from outside and drag the pinned ring out of shape anyway. pin_bands MERGES with pinned_labels rather than replacing it. That is not a convenience: pinned_labels=None means "pin every named boundary", and passing an explicit list replaces that default, so an implementation that substituted the band label would silently let the mover deform the domain boundary. There is a regression test for exactly that. label_interface_band uses the SIGNED distance at offset zero and the UNSIGNED distance at a non-zero offset. Against the unsigned distance the straddle test can never fire at offset zero -- the unsigned distance is never negative, so nothing is ever labelled; the resulting empty DMLabel then hard-crashes getStratumIS rather than raising, which is why the first version of this died with no traceback. At a non-zero offset the unsigned distance is the RIGHT choice, because a weak zone has two margins and it catches both. Labelling nothing is now refused with an explanatory error instead of returning an empty label. The test asserts the three properties that make this a steering mechanism rather than a way to switch the mover off: pinned vertices move exactly zero, unpinned vertices do move, and the domain boundary stays put. Underworld development team with AI support from Claude Code
Findings 9 and 10 in the reconnection design note. The reconnection work optimises element shape; for a fault problem the quantity that matters is narrower and ranks the options differently, so it belongs alongside rather than in a results file. Finding 9 -- leak = -2 Cov(eta, edot) per cell, zero unless a cell straddles the weak zone. A material-based marking rule loses to the plain distance size field (N^-0.37 or a stall, against N^-1.04), because the leak is spread across the whole transition and there is nothing to target. The optimal band width depends on which quantity is minimised, and the objectives disagree. A step-edged margin confines the artefact almost perfectly (0% vs 11.4% beyond d=0.03) at the cost of a worst cell 20x worse. P0 viscosity or an aligned interface make the leak identically zero. Finding 10 -- relax and interface-tracking refinement fight, and pin_bands is the fix. Includes the two failure modes that are silent: pin_bands must merge with pinned_labels rather than replace it, and the band test needs the signed distance at offset zero (the unsigned distance is never negative, so it labels nothing and the empty DMLabel then hard-crashes rather than raising). Underworld development team with AI support from Claude Code
… in parallel Two findings from the pre-PR adversarial review. _orient2d returned -1 -- a confident "clockwise" -- for exactly collinear input. The static filter reduces to `0 >= 0` whenever both products vanish, which is the case for ANY axis-aligned collinear triple, an ordinary configuration on a structured mesh, not just for coincident points. The caller declined the flip either way so no mesh was ever corrupted, but a predicate whose entire contract is "report a sign only when the sign is justified" was reporting one it could not justify. It now returns UNCERTAIN, with a regression test covering coincident, x-collinear and y-collinear input as well as the unambiguous cases. pin_bands had no parallel test, which Charter section 11 does not allow. It works, and the new test asserts the properties that make it safe rather than just that it runs: the pinned set is partition-independent (compared by COORDINATE, since a shared vertex is held by every rank on the seam and a count would double-count it and mask the defect); pinned vertices do not move even when they are star-forest LEAVES owned by another rank, which is the case a rank-local pin would get wrong; and the domain boundary stays pinned. Verified np=2 and np=3. Underworld development team with AI support from Claude Code
Adversarial reviewWe reviewed this branch against itself before marking it ready. Two findings were Fixed in this branch1. 2. Accepted, with reasons
Under a frozen seam the 99th-percentile angle recovers but the absolute maximum
Performance.
Where the numbers came from, and what we got wrong getting themSeveral figures in the description are revisions of earlier, wrong ones. Recording
The controls that caught these are in the tests: parallel confluence, volume |
Adversarial review — live probes at head 4b041f7 (isolated worktree, own build)Three findings block merge:
Minor: (4) Attacks that failed: serial suites 18/18 in 25 s; ptest_0843 np=2 and np=3 pass against the serial 412-cell/7-pass reference; ptest_0844 np=2/np=3 4 passed each (shared-point cones untouched, chart/area invariant, solve on repaired mesh converges); ptest_0845 np=2 3 passed; no bare exception swallows (the four narrow SF-getGraph guards carry rationale comments); all Holding both this and #489 until findings 1–3 are fixed; 1 is a hang on the exact workflow the #489 skills then recommend. |
Three capabilities for locally-refined (adapt-on-top) meshes, plus the design note
recording what was measured. All parallel, all tested at np=2/3/4.
What lands
mesh.adapt(engine="edge_split")— longest-edge refinement with no conformingclosure. Splitting an edge divides every incident cell at the same new vertex, so
there is no hanging node to repair and refinement cannot escape the marked region:
the band hugs the feature instead of a bounded halo around it. No new topology code
— it drives the compiled
uwnvb_bisecttransform that NVB already uses, so itinherits star-forest propagation, co-partitioning, labels and coordinates.
Bit-confluent: identical mesh at np=1/2/3/4, verified to 56k cells.
Marking is on the cell diameter, not
(d! V)^(1/d). The volume proxy reportedthe target met while the mesh was 3.2x coarser across the feature.
mesh.adapt(..., repair=True)— a reconnection (Lawson flip) pass after eachgeneration. 2-D only; raises rather than silently doing nothing elsewhere.
The acceptance criterion is not Delaunay, although these are Lawson flips.
Delaunay maximises the minimum angle and says nothing about the maximum, while
the P1 interpolation bound depends on the maximum (Babuska-Aziz). Measured:
flipping a gmsh-generated mesh towards Delaunay raised the 99th-percentile
maximum angle from 126.8 to 129.3 degrees, because gmsh optimises element shape
rather than the empty-circle property. Since every UW3 mesh starts from gmsh, a
pass that can degrade one is unusable. Gating on the angle makes it monotone.
Worth it on a poor base — 99th-pct max angle 156 -> 115 degrees on an
aspect-ratio-4 grid, slivers below q=0.1 3.84% -> 0.00% on a non-Delaunay one. On
a clean gmsh base it moves 124.7 -> 120.5 and the error not at all.
Opt-in, because it gives up the one property
edge_splithas: which cavities maybe flipped depends on where the partitioner cut, so the repaired mesh is not
partition-independent. Conformity, orientation, volume, labels and the
star-forest stay exact at every rank count.
mesh.relax(pin_bands=...)— hold an interface while relaxing everything else.Relaxation and interface-tracking refinement fight: the MMPDE mover optimises shape
against an equilateral reference and knows nothing about where the material
changes, so it slides the small cells refinement placed on an interface off it.
Measured on a step-edged fault, manufactured stress across the interface +77%.
Pinning the band leaves that unchanged to five decimal places while the mover still
reshapes the rest of the domain.
Implementation notes worth a reviewer's attention
A flip is not expressible as a
DMPlexTransform— a child's cone may onlyreference its own parent's closure, and a flip's output cells use the other
parent's apex. So the DM is rebuilt. It is rebuilt on the same point chart: a
2-D flip adds and removes no points, so preserving the numbering lets the point
star-forest transfer verbatim, labels transfer by point id and coordinates transfer
unchanged. That removes the whole reconstruct-the-SF-by-matching-seam-coordinates
stage. Surgery on the source DM is impossible:
DMPlexSymmetrizerefuses to run ona plex that already has supports and nothing outside
DMDestroyfrees them.Exact predicates are a static filter that declines when it cannot resolve a
sign, not adaptive precision. Declining is always safe because a flip is an
optimisation, never a requirement.
Tests
22 serial, 9 parallel at np=2/3/4. 149 adapt/relax/smoothing regression tests pass;
style gates clean.
The load-bearing ones are the controls: parallel confluence (three of the four
edge_splitdefects during development were invisible serially), volumeconservation rather than orientation (checking that new cells are positively
oriented is worthless when they were built anticlockwise by construction), and the
maximum-angle assertion, which is what caught the Delaunay criterion.
Companion
Skills documentation is split out to #489 — it documents this API, so it should
merge after this.
Underworld development team with AI support from Claude Code