Fault networks: the network mesher, one specification with two realisations, and the junction glue

.claude/skills/nonlinear-solver/SKILL.md

@@ -77,7 +77,7 @@ case felt like whack-a-mole. Check them first.
 
 | Trap | Symptom | Fix |
 |---|---|---|
-| Consistent Newton makes the velocity block **non-symmetric**; a Chebyshev/Richardson MG smoother assumes SPD | smoother diverges / stalls → `DIVERGED_LINEAR_SOLVE` or an endless grind | **Now the default** — the FMG bundle ships `mg_levels_ksp_type=gmres` + `pc_type=sor` + `norm_type=none`. Only an issue if you override it, or on GAMG (which uses PETSc's chebyshev default) |
+| **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: a small δ soft-min (`yield_mode="softmin"`, powermean, `yield_anchor="yield"`), a rounded viscosity floor, or rate-strengthening ξ; Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way |
 | `preconditioner="fmg"` (vs explicit `pc_type=mg` + manual mg opts) | outer KSP "converges" in **1 iteration** → no real Newton correction → stall → `DIVERGED_LINE_SEARCH` | use explicit `pc_type=mg` with the smoother opts above; bound the outer KSP (`ksp_max_it`~80) so a hostile step fails fast |
 | Cold plastic start `v=0`, or any rigid/unyielded point | `DIVERGED_FNORM_NAN` at iteration 0 | **Not** a div/0: `ε̇=0` gives `η_pl=+inf`, which `Min` and the sqrt soft-min carry correctly to the viscous branch. Only a soft-min form that computes `η_ve·η_pl/(η_ve+η_pl)` breaks (`inf/inf`). Fixed in the power-mean; if you hand-roll a blend, write the harmonic mean as `η_ve/(1+η_ve/η_pl)`. **Do not reach for a strain-rate floor** — it hides this rather than fixing it |
 | LU velocity block with all-Dirichlet-ish BC | pressure nullspace singular | attach the Stokes nullspace / avoid a bare LU there |

docs/advanced/fault-networks.md

@@ -15,24 +15,107 @@ net = uw.meshing.FaultNetwork(
     [("Main", main_pts), ("Splay", splay_pts), ("Cross", cross_pts)],
     hierarchy=["Main", "Splay", "Cross"])   # seniority order
 
-mesh = net.prepare(h=0.006).build()          # junctions -> mesh -> split
+mesh = net.prepare(h=0.006).build(width=0.01)  # junctions -> mesh -> split
 
 v = uw.discretisation.MeshVariable("V", mesh, 2, degree=2)
 p = uw.discretisation.MeshVariable("P", mesh, 1, degree=0,
                                    continuous=False)
 stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
-stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
-stokes.constitutive_model.yield_mode = "min"
-stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
-stokes.constitutive_model.Parameters.yield_stress = \
-    net.damage_yield(v, dial=0.05)           # the junction glue
-stokes.consistent_jacobian = True
-net.apply_contact(stokes)                    # no-opening pairs, all pieces
+stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = \
+    net.junction_patch(eta_0=1.0)            # the junction glue (linear)
+net.apply(stokes)                            # no-opening pairs, all pieces
 # ... wall boundary conditions ...
 info = net.solve(stokes)
 print(net.slips(stokes))                     # peak slip per piece
 ```
 
+## One specification, two realisations
+
+A fault is specified once — a trace, its rank in the hierarchy, and the
+properties it carries — and then *realised*. Which realisation you get
+is a keyword on `build`, not a different set of calls:
+
+```python
+net.prepare(h=0.006)
+mesh = net.build(width=0.01)                      # cut, node-pair contact
+mesh = net.build(width=0.002, realisation="ti")   # volumetric weak plane
+net.apply(stokes, eta_1=0.01)                     # eta_1: TI only
+```
+
+Both realisations place the same ribbon band along the same prepared
+pieces, so the cells are identical and results from the two may be
+compared directly. The band is meshed around the trace's own points and
+segments, which become mesh vertices and edges, so **the mesh can be cut
+whatever the width is** — measured complete, with exact vertex
+coincidence, down to a band a tenth of the background element size. The
+realisation is a free choice, not something the mesh grants or refuses.
+
+What differs is what `width` *means*. For the split it is a resolution
+parameter: the band exists to give the cut its own vertices, and its
+thickness is not physics. For the weak plane it is constitutive — the
+layer thickness that sets the slip rate through `V = 2 e_nt w` — so it
+wants two or three elements across it. That is the whole asymmetry, and
+it is about the rheology rather than the mesh.
+
+`slips()` reports each realisation in its own quantity: the tangential
+jump between the two nodes of a cut pair, or the jump in tangential
+velocity across the layer, sampled one half-width plus a cell either
+side of the spine. Both are the fault's own throughput; a probe placed
+further out reads the surrounding flow as well and over-reads short
+strands.
+
+`build(width=None)` keeps the older no-band path — graded refinement
+cut directly. It is split-only, and its mesh is not the one a weak
+plane would use, so do not compare across that choice.
+
+**The band is material, not scaffolding.** It is easy to read the band
+as something the weak plane needs and the split merely tolerates. It is
+not: the band is a meshed region of material *around* the fault, and a
+segmented fault does its interesting work exactly there — at the strand
+tips, and in the ligaments where one cut stops short of the next. Damage
+in those places needs cells to live in, and the band is where they are.
+`net.band` is the mask, `net.footprints` the per-strand ones, in either
+realisation.
+
+`net.band_yield(tau_y)` gives a rheology for the whole band: von
+Mises yield confined to it, everything outside set far too strong to
+yield. Read the next paragraph before using it as glue.
+
+```python
+stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
+stokes.constitutive_model.Parameters.yield_stress = net.band_yield(4.0)
+stokes.consistent_jacobian = True
+```
+
+A released fault flank sits far below `tau_y` and is untouched, so
+the breakdown appears where the mechanics puts it — but that is tips
+and bends as much as joints. Measured on the S-fault rig: at the
+strength that repairs a stepover, more than half the yielded band cells
+were on strand flanks and free tips, and one step weaker the whole
+main strand had become a weak fault. A uniform threshold cannot pick
+out the welds alone, because the stress concentration at a weld is not
+far enough above the tip and bend concentrations. `band_yield` is a
+damage model for the band; the junction glue is `junction_patch`.
+
+The two realisations' interface parameters correspond, which is worth
+keeping in view when comparing them: the zero-thickness limit of a band
+of viscosity `eta_band` and width `w` is an interface viscosity
+`eta_f = eta_band / w`, which is the `conds` argument of
+`add_fault_bc`. The weak plane's `V = 2 e_nt w` is precisely what the
+contact replaces with a genuine slip rate.
+
+**Properties belong to the fault.** `net.surface(name)` returns the
+retained {class}`~underworld3.meshing.surfaces.Surface` for a piece.
+Friction, accumulated slip, a damage state live there, on the fault,
+and outlive any one realisation of it:
+
+```python
+main = net.surface("Main")
+friction = main.add_variable("mu", size=1)
+friction.data[:] = 0.6
+```
+
 ## The recipe, and why each piece is the way it is
 
 **Hierarchy.** At an X crossing the senior fault runs through and the
@@ -48,17 +131,70 @@ same answer when the junction patch is refined 2x. Make the join as
 small as the mesh allows; buy fidelity with elements, not physical
 size.
 
-**The glue.** `damage_yield` places a compact viscoplastic plug at
-each junction: yield `dial * (1 + 2 * edot_II)` inside, effectively
-rigid outside, sharp `Piecewise` boundaries. The strength and the
-rate-regularisation move together on ONE dial (separating them makes
-the solve harsh without making the zone weaker). Zone stress is
-proportional to the dial down to a ~100x viscosity contrast with
-Newton-from-cold still converging — the compact plug conditions like a
-hole, not like a thin weak layer, so the classic thin-inclusion
-Schur breakdown never appears. `dial=0.05` is near-invisible in the
-stress field at unchanged cost; `dial=0.01` reaches the transmission
-ceiling of an inviscid plug at roughly double cost.
+**The glue, and where it goes.** The split only goes wrong at the
+joints: away from them the cut *is* the target every volumetric
+representation converges to, and adding weakness along a whole strand
+makes the fault over-weak in a way that depends on the band width. So
+the glue is placed, not found. `junction_cells()` reads the places off
+the mesh itself: the ribbon (the band with its extrapolated margins)
+is everything the weak-plane realisation would treat as fault, the cut
+chains are what the split sliced, and a band cell whose nearest spine
+point lies in a piece's margin *and* which sits inside a second
+piece's ribbon is where two pieces meet without being joined — a
+kissing branch, an abutting pair, the intact bridge of a stepover.
+Free tips are excluded on purpose: a margin that runs into intact
+material is a tip, and damage there lengthens the fault instead of
+joining it (measured: with the free tips included, nearly every
+yielded cell was at a tip and the main strand grew 1-6% longer in
+slip). The cells are dilated by one vertex ring, and that ring is
+not optional: the weld's stiffness lives in the intact material
+around the two tips, and the bare junction cells recover only a
+fifth to a quarter of the deficit even when fully plastic.
+
+Two pieces that continue one another along a line — an abutting pair,
+a stepover's continuation — are placed on **one spine**: two ribbons
+laid along the same line interleave their vertices into sliver cells
+(measured: 7800 cells below 1e-6 in area, and the velocity solve
+five times slower). `build()` groups such pieces (end tangents within
+25 degrees, the far start within half a width of the line, within the
+margins' reach), bridges the gap with spine vertices at the local rung,
+and cuts each piece at its own ends; the gap is spine the split does
+not cut, which is exactly what the junction rule reads.
+
+For the rule to see a joint, the ribbons have to meet across it.
+`build()` sees to that: at an end that sits on a prepared junction the
+tip margin is extended until the ribbon reaches the other piece's cut,
+so the whole ligament lies in both ribbons; free tips keep the default
+margin. An abutting pair that `prepare()` did not record as a junction
+(a gap wider than the ligament) is covered as far as the default
+margins overlap — a gap wider than that is two faults, and stays
+welded, which is what a gap of intact rock means.
+
+`junction_patch(eta_0, ratio=0.01)` then makes those cells weak
+isotropic material, `eta = ratio * eta_0`. A viscosity ratio rather
+than a yield stress, because the joint only has to be broken and a
+ratio needs no stress scale — nothing about the block or the loading
+has to be known to set it. Measured against the two end members on
+the S-fault rig (the fault *longer*, one continuous cut, and the fault
+*cut*, abutting cuts, at two resolutions): the patch recovers
+0.8-0.97 of the continuous fault's transmission across the joint; the
+slip crosses on the cut itself (the segment's pair jump reaches the
+continuous fault's); the rest of the network keeps the split's answer
+(main strand within 2.5%); the weak patch reproduces a fully plastic
+patch on the same cells to 1-2% and is insensitive to the ratio from
+0.01 to 0.001; the solve is linear and costs the split's velocity
+iterations, with only the pressure block noticing the contrast
+(hence 0.01, not smaller). Gluing a joint does change the partition
+between the strands that meet there — a reconnected main line takes
+back slip a through-going branch was carrying past the weld — which
+is the junction working, not the patch leaking.
+
+`damage_yield` is the older glue: a viscoplastic plug of radius
+`max(2.5 h, 1.2 pull)` at each *prepared* junction point, yield
+`dial * (1 + 2 * edot_II)` inside, strength and rate-regularisation on
+one dial, sharp `Piecewise` boundaries. It stays available for
+studies of the glue itself, and it does not see stepover bridges,
+which are not prepared junctions.
 
 **No prescribed reconnection.** Nothing tells the network how to link
 up: the stress lobes of the abutting tips decide. A collinear gap
@@ -128,7 +264,9 @@ through redistribution — single faults are parallel-validated).
 
 ## Limitations
 
-- 3-D: planar convex patches, X crossings only, serial (above).
+- 3-D: planar convex patches, X crossings only, serial (above); the
+  weak-plane realisation is 2-D for now — place 3-D zones with
+  `place_thin_volume` directly.
 - One damage dial per network in `damage_yield` (per-junction values:
   build the expression with `uw.meshing.damage_zone_yield` directly).
 - Time-dependent damage (wear-in/healing) is study-level for now: see

docs/developer/design/fault-zone-hybrid-architecture.md

@@ -196,15 +196,35 @@ net.apply_contact(stokes)                          # only on the sliced pieces
 finite-width model. The two share one mesh, which is what makes the comparison
 between them clean.
 
+### What landed (2026-08)
+
+The whole-network end of this arrived first, in a slightly different
+shape. `build(width=..., realisation="split"|"ti")` places one ribbon
+band along every prepared piece and either cuts it or leaves it whole,
+so the two representations share one mesh as intended; `apply(solver,
+...)` imposes whichever was built, and `ti_fields` paints the weak-plane
+viscosity and director. The realisation is a property of the whole
+network, not yet of an individual piece — there is no per-fault
+`slice="auto"` criterion, so a model that is sliced *here* and finite-width
+*there* still has to be assembled by hand.
+
+The director question moved rather than closed. Within one strand's
+footprint the director is the unit normal of the nearest **segment** of
+that strand's trace, so it no longer moves when a trace is re-sampled.
+Which strand *owns* a cell is still nearest-sample over the concatenated
+spines, so the partition boundary and the high-angle-junction objection
+above are unchanged.
+
 ## What still has to be built
 
 In dependency order, for 2-D:
 
 1. Zone-boundary distance exposed so the criterion can be evaluated.
-2. The nearest-fault director, needed by the TI rheology and currently
-   hand-rolled in every script (issue #540 — broken three ways in 2-D — and
-   issue #544 for the ownership question above).
-3. The criterion itself, and the `slice="auto"` wiring.
+2. ~~The nearest-fault director~~ — landed as `FaultNetwork.ti_fields`
+   (nearest segment within a footprint); the ownership question (issue
+   #544) is untouched.
+3. The criterion itself, and the `slice="auto"` wiring — i.e. a
+   per-piece rather than per-network realisation.
 
 Placing a contact tip inside a zone works: `add_fault` puts a vertex on every
 control point of the trace, so a truncated trace terminates cleanly. That path

src/underworld3/discretisation/discretisation_mesh.py

@@ -8137,7 +8137,10 @@ def add_fault(self, faults, verbose=False):
         fault is one open polyline with both tips strictly inside the
         domain; segments must not share vertices, so branches and
         crossings are represented as OFFSET segments (a one-to-two-cell
-        ligament). The result is standalone — no geometric-MG tail, since
+        ligament). The result inherits a MESH-OWNED geometric-MG tail (the
+        parent's coarse levels, the cut mesh finest — a cut is the same
+        grid re-represented); a parent without one yields a standalone
+        mesh — no tail, since
         the coarse levels do not carry the fault (see
         :meth:`add_conforming_surface`); solvers take their
         algebraic-multigrid defaults.
@@ -8146,7 +8149,23 @@ def add_fault(self, faults, verbose=False):
         and ``docs/developer/design/FAULT_CONTACT_DEPLOYMENT_2026-08.md``.
         """
         from underworld3.utilities.fault_split import add_fault
-        return add_fault(self, faults, verbose=verbose)
+        child = add_fault(self, faults, verbose=verbose)
+        # The split mesh INHERITS a mesh-owned geometric-MG tail: a cut
+        # re-represents the same grid with the surface conformed (finer only
+        # by the duplicated vertices), so the parent's coarse levels serve
+        # unchanged with the cut mesh as the finest level — the coarse
+        # levels do not need the fault (#620/#629). Without this every
+        # solver on a split mesh fell to GAMG unless it called
+        # set_custom_fmg by hand. The FAC zone is NOT inherited: a split
+        # fault needs no patch (the keying ruling).
+        own_tail = getattr(self, "_custom_mg_coarse_meshes", None)
+        if (own_tail is not None
+                and getattr(child, "_custom_mg_coarse_meshes", None) is None):
+            child._custom_mg_coarse_meshes = list(own_tail)
+            child._custom_mg_builder = getattr(self, "_custom_mg_builder",
+                                               "barycentric")
+            child._custom_mg_fac_zone = None
+        return child
 
 
     def adapt(self, metric_field, max_levels=None, node_budget=None,