Split-node faults: zero-thickness fault contacts in 2-D and 3-D, parallel, with interface constitutive laws - #502
Open
lmoresi wants to merge 82 commits into
Open
Split-node faults: zero-thickness fault contacts in 2-D and 3-D, parallel, with interface constitutive laws#502lmoresi wants to merge 82 commits into
lmoresi wants to merge 82 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
…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
… 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
…genuine discontinuity split_along_label grows the point chart the way rebuild_without_vertices compacts it: replica vertices for the chain interior, doubled fault facets (<name>Plus keeps the originals' labels, <name>Minus clones them via the clone map), cells kept in source order, edges re-derived so plus-side facets survive with their labels, coordinates repeated, and the star-forest renumbered by the same single broadcast. Tips stay unsplit, so a slip datum tapers to zero there. split_fault wraps the result as a standalone Mesh (no MG tail: the coarse levels do not carry the fault, so an essential condition there would leave the custom-P coarse operator singular) and records the Minus->Plus point pairing that any interface condition needs -- the sides are geometrically coincident, so no coordinate query can ever recover it. Refused loudly rather than mishandled: junctions, loops, single-facet chains, boundary-touching faults, and any fault whose cell fans touch the partition seam. The seam verdict outranks the chain-fragment symptoms each rank sees locally and is allgathered so every rank raises together. Serial tests: exact chart arithmetic (Euler characteristic 1 -> 0), conserved geometry with per-cell data alignment, DOF independence proven through a solve, label pairing into UW_Boundaries, five refusals, and re-application at a second fault position from the same base (the fault moves). Parallel (np=2/3): a rank-interior fault splits with zero star-forest coordinate drift; a seam-crossing fault is refused on every rank. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
Essential-BC DOFs are absent from the global vector, so the rotated path's per-field globalToLocal scatter left them at ZERO in the output fields wherever the Dirichlet datum was non-zero. The solve itself was always right -- the residual path inserts the values internally -- but every field-based diagnostic (projection, Integral, evaluate, renders) read a garbage boundary strip. Measured: far-field sigma_xy of 0.79 for a true 1.00, and a projected stress row off by 3x one fault length out. Invisible to every earlier rotated test because their Dirichlet walls were homogeneous: zero was accidentally the datum. Fix: petsc_dm_insert_boundary_values in cython/petsc_discretisation.pyx -- wrapping DMPlexInsertBoundaryValues, which petsc4py does not expose; the same insertion the consistent-boundary-flux paths already use (issues #407/#411) -- called per field after the scatter in _finalize_rotated_solution. The regression test drives an inhomogeneous lid with rotated free-slip sides and asserts the FIELD carries the wall values. The tell that found it: the divergence theorem, Integral(2 e_xy) against its own boundary term. The rotated boundary's prescribed datum and the sigma_nn reaction recovery were never affected (both live in the global vector). Fixes #497. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…inery For each coincident DOF pair of a split fault, build_rotation now writes an orthogonal 2*dim mean/jump block rotated into the fault (n,t) frame: mean rows on the Plus point, jump rows on the Minus point. The jump-normal row joins the strongly-constrained set (no opening, datum 0) and the slip row stays FREE, which is the zero-shear-traction condition -- its conjugate reaction is the shear traction and an unconstrained row carries none. Both points of a pair are rank-local by the split's seam refusal, so every block sits inside one rank's diagonal portion of Q exactly like the single-node wall blocks, and the whole rotated Newton loop, feasibility projection and line search carry over unchanged. A future datum on the slip row is jump-only prescribed slip; a nonlinear relation on it is friction. fault_contact.py holds the registration (against the pairing recorded by split_fault -- the sides are coincident, so the clone map is the only route to the pairs), the solve driver, and the fault_slip diagnostic that reads the jump through the pairing from the global vector. Validated against the analytic mode-II crack at the gate study's resolution: leak max |[v].n| = 5e-18; emergent slip elliptical to 1.35% RMS shape error with the residual confined to the under-resolved sqrt(r) tip zone; peak 92% of the infinite-medium value Delta-tau a / eta (finite box); Coulomb stress lobes in the same decay family as the inclusion and kinematic-slip representations, with amplitudes ordered exactly as their slips; converged in one Newton increment. Study and figures: ~/+Simulations/fault_split_gate/. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…esign add_viscous_fault_bc(solver, conds, boundary): tau = eta_f * V on the split fault, entering the rotated system as 2*eta_f*M on the slip rows (M the consistent 1-D P2 trace mass per fault facet; the tips drop out topologically because the jump space vanishes where the sides share a point, so the interface law needs no boundary condition of its own). The solve driver gains the two hooks every interface law uses: K.(Qu) added to the rotated residual, K added to a COPY of the ptap-refreshed operator (never injected -- the refresh owns its structure; the null space is re-attached on each rebuilt copy). Measured at the study resolution: the family bridges welded to free monotonically over three decades, following the crack compliance V/V_free = 1/(1 + 0.91 eta_f a / eta) with half-slip at eta_f = eta/a, and machine-zero opening at every member. The welded limit REMOVES the fault (recovers the uncut continuum) rather than stiffening it: only the jump is penalised, the mean velocity is never touched. A friction law replaces the constant with 2*(d tau/d V)*M in the same hook, and the measured compliance curve is its acceptance oracle at the secant eta_f = tau/V. docs/developer/design/FAULT_CONTACT_DEPLOYMENT_2026-08.md records the agreed deployment architecture: the three layers (persistent replicated manifold / ephemeral split mesh with the trace mapping / pair-transform constitutive layer), the two entry paths over one backend, and the parallel seam-vertex crossing rule (pin a vertex at each fault-seam crossing; one shared replica pair per crossing, keyed by root point and side). Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…-sector Three levels, each a legitimate model on its own. J0: disjoint segments with a ligament of one or two local h -- available today under the existing refusals, validated by the King two-fault interaction pattern; a crossing is two offset abutments. J1: true tip-on-fault abutment -- split the through-going fault first, then the abutter with its unsplit tip on the master's slit; the only code change is the boundary-touching refusal distinguishing, for tips only, a domain boundary from a prior fault's slit. J2: the exact degree-d sector split, where branch compatibility is a telescoping identity of sector differences (no cycle constraint, no multiplier) and the acceptance test is machine-precision closure of the branch slips at the junction. Junctions stay off partition seams at every level. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…plain solve() The user-facing layer over the split-node backend, per the deployment design. Mesh.add_fault(faults) runs the whole pipeline in one call — tips placed onto vertices, conforming cut, split with the coincident pairing recorded — accepting a Surface, a (name, points) pair, or a sequence of either. A sequence is a NETWORK: every fault is cut first, then every fault is split, which makes the offset-junction (J0) pattern a one-liner. Splitting renumbers the whole chart, so prior faults' pairings are carried through the new split's point_map rather than copied verbatim — verbatim ids silently index the wrong points, which is exactly how the network case failed first. solver.add_fault_bc(conds, boundary) is the value-first BC: conds = 0 the frictionless contact, conds > 0 the viscous interface tau = eta_f V. With a fault registered, an ORDINARY solve() dispatches to the rotated fault-contact path; guard() and estimate_difficulty() refuse loudly there, exactly as for rotated free-slip. End-to-end tests: one-call fault + plain solve reproducing the crack behaviour, the viscous law through the solver method, Surface-object input, and a two-segment offset network slipping on both faults with machine-zero opening. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…named The blanket refusal (any fault-fan point shared) was far stricter than the invariant requires. What the rank-local split actually needs is that every DUPLICATED vertex owns its whole cell fan — and an UNSHARED vertex always does (any cell touching it from another rank would make it shared), which also keeps every re-homed spoke rank-local. So the rules become: a support-1 facet in a fault vertex's star is the domain boundary only if it is UNSHARED (on a seam every facet looks one-sided locally); and the split refuses only a CHAIN VERTEX on the seam — a seam crossing — with a message naming the crossing milestone. Faults may now run arbitrarily close to a partition seam. np=2: a fault reaching toward the box centre (refused by the old rule) splits cleanly with zero star-forest drift; np=3 lands a chain vertex on a seam and gets the collective crossing refusal. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
The inventory now matches what ships: the en echelon entry carries the tip-lobe receiver, envelope crossings, the P0/C dressing note and the measured per-solve costs; the harness contracts replace the retired symmetric-log recommendation with the standing rendering rules (P0 cells on true connectivity, linear +-1 colour, straight segments over sampled curves — each with the measurement that decided it); the regeneration section quotes today's measured timings; figure count corrected to nine + three animations. Underworld development team with AI support from Claude Code
… traces exact The roughness of curved (polyline-sampled) fault traces is neither the meshing nor the constraint integration: the default per-node normal averages the adjacent facet normals, which zig-zags at the sampling kinks, and the no-opening constraint then forbids smooth slip past each kink — slip notches and normal-traction sawteeth that GROW under mesh refinement. Controlled experiment (circular arc, 78 deg total turn, N = 2..24 segments vs straight control, h and h/2): near-trace cell quality is uncorrelated with the roughness (the straight cut is the worst-quality case yet perfectly smooth), leak <= 1.4e-17 throughout, and supplying the smooth curve's analytic normal on the SAME kinked mesh cuts the traction sawtooth 7-17x, makes it h-independent, and collapses the slip profiles of different samplings onto one curve. - add_fault_bc(conds, boundary, normal=...) with the rotated-free-slip conventions (sympy 1xdim Matrix in mesh.X, or a constant array); every law variant (frictionless/viscous/Coulomb/rate-state) accepts it, and _fault_pair_nodes remains the single frame authority, so the pair blocks, the law rows and every diagnostic see one frame. - Validation before storage; a rejected normal cannot corrupt an already-registered override. - unwrap canonicalises coordinate BaseScalars onto the FIRST mesh's frame, so the stray-symbol check would reject any later mesh's own coordinates: re-tag by name before checking/lambdifying. The same latent bug in rotated_bc's analytic free-slip normal is issue #501. - Regression test: sampled arc, averaged vs analytic normal (the averaged case is the negative control), foreign symbols still refused. - California example rewritten: the San Andreas as ONE continuous dextral trace with a smooth tanh S-bend (the Big Bend, the smoothed stepover) carrying its analytic normal. No beading along the trace; the restraining bend fills with a dCFF compression bowtie exactly where the Transverse Ranges belong. Measured slip +0.242 dextral. - User guide, method write-up and the curation handoff updated (the "straight segments by design" rule is replaced by the mechanism and the capability; deliberately kinked faults should NOT be smoothed). Experiment record: ~/+Simulations/curved_fault_roughness/. Underworld development team with AI support from Claude Code
…study
Two capability increments toward 3-D examples and the parallel
instantaneous-elastic workflow:
- BoxInternalPatch grows patch_cellSize / grading_distance: a gmsh
Distance/Threshold size field grades the mesh from the fault outward
(resolution belongs at the fault, as in 2-D practice). The King L/W
study runs at 36k cells graded where uniform would need ~130k.
- add_fault_bc(..., normal="trace"): the smoothed normal built from the
fault's OWN control polyline (central-difference tangents at control
points, tangent angle interpolated along each segment), for digitized
traces with no analytic formula. add_fault stores each fault's trace
on the split mesh. Validation moved ahead of storage so a rejected
spec cannot corrupt a registered override; regression test extended
with the trace mode (the averaged-normal case remains the negative
control).
Measurements recorded outside the repo (~/+Simulations):
- king_LW_3d: buried rectangular strike-slip rupture, L = 0.5, W in
{0.5, 0.25, 0.125} — leak 0 in all cases; W throttles both the slip
(0.109/0.177/0.225 vs the 0.25 2-D limit) and the Delta CFF reach
(~4-5x at r = 0.3 between L/W = 1 and 4).
- fault_network_parallel: the 7-fault California network in one
add_fault call at np = 1/2/4/8 — all split, all converge, leaks
<= 7e-17, per-fault peak slip rank-count-independent to ~1e-5.
Underworld development team with AI support from Claude Code
…heckpoint fix Two changes that together make the 3-D split-node workflow run-parallel / render-serial, end to end: - split_fault now REDISTRIBUTES the cut mesh before a parallel 3-D split: the patch's cell star, plus one growth layer, is gathered onto the rank that already owns most of it via a shell partitioner; everything else keeps the load-balanced partition. The star marking is SF-propagated (a cell can touch a patch vertex without holding a labelled face in its own closure, so local label reading under-marks near the seam), and one growth layer is exactly sufficient because a point is shared iff its incident cells span ranks. The default partitioner's balance cuts are attracted to the refined patch region — graded fault meshes refused at EVERY np before this. Now: split at np = 2-8 with pair topology identical to serial; measured imbalance 1.8x at np = 8 (the thin star, not the refined band). Several 3-D faults on one parallel mesh are refused loudly (pairing does not yet migrate through the redistribution). - _rebuild_point_sf now installs the rebuilt SF on the COORDINATE DM too (new helper _install_point_sf, all three exits). The coordinate DM is created when the rebuilt chart's coordinates are written — before any SF exists — and its empty snapshot made a parallel HDF5 save write shared vertices as owned on every rank (split mesh at np = 4: 7121 vertex rows vs 6334 true owned); the serial reload of such a checkpoint died in coordinatesLoad. Solves never noticed (field sections are created after the SF). NOTE for the merge with the placed-surface branch: reconnect.py changed here. ptest_0848 rewritten to the new contract: the seam-straddling patch — the case that used to refuse — must now split, single-owner, zero SF drift; np 2/3/4 pass. Serial fault suites unregressed (11 tests). Proof of the workflow: the King L/W W = 0.25 case solved at np = 4 (peak slip 0.1772, identical to serial; leak 0), checkpointed with write_timestep, and rendered serially from the parallel-written checkpoint (~/+Simulations/king_LW_3d/). Underworld development team with AI support from Claude Code
Member
Author
Adversarial reviewWe reviewed this branch as adversaries before asking anyone else to. Six 1. Nodal 2. The held-shut regime has no runtime tell. A zero-strength fault 3. Welded-probe slip shows a small rank dependence. In the 4. The redistribution target keeps its full far-field share. The 5. The unwrap coordinate-retag fix exists only on the fault path. 6. Diff size. 90 files, +16.6k lines, of which roughly a third is 7. The Surface idiom is plumbed but not practised, and crossings are a None of these are correctness defects in what the PR claims; 1 and 2 are Underworld development team with AI support from Claude Code |
Contributor
There was a problem hiding this comment.
Pull request overview
This PR introduces a user-facing split-node fault capability (2-D and 3-D, including parallel workflows) by extending the rotated strong-BC machinery to support coincident DOF-pair contact constraints and interface constitutive laws, alongside supporting meshing, visualisation, tests, and documentation.
Changes:
- Extend rotated strong-BC solving to support split-fault contact pairs and interface-law Jacobian contributions, plus correct field copy-back for inhomogeneous essential BCs.
- Add 3-D meshing support for an internal planar patch (
BoxInternalPatch) and expand visualisation helpers for hierarchies/fault facets and native P1 plotting. - Add extensive serial/parallel tests and developer/user documentation (including figures and teaching example scripts).
Reviewed changes
Copilot reviewed 67 out of 90 changed files in this pull request and generated 1 comment.
Show a summary per file
| File | Description |
|---|---|
| tests/test_1018_rotated_freeslip.py | Adds regression test for rotated-path field copy-back with inhomogeneous Dirichlet walls. |
| tests/test_0847_mesh_hierarchy_plot.py | Adds rendering-behaviour tests for plot_mesh_hierarchy (actors, faults, glyphs, legend). |
| tests/test_0847_fault_api.py | Adds end-to-end tests for the user fault API (add_fault, add_fault_bc, normals). |
| tests/test_0846_visualisation_native_mesh.py | Tests that continuous P1 fields render on native mesh connectivity (not Delaunay). |
| tests/test_0845_relax_pinned_band.py | Adds serial tests for relax(pin_bands=...) behaviour and boundary pinning. |
| tests/test_0843_edge_split_adapt.py | Adds serial tests for engine="edge_split" refinement properties and confluence contract. |
| tests/test_0840_nvb_3d_serial_adapt.py | Refactors/strengthens 3-D adapt MG-level assertions using shared ladder check. |
| tests/test_0836_nvb_graded_adapt.py | Refactors MG-level assertions and adds tests for mg_coarsening_ratio effect. |
| tests/parallel/ptest_0848_fault_split_3d_parallel.py | Adds 3-D parallel fault split tests including redistribution and SF drift checks. |
| tests/parallel/ptest_0845_relax_pinned_band_parallel.py | Adds parallel tests for partition-independent pin sets and pinned-vertex immobility. |
| tests/parallel/ptest_0845_fault_split_parallel.py | Adds parallel tests for split-node faults (collective refusal, seam/crossing cases). |
| tests/parallel/ptest_0843_edge_split_parallel.py | Adds parallel confluence tests for engine="edge_split" and adapt-tail checks. |
| tests/_mg_ladder.py | Introduces shared MG coarsening “ladder” assertion helper for adapt children. |
| src/underworld3/visualisation/init.py | Exposes new/updated visualisation helpers and constants in the public API. |
| src/underworld3/utilities/rotated_bc.py | Extends rotation assembly for fault pair blocks; adds interface assembly/tangent; fixes copy-back gap. |
| src/underworld3/utilities/nvb.py | Adds explicit TODO note about known 3-D prolongation inaccuracy. |
| src/underworld3/utilities/init.py | Exports new utilities modules (edge_split / fault_contact / fault_split / etc.). |
| src/underworld3/meshing/cartesian.py | Adds BoxInternalPatch mesh factory for embedded internal planar patches in 3-D. |
| src/underworld3/meshing/init.py | Re-exports BoxInternalPatch in meshing public namespace. |
| src/underworld3/cython/petsc_generic_snes_solvers.pyx | Adds add_fault_bc and routes fault-contact solves through rotated solve path; updates guard probes. |
| src/underworld3/cython/petsc_discretisation.pyx | Adds Cython shim petsc_dm_insert_boundary_values wrapping DMPlexInsertBoundaryValues. |
| docs/developer/subsystems/conforming-surfaces-and-fault-zones.md | New subsystem doc explaining conforming surfaces, zones, and implications. |
| docs/developer/index.md | Adds new subsystem doc to developer documentation index. |
| docs/developer/design/figures/split-node-faults/stack-progression.typ | Adds Typst/cetz source for split-node fault method figure. |
| docs/developer/design/figures/split-node-faults/split-anatomy.typ | Adds Typst/cetz source for split anatomy figure. |
| docs/developer/design/figures/split-node-faults/split-anatomy-data.json | Adds generated geometry data for split anatomy figure. |
| docs/developer/design/figures/split-node-faults/render-components-3d.py | Adds PyVista script to render 3-D split components figure from real mesh. |
| docs/developer/design/figures/split-node-faults/pair-transform.typ | Adds Typst/cetz source for pair-transform schematic. |
| docs/developer/design/figures/split-node-faults/grid-hierarchy.typ | Adds Typst/cetz source for grid hierarchy schematic. |
| docs/developer/design/figures/split-node-faults/grid-hierarchy-data.json | Adds generated geometry data for grid hierarchy schematic. |
| docs/developer/design/figures/split-node-faults/generate-stack-progression-data.py | Adds generator for stack progression figure geometry. |
| docs/developer/design/figures/split-node-faults/generate-split-anatomy-data.py | Adds generator for split anatomy figure geometry. |
| docs/developer/design/figures/split-node-faults/generate-grid-hierarchy-data.py | Adds generator for grid hierarchy figure geometry. |
| docs/developer/design/FAULT_CONTACT_DEPLOYMENT_2026-08.md | Adds deployment/architecture design note for split-node fault contacts. |
| docs/developer/ai-notes/fault-teaching-examples-handoff-2026-08.md | Adds handoff note for teaching examples and figure regeneration workflow. |
| docs/advanced/index.md | Adds new advanced-doc entries for split-node faults and teaching examples. |
| docs/advanced/figures/fault-examples/orientations.py | Adds teaching script: orientation sweep vs slip behaviour. |
| docs/advanced/figures/fault-examples/mohr_graded.py | Adds teaching script: hydrostatic loading sampled along welded faults. |
| docs/advanced/figures/fault-examples/mohr_friction.py | Adds teaching script + animation: Coulomb friction envelope vs Mohr circle. |
| docs/advanced/figures/fault-examples/mohr_cohesion.py | Adds teaching script + animation: cohesive Mohr-Coulomb behaviour. |
| docs/advanced/figures/fault-examples/mohr_circle.py | Adds teaching script: Mohr circle measured by welded faults. |
| docs/advanced/figures/fault-examples/mohr_animate.py | Adds teaching script + animation: Mohr circle build as fault rotates. |
| docs/advanced/figures/fault-examples/ladder.py | Adds teaching script: “fault-strength ladder” comparing multiple laws. |
| docs/advanced/figures/fault-examples/.gitignore | Ignores generated frame/log outputs for teaching example renders. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+451
to
+457
| from underworld3.cython.petsc_discretisation import \ | ||
| petsc_dm_insert_boundary_values | ||
| for name, var in solver.fields.items(): | ||
| sg = U.getSubVector(solver._subdict[name][0]) | ||
| solver._subdict[name][1].globalToLocal(sg, var.vec) | ||
| U.restoreSubVector(solver._subdict[name][0], sg) | ||
| petsc_dm_insert_boundary_values(solver._subdict[name][1], var.vec) |
…onverted The two pre-release follow-ups from PR #502, delivered: - The 3-D pipeline now consumes uw.meshing.FaultSurface end to end: FaultSurface.rim_polygon() extracts the ordered rim of a planar triangulated surface (boundary-edge chain; non-planar refuses toward the discrete-entity embed); BoxInternalPatch(patch_points=fault_surface) takes the object, adopts its name and stores it on the mesh; split_fault carries it (and the 2-D stored traces) onto the split child; add_fault_bc(..., normal="surface") builds the constraint frame from the surface's OWN face normals (nearest-face lookup — exact for a planar patch, unchanged for future curved sheets). - uw.meshing.prepare_fault_network converts an imported 2-D trace set into splittable offset-junction form, loudly: X crossings cut both traces, T abutments trim the abutting end, near-miss endpoints are pulled to clearance. Pull-backs are angle-corrected (ligament / sin theta — the Euclidean-clearance check caught the along-trace version being short at oblique junctions) and pieces shorter than two ligaments are dropped and reported. - Ligament sensitivity measured (T-junction, 35 degrees, both faults frictionless, ~/+Simulations/fault_junction_ligament/): branch peak slip is 1.01x / 0.96x / 0.79x its isolated value at ligament = 1h / 2h / 4h, and at FIXED ligament 0.04 the whole branch slip profile is IDENTICAL at h = 0.02 and h = 0.01 — the offset junction is a converged physical answer for a given plug size, not a resolution artifact. Keep ligaments at 1-2 local cell sizes; a true shared-vertex junction design is an accuracy refinement, not a correctness necessity. Tests: test_0849_fault_network_prep (X/T/near-miss conversion, Euclidean clearance, end-to-end add_fault) and test_0848::test_fault_surface_route (rim extraction, name adoption, surface-normal compile, non-planar refusal). User guide updated to lead 3-D with the FaultSurface idiom. Underworld development team with AI support from Claude Code
Member
Author
|
Two of the pre-release follow-ups have now landed on the branch (and in this PR): the FaultSurface object is the 3-D input idiom end to end (rim embed, stored on the mesh, The ligament approximation is now measured, not asserted: at a 35-degree T-junction with both faults slipping, the branch's peak slip is 1.01x / 0.96x / 0.79x its isolated value for ligaments of 1 / 2 / 4 local cell sizes, and at fixed ligament 0.04 the entire branch slip profile is identical at h = 0.02 and h = 0.01. The offset junction is a converged physical answer for a given plug size; keep ligaments at 1-2 local cells. A true shared-vertex junction design is therefore an accuracy refinement, staged, not a correctness gap. Underworld development team with AI support from Claude Code |
Mesh.add_fault now stores Surface OBJECTS on the split mesh (alongside the raw traces), and normal="surface" resolves in 2-D to the smoothed normal of the Surface's own control polyline — the same spelling and the same object-first idiom as the 3-D FaultSurface route. Deliberately NOT the Surface's director (signed-distance gradient): that field is discontinuous exactly ON the trace, where the fault nodes live; the director stays the right tool for TI weak zones, the polyline normal for the contact frame. Test extended with the frame-equality check. Underworld development team with AI support from Claude Code
The prepare_fault_network capability shown end to end: raw traces that genuinely intersect (a dextral trunk, a splay abutting its midpoint at 30 degrees, a conjugate crossing outright) are auto-converted to offset junctions and split in one add_fault call; the trunk and splay rupture together, the conjugate is welded. The Delta CFF map's brightest features are the junction plugs themselves — the intact ligaments take the concentrated load shed by their slipped neighbours, which is where a through-going event would break next — shown at map scale and in ligament-scale zooms at the T and X junctions. Page section added with the measured ligament-convergence statement. Underworld development team with AI support from Claude Code
prepare_fault_network was cutting the through-going trace at T abutments, contradicting its own documentation (the trunk in the branching example was interrupted twice). Now: a T abutment trims only the abutting end; X crossings cut both traces by default; and through=[names] declares MASTER faults that stay continuous even at crossings (the other trace yields on both sides; two masters crossing is a hard error). Tests updated to the corrected contract. The cost of interrupting a fault is now measured, not guessed (branching_compare.py, on the teaching page): with the trunk through-going its peak slip is 0.2590 — slightly MORE than the isolated trunk's 0.2558, because the slipping splay feeds it; cut at the crossing it drops to 0.1963 (about a quarter forfeited) and is pinned at the junction plug, but the profiles converge away from it — an interruption's reach is the segment scale, not the system scale. The both-cut state the old defect produced was the worst of all (0.144). branching.png regenerated under the corrected default so the committed figure matches the code's behaviour. Underworld development team with AI support from Claude Code
Does a near-miss tributary reproduce a genuine branch? Two findings, one per axis. The GAP is second-order: with the trunk continuous and the splay abutting, ligaments of 1h/2h/4h give the same slip on every arm except the splay's near-junction toe, which moves with the trimmed tip. The DECOMPOSITION is first-order: routing continuity through the bend instead (west arm + splay as one kinked fault) locks at the 33-degree kink — slip through a kink is geometrically incompatible with no-opening, the curved-fault mechanism appearing as genuine physics — costing both trunk arms a third of their slip. The rule for imported networks: declare the straightest path continuous; never route continuity around a corner. The continuous-trunk decomposition is admissible for the true junction and differs from it only in the ligament-scale toe, which is all a true shared-vertex junction design would buy. Underworld development team with AI support from Claude Code
Delta CFF maps for both decompositions of the Y-branch, map scale plus branch-point zooms at one colour scale. Trunk-continuous: the stress shadow runs straight through the junction, which carries no feature at all — mechanically invisible, as the declared through-going pair should be. Bend-continuous: the locked kink is a barrier — both arms' shadows terminate against it and dump load into the corner, the rupture-arresting-bend signature. Noted in the log: decomposition B's welded reference landed 139 stress units away in pressure gauge from its slipping solve — the far-field anchor removed it (printed), the strongest case yet for the anchoring discipline. Underworld development team with AI support from Claude Code
Three chains terminating at one shared vertex are refused when the second split's tip lands on the first fault's slit — but the message said "touches the domain boundary", which is what the predicate sees (one-sided facets in the tip's star), not what the user did. The refusal now detects that the offending facet carries a Plus/Minus side label and says so: a shared point would clamp every arm's slip to zero there — the kink-lock in its purest form, STIFFER than a true sector junction — and points at prepare_fault_network, whose offset form brackets the true branch (the true-branch teaching example). Test added; deliberately NOT relaxed to allow the touching case, since its mechanics is reproduced to ligament accuracy by the continuous-bend decomposition already measured. Underworld development team with AI support from Claude Code
Member
Author
Adversarial review: PR #502 — split-node faults (fault delta over #488)Scope: the 53-commit fault delta VerdictThe core mechanism is sound and well-tested: the split is a clean rebuild with Measurements
MERGE-BLOCKERS
Findings (non-blocking, severity order)
Contract checks that PASSED scrutiny
Rebase notes (08f8603 supersession)
|
…rossing machinery The shipped 2-D seam-crossing capability did not survive measurement: the cut drags partition-seam vertices onto the fault line, so the crossing design's premise (one crossing vertex, both flanking facets rank-local) fails and the shipped ptest geometry died at np=2 with a chain-fragment error and skipped at np=4 (PR #502 review, blocker 1). Rather than repair the most intricate code in the delta, give 2-D the strategy 3-D already ships and measures healthy at np=2-8: BEFORE the split, move the fault's cell star plus one growth layer to the rank that already owns most of it (shell partitioner; everything else keeps the balanced partition), so the chain is rank-interior and the split runs with serial topology. _redistribute_fault_interior is now dimension-generic and takes a list of (label, value) pairs; add_fault redistributes ONCE for a whole network, keyed on every fault together, which is also what keeps prior pairings valid (a pairing cannot yet migrate through a redistribution — split_fault still refuses that per-split case, now in both dimensions). In 2-D the move only happens when the chain actually touches the seam (_fault_labels_touch_seam, collective verdict); 3-D stays unconditional. The crossing apparatus this makes unreachable is deleted: the crossings fan walk in _take_sides, the _sf_sum census and third-rank corner rule, _pull_seam_vertices_onto_crossings, and the dup_new keyed star-forest exchange in reconnect._rebuild_point_sf (back to the simple renumbering broadcast; _install_point_sf stays). The low-level splitter now refuses ANY shared chain vertex, checked FIRST so a straddling fault always gets the honest seam verdict instead of a fragment symptom — the misdiagnosis the review measured. ptest_0845's crossing test now asserts SUCCESS (split + slip, machine-zero leak) with no tolerated refusals — measured 4 passed at np=2 and np=4, 3 passed + 1 documented low-level skip at np=3, peak slip rank-independent to ~1e-8 against serial. The organic straddling cut keeps its collective low-level refusal and now also asserts that add_fault splits the same geometry cleanly. 3-D ptest_0848 stays 9/9 at np=2 and np=4 with the generalised helper. Underworld development team with AI support from Claude Code
fault_normal_traction(solver, "A", info) on a two-fault network returned the UNION of both faults' pair nodes (measured: 44 rows where fault A has 23) — fault B's nodes interleaved into A's traction profile as bogus arc-length samples — and update_fault_state mixed the faults the same way in its slip vector, its theta store (foreign point keys written into this fault's by_point) and its convergence monitor (PR #502 review, blocker 2). The cause: _InterfaceAssembler always walked EVERY registered law fault plus the include fault. When include is given the assembler now restricts its walk to exactly those faults — a per-fault view. The SOLVE path passes no include and still assembles every law-carrying fault, untouched. Single-fault runs are bit-identical, which is why the suite never saw the bug. The new regression test (two-fault J0 network, viscous laws on both) asserts the traction profile's row count against _fault_pair_nodes — the per-fault oracle that was always correct — the arc-length span, and that a theta update for fault A writes only fault A's own points. Validated as a negative control: it fails on the unfixed assembler with exactly the measured union count (44 vs 23) and passes after the fix. Also freshens comments that still described the retired 2-D seam-crossing machinery. Underworld development team with AI support from Claude Code
Three raises inside collective paths were rank-local: the pair-straddles-ranks RuntimeError and the fault-node-on-a-wall- boundary ValueError in build_rotation's pair-block loop, and the pairing/closure-mismatch and vanishing-analytic-normal raises inside _fault_pair_nodes (reached from both build_rotation and the interface assembler). A rank that raised alone left its peers blocked in the next collective — Q assembly, the trace-mass exchange, the Newton loop — turning an error into a deadlock at np > 2 (PR #502 review, finding 3). Funnel them through the allgather-verdict pattern fault_split already uses: _collective_raise gathers each rank's local error string and raises the SAME error on every rank if any is non-empty. build_rotation wraps its whole fault block (fault_names is registration state, identical across ranks, so every rank reaches the verdict exchange); the assembler collects _fault_pair_nodes failures per fault and verdicts after its walk, before the collective mass exchange. The straddle check should now be structurally impossible — the redistribution makes every pair rank-local — but checks that "cannot fire" are exactly where this PR's blocker 1 lived, so it stays, and its verdict is now collective like the rest. Underworld development team with AI support from Claude Code
Initialised in the solver setup, never written or read anywhere in the repository — the interface laws actually live in solver._fault_interface_laws, set lazily by utilities/fault_contact.py (PR #502 review, finding 6). Deleted, and the neighbouring comment now points at the real home of the laws. Underworld development team with AI support from Claude Code
Member
Author
|
Response commits 1d48731..c0d64cd — both blockers fixed, plus findings 3 and 6:
Serial suite 28/28 (the new regression test included). Remaining known work before merge: this branch is stacked on #488 (which owns the two failing CI tests) and carries the superseded 08f8603 form of the #500 fix — both resolved at the rebase once #488 lands. Underworld development team with AI support from Claude Code |
The junction policy Louis set: declared master faults keep geometric continuity; every other junction is left as an offset gap carrying DAMAGE-ZONE material, and the stress lobes decide how the faults link. Two pieces land here: - prepare_fault_network(..., return_junctions=True) hands out every junction's kind, point, pull-back and fault pair. - uw.meshing.damage_zone_yield(mesh, junctions, tau_damage, radius) builds the composite yield-stress expression over all plugs (SHARP Piecewise regions folded with Min against a sane finite far cap). The recipe behind it is measured (~/+Simulations/fault_junction_rheology): a von Mises plug's tau_J dial is smooth over a decade and interpolates the geometric-decomposition brackets; a collinear gap grows a straight, sharp-edged yielded band tip-to-tip (verified at 6 cells across the gap with adapt-on-top refinement through the gap); TWO Picard passes match picard=8 to 0.01% (find the lobe, respond); blending a rheological parameter against a huge sentinel through any smooth mask tail contaminates the plug — hence the sharp regions; the fully-weakened band (~100x) locally recreates the thin-weak-inclusion conditioning cost, so production plugs should cap weakening near 10x. Newton via the consistent tangent on the rotated path is issue #507 (pre-existing, wall-only reproducer, NOT a split-node defect). Underworld development team with AI support from Claude Code
…eature/fault-split-node
Measured side by side on the collinear-bridge strain rate (the maintainer's challenge): continuous-P1 nodal projection renders smoothly with no ringing — the old "never P1" rule was calibrated against the kinked-normal sawtooth fields later cured at the source, and the split mesh's P1 space is naturally discontinuous across faults (doubled nodes), so the jump is representable. P0 remains right where the field carries genuine cell-scale structure: the plastic yield-zone boundary is sharp in P0 and blurred a node spacing by P1. Rule now: choose by physics; components recovered separately, invariants in numpy; one-sided Oranges on white for strain rate with mesh edges visible (the maintainer's convention, restored). Underworld development team with AI support from Claude Code
Underworld development team with AI support from Claude Code
…split-node # Conflicts: # src/underworld3/utilities/rotated_bc.py
…#507) Differentiating an unwrapped strain-rate invariant produces half-integer powers of (grad v : grad v) whose value or derivative is 0/0 at a state of rest, so EVERY consistent-tangent assembly with edot = 0 anywhere filled the operator with NaN — measured J(0) norm = nan for a ViscoPlastic model at ANY yield stress, surfacing as GAMG's "Computed maximum singular value as zero" (error 77) on the standard, rotated free-slip, and split-node fault paths alike. Mainline, since the consistent-tangent feature existed. Why it hid, and why detection cannot fix it (both measured): - Layer 1's automatic cold-start Picard injection fires correctly and STILL fails on boundary-driven problems: one nrichardson sweep propagates boundary data a single element layer, leaving the deep interior at exactly zero strain rate. Body-force-driven problems (the yield campaigns) fill F(0) everywhere — the only class the protection ever actually covered. - The "continuation" tangent's alpha = 0 phase does not protect: the blended kernel evaluates the Newton branch pointwise and IEEE 0*NaN = NaN. - A rigidly-translating stuck region has edot = 0 at the CONVERGED solution: the singular state is physics, not a start-up artifact, so no warm-start policy can make it unreachable. The fix implements the derivative's removable-singularity limit: _jacobian_unwrap now adds 1e-36 under every half-integer power whose argument carries unknowns (+1/2 the invariant, -1/2 its reciprocal in eta_pl, -3/2 their derivatives). Jacobian sources only — the residual never passes through this function and the default Picard tangent never calls it, so both stay bit-identical by construction. The perturbation at any resolvable strain rate is under one part in 1e24; the cold Jacobian of a never-yielding probe now equals its smooth-state Jacobian to all digits. The Layer-1 warm-up remains as the convergence aid it was designed to be; its comment now states the measured facts. Gates: 53 passed across yield-homotopy, rotated free-slip, fault contact/API, plus the new cold-start regression test (test_1067_newton_cold_start: True and "continuation", solve from cold and assert J(0) finite). Every previously-failing reproducer (standard / rotated wall / fault, True / continuation) now solves. 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.
What this is
Zero-thickness split-node faults for UW3, in 2-D and 3-D, as a user-facing
capability. A fault becomes a genuine velocity discontinuity in the mesh:
the labelled facet chain (2-D) or triangulated patch (3-D) is split by
duplicating its interior nodes, a strong no-opening constraint ties the
coincident pairs together, and interface constitutive laws set the
tangential condition. There is no thin weak band and no viscosity contrast
for the solver to fight: the conditioning benchmark measures 10 outer
Schur iterations against 147 for the equivalent one-element weak inclusion
at 1e-4 contrast on the same mesh.
The user surface
Mesh.add_faultplaces tips and kinks onto mesh vertices, cuts aconforming chain, splits it, and records the coincident DOF pairing.
Lists of faults build offset networks in one call.
add_fault_bc(conds, boundary):0= frictionless,eta_f= viscousinterface;
add_coulomb_fault_bcandadd_rate_state_fault_bcfor thefriction family. Laws are sympy expressions in canonical symbols
(slip rate, effective normal stress, state); consistent tangents come
from
sympy.diff, never hand-coded. The effective normal stress isSIGNED and reaction-fed; each law clamps its own strength (a cohesive
law keeps strength into mild tension, a cohesionless one does not).
normal=on every law variant supplies a smooth fault normal — ananalytic sympy expression, or
"trace"to build it from the fault's ownpolyline. This matters for sampled curved traces: the default
facet-averaged normal zig-zags at sampling kinks and produces slip
notches and traction sawteeth that grow under refinement; the smooth
normal removes them (measured 7-17x) and restores h-convergence.
(
fault_slip,fault_pair_jumps,fault_normal_traction) — the sidesare geometrically coincident, so coordinate queries cannot see them.
BoxInternalPatchembeds a planar interior patch,now with
patch_cellSize/grading_distanceto grade resolution fromthe fault outward.
Parallel
vertex, manufactured automatically; a fault along a seam refuses
loudly and collectively. Swept at np 2-5. A seven-fault network splits
and solves at np 1-8 with per-fault slip rank-independent to ~1e-5.
split_faultredistributes first — the patch's cell star plus onegrowth layer moves to the rank that already owns most of it (shell
partitioner, SF-propagated marking), everything else keeps the balanced
partition. Splits at np 2-8 with pair topology identical to serial;
measured imbalance 1.8x at np 8 (the thin star, not the refined band).
Several 3-D faults on one parallel mesh are refused (pairing migration
through the redistribution is future work).
coordinate DM as well (
_install_point_sf) — previously a parallelHDF5 save of any rebuilt mesh wrote shared vertices as owned on every
rank and the file would not reload. The run-parallel / render-serial
workflow (solve at any np,
write_timestep, read checkpoints seriallyfor analysis) is exercised end to end by the King study scripts.
Validation
converging to the (8/3 pi)(d tau/eta) sqrt(a^2-r^2) oracle from below;
leak at machine zero (1e-17 or exactly 0) in every benchmark.
stick/slide switching with reaction-fed normal stress (2-D and 3-D);
rate-state smoke tests with exact ageing updates.
test_0845(split topology),test_0846(contact physics),test_0847(user API, analytic/trace normals with negative control),test_0848(3-D), plus parallelptest_0845/ptest_0848.docs/advanced/split-node-faults.md),method and benchmark write-up
(
docs/developer/design/SPLIT_NODE_FAULT_METHOD_2026-08.md), afault-mechanics teaching page with nine committed figures and three
animations, and the deployment design doc.
Known limitations (all refuse loudly)
Closed-loop faults; daylighting (faults reaching the domain boundary);
junctions sharing vertices; multiple 3-D faults on one parallel mesh.
Split meshes carry no geometric-MG tail (coarse levels do not contain
the fault) — solvers take algebraic-multigrid defaults. The bilateral
no-opening constraint holds a fault shut where it would physically open;
the tell is a tensile fault normal traction, documented in the guide.
Coordination
This merge brings
reconnect.py/line_cut.py(the cut and rebuildmachinery) to development for the first time. The placed-surface branch
also modifies
reconnect.py; per the agreed protocol, whichever branchmerges second reconciles — this PR merging first means the
placed-surface session rebases onto development and folds
_install_point_sfinto its rebuild paths.An adversarial review is posted as a comment below, per our review
practice.
Pre-release follow-ups (agreed direction, not in this PR)
Before importing generic fault models: (a) lead with the
uw.meshing.Surfaceidiom in examples/docs (accepted and tested today;the examples use the
(name, points)shorthand) and wireFaultSurfaceinto the 3-D split path; (b) source
normal=from the surface objects'own geometry rather than user formulae; (c) an importer that detects
crossing/abutting traces and auto-converts them to the offset-junction
(ligament) form, loudly — true shared-vertex junctions need their own
design (a crossing vertex takes four coincident copies and a non-binary
pairing) and are staged behind a ligament-sensitivity measurement.
Underworld development team with AI support from Claude Code