diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md index 04ea6a27..297b3df8 100644 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ b/.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 | diff --git a/docs/advanced/fault-networks.md b/docs/advanced/fault-networks.md index a4fe0b12..9d33a4c7 100644 --- a/docs/advanced/fault-networks.md +++ b/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 diff --git a/docs/developer/design/fault-zone-hybrid-architecture.md b/docs/developer/design/fault-zone-hybrid-architecture.md index cd92d84f..ee59b9fd 100644 --- a/docs/developer/design/fault-zone-hybrid-architecture.md +++ b/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 diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 769dd6fa..d3c67611 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/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, diff --git a/src/underworld3/meshing/fault_network.py b/src/underworld3/meshing/fault_network.py index d3d888f9..27151b05 100644 --- a/src/underworld3/meshing/fault_network.py +++ b/src/underworld3/meshing/fault_network.py @@ -1,17 +1,36 @@ """The user-facing 2-D fault-network toolkit. +A fault is specified ONCE — a trace, its place in the hierarchy, and +the properties it carries — and then realised. The realisation is a +keyword, not a different subsystem: the same specification, the same +prepared pieces and the same meshed band become either a cut with +node-pair contact (``realisation="split"``) or a volumetric weak plane +(``realisation="ti"``). The band carries the fault's own points and +segments as mesh vertices and edges, so it can be cut whatever its +width — the choice of realisation is not constrained by the mesh. +What differs is what the width MEANS: for the split it is a resolution +parameter that gives the cut its vertices, while for the weak plane it +is constitutive (``V = 2 e_nt w``) and wants two or three elements +across it. + One object carries the validated network recipe end to end: 1. **Hierarchy-respecting junction preparation** — where faults cross or abut, the junior trace is severed and pulled back a short ligament; the senior runs through (:func:`prepare_fault_network`). 2. **Network-refined meshing** — a graded mesh following every trace, - then split-node faults cut along the prepared pieces. -3. **Contact** — the no-opening pair constraint on every piece. + then one ribbon band placed along every prepared piece, cut or left + whole according to the realisation. +3. **Imposition** — :meth:`FaultNetwork.apply` gives the solver the + no-opening pair constraint, or the weak-plane rheology. 4. **Damage-zone glue** — small viscoplastic plugs at the junctions connect the network mechanically; the stress lobes of the abutting tips decide how slip transfers, no reconnection geometry is ever prescribed. +5. **Fault-attached properties** — :meth:`FaultNetwork.surface` returns + the retained :class:`~underworld3.meshing.surfaces.Surface` for a + piece; friction, accumulated slip and damage live there, on the + fault, and outlive any one realisation of it. The dials encode the measured rulings (2026-08): junction gaps of one or two elements transmit slip within a few percent of a continuous @@ -29,7 +48,7 @@ >>> net = uw.meshing.FaultNetwork( ... [("Main", main_pts), ("Splay", splay_pts)], ... hierarchy=["Main", "Splay"]) # Main severs Splay ->>> mesh = net.prepare(h=0.006).build() +>>> mesh = net.prepare(h=0.006).build(width=0.01) # realisation="split" >>> stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) >>> stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel >>> stokes.constitutive_model.yield_mode = "min" @@ -37,11 +56,18 @@ >>> stokes.constitutive_model.Parameters.yield_stress = \\ ... net.damage_yield(v, dial=0.05) >>> stokes.consistent_jacobian = True ->>> net.apply_contact(stokes) +>>> net.apply(stokes) >>> # ... boundary conditions ... >>> info = net.solve(stokes) >>> net.slips(stokes) {'Main': 0.14, 'Splay_1': 0.05, 'Splay_2': 0.04} + +The same specification as a weak plane — one keyword, one more +constitutive number, everything else unchanged:: + +>>> mesh = net.prepare(h=0.006).build(width=0.01, realisation="ti") +>>> net.apply(stokes, eta_1=0.01) +>>> net.slips(stokes) # the layer's own throughput """ import numpy as np @@ -51,6 +77,41 @@ from .faults import FaultSurface +def _nearest_segment_normals(P, X): + """Unit normal of the polyline segment nearest to each point of ``X``. + + The director of a weak plane, cell by cell: a curved trace has no one + orientation, and the nearest SEGMENT is the piece of fault the cell + actually lies against. 2-D; ``P`` is ``(n, 2)`` and ``X`` ``(m, 2)``. + """ + P = np.asarray(P, dtype=float)[:, :2] + X = np.asarray(X, dtype=float)[:, :2] + A, D = P[:-1], np.diff(P, axis=0) + L2 = np.maximum(np.einsum("sj,sj->s", D, D), 1e-300) + W = X[:, None, :] - A[None, :, :] + t = np.clip(np.einsum("psj,sj->ps", W, D) / L2[None, :], 0.0, 1.0) + R = W - t[:, :, None] * D[None, :, :] + k = np.argmin(np.einsum("psj,psj->ps", R, R), axis=1) + T = D[k] / np.sqrt(L2[k])[:, None] + return np.column_stack([-T[:, 1], T[:, 0]]) + + +def _densify_polyline(E, piece, per_segment=4): + """Points along an extended spine, ``per_segment`` per edge, each + flagged as lying on a CUT — an edge whose two vertices belong to the + same cut piece (``piece`` is the piece index per vertex of ``E``; + -1 for margin and gap vertices).""" + E = np.asarray(E, dtype=float) + piece = np.asarray(piece) + f = np.linspace(0.0, 1.0, per_segment, endpoint=False) + Q = (E[:-1, None, :] + f[None, :, None] * np.diff(E, axis=0)[:, None, :]) + Q = np.vstack([Q.reshape(-1, E.shape[1]), E[-1:]]) + edge_on_cut = (piece[:-1] == piece[1:]) & (piece[:-1] >= 0) + on_cut = np.concatenate([np.repeat(edge_on_cut, per_segment), + [False]]) + return Q, on_cut + + class FaultNetwork: """A hierarchy of fault traces (2-D) or planar patches (3-D) meshed, split, and glued. @@ -99,6 +160,17 @@ def __init__(self, faults, hierarchy=None): self.junctions = None self.report = None self.mesh = None + # the realisation state, set by build() + self.realisation = None + self.width = None + self.info = None + self.fault_surfaces = {} + self.ti = None + self.glue = None + self.margin_rings = None + self.spines = None # (name, polyline, piece index per vertex) + self._eta0_var = None + self._band_yield_var = None # ------------------------------------------------------------------ def prepare(self, h, ligament=2.0, through=None, verbose=True): @@ -132,18 +204,82 @@ def prepare(self, h, ligament=2.0, through=None, verbose=True): # ------------------------------------------------------------------ def build(self, base=None, h_far=None, band=0.03, ramp=0.08, - max_levels=2, qdegree=2, mesher="embed"): + max_levels=2, qdegree=2, mesher=None, + width=None, realisation="split", + margin_rings=2, carve_clearance=0.3): """Mesh the network: graded refinement along every RAW trace, - then split-node faults along the PREPARED pieces. + then the chosen REALISATION of the faults on that mesh. ``base`` is an existing coarse mesh to adapt (default: a unit ``UnstructuredSimplexBox`` at ``h_far = 4 h``); the refinement holds ``h`` within ``band`` of any trace and grades to ``h_far`` over ``ramp``. + + ``width`` is the fault band's thickness. Give it and the network + is placed as a ribbon band (2-D) that both realisations share: + the SAME mesh is cut and split (``realisation="split"``) or left + whole for a volumetric weak-plane rheology + (``realisation="ti"``), which is what makes the two comparable. + What ``width`` MEANS differs: for the split it is a resolution + parameter — the band exists to give the cut its own vertices — + while for TI it is constitutive, the layer thickness that sets + the slip rate ``V = 2 e_nt w``, so it wants two or three elements + across it. The width does NOT decide whether the fault can be + cut: the band is meshed around the trace's own points and + segments, so the cut chain is part of the mesh by construction at + any width (measured to ``w = h_far / 10``; see + ``~/+Simulations/fault_split_at_ti_width``). + + ``carve_clearance`` sizes the cells cleared around the ribbons + before they are meshed in. On a graded base (``max_levels > 0`` + refining toward the traces) the default can fail to leave one + simple cavity at fine widths; raise it (1.0 worked where 0.3 did + not) or build on the uniform base with ``max_levels=0``. + + ``width=None`` keeps the original no-band path: graded refinement + cut directly. It is split-only, and its mesh is NOT the one a TI + run would use, so do not compare the two across that choice. + + ``mesher`` picks how the fault is meshed into the mesh, from + the choices that dimension offers. In 2-D: ``"network"`` + (default) places every strand in one fused call, so strands may + touch, and ``"ladder"`` places them sequentially and lets placed + levels nest (see + :func:`~underworld3.utilities.place_surface.place_fault_ribbon_2d`). + In 3-D: ``"embed"`` (default) and ``"place"``, described in + :meth:`_build_3d`. """ if self.prepared is None: raise RuntimeError("call prepare(h=...) first") + meshers = {2: ("network", "ladder"), 3: ("embed", "place")}[self.dim] + if mesher is None: + mesher = meshers[0] + if mesher not in meshers: + raise ValueError( + f"mesher must be one of {meshers} in {self.dim}-D, not " + f"{mesher!r}") + if realisation not in ("split", "ti"): + raise ValueError( + f"realisation must be 'split' or 'ti', not {realisation!r}") + if realisation == "ti" and width is None: + raise ValueError( + "realisation='ti' needs width=: the weak plane is a LAYER, " + "and its thickness is constitutive (V = 2 e_nt w). Pass the " + "band width you intend to resolve.") + self.realisation = realisation + self.width = None if width is None else float(width) if self.dim == 3: + if realisation != "split": + raise NotImplementedError( + "the 3-D network builds the split realisation only; " + "place the patches with place_thin_volume for a " + "volumetric zone.") + if width is not None: + raise NotImplementedError( + "the 3-D network does not place a band: its patches " + "are meshed conforming (mesher='embed') or placed as " + "sheets (mesher='place'), both of zero thickness. For " + "a finite-width 3-D zone call place_thin_volume.") return self._build_3d(h_far=h_far, qdegree=qdegree, mesher=mesher) from .cartesian import UnstructuredSimplexBox @@ -168,10 +304,127 @@ def metric(pts_, _ss=surfs, _h=h, _hf=h_far, _b=band, _r=ramp): return 1.0 / hh ** 2 child = base.adapt(metric, max_levels=max_levels) - self.mesh = child.add_fault( - [(n, p.copy()) for n, p in self.prepared]) + if width is None: + self.mesh = child.add_fault( + [(n, p.copy()) for n, p in self.prepared]) + self.info = None + else: + from underworld3.utilities.place_surface import ( + place_fault_ribbon_2d) + # The ribbons must MEET across every junction ligament: the + # ribbon is the fault as the weak plane sees it, continuous + # through a junction where the cut stops short, and that + # difference is what junction_cells() reads. So the tip + # margin reaches at least as far as the longest pull-back. + # Collinear abutting pieces share ONE spine (two ribbons + # overlapping along a line interleave their vertices into + # slivers); the cut still stops at each piece's own ends, and + # the gap between them is spine the split does not cut. + self.spines = self._shared_spines(margin_rings) + margin_rings = self._junction_margins(int(margin_rings)) + self.margin_rings = margin_rings + self.mesh, self.info = place_fault_ribbon_2d( + child, [(n, p.copy()) for n, p in self.prepared], + self.width, margin_rings=margin_rings, + clearance=carve_clearance, + split=(realisation == "split"), mesher=mesher, + spines=[(n, S) for n, S, _idx in self.spines]) + self._make_surfaces() return self.mesh + def _shared_spines(self, rings): + """Group the prepared pieces into spines: a piece whose end faces + the start of another within the margins' reach, along the same + line (end tangents within ~25 degrees, the far start within half + a width of the line), continues onto that piece's spine. Returns + ``(name, polyline, piece_index_per_vertex)`` per spine; vertices + inserted across a gap carry index -1 (spine the split does not + cut).""" + pieces = [np.asarray(P, dtype=float) for _n, P in self.prepared] + n = len(pieces) + follower = {} + for i, A in enumerate(pieces): + tA = A[-1] - A[-2] + tA /= np.linalg.norm(tA) + sA = float(np.linalg.norm(np.diff(A, axis=0), axis=1).mean()) + best = None + for j, B in enumerate(pieces): + if j == i: + continue + gap = B[0] - A[-1] + dist = float(np.linalg.norm(gap)) + sB = float(np.linalg.norm(np.diff(B, axis=0), axis=1).mean()) + if dist > rings * (sA + sB) + self.width: + continue + tB = B[1] - B[0] + tB /= np.linalg.norm(tB) + along = float(gap @ tA) + lateral = float(abs(gap[0] * tA[1] - gap[1] * tA[0])) + if (tA @ tB < np.cos(np.radians(25)) or along < 0.0 + or lateral > 0.5 * self.width): + continue + if best is None or dist < best[1]: + best = (j, dist) + if best is not None and best[0] not in follower.values(): + follower[i] = best[0] + heads = [i for i in range(n) if i not in follower.values()] + spines = [] + for h in heads: + chain = [h] + while chain[-1] in follower: + chain.append(follower[chain[-1]]) + S, idx = [pieces[chain[0]]], [np.full(len(pieces[chain[0]]), chain[0])] + for a, b in zip(chain[:-1], chain[1:]): + A, B = pieces[a], pieces[b] + s = 0.5 * (np.linalg.norm(A[-1] - A[-2]) + + np.linalg.norm(B[1] - B[0])) + gap = float(np.linalg.norm(B[0] - A[-1])) + # bridge the gap at the local rung; a gap of one rung is + # one edge between the two cut ends, which is already + # spine no cut owns (its ends belong to different pieces) + k = max(0, int(round(gap / s)) - 1) + f = np.linspace(0.0, 1.0, k + 2)[1:-1] + S.append(A[-1] + f[:, None] * (B[0] - A[-1])) + idx.append(np.full(k, -1)) + S.append(B) + idx.append(np.full(len(B), b)) + name = "+".join(self.prepared[c][0] for c in chain) + spines.append((name, np.vstack(S), np.concatenate(idx))) + return spines + + def _junction_margins(self, rings): + """Per-end margin counts per SPINE: ``rings`` at a free tip; at + an end that sits on a prepared junction, enough rings for the + ribbon to reach the OTHER piece's cut (plus half a width) — the + whole ligament lies in both ribbons, so the junction cells cover + it end to end. The margin is built by continuing the end + segment, so the count depends on that segment's length.""" + pieces = {n: np.asarray(P, dtype=float) for n, P in self.prepared} + out = [] + for _sname, S, idx in self.spines: + ends = [] + for end, seg, piece in ((S[0], S[1] - S[0], idx[0]), + (S[-1], S[-1] - S[-2], idx[-1])): + name = self.prepared[int(piece)][0] + reach = 0.0 + for j in (self.junctions or []): + if name not in j["faults"]: + continue + pull = float(j["pull"]) + near = (np.linalg.norm(np.asarray(j["point"]) - end) + <= pull + self.width) + if not near: + continue + other = [n for n in j["faults"] if n != name][0] + gap = float(np.linalg.norm(pieces[other] - end, + axis=1).min()) + reach = max(reach, gap + 0.5 * self.width) + seg_len = float(np.linalg.norm(seg)) + ends.append(max(rings, int(np.ceil(reach / seg_len)) + 1) + if reach > 0 else rings) + out.append((ends[0], ends[1])) + return out + def _build_3d(self, h_far=None, qdegree=2, mesher="embed", minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), band=None, ramp=None, max_levels=2, clearance=0.8): @@ -216,6 +469,7 @@ def _build_3d(self, h_far=None, qdegree=2, mesher="embed", for name, _p in self.prepared: mesh = split_fault(mesh, name) self.mesh = mesh + self._make_surfaces() return self.mesh if mesher != "place": raise ValueError(f"mesher must be 'place' or 'embed', " @@ -256,6 +510,7 @@ def metric(pts_, _sheets=sheets, _h=h, _hf=h_far, _b=b, _r=r): for n, _sp, _st in sheets: mesh = split_fault(mesh, n) self.mesh = mesh + self._make_surfaces() return self.mesh @staticmethod @@ -383,14 +638,355 @@ def area(p, q, r): "degenerate; drop it or refine h.") # ------------------------------------------------------------------ - def apply_contact(self, solver, conds=0): - """Register the split-node contact on every prepared piece.""" + def _make_surfaces(self): + """Retain a :class:`~underworld3.meshing.surfaces.Surface` per + prepared piece, on the BUILT mesh. + + The fault is specified once and then realised; the surface object + is what survives both realisations, and it is where properties + that belong to the FAULT rather than to the mesh live — + ``add_variable("friction")``, accumulated slip, a damage state. + The realisation reads them; it does not own them. In 3-D the + input :class:`FaultSurface` objects play the same role. + """ + if self.dim == 3: + self.fault_surfaces = {s.name: s for s in (self.surfaces or [])} + return self.fault_surfaces + from .surfaces import Surface + + self.fault_surfaces = {} + for name, P in self.prepared: + pts = np.asarray(P, dtype=float) + if pts.shape[1] == 2: + pts = np.column_stack([pts, np.zeros(len(pts))]) + self.fault_surfaces[name] = Surface(name, self.mesh, pts) + return self.fault_surfaces + + def surface(self, name): + """The retained fault surface for one prepared piece.""" + if not self.fault_surfaces: + raise RuntimeError("call build() first") + if name not in self.fault_surfaces: + raise KeyError(f"no fault surface {name!r}; the network holds " + f"{sorted(self.fault_surfaces)}") + return self.fault_surfaces[name] + + # ------------------------------------------------------------------ + def apply(self, solver, conds=0, eta_1=None, eta_0=1.0, tag="", + normal=None): + """Impose the network on ``solver``, in whichever realisation + ``build()`` made. + + ``realisation="split"`` registers the no-opening contact pair on + every prepared piece (``conds`` is the datum, 0 for free slip). + ``realisation="ti"`` paints the weak-plane fields on the honoured + footprints and hands the solver a + :class:`~underworld3.constitutive_models.TransverseIsotropicFlowModel`: + ``eta_1`` (required) is the weak-plane viscosity, ``eta_0`` the + background — a float, or a per-cell array if the background is + itself painted (a terrane, say). ``tag`` disambiguates the field + names when more than one network is applied to one mesh. + + ``normal`` is the split's fault normal (see + :meth:`~underworld3.systems.Stokes.add_fault_bc`). Pass + ``"trace"`` when the traces are SAMPLED SMOOTH CURVES — the + default per-node normal zig-zags at the sampling kinks, and + the no-opening constraint then notches the slip. + """ + if self.mesh is None: + raise RuntimeError("call build() first") + if self.realisation == "split": + for name, _p in self.prepared: + solver.add_fault_bc(conds, boundary=name, normal=normal) + return self + if eta_1 is None: + raise ValueError( + "realisation='ti' needs eta_1=: the weak-plane viscosity " + "is the other half of the constitutive pair (with width).") + import underworld3 as uw + + eta1, ndir, foot = self.ti_fields(eta_1, eta_0=eta_0, tag=tag) + solver.constitutive_model = \ + uw.constitutive_models.TransverseIsotropicFlowModel + params = solver.constitutive_model.Parameters + params.shear_viscosity_0 = ( + float(eta_0) if np.ndim(eta_0) == 0 else self._eta0_var.sym[0]) + params.shear_viscosity_1 = eta1.sym[0] + params.director = ndir.sym + return self + + def apply_contact(self, solver, conds=0, normal=None): + """Register the split-node contact on every prepared piece. + + The split realisation's half of :meth:`apply`, kept under its own + name for callers that only ever want the contact. + """ + if self.realisation not in (None, "split"): + raise RuntimeError( + f"this network was built as {self.realisation!r}; there " + f"are no fault pairs to constrain. Use apply().") if self.mesh is None: raise RuntimeError("call build() first") for name, _p in self.prepared: - solver.add_fault_bc(conds, boundary=name) + solver.add_fault_bc(conds, boundary=name, normal=normal) return self + # ------------------------------------------------------------------ + def ti_fields(self, eta_1, eta_0=1.0, tag=""): + """The weak-plane (TI) realisation's painted P0 fields. + + ``eta_1`` inside each fault's HONOURED footprint — the band cells + whose nearest spine sample is a USER point, never the whole band, + whose margin is extrapolated surround — and the background + elsewhere; the director is the unit normal of the nearest segment + of the strand that owns the cell, so a curved trace carries its + own orientation cell by cell. + + Returns ``(eta_1_var, director_var, footprint_mask)``. The mask is + also the right ``fac_zone`` key for a multigrid patch (#629). + """ + if self.info is None: + raise RuntimeError( + "no band on this mesh: build(width=...) first (the weak " + "plane is a layer, and the layer has to be meshed).") + import underworld3 as uw + + foots = self.info["footprints"] + foot = np.zeros_like(next(iter(foots.values()))) + for m_ in foots.values(): + foot = foot | m_ + + eta1 = uw.discretisation.MeshVariable( + f"fnEta1{tag}", self.mesh, 1, degree=0) + eta0_vals = np.broadcast_to(np.asarray(eta_0, dtype=float), + (len(eta1.coords),)) + eta1.array[:, 0, 0] = np.where(foot, float(eta_1), eta0_vals) + self._eta0_var = None + if np.ndim(eta_0) != 0: + self._eta0_var = uw.discretisation.MeshVariable( + f"fnEta0{tag}", self.mesh, 1, degree=0) + self._eta0_var.array[:, 0, 0] = eta0_vals + + dim = self.mesh.dim + ndir = uw.discretisation.MeshVariable( + f"fnDir{tag}", self.mesh, dim, degree=0, continuous=False) + cen = np.asarray(ndir.coords)[:, :dim] + dvals = np.zeros((len(cen), dim)) + dvals[:, -1] = 1.0 # any unit vector outside the + for name, P in self.prepared: # footprints: eta_1 == eta_0 + m_ = foots[name] # there, so TI is isotropic + if not m_.any(): + continue + dvals[m_] = _nearest_segment_normals(P, cen[m_]) + ndir.array[...] = dvals.reshape(ndir.array.shape) + self.ti = {"eta_1": eta1, "director": ndir, "footprint": foot} + return eta1, ndir, foot + + # ------------------------------------------------------------------ + @property + def band(self): + """The band's cell mask — the material the fault is embedded in. + + The band is not scaffolding for the weak plane. It is a meshed + region of material AROUND the fault, and the split wants it as + much: a segmented fault does its interesting work at the tips and + in the ligaments between strands, and damage there needs cells to + live in. This mask (and :attr:`footprints`, per strand) is how a + rheology addresses that region in either realisation. + """ + if self.info is None: + raise RuntimeError( + "no band on this mesh: build(width=...) first") + return self.info["band"] + + @property + def footprints(self): + """Per-strand FAULT footprints — the band cells whose nearest + spine sample is a USER point, never the extrapolated margin.""" + if self.info is None: + raise RuntimeError( + "no band on this mesh: build(width=...) first") + return self.info["footprints"] + + def band_yield(self, tau_y, tau_far=1.0e8, tag=""): + """Von Mises yield confined to the band, as an expression. + + The band's cells yield at ``tau_y``; everything else is given + ``tau_far``, high enough never to yield. Pair it with a + :class:`~underworld3.constitutive_models.ViscoPlasticFlowModel` + and ``consistent_jacobian = True``:: + + stokes.constitutive_model.Parameters.yield_stress = \ + net.band_yield(tau_y=4.0) + + This is the damage the SPLIT realisation wants. A released fault + flank sits far below ``tau_y`` and is untouched; the places that + sit far above it — a strand's tips, the weld where a cut stops + short, the sliver at a junction — yield by themselves, so the + breakdown appears where the mechanics puts it rather than where a + geometric plug was placed. Compare :meth:`damage_yield`, which + places plugs at the junctions by construction and is the right + tool when the junction glue itself is the object of study. + + The region is a sharp mask, not a blend: never taper a + rheological parameter towards a large sentinel. + """ + if self.info is None: + raise RuntimeError( + "no band on this mesh: build(width=...) first (the " + "damage needs cells to live in)") + import underworld3 as uw + + ybar = uw.discretisation.MeshVariable( + f"fnTauY{tag}", self.mesh, 1, degree=0) + ybar.array[:, 0, 0] = np.where(self.info["band"], float(tau_y), + float(tau_far)) + self._band_yield_var = ybar + return ybar.sym[0] + + # ------------------------------------------------------------------ + def junction_cells(self, ring=1): + """The cells where the split cannot join: the ribbon minus the cut. + + The ribbon (the band, with its extrapolated tip margins) is + everything the weak-plane realisation treats as fault; the cut + chains are what the split actually sliced. A band cell whose + nearest point on a piece's extended spine lies in that piece's + margin (or on the gap of a shared spine) is fault the split did + not cut. Where such a cell also lies inside a SECOND piece's + ribbon, two pieces meet there and the split has left them welded + — a stop-short abutment, a kissing branch, the intact bridge of + a stepover. Those cells, dilated by ``ring`` vertex rings within + the band, are the junction cells. + + The rule is geometric and comes entirely from the placement: no + stress threshold, and the free tips are excluded on purpose (a + margin that runs into intact material rather than another + ribbon is a tip, not a junction; damage placed there lengthens + the fault instead of joining it). + + ``ring=1`` is not optional in practice. Measured on the S-fault + rig (2026-08-27): the bare junction cells cannot repair the weld + even when fully plastic (a fifth to a quarter of the deficit), + because the weld's stiffness lives in the ring of intact material + around the two tips; one ring restores 0.8-0.97 of a continuous + fault's transmission at both resolutions tested. + + Returns a boolean cell mask over the mesh's cells. + """ + if self.info is None: + raise RuntimeError( + "no band on this mesh: build(width=...) first (the " + "junction cells are read off the ribbon)") + if self.realisation != "split": + raise RuntimeError( + "junction cells are the SPLIT realisation's joints; the " + "weak plane has no cut to fall off") + from underworld3.utilities.place_surface import _cell_centroids_of + + band = self.info["band"] + ids, cen = _cell_centroids_of(self.mesh.dm, band) + mask = np.zeros_like(band) + if len(ids) == 0: + return mask + half_width = 0.5 * self.width + from underworld3.utilities.place_surface import _extend_polyline_2d + + def nearest(Q): + d = cen[:, None, :] - Q[None, :, :] + d2 = np.einsum("ijk,ijk->ij", d, d) + j = np.argmin(d2, axis=1) + return j, np.sqrt(d2[np.arange(len(ids)), j]) + + # off the cut, read spine by spine: a cell in a spine's ribbon + # whose nearest spine point is a gap edge or a tip margin is fault + # that spine's split did not cut. (It may well sit on ANOTHER + # spine's cut — the senior's flank beside a kissing tip is exactly + # where the glue belongs.) + off_cut = np.zeros(len(ids), dtype=bool) + for k, (_sname, _S, idx) in enumerate(self.spines): + E = np.asarray(self.info["extended"][k], dtype=float) + m0, m1 = self.info["margin_rings"][k] + piece = np.concatenate([np.full(m0, -1), idx, np.full(m1, -1)]) + Q, on_cut = _densify_polyline(E, piece) + j, dist = nearest(Q) + # a sample spacing's worth of tolerance: the densified spine + # is a polyline, the cell centroid a point beside it + inside = dist <= half_width + 0.35 * float(self.info["spacing"][k]) + off_cut |= inside & ~on_cut[j] + # two pieces meet: the cell lies in the ribbons of two CUT pieces, + # each continued by the default margin + ribbons = np.zeros(len(ids), dtype=int) + for name, P in self.prepared: + P = np.asarray(P, dtype=float) + E = _extend_polyline_2d(P, 2) + spacing = float(np.linalg.norm(np.diff(P, axis=0), axis=1).mean()) + _j, dist = nearest(E) + ribbons += dist <= half_width + 0.5 * spacing + # unjoined = off some cut, where two pieces' ribbons meet + mask[ids[off_cut & (ribbons >= 2)]] = True + for _ in range(int(ring)): + mask = self._vertex_ring(mask) & band + return mask + + def _vertex_ring(self, mask): + """Cells sharing a vertex with a masked cell (rank-local).""" + dm = self.mesh.dm + cS, cE = dm.getHeightStratum(0) + vS, vE = dm.getDepthStratum(0) + out = mask.copy() + for c in np.flatnonzero(mask): + for q in dm.getTransitiveClosure(int(c) + cS)[0]: + if vS <= int(q) < vE: + for c2 in dm.getTransitiveClosure(int(q), useCone=False)[0]: + if cS <= int(c2) < cE: + out[int(c2) - cS] = True + return out + + def junction_patch(self, eta_0=1.0, ratio=0.01, ring=1, tag=""): + """The junction glue: a weak isotropic patch on the junction cells. + + Returns the viscosity to give an isotropic flow model, as a P0 + field expression: ``eta_0`` everywhere, ``ratio * eta_0`` on + :meth:`junction_cells`. ``eta_0`` is a float or a per-cell array + (a painted background). Use it as the split realisation's + background viscosity:: + + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = \ + net.junction_patch(eta_0=1.0) + net.apply(stokes) + + The glue is a viscosity RATIO rather than a yield stress because + the joint only has to be broken, and a ratio needs no stress + scale. Measured on the S-fault rig (2026-08-27, coarse and fine): + the weak patch reproduces a fully plastic patch on the same cells + to 1-2%, is insensitive to ``ratio`` from 0.01 to 0.001, carries + the segment's slip through the cut itself (the segment's pair + jump reaches a continuous fault's), leaves the rest of the + network the split's own answer (main strand within 2.5%), and + costs the split's velocity iterations. The pressure block alone + notices the contrast, which is why 0.01 is the default rather + than something smaller. The solve stays linear. + + Isotropic on purpose: inside a junction there is no single plane + to be weak on (see :meth:`damage_yield`, whose plugs are placed + at the prepared junction POINTS by radius; this patch is read off + the mesh instead and also catches stepover bridges, which are not + prepared junctions). + """ + import underworld3 as uw + + cells = self.junction_cells(ring=ring) + eta = uw.discretisation.MeshVariable( + f"fnGlue{tag}", self.mesh, 1, degree=0) + eta0_vals = np.broadcast_to(np.asarray(eta_0, dtype=float), + (len(eta.coords),)) + eta.array[:, 0, 0] = np.where(cells, float(ratio) * eta0_vals, + eta0_vals) + self.glue = {"cells": cells, "viscosity": eta, "ratio": float(ratio)} + return eta.sym[0] + # ------------------------------------------------------------------ def damage_yield(self, velocity, dial=0.05, radius=None, tau_far=1.0e3): @@ -455,7 +1051,21 @@ def solve(self, solver, **kwargs): # ------------------------------------------------------------------ def slips(self, solver): - """Peak tangential slip per prepared piece (rank-local).""" + """Peak tangential slip per prepared piece, in each realisation's + OWN quantity (rank-local). + + The split's slip is the tangential jump between the two nodes of + a cut pair. The weak plane has no pair: its slip is the jump in + tangential velocity across the layer, sampled one band half-width + plus a cell either side of the spine — read from the velocity + field itself rather than integrated from the in-band strain rate, + which is vertex-phase sensitive once ``w`` approaches ``h``. + Both are the layer's own throughput, so the two numbers may be + compared; a probe placed further out than this reads the + surrounding flow as well and over-reads short strands. + """ + if self.realisation == "ti": + return self._slips_ti(solver) from underworld3.utilities.fault_contact import fault_pair_jumps info = getattr(solver, "_rotated_freeslip_info", None) if info is None: @@ -472,12 +1082,38 @@ def slips(self, solver): out[name] = float(np.linalg.norm(tangential, axis=1).max()) return out + def _slips_ti(self, solver): + """The weak plane's slip: the tangential velocity jump across the + band, one half-width plus a cell either side of each spine.""" + import underworld3 as uw + + if self.info is None: + raise RuntimeError("no band on this mesh: build(width=...)") + out = {} + for k, (name, P) in enumerate(self.prepared): + P = np.asarray(P, dtype=float)[:, :2] + t = np.gradient(P, axis=0) + t /= np.linalg.norm(t, axis=1)[:, None] + n = np.column_stack([-t[:, 1], t[:, 0]]) + skirt = 0.5 * self.width + float(self.info["spacing"][k]) + vp = np.asarray(uw.function.evaluate( + solver.u.sym, P + skirt * n)).reshape(len(P), -1)[:, :2] + vm = np.asarray(uw.function.evaluate( + solver.u.sym, P - skirt * n)).reshape(len(P), -1)[:, :2] + out[name] = float( + np.abs(np.einsum("ij,ij->i", vp - vm, t)).max()) + return out + # ------------------------------------------------------------------ def __repr__(self): n_j = len(self.junctions) if self.junctions is not None else "?" n_p = len(self.prepared) if self.prepared is not None else "?" state = ("meshed" if self.mesh is not None else "prepared" if self.prepared is not None else "raw") + if self.mesh is not None: + state += f" as {self.realisation}" + if self.width is not None: + state += f", w={self.width:g}" return (f"FaultNetwork({len(self.faults)} faults -> {n_p} " f"pieces, {n_j} junctions, {state}; " f"hierarchy={self.hierarchy})") diff --git a/src/underworld3/meshing/surfaces.py b/src/underworld3/meshing/surfaces.py index 690bdfc5..d62027e2 100644 --- a/src/underworld3/meshing/surfaces.py +++ b/src/underworld3/meshing/surfaces.py @@ -2984,15 +2984,25 @@ def dist_to(P, Q): for j, (nj, pj, _cj) in enumerate(traces): if j == i: continue - if dist_to(P, pj) >= lig: + d0 = dist_to(P, pj) + if d0 >= lig: continue if any(abs(arc - arc_end) < 4.0 * lig for arc, _ in ci): continue # already handled above - pull = lig - for _ in range(5): + # Pull back by the clearance DEFICIT, not a whole ligament: + # the join should be as small as the mesh allows. When the + # other trace's END is the near part (two ends facing each + # other), each side yields half — the other end is pulled + # by its own pass below. + facing = (dist_to(pj[0], pi) < lig + or dist_to(pj[-1], pi) < lig) + share = 0.5 if facing else 1.0 + target = d0 + share * (lig - d0) + pull = share * (lig - d0) + for _ in range(6): s_q = arc_end + pull if arc_end < 1e-12 \ else arc_end - pull - if dist_to(_point_at_arc(pi, s_q), pj) >= lig: + if dist_to(_point_at_arc(pi, s_q), pj) >= target: break pull *= 1.6 ci.append((arc_end, pull)) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 872194c4..3c184d4a 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -1760,8 +1760,12 @@ def build_transfers(solver, field_id=None): _attempts = [builder] + (["rbf"] if builder != "rbf" else []) h = Ps = None for _i, _b in enumerate(_attempts): - h = CustomMGHierarchy(level_tail + [solver.mesh], builder=_b, - field_id=field_id) + h = CustomMGHierarchy( + level_tail + [solver.mesh], builder=_b, field_id=field_id, + # a MESH-OWNED FAC zone (set by the placement that built the + # mesh — the fault band): the finest level's strong patch + # smoother keys on it with no per-solver set_custom_fmg call + fac_zone=getattr(solver.mesh, "_custom_mg_fac_zone", None)) try: Ps = h.build(solver) break @@ -1951,3 +1955,33 @@ def inject_custom_mg(solver): levels.append(fine) Ps = [_to_petsc_aij(builder(levels[l - 1], levels[l])) for l in range(1, len(levels))] _install_transfers(solver, Ps, verbose=cfg.get("verbose", False)) + + +def adopt_hierarchy(mesh, base_mesh, fac_zone=None, builder=None): + """Make ``mesh`` OWN the multigrid hierarchy of ``base_mesh`` — the + static coarse tail every solver built on ``mesh`` then drives + automatically (standard and rotated paths alike, through + :func:`build_transfers`'s mesh-owned route), with ``fac_zone`` as the + finest level's FAC patch key. For a mesh produced by SURGERY on a + refined base (a placed fault band, a network of ribbons): the coarse + levels do not need the fault — the finest level inherits the tail + (#620/#629). Without this, a fresh Mesh from a surgery DM owns no + hierarchy and every solver on it silently falls back to GAMG (whose + Chebyshev smoother is the configuration a nonlinear solve should not + be handed by default). A later ``mesh.add_fault`` child inherits the + tail itself; the FAC zone is NOT inherited by the cut child (a split + fault needs no patch — the keying ruling). + """ + # a base that is itself an adapt()/cut child owns its tail: extend it + # with the base (the child's finest level) — the same rule + # Mesh._adopt_cut_child applies; a plain refined base contributes its + # static level wraps (coarsest .. base-finest) + own = getattr(base_mesh, "_custom_mg_coarse_meshes", None) + mesh._custom_mg_coarse_meshes = (list(own) + [base_mesh] if own is not None + else list(base_mesh._coarse_level_meshes())) + mesh._custom_mg_builder = (builder if builder is not None + else getattr(base_mesh, "_custom_mg_builder", + "barycentric")) + mesh._custom_mg_fac_zone = (None if fac_zone is None + else np.asarray(fac_zone, dtype=bool)) + return mesh diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index c656e786..6dcfb949 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -531,6 +531,12 @@ def discrete_loop(loop_verts): gmsh.option.setNumber("Mesh.MeshSizeMin", 0.5 * float(lengths.min())) gmsh.option.setNumber("Mesh.MeshSizeMax", 2.0 * float(lengths.max())) if size_of is not None: + # the callback is AUTHORITATIVE for the interior: without these + # gmsh extends the constrained curves' segmentation inward and + # takes the minimum, and a graded callback has no effect + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) gmsh.model.mesh.setSizeCallback( lambda dim, tag, x, y, z, lc: float(size_of(x, y))) gmsh.model.mesh.generate(2) @@ -5008,10 +5014,15 @@ def _footprint_from_samples(dm, band_mask, samples_ext, is_user_sample): def _extend_polyline_2d(P, rings): """Continue a polyline ``rings`` points outward at both ends, linearly — the 2-D tip-margin builder (:func:`_extend_grid` one dimension - down): end tangents at the local spacing, no invented curvature.""" + down): end tangents at the local spacing, no invented curvature. + ``rings`` is one count for both ends or a ``(start, end)`` pair (a + junction end wants a longer reach than a free tip).""" P = np.asarray(P, dtype=float) - for _ in range(int(rings)): - P = np.vstack([2.0 * P[0] - P[1], P, 2.0 * P[-1] - P[-2]]) + r0, r1 = (rings, rings) if np.ndim(rings) == 0 else rings + for _ in range(int(r0)): + P = np.vstack([2.0 * P[0] - P[1], P]) + for _ in range(int(r1)): + P = np.vstack([P, 2.0 * P[-1] - P[-2]]) return P @@ -5232,7 +5243,8 @@ def _occ_ladder_assembly_2d(polylines, width, size, assembly="fuse", gmsh.finalize() -def _occ_assembly_2d(polylines, width, size, assembly="fuse", domain=None): +def _occ_assembly_2d(polylines, width, size, assembly="fuse", domain=None, + embed=None): """Thicken each polyline into a ribbon, resolve overlaps, mesh. The 2-D thin volume: a ribbon is the mitred outline of one polyline, and @@ -5251,6 +5263,16 @@ def _occ_assembly_2d(polylines, width, size, assembly="fuse", domain=None): on the boundary's own facets (snapped onto them after meshing, defensively). ``cad_area`` is then the clipped area, so the meshed-vs-CAD gate holds unchanged. + + ``embed`` — polylines (each ``(n, 2)``, strictly inside the resolved + faces) whose points and segments are EMBEDDED in the face before + meshing, so they are vertices and edges of the mesh exactly. This is + the NETWORK path (``mesher="network"``): the fuse resolves junctions + between touching ribbons as ordinary cells, and the embedded spines + give a split cut its own vertices to walk (#595: nothing snaps) — + the two properties the sequential ladder and the plain fuse each + had only one of. Each embedded point carries its local segment + length as mesh size so gmsh does not subdivide the spine. """ import gmsh @@ -5331,6 +5353,31 @@ def outline(P): faces = gmsh.model.getEntities(2) cad_area = sum(occ.getMass(2, t) for _d, t in faces) + if embed: + face_tags = [t for _d, t in faces] + for P in embed: + P = np.asarray(P, dtype=float)[:, :2] + if len(P) < 2: + continue + seg = np.linalg.norm(np.diff(P, axis=0), axis=1) + hp = np.concatenate([[seg[0]], 0.5 * (seg[1:] + seg[:-1]), + [seg[-1]]]) + ptags = [occ.addPoint(q[0], q[1], 0.0, meshSize=float(h)) + for q, h in zip(P, hp)] + ltags = [occ.addLine(ptags[i], ptags[i + 1]) + for i in range(len(ptags) - 1)] + occ.synchronize() + mid = P[len(P) // 2] + host = [t for t in face_tags + if gmsh.model.isInside(2, t, [mid[0], mid[1], 0.0])] + if not host: + raise RuntimeError( + "network assembly: an embedded spine lies outside " + "the resolved ribbon faces (trim its ends inside " + "the caps).") + gmsh.model.mesh.embed(0, ptags, 2, host[0]) + gmsh.model.mesh.embed(1, ltags, 2, host[0]) + gmsh.option.setNumber("Mesh.MeshSizeMin", 0.7 * size) gmsh.option.setNumber("Mesh.MeshSizeMax", 1.3 * size) gmsh.model.mesh.generate(2) @@ -5537,6 +5584,12 @@ def _ring_growing(cells, drop, held_mask): if not pinch: ring = _cavity_ring(cells, np.flatnonzero(drop)) if ring is None: + # TODO(BUG): on a GRADED base (adapt() toward the traces) the cleared + # cells stop forming one simple hole at the default clearance: the + # S-fault rig at fine width (w = 0.01) builds on the uniform base at + # clearance 0.3 but on the one-level graded base only at 1.0. The + # clearing should be sized by the LOCAL cell size, not one factor. + # Found 2026-08-27 running FaultNetwork.build on the rig. raise RuntimeError( "the cells cleared for the thin volume do not leave one " "simple hole. Raise `clearance`.") @@ -6196,10 +6249,21 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, payload = None if comm.rank == 0: try: - assemble = (_occ_ladder_assembly_2d if mesher == "ladder" - else _occ_assembly_2d) - asm_pts, asm_tris, cad_area = assemble( - polylines, width, size, assembly, domain=domain_loops) + if mesher == "ladder": + asm_pts, asm_tris, cad_area = _occ_ladder_assembly_2d( + polylines, width, size, assembly, domain=domain_loops) + elif mesher == "network": + # fuse + embedded spines: each polyline's interior points + # (the end points sit ON the caps) become mesh vertices; + # junctions between touching ribbons are free + spines = [np.asarray(P, dtype=float)[1:-1] + for P in polylines] + asm_pts, asm_tris, cad_area = _occ_assembly_2d( + polylines, width, size, assembly, domain=domain_loops, + embed=spines) + else: + asm_pts, asm_tris, cad_area = _occ_assembly_2d( + polylines, width, size, assembly, domain=domain_loops) P = asm_pts[asm_tris] twice = ((P[:, 1, 0] - P[:, 0, 0]) * (P[:, 2, 1] - P[:, 0, 1]) - (P[:, 1, 1] - P[:, 0, 1]) * (P[:, 2, 0] - P[:, 0, 0])) @@ -6373,7 +6437,34 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, Xall = np.vstack([X, asm_pts]) holes = [[len(X) + int(v) for v in loop] for loop in hole_loops] - gap_tris, extra = _gmsh_fill_2d(Xall, ring, None, holes=holes) + # GRADED fill: the annulus between the assembly's skin and the + # cavity ring is meshed from the skin's own size out to the + # ring's edge length, interpolated by relative distance. Without + # this the fill inherits the skin segmentation throughout, and + # where several ribbons sit within a cavity of each other the + # merged cavity fills at band resolution end to end (measured on + # the S-fault rig: the fill cost as many cells as the bands — + # the #629 "fill shell was the fat" finding on the network path). + size_of = None + if hole_loops: + from scipy.spatial import cKDTree as _KDT + _skin = asm_pts[np.unique(np.concatenate( + [np.asarray(l, dtype=int) for l in hole_loops]))] + _ring_pts = Xall[np.asarray(ring, dtype=int)] + _rl = np.linalg.norm(np.diff( + np.vstack([_ring_pts, _ring_pts[:1]]), axis=0), axis=1) + _h_ring, _h_skin = float(np.median(_rl)), float(size) + if _h_ring > 1.2 * _h_skin: + _kd_s, _kd_r = _KDT(_skin), _KDT(_ring_pts) + + def size_of(x, y, _s=_h_skin, _h=_h_ring, + _ks=_kd_s, _kr=_kd_r): + q = np.array([[x, y]]) + ds = float(_ks.query(q)[0][0]) + dr = float(_kr.query(q)[0][0]) + return _s + (_h - _s) * ds / (ds + dr + 1e-30) + gap_tris, extra = _gmsh_fill_2d(Xall, ring, None, holes=holes, + size_of=size_of) placed = np.vstack([asm_pts, extra]) if len(extra) else asm_pts def mixed(v): @@ -6644,9 +6735,14 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, boundaries are themselves of interest. Two zones converging at a shallow angle make the overlap a spike, and its fragmented tip meshes to arbitrarily bad angles; the fused union has no such tip. - mesher : {None, "ladder"}, keyword-only + mesher : {None, "ladder", "network"}, keyword-only How the band itself is meshed. ``None`` (default) is the CAD-built - band (frontal fill). ``"ladder"`` is the STRUCTURED band with an + band (frontal fill). ``"network"`` (2-D) is the CAD-built band of + a whole NETWORK of polylines in one call — ribbons fused (touching + strands, junctions free) with every spine EMBEDDED so cuts walk + exact vertices at any resolution; the one path for kissing joins, + shared-band stepovers and plain strands alike. ``"ladder"`` is + the STRUCTURED band with an exact mid-surface vertex sheet — the mandatory choice when the band's spine/mid-surface is to be cut/split (#595: a cut through a remeshed band snaps rail vertices), and the choice that makes @@ -6685,8 +6781,9 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, if assembly not in ("fuse", "fragment"): raise ValueError( f"assembly must be 'fuse' or 'fragment', not {assembly!r}") - if mesher not in (None, "ladder"): - raise ValueError(f"mesher must be None or 'ladder', not {mesher!r}") + if mesher not in (None, "ladder", "network"): + raise ValueError( + f"mesher must be None, 'ladder' or 'network', not {mesher!r}") if dm.getDimension() == 2: return _place_thin_volume_2d(dm, patches, width, label, label_value, @@ -7440,7 +7537,8 @@ def place_fault_ribbon(base_mesh, sheet, width, *, normals=None, def place_fault_ribbon_2d(base_mesh, traces, width, *, margin_rings=2, band_label="Band", band_value=71, - clearance=0.3, split=True, verbose=False): + clearance=0.3, split=True, mesher="ladder", + spines=None, verbose=False): """Split-ready 2-D fault ribbons from the traces' OWN sampling (#629). The 2-D production fault-prep path, honouring the same contract set @@ -7465,15 +7563,37 @@ def place_fault_ribbon_2d(base_mesh, traces, width, *, margin_rings=2, boundary for ``add_fault_bc``. width : float Band thickness (split-node models: a resolution parameter). - margin_rings : int + margin_rings : int or sequence of (int, int) The band extends this many points beyond each fault end, by - linear tangent continuation. Must be >= 1 (the tip rule). + linear tangent continuation. Must be >= 1 (the tip rule). One + count for every end, or a ``(start, end)`` pair per trace. band_label, band_value : str, int - Cell label of each band; trace ``k`` gets ``band_value + k`` so - per-fault zones stay distinguishable - (``mesh.cells_labelled(band_label)`` unions them). + Cell label of the band. With ``mesher="ladder"`` trace ``k`` gets + ``band_value + k`` so per-fault zones stay distinguishable + (``mesh.cells_labelled(band_label)`` unions them); with + ``mesher="network"`` the fused band is ONE region carrying + ``band_value``. Either way ``info["band"]`` is the union mask and + ``info["footprints"]`` the per-fault ones. clearance : float Carve clearance (the measured thin-shell default 0.3). + spines : list of (label, polyline), optional + The polylines to PLACE, when they differ from the traces to cut: + a collinear abutting pair must share one spine (two ribbons + overlapping along one line interleave their vertices into + slivers), with the cut stopping at each piece's own ends. Every + trace vertex must be a vertex of some spine. Default: the + traces themselves. + mesher : {"ladder", "network"} + How the band is meshed. ``"ladder"`` (default) places each trace + SEQUENTIALLY as a structured band: bands may not touch, and level + pairs NEST (a 2:1 sub-sampled ladder shares every vertex — the + composed-hierarchy economics, #629). ``"network"`` places the + whole set in ONE call: the ribbons are fused in CAD so strands may + touch — a kissing junction, a shared stepover band — and every + spine is embedded, so the cuts still walk exact vertices at any + resolution. Choose ``"network"`` whenever the traces come within a + band width of one another; choose ``"ladder"`` when the placed + levels must nest. split : bool Cut + split each trace (``mesh.add_fault``). ``False`` returns the placed, unlabelled-fault mesh for painted (volumetric) @@ -7487,10 +7607,13 @@ def place_fault_ribbon_2d(base_mesh, traces, width, *, margin_rings=2, mesh : uw.discretisation.Mesh The fault-resolving mesh (split when ``split=True``). info : dict - ``n_cells``, ``spacing`` / ``n_rungs`` (per trace), and - ``footprints`` — per-label FAULT-footprint cell masks (the - honoured-paint rule): what painted rheology / ``fac_zone`` - keys should use, never the whole band. + ``n_cells``, ``spacing`` / ``n_rungs`` (per trace), the + ``mesher`` used, the ``extended`` spines (user samples plus the + margin) and the band ``width``; ``footprints`` — per-label + FAULT-footprint cell masks (the honoured-paint rule): what + painted rheology / ``fac_zone`` keys should use, never the whole + band; and ``band`` — the union band mask, for a band-confined + yield or a structural patch key. Notes ----- @@ -7501,61 +7624,107 @@ def place_fault_ribbon_2d(base_mesh, traces, width, *, margin_rings=2, """ from underworld3 import discretisation - if margin_rings < 1: + if spines is None: + spines = [(label, np.asarray(P, dtype=float)) for label, P in traces] + if np.ndim(margin_rings) == 0: + margin_rings = [(int(margin_rings), int(margin_rings))] * len(spines) + margin_rings = [(int(a), int(b)) for a, b in margin_rings] + if len(margin_rings) != len(spines) or min(min(m) for m in margin_rings) < 1: raise ValueError( - "margin_rings must be >= 1: the split cannot reach the band " - "rim (the tip rule); the margin is extrapolated surround.") + "margin_rings must be >= 1 at every end of every trace: the " + "split cannot reach the band rim (the tip rule); the margin is " + "extrapolated surround.") + if mesher not in ("ladder", "network"): + raise ValueError( + f"mesher must be 'ladder' or 'network', not {mesher!r}") labels = [label for label, _P in traces] if len(set(labels)) != len(labels): raise ValueError( f"trace labels must be unique (each becomes a boundary); got " f"{labels}") + if spines is None: + spines = [(label, np.asarray(P, dtype=float)) for label, P in traces] dm = base_mesh.dm spacing_all, rungs_all, extended = [], [], [] - for k, (label, P) in enumerate(traces): + for label, P in spines: P = np.asarray(P, dtype=float) if P.ndim != 2 or P.shape[1] != 2 or len(P) < 3: raise ValueError( f"trace {label!r}: expected an (n, 2) polyline with " f"n >= 3, got shape {P.shape}") - S = _extend_polyline_2d(P, margin_rings) - R = _mitred_reach_2d(S) - spacing = float(np.linalg.norm(np.diff(P, axis=0), axis=1).mean()) - dm, _info = place_thin_volume( - dm, [(S, R)], width, label=band_label, - label_value=band_value + k, clearance=clearance, - size=spacing, mesher="ladder", verbose=verbose) - spacing_all.append(spacing) + extended.append(_extend_polyline_2d(P, margin_rings[len(extended)])) + spacing_all.append( + float(np.linalg.norm(np.diff(P, axis=0), axis=1).mean())) rungs_all.append(len(P)) - extended.append(S) + + if mesher == "network": + # ONE placement call for the whole network: the ribbons are fused + # in CAD, so touching strands and shared bands are ordinary cells + # of the union, and every spine is EMBEDDED, so the cut below + # walks its own vertices (#595). The fused band is one region and + # therefore carries one label value. + dm, _info = place_thin_volume( + dm, extended, width, label=band_label, label_value=band_value, + clearance=clearance, size=float(np.mean(spacing_all)), + mesher="network", verbose=verbose) + else: + for k, S in enumerate(extended): + dm, _info = place_thin_volume( + dm, [(S, _mitred_reach_2d(S))], width, label=band_label, + label_value=band_value + k, clearance=clearance, + size=spacing_all[k], mesher="ladder", verbose=verbose) mesh = discretisation.Mesh( dm, simplex=True, qdegree=base_mesh.qdegree, coordinate_system_type=base_mesh.CoordinateSystem.coordinate_type, boundaries=base_mesh.boundaries, verbose=False) + # The placed mesh OWNS the base's multigrid tail (the coarse levels do + # not need the fault): every solver on it drives FMG automatically, + # and the band label is its FAC patch key for volumetric rheologies. + from underworld3.utilities.custom_mg import adopt_hierarchy + band_all = np.zeros(int(mesh.dm.getHeightStratum(0)[1]), dtype=bool) + for k in range(len(spines)): + band_all |= mesh.cells_labelled(band_label, band_value + k) + adopt_hierarchy(mesh, base_mesh, fac_zone=band_all) if split: # ONE network call — cut all, then split all (chained add_fault # calls do not compose: each split re-derives the pairing records # and drops the earlier fault's). mesh = mesh.add_fault([(label, np.asarray(P, dtype=float)) for label, P in traces]) + mesh._custom_mg_fac_zone = None # a split fault needs no patch # Per-strand FAULT FOOTPRINT masks (the honoured-paint rule): band - # cells whose nearest extended sample is a USER point. This is the - # mask a volumetric rheology (or a fac_zone key) should use — never - # the whole band, whose margin is extrapolated surround. + # cells whose nearest extended sample is one of THAT trace's own + # vertices. This is the mask a volumetric rheology (or a fac_zone + # key) should use — never the whole band, whose margin is + # extrapolated surround (and, on a shared spine, whose gap edges + # belong to no trace). A band cell belongs to the trace whose vertex + # is nearest to it, read off the CONCATENATED spine samples. + band = np.zeros_like(band_all) + for k in range(len(spines)): + band |= mesh.cells_labelled(band_label, band_value + k) + S_all = np.vstack(extended) + scale = 1e-9 * float(np.mean(spacing_all)) footprints = {} - for k, (label, P) in enumerate(traces): - band_k = mesh.cells_labelled(band_label, band_value + k) - m = margin_rings - is_user = np.zeros(len(extended[k]), dtype=bool) - is_user[m:len(extended[k]) - m] = True + for label, P in traces: + P = np.asarray(P, dtype=float) + d = S_all[:, None, :] - P[None, :, :] + is_user = (np.einsum("ijk,ijk->ij", d, d).min(axis=1) < scale ** 2) + if is_user.sum() != len(P): + raise ValueError( + f"trace {label!r}: {int(is_user.sum())} of its {len(P)} " + f"vertices lie on a spine; every trace vertex must be a " + f"spine vertex") footprints[label] = _footprint_from_samples( - mesh.dm, band_k, extended[k], is_user) + mesh.dm, band, S_all, is_user) info = {"n_cells": int(mesh.dm.getHeightStratum(0)[1]), "spacing": spacing_all, "n_rungs": rungs_all, - "footprints": footprints} + "footprints": footprints, "band": band, "mesher": mesher, + "extended": extended, "width": float(width), + "margin_rings": margin_rings, + "spines": [label for label, _P in spines]} if verbose: import underworld3 as _uw _uw.pprint(f"[place_fault_ribbon_2d] {info['n_cells']} cells, " diff --git a/tests/test_0857_network_mesher.py b/tests/test_0857_network_mesher.py new file mode 100644 index 00000000..c0630fae --- /dev/null +++ b/tests/test_0857_network_mesher.py @@ -0,0 +1,155 @@ +"""The NETWORK mesher (:func:`place_surface.place_thin_volume`, +``mesher="network"``): fused ribbons with EMBEDDED spines. + +Two properties the earlier paths had only one of each: touching +strands are resolved as one fused region (junctions free — the plain +fuse), and every polyline's interior sample points are vertices of the +mesh so a split cut walks them exactly (#595 — the sequential ladder, +which could not touch). Asserted on a KISSING pair: a main line and a +splay leaving its side at a gap below the band width. + +The plain fuse is the negative control: it must NOT carry the spine +vertices, or this test is not measuring the embedding it claims to. +""" +import numpy as np +import pytest +from scipy.spatial import cKDTree + +import underworld3 as uw +from underworld3.utilities.place_surface import place_thin_volume + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, + pytest.mark.skipif(uw.mpi.size > 1, + reason="serial suite")] + +WIDTH = 0.02 + + +def _box2(cell): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cell, + regular=False, qdegree=2) + + +def _kissing_pair(spacing=0.01): + """A straight main line and a splay leaving its side at gap WIDTH/2 + (kissing: inside the band), both sampled at the rung scale.""" + s = np.arange(0.0, 0.5 + 1e-9, spacing) + main = np.column_stack([0.25 + s, 0.4 + 0.3 * s]) + t = np.array([1.0, 0.3]) / np.hypot(1.0, 0.3) + n = np.array([-t[1], t[0]]) + d = np.cos(np.deg2rad(25)) * t + np.sin(np.deg2rad(25)) * n + start = main[len(main) // 2] + 0.5 * WIDTH * n + u = np.arange(0.0, 0.15 + 1e-9, spacing) + splay = start + u[:, None] * d + return [main, splay] + + +def _vertex_distance(dm, points): + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + return cKDTree(X).query(points)[0] + + +def test_a_kissing_pair_embeds_its_spines_as_vertices(): + from underworld3.utilities.line_cut import cell_areas + + mesh = _box2(0.05) + before = float(cell_areas(mesh.dm).sum()) + pair = _kissing_pair() + new, info = place_thin_volume(mesh.dm, pair, width=WIDTH, + label="Net", label_value=9, + size=0.01, mesher="network") + assert info["n_zone_cells"] > 0 + assert float(cell_areas(new).sum()) == pytest.approx(before, rel=1e-12) + assert info["min_angle"] > 10.0 + # every INTERIOR spine point is a mesh vertex (the end points sit on + # the caps and are not promised) + for P in pair: + assert _vertex_distance(new, P[1:-1]).max() < 1e-9 + + # negative control: the plain fuse places the same pair but does NOT + # carry the spine vertices + plain, _ = place_thin_volume(mesh.dm, pair, width=WIDTH, label="Net", + label_value=9, size=0.01) + assert _vertex_distance(plain, pair[0][1:-1]).max() > 1e-4, ( + "the plain fuse now carries the spine vertices; the control no " + "longer distinguishes the network path") + + +def test_a_kissing_pair_can_be_cut_along_both_spines(): + """The point of the embedding: a split cut along each strand of a + touching pair, one network add_fault call, on the fused mesh.""" + base = _box2(0.05) + pair = _kissing_pair() + dm, _ = place_thin_volume(base.dm, pair, width=WIDTH, label="Net", + label_value=9, size=0.01, mesher="network") + mesh = uw.discretisation.Mesh( + dm, simplex=True, qdegree=2, + coordinate_system_type=base.CoordinateSystem.coordinate_type, + boundaries=base.boundaries, verbose=False) + # cut the interior chains (the cap points are not vertices) + cut = mesh.add_fault([("Main", pair[0][2:-2]), + ("Splay", pair[1][2:-2])]) + assert cut.dm.getDepthStratum(0)[1] > mesh.dm.getDepthStratum(0)[1], ( + "the cut duplicated no vertices") + + +def test_an_unknown_mesher_is_refused(): + mesh = _box2(0.2) + with pytest.raises(ValueError, match="network"): + place_thin_volume(mesh.dm, _kissing_pair(), width=WIDTH, + mesher="lattice") + + +def test_a_placed_mesh_owns_the_base_hierarchy_and_the_cut_inherits_it(): + """A placement built on a refined base OWNS the base's geometric-MG + tail with the band as FAC zone; the split child inherits the tail + (the cut is the same grid) but not the zone (a split fault needs no + patch). A Stokes solve on either drives custom-P multigrid with no + set_custom_fmg — the configuration a nonlinear solve must get by + default, never GAMG by silent fallback.""" + from underworld3.utilities.place_surface import place_fault_ribbon_2d + + def fresh_base(): # a placement marks its base DM + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.125, + regular=False, qdegree=2, refinement=1) + + pair = _kissing_pair() + traces = [("Main", pair[0][2:-2]), ("Splay", pair[1][2:-2])] + base = fresh_base() + n_base = len(base._coarse_level_meshes()) + assert n_base >= 2 + placed, _info = place_fault_ribbon_2d(base, traces, WIDTH, split=False, + mesher="network") + assert placed._custom_mg_coarse_meshes is not None + assert len(placed._custom_mg_coarse_meshes) == n_base + assert placed._custom_mg_fac_zone is not None + assert int(placed._custom_mg_fac_zone.sum()) > 0 + + cut, _info = place_fault_ribbon_2d(fresh_base(), traces, WIDTH, + split=True, mesher="network") + assert cut._custom_mg_coarse_meshes is not None + assert len(cut._custom_mg_coarse_meshes) == n_base + assert cut._custom_mg_fac_zone is None + + for mesh in (placed, cut): + v = uw.discretisation.MeshVariable("v", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1, + continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + x, y = mesh.X + for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc((y - 0.5, 0.0), wall) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + stokes.solve() + assert stokes.snes.getConvergedReason() > 0 + velpc = stokes.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC() + assert velpc.getType() == "mg", ( + f"velocity block runs {velpc.getType()}: the mesh-owned " + "hierarchy was not picked up") + assert velpc.getMGLevels() == n_base + 1 + assert not getattr(stokes, "pc_fallbacks", {}) diff --git a/tests/test_0858_fault_network_realisations.py b/tests/test_0858_fault_network_realisations.py new file mode 100644 index 00000000..7a0fce1f --- /dev/null +++ b/tests/test_0858_fault_network_realisations.py @@ -0,0 +1,207 @@ +"""One fault specification, two realisations (:class:`FaultNetwork`). + +A fault is specified once — trace, hierarchy, properties — and the +realisation is a keyword. The split cuts the band and constrains the +node pairs; the weak plane leaves the same band whole and paints a +transversely-isotropic rheology on it. The property under test is that +these are realisations of ONE specification: the same call, the same +prepared pieces, and a mesh whose CELLS are identical (only the split's +duplicated vertices differ), which is what makes the two comparable. + +The negative controls are the two ways the pair could be faked: a weak +plane without a meshed layer (there is nothing to be weak), and a +directorless painted zone (an isotropic blob, not a fault). +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, + pytest.mark.skipif(uw.mpi.size > 1, + reason="serial suite")] + +H = 0.03 +WIDTH = 0.02 + + +def _trace(): + """A gently curved trace, well inside a unit box.""" + s = np.linspace(0.3, 0.7, 15) + return np.column_stack([s, 0.5 + 0.15 * (s - 0.5)]) + + +def _base(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=8 * H, + regular=False, refinement=1, qdegree=2) + + +def _network(realisation, base, width=WIDTH): + net = uw.meshing.FaultNetwork([("Main", _trace())]) + net.prepare(h=H, verbose=False) + net.build(base=base, width=width, realisation=realisation, + max_levels=1) + return net + + +def test_the_two_realisations_share_one_mesh(): + base = _base() + split = _network("split", base) + ti = _network("ti", base) + + assert split.realisation == "split" and ti.realisation == "ti" + assert [n for n, _p in split.prepared] == [n for n, _p in ti.prepared] + assert split.info["n_cells"] == ti.info["n_cells"], ( + "the realisations no longer share a mesh; the comparison between " + "them is then confounded by the discretisation") + # the ONLY difference: the split duplicated the cut's nodes + assert (split.mesh.dm.getDepthStratum(0)[1] + > ti.mesh.dm.getDepthStratum(0)[1]) + + +def test_the_band_can_be_cut_at_any_width(): + """The realisation is a free choice, not something the mesh grants: + the band is meshed around the trace's own points and segments, so a + band a tenth of the background element still holds the whole cut + chain. Measured over w/h from 1 to 1/10 in + ``~/+Simulations/fault_split_at_ti_width``; two widths here. + """ + from scipy.spatial import cKDTree + from underworld3.utilities.place_surface import place_fault_ribbon_2d + + h_far = 4 * H + for width in (h_far / 2, h_far / 10): + s = np.arange(0.3, 0.7 + 1e-12, width / 2) + P = np.column_stack([s, np.full_like(s, 0.5)]) + mesh, _info = place_fault_ribbon_2d( + _base(), [("F", P)], width, mesher="network", split=False) + X = np.asarray(mesh.dm.getCoordinatesLocal().array).reshape(-1, 2) + assert cKDTree(X).query(P)[0].max() < 1e-12, ( + f"w={width:g}: the trace is not on the mesh") + n_before = mesh.dm.getDepthStratum(0)[1] + cut = mesh.add_fault([("F", P)]) + # every interior node duplicated; the two tips stay welded + assert cut.dm.getDepthStratum(0)[1] - n_before == len(P) - 2, ( + f"w={width:g}: the cut chain is incomplete") + + +def test_the_weak_plane_is_painted_on_the_fault_footprint(): + ti = _network("ti", _base()) + eta1, ndir, foot = ti.ti_fields(eta_1=0.01, eta_0=1.0) + + assert foot.sum() > 0, "no footprint cells: nothing would be weak" + assert foot.sum() < ti.info["band"].sum() + 1 + vals = eta1.array[:, 0, 0] + assert np.allclose(vals[foot], 0.01) + assert np.allclose(vals[~foot], 1.0) + + # the director is a unit normal of the trace, cell by cell — a + # painted zone WITHOUT this is an isotropic blob, not a fault + d = np.asarray(ndir.array).reshape(-1, 2)[foot] + assert np.allclose(np.linalg.norm(d, axis=1), 1.0) + P = _trace() + t = (P[-1] - P[0]) / np.linalg.norm(P[-1] - P[0]) + assert np.abs(d @ t).max() < 0.05, ( + "the directors are not perpendicular to the trace") + + +def test_the_band_is_damage_material_in_either_realisation(): + """The band is a meshed region of material around the fault, not + scaffolding for the weak plane: the split gets the same mask and the + same band-confined yield.""" + for realisation in ("split", "ti"): + net = _network(realisation, _base()) + assert net.band.sum() > 0 + assert set(net.footprints) == {n for n, _p in net.prepared} + assert net.footprints["Main"].sum() <= net.band.sum() + + tau = net.band_yield(tau_y=4.0) + painted = net._band_yield_var.array[:, 0, 0] + assert np.allclose(painted[net.band], 4.0) + assert np.allclose(painted[~net.band], 1.0e8) + assert tau.free_symbols # a usable expression + + +def test_a_weak_plane_needs_a_layer_to_be_weak_in(): + net = uw.meshing.FaultNetwork([("Main", _trace())]) + net.prepare(h=H, verbose=False) + with pytest.raises(ValueError, match="width"): + net.build(base=_base(), realisation="ti", max_levels=1) + with pytest.raises(ValueError, match="realisation"): + net.build(base=_base(), width=WIDTH, realisation="smeared", + max_levels=1) + # the mesher choices are the ones THIS dimension offers + with pytest.raises(ValueError, match="2-D"): + net.build(base=_base(), width=WIDTH, mesher="embed", max_levels=1) + + +def test_the_solver_is_given_whichever_realisation_was_built(): + """One call imposes the network; what it imposes follows the + realisation, and asking for the other one is refused rather than + silently ignored.""" + for realisation in ("split", "ti"): + net = _network(realisation, _base()) + mesh = net.mesh + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + q = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, + continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=q) + net.apply(stokes, eta_1=0.01) + if realisation == "split": + assert getattr(stokes, "_fault_contact_faults", []) == \ + [n for n, _p in net.prepared] + else: + assert net.ti["footprint"].sum() > 0 + params = stokes.constitutive_model.Parameters + assert params.director.shape == (1, mesh.dim) + with pytest.raises(RuntimeError, match="no fault pairs"): + net.apply_contact(stokes) + + +def test_a_weak_plane_cannot_be_gauged_without_a_layer(): + net = uw.meshing.FaultNetwork([("Main", _trace())]) + net.prepare(h=H, verbose=False) + net.build(base=_base(), max_levels=1) # no width: no band + assert net.info is None + with pytest.raises(RuntimeError, match="width"): + net.ti_fields(eta_1=0.01) + + +def test_the_weak_plane_gauges_the_jump_across_its_layer(): + """The weak plane has no node pair, so its slip is the tangential + velocity jump across the layer. Read on a PRESCRIBED linear shear, + where the jump is known exactly: no solve, no solver behaviour, just + the gauge's own arithmetic and sampling. + """ + import types + + s = np.linspace(0.3, 0.7, 15) + trace = np.column_stack([s, np.full_like(s, 0.5)]) # horizontal + net = uw.meshing.FaultNetwork([("Main", trace)]) + net.prepare(h=H, verbose=False) + net.build(base=_base(), width=WIDTH, realisation="ti", max_levels=1) + + a = 3.0 # v = (a y, 0): t = x-hat + v = uw.discretisation.MeshVariable("Ug", net.mesh, net.mesh.dim, + degree=2) + v.array[:, 0, 0] = a * np.asarray(v.coords)[:, 1] + v.array[:, 0, 1] = 0.0 + + skirt = 0.5 * WIDTH + float(net.info["spacing"][0]) + got = net.slips(types.SimpleNamespace(u=v)) + assert got["Main"] == pytest.approx(a * 2 * skirt, rel=1e-6) + + +def test_the_fault_carries_its_own_properties(): + """The surface object survives the realisation and is where + fault-attached (not mesh-attached) properties live.""" + for realisation in ("split", "ti"): + net = _network(realisation, _base()) + name = net.prepared[0][0] + surf = net.surface(name) + friction = surf.add_variable("mu", size=1) + friction.data[:] = 0.6 + assert np.allclose(np.asarray(friction.data), 0.6) + with pytest.raises(KeyError): + net.surface("NotAFault") diff --git a/tests/test_0859_fault_network_junction_glue.py b/tests/test_0859_fault_network_junction_glue.py new file mode 100644 index 00000000..dfc84a28 --- /dev/null +++ b/tests/test_0859_fault_network_junction_glue.py @@ -0,0 +1,189 @@ +"""The junction glue is read off the mesh: ribbon minus cut. + +A split network is welded wherever one cut stops short of another — +a kissing branch, an abutting pair. The band (the ribbon, with its +extrapolated margins) covers those places; the cut chains do not. The +difference, restricted to where two ribbons meet, is the set of cells +the glue belongs in (:meth:`FaultNetwork.junction_cells`), and a weak +isotropic patch on them (:meth:`FaultNetwork.junction_patch`) lets the +cut carry slip through the joint. Measured on the S-fault rig +(2026-08-27); here on the smallest network that has both joint kinds. + +The negative controls: free tips are NOT junctions (damage there would +lengthen the fault rather than join it), and the patch must leave the +strand away from the joints as the split alone had it. +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, + pytest.mark.skipif(uw.mpi.size > 1, + reason="serial suite")] + +H = 0.03 +WIDTH = 0.04 # two rungs across: a band, not a line of cells + + +def _pieces(): + """A main strand, a collinear continuation abutting it across a + gap of one and a half cells (wider than the ligament, so prepare() + leaves it as it is), and a splay ending on the main (a T, trimmed + back by prepare() into a kissing junction).""" + main = np.column_stack([np.linspace(0.25, 0.50, 12), + np.full(12, 0.5)]) + cont = np.column_stack([np.linspace(0.55, 0.75, 9), + np.full(9, 0.5)]) + s = np.linspace(0.0, 1.0, 8) + splay = np.column_stack([0.38 + 0.12 * s, 0.5 + 0.18 * s]) + return [("Main", main), ("Cont", cont), ("Splay", splay)] + + +def _longer(): + """The other end member: the same fault, just LONGER — Main and its + continuation as one trace, no joint to weld.""" + main, cont, splay = (p for _n, p in _pieces()) + line = np.column_stack([np.linspace(main[0, 0], cont[-1, 0], 22), + np.full(22, 0.5)]) + return [("Main", line), ("Splay", splay)] + + +def _network(pieces=None): + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=8 * H, + regular=False, refinement=1, qdegree=2) + pieces = _pieces() if pieces is None else pieces + net = uw.meshing.FaultNetwork(pieces, + hierarchy=[n for n, _p in pieces]) + net.prepare(h=H, ligament=1.0, verbose=False) + net.build(base=base, width=WIDTH, realisation="split", max_levels=1) + return net + + +def _centroids(net): + from underworld3.utilities.place_surface import _cell_centroids_of + band = net.info["band"] + ids, cen = _cell_centroids_of(net.mesh.dm, np.ones_like(band)) + out = np.zeros((len(band), 2)) + out[ids] = cen + return out + + +def test_junction_cells_sit_at_the_joints_and_not_at_the_free_tips(): + net = _network() + cells = net.junction_cells() + cen = _centroids(net) + + assert cells.any() + assert not (cells & ~net.info["band"]).any(), "glue outside the ribbon" + + joints = [np.array([0.525, 0.5]), np.array([0.38, 0.5])] # gap, T + near_joint = np.min([np.linalg.norm(cen[cells] - j, axis=1) + for j in joints], axis=0) + assert near_joint.max() < 5 * H, "a junction cell far from any joint" + for j in joints: + assert (np.linalg.norm(cen[cells] - j, axis=1) < 3 * H).any(), ( + f"no junction cells at the joint {j}") + + free_tips = [np.array([0.25, 0.5]), np.array([0.75, 0.5]), + np.array([0.5, 0.68])] + for tip in free_tips: + assert not (np.linalg.norm(cen[cells] - tip, axis=1) < 2 * H).any(), ( + f"a free tip {tip} was treated as a junction") + + +def test_a_near_miss_is_pulled_back_to_one_ligament_not_three(): + """Two collinear pieces closer than the ligament are an offset + junction: the join opens to ONE ligament, shared between the two + ends — not a ligament on each side on top of the gap it had.""" + main = np.column_stack([np.linspace(0.25, 0.50, 12), np.full(12, 0.5)]) + cont = np.column_stack([np.linspace(0.53, 0.75, 9), np.full(9, 0.5)]) + net = uw.meshing.FaultNetwork([("Main", main), ("Cont", cont)]) + net.prepare(h=H, ligament=2.0, verbose=False) + ends = dict((n, P) for n, P in net.prepared) + gap = ends["Cont"][0, 0] - ends["Main"][-1, 0] + assert gap == pytest.approx(2.0 * H, rel=0.15), ( + f"the join opened to {gap:.3f}, not the ligament {2 * H:.3f}") + assert [j["kind"] for j in net.junctions] == ["near-miss", "near-miss"] + + +def test_a_collinear_pair_shares_one_spine_and_makes_no_slivers(): + """Two ribbons placed along one line interleave their vertices into + sliver cells (measured: 7800 cells below 1e-6 in area on the rig). + The pieces of a collinear abutting pair are placed on ONE spine, cut + at their own ends, with the gap as spine the split does not cut.""" + net = _network() + assert [n for n, _S, _i in net.spines] == ["Main+Cont", "Splay"] + main_cont = net.spines[0] + assert (main_cont[2] == -1).sum() >= 1, "no gap vertex on the shared spine" + + dm = net.mesh.dm + cS, _cE = dm.getHeightStratum(0) + ids = np.flatnonzero(net.info["band"]) + areas = np.array([dm.computeCellGeometryFVM(int(c) + cS)[0] for c in ids]) + assert areas.min() > 1e-3 * np.median(areas), ( + f"sliver cells in the band: min area {areas.min():.2e} against a " + f"median of {np.median(areas):.2e}") + + +def test_the_weak_plane_has_no_junction_cells(): + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=8 * H, + regular=False, refinement=1, qdegree=2) + net = uw.meshing.FaultNetwork(_pieces()) + net.prepare(h=H, ligament=1.0, verbose=False) + net.build(base=base, width=WIDTH, realisation="ti", max_levels=1) + with pytest.raises(RuntimeError, match="SPLIT"): + net.junction_cells() + + +def _slip_at(solver, name, x0): + """The cut's tangential jump at the pair nearest ``x = x0``.""" + from underworld3.utilities.fault_contact import fault_pair_jumps + coords, jumps, normals = fault_pair_jumps( + solver, name, solver._rotated_freeslip_info) + k = np.argmin(np.abs(coords[:, 0] - x0)) + jn = float(jumps[k] @ normals[k]) + return float(np.linalg.norm(jumps[k] - jn * normals[k])) + + +def _shear_solve(net, glue, tag): + mesh = net.mesh + x, y = mesh.X + v = uw.discretisation.MeshVariable(f"U{tag}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"P{tag}", mesh, 1, degree=1, + continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ( + net.junction_patch(eta_0=1.0) if glue else 1.0) + for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc((2.0 * (y - 0.5), 0.0), wall) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + net.apply(stokes) + net.solve(stokes) + return stokes + + +def test_the_patch_lands_between_the_cut_and_the_longer_fault(): + """The two end members: the fault CUT (abutting pieces, welded by + the intact bridge) and the fault LONGER (one continuous trace). The + glue must move the continuation's slip from the first toward the + second and not beyond it. Read at one station on the continuation, + under a shear along the fault.""" + X0 = 0.70 + cut = _network() + welded = _slip_at(_shear_solve(cut, glue=False, tag="w"), "Cont", X0) + glued = _slip_at(_shear_solve(cut, glue=True, tag="g"), "Cont", X0) + longer = _slip_at(_shear_solve(_network(_longer()), glue=False, + tag="l"), "Main", X0) + + assert longer > 1.2 * welded, ( + f"no weld to repair: cut {welded:.3f} vs longer {longer:.3f}") + assert glued > welded + 0.5 * (longer - welded), ( + f"the glue repaired less than half the weld: {welded:.3f} -> " + f"{glued:.3f} against {longer:.3f}") + assert glued < 1.1 * longer, ( + f"the glue overshot the longer fault: {glued:.3f} vs {longer:.3f}")