Adapt-on-top: closure-free edge_split engine, reconnection repair, and interface-pinned relaxation - #488
Open
lmoresi wants to merge 23 commits into
Open
Adapt-on-top: closure-free edge_split engine, reconnection repair, and interface-pinned relaxation#488lmoresi wants to merge 23 commits into
lmoresi wants to merge 23 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
Member
Author
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 |
…dition mesh.add_conforming_surface(points, name) splits every edge the surface crosses at the crossing point, so the surface becomes a chain of element edges. No element straddles it, a material property can be assigned per CELL and be exactly right, and the surface becomes a named boundary that a solver can apply conditions on. The point of adding it on top of an existing mesh, rather than building it into the mesh generator, is that its position need not be known when the mesh is made: the base mesh and its multigrid hierarchy stay fixed while the surface moves, which is what an outer optimisation over its position needs. Why the straddling matters: a linear element forms stress from the interpolated viscosity times the interpolated strain rate, so it carries mean(eta)*mean(edot) where the honest average is mean(eta*edot). The difference is -2 Cov(eta, edot) per cell, zero for any element wholly inside or outside the zone and positive only across the transition. Refinement shrinks the straddling band but never empties it. Measured on a step viscosity 1 -> 1e4: the leak is 285 on an uncut mesh and EXACTLY zero on a cut one with a cell-wise viscosity. A continuous P1 viscosity still leaks on a cut mesh (227 against 240 uncut) because the nodes ON the surface are shared by both sides -- the cut is what makes a per-cell assignment correct, not smooth. SolCx, the acceptance test, eta 1 -> 1e6 on an irregular mesh at matched cell count: a regular mesh that already conforms takes 14.1 s for a relative L2 error of 2.5e-05; the cut irregular mesh takes 17.9 s for 1.3e-05; the same mesh uncut takes 271.8 s for 4.5e-02. So the cut costs about 27 % over the ideal and is 15x faster and 3600x more accurate than leaving the mesh unaligned. No new C. The compiled uwnvb_bisect transform already inserts a vertex per marked edge; only the coordinate needed overriding, and the topology follows. Implementation notes worth keeping: * PASSES OF PAIRWISE-INDEPENDENT EDGES, not one pass. The transform can split two edges of a triangle at once and emit the joining segment -- the whole cut in a single pass. That is correct in serial and WRONG IN PARALLEL: the double-split path leaves the child point star-forest inconsistent and wrapping the result as a Mesh dies in PetscSectionCreateGlobalSection at np>=3. Its own source calls those tables "a safety net"; nothing had exercised them across a partition. Independent single splits still build the cut, because the second pass joins its new vertex to the opposite vertex of the cell, which is the first pass's new vertex. * SNAP OR CUT, measured ALONG THE EDGE. A crossing landing near a vertex leaves a sliver -- in the worst case an area of 1e-24 and a zero angle. A crossing within snap_frac of an edge's end moves that vertex onto the surface instead. The along-edge measure is the short side of the sliver that would otherwise be created and carries no length scale. GAMG on a Poisson solve, which is sensitive to element shape where the geometric hierarchy deliberately is not: uncut 20 iterations, snap_frac 0.00 32, 0.05 28, 0.10 23, 0.20 21. Hence the 0.10 default. A Lawson flip pass helps less (32 -> 29, 28 -> 25), so snapping is the better lever and repair is a touch-up rather than a requirement. * EVERY rank-local decision is reconciled. Four collective bugs, all the same shape -- a rank-local branch around a collective -- and all invisible at np=1 and np=2, because a two-way split happens to give every rank a piece of the surface. np=3 exposed all four: the tip / triple-crossing / multiply-crossed validations, the "nothing to cut" guard, the snap-set reconcile itself, and the substantive one -- the snap decision is read off an EDGE, so a rank holding one side of a shared vertex could decide differently from its neighbour, leaving the ranks disagreeing about which edges were crossed and the split loop never emptying. * cut_hierarchy is OFF by default. It is tempting to argue a surface-free coarse level "solves a different problem", but custom-P sets pc_mg_galerkin=both, so every coarse operator is PtAP from the FINE operator and inherits the contrast whatever the coarse mesh looks like. What a coarse cut would buy is a coarse SPACE able to represent the kink; measured on SolCx at contrasts of 1e2 and 1e6, cutting the coarse levels moved the error in the fifth significant figure and the solve time not at all. Scope: two dimensions, and surfaces crossing the mesh from boundary to boundary. A surface ending inside the mesh (a fault tip) is refused rather than silently mis-meshed, as is a triangle crossed three times. Tests: 18 serial, 8 parallel passing at np=2/3/4. The parallel file asserts the mesh by sorted owned-vertex COORDINATES and a hash rather than counts (derived counters lie in parallel), and solves a Dirichlet problem on the surface, matching the serial domain integral to 4e-17. Both solves are driven to a tight tolerance so that can be asserted strictly: at default tolerance the two differ by 1.5e-8, which is two iterative solves converging within their own rtol rather than a partition effect, and a loosened bound would have hidden the question. Underworld development team with AI support from Claude Code
A refinement engine takes as many passes as it needs to reach the size the
metric asks for: independence caps how many edges one pass may split, and a
conforming closure cascades. So a pass is how the engine REACHES a size, while
a multigrid level is a COARSENING RATIO. adapt() conflated them by recording
every pass as a level, and nothing connected the two numbers:
edge_split n_pass = 8*dim*max_levels is only a CAP; the loop runs to metric
satisfaction, so max_levels 1/2/3 returned byte-identical meshes
and 10 passes became 10 levels;
nvb n_gen = dim*max_levels, and a bisection is a 2^(1/dim) step in h,
so `dim` generations make ONE h-halving -- you got dim times as
many levels as isotropic-equivalent ones.
Both then degenerate: once the metric is nearly met the passes coarsen nothing
(measured ratios 1.06, 1.02, 1.007) and each such level still costs a full
Galerkin RAP and smoother sweep. That hierarchy stopped SolCx converging at all.
adapt() now takes mg_coarsening_ratio (default 2.0, applied identically by both
engines) and keeps one level per that much coarsening in h.
THE MEASURE IS RESOLUTION, NOT ELEMENT COUNT. Under adapt-on-top the mesh only
grows where the feature is, so a genuine halving of h shows up as a global cell
ratio near 1: on a thin band NVB grew the mesh 1.06-1.11x per generation while
the in-band h went 0.125 -> 0.0626 -> 0.0313 -> 0.0157. A count-based rule keeps
nothing and collapses the hierarchy; the whole-mesh median h is flat and equally
useless. The selector uses a low percentile of cell diameter, reduced with MIN
across ranks, and replaces rather than appends when the level below the finest is
within the ratio -- appending reintroduces the near-duplicate pair it exists to
remove.
Measured on SolCx with the interface CONFORMING at the finest level, so the
discretisation pathology of an unaligned jump does not swamp the comparison
(uncut, SolCx at 1e6 does not finish at all):
engine hierarchy levels vel its seconds
nvb per-pass 7 4 19.67
nvb doubling 5 5 6.95
edge_split per-pass 11 5 161.04
edge_split doubling 6 6 22.16
2.3x to 7.3x faster for +0 to +1 iterations, at errors identical to four
significant figures, and contrast-independent (iterations barely move from 1e4
to 1e6). The extra levels were overhead. A ratio sweep at np=1/2/4 shows the
ranking is stable and that cost keeps falling to ratio 3 before saturating; the
default stays at the conservative 2.0 and the knob is exposed.
Prolongations are COMPOSED across the passes a level spans, so the recorded
transfer stays exact instead of falling back to the geometric builder. Composed
in numpy: each row of a bisection prolongation holds one or two entries, so
expanding the fine map through the coarse rows and summing duplicates is the
whole operation. Validated against a dense oracle on 200 random cases, exact and
a partition of unity.
Tests updated to the new contract:
* test_0753 asserted two SINGLE-GENERATION properties -- every fine vertex lies
on a coarse edge, and at most 2 nonzeros per row. Neither survives composition
and neither should: a composed span can place a vertex strictly INSIDE a coarse
cell, where it depends on that cell's dim+1 vertices. The reference is now
barycentric-in-cell, which covers every fine vertex instead of the ~64 % that
lie on an edge, so the test checks MORE than it did; the sparsity bound becomes
dim+1.
* test_0836 / test_0840 tied the level count to the generation count. They now
assert the property that defines the contract: no level is a near-duplicate of
its neighbour, and interior adapted steps reach the requested ratio. The step
INTO the finest level is exempt -- the finest is the child and is mandatory, so
when the whole adapt is less than one doubling its single step is whatever the
metric asked for (1.74 measured in 3-D).
FOUND ON THE WAY, NOT FIXED: nvb.nested_prolongation is wrong in 3-D for vertices
a closure cascade places strictly inside a coarse tet -- worst |P.u - P1(x)| =
1.19, measured PER GENERATION with no composition involved, against 1.9e-15 in
2-D. It was masked because the old reference was edge-based and skipped exactly
those vertices. Marked with TODO(BUG) at the source and xfailed (strict) in
test_0753; it predates this change and is not caused by it.
Underworld development team with AI support from Claude Code
…face count `_resolve_snapping` initialised its on-surface set to all-False and only added vertices it decided to SNAP. A vertex ALREADY lying on the surface was therefore invisible to it -- the edges radiating from such a vertex have signed distance exactly zero and register no strict sign change, so nothing ever proposes them. That is fine for a surface crossing open mesh, and wrong for a fault NETWORK. A junction (or a tip) is placed by pulling a mesh vertex onto it, so it lies exactly on every branch that meets there. The validation then read the cell beyond it as "entered but not left" and refused a legal branch. Seeding the set with vertices already on the surface fixes it. Measured, on a 1/20 box with the junction pulled onto a vertex: Y three arms from one junction 3 branches, zone 116 cells, 0 inverted T one fault abutting another 2 branches, zone 114 cells, 0 inverted X two faults crossing 2 branches, zone 166 cells, 0 inverted all branches labelled chains of mesh edges, in every case. Y previously failed; T and X already worked, which is what made the cause specific -- both of those have a branch passing THROUGH the junction, so an ordinary crossing marked the vertex as a side effect. This is the "crossings computed twice from different sources" smell already recorded in the design review, producing a false refusal. The pass loop derives its on-surface set correctly (`distance < 1e-12 * scale`); only the validation path did not. The single-source-of-truth refactor should absorb this. Why networks matter here: a one-element fault zone taken as the cells in the SUPPORT of the labelled facets makes a network's zone the UNION of its branch zones -- no geometry to reconcile where branches meet, in any dimension. The alternative (offset surfaces plus end caps) has to mesh T- and X-junctions conformally, and for a one-element-wide fault that is self-contradictory: the cap has extent equal to the thickness, so resolving it needs h << h. Maintainer ruling 2026-08-02, recorded because it scopes the work: intersecting faults are transient -- if they slip they change the geometry -- so an approximation to the fault volume is fine, and junction geometry need not be resolved exactly. The union-of-cells zone bulges where branches meet, since the fan around the shared vertex is picked up by each branch. That is an accepted characteristic, not a defect to engineer away. Underworld development team with AI support from Claude Code
Member
Author
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. |
… no cut below the child Adversarial review of this branch found six correctness defects and a test suite several of whose tests passed with the feature removed. This is sections A and B of that triage, plus a maintainer ruling that removes a whole path. THE SURFACE EXISTS ON THE FINEST LEVEL ONLY. `cut_hierarchy=` is gone, along with `_cut_coarse_levels`. Cutting the coarse multigrid levels produced a hierarchy of cut copies of the base levels, which defeats the point of the stack-on formulation: the surface's position is a design variable in an outer optimisation, so the base and the hierarchy resting on it have to stay fixed while the surface moves. It bought nothing either — custom-P sets pc_mg_galerkin=both, so every coarse operator is PtAP from the FINE operator and carries the contrast whatever the coarse mesh looks like (SolCx at 1e2 and 1e6: fifth significant figure, no time difference). It was also the path with zero tests and the one where two of the defects below bite. EVERY REFUSAL IS NOW GLOBAL. A rank-local raise aborts one rank while its peers walk into the next collective and block there, so the error becomes a hang. Nine defects of this shape have now been found in this module, and the parallel suite could not see any of them because it only ever took the happy path. Audited as a class rather than fixing the five named: * the cell-inversion raise, the `_child_vertex_of` raise (which sat inside a rank-local "did this rank split anything?" guard as well), and the guard around the coordinate write are all gone or reduced first; * `_global_extent` replaces five rank-local `np.ptp(...).max()` calls. Those raised outright on a rank owning no vertices, and one of them fed the crossing tolerance — so the module's central invariant, that every rank computes the same crossing from the coordinates alone, was false (measured spread 0.58-0.67 against 1.0 serial); * every number in `info` is reduced, counted over owned points, so the documented identity between them can hold at np>1. `n_snapped` becomes `n_on_surface`, which is what it has counted since junctions were seeded into it. Measured negative control: restoring the rank-local form of the inversion test HANGS at np=3 on exactly that case while the three refusals before it pass. THE STRESS LEAK IS ASSERTED. It is the claim every docstring and commit message on this branch rests on and it was tested nowhere. On a 1/16 box at contrast 1e4: uncut leaks 285.4 with a cell-wise viscosity, cut leaks exactly 0.0, and a continuous P1 viscosity leaks 298.7 even when cut — so the feature is "cut AND assign per cell", not "cut". Stubbing add_conforming_surface to return the mesh unchanged fails it. Tests that passed with the feature stubbed out, and now do not: * the parallel snap test selected vertices within 1e-6 of the surface and asserted the worst was under 1e-12. On the uncut base that set is EMPTY (nearest vertex 5.5e-3), so it held with the feature removed. Now the count and identity of on-surface vertices against serial; * `no_inverted_cells` was true by construction twice over — cut_along_lines already raises on the same areas, and min_angles is arccos of a clipped value. Now the documented angle table (1.60/3.88/6.56/13.93 deg); * the coarsening-ratio knob passed with the ratio hard-coded ([3,3,3] is still non-increasing). Now strict decrease; * `_assert_coarsening_ladder` re-derived the implementation's own level selector and passed ratio 2.0 at 1.817 against 1.800. Now an INDEPENDENT estimator (mean edge length in the refined band), shared between the 2-D and 3-D suites instead of duplicated verbatim, asserting the adapted SPAN rather than a per-step number the engine never promised. That same step measures 1.401 independently; * the parent-cell map was discarded unconditionally after subsampling, which tautologised the repair test. Kept per level when the level is one generation. test_0753 (tier_a): the barycentric reference REPLACED an edge-membership one on the grounds that it covered every fine vertex rather than 64 %. That 64 % is the 3-D case, which is xfailed; in 2-D nothing composes and the old reference already covered 100 %, so it was a loosening. Both references are kept now — edge membership catches a PHANTOM parent edge, which is the 3-D defect and which barycentric position and linear-field reproduction are both blind to. Added a 2-D case that genuinely composes, so the docstring's claim is exercised somewhere that runs. Sparsity is bounded PER ROW, not on the mean, since dim+1 IS point-location density. The 3-D defect is asserted positively instead of by strict xfail on one row in 2336. Smaller: _boundaries_with could land a surface on Null_Boundary(666); _cut_coarse_levels caught only ValueError when two of three failures are RuntimeError; the cut child is marked as not having coincident DOFs, so _refine_restrict interpolates rather than injecting from a displaced node; uw.pprint(0, ...) printed a literal 0; a malformed RST table would have broken the Sphinx build. A2 (coarse levels carry no boundary, so an essential BC on the surface is unsound) is DEFERRED. The docstring no longer claims otherwise. Underworld development team with AI support from Claude Code
… parallel The fault is a one-element-wide zone defined at the FINEST level of the adapt-on-top, and the zone is the cells in the SUPPORT of the labelled facets — not a geometrically bounded region. `mesh.cells_supporting(name)` is that zone. It needs no end cap, no edge band and no rim; it terminates automatically where the chain of facets ends, it says nothing about dimension, and the zone of a network is the union of its branches' zones with no geometry to reconcile where they meet. Bounding it geometrically is self-contradictory for a one-element fault anyway: the cap has extent equal to the thickness, so resolving it would need h much smaller than h. Measured, and asserted: * the zone is EXACTLY 2 x facets at every resolution tried. A cell carrying two labelled edges would have been cut in two, so no cell is double-counted and every facet contributes both neighbours — one element each side, by construction; * thickness tracks the LOCAL h: 0.189 / 0.183 / 0.184 across a 4x uniform refinement, and 0.195 / 0.202 / 0.198 under the adapt metric. So width is a REFINEMENT parameter — the surface lives at the finest level and the metric decides how wide one element is, controlled locally and at bounded cost; * adapt THEN cut composes, and the child keeps its multigrid tail. That is the order the design needs. (adapt refusing to chain ON a cut child is the other direction and is not what the fault requires.) The max centroid distance will NOT do as the thickness statistic: it is one outlier cell and it came out bit-identical at two different adapt resolutions, reporting no scaling where the mean shows it cleanly. add_conforming_surface takes a Surface, not (points, name). It is what fault_metric, fault_metric_tensor and refinement_metric_function already take, so one object drives the refinement metric AND the cut instead of being unpacked and its name re-stated, and it carries signed_distance and director for the weak-plane model afterwards. Control points are read in MODEL space via the machinery's own _fault_collect_polylines — surface.control_points is the dimensionalised gateway and would be the wrong space under an active units system. pull_vertex_onto() is promoted out of the test file into the library, because a TIP and a JUNCTION are the same problem — a distinguished point that must coincide with a mesh vertex, after which every branch meeting there arrives at the already-legal "one crossed edge, one on-surface corner" case. It is now COLLECTIVE: the test helper took a rank-local nearest vertex, which moves a DIFFERENT vertex on each rank so the branches meet at different places either side of a seam. Reduced as (distance, x, y) so the tie-break rides along in the same reduction, and the move is applied by POSITION so a ghost copy lands in the same place without a star-forest exchange. Fault NETWORKS now run in parallel — Y, T and X at np=2/3/4, previously untested. Negative control: restoring the rank-local vertex choice fails all three at np=3. The fault zone is checked across the partition too, by owned count AND by a hash of the sorted zone centroids, since a count alone can agree between two different sets of cells. Also asserted, because the docstring tells users to rely on it: degree-0 DOF order IS plex cell order, so cells_supporting can be assigned straight into a P0 viscosity. Were that untrue the contrast would land on the wrong cells and every downstream result would be quietly wrong while looking plausible. Underworld development team with AI support from Claude Code
`vis.labelled_facets_to_pv_mesh(mesh, name)` returns the facets carrying a boundary label as a PolyData of their own — lines in 2-D, triangles in 3-D, since a labelled facet's closure gives its vertices whatever the dimension. An embedded surface drawn WITH the mesh is a few lines among thousands in 2-D and completely occluded in 3-D, so it has to be separable to be looked at. It saves to `.vtp` for interactive viewing, which is how the 3-D version will have to be inspected. `docs/developer/subsystems/conforming-surfaces-and-fault-zones.md` is the design note the branch was missing entirely: why straddling elements are a representation problem rather than a resolution one, the leak table, why the zone is the facet support and not a bounded region, the thickness-tracks-h measurements, the snap_frac trade, tips and junctions, and the limitations. The GAMG table moves out of the line_cut docstring into it, leaving the sentence that justifies the default — which also removes the malformed RST that would have broken the Sphinx build. `line_cut` is exported from `utilities/__init__` alongside `edge_split` and `reconnect`, so it is not deep-import-only and its cross-references resolve. Underworld development team with AI support from Claude Code
… reach Two additions to `cut_along_lines`, both driven by the same measured fact: the cut's slivers are made by the SPLITS, so anything that replaces a split with a vertex move helps and anything that turns a move back into a split hurts. `snap_quality` — a triangle-quality floor on snapping. A cell thinner than the tolerance band has every corner pulled onto the line from both sides and is flattened; measured on a graded mesh, every collapsed cell at snap_frac 0.4 had all three corners snapped, and the cut was refused outright. A proposed move that would take an incident cell below the floor is now vetoed and that crossing is split instead. The floor is absolute and monotone (never below it, never worse if already below), because a floor expressed as a fraction of the CURRENT quality compounds when the routine is applied repeatedly — 0.5 over six rounds licenses 0.5**6, and the worst angle duly fell 15.4 -> 2.3 degrees with every individual round looking well behaved. The guard must test QUALITY, not inversion: a flattened cell lands at ~1e-16 of either sign, so half survive an inversion test, the worst angle still reaches zero, and the returned mesh looks fine while the cut chain has silently broken. Guarding on inversion alone was measured doing exactly that. The default is deliberately LOW (0.15). The guard protects the snapped mesh, which is not the mesh that comes back. Raising the floor from 0.15 to 0.55 held the snapped mesh's worst angle up (15.6 -> 24.5 degrees) while driving the CUT's down (10.9 -> 0.16) and the split count up (139 -> 359). It is a backstop against flattening, not a quality target. `None` removes it entirely, restoring the pre-guard behaviour and its refusal. `snap_dist` — snap any vertex within that multiple of its own local h of the line, whatever the crossings on its edges look like. `snap_frac` is measured ALONG an edge and is blind to a vertex sitting close to the line while every edge meeting it is crossed near its midpoint. That vertex becomes the apex of a cell with one edge on the cut, which is the characteristic sliver: of the sixty cells below 15 degrees in a box-fault cut, ALL sixty had two corners on a cut and ALL sixty were elongated along it, apex about 0.45 W away. Five separate knobs (snap tolerance, quality floor, staged refinement, metric ramp slope, metric core width) each returned a worst angle of 10.80 degrees and ~59 poor cells, to the digit, because none of them can reach that configuration. `snap_dist` 0.30 halves the population (60 -> 35 on the box, 26 -> 13 on a single cut). It is OFF by default: it also makes the worst single cell worse (10.8 -> 1.4 degrees), because the splits it leaves behind sit in harder places and nothing guards the splits. That gap is the next piece of work, not something to enable by default ahead of it. Also: `add_conforming_surface` forwards both, and its `snap_frac` docstring now records that 0.10 is not the right value on a graded mesh (0.30 took the worst angle from 4.96 to 10.81 degrees and cells below 15 from 231 to 31) without changing a default chosen on a uniform one. Tests: two serial tests — the guard turns the flattening refusal into a valid cut, and `snap_dist` finds vertices the along-edge test does not. The parallel collective-refusal case now passes `snap_quality=None` so the refusal path it exists to protect is still reachable. 36 serial, 16 parallel at np=2/3/4. Underworld development team with AI support from Claude Code
… label `_cell_regions` builds a per-cell signature from every non-topology label and locks any edge whose two cells disagree, on the reasoning that such an edge is a material interface even when unlabelled. `uwnvb_refedge` is not a material label: it records which of a triangle's edges is its refinement edge, and it takes values 0/1/2 across any NVB-adapted mesh. Measured on an adapted fault mesh, that read as three regions of 2230/2184/134 cells, and every edge between them was locked. The effect was not marginal. Of the edges around a sub-15-degree cell in a cut mesh, 113 were declined as a "region interface" against 54 genuinely locked on the fault. Excluding the label takes the pass from 101 flips to 483, and cells below 15 degrees from 60 -> 18 rather than 60 -> 59; cells below 25 degrees go 420 -> 239 and the 1st-percentile angle 14.4 -> 18.1 degrees. This is the same trap `_labelled_points` already documents for `Elements`, one level along. That fix — ignore a label carried by CELLS — cured `Elements` because `Elements` is uniform, so it never reaches `_cell_regions`' final "are all signatures equal" test. A bookkeeping label that VARIES over cells does. The gate itself was never the problem, and is unchanged: of the edges around a sliver, the 44 with a minimum-angle gain are exactly the 44 with a maximum-angle gain, so a Delaunay-style gate would have flipped the same set. Only the lock differed. The fault is untouched, as it must be — flips are locked on labelled edges. Cut and cut+flip agree to the digit on both flanks: 317 and 318 facets, every chain vertex within 1.3e-16 of the line, zero straddling cells, zero inverted, and the minimum cell area rises 4.7e-7 -> 6.7e-7. Test: a regression with its own `adapt`-built fixture, since the file's shared `_refined_dm` goes through `bisect_longest_edges` and never carries the slot label — which is why the defect survived this suite. It asserts the label is present AND that it takes more than one value on cells, so a fixture that could not expose the defect fails loudly rather than passing vacuously. Underworld development team with AI support from Claude Code
… applies `add_conforming_surface` appended the mesh it was cutting to the child's coarse tail unconditionally, on the stated reasoning that "adding a surface refines this mesh, so this mesh plus everything below it is a valid coarse tail". The premise is wrong. A cut re-represents the same grid with the surface conformed; it adds no resolution. Measured on a box fault, the two cuts produced two levels that coarsened h by 1.11x and 1.17x on the 5th-percentile measure, against a threshold of 1.8 — each one a full Galerkin RAP and a smoother sweep for no correction. `_subsample_mg_levels` already decides exactly this question for an engine pass, including the "replace the level below rather than append to it" case, and it is the committed answer to it. So the cut path now calls it, handing it the pair (self, child) measured against the level beneath them, rather than carrying a second rule that could drift from the first. `mg_coarsening_ratio` is exposed to match `adapt`. Box fault: 9 levels -> 7, and the top transition goes from 1.06x to 1.96x in mean h. One cut: 8 -> 7. The hierarchy is now the same depth as the adapted mesh it was cut from, which is the point — cutting is not refining. Two things fall out, both measured on the same shear solve: * the barycentric transfer stops failing. Transfer 7->8 ran BETWEEN the two near-duplicate cut levels, and it was there that the builder ran out of coarse DOFs with a fine image and fell back to the dense-RBF one (#424) — dense Galerkin coarse operators, and a measured 94s/0.6GB turning into >21min/12GB when the mesh was also relaxed. The fallback no longer fires, in this solve or anywhere in the two test suites. * the solve is 1.87x faster for the same answer: 93.7s -> 50.0s, strain-rate ratio 133 either way and the fault strain rate 58.04 -> 58.03. The reported V-cycle count went 8 -> 15, which is NOT a regression and should not be read as one: it counts the last inner solve only, and the hierarchy under it changed. Time the solve. Test: `test_the_surface_exists_on_the_finest_level_only` asserted the tail keeps its length and that its finest level is the base finest. Both described the old contract. Its substance — coarse levels carry no surface label, the base is not mutated, the tail is built from the base's own uncut level objects — is unchanged and still asserted; the count and the identity of the finest level now say that the cut REPLACED the base finest. 45 serial, 16 parallel at np=2/3/4. Underworld development team with AI support from Claude Code
A conforming cut has only two primitives — snap a vertex onto the surface, or split an edge it crosses — and every sliver it leaves follows from that. A crossing falling near a vertex must either drag the vertex to it or carve a thin cell beside it, and tightening the snap tolerance only trades one for the other. Delete is the missing third: it dissolves the case, and it is the only one of the three that removes work rather than adding it. On a box fault cut into an adapted mesh, counting cells under 15 degrees: the cut leaves 60, flipping takes that to 18, and deleting afterwards to 4 while removing 242 cells. The order is not symmetric — deleting first leaves the count at 60, because a cavity, once ear-clipped, no longer presents the quad the flip pass was looking for. The pair then converges: a second round of each finds nothing. The fault itself is bit-identical through both passes, at every rank count. The acceptance test needs both shape measures, unlike the flip pass. Gating on the largest angle alone — correct for flipping, since the P1 interpolation bound depends on it — let the minimum angle fall from 10.80 to 10.23 degrees and RAISED the sliver count from 60 to 61, because a needle has one tiny angle and two close to 90 and never registers as obtuse. Hence gate="both". Parallel is one exchange, not a redistribution. Deletion compacts the point chart, so unlike a flip it cannot hand the star-forest across verbatim: every point after a deleted one shifts, and each leaf's remote index is a number only its owner holds. rebuild_without_vertices renumbers locally and broadcasts the new numbering root-to-leaf once. Freezing the seam is what keeps the leaf set itself unchanged, so the forest is renumbered and never rebuilt; it costs 113-115 deletions against 121 serial at np=2..4. Also fixes the third instance of one labelling trap. Null_Boundary marks every vertex of every UW3 mesh with the reserved value 666, and UW_Boundaries re-packs every per-boundary stratum, sentinel included, into one stacked label — so reading labelled POINTS as interfaces flags the entire vertex stratum. That costs the flip pass nothing, since it asks only about edges, and it refused 1114 of 1114 candidates the first time the removal pass met a cut mesh. _labelled_points is now _interface_edges and reads edges only, which is the right reading anyway: in 2-D an interface is a curve. It is also the only reading that protects a fault, since cut_along_lines labels the cut's edges and not its vertices. Underworld development team with AI support from Claude Code
The cut can only snap a vertex onto the surface or split an edge it crosses, so a crossing landing near a vertex either drags the vertex to it or carves a thin cell beside it, and tightening snap_frac only trades one for the other. repair=True runs the two operations the cut does not have: flip, then delete. On a box fault, cells under 15 degrees go 60 -> 4 while 242 cells are removed. The surface's own facet count is unchanged, since both passes refuse to act on a labelled edge. Deletion is offered only the vertices within repair_reach * h of the surface. It removes degrees of freedom, and the cut is what justifies removing these particular ones; a pass turned loose on the whole mesh would coarsen it wherever the shape happened to be poor. Flipping is offered everything, because it conserves the point set. Off by default, like adapt(repair=...), because the cut alone gives the same mesh at any rank count and repair gives that up. Also fixes needle blindness in the flip pass. It gated only on the pair's largest angle — right as an OBJECTIVE, since the P1 interpolation bound depends on it and Delaunay is the wrong criterion here — but nothing stopped it buying that gain by making a thin cell, whose largest angle is unremarkable and so never registers. Measured: on a cut graded mesh, flipping alone took the smallest angle in the mesh DOWN. The objective is unchanged; this adds a floor under the other end, which is the same correction the deletion gate already carries. Found by composing the two passes, which is the only place it shows. Underworld development team with AI support from Claude Code
… them meshVariable_to_pv_mesh_object triangulates a variable's nodal points with delaunay_2d. That exists so higher-order fields can be plotted at all -- the base mesh does not carry their DOFs. For a CONTINUOUS P1 field it is the wrong thing to do: the DOFs are the vertices, so the triangulation is already in the DM. And it is lossy, not merely redundant. delaunay_2d takes one alpha for the whole domain and discards triangles whose circumradius exceeds it, so on a graded mesh it deletes the COARSE cells. Measured on a fault mesh graded 8:1, 361 of 11610 cells were dropped, and they render as blank holes in the middle of the field -- which reads as missing data and was in fact mistaken for one. meshVariable_to_native_pv_mesh returns the DM's own cells, renumbered so that point i is the variable's DOF i, and mesh_to_pv_mesh already did the hard half of that. The renumbering is the load-bearing detail: the documented usage attaches values by DOF index, so handing back the right cells in the DM's vertex order would draw a plausible field with the values shuffled. The permutation is found by coordinate match and asserted, not assumed, and the helper returns None -- falling back to Delaunay -- whenever the DOFs are not one-per-vertex. Automatic, so every existing call site is fixed without change. Passing an explicit alpha keeps the old path. The test fixture is deliberately GRADED, with a control asserting that the Delaunay route really does lose cells on it: on a uniform mesh the two agree and the regression is invisible. Underworld development team with AI support from Claude Code
plot_mesh_hierarchy draws a mesh, its multigrid tail and its faults in one figure -- one colour per level, coarsest palest and thickest, fault zones filled in a contrasting red. It answers the three questions that come up every time a mesh is built this way: did the hierarchy come out with the levels expected, is the refinement where the fault is, and did the fault survive the repair passes. Written for 3-D rather than adapted to it later. Nothing reads the dimension except the defaults: in 3-D the wireframes come from each level's SURFACE, because extracting every interior edge of a tetrahedral hierarchy is an unreadable haze, and `clip` cuts the model open so the interior levels and the fault can be seen at all. The fault selector is cells_supporting, which is already dimension-general -- a fault zone is the support of its labelled facets whether those are segments or triangles. The colour taper is load-bearing, not decoration: drawn at one width the finest level's edges cover every level beneath it and the hierarchy cannot be read at all. Tests assert what was DRAWN -- an actor per level and one per fault -- since that is how this can silently mislead. A hierarchy missing a level reads as a shallower mesh; a fault that contributed no actor reads as a mesh with no fault in it. Both would look like perfectly good figures. Underworld development team with AI support from Claude Code
plot_mesh_hierarchy filled cells_supporting(name) in red. That is the fault ZONE -- every cell with a labelled facet, which is one element on EACH side -- so a one-element-wide fault came out two or three elements thick and looked like something the mesh does not contain. The facets are the fault as the mesh represents it, and labelled_facets_to_pv_mesh already returns them, dimension-general: segments in 2-D, triangles in 3-D. That is now the default. fault_style="cells" keeps the zone fill for the question it does answer -- which cells carry the weak viscosity -- and an unrecognised style is refused rather than silently drawing nothing. The test now asserts the two sets DIFFER, so the default cannot quietly revert to the fat one and still pass. It counts n_lines + n_faces_strict, not n_cells: `pv.PolyData(points)` gives every point its own vertex cell, so n_cells is n_points plus the lines and reads as a wildly wrong facet count -- 127 for a 63-segment chain. Underworld development team with AI support from Claude Code
…angles for faults Shape carries the distinction as well as colour, so the figure survives being printed in grey and does not ask anyone to tell four blues apart. In 3-D the same three roles become sphere, cube and cone, and that branch is exercised by the tests rather than left until there is a 3-D fault to look at. Sizing them took two goes and both failures are worth recording, because they are the same mistake at different scales. Scaling each level's glyphs by ITS OWN cell size seemed natural -- coarse level, coarse marks. Zoomed in on the fault it is a disaster: the coarse level's marks are drawn at the coarse spacing and blanket the fine mesh completely, which is precisely the view the figure exists for. Sizing them all by the finest level's MEAN cell size then failed for the reason a graded mesh always breaks a mean: the fault meshes here average h = 0.017 while h at the fault is 0.002, so every mark came out several times larger than the cell it stood on. The 5th percentile is what is wanted -- the size of the cells that actually need marking. This is the same trap as judging a multigrid level by its mean h. Node actors are unlabelled: a legend line per level per glyph doubles its length to say nothing the shapes do not. Underworld development team with AI support from Claude Code
PyVista's default legend face is a triangle for every entry, so the key showed triangles beside wireframes and beside square nodes -- a key that contradicts the figure it is keying, which is worse than no key. Since plot_mesh_hierarchy chose the shapes, it is the thing that can label them, so it now builds its own. Named faces are only triangle / circle / rectangle / none, and a wireframe is none of those: without supplying line geometry a mesh level and a square node key identically and the distinction the figure makes is lost in its own legend. Wireframe and fault-facet entries therefore carry a pv.Line. The key is exposed as plotter._uw_legend_key so it can be INSPECTED. A legend disagreeing with its figure is invisible to any check that counts actors, which is all the previous tests did. Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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