From e71cbd450c2abb11fff3816cf9705779d76eeed3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 13:51:56 +1000 Subject: [PATCH 1/4] Fuse overlapping fault zones instead of fragmenting them MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit place_thin_volume gains assembly={"fuse","fragment"}, defaulting to "fuse". Where two zones overlap, fragment cuts along the boundary of the overlap; for a tangential merge — a ribbon soling into another, the listric case — that boundary is a lens closing at the convergence angle, and the mesher resolves it as a chain of slivers. Measured on the sole geometry now in test_0855: minimum angle 0 degrees and 22 cells under 5, against 37 degrees and none for the fused union. Nothing downstream needs those internal boundaries. The zone carries one label, and a cell's fault properties are read from the Surface objects by proximity, not from the CAD piece it was meshed in — so the seams are conditioning cost with no physics on them. The same swap applies in 3-D (_occ_assembly_3d) and is the default there too. The CAD area/volume gate is unaffected: it is computed from the faces that survive the boolean, which are the union either way, and matches the meshed area to 1e-16 relative under both. The fragment branch is kept and tested as the negative control: a test that only asserted the fused mesh is clean could not tell whether the geometry still exercises the defect. Underworld development team with AI support from Claude Code --- .../conforming-surfaces-and-fault-zones.md | 29 +++++-- src/underworld3/utilities/place_surface.py | 75 ++++++++++++++----- tests/test_0855_place_thin_volume.py | 45 ++++++++++- 3 files changed, 121 insertions(+), 28 deletions(-) diff --git a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md index 98d2b2ab0..41b0813c5 100644 --- a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md +++ b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md @@ -300,15 +300,28 @@ representation, where the width `w` is a mesh parameter and constitutive and measured (`w = h/4` passes every gate). The construction is "mesh the whole lot, then embed", and the two stages are -forced by kernel: OCC `fragment` — the only operation that resolves -fault–fault intersections — sees only CAD entities, and the cavity fill +forced by kernel: the OCC booleans — the only operations that resolve +fault–fault intersections — see only CAD entities, and the cavity fill honours only discrete ones. So the network's patches are thickened by -`±w/2` and fragmented **together** in OCC (a junction becomes ordinary -volumes of the union — no geometric junction treatment, the rheology -decides), the assembly is meshed standalone at layer scale, its boundary -skin is extracted, and the meshed assembly is embedded whole: a cavity is -carved around it and gmsh fills the annular gap with the skin as an interior -**hole** in the fill volume, both surfaces verbatim. +`±w/2` and resolved against one another **together** in OCC (a junction +becomes ordinary volumes of the union — no geometric junction treatment, the +rheology decides), the assembly is meshed standalone at layer scale, its +boundary skin is extracted, and the meshed assembly is embedded whole: a +cavity is carved around it and gmsh fills the annular gap with the skin as an +interior **hole** in the fill volume, both surfaces verbatim. + +Which boolean resolves the overlaps is the `assembly` argument, and the +default is `"fuse"` — the union as **one** region. The alternative, +`"fragment"`, keeps every overlap piece as its own region, so the mesh must +conform to the boundary of the overlap; where two zones converge +tangentially that boundary is a lens closing at the convergence angle, and +the mesh resolves it as a chain of slivers (measured on a ribbon soling into +another: minimum angle 0°, 22 cells under 5°, against 37° and none for the +fused union). Nothing downstream needs those internal boundaries: the zone +carries a single label, and a cell's fault properties are read from the +`Surface` objects by proximity, not from the piece it was meshed in. Ask for +`"fragment"` only when the boundaries between overlapping zones are +themselves the object of interest. In the result the layer's **cells** carry `(label, value)` — the zone exists to hand cells to the rheology — and the skin's faces carry diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 31c25742d..201a17071 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -2693,8 +2693,11 @@ def _patch_frame(patch): return n -def _occ_assembly_3d(patches, width, size, box=None): - """Thicken each planar patch by ±width/2, fragment together, mesh. +def _occ_assembly_3d(patches, width, size, box=None, assembly="fuse"): + """Thicken each planar patch by ±width/2, resolve overlaps, mesh. + + ``assembly`` is :func:`place_thin_volume`'s: ``"fuse"`` returns the union + as one solid, ``"fragment"`` keeps every overlap piece. ``box = (lo, hi)`` applies the specify-long contract: the thickened solids are INTERSECTED with the domain box, so patches may protrude — @@ -2725,7 +2728,14 @@ def _occ_assembly_3d(patches, width, size, box=None): out = occ.extrude([(2, surf)], *(width * n)) solids += [t for d, t in out if d == 3] if len(solids) > 1: - occ.fragment([(3, solids[0])], [(3, t) for t in solids[1:]]) + # The 2-D lesson one dimension up: the seams fragment leaves at + # an overlap have no physics on them (properties reach the cells + # from the Surface objects), and a shallow-angle overlap gives + # them a spike to mesh. See :func:`_occ_assembly_2d`. + if assembly == "fuse": + occ.fuse([(3, solids[0])], [(3, t) for t in solids[1:]]) + else: + occ.fragment([(3, solids[0])], [(3, t) for t in solids[1:]]) if box is not None: lo, hi = (np.asarray(b, dtype=float) for b in box) occ.synchronize() @@ -3028,12 +3038,16 @@ def cap_tag_of(k): gmsh.finalize() -def _occ_assembly_2d(polylines, width, size): - """Thicken each polyline segment into a quad, fragment together, mesh. +def _occ_assembly_2d(polylines, width, size, assembly="fuse"): + """Thicken each polyline into a ribbon, resolve overlaps, mesh. - The 2-D thin volume: a ribbon is the union of one quad per polyline - segment, kinks and crossings resolved by ``fragment`` exactly as the 3-D - junctions are. Returns ``(points, triangles, cad_area)``. + The 2-D thin volume: a ribbon is the mitred outline of one polyline, and + the ribbons of a network are resolved against one another in CAD. + ``assembly`` chooses that resolution: ``"fuse"`` returns the union as ONE + face, ``"fragment"`` keeps every overlap piece as its own face. Both mesh + the same region; they differ in the internal seams the mesher must + honour. Returns ``(points, triangles, cad_area)``, the area being that of + the resolved faces — the union — under either choice. """ import gmsh @@ -3088,7 +3102,19 @@ def outline(P): if not surfs: raise ValueError("the polylines contain no segment to thicken") if len(surfs) > 1: - occ.fragment([(2, surfs[0])], [(2, t) for t in surfs[1:]]) + # Where ribbons overlap, fragment's seams are the boundary of the + # overlap, and for a tangential merge that boundary is a lens + # closing at the convergence angle — a chain of slivers the + # mesher must resolve (measured on the sole geometry of + # test_0855: minimum angle 0 degrees and 22 cells under 5, + # against 37 degrees and none fused). The seams carry no physics + # to preserve: the zone is one label, and a cell's fault + # properties are read from the Surface objects by proximity, not + # from the piece it was meshed in. + if assembly == "fuse": + occ.fuse([(2, surfs[0])], [(2, t) for t in surfs[1:]]) + else: + occ.fragment([(2, surfs[0])], [(2, t) for t in surfs[1:]]) occ.synchronize() faces = gmsh.model.getEntities(2) @@ -3807,7 +3833,7 @@ def remove_embedded(dm, label, label_value=1, clearance=0.6, verbose=False): def _place_thin_volume_2d(dm, polylines, width, label, label_value, - clearance, size, verbose): + clearance, size, assembly, verbose): """The ribbon: the identical construction one dimension down. Serial AND parallel through the same gather-first mechanism as the 3-D @@ -3823,7 +3849,7 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, if comm.rank == 0: try: asm_pts, asm_tris, cad_area = _occ_assembly_2d(polylines, width, - size) + size, assembly) 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])) @@ -4057,14 +4083,15 @@ def mixed(v): def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, - clearance=0.7, size=None, verbose=False): + clearance=0.7, size=None, assembly="fuse", + verbose=False): """Embed a THIN VOLUME of the given width around each patch, junctions free. The finite-width fault representation: each planar patch is thickened by ``±width/2``, the thickened volumes of the whole network are resolved - against one another with OCC ``fragment`` — a junction becomes ordinary - cells of the union, no geometric treatment, the rheology decides — the - assembly is meshed standalone at layer scale (sub-``h`` widths are the + against one another in OCC — a junction becomes ordinary cells of the + union, no geometric treatment, the rheology decides — the assembly is + meshed standalone at layer scale (sub-``h`` widths are the point: ``V = 2 ε̇ w`` makes the width constitutive), and the meshed assembly is embedded whole into the existing mesh: a cavity is carved around it and gmsh fills the annular gap with the assembly's boundary @@ -4102,6 +4129,17 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, skin. size : float or None The layer's own mesh size; ``None`` takes ``0.9 * width``. + assembly : {"fuse", "fragment"} + How overlapping zones are resolved in CAD before meshing. ``"fuse"`` + (the default) returns the union as one region with no internal seam; + ``"fragment"`` keeps each overlap piece as its own region, so the + mesh conforms to the boundaries of the overlap. The zone mesh carries + one label either way — a cell's fault properties come from the + :class:`Surface` objects, not from the piece it was meshed in — so + ``"fragment"`` is worth asking for only when those internal + 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. verbose : bool Report the counts. @@ -4125,10 +4163,13 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, if width <= 0.0: raise ValueError("width must be positive") size = 0.9 * width if size is None else float(size) + if assembly not in ("fuse", "fragment"): + raise ValueError( + f"assembly must be 'fuse' or 'fragment', not {assembly!r}") if dm.getDimension() == 2: return _place_thin_volume_2d(dm, patches, width, label, label_value, - clearance, size, verbose) + clearance, size, assembly, verbose) if dm.getDimension() != 3: raise NotImplementedError( f"place_thin_volume takes a 2-D or 3-D simplex mesh; this mesh " @@ -4151,7 +4192,7 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, if comm.rank == 0: try: asm_pts, asm_tets, cad_vol = _occ_assembly_3d( - patches, width, size, box=(box_lo, box_hi)) + patches, width, size, box=(box_lo, box_hi), assembly=assembly) v6 = np.einsum( "ij,ij->i", np.cross(asm_pts[asm_tets][:, 1] - asm_pts[asm_tets][:, 0], diff --git a/tests/test_0855_place_thin_volume.py b/tests/test_0855_place_thin_volume.py index 803fe2ec9..2e937c1b4 100644 --- a/tests/test_0855_place_thin_volume.py +++ b/tests/test_0855_place_thin_volume.py @@ -1,8 +1,9 @@ """The embedded thin-volume fault mesh (:func:`place_surface.place_thin_volume`). -The finite-width representation: patches are thickened by ±width/2 in OCC, -``fragment`` resolves the network's junctions in CAD — the only kernel that -can — the assembly is meshed standalone at layer scale, and the meshed +The finite-width representation: patches are thickened by ±width/2 in OCC and +resolved against one another in CAD — the only kernel that can, by ``fuse`` +into one region or ``fragment`` into the overlap pieces — the assembly is +meshed standalone at layer scale, and the meshed assembly is embedded whole into the existing mesh: cavity carved, annular gap filled by gmsh with the assembly's skin as an interior HOLE, both constraint surfaces verbatim. Junctions need no geometric treatment: they are ordinary @@ -174,6 +175,44 @@ def test_a_kinked_ribbon_does_not_sliver(): assert info["min_angle"] > 10.0 +def test_a_tangential_merge_fuses_instead_of_slivering(): + """The sole: a ribbon converging onto another until the two coincide. + + Where two zones overlap, ``fragment`` cuts along the boundary of the + overlap, and for a tangential merge that boundary is a lens whose tips + close at the convergence angle — the mesh must resolve a chain of + slivers, and does (measured: 22 cells under 5 degrees, one of them + degenerate). ``fuse`` returns the union as one face with no such + boundary. The zone carries one label either way, so nothing downstream + can tell the difference except the conditioning. + + The ``fragment`` branch is the negative control: it must show the + slivers, or this test is not measuring what it claims to. + """ + from underworld3.utilities.line_cut import min_angles + + sole = [np.array([[0.20, 0.50], [0.80, 0.50]]), + np.array([[0.20, 0.56], [0.50, 0.505], [0.80, 0.50]])] + + mesh = _box2(0.05) + fused, _ = place_thin_volume(mesh.dm, sole, width=0.02, label="Sole", + label_value=7) + assert float(min_angles(fused).min()) > 5.0 + + torn, _ = place_thin_volume(mesh.dm, sole, width=0.02, label="Sole", + label_value=7, assembly="fragment") + assert int((min_angles(torn) < 5.0).sum()) > 5, ( + "the fragmented merge did not sliver; the geometry no longer " + "exercises the defect this test exists for") + + +def test_an_unknown_assembly_boolean_is_refused(): + mesh = _box2(0.2) + with pytest.raises(ValueError, match="fuse.*fragment"): + place_thin_volume(mesh.dm, [np.array([[0.3, 0.5], [0.7, 0.5]])], + width=0.02, assembly="union") + + def test_a_second_zone_leaves_the_first_intact(): mesh = _box2(0.05) l1 = np.array([[0.3, 0.35], [0.7, 0.65]]) From b8fe93719ae0b36e4f697f2477c70daf5b20cdfd Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 14:45:08 +1000 Subject: [PATCH 2/4] Record the hybrid fault-zone architecture as a design note MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The plan for this work lived outside the repository. The architecture belongs in docs/developer/design/ where collaborators can read it: the two representations and what separates them, the tip-pinning mechanism behind the seam defect, the fuse ruling landed here, and the geometric slice criterion that is proposed and not yet built. Marked throughout as to what is implemented and what is not, and it records two things that are decided but not obvious: the CAD area gate did not need the change the plan expected, and the question of what owns a cell where two zones overlap is still open — contacts refuse to overlap and pay a gap, zones overlap freely and resolve ownership by nearest surface vertex, which is neither continuous nor independent of how the traces were sampled. Underworld development team with AI support from Claude Code --- .../design/fault-zone-hybrid-architecture.md | 221 ++++++++++++++++++ docs/developer/index.md | 1 + 2 files changed, 222 insertions(+) create mode 100644 docs/developer/design/fault-zone-hybrid-architecture.md diff --git a/docs/developer/design/fault-zone-hybrid-architecture.md b/docs/developer/design/fault-zone-hybrid-architecture.md new file mode 100644 index 000000000..827cd2819 --- /dev/null +++ b/docs/developer/design/fault-zone-hybrid-architecture.md @@ -0,0 +1,221 @@ +# Hybrid fault zones: contacts along the fault, TI zones at the junctions + +Underworld3 represents a fault two ways, and neither is right everywhere. +This note records the design that uses both in one model, with the handover +between them chosen from the geometry rather than by the modeller. + +Status: the zone assembly change described under +[Fusing the zone](#fusing-the-zone) is implemented. Everything under +[Choosing the handover](#choosing-the-handover-from-the-geometry) is proposed +and not yet built. Numbers quoted as measured come from the exploratory +campaign in `~/+Simulations/listric_extension/`. + +## The two representations + +A **zero-thickness contact** ([split-node](../../advanced/split-node-faults.md)) +is a surface across which the velocity jumps. It is cheap, and it matches the +fault assumptions exactly: a fault of no width, carrying a slip discontinuity +and an interface law. + +A **finite-width zone** ([`place_thin_volume`](../subsystems/conforming-surfaces-and-fault-zones.md)) +is a band of weak material of a width you choose. It is more expensive, and it +applies where the contact assumptions fail — at junctions, at merges, anywhere +the fault has a width or the geometry has no unique centreline. + +On a single straight fault the two agree to 3.5%. Their cost per solve, on the +same mesh, was measured at 3.2 s for the contact, 13.5 s for a +transversely-isotropic (TI) band, and about 296 s for an isotropic band of the +same width. The isotropic band is the one to avoid; the TI band buys back most +of the cost by making the weakness directional instead of resolving a +viscosity contrast in every direction. + +## The seam defect + +Abutting the two representations end to end does not work. A hybrid built that +way returned 74.6% of the continuous-contact reference — worse than either +pure representation on its own. + +The mechanism is **tip pinning**. A contact that stops at the edge of a zone +terminates as a free crack tip, and slip is forced to zero there. The zone +then has to re-accelerate the material it receives, and the model loses the +slip that the pinned tip refused to carry. + +A tip that ends *inside* weak material is not pinned. Running each contact one +cell **into** the zone recovers 96.3% of the reference, saturating near 99% +with further penetration. One mechanism accounts for the 74.6%, for the cure, +and for the residual seen in the earlier junction work. + +This is the finding that makes a hybrid worth building: the handover is not +inherently lossy, it was being built at the wrong place. + +## Fusing the zone + +When a network's zones are thickened and resolved against one another in CAD, +the boolean that resolves them decides what the mesher has to honour. +`fragment` keeps every overlap piece as its own region, so the mesh conforms to +the boundary of the overlap. Where two faults converge tangentially — a splay +soling into a detachment, the listric case — that boundary is a lens closing +at the convergence angle, and the mesh resolves it as a chain of slivers. On +one measured pair the fragmented assembly meshed to a minimum angle of 0.13 +degrees with 52 cells under 5 degrees, and the Stokes SNES diverged on it. + +`fuse` returns the union as one region with no internal seam. The same +geometry meshed at 26.5 degrees minimum with no cell under 15, in fewer +triangles, and solved. + +Nothing downstream needs those seams. A zone carries one label, and a cell's +fault properties are read from the `Surface` objects by proximity, not from +the CAD piece the cell was meshed in. The zone mesh never has to know which +branch a cell came from, and its not knowing is what lets the connection +region select itself. + +`place_thin_volume` therefore takes `assembly={"fuse", "fragment"}` and +defaults to `"fuse"`, in both 2-D and 3-D. + +```{note} +The CAD area gate did not need changing, contrary to the expectation this +work started from. A fused Y does have less area than the sum of its ribbons, +but `cad_area` is computed from the faces that *survive* the boolean, and +those are the union under either choice. The measured CAD-against-meshed +relative difference is `1e-16` for both. +``` + +## Choosing the handover from the geometry + +After fusing, the zone is one region and there is no marker saying which part +of it is a well-defined fault and which part is a junction. We propose to read +that off the geometry. + +For a point on a fault trace, let $d(x)$ be the distance from the trace to the +**zone boundary**, and let $w$ be the zone width. + +- $d \approx w/2$ — the trace is the local medial axis of a single band. The + slip surface is well defined, so **cut a contact here**. +- $d > w/2 + \text{tol}$ — the stem is wider than one band, so two or more + faults overlap and no unique centreline exists. **Leave it to the TI zone.** + +The rule imposes no hierarchy. It does not need to know which fault is senior, +it is symmetric in the branches, and it is computable directly from the fused +outline with machinery `Surface` already has in 2-D. It is the formal version +of the observation that a fat Y's centreline departs from each trace somewhere +along the stem. + +Two other formulations are worth testing against it: the true medial axis of +the fused polygon, with the criterion becoming "the medial axis departs from +the trace by more than tol"; and the local zone width measured normal to the +trace. We prefer the distance-to-boundary form unless it misbehaves at the +tangency corner. + +Cutting the contact inside the zone needs no new meshing primitive. `add_fault` +is a *cut*, not a placement, so it slices the zone's own cells and the zone +survives — cell counts went from 72 to 74 in the measured case, the increase +being duplication at the split. + +### What owns a cell where two zones overlap + +The criterion above says where to stop cutting contacts. It does not say what +the rheology should be in the stem, and the two representations answer that +question in opposite ways today. + +A **contact** cannot overlap another contact. `fault_split` refuses a fault +that terminates on an already-split fault's slit, and gives the reason: a +shared point would clamp every arm's slip to zero, which is stiffer than a +true junction. So no node ever carries two interface laws — the network forces +the faults apart into the offset form first, and pays the gap. That gap is the +tip pinning described above. + +A **zone** may overlap freely, and carries a single label, so the question +moves into the property lookup. `SurfaceCollection.compute_nearest_fields` +gives each node the normal and identifier of the nearest surface **vertex**, +so the nearer trace wins node by node. Two consequences worth deciding on +before the criterion lands: + +- the partition boundary is the medial axis between the traces, which carries + no physics, and the director flips discontinuously across it; +- it is nearest *vertex* rather than nearest *surface*, so the boundary moves + when a trace is re-sampled without being moved (issue #544). + +It is benign for a near-tangential merge, where both directors nearly agree. +It is not benign at a high-angle junction, where two faults produce a region +weak in **two** directions. That is orthotropic rather than transversely +isotropic, and no single director represents it, so choosing a nearest fault +chooses which of the two weaknesses to discard. + +### The intended API + +```python +net = uw.meshing.FaultNetwork([...]) # Surface objects +mesh = net.prepare(h=...).build(zone="fuse", slice="auto") +constitutive_model, director = net.ti_rheology(v, eta_1=...) +net.apply_contact(stokes) # only on the sliced pieces +``` + +`slice="auto"` applies the criterion above; `slice=None` gives the pure +finite-width model. The two share one mesh, which is what makes the comparison +between them clean. + +## 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. + +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 +was blocked by a defect in which the placed vertex could be taken from a +neighbouring cell, leaving the cell that held the tip without a corner on the +line (issue #542, fixed). It appeared inside zones rather than on plain meshes +because a ribbon's vertex rows sit at exactly the half-width and manufacture +near-ties. + +3-D is the same argument with more mechanical work. The `fuse` change already +applies. The centreline becomes a medial *surface*, and the criterion is +unchanged — distance from the patch to the zone boundary against $w/2$. +Cutting is `place_sheet` plus `split_fault` rather than `cut_along_lines`, and +the tip condition becomes a rim condition, which is likely harder. +`SurfaceCollection.transfer_normals` already works in 3-D, since triangulated +surfaces have genuine cell normals, so the director may come free there. + +We do 2-D end to end before starting 3-D. + +## How we will know it works + +- **Straight-fault control.** One fault, a meshed zone in the middle, + contacts either side. The answer is known: 100% of the continuous-contact + reference, with no notch at either seam. Any implementation must reproduce + it. +- **Self-selection sweep.** Vary the fault separation so that the correct + partition of slip between the two faults changes, and require one fixed + geometric rule to track a gmsh-union control across the range. Agreement at + a single separation proves little, because that geometry may admit only one + sensible partition. +- **Traction continuity.** $\tau_{nt}$ and $\sigma_{nn}$ are continuous across + the zone boundary; $\tau_{tt}$ and $p$ may jump. Checking the jump at + stations along the fault is an oracle that does not depend on the + implementation. + +```{warning} +The ratio of $\tau_{nt}$ inside the zone to outside it is **not** a valid +metric, and cost real time in the exploratory work. Traction is continuous +across the zone boundary, so that ratio compares different regions of the +domain and resolves onto a meaningless plane away from the trace. It also +changes with the discretisation: the same ratio read 10x at cell centroids and +2.7x projected to P1, because continuity blends the value across one to two +cells. Quote cell values for metrics and P1 for figures. +``` + +## Related + +- [Conforming surfaces and fault zones](../subsystems/conforming-surfaces-and-fault-zones.md) + — the placement family, including `place_thin_volume` and the `assembly` + argument. +- [Fault contact deployment](FAULT_CONTACT_DEPLOYMENT_2026-08.md) — the + offset-junction convention this design replaces at junctions. +- Issue #539 — the TI compliance tensor (the inverse has the same structure + with $\eta \to 1/(4\eta)$), needed for plasticity's compliance-to-stiffness + map and to infer strain rate from recovered stress. +- Issue #540 — the nearest-fault director in 2-D. diff --git a/docs/developer/index.md b/docs/developer/index.md index 29c6cb28c..356a4f74c 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -163,6 +163,7 @@ design/PROJECTED_NORMALS_API_DESIGN design/TURBULENCE_MODEL_DESIGN design/declined-coord-units-proposal design/nonlinear-solver-homotopy-warmstart +design/fault-zone-hybrid-architecture ``` ```{toctree} From 2ee633fdf3a1237a674ead40613122ab56289698 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 18:22:33 +1000 Subject: [PATCH 3/4] Record what slicing costs and what it does not MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Four settings on one mesh, measured. Keeping the weak zone under the contacts reproduces the pure-zone answer to 99% and solves about three times faster; stripping the zone back to the handover buys a further factor of two and loses a quarter of the slip on the fault that receives through the merge. Those are the same mechanism. With the zone intact a contact tip ends inside weak material and is never pinned; stripped, it ends at the weak material's edge, which is the abutting case measured at 74.6% on the straight-fault control. So the criterion only has to be geometrically right for the stripped variant. With the zone kept, slicing less is slower, slicing more is faster, and slicing past what the split machinery supports is an explicit refusal — never a wrong answer. Also records a negative result, so it is not re-derived: the three-times speedup is NOT contrast-driven conditioning relief. It is flat across three decades of viscosity contrast and the zone-only cost barely moves, because transverse isotropy had already removed that penalty. The reason is still open, and the rotated path's iteration counters read zero, which is what has to be fixed to settle it. Underworld development team with AI support from Claude Code --- .../design/fault-zone-hybrid-architecture.md | 42 +++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/docs/developer/design/fault-zone-hybrid-architecture.md b/docs/developer/design/fault-zone-hybrid-architecture.md index 827cd2819..cd92d84f1 100644 --- a/docs/developer/design/fault-zone-hybrid-architecture.md +++ b/docs/developer/design/fault-zone-hybrid-architecture.md @@ -48,6 +48,48 @@ and for the residual seen in the earlier junction work. This is the finding that makes a hybrid worth building: the handover is not inherently lossy, it was being built at the wrong place. +### Two ways to slice, and only one of them makes a seam + +Measured on the listric pair at a trace separation of 0.30, four settings on +one mesh (5110 cells in every case; slicing adds two, from duplication at the +split): + +| setting | contacts | weak zone | solve | slip against control | +|---|---|---|---|---| +| control | none | everywhere | 20.4 s | — | +| sliced | to the handover | everywhere | 4.7 s | 99.1% / 99.2% | +| stripped | to the handover | only across the merge | 2.7 s | 92.0% / 75.4% | +| stripped, penetrating | one cell past it | only across the merge | 2.5 s | 92.6% / 86.7% | + +The three sliced rows are one mechanism seen from three sides. Where the zone +is left intact, a contact's tip ends *inside* weak material by construction, +so it is never pinned and the answer is the control's. Where the zone is +stripped back to the handover, the tip ends exactly at the weak material's +edge — the abutting case — and slip falls to 75%. Running the tip one cell +further in recovers part of it. + +This matters for what the slice criterion IS. If the zone is kept, slicing +cannot produce a wrong answer: slicing less is slower, slicing more is faster, +and slicing where the split machinery cannot go is an explicit refusal. The +criterion is a performance knob with a hard stop. Only the stripped variant, +which buys a further factor of about two, needs the handover to be +geometrically right, and needs penetration. + +```{note} +Keeping the contacts *and* the zone is roughly three times faster than the +zone alone, and we do not yet know why. The natural explanation — that the +contact supplies the velocity discontinuity the solver would otherwise build +out of the viscosity contrast — is measurably wrong: the advantage is flat +across three decades of contrast (2.8x at 1e-3, 3.1x at 1e-5) and the +zone-only cost barely moves with it. Transverse isotropy had already removed +that conditioning penalty; the ill-conditioned form was the *isotropic* band, +at 296 s against the TI band's 13.5 s. The likely candidate is that the two +configurations take different solver paths with different outer iteration +counts. Settling it needs the rotated solve's own KSP instrumented — +``snes.getLinearSolveIterations()`` reads zero for it, because the rotated +path builds its own prefixed KSP. +``` + ## Fusing the zone When a network's zones are thickened and resolved against one another in CAD, From 92d7f5ba25f02896e1d0e6358d7ed81fa3980f0b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 13:48:24 +1000 Subject: [PATCH 4/4] Make the assembly argument keyword-only MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Review finding on #541. `assembly` was inserted ahead of `verbose`, so a caller passing `verbose` positionally would have had it bound to `assembly` — silently, and with a value that fails the new validation only if it happens not to be "fuse" or "fragment". Both are now keyword-only, which is what they should have been. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/place_surface.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 201a17071..672961ac5 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -4083,7 +4083,7 @@ def mixed(v): def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, - clearance=0.7, size=None, assembly="fuse", + clearance=0.7, size=None, *, assembly="fuse", verbose=False): """Embed a THIN VOLUME of the given width around each patch, junctions free. @@ -4129,8 +4129,10 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, skin. size : float or None The layer's own mesh size; ``None`` takes ``0.9 * width``. - assembly : {"fuse", "fragment"} - How overlapping zones are resolved in CAD before meshing. ``"fuse"`` + assembly : {"fuse", "fragment"}, keyword-only + How overlapping zones are resolved in CAD before meshing. Keyword-only + so that inserting it ahead of ``verbose`` cannot rebind a positional + ``verbose`` from an existing caller. ``"fuse"`` (the default) returns the union as one region with no internal seam; ``"fragment"`` keeps each overlap piece as its own region, so the mesh conforms to the boundaries of the overlap. The zone mesh carries