diff --git a/articles/putting-a-fault-in-a-mesh/examples/anatomy_mesh.py b/articles/putting-a-fault-in-a-mesh/examples/anatomy_mesh.py new file mode 100644 index 0000000..aaf2797 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/anatomy_mesh.py @@ -0,0 +1,90 @@ +"""The one mesh behind the anatomy figures (band-anatomy, split-anatomy). + +A structured triangulation of [0,3] x [0,1], twelve cells by four, and a +gently curved fault trace across its middle from x = 0.5 to x = 2.5. The +CONFORMING mesh is the grid bent so that the trace runs along its edges +(the row y = 0.5 lifted onto the curve, the displacement tapering to +nothing at the top and bottom walls) — the remeshing, in miniature. The +non-conforming panel keeps the grid flat and lets the same curve cross it. +""" +import math + +H = 0.25 +NX, NY = 12, 4 +AMP = 0.12 # the curve's rise at mid-span +X0, X1 = 0.5, 2.5 # the trace's span +W = 0.5 # the band width: one cell either side + + +def rise(x): + """The trace's offset from y = 0.5: a half sine over its span, zero + at both tips and beyond them.""" + if x <= X0 or x >= X1: + return 0.0 + return AMP * math.sin(math.pi * (x - X0) / (X1 - X0)) + + +def curve(x): + return [x, 0.5 + rise(x)] + + +def normal(x): + """Unit normal to the trace at x (pointing +y).""" + if x <= X0 or x >= X1: + return [0.0, 1.0] + dy = AMP * math.pi / (X1 - X0) * math.cos(math.pi * (x - X0) / (X1 - X0)) + n = math.hypot(1.0, dy) + return [-dy / n, 1.0 / n] + + +def grid(): + """The flat grid: (coords, tris, vid) with tris as vertex triples.""" + verts, coords = {}, [] + + def vid(i, j): + if (i, j) not in verts: + verts[(i, j)] = len(coords) + coords.append([i * H, j * H]) + return verts[(i, j)] + + tris = [] + for i in range(NX): + for j in range(NY): + a, b = vid(i, j), vid(i + 1, j) + c, d = vid(i + 1, j + 1), vid(i, j + 1) + tris.append([a, b, c]) + tris.append([a, c, d]) + return coords, tris, vid + + +def conforming(coords): + """Bend the grid onto the curve: each vertex rises by the trace's + offset at its x, scaled by (1 - |y - 0.5| / 0.5) so the walls stay + put. The row y = 0.5 lands exactly on the curve.""" + out = [] + for x, y in coords: + out.append([x, y + rise(x) * (1.0 - abs(y - 0.5) / 0.5)]) + return out + + +def centroids(coords, tris): + return [[sum(coords[v][0] for v in t) / 3.0, + sum(coords[v][1] for v in t) / 3.0] for t in tris] + + +def chain(vid): + """The trace's vertices on the conforming mesh: the middle row over + the span, first tip to second.""" + return [vid(i, NY // 2) for i in range(int(X0 / H), int(X1 / H) + 1)] + + +def distance_to_curve(p, n_samp=400): + """Distance from p to the trace, and the x of the nearest sample.""" + best, bx = float("inf"), X0 + for k in range(n_samp + 1): + x = X0 + (X1 - X0) * k / n_samp + c = curve(x) + d = math.hypot(p[0] - c[0], p[1] - c[1]) + if d < best: + best, bx = d, x + return best, bx diff --git a/articles/putting-a-fault-in-a-mesh/examples/band-anatomy-data.json b/articles/putting-a-fault-in-a-mesh/examples/band-anatomy-data.json new file mode 100644 index 0000000..671258b --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/band-anatomy-data.json @@ -0,0 +1 @@ +{"flat": [[0.0, 0.0], [0.25, 0.0], [0.25, 0.25], [0.0, 0.25], [0.25, 0.5], [0.0, 0.5], [0.25, 0.75], [0.0, 0.75], [0.25, 1.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.25], [0.5, 0.5], [0.5, 0.75], [0.5, 1.0], [0.75, 0.0], [0.75, 0.25], [0.75, 0.5], [0.75, 0.75], [0.75, 1.0], [1.0, 0.0], [1.0, 0.25], [1.0, 0.5], [1.0, 0.75], [1.0, 1.0], [1.25, 0.0], [1.25, 0.25], [1.25, 0.5], [1.25, 0.75], [1.25, 1.0], [1.5, 0.0], [1.5, 0.25], [1.5, 0.5], [1.5, 0.75], [1.5, 1.0], [1.75, 0.0], [1.75, 0.25], [1.75, 0.5], [1.75, 0.75], [1.75, 1.0], [2.0, 0.0], [2.0, 0.25], [2.0, 0.5], [2.0, 0.75], [2.0, 1.0], [2.25, 0.0], [2.25, 0.25], [2.25, 0.5], [2.25, 0.75], [2.25, 1.0], [2.5, 0.0], [2.5, 0.25], [2.5, 0.5], [2.5, 0.75], [2.5, 1.0], [2.75, 0.0], [2.75, 0.25], [2.75, 0.5], [2.75, 0.75], [2.75, 1.0], [3.0, 0.0], [3.0, 0.25], [3.0, 0.5], [3.0, 0.75], [3.0, 1.0]], "bent": [[0.0, 0.0], [0.25, 0.0], [0.25, 0.25], [0.0, 0.25], [0.25, 0.5], [0.0, 0.5], [0.25, 0.75], [0.0, 0.75], [0.25, 1.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.25], [0.5, 0.5], [0.5, 0.75], [0.5, 1.0], [0.75, 0.0], [0.75, 0.27296100594190537], [0.75, 0.5459220118838107], [0.75, 0.7729610059419054], [0.75, 1.0], [1.0, 0.0], [1.0, 0.29242640687119287], [1.0, 0.5848528137423857], [1.0, 0.7924264068711928], [1.0, 1.0], [1.25, 0.0], [1.25, 0.3054327719506772], [1.25, 0.6108655439013544], [1.25, 0.8054327719506772], [1.25, 1.0], [1.5, 0.0], [1.5, 0.31], [1.5, 0.62], [1.5, 0.81], [1.5, 1.0], [1.75, 0.0], [1.75, 0.3054327719506772], [1.75, 0.6108655439013544], [1.75, 0.8054327719506772], [1.75, 1.0], [2.0, 0.0], [2.0, 0.29242640687119287], [2.0, 0.5848528137423857], [2.0, 0.7924264068711928], [2.0, 1.0], [2.25, 0.0], [2.25, 0.2729610059419054], [2.25, 0.5459220118838108], [2.25, 0.7729610059419054], [2.25, 1.0], [2.5, 0.0], [2.5, 0.25], [2.5, 0.5], [2.5, 0.75], [2.5, 1.0], [2.75, 0.0], [2.75, 0.25], [2.75, 0.5], [2.75, 0.75], [2.75, 1.0], [3.0, 0.0], [3.0, 0.25], [3.0, 0.5], [3.0, 0.75], [3.0, 1.0]], "tris": [[0, 1, 2], [0, 2, 3], [3, 2, 4], [3, 4, 5], [5, 4, 6], [5, 6, 7], [7, 6, 8], [7, 8, 9], [1, 10, 11], [1, 11, 2], [2, 11, 12], [2, 12, 4], [4, 12, 13], [4, 13, 6], [6, 13, 14], [6, 14, 8], [10, 15, 16], [10, 16, 11], [11, 16, 17], [11, 17, 12], [12, 17, 18], [12, 18, 13], [13, 18, 19], [13, 19, 14], [15, 20, 21], [15, 21, 16], [16, 21, 22], [16, 22, 17], [17, 22, 23], [17, 23, 18], [18, 23, 24], [18, 24, 19], [20, 25, 26], [20, 26, 21], [21, 26, 27], [21, 27, 22], [22, 27, 28], [22, 28, 23], [23, 28, 29], [23, 29, 24], [25, 30, 31], [25, 31, 26], [26, 31, 32], [26, 32, 27], [27, 32, 33], [27, 33, 28], [28, 33, 34], [28, 34, 29], [30, 35, 36], [30, 36, 31], [31, 36, 37], [31, 37, 32], [32, 37, 38], [32, 38, 33], [33, 38, 39], [33, 39, 34], [35, 40, 41], [35, 41, 36], [36, 41, 42], [36, 42, 37], [37, 42, 43], [37, 43, 38], [38, 43, 44], [38, 44, 39], [40, 45, 46], [40, 46, 41], [41, 46, 47], [41, 47, 42], [42, 47, 48], [42, 48, 43], [43, 48, 49], [43, 49, 44], [45, 50, 51], [45, 51, 46], [46, 51, 52], [46, 52, 47], [47, 52, 53], [47, 53, 48], [48, 53, 54], [48, 54, 49], [50, 55, 56], [50, 56, 51], [51, 56, 57], [51, 57, 52], [52, 57, 58], [52, 58, 53], [53, 58, 59], [53, 59, 54], [55, 60, 61], [55, 61, 56], [56, 61, 62], [56, 62, 57], [57, 62, 63], [57, 63, 58], [58, 63, 64], [58, 64, 59]], "cent_a": [[0.16666666666666666, 0.08333333333333333], [0.08333333333333333, 0.16666666666666666], [0.16666666666666666, 0.3333333333333333], [0.08333333333333333, 0.4166666666666667], [0.16666666666666666, 0.5833333333333334], [0.08333333333333333, 0.6666666666666666], [0.16666666666666666, 0.8333333333333334], [0.08333333333333333, 0.9166666666666666], [0.4166666666666667, 0.08333333333333333], [0.3333333333333333, 0.16666666666666666], [0.4166666666666667, 0.3333333333333333], [0.3333333333333333, 0.4166666666666667], [0.4166666666666667, 0.5833333333333334], [0.3333333333333333, 0.6666666666666666], [0.4166666666666667, 0.8333333333333334], [0.3333333333333333, 0.9166666666666666], [0.6666666666666666, 0.09098700198063513], [0.5833333333333334, 0.17432033531396848], [0.6666666666666666, 0.3562943392752387], [0.5833333333333334, 0.4319740039612703], [0.6666666666666666, 0.6062943392752388], [0.5833333333333334, 0.6743203353139684], [0.6666666666666666, 0.840987001980635], [0.5833333333333334, 0.9166666666666666], [0.9166666666666666, 0.09747546895706428], [0.8333333333333334, 0.1884624709376994], [0.9166666666666666, 0.38341340885182795], [0.8333333333333334, 0.4679119438560339], [0.9166666666666666, 0.6410670774991297], [0.8333333333333334, 0.7037698082323031], [0.9166666666666666, 0.855129137604366], [0.8333333333333334, 0.9243203353139684], [1.1666666666666667, 0.10181092398355907], [1.0833333333333333, 0.19928639294062336], [1.1666666666666667, 0.4029082409077415], [1.0833333333333333, 0.49604825483831094], [1.1666666666666667, 0.6670503765314725], [1.0833333333333333, 0.7275706641880854], [1.1666666666666667, 0.86595305960729], [1.0833333333333333, 0.9308088022903975], [1.4166666666666667, 0.10333333333333333], [1.3333333333333333, 0.20514425731689237], [1.4166666666666667, 0.411810923983559], [1.3333333333333333, 0.5120994386173439], [1.4166666666666667, 0.6802885146337848], [1.3333333333333333, 0.7420994386173438], [1.4166666666666667, 0.8718109239835591], [1.3333333333333333, 0.9351442573168924], [1.6666666666666667, 0.10181092398355907], [1.5833333333333333, 0.20514425731689237], [1.6666666666666667, 0.40876610528401053], [1.5833333333333333, 0.5136218479671181], [1.6666666666666667, 0.6787661052840105], [1.5833333333333333, 0.7451442573168924], [1.6666666666666667, 0.8718109239835591], [1.5833333333333333, 0.9366666666666666], [1.9166666666666667, 0.09747546895706428], [1.8333333333333333, 0.19928639294062336], [1.9166666666666667, 0.394237330854752], [1.8333333333333333, 0.5003837098648057], [1.9166666666666667, 0.6627149215049777], [1.8333333333333333, 0.7362415742410748], [1.9166666666666667, 0.86595305960729], [1.8333333333333333, 0.9351442573168924], [2.1666666666666665, 0.09098700198063514], [2.0833333333333335, 0.1884624709376994], [2.1666666666666665, 0.3704364748989697], [2.0833333333333335, 0.47440041083246315], [2.1666666666666665, 0.6345786105227006], [2.0833333333333335, 0.7167467421851613], [2.1666666666666665, 0.855129137604366], [2.0833333333333335, 0.9308088022903975], [2.4166666666666665, 0.08333333333333333], [2.3333333333333335, 0.17432033531396848], [2.4166666666666665, 0.34098700198063514], [2.3333333333333335, 0.4396276726085721], [2.4166666666666665, 0.5986406706279369], [2.3333333333333335, 0.6896276726085722], [2.4166666666666665, 0.840987001980635], [2.3333333333333335, 0.9243203353139684], [2.6666666666666665, 0.08333333333333333], [2.5833333333333335, 0.16666666666666666], [2.6666666666666665, 0.3333333333333333], [2.5833333333333335, 0.4166666666666667], [2.6666666666666665, 0.5833333333333334], [2.5833333333333335, 0.6666666666666666], [2.6666666666666665, 0.8333333333333334], [2.5833333333333335, 0.9166666666666666], [2.9166666666666665, 0.08333333333333333], [2.8333333333333335, 0.16666666666666666], [2.9166666666666665, 0.3333333333333333], [2.8333333333333335, 0.4166666666666667], [2.9166666666666665, 0.5833333333333334], [2.8333333333333335, 0.6666666666666666], [2.9166666666666665, 0.8333333333333334], [2.8333333333333335, 0.9166666666666666]], "cent_b": [[0.16666666666666666, 0.08333333333333333], [0.08333333333333333, 0.16666666666666666], [0.16666666666666666, 0.3333333333333333], [0.08333333333333333, 0.4166666666666667], [0.16666666666666666, 0.5833333333333334], [0.08333333333333333, 0.6666666666666666], [0.16666666666666666, 0.8333333333333334], [0.08333333333333333, 0.9166666666666666], [0.4166666666666667, 0.08333333333333333], [0.3333333333333333, 0.16666666666666666], [0.4166666666666667, 0.3333333333333333], [0.3333333333333333, 0.4166666666666667], [0.4166666666666667, 0.5833333333333334], [0.3333333333333333, 0.6666666666666666], [0.4166666666666667, 0.8333333333333334], [0.3333333333333333, 0.9166666666666666], [0.6666666666666666, 0.08333333333333333], [0.5833333333333334, 0.16666666666666666], [0.6666666666666666, 0.3333333333333333], [0.5833333333333334, 0.4166666666666667], [0.6666666666666666, 0.5833333333333334], [0.5833333333333334, 0.6666666666666666], [0.6666666666666666, 0.8333333333333334], [0.5833333333333334, 0.9166666666666666], [0.9166666666666666, 0.08333333333333333], [0.8333333333333334, 0.16666666666666666], [0.9166666666666666, 0.3333333333333333], [0.8333333333333334, 0.4166666666666667], [0.9166666666666666, 0.5833333333333334], [0.8333333333333334, 0.6666666666666666], [0.9166666666666666, 0.8333333333333334], [0.8333333333333334, 0.9166666666666666], [1.1666666666666667, 0.08333333333333333], [1.0833333333333333, 0.16666666666666666], [1.1666666666666667, 0.3333333333333333], [1.0833333333333333, 0.4166666666666667], [1.1666666666666667, 0.5833333333333334], [1.0833333333333333, 0.6666666666666666], [1.1666666666666667, 0.8333333333333334], [1.0833333333333333, 0.9166666666666666], [1.4166666666666667, 0.08333333333333333], [1.3333333333333333, 0.16666666666666666], [1.4166666666666667, 0.3333333333333333], [1.3333333333333333, 0.4166666666666667], [1.4166666666666667, 0.5833333333333334], [1.3333333333333333, 0.6666666666666666], [1.4166666666666667, 0.8333333333333334], [1.3333333333333333, 0.9166666666666666], [1.6666666666666667, 0.08333333333333333], [1.5833333333333333, 0.16666666666666666], [1.6666666666666667, 0.3333333333333333], [1.5833333333333333, 0.4166666666666667], [1.6666666666666667, 0.5833333333333334], [1.5833333333333333, 0.6666666666666666], [1.6666666666666667, 0.8333333333333334], [1.5833333333333333, 0.9166666666666666], [1.9166666666666667, 0.08333333333333333], [1.8333333333333333, 0.16666666666666666], [1.9166666666666667, 0.3333333333333333], [1.8333333333333333, 0.4166666666666667], [1.9166666666666667, 0.5833333333333334], [1.8333333333333333, 0.6666666666666666], [1.9166666666666667, 0.8333333333333334], [1.8333333333333333, 0.9166666666666666], [2.1666666666666665, 0.08333333333333333], [2.0833333333333335, 0.16666666666666666], [2.1666666666666665, 0.3333333333333333], [2.0833333333333335, 0.4166666666666667], [2.1666666666666665, 0.5833333333333334], [2.0833333333333335, 0.6666666666666666], [2.1666666666666665, 0.8333333333333334], [2.0833333333333335, 0.9166666666666666], [2.4166666666666665, 0.08333333333333333], [2.3333333333333335, 0.16666666666666666], [2.4166666666666665, 0.3333333333333333], [2.3333333333333335, 0.4166666666666667], [2.4166666666666665, 0.5833333333333334], [2.3333333333333335, 0.6666666666666666], [2.4166666666666665, 0.8333333333333334], [2.3333333333333335, 0.9166666666666666], [2.6666666666666665, 0.08333333333333333], [2.5833333333333335, 0.16666666666666666], [2.6666666666666665, 0.3333333333333333], [2.5833333333333335, 0.4166666666666667], [2.6666666666666665, 0.5833333333333334], [2.5833333333333335, 0.6666666666666666], [2.6666666666666665, 0.8333333333333334], [2.5833333333333335, 0.9166666666666666], [2.9166666666666665, 0.08333333333333333], [2.8333333333333335, 0.16666666666666666], [2.9166666666666665, 0.3333333333333333], [2.8333333333333335, 0.4166666666666667], [2.9166666666666665, 0.5833333333333334], [2.8333333333333335, 0.6666666666666666], [2.9166666666666665, 0.8333333333333334], [2.8333333333333335, 0.9166666666666666]], "chain": [12, 17, 22, 27, 32, 37, 42, 47, 52], "band_a": [18, 19, 20, 21, 26, 27, 28, 29, 34, 35, 36, 37, 42, 43, 44, 45, 50, 51, 52, 53, 58, 59, 60, 61, 66, 67, 68, 69, 74, 75, 76, 77], "dir_a": {"18": [-0.1791278488903745, 0.9838258045771656], "19": [-0.1837025565064759, 0.9829818771131975], "20": [-0.1791278488903745, 0.9838258045771656], "21": [-0.1837025565064759, 0.9829818771131975], "26": [-0.1478989721424039, 0.9890024742331136], "27": [-0.16110943939460595, 0.9869365473716917], "28": [-0.1478989721424039, 0.9890024742331136], "29": [-0.16110943939460595, 0.9869365473716917], "34": [-0.09383196299423324, 0.9955880487032018], "35": [-0.11400073992727823, 0.9934806647821752], "36": [-0.09383196299423324, 0.9955880487032018], "37": [-0.11400073992727823, 0.9934806647821752], "42": [-0.024596164230675614, 0.9996974685899418], "43": [-0.04872828607559788, 0.998812071480984], "44": [-0.024596164230675614, 0.9996974685899418], "45": [-0.04872828607559788, 0.998812071480984], "50": [0.04872828607559786, 0.998812071480984], "51": [0.02459616423067559, 0.9996974685899418], "52": [0.04872828607559786, 0.998812071480984], "53": [0.02459616423067559, 0.9996974685899418], "58": [0.1140007399272782, 0.9934806647821752], "59": [0.09383196299423324, 0.9955880487032018], "60": [0.1140007399272782, 0.9934806647821752], "61": [0.09383196299423324, 0.9955880487032018], "66": [0.16110943939460592, 0.9869365473716917], "67": [0.1478989721424039, 0.9890024742331136], "68": [0.16110943939460592, 0.9869365473716917], "69": [0.1478989721424039, 0.9890024742331136], "74": [0.1837025565064759, 0.9829818771131975], "75": [0.1791278488903745, 0.9838258045771656], "76": [0.1837025565064759, 0.9829818771131975], "77": [0.1791278488903745, 0.9838258045771656]}, "band_b": [18, 19, 20, 21, 26, 27, 28, 29, 35, 36, 37, 38, 43, 44, 45, 46, 51, 52, 53, 54, 59, 60, 61, 62, 66, 67, 68, 69, 74, 75, 76, 77], "dir_b": {"18": [-0.18151302229643904, 0.9833885410847598], "19": [-0.18430173394456653, 0.9828697120498862], "20": [-0.17850478287533372, 0.9839390440929916], "21": [-0.18256790844393192, 0.983193245911712], "26": [-0.1540402697713341, 0.9880645704045735], "27": [-0.16432146800919897, 0.9864068405841993], "28": [-0.14731569679664494, 0.98908952348982], "29": [-0.15868252166126315, 0.9873296599004928], "35": [-0.11929359162070788, 0.9928590227208653], "36": [-0.09425335329998778, 0.995548243628458], "37": [-0.1124586006450129, 0.9936564110098447], "38": [-0.08915190230652295, 0.9960180411594602], "43": [-0.05109847801328024, 0.9986936194573021], "44": [-0.02508488507213821, 0.9996853247602054], "45": [-0.04825313496509566, 0.9988351390324834], "46": [-0.02361816787394588, 0.9997210521671922], "51": [0.026549894089562934, 0.9996474894300655], "52": [0.04967748606765285, 0.998765311461105], "53": [0.02361816787394586, 0.9997210521671922], "54": [0.04539473394090614, 0.9989691277163847], "59": [0.09801443556848749, 0.9951849930641995], "60": [0.11476688599598552, 0.9933924510880806], "61": [0.09298712137823809, 0.9956673115342236], "62": [0.10775485881489509, 0.9941774944152486], "66": [0.16627056238462376, 0.9860801691973634], "67": [0.1515892675353472, 0.9884435714637921], "68": [0.16087169181916747, 0.9869753283498219], "69": [0.14554128987838752, 0.9893521784180469], "74": [0.18468204578253408, 0.9827983221218777], "75": [0.18060764653052153, 0.9835552236731329], "76": [0.1832444605666451, 0.9830673769745586], "77": [0.17730813005971677, 0.9841553876368947]}, "curve": [[0.5, 0.5], [0.5333333333333333, 0.5062803147491532], [0.5666666666666667, 0.5125434155921185], [0.6, 0.5187721358048277], [0.6333333333333333, 0.5249494028981311], [0.6666666666666666, 0.5310582854123025], [0.7, 0.5370820393249937], [0.7333333333333334, 0.543004153945436], [0.7666666666666666, 0.548808397169096], [0.8, 0.5544788599687456], [0.8333333333333333, 0.5599999999999999], [0.8666666666666667, 0.5653566842018033], [0.9, 0.5705342302750968], [0.9333333333333333, 0.5755184469259805], [0.9666666666666667, 0.580295672763063], [1.0, 0.5848528137423857], [1.0333333333333332, 0.5891773790572873], [1.0666666666666667, 0.5932575153748365], [1.1, 0.5970820393249937], [1.1333333333333333, 0.6006404681534508], [1.1666666666666665, 0.6039230484541326], [1.2, 0.6069207829026042], [1.2333333333333334, 0.6096254549171121], [1.2666666666666666, 0.6120296511796642], [1.3, 0.6141267819554184], [1.3333333333333335, 0.6159110991546882], [1.3666666666666667, 0.6173777120880567], [1.4, 0.6185226008714165], [1.4333333333333333, 0.6193426274441928], [1.4666666666666668, 0.6198355441705489], [1.5, 0.62], [1.5333333333333334, 0.6198355441705489], [1.5666666666666667, 0.6193426274441928], [1.6, 0.6185226008714165], [1.6333333333333333, 0.6173777120880567], [1.6666666666666667, 0.6159110991546882], [1.7, 0.6141267819554185], [1.7333333333333334, 0.6120296511796642], [1.7666666666666666, 0.6096254549171121], [1.8, 0.6069207829026042], [1.8333333333333333, 0.6039230484541327], [1.8666666666666667, 0.6006404681534508], [1.9, 0.5970820393249937], [1.9333333333333333, 0.5932575153748365], [1.9666666666666666, 0.5891773790572873], [2.0, 0.5848528137423857], [2.033333333333333, 0.580295672763063], [2.0666666666666664, 0.5755184469259805], [2.1, 0.5705342302750968], [2.1333333333333333, 0.5653566842018033], [2.166666666666667, 0.5599999999999999], [2.2, 0.5544788599687456], [2.2333333333333334, 0.548808397169096], [2.2666666666666666, 0.5430041539454361], [2.3, 0.5370820393249938], [2.333333333333333, 0.5310582854123026], [2.3666666666666667, 0.5249494028981311], [2.4, 0.5187721358048277], [2.4333333333333336, 0.5125434155921184], [2.466666666666667, 0.5062803147491532], [2.5, 0.5]], "width": 0.5} \ No newline at end of file diff --git a/articles/putting-a-fault-in-a-mesh/examples/fault-anatomy.typ b/articles/putting-a-fault-in-a-mesh/examples/fault-anatomy.typ new file mode 100644 index 0000000..c768728 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/fault-anatomy.typ @@ -0,0 +1,119 @@ +// Three representations of one fault on one mesh (cetz draws; the +// geometry comes from generate-band-anatomy-data.py and +// generate-split-anatomy-data.py over the shared anatomy_mesh.py). +// (a) the non-conforming paint, the trace across the flat grid +// (b) the ribbon, on the mesh bent so the trace runs along its edges +// (c) the split, on the bent mesh, exploded for display +// No captions in the drawing: they live in the figure caption. +#import "@preview/cetz:0.3.4" + +#set page(width: auto, height: auto, margin: 10pt) +#set text(font: ("Noto Sans", "Helvetica", "Arial"), size: 5pt) + +#let band = json("band-anatomy-data.json") +#let split = json("split-anatomy-data.json") + +#let mesh-fill = rgb("#f4f6f9") +#let band-fill = rgb("#c9e6c9") +#let fill-plus = rgb("#dce8fc") +#let fill-minus = rgb("#fce4ec") +#let mesh-stroke = rgb("#8899aa") +#let fault-col = rgb("#c62828") +#let dir-col = rgb("#1b5e20") +#let tip-col = rgb("#1a1a1a") +#let note-col = rgb("#555555") + +#let cells(oy, coords, tris, fill-of) = { + import cetz.draw: * + let P(v) = (coords.at(v).at(0), coords.at(v).at(1) + oy) + for (k, t) in tris.enumerate() { + line(P(t.at(0)), P(t.at(1)), P(t.at(2)), close: true, + fill: fill-of(k), stroke: (paint: mesh-stroke, thickness: 0.3pt)) + } +} + +#let directors(oy, cent, keys, dirs) = { + import cetz.draw: * + for k in keys { + let c = (cent.at(k).at(0), cent.at(k).at(1) + oy) + let d = dirs.at(str(k)) + let h = 0.055 + line((c.at(0) - h * d.at(0), c.at(1) - h * d.at(1)), + (c.at(0) + h * d.at(0), c.at(1) + h * d.at(1)), + stroke: (paint: dir-col, thickness: 0.7pt)) + } +} + +#let GAP = 1.35 +#let YA = GAP +#let YB = 2 * GAP +#let YC = 0.0 + +#cetz.canvas(length: 1.55cm, { + import cetz.draw: * + + // ---- (b) the ribbon on the bent mesh ------------------------------------- + cells(YA, band.bent, band.tris, + k => if band.band_a.contains(k) { band-fill } else { mesh-fill }) + directors(YA, band.cent_a, band.band_a, band.dir_a) + let Pa(v) = (band.bent.at(v).at(0), band.bent.at(v).at(1) + YA) + for i in range(band.chain.len() - 1) { + line(Pa(band.chain.at(i)), Pa(band.chain.at(i + 1)), + stroke: (paint: fault-col, thickness: 0.9pt)) + } + for v in band.chain { + circle(Pa(v), radius: 0.028, fill: fault-col, stroke: none) + } + content((-0.3, YA + 0.95), [(b)]) + content((3.5, YA + 0.5), text(fill: fault-col)[$Gamma$]) + line((3.12, YA + 0.25), (3.12, YA + 0.75), stroke: (paint: dir-col, thickness: 0.5pt)) + line((3.07, YA + 0.25), (3.17, YA + 0.25), stroke: (paint: dir-col, thickness: 0.5pt)) + line((3.07, YA + 0.75), (3.17, YA + 0.75), stroke: (paint: dir-col, thickness: 0.5pt)) + content((3.28, YA + 0.5), text(fill: dir-col, size: 4.5pt)[$w$]) + + // ---- (a) the same trace on the flat grid --------------------------------- + cells(YB, band.flat, band.tris, + k => if band.band_b.contains(k) { band-fill } else { mesh-fill }) + directors(YB, band.cent_b, band.band_b, band.dir_b) + for i in range(band.curve.len() - 1) { + line((band.curve.at(i).at(0), band.curve.at(i).at(1) + YB), + (band.curve.at(i + 1).at(0), band.curve.at(i + 1).at(1) + YB), + stroke: (paint: fault-col, thickness: 0.9pt)) + } + content((-0.3, YB + 0.95), [(a)]) + content((3.5, YB + 0.5), text(fill: fault-col)[$Gamma$]) + + // ---- (c) the split on the bent mesh, exploded ---------------------------- + let touches(t) = t.any(v => split.chain.contains(v)) + cells(YC, split.exploded, split.moved_tris, + k => if split.side.at(k) < 0 { fill-minus } + else if touches(split.tris.at(k)) { fill-plus } else { mesh-fill }) + let Pc(v) = (split.exploded.at(v).at(0), split.exploded.at(v).at(1) + YC) + let lower(v) = if str(v) in split.replicas { split.replicas.at(str(v)) } else { v } + for i in range(split.chain.len() - 1) { + line(Pc(split.chain.at(i)), Pc(split.chain.at(i + 1)), + stroke: (paint: fault-col, thickness: 0.9pt)) + line(Pc(lower(split.chain.at(i))), Pc(lower(split.chain.at(i + 1))), + stroke: (paint: fault-col, thickness: 0.9pt, dash: "densely-dashed")) + } + for v in split.interior { + circle(Pc(v), radius: 0.028, fill: fault-col, stroke: none) + } + for (orig, rep) in split.replicas { + circle(Pc(rep), radius: 0.028, fill: white, + stroke: (paint: fault-col, thickness: 0.8pt)) + } + for v in split.tips { + circle(Pc(v), radius: 0.04, fill: white, stroke: (paint: tip-col, thickness: 0.9pt)) + circle(Pc(v), radius: 0.014, fill: tip-col, stroke: none) + } + content((-0.3, YC + 0.95), [(c)]) + content((0.42, YC + 0.38), text(size: 4pt)[tip]) + content((2.58, YC + 0.38), text(size: 4pt)[tip]) + content((3.5, YC + 0.66), text(fill: fault-col, size: 4.5pt)[$Gamma^+$]) + content((3.5, YC + 0.34), text(fill: fault-col, size: 4.5pt)[$Gamma^-$]) + line((1.5, YC + 0.66), (1.85, YC + 0.95), stroke: (paint: note-col, thickness: 0.3pt)) + content((2.25, YC + 1.02), text(size: 4pt)[$v^+$ (original)]) + line((1.5, YC + 0.49), (1.85, YC + 0.12), stroke: (paint: note-col, thickness: 0.3pt)) + content((2.25, YC + 0.05), text(size: 4pt)[$v^-$ (replica)]) +}) diff --git a/articles/putting-a-fault-in-a-mesh/examples/generate-band-anatomy-data.py b/articles/putting-a-fault-in-a-mesh/examples/generate-band-anatomy-data.py new file mode 100644 index 0000000..08a1a45 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/generate-band-anatomy-data.py @@ -0,0 +1,42 @@ +"""Geometry for the band-anatomy figure (cetz draws, Python computes). + +(a) the ribbon on the CONFORMING mesh: the trace along the bent row of + edges, the band the cells either side of it, each with a director + (the trace normal); no vertex duplicated, every field continuous. +(b) the non-conforming paint on the FLAT grid: the same curve and the + same width, but the band is whichever cells fall within w/2 of it. +""" +import json +import os + +import anatomy_mesh as M + +coords, tris, vid = M.grid() +bent = M.conforming(coords) +ch = M.chain(vid) + +# (a) the ribbon: the cells that touch the chain, over its span +cent_a = M.centroids(bent, tris) +band_a = [k for k, t in enumerate(tris) + if any(v in ch for v in t) and M.X0 < cent_a[k][0] < M.X1] +dir_a = {k: M.normal(cent_a[k][0]) for k in band_a} + +# (b) the paint: flat grid, cells within w/2 of the curve +cent_b = M.centroids(coords, tris) +band_b, dir_b = [], {} +for k, c in enumerate(cent_b): + d, x = M.distance_to_curve(c) + if d <= 0.5 * M.W and M.X0 <= c[0] <= M.X1: + band_b.append(k) + dir_b[k] = M.normal(x) + +samples = [M.curve(M.X0 + (M.X1 - M.X0) * k / 60) for k in range(61)] +out = dict(flat=coords, bent=bent, tris=tris, cent_a=cent_a, cent_b=cent_b, + chain=ch, band_a=band_a, dir_a={str(k): v for k, v in dir_a.items()}, + band_b=band_b, dir_b={str(k): v for k, v in dir_b.items()}, + curve=samples, width=M.W) +here = os.path.dirname(os.path.abspath(__file__)) +with open(os.path.join(here, "band-anatomy-data.json"), "w") as f: + json.dump(out, f) +print(f"wrote band-anatomy-data.json: band (a) {len(band_a)} cells, " + f"band (b) {len(band_b)} cells") diff --git a/articles/putting-a-fault-in-a-mesh/examples/generate-split-anatomy-data.py b/articles/putting-a-fault-in-a-mesh/examples/generate-split-anatomy-data.py new file mode 100644 index 0000000..9dbd23e --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/generate-split-anatomy-data.py @@ -0,0 +1,53 @@ +"""Geometry for the split-anatomy figure (cetz draws, Python computes). + +The same CONFORMING mesh as band-anatomy (a): the curved trace runs +along a row of edges. The chain has nine vertices: two tips (shared) +and seven interior (duplicated by the split). The "after" panel is +EXPLODED for display — the Minus block under the trace translated down +— because the real copies are geometrically coincident. +""" +import json +import os + +import anatomy_mesh as M + +DELTA = 0.09 # display-only explosion offset + +coords, tris, vid = M.grid() +bent = M.conforming(coords) +ch = M.chain(vid) +tips = [ch[0], ch[-1]] +interior = ch[1:-1] +cent = M.centroids(bent, tris) + +# Minus side = cells touching the chain from below, within the span; +# everything else, including the cells beyond the tips, stays welded +side = [] +for k, t in enumerate(tris): + touches = any(v in ch for v in t) + below = cent[k][1] < M.curve(cent[k][0])[1] + side.append(-1 if (below and touches and M.X0 < cent[k][0] < M.X1) else +1) + +# Exploded: the replicas and the whole Minus block strictly inside the +# span move down together; the tips and the columns at the tips stay, +# so the cells at the two ends shear — the pinned tip +exploded = [list(c) for c in bent] +replicas = {} +for v in interior: + replicas[v] = len(exploded) + exploded.append([bent[v][0], bent[v][1] - DELTA]) +for v, (x, y) in enumerate(coords): + if y < 0.5 - 1e-9 and M.X0 + 1e-9 < x < M.X1 - 1e-9: + exploded[v][1] -= DELTA +moved_tris = [[replicas.get(v, v) for v in t] if s < 0 else list(t) + for t, s in zip(tris, side)] + +out = dict(coords=bent, exploded=exploded, tris=tris, moved_tris=moved_tris, + side=side, chain=ch, tips=tips, interior=interior, + replicas=replicas, delta=DELTA) +here = os.path.dirname(os.path.abspath(__file__)) +with open(os.path.join(here, "split-anatomy-data.json"), "w") as f: + json.dump(out, f) +print("wrote split-anatomy-data.json:", + f"{len(bent)} verts, {len(tris)} tris, chain {len(ch)},", + f"replicas {len(replicas)}") diff --git a/articles/putting-a-fault-in-a-mesh/examples/rig_geometry.py b/articles/putting-a-fault-in-a-mesh/examples/rig_geometry.py new file mode 100644 index 0000000..b1d5630 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/rig_geometry.py @@ -0,0 +1,340 @@ +"""The S-fault rig's geometry of record — ONE source for every script. + +Adjust the dials here; geometry_preview.py renders them, s_fault.py +solves them, sf_viz.py / sf_compare_viz.py draw them. + +Strands: + Main — the tanh S on the 40° diagonal (the San Andreas analogue); + also the TERRANE boundary (strong above/left of its line). + Branch — the through-line alternative in the WEAK block (tip-to-tip + near-join with the bend; the San Jacinto analogue). + Splay — the Y near-join in the TOUGH layer: a short strand leaving + the outgoing limb's side at SPLAY_PHI into the strong block, + stopped short (tip-to-SIDE junction). The trace gap is the + junction physics (the intact-gap linkage); it KISSES the main + (gap w/2) so the mesh surround must interact — placed as one + fused network with the main line (see SPLAY_GAP note). In TI + the two weak zones merge at the Y (one director per cell: + paint order decides — the #544 caveat); split keeps two cuts + and an intact sliver. + Cont — a near-touching SEGMENT of the main fault along its + continuation line beyond the upper tip (tip-to-tip stepover, + the California fault-database segmentation pattern). Same + near-touching gap as the Y; same test from the other + direction: collinear instead of oblique. +""" +import numpy as np + +# ----------------------------------------------------------------- dials +THETA = np.deg2rad(40.0) # strike of the diagonal baseline +C = np.array([0.5, 0.5]) # baseline centre +A = 0.05 # S half-step +LAM = 0.10 # bend length scale (tanh) +S_MAIN = (-0.42, 0.42) # main-trace arc range (blind tips) +S_BRANCH = (0.06, 0.42) # branch: through-line at t = -A +SPLAY_S0 = 0.06 # splay leaves the main line near this s — + # just past the bend's centre, so at +22 deg it + # CONTINUES the bend's turning (Louis, final + # pass 2026-08-26: "more of a continuation of + # the curved segment; shift it left") +SPLAY_GAP_RUNGS = 1.0 # trace-to-trace gap at the Y in RUNGS (s_f): + # KISSING = one rung at any resolution + # (a fixed 0.015 was w/2 at coarse but 1.5 w + # at fine). traces(s_f) sets SPLAY_GAP. +SPLAY_GAP = 0.010 # (derived; the coarse value) KISSING (w/2 + # at coarse; Louis's geometry of record). + # Any gap: the whole rig is ONE network + # placement (fused ribbons + embedded + # spines, place_thin_volume mesher= + # "network"); junctions are free and the + # cuts walk exact vertices. +SPLAY_PHI = np.deg2rad(22.0) # splay divergence from local strike, +t side +SPLAY_LEN = 0.18 # splay length +EE_N = 3 # EN-ECHELON: the splay + 2 more parallel +EE_LATERAL = 0.08 # segments: in the SPLAY's own frame, each +EE_ALONG = 0.05 # start = previous + EE_LATERAL n_sp (across + # the strands, up-left into the strong block) + # + EE_ALONG d_sp (forward along the strand). + # Louis 2026-08-26 (3rd pass): the array + # climbs UP as it goes LEFT; spacing 0.058; + # only #1 kisses the main +STEP_S = (-0.53, -0.37) # STEPOVER at the main's LOWER (clean) tip: a +STEP_OFF = 0.035 # segment parallel to the lower limb, offset + # STEP_OFF (ABSOLUTE — the same geometry at + # both resolutions; Louis: "close, hard to + # separate") into the WEAK block, overlapping + # the tip (s=-0.42) by 0.05 and running 0.11 + # toward the wall (wall clearance limits it) +CONT_GAP_RUNGS = 1.0 # tip-to-tip stepover gap beyond the main's + # upper tip, in RUNGS (s_f), as the splay's: + # 1 = the split's floor (chains must be + # disjoint, so one uncut edge is the least + # it can leave); 0 = the CONTINUOUS control, + # one cut the whole line. traces(s_f) sets + # CONT_GAP. It was an absolute 0.015, which + # refinement walked AWAY from the floor (2 + # edges coarse, 6 superfine) — the test the + # segment exists for is that the band closes + # a gap the cut cannot. None = leave CONT_GAP + # as set (the absolute dial of older probes). +CONT_GAP = 0.015 # (derived; the coarse value) tip-to-tip + # stepover, near-touching, California-DB + # style. Placeable at ANY value: the + # Main-line family shares ONE band + # (main_line_union), so the gap lives only + # in the cuts + paint. +CONT_LEN = 0.10 # continuation-segment length along the line + # (0.18 put the extended band's cavity into + # the top-right wall — interior clearance) +# ----------------------------------------------------------------------- + +ES = np.array([np.cos(THETA), np.sin(THETA)]) +ET = np.array([-np.sin(THETA), np.cos(THETA)]) +RES = {"coarse": (0.03, 0.01), "fine": (0.01, 0.005), + "superfine": (0.005, 0.0025)} # (w, s_f) +W_OF_SF = {sf: w for w, sf in RES.values()} # s_f -> w + + +def main_curve(s): + s = np.atleast_1d(s) + return C + s[:, None] * ES + (A * np.tanh(s / LAM))[:, None] * ET + + +def resample_arclength(P, spacing): + """Equispaced-in-arclength, ends kept, interval count EVEN.""" + seg = np.linalg.norm(np.diff(P, axis=0), axis=1) + arc = np.concatenate([[0.0], np.cumsum(seg)]) + n = max(2, int(round(arc[-1] / spacing))) + n += n % 2 + si = np.linspace(0.0, arc[-1], n + 1) + return np.column_stack([np.interp(si, arc, P[:, 0]), + np.interp(si, arc, P[:, 1])]) + + +def splay_line(k=0): + """En-echelon strand k (k=0 is THE splay, kissing the main): from a + stop-short start beside the outgoing limb, diverging into the + strong block; strand k starts EE_STEP back along strike and out.""" + t0 = float(A * np.tanh(SPLAY_S0 / LAM)) + # local strike of the main line at SPLAY_S0 + sech2 = 1.0 / np.cosh(SPLAY_S0 / LAM) ** 2 + tan_l = ES + (A / LAM) * sech2 * ET + tan_l /= np.linalg.norm(tan_l) + nrm_l = np.array([-tan_l[1], tan_l[0]]) # +t side + if nrm_l @ ET < 0: + nrm_l = -nrm_l + d = np.cos(SPLAY_PHI) * tan_l + np.sin(SPLAY_PHI) * nrm_l + n_sp = np.array([-d[1], d[0]]) # up-left of the strand + start = (C + SPLAY_S0 * ES + t0 * ET + SPLAY_GAP * nrm_l + + k * (EE_LATERAL * n_sp + EE_ALONG * d)) + return np.array([start, start + SPLAY_LEN * d]) + + +def stepover_line(w=None): + """The stepover at the main's lower tip: parallel to the lower limb + (t = -A there), offset STEP_OFF (absolute) into the weak block.""" + ss = np.linspace(*STEP_S, 500) + return C + ss[:, None] * ES - (A + STEP_OFF) * ET + + +def main_line_union(s_f): + """The Main-line family as ONE parametrisation — the library's rule + for close pairs ('surfaces must be separated by at least a cell; + place close pairs as one thin volume instead'): one ladder polyline + over the WHOLE line, with the Main and Cont cut chains as vertex + SUBCHAINS of it (#595: nothing snaps — cuts consume the spine's own + vertices, so chains and band must share them). The stepover gap + therefore lives only in the cuts and the paint, and can close below + the placement floor. Returns (union, main_chain, cont_chain); + chain ends are vertex-quantised, so the realised gap is the nearest + ≥ CONT_GAP the sampling allows.""" + s_end = S_MAIN[1] + CONT_GAP + CONT_LEN + sm = np.linspace(S_MAIN[0], s_end, 4000) + union = resample_arclength(main_curve(sm), s_f) + if CONT_GAP <= 0.0: + # the CONTINUOUS control: no stepover — one cut, the whole line + return union, union, union[:0] + s_along = (union - C) @ ES + # the gap in SPINE EDGES (distance thresholds are at the mercy of + # vertex phase): n_edges=1 -> abutting chains, one uncut edge; + # n_edges=2 -> one free bridge vertex; realised gap = n_edges rungs + i_tip = int(np.searchsorted(s_along, S_MAIN[1] + 1e-9)) - 1 + rung = float(np.linalg.norm(np.diff(union, axis=0), axis=1).mean()) + n_edges = max(1, int(round(CONT_GAP / rung))) + return union, union[:i_tip + 1], union[i_tip + n_edges:] + + +def traces(s_f, w=None): + """The rig's strands, sampled at the rung scale — [(label, points)]: + Main, Branch, Splay (+ Splay2, Splay3 en echelon), Step (the lower + stepover), Cont (the upper stepover segment, if the gap is open). + Main and Cont are subchains of the shared main_line_union spine.""" + global SPLAY_GAP, CONT_GAP + SPLAY_GAP = SPLAY_GAP_RUNGS * s_f + if CONT_GAP_RUNGS is not None: + CONT_GAP = CONT_GAP_RUNGS * s_f + w = W_OF_SF[s_f] if w is None else w + _union, main, cont = main_line_union(s_f) + sb = np.linspace(*S_BRANCH, 2000) + branch = resample_arclength(C + sb[:, None] * ES - A * ET, s_f) + out = [("Main", main), ("Branch", branch)] + for k in range(EE_N): + sp = splay_line(k) + seg = resample_arclength( + np.linspace(0, 1, 200)[:, None] * (sp[1] - sp[0]) + sp[0], s_f) + out.append(("Splay" if k == 0 else f"Splay{k + 1}", seg)) + out.append(("Step", resample_arclength(stepover_line(w), s_f))) + if len(cont): + out.append(("Cont", cont)) + return out + + +def strong_mask(points): + """The terrane rule: True on the +t (strong) side of the fault line.""" + rel = np.asarray(points) - C + return (rel @ ET) > A * np.tanh((rel @ ES) / LAM) + + +def receiver_normals(points): + """Local fault-line orientation, continued off the line (for dCFF).""" + rel = np.asarray(points) - C + s_co = rel @ ES + sech2 = 1.0 / np.cosh(s_co / LAM) ** 2 + n = ET[None, :] - (A / LAM) * sech2[:, None] * ES[None, :] + return n / np.linalg.norm(n, axis=1, keepdims=True) + + +def ti_fields(mesh, foots, eta1_val, eta0_vals, tag, w=None, traces_=None): + """The TI representation's P0 fields, painted on the wrapper's + HONOURED footprints (the review-1 API): eta_1 = weak inside the + fault footprints only, background elsewhere; per-cell directors + following each strand's own orientation (curved Main and its Cont + segment, ET-normal Branch, splay-normal Splay). Returns + (eta1_var, dir_var, foot). + + With ``w`` and ``traces_`` given, the stepover corridor between the + Main's tip and the Cont's start is painted too: a footprint stops + half a rung past its tip, so a one-edge gap leaves a plug of a few + intact cells in the band (4 at superfine) and the TI still has a + gap. Two collinear segments of ONE line closer than the band is + wide are one weak zone — that is the band's whole claim on the + stepover, and the paint has to say it.""" + import underworld3 as uw + + foot = np.zeros_like(next(iter(foots.values()))) + for m_ in foots.values(): + foot = foot | m_ + eta1 = uw.discretisation.MeshVariable(f"etaT{tag}", mesh, 1, degree=0) + cen_abs = np.asarray(eta1.coords) + if w is not None and traces_ is not None: + tr = dict(traces_) + if "Main" in tr and "Cont" in tr: + a, b = np.asarray(tr["Main"][-1]), np.asarray(tr["Cont"][0]) + d = b - a + L = float(np.linalg.norm(d)) + d = d / L + rel = cen_abs - a + along = rel @ d + across = np.abs(rel @ np.array([-d[1], d[0]])) + corridor = (along > -0.5 * L) & (along < 1.5 * L) & (across < 0.5 * w) + n_new = int((corridor & ~foot).sum()) + foot = foot | corridor + foots = dict(foots) + foots["Main"] = foots["Main"] | corridor + uw.pprint(f"[ti_fields] stepover corridor: {n_new} band cell(s) " + f"painted between the Main's tip and the Cont's start") + eta1.array[:, 0, 0] = np.where(foot, eta1_val, eta0_vals) + ndir = uw.discretisation.MeshVariable(f"dirT{tag}", mesh, 2, degree=0, + continuous=False) + cen = cen_abs - C + s_co = cen @ ES + dvals = np.tile(ET, (len(s_co), 1)) + sech2 = 1.0 / np.cosh(s_co / LAM) ** 2 + n_main = ET[None, :] - (A / LAM) * sech2[:, None] * ES[None, :] + n_main /= np.linalg.norm(n_main, axis=1, keepdims=True) + if "Main" in foots: + dvals[foots["Main"]] = n_main[foots["Main"]] + if "Cont" in foots: + dvals[foots["Cont"]] = n_main[foots["Cont"]] + sp = splay_line() + d_sp = (sp[1] - sp[0]) / np.linalg.norm(sp[1] - sp[0]) + for lbl, m_ in foots.items(): + if lbl.startswith("Splay"): # all en-echelon strands + dvals[m_] = np.array([-d_sp[1], d_sp[0]]) + # "Step" and "Branch": strike-parallel -> the ET default + ndir.array[...] = dvals.reshape(ndir.array.shape) + return eta1, ndir, foot + + +# The drive is built by callers from mesh coordinates: +# t_sym = (x - 0.5) * ET[0] + (y - 0.5) * ET[1] +# v = sense * 2 * t_sym * ES (plate-motion-frame shear) + + +def place_rig(base_mesh, s_f, width, margin_rings=2, split=True, uncut=()): + """ONE placement path for the whole rig (Louis, 2026-08-25: no more + parallel paths): every strand's extended spine goes into a single + ``place_thin_volume(..., mesher="network")`` call — the ribbons are + fused in CAD (junctions free, so the Y may kiss and the stepover + shares its band) and every spine is EMBEDDED in the fused face, so + the cuts walk exact vertices (#595) at any resolution. Cuts: ONE + network add_fault call. Honoured footprints per strand by nearest + sample over the concatenated spines. Returns (mesh, info).""" + import underworld3 as uw + from underworld3.utilities.place_surface import ( + place_thin_volume, _extend_polyline_2d, _footprint_from_samples) + + tr = traces(s_f, width) + union, chain_main, chain_cont = main_line_union(s_f) + m = margin_rings + spines = [("MainLine", union)] + [(lbl, P) for lbl, P in tr + if lbl not in ("Main", "Cont")] + ext = [_extend_polyline_2d(np.asarray(P, dtype=float), m) + for _n, P in spines] + spacing = float(np.linalg.norm(np.diff(union, axis=0), axis=1).mean()) + dm, _ = place_thin_volume(base_mesh.dm, ext, width, label="Band", + label_value=71, clearance=0.3, size=spacing, + mesher="network") + mesh = uw.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 MG tail + the band as FAC zone: every + # solver on it (Stokes, rotated, projections, probes, renders) drives + # FMG/FAC automatically — no per-solver set_custom_fmg, no silent GAMG + from underworld3.utilities.custom_mg import adopt_hierarchy + adopt_hierarchy(mesh, base_mesh, + fac_zone=mesh.cells_labelled("Band", 71)) + if split: + # ``uncut``: strands left WITHOUT a cut (band only — the fused + # control: no jump, no weld, by construction) + mesh = mesh.add_fault([(lbl, np.asarray(P, dtype=float)) + for lbl, P in tr if lbl not in uncut]) + mesh._custom_mg_fac_zone = None # a split fault needs no patch + band = mesh.cells_labelled("Band", 71) + S_all = np.vstack(ext) + off = np.cumsum([0] + [len(S) for S in ext]) # spine offsets + def user(k, lo, hi): + u = np.zeros(len(S_all), dtype=bool) + u[off[k] + lo:off[k] + hi] = True + return u + foots = {"Main": _footprint_from_samples( + mesh.dm, band, S_all, user(0, m, m + len(chain_main)))} + if len(chain_cont): + foots["Cont"] = _footprint_from_samples( + mesh.dm, band, S_all, + user(0, m + len(union) - len(chain_cont), m + len(union))) + for k, (lbl, _P) in enumerate(spines[1:], start=1): + foots[lbl] = _footprint_from_samples( + mesh.dm, band, S_all, user(k, m, len(ext[k]) - m)) + # the extended spines (what the RIBBON covers) and the cut chains on + # each (what the SPLIT sliced): their difference along the spine is + # where the surface representation cannot join (Louis, 2026-08-27) + cut_chains = {"MainLine": [c for c in (chain_main, chain_cont) if len(c)]} + for lbl, P in spines[1:]: + cut_chains[lbl] = [np.asarray(P, dtype=float)] + return mesh, {"n_cells": int(mesh.dm.getHeightStratum(0)[1]), + "n_rungs": [len(P) for _n, P in spines], + "footprints": foots, "band": band, + "spines_ext": [(n_, E) for (n_, _P), E in zip(spines, ext)], + "cut_chains": cut_chains} diff --git a/articles/putting-a-fault-in-a-mesh/examples/sf_fields.py b/articles/putting-a-fault-in-a-mesh/examples/sf_fields.py new file mode 100644 index 0000000..ab24b5f --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/sf_fields.py @@ -0,0 +1,396 @@ +"""Shared IN-SESSION solve + nodal recovery for the rig's field renders +(in-session because of #640: split-mesh checkpoint reload mixes the +cut DOFs). One solver copy for every render script. + +Recovery follows the standing rulings (reference_stress_recovery_across +_a_material_jump; fault-teaching handoff): + * project the stress COMPONENTS to continuous P1 — never an + invariant, never evaluate() a composite sqrt; + * strain rate across a viscosity jump = recovered stress divided by + each material's OWN viscosity (the jump re-enters through eta, not + through a smeared gradient). For TI, in the (n, t) frame: + e_nn, e_tt = tau/(2 eta0); e_nt = tau_nt/(2 eta1). + * invariants in numpy. +Nodal eta / director are taken from the NEAREST cell centroid (the +honoured paint, cell-based; the one-ring at the band edge is the only +ambiguity). +""" +import numpy as np +import sympy +from scipy.spatial import cKDTree + +import underworld3 as uw +from underworld3.utilities import fault_contact + +import rig_geometry as G + + +def solve_rep(rep, contrast=1.0, res="coarse", cont_gap=0.004, + eta1_val=1e-3, sense=1.0, tol=1e-5, tag="", yield_stress=0.0, + eta1_frac=0.1, uncut=(), uncut_eta1=-1.0): + """Solve the rig for one representation; returns the live state. + ``yield_stress > 0`` confines a von Mises yield to the Band label + (split: cuts + yielding band; ``rep="vp"``: the isotropic yielding + band IS the fault — no cuts, no director). ``rep="hybrid"`` = the + recipe (fault_interface_equivalence/HYBRID_RECIPE.md): frictionless + cuts PLUS a mild TI band, eta_1 = eta1_frac x eta_0 on the + footprints.""" + G.CONT_GAP_RUNGS = None # the absolute dial + G.CONT_GAP = float(cont_gap) + W, S_F = G.RES[res] + TRACES = G.traces(S_F) + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1 / 8, + regular=False, qdegree=2, refinement=1) + mesh, pinfo = G.place_rig(base, S_F, W, + split=(rep in ("split", "hybrid")), + uncut=tuple(uncut)) + x, y = mesh.X + t = f"{rep[0]}{tag}" + v = uw.discretisation.MeshVariable(f"v{t}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"p{t}", mesh, 1, degree=1, + continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + eta0 = uw.discretisation.MeshVariable(f"eta0{t}", mesh, 1, degree=0) + cen = np.asarray(eta0.coords) + strong = G.strong_mask(cen) + eta0_vals = np.where(strong, float(contrast), 1.0) + eta0.array[:, 0, 0] = eta0_vals + st = dict(rep=rep, contrast=float(contrast), mesh=mesh, v=v, p=p, + stokes=stokes, eta0=eta0, eta0_vals=eta0_vals, pinfo=pinfo, + traces=TRACES, W=W, s_f=S_F, eta1=None, ndir=None, + foot=None) + if rep in ("ti", "hybrid") or uncut: + eta1, ndir, foot = G.ti_fields(mesh, pinfo["footprints"], + float(eta1_val), eta0_vals, t) + if rep == "hybrid": + eta1.array[:, 0, 0] = np.where(foot, float(eta1_frac) * eta0_vals, + eta0_vals) + elif rep == "split": + eta1.array[:, 0, 0] = eta0_vals + if uncut and float(uncut_eta1) >= 0.0: + ufoot = np.zeros_like(foot) + for u in uncut: + ufoot |= pinfo["footprints"][u] + eta1.array[:, 0, 0] = np.where(ufoot, float(uncut_eta1) * eta0_vals, + eta1.array[:, 0, 0]) + st["uncut"] = tuple(uncut) + stokes.constitutive_model = \ + uw.constitutive_models.TransverseIsotropicFlowModel + P = stokes.constitutive_model.Parameters + P.shear_viscosity_0 = eta0.sym[0] + P.shear_viscosity_1 = eta1.sym[0] + P.director = ndir.sym + st.update(eta1=eta1, ndir=ndir, foot=foot) + elif yield_stress > 0.0: + ybar = uw.discretisation.MeshVariable(f"y{t}", mesh, 1, degree=0) + ybar.array[:, 0, 0] = np.where(pinfo["band"], float(yield_stress), + 1e8) + stokes.constitutive_model = \ + uw.constitutive_models.ViscoPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta0.sym[0] + stokes.constitutive_model.Parameters.yield_stress = ybar.sym[0] + # Newton on a FIXED rounded yield surface (Louis, 2026-08-26): the + # hard-Min perfect-plastic tangent is 0 along the flow in every + # yielded cell (symmetric, semi-definite) and stalls Newton when + # many cells yield; delta=0.1 (yield_anchor "yield": exact yield + # point, rounded corner) restores it. Not the retired in-solve + # delta-march. + cm_ = stokes.constitutive_model + cm_.yield_mode = "softmin" + cm_.yield_smoother = "powermean" + cm_.yield_anchor = "yield" + cm_.yield_softness = 0.1 + stokes.consistent_jacobian = True + stokes.petsc_options["snes_max_it"] = 60 + st["tau_y"] = float(yield_stress) + else: + if rep == "vp": + raise ValueError("rep='vp' needs yield_stress > 0") + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta0.sym[0] + stokes.bodyforce = [0.0, 0.0] + ES, ET = G.ES, G.ET + t_sym = (x - 0.5) * float(ET[0]) + (y - 0.5) * float(ET[1]) + for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc((sense * 2.0 * t_sym * float(ES[0]), + sense * 2.0 * t_sym * float(ES[1])), wall) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = float(tol) + if rep in ("split", "hybrid"): + for lbl, _tr in TRACES: + if lbl in uncut: + continue + stokes.add_fault_bc(0, boundary=lbl) + info = fault_contact.solve_with_fault(stokes) + st["converged"] = bool(info.get("converged")) + else: + stokes.solve() + st["converged"] = stokes.snes.getConvergedReason() > 0 + return st + + +def stress_expression(st): + """Deviatoric stress components (sympy 2x2) for this rep's law.""" + E = st["stokes"].strainrate + E = E - sympy.eye(2) * (E[0, 0] + E[1, 1]) / 2 + eta0 = st["eta0"].sym[0] + if st["rep"] in ("ti", "hybrid"): + nd = st["ndir"].sym + nv = sympy.Matrix([nd[0], nd[1]]) + tv = sympy.Matrix([-nd[1], nd[0]]) + e_nt = (nv.T * E * tv)[0, 0] + return (2 * eta0 * E + - 2 * (eta0 - st["eta1"].sym[0]) * e_nt + * (nv * tv.T + tv * nv.T)) + return 2 * eta0 * E + + +def nodal_stress(st, name="tP1"): + """Project the stress COMPONENTS to continuous P1; returns + (P1 variable, txx, tyy, txy) — deviatoric, nodal.""" + mesh = st["mesh"] + S = uw.discretisation.MeshVariable(f"{name}{st['rep'][0]}", mesh, 1, + degree=1) + proj = uw.systems.Projection(mesh, S) + tau = stress_expression(st) + comps = {} + for key, expr in (("xx", tau[0, 0]), ("yy", tau[1, 1]), + ("xy", tau[0, 1])): + proj.uw_function = expr + proj.solve(zero_init_guess=True) + comps[key] = np.asarray(S.data[:, 0]).copy() + d = 0.5 * (comps["xx"] + comps["yy"]) + return S, comps["xx"] - d, comps["yy"] - d, comps["xy"] + + +def nodal_material(st, S): + """eta0, eta1, director at the P1 nodes from the NEAREST cell + centroid (the honoured, cell-based paint).""" + cen = np.asarray(st["eta0"].coords) + idx = cKDTree(cen).query(np.asarray(S.coords))[1] + eta0 = st["eta0_vals"][idx] + if st["rep"] in ("ti", "hybrid"): + eta1 = np.asarray(st["eta1"].array).reshape(len(cen))[idx] + nv = np.asarray(st["ndir"].array).reshape(len(cen), 2)[idx] + else: + eta1 = eta0.copy() + nv = np.tile(G.ET, (len(idx), 1)) + return eta0, eta1, nv + + +def nodal_strainrate(st, name="eP1"): + """THE reliable strain-rate map (measured 2026-08-25, TI band, coarse, + eta ratio 1e3): project the strain-rate COMPONENTS to continuous P1 + and form the invariant in numpy. In-band median 5.3 / p90 10.2 / + max 15.9 against the direct in-cell truth 5.6 / 9.7 / 13; far field + 0.98 = the drive. The one-ring smear at the band edge is the only + cost. Returns (P1 variable, edot_II).""" + S, exx, eyy, exy = nodal_strainrate_components(st, name=name) + return S, np.sqrt(0.5 * (exx ** 2 + eyy ** 2) + exy ** 2) + + +def nodal_strainrate_components(st, name="eP1"): + """The deviatoric strain-rate COMPONENTS projected to continuous P1; + returns (P1 variable, exx, eyy, exy) in the variable's DOF order — + the recipe of :func:`nodal_strainrate`, components kept.""" + mesh = st["mesh"] + S = uw.discretisation.MeshVariable(f"{name}{st['rep'][0]}", mesh, 1, + degree=1) + proj = uw.systems.Projection(mesh, S) + E = st["stokes"].strainrate + comps = {} + for key, expr in (("xx", E[0, 0]), ("yy", E[1, 1]), ("xy", E[0, 1])): + proj.uw_function = expr + proj.solve(zero_init_guess=True) + comps[key] = np.asarray(S.data[:, 0]).copy() + d = 0.5 * (comps["xx"] + comps["yy"]) + return S, comps["xx"] - d, comps["yy"] - d, comps["xy"] + + +def nodal_strainrate_from_stress(st, S, txx, tyy, txy): + """The material-jump recipe (recovered stress / own eta; TI in the + (n, t) frame, e_nt = tau_nt / (2 eta1)). CAVEAT (measured): only + reliable when the recovered interface stress is accurate to better + than eta0/eta1 — at coarse w with a 1e3 ratio it is NOT (in-band + p90 33, max 1037 vs truth 9.7 / 13): the band-edge recovery error is + amplified 1000x. Kept for the single-cell-layer case it was measured + on; use nodal_strainrate for a multi-cell band.""" + eta0, eta1, nv = nodal_material(st, S) + tv = np.column_stack([-nv[:, 1], nv[:, 0]]) + # rotate tau into (n, t) + tnn = (nv[:, 0] ** 2 * txx + 2 * nv[:, 0] * nv[:, 1] * txy + + nv[:, 1] ** 2 * tyy) + ttt = (tv[:, 0] ** 2 * txx + 2 * tv[:, 0] * tv[:, 1] * txy + + tv[:, 1] ** 2 * tyy) + tnt = (nv[:, 0] * tv[:, 0] * txx + + (nv[:, 0] * tv[:, 1] + nv[:, 1] * tv[:, 0]) * txy + + nv[:, 1] * tv[:, 1] * tyy) + enn, ett, ent = tnn / (2 * eta0), ttt / (2 * eta0), tnt / (2 * eta1) + return np.sqrt(0.5 * (enn ** 2 + ett ** 2) + ent ** 2) + + +def stress_invariant(txx, tyy, txy): + return np.sqrt(0.5 * (txx ** 2 + tyy ** 2) + txy ** 2) + + +def native_slip_profile(st, label, mask=True): + """The NATIVE slip gauge (Louis's ruling, 2026-08-25) on a live + state: split = tangential pair jump along the cut; TI = band + integral of the director-plane shear strain rate along the + mid-line (5 stations across the band). Samples inside another + strand's band are dropped when ``mask``. Returns (points, slip) + with points ON the trace (for a coloured-trace overlay).""" + from underworld3.utilities import fault_contact + from scipy.spatial import cKDTree + TR = dict(st["traces"]) + trace = np.asarray(TR[label]) + W = st["W"] + others = (np.vstack([q for l2, q in st["traces"] if l2 != label]) + if mask else None) + + def own(X): + if others is None: + return np.ones(len(X), dtype=bool) + k = cKDTree(others).query(X)[0] > W + return k if k.sum() >= 3 else np.ones(len(X), dtype=bool) + + if st["rep"] in ("split", "hybrid") and label not in st.get("uncut", ()): + coords, jumps, normals = fault_contact.fault_pair_jumps( + st["stokes"], label, st["stokes"]._rotated_freeslip_info) + jn = np.einsum("ij,ij->i", jumps, normals) + tang = np.linalg.norm(jumps - jn[:, None] * normals, axis=1) + keep = own(coords) + return coords[keep], tang[keep] + P = trace[::2] + t = np.gradient(P, axis=0) + t /= np.linalg.norm(t, axis=1)[:, None] + n = np.column_stack([-t[:, 1], t[:, 0]]) + E = st["stokes"].strainrate + # stations strictly INSIDE the band, off the mid-line and rail + # edges (points ON an edge locate into either neighbour, host + # cells included): the band-integral is w x the mean of 2 e_nt + offs = W * np.array([-0.45, -0.3, -0.1, 0.1, 0.3, 0.45]) + ent = np.zeros((len(offs), len(P))) + for k, o in enumerate(offs): + Q = P + o * n + exx = np.asarray(uw.function.evaluate(E[0, 0], Q)).ravel() + eyy = np.asarray(uw.function.evaluate(E[1, 1], Q)).ravel() + exy = np.asarray(uw.function.evaluate(E[0, 1], Q)).ravel() + ent[k] = (n[:, 0] * t[:, 0] * exx + + (n[:, 0] * t[:, 1] + n[:, 1] * t[:, 0]) * exy + + n[:, 1] * t[:, 1] * eyy) + # band slip = v_t difference across the band + one-cell skirt + # (continuous field -> robust; in-cell gradient integrals are + # vertex-phase sensitive at fine w) + skirt = W / 2 + st["s_f"] + v = st["v"] + vp = np.asarray(uw.function.evaluate(v.sym, P + skirt * n) + ).reshape(len(P), -1)[:, :2] + vm = np.asarray(uw.function.evaluate(v.sym, P - skirt * n) + ).reshape(len(P), -1)[:, :2] + slip = np.abs(np.einsum("ij,ij->i", vp - vm, t)) + keep = own(P) + return P[keep], slip[keep] + + +def midline_shear_rate(st, trace, eps=0.02): + """The TI band's slip measure of record (Louis, 2026-08-25, reaffirmed + 2026-09-05): the shear strain rate resolved on the director plane + ALONG THE MID-LINE, 2 e_nt = 2 n.E.t at the spine vertices. + + Computed as the NODAL RECOVERY of the P2 velocity's own gradient: in + every cell of the vertex's fan the strain rate is evaluated at the + vertex (a point stepped ``eps`` of the way toward the cell's + centroid, so the locator lands in THAT cell — the P2 gradient is + linear within a cell, so the step costs O(eps)) and the fan is + averaged with the cells' angles at the vertex as weights. No span, + so a neighbouring band cannot enter; no L2 projection, so no + overshoot at the peak (the P1 projection read 91 on the fine Main + where the in-cell truth is ~55); defined through a junction. The + fan is summed over ranks (a shared vertex's fan is split). + Returns (points, 2 e_nt) for the trace vertices that are mesh + vertices (every placed trace vertex is a spine vertex). + """ + from mpi4py import MPI + from scipy.spatial import cKDTree + from underworld3.utilities.reconnect import _coords + mesh = st["mesh"] + dm = mesh.dm + comm = uw.mpi.comm + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + X = _coords(dm)[: vE - vS] + P = np.asarray(trace, dtype=float) + d, idx = cKDTree(X).query(P) + on = d < 1e-6 + t = np.gradient(P, axis=0) + t /= np.linalg.norm(t, axis=1)[:, None] + n = np.column_stack([-t[:, 1], t[:, 0]]) + E = st["stokes"].strainrate + # one batch of evaluation points: every (vertex, fan cell) pair + pts, owner, weight = [], [], [] + for k in np.flatnonzero(on): + iv = int(idx[k]) + for c in dm.getTransitiveClosure(iv + vS, useCone=False)[0]: + c = int(c) + if not (cS <= c < cE): + continue + verts = [int(q) - vS for q in dm.getTransitiveClosure(c)[0] + if vS <= int(q) < vE] + cen = X[verts].mean(axis=0) + others = [q for q in verts if q != iv] + a, b = X[others[0]] - X[iv], X[others[1]] - X[iv] + ang = abs(np.arctan2(a[0] * b[1] - a[1] * b[0], a @ b)) + pts.append(X[iv] + eps * (cen - X[iv])) + owner.append(k) + weight.append(ang) + num = np.zeros(len(P)) + den = np.zeros(len(P)) + if pts: + Q = np.asarray(pts, dtype=float) + exx = np.asarray(uw.function.evaluate(E[0, 0], Q)).ravel() + eyy = np.asarray(uw.function.evaluate(E[1, 1], Q)).ravel() + exy = np.asarray(uw.function.evaluate(E[0, 1], Q)).ravel() + owner = np.asarray(owner) + weight = np.asarray(weight) + ent = (n[owner, 0] * t[owner, 0] * exx + + (n[owner, 0] * t[owner, 1] + n[owner, 1] * t[owner, 0]) * exy + + n[owner, 1] * t[owner, 1] * eyy) + np.add.at(num, owner, weight * ent) + np.add.at(den, owner, weight) + if comm.size > 1: + for arr in (num, den): + comm.Allreduce(MPI.IN_PLACE, arr, op=MPI.SUM) + rate = np.full(len(P), np.nan) + have = den > 0 + rate[have] = 2.0 * np.abs(num[have] / den[have]) + return P, rate + + +def coloured_trace(points, values, z=0.003): + """A pv polyline through ``points`` carrying ``values`` as point + scalars — the coloured-trace overlay (sorted along the trace by + the order given).""" + import pyvista as pv + line = pv.lines_from_points( + np.column_stack([points, np.full(len(points), z)])) + line.point_data["slip"] = np.asarray(values) + return line + + +def yielded_cells(st): + """Band cells on the plastic branch (2 eta0 edot_II > tau_y) from the + in-cell strain rate at the band centroids. Returns (band_mask, + yielded_within_band) — both over ALL cells for cell-data rendering.""" + band = st["pinfo"]["band"] + cen = np.asarray(st["eta0"].coords)[band] + E = st["stokes"].strainrate + c = [np.asarray(uw.function.evaluate(E[i, j], cen)).ravel() + for i, j in ((0, 0), (1, 1), (0, 1))] + d = 0.5 * (c[0] + c[1]) + e2 = np.sqrt(0.5 * ((c[0] - d) ** 2 + (c[1] - d) ** 2) + c[2] ** 2) + yielded = np.zeros(len(band), dtype=bool) + yielded[np.flatnonzero(band)] = (2.0 * st["eta0_vals"][band] * e2 + > st["tau_y"]) + return band, yielded diff --git a/articles/putting-a-fault-in-a-mesh/examples/sf_note_stress_slip.py b/articles/putting-a-fault-in-a-mesh/examples/sf_note_stress_slip.py new file mode 100644 index 0000000..8651fb9 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/sf_note_stress_slip.py @@ -0,0 +1,142 @@ +"""Figure for UWTN 2026-017: the cut and the band on the same network, +stress and slip, at the width where they agree. + +Both columns are the rig at w = 0.005 (superfine): the split node, and +the TI band at eta_1/w = 0.1, the ratio the convergence test holds +fixed. Background is log10 tau_II on one scale; every trace is drawn as +a tube coloured by the SIGNED slip rate in each representation's own +currency -- the pair jump on the cut, the velocity jump across the +band's own edges -- on one two-slope map with grey at zero. The lower +row zooms the Y junction, which is the one place the two disagree. + +Dumps come from sf_split_seam_viz.py: + + ../../uwrun sf_split_seam_viz.py -uw_rep split -uw_res superfine + ../../uwrun sf_split_seam_viz.py -uw_rep ti -uw_res superfine -uw_eta1 5e-4 + + ../../uwrun sf_note_stress_slip.py -> sf_note_stress_slip.png +""" +import numpy as np +import pyvista as pv +from matplotlib.colors import (LinearSegmentedColormap, Normalize, + TwoSlopeNorm) + +import rig_geometry as G + +pv.OFF_SCREEN = True +W, S_F = G.RES["superfine"] +ETA1 = 5e-4 + +# the upper stepover closed to the split's floor (one uncut edge): +# ../../uwrun sf_split_seam_viz.py -uw_rep split -uw_res superfine -uw_cont_gap_rungs 1 +# ../../uwrun sf_split_seam_viz.py -uw_rep ti -uw_res superfine -uw_eta1 5e-4 -uw_cont_gap_rungs 1 +CASES = [("split node", "fields_split_superfine_gather_np1_rank0_cont1.npz"), + (f"TI band, eta_1/w = {ETA1 / W:g}", + "fields_ti_superfine_gather_np1_rank0_eta0.0005_cont1.npz")] + +G.CONT_GAP_RUNGS = 1.0 +TR = dict(G.traces(S_F, W)) +# the Y: the Splay's stop-short tip beside the Main's outgoing limb +CY = TR["Splay"][0] +# the stepover: the Main's upper tip, the Cont segment one edge beyond it +CS = 0.5 * (TR["Main"][-1] + TR["Cont"][0]) +# each camera fills its 1000x1000 viewport exactly, so all four +# panels carry the same pixels; the gutter is added afterwards +PANEL, GUTTER = 1000, 22 +VIEWS = [((0.5, 0.5), 0.5, ""), + ((CY[0], CY[1]), 0.085, "the Y junction"), + ((CS[0], CS[1]), 0.06, "the stepover")] + +dumps = [np.load(f, allow_pickle=True) for _t, f in CASES] +allv = np.concatenate([d["logtau"] for d in dumps]) +clim = (float(np.percentile(allv, 2)), float(np.percentile(allv, 99.5))) + +# the trace colour is the SIGNED slip in each representation's own +# currency, on one map with grey at zero: blue sinistral, green dextral, +# each sense scaled to its own extreme so the small one stays visible. +# With no sinistral slip anywhere (this network, at this width) the blue +# half would be a range that does not exist, so the map drops to its +# dextral half and the bar runs 0 to the peak. +SENSE = ["#1f4fbf", "#7aa6ff", "#d9d9d9", "#7fe07a", "#118a2e"] +alls = np.concatenate([d["slips"][:, 2] for d in dumps]) +alls = alls[np.isfinite(alls)] +lo, hi = float(alls.min()), float(alls.max()) +TWO_SIDED = lo < -1e-3 +if TWO_SIDED: + SLIP_CMAP = LinearSegmentedColormap.from_list("slip_sense", SENSE) + NORM = TwoSlopeNorm(vcenter=0.0, vmin=lo, vmax=hi) + TICKS = {0.0: f"{lo:.3f}", 0.25: f"{lo / 2:.3f}", 0.5: "0", + 0.75: f"{hi / 2:.2f}", 1.0: f"{hi:.2f}"} +else: + SLIP_CMAP = LinearSegmentedColormap.from_list("slip_dextral", SENSE[2:]) + NORM = Normalize(vmin=0.0, vmax=hi) + TICKS = {0.0: "0", 0.5: f"{hi / 2:.2f}", 1.0: f"{hi:.2f}"} +KEY = "slip rate (+ dextral)" + +pl = pv.Plotter(off_screen=True, shape=(len(VIEWS), len(CASES)), + window_size=(PANEL * len(CASES), PANEL * len(VIEWS)), + border=False) +for j, ((title, _f), d) in enumerate(zip(CASES, dumps)): + for i, ((cx, cy), scale, ztag) in enumerate(VIEWS): + pl.subplot(i, j) + pl.set_background("white") + pvm = pv.UnstructuredGrid(d["cells"], d["celltypes"], d["points"]) + pvm.point_data["logtau"] = d["logtau"] + pl.add_mesh(pvm, scalars="logtau", cmap="magma", clim=clim, + show_edges=False, lighting=False, + show_scalar_bar=(i == 0 and j == 0), + scalar_bar_args=dict(title="log10 tau_II", n_labels=4, + color="black", vertical=True, + position_x=0.84, position_y=0.06)) + if i > 0: + pl.add_mesh(pvm.extract_all_edges(), color="white", + line_width=0.4, opacity=0.3, lighting=False) + arr = d["slips"] + for seg in np.split(arr, np.flatnonzero(np.isnan(arr[:, 0]))): + seg = seg[~np.isnan(seg[:, 0])] + if len(seg) < 2: + continue + t = seg[:, :2] - seg[:, :2].mean(axis=0) + _u, _s, vt = np.linalg.svd(t, full_matrices=False) + seg = seg[np.argsort(t @ vt[0])] + vals = np.where(np.isfinite(seg[:, 2]), seg[:, 2], 0.0) + if seg.shape[1] > 3: + vals = np.where(seg[:, 3] > 0, vals, 0.0) + line = pv.lines_from_points(np.column_stack( + [seg[:, :2], np.full(len(seg), 0.004)])) + line.point_data[KEY] = np.asarray(NORM(vals), dtype=float) + pl.add_mesh(line, scalars=KEY, cmap=SLIP_CMAP, clim=(0.0, 1.0), + line_width=6 if i == 0 else 10, lighting=False, + render_lines_as_tubes=True, + show_scalar_bar=(i == 0 and j == 1), + annotations=TICKS, + scalar_bar_args=dict(title=KEY, n_labels=0, + color="black", vertical=True, + position_x=0.84, + position_y=0.45)) + txt = f"{title}, w = {W}" if i == 0 else f"{ztag} — {title}" + pl.add_text(txt, font_size=13, color="black", + position=(18.0, 960.0)) + pl.view_xy() + pl.camera.parallel_projection = True + pl.camera.focal_point = (cx, cy, 0.0) + pl.camera.parallel_scale = scale +img = pl.screenshot(None, return_img=True) +pl.close() + +# equal panel size and equal margins: cut the render into its four +# viewports and lay them out on white with one gutter width everywhere +nr, nc = len(VIEWS), len(CASES) +H = nr * PANEL + (nr + 1) * GUTTER +Wd = nc * PANEL + (nc + 1) * GUTTER +canvas = np.full((H, Wd, img.shape[2]), 255, dtype=img.dtype) +for i in range(nr): + for j in range(nc): + tile = img[i * PANEL:(i + 1) * PANEL, j * PANEL:(j + 1) * PANEL] + y0 = GUTTER + i * (PANEL + GUTTER) + x0 = GUTTER + j * (PANEL + GUTTER) + canvas[y0:y0 + PANEL, x0:x0 + PANEL] = tile +out = "sf_note_stress_slip.png" +from PIL import Image +Image.fromarray(canvas).save(out) +print(f"wrote {out} slip range [{lo:.4f}, {hi:.4f}] two-sided {TWO_SIDED} logtau clim {clim}") diff --git a/articles/putting-a-fault-in-a-mesh/examples/sf_split_seam_viz.py b/articles/putting-a-fault-in-a-mesh/examples/sf_split_seam_viz.py new file mode 100644 index 0000000..328ced2 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/sf_split_seam_viz.py @@ -0,0 +1,466 @@ +"""The SPLIT rig solved through the partition seam, against serial, the +fine grid and the TI model: stress and slip. + +Two stages, as the other parallel figures. Under mpirun (or serial) the +rig is built in the realisation and seam mode asked for, solved (the +frictionless contact solve for the split, the weak plane for TI), the +stress COMPONENTS are recovered to continuous P1 (sf_fields' recipe, +invariant formed nodally) and every rank dumps its P1 mesh, the stress, +its seam vertices and the NATIVE slip along each strand: the tangential +pair jump on the cut (split), the tangential velocity jump across the +band (TI). In serial the dumps are rendered side by side. + + ../../uwrun sf_split_seam_viz.py -uw_rep split -uw_res coarse # serial + ../../uwrun sf_split_seam_viz.py -uw_rep split -uw_res fine + ../../uwrun -np 2 sf_split_seam_viz.py -uw_rep split -uw_res fine -uw_seams conform + ../../uwrun sf_split_seam_viz.py -uw_rep ti -uw_res fine + ../../uwrun sf_split_seam_viz.py -uw_render 1 + -> sf_split_seam_logtau_slip.png (fields, trace coloured by slip) + sf_split_seam_slip_profiles.png (slip along the Main, four cases) +""" +import glob + +import numpy as np +import underworld3 as uw + +params = uw.Params( + uw_rep=uw.Param("split", "split | ti"), + uw_res=uw.Param("coarse", "coarse | fine | superfine"), + uw_seams=uw.Param("gather", "gather | ligament | conform"), + uw_render=uw.Param(0.0, "1 = render the dumps (serial)"), + uw_np=uw.Param(2.0, "render: the np of the parallel split dump"), + uw_carve_clearance=uw.Param(1.0, "build() carve_clearance (fine: 1.0)"), + uw_eta1=uw.Param(1e-3, "TI weak shear viscosity in the footprints"), + uw_tol=uw.Param(1e-5, "stokes tolerance"), + uw_splay_gap_rungs=uw.Param(-1.0, "the Splay's stop-short gap from the " + "Main in rungs (rig default when < 0); 0 = " + "a genuine FORK, the band's own junction " + "(the cut cannot take touching strands)"), + uw_cont_gap_rungs=uw.Param(-1.0, "the upper stepover's tip-to-tip gap " + "in rungs (rig default when < 0); 1 = the " + "split's floor, one uncut edge; 0 = the " + "CONTINUOUS control, one cut the whole " + "line, no Cont strand"), + uw_strands=uw.Param("all", "all | main | main+splay (the Main alone: no " + "fork, so the band and the cut represent the SAME " + "structure; with the Splay: the junction and its " + "partition; sf_main_only_convergence.py and " + "sf_main_splay_junction.py draw those dumps)"), +) +import rig_geometry as G +import sf_fields as F + +REP, RES, SEAMS = str(params.uw_rep), str(params.uw_res), str(params.uw_seams) +STRANDS = str(params.uw_strands) +GAP = float(params.uw_splay_gap_rungs) +if GAP >= 0.0: + G.SPLAY_GAP_RUNGS = GAP +CONT = float(params.uw_cont_gap_rungs) +if CONT >= 0.0: + G.CONT_GAP_RUNGS = CONT +H_BG = 1 / 8 + + +def dump_name(rep, res, seams, np_, rank, eta1=None): + """TI dumps carry the weak viscosity when it is not the default: the + band's interface strength is eta_1 / w, so a width sweep at fixed + eta_1 doubles the strength each halving — a convergence test holds + the ratio.""" + tag = "" if (eta1 is None or rep != "ti" or abs(eta1 - 1e-3) < 1e-12) \ + else f"_eta{eta1:g}" + if STRANDS == "main": + tag += "_main" + elif STRANDS == "main+splay": + tag += "_mainsplay" + if GAP >= 0.0: + tag += f"_gap{GAP:g}" + if CONT >= 0.0: + tag += f"_cont{CONT:g}" + return f"fields_{rep}_{res}_{seams}_np{np_}_rank{rank}{tag}.npz" + + +if float(params.uw_render) < 0.5: + from mpi4py import MPI + import underworld3.visualisation as vis + from underworld3.utilities import fault_contact + from underworld3.utilities.place_surface import _shared_point_flags + comm = MPI.COMM_WORLD + W, S_F = G.RES[RES] + traces = G.traces(S_F, W) + keep = {"main": ("Main",), "main+splay": ("Main", "Splay")}.get(STRANDS) + if keep is not None: + traces = [tr for tr in traces if tr[0] in keep] + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=H_BG, + regular=False, qdegree=2, refinement=1) + net = uw.meshing.FaultNetwork( + [(lbl, np.asarray(P, dtype=float)) for lbl, P in traces], + hierarchy=[lbl for lbl, _ in traces]) + net.prepare(h=S_F, ligament=1.0, verbose=False) + mesh = net.build(base=base, h_far=H_BG / 2, width=W, realisation=REP, + max_levels=1, margin_rings=2, seams=SEAMS, + carve_clearance=float(params.uw_carve_clearance)) + x, y = mesh.X + tag = f"{REP[0]}{RES[0]}" + v = uw.discretisation.MeshVariable(f"v{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) + eta0 = uw.discretisation.MeshVariable(f"eta0{tag}", mesh, 1, degree=0) + eta0_vals = np.ones(len(np.asarray(eta0.coords))) + eta0.array[:, 0, 0] = eta0_vals + st = dict(rep=REP, mesh=mesh, stokes=stokes, v=v, eta0=eta0, + eta0_vals=eta0_vals, eta1=None, ndir=None, foot=None, + traces=traces, W=W, s_f=S_F) + if REP == "ti": + eta1, ndir, foot = G.ti_fields(mesh, net.info["footprints"], + float(params.uw_eta1), eta0_vals, tag, + w=W, traces_=traces) + stokes.constitutive_model = \ + uw.constitutive_models.TransverseIsotropicFlowModel + P = stokes.constitutive_model.Parameters + P.shear_viscosity_0 = eta0.sym[0] + P.shear_viscosity_1 = eta1.sym[0] + P.director = ndir.sym + st.update(eta1=eta1, ndir=ndir, foot=foot) + else: + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta0.sym[0] + stokes.bodyforce = [0.0, 0.0] + ES, ET = G.ES, G.ET + t_sym = (x - 0.5) * float(ET[0]) + (y - 0.5) * float(ET[1]) + for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc((2.0 * t_sym * float(ES[0]), + 2.0 * t_sym * float(ES[1])), wall) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = float(params.uw_tol) + if REP == "split": + for lbl, _tr in traces: + stokes.add_fault_bc(0, boundary=lbl) + info = fault_contact.solve_with_fault(stokes) + converged = bool(info.get("converged")) + else: + stokes.solve() + info = None + converged = stokes.snes.getConvergedReason() > 0 + S_t, txx, tyy, txy = F.nodal_stress(st, name="tP1") + tau = F.stress_invariant(txx, tyy, txy) + + # the NATIVE slip along each strand, made global. Split: every + # coincident pair's tangential jump, all-gathered (rank-owned pairs + # only, so the seam pair is counted once). TI: the tangential velocity + # jump across band + one-cell skirt, counted where this rank owns both + # probes (evaluate answers for unowned points by extrapolation). + slip_profiles = {} + if REP == "split": + # SIGNED: the tangential jump v_plus - v_minus projected on the + # trace's own direction at the nearest trace vertex; Plus is the + # left of that direction (add_fault orients the chain by the + # polyline), the same convention as the band's left-minus-right. + # POSITIVE IS RIGHT-LATERAL (dextral), negative left-lateral, and + # the sign does not depend on the way the trace is walked + # (reversing it swaps left and right and flips t together). + from scipy.spatial import cKDTree + for lbl, tr in traces: + coords, jumps, normals = fault_contact.fault_pair_jumps( + stokes, lbl, stokes._rotated_freeslip_info, gather=True) + if len(coords) == 0: + continue + P = np.asarray(tr, dtype=float) + t = np.gradient(P, axis=0) + t /= np.linalg.norm(t, axis=1)[:, None] + near = cKDTree(P).query(coords)[1] + tang = np.einsum("ij,ij->i", jumps, t[near]) + slip_profiles[lbl] = np.column_stack( + [coords, tang, np.ones(len(tang))]) + else: + from scipy.spatial import cKDTree + # the TI record: the shear strain rate resolved on the director + # plane along the mid-line, at the spine vertices (P1-recovered + # components) — a strain rate, no span, continuous through a + # junction. The velocity-jump gauge across band + skirt is kept + # for the peak numbers only (a span one cell wide wanders in + # location and, past a junction, into the neighbouring band). + midline = {} + for lbl, tr in traces: + Pm, rate = F.midline_shear_rate(st, tr) # collective + midline[lbl] = np.column_stack([Pm, rate]) + + def global_slip(P): + t = np.gradient(P, axis=0) + t /= np.linalg.norm(t, axis=1)[:, None] + n = np.column_stack([-t[:, 1], t[:, 0]]) + # the band's OWN edges, no skirt: the layer's shear is a triangle + # peaked at the mid-line that reaches the rails (the transect + # probe, 2026-09-05: the jump across +-W/2 is within 3-6% of the + # jump across +-(W/2 + s_f)), and a probe on the rail cannot + # wander into a neighbouring band the mask keeps W away + skirt = W / 2 + Qp, Qm = P + skirt * n, P - skirt * n + vp = np.asarray(uw.function.evaluate(v.sym, Qp)).reshape(len(P), -1)[:, :2] + vm = np.asarray(uw.function.evaluate(v.sym, Qm)).reshape(len(P), -1)[:, :2] + slip = np.einsum("ij,ij->i", vp - vm, t) # SIGNED: left - right + own = ((mesh._robust_owning_cells(np.ascontiguousarray(Qp)) >= 0) + & (mesh._robust_owning_cells(np.ascontiguousarray(Qm)) >= 0)) + # gathered by the owning rank: the value where owned, and a + # count, so the sign survives (a MAX would not keep it) + val = np.where(own, slip, 0.0) + cnt = own.astype(float) + comm.Allreduce(MPI.IN_PLACE, val, op=MPI.SUM) + comm.Allreduce(MPI.IN_PLACE, cnt, op=MPI.SUM) + out = np.where(cnt > 0, val / np.maximum(cnt, 1.0), np.nan) + return out + + from scipy.spatial import cKDTree + all_traces = [np.asarray(tr, dtype=float) for _l, tr in traces] + # the Main and its Cont segment are ONE line: the gauge must not + # mask the Main's samples as "inside the Cont's band" at a closed + # stepover, or the tube shows a nick where the band has none + ONE_LINE = {"Main", "Cont"} + for k, (lbl, tr) in enumerate(traces): + P = np.asarray(tr, dtype=float)[::2] + others = [q for j, (l2, q) in enumerate(zip([l for l, _ in traces], + all_traces)) + if j != k and not (lbl in ONE_LINE and l2 in ONE_LINE)] + clear = (cKDTree(np.vstack(others)).query(P)[0] > W if others + else np.ones(len(P), dtype=bool)) + # the raw jump everywhere, and whether the trace point is + # clear of every other strand's band (column 3): the render + # masks with it, the junction studies show it + slip = global_slip(P) + slip_profiles[lbl] = np.column_stack([P, slip, clear.astype(float)]) + + # the P1 stress on the mesh's OWN cells (a split mesh must never be + # re-triangulated through its DOF cloud: the slit would close) + pvm = vis.mesh_to_pv_mesh(mesh).cast_to_unstructured_grid() + dm = mesh.dm + vS, vE = dm.getDepthStratum(0) + pStart, _ = dm.getChart() + shared = _shared_point_flags(dm).astype(bool)[vS - pStart: vE - pStart] + from underworld3.utilities.reconnect import _coords + Xv = _coords(dm)[: vE - vS] + # S_t is P1: one DOF per vertex in the DM's vertex order, which is + # also the pv mesh's point order — the split's coincident vertices + # keep their own DOFs, one per side + tau = np.asarray(tau).ravel() + if len(tau) != vE - vS: + raise RuntimeError("the P1 stress does not have one DOF per vertex") + logtau = np.log10(tau + 1e-12) + slips = np.vstack([np.vstack([sp, [np.nan] * 4]) + for sp in slip_profiles.values()]) \ + if slip_profiles else np.zeros((0, 4)) + labels = list(slip_profiles) + if REP == "ti": + midline_arr = np.vstack([np.vstack([mp, [np.nan] * 3]) + for mp in midline.values()]) + midline_labels = list(midline) + else: + midline_arr, midline_labels = np.zeros((0, 3)), [] + np.savez(dump_name(REP, RES, SEAMS, comm.size, comm.rank, + eta1=float(params.uw_eta1)), + points=np.asarray(pvm.points), + cells=np.asarray(pvm.cells, dtype=np.int64), + celltypes=np.asarray(pvm.celltypes, dtype=np.uint8), + logtau=logtau, + shared_xy=Xv[shared], + traces=np.vstack([np.vstack([tr, [np.nan, np.nan]]) + for _l, tr in traces]), + slips=slips, slip_labels=np.asarray(labels), + midline=midline_arr, midline_labels=np.asarray(midline_labels)) + peaks = {lbl: float(np.nanmax(np.where(sp[:, 3] > 0, np.abs(sp[:, 2]), np.nan))) + for lbl, sp in slip_profiles.items()} + if REP == "ti": + uw.pprint("mid-line 2 e_nt peaks: " + ", ".join( + f"{lbl}: {float(np.nanmax(mp[:, 2])):.2f}" for lbl, mp in midline.items())) + uw.pprint(f"[{REP} {RES} {SEAMS} np{comm.size}] converged {converged}; " + f"cells {comm.gather(int(dm.getHeightStratum(0)[1]))}; " + f"peak slip {{ {', '.join(f'{k}: {v_:.4f}' for k, v_ in peaks.items())} }}; dumped") +else: + import pyvista as pv + pv.OFF_SCREEN = True + NP = int(float(params.uw_np)) + import re + + def ranks_of(pattern): + # the rank wildcard must not pick up the Main-only / Main+Splay / + # other-eta dumps, whose names carry a suffix after the rank and + # whose MESHES differ (three meshes drawn over each other dappled + # the zoomed field) + return sorted(f for f in glob.glob(pattern) + if re.search(r"rank\d+\.npz$", f)) + + cases = [("split, coarse, serial", + ranks_of(dump_name("split", "coarse", "gather", 1, "*"))), + ("split, fine, serial", + ranks_of(dump_name("split", "fine", "gather", 1, "*"))), + (f"split, fine, np={NP} through the seam", + ranks_of(dump_name("split", "fine", "conform", NP, "*"))), + ("TI weak plane, fine, serial", + ranks_of(dump_name("ti", "fine", "gather", 1, "*")))] + for title, files in cases: + if not files: + raise SystemExit(f"no dump for {title}: run the dump stage first") + seam = np.vstack([np.load(f, allow_pickle=True)["shared_xy"] + for f in cases[2][1]]) + seam = seam[np.linalg.norm(seam - 0.5, axis=1) < 0.2] + czoom = seam.mean(axis=0) if len(seam) else np.array([0.5, 0.5]) + VIEWS = [((0.5, 0.5), 0.54, ""), + ((czoom[0], czoom[1]), 0.12, " — the Main's seam crossing")] + allv = np.concatenate([np.load(f, allow_pickle=True)["logtau"] + for _t, fs in cases for f in fs]) + clim = (float(np.percentile(allv, 2)), float(np.percentile(allv, 99.5))) + # the trace colour is the SIGNED slip (positive dextral) in the one + # currency for both representations — the pair jump on the cut, the + # jump across the band's own edges — on a diverging map, neutral grey + # at zero, blue sinistral and green dextral (both stand off the magma + # stress), symmetric about zero + # two-slope: zero stays at the map's centre, each sense scaled to its + # own extreme (a symmetric range left the sinistral lobes, a twentieth + # of the Main's dextral slip, invisible) + from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm + SLIP_CMAP = LinearSegmentedColormap.from_list( + "slip_sense", ["#1f4fbf", "#7aa6ff", "#d9d9d9", "#7fe07a", "#118a2e"]) + allslip = np.concatenate([np.load(fs[0], allow_pickle=True)["slips"][:, 2] + for _t, fs in cases]) + allslip = allslip[np.isfinite(allslip)] + slip_lo = min(float(allslip.min()), -1e-6) + slip_hi = max(float(allslip.max()), 1e-6) + SLIP_NORM = TwoSlopeNorm(vcenter=0.0, vmin=slip_lo, vmax=slip_hi) + slip_ticks = {0.0: f"{slip_lo:.3f}", 0.25: f"{slip_lo / 2:.3f}", 0.5: "0", + 0.75: f"{slip_hi / 2:.2f}", 1.0: f"{slip_hi:.2f}"} + pl = pv.Plotter(off_screen=True, shape=(len(VIEWS), len(cases)), + window_size=(1000 * len(cases), 1000 * len(VIEWS)), + border=False) + for j, (title, files) in enumerate(cases): + for i, ((cx, cy), scale, ztag) in enumerate(VIEWS): + pl.subplot(i, j) + pl.set_background("white") + for f in files: + d = np.load(f, allow_pickle=True) + pvm = pv.UnstructuredGrid(d["cells"], d["celltypes"], d["points"]) + pvm.point_data["logtau"] = d["logtau"] + pl.add_mesh(pvm, scalars="logtau", cmap="magma", clim=clim, + show_edges=False, lighting=False, + show_scalar_bar=(i == 0 and j == 0), + scalar_bar_args=dict(title="log10 tau_II", n_labels=4, + color="black", vertical=True, + position_x=0.86, position_y=0.05)) + if i > 0: + pl.add_mesh(pvm.extract_all_edges(), color="white", + line_width=0.4, opacity=0.35, lighting=False) + if len(d["shared_xy"]): + pl.add_points(np.column_stack( + [d["shared_xy"], np.full(len(d["shared_xy"]), 0.003)]), + color="black", point_size=4 if i == 0 else 8, + render_points_as_spheres=True) + # the trace colour is each representation's OWN currency: the + # pair jump on the cut (split), the mid-line resolved shear + # strain rate 2 e_nt (TI) — two scalar bars. Rank 0's dump + # carries the profiles. Split pairs are unordered along the + # cut: order each strand's samples by their projection onto + # the trace's mean tangent. + d0 = np.load(files[0], allow_pickle=True) + sp_all = d0["slips"] + for seg in np.split(sp_all, np.flatnonzero(np.isnan(sp_all[:, 0]))): + seg = seg[~np.isnan(seg[:, 0])] + if len(seg) < 2: + continue + t = seg[:, :2] - seg[:, :2].mean(axis=0) + _u, _s, vt = np.linalg.svd(t, full_matrices=False) + order = np.argsort(t @ vt[0]) + seg = seg[order] + vals = np.where(np.isfinite(seg[:, 2]), seg[:, 2], 0.0) + if seg.shape[1] > 3: + vals = np.where(seg[:, 3] > 0, vals, 0.0) + line = pv.lines_from_points(np.column_stack( + [seg[:, :2], np.full(len(seg), 0.004)])) + key = "slip rate (signed; + dextral)" + line.point_data[key] = np.asarray(SLIP_NORM(vals), dtype=float) + pl.add_mesh(line, scalars=key, cmap=SLIP_CMAP, clim=(0.0, 1.0), + line_width=6, lighting=False, + render_lines_as_tubes=True, + show_scalar_bar=(i == 0 and j == 1), + annotations=slip_ticks, + scalar_bar_args=dict(title=key, n_labels=0, + color="black", vertical=True, + position_x=0.86, + position_y=0.45)) + pl.add_text(f"{title}{ztag}", font_size=11, color="black", + position="upper_left") + pl.view_xy() + pl.camera.parallel_projection = True + pl.camera.focal_point = (cx, cy, 0.0) + pl.camera.parallel_scale = scale + out = "sf_split_seam_logtau_slip.png" + pl.screenshot(out) + pl.close() + print(f"wrote {out}") + + # CONVERGENCE of the band to the cut, along the Main, in ONE + # currency — the velocity jump across the band's own edges for TI, + # the pair jump on the cut for the split (the integral of the + # resolved strain rate across the layer, either way): every split + # and TI dump on disk, by width. + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + def along(seg): + t = seg[:, :2] - np.asarray(G.C) + s_ = t @ G.ES + o = np.argsort(s_) + return s_[o], np.abs(seg[o, 2]) + + def strand(d0, name): + labels = list(d0["slip_labels"]) + arr = d0["slips"] + segs = [q[~np.isnan(q[:, 0])] for q in + np.split(arr, np.flatnonzero(np.isnan(arr[:, 0])))] + segs = [q for q in segs if len(q)] + return segs[labels.index(name)] + + fig, ax = plt.subplots(figsize=(9, 4.8)) + widths = {"coarse": G.RES["coarse"][0], "fine": G.RES["fine"][0], + "superfine": G.RES["superfine"][0]} + conv = [("split", "coarse", None, "C0-", 1.2), + ("split", "fine", None, "C0-", 1.8), + ("split", "superfine", None, "k-", 2.4), + ("ti", "fine", 1e-3, "C2--", 1.6), + ("ti", "superfine", 1e-3, "C3--", 1.6), + ("ti", "fine", 5e-4, "C2-.", 1.6), + ("ti", "superfine", 5e-4, "C3-.", 2.2)] + peaks = [] + for rep, res, eta1, sty, lw in conv: + files = sorted(f for f in glob.glob(dump_name(rep, res, "gather", 1, "*", + eta1=eta1)) + if re.search(r"rank\d+\.npz$", f)) + if not files: + continue + d0 = np.load(files[0], allow_pickle=True) + s_, v_ = along(strand(d0, "Main")) + if rep == "split": + name = f"split, w = {widths[res]}" + else: + name = (f"TI band, w = {widths[res]}, eta_1 = {eta1:g} " + f"(eta_1/w = {eta1 / widths[res]:g})") + ax.plot(s_, v_, sty, lw=lw, label=name) + peaks.append((name, float(np.nanmax(v_)))) + files = sorted(glob.glob(dump_name("split", "fine", "conform", NP, "*"))) + if files: + d0 = np.load(files[0], allow_pickle=True) + s_, v_ = along(strand(d0, "Main")) + ax.plot(s_, v_, "C0:", lw=1.4, + label=f"split, w = {widths['fine']}, np={NP} through the seam") + ax.set_xlabel("along the Main (from its centre)") + ax.set_ylabel("slip rate: pair jump (split) / v_t jump across the band (TI)") + ax.set_ylim(0, None) + ax.set_title("The TI band against the split node along the Main: " + "the width halved at fixed eta_1 and at fixed eta_1/w", + fontsize=10) + ax.legend(frameon=False, fontsize=7, loc="lower center") + fig.tight_layout() + for name, pk in peaks: + print(f"peak on the Main {name}: {pk:.4f}") + out2 = "sf_split_seam_slip_profiles.png" + fig.savefig(out2, dpi=160) + print(f"wrote {out2}") diff --git a/articles/putting-a-fault-in-a-mesh/examples/sf_stepover_closure.py b/articles/putting-a-fault-in-a-mesh/examples/sf_stepover_closure.py new file mode 100644 index 0000000..4a39ddf --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/sf_stepover_closure.py @@ -0,0 +1,67 @@ +"""The stepover closed to the split's floor, against the continuous +control, at superfine: how much of the through-going slip each +representation carries onto the segment beyond the gap. + +Four dumps from sf_split_seam_viz.py (-uw_cont_gap_rungs 1 and 0, for +-uw_rep split and ti -uw_eta1 5e-4). The Cont region is the arc beyond +the Main's design tip, S_MAIN[1]; in the continuous control it is the +same arc of the one uncut line. + + ../../uwrun sf_stepover_closure.py +""" +import numpy as np + +import rig_geometry as G + +RES = "superfine" +W, S_F = G.RES[RES] + + +def strands(f): + d = np.load(f, allow_pickle=True) + labels = list(d["slip_labels"]) + arr = d["slips"] + segs = [q[~np.isnan(q[:, 0])] for q in + np.split(arr, np.flatnonzero(np.isnan(arr[:, 0])))] + segs = [q for q in segs if len(q)] + return dict(zip(labels, segs)) + + +def line_profile(st): + """The Main-line family on one arc coordinate: s along ES from C.""" + parts = [st["Main"]] + ([st["Cont"]] if "Cont" in st else []) + P = np.vstack([q[:, :2] for q in parts]) + v = np.concatenate([q[:, 2] for q in parts]) + ok = np.concatenate([q[:, 3] > 0 if q.shape[1] > 3 + else np.ones(len(q), bool) for q in parts]) + s_ = (P - np.asarray(G.C)) @ G.ES + o = np.argsort(s_) + return s_[o], np.where(ok[o], v[o], np.nan) + + +def region_mean(s_, v, lo, hi): + m = (s_ > lo) & (s_ < hi) & np.isfinite(v) + return float(np.mean(np.abs(v[m]))) if m.any() else float("nan") + + +s_tip = G.S_MAIN[1] +s_end = s_tip + G.CONT_LEN +rows = [] +for rep, eta in (("split", ""), ("ti", "_eta0.0005")): + out = {} + for gap in (1, 0): + f = f"fields_{rep}_{RES}_gather_np1_rank0{eta}_cont{gap}.npz" + s_, v = line_profile(strands(f)) + out[gap] = (region_mean(s_, v, s_tip + 0.01, s_end - 0.01), + region_mean(s_, v, -0.2, 0.2), + float(np.nanmax(np.abs(v)))) + seg_closed, main_closed, pk_closed = out[1] + seg_cont, main_cont, pk_cont = out[0] + rows.append((rep, seg_closed, seg_cont, seg_closed / seg_cont, + main_closed, main_cont, pk_closed, pk_cont)) + +print(f"{'rep':<6}{'seg closed':>12}{'seg cont':>10}{'ratio':>8}" + f"{'main closed':>13}{'main cont':>11}{'peak closed':>13}{'peak cont':>11}") +for r in rows: + print(f"{r[0]:<6}{r[1]:>12.4f}{r[2]:>10.4f}{r[3]:>8.3f}" + f"{r[4]:>13.4f}{r[5]:>11.4f}{r[6]:>13.4f}{r[7]:>11.4f}") diff --git a/articles/putting-a-fault-in-a-mesh/examples/split-anatomy-data.json b/articles/putting-a-fault-in-a-mesh/examples/split-anatomy-data.json new file mode 100644 index 0000000..4de16bf --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/examples/split-anatomy-data.json @@ -0,0 +1 @@ +{"coords": [[0.0, 0.0], [0.25, 0.0], [0.25, 0.25], [0.0, 0.25], [0.25, 0.5], [0.0, 0.5], [0.25, 0.75], [0.0, 0.75], [0.25, 1.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.25], [0.5, 0.5], [0.5, 0.75], [0.5, 1.0], [0.75, 0.0], [0.75, 0.27296100594190537], [0.75, 0.5459220118838107], [0.75, 0.7729610059419054], [0.75, 1.0], [1.0, 0.0], [1.0, 0.29242640687119287], [1.0, 0.5848528137423857], [1.0, 0.7924264068711928], [1.0, 1.0], [1.25, 0.0], [1.25, 0.3054327719506772], [1.25, 0.6108655439013544], [1.25, 0.8054327719506772], [1.25, 1.0], [1.5, 0.0], [1.5, 0.31], [1.5, 0.62], [1.5, 0.81], [1.5, 1.0], [1.75, 0.0], [1.75, 0.3054327719506772], [1.75, 0.6108655439013544], [1.75, 0.8054327719506772], [1.75, 1.0], [2.0, 0.0], [2.0, 0.29242640687119287], [2.0, 0.5848528137423857], [2.0, 0.7924264068711928], [2.0, 1.0], [2.25, 0.0], [2.25, 0.2729610059419054], [2.25, 0.5459220118838108], [2.25, 0.7729610059419054], [2.25, 1.0], [2.5, 0.0], [2.5, 0.25], [2.5, 0.5], [2.5, 0.75], [2.5, 1.0], [2.75, 0.0], [2.75, 0.25], [2.75, 0.5], [2.75, 0.75], [2.75, 1.0], [3.0, 0.0], [3.0, 0.25], [3.0, 0.5], [3.0, 0.75], [3.0, 1.0]], "exploded": [[0.0, 0.0], [0.25, 0.0], [0.25, 0.25], [0.0, 0.25], [0.25, 0.5], [0.0, 0.5], [0.25, 0.75], [0.0, 0.75], [0.25, 1.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.25], [0.5, 0.5], [0.5, 0.75], [0.5, 1.0], [0.75, -0.09], [0.75, 0.18296100594190537], [0.75, 0.5459220118838107], [0.75, 0.7729610059419054], [0.75, 1.0], [1.0, -0.09], [1.0, 0.20242640687119287], [1.0, 0.5848528137423857], [1.0, 0.7924264068711928], [1.0, 1.0], [1.25, -0.09], [1.25, 0.2154327719506772], [1.25, 0.6108655439013544], [1.25, 0.8054327719506772], [1.25, 1.0], [1.5, -0.09], [1.5, 0.22], [1.5, 0.62], [1.5, 0.81], [1.5, 1.0], [1.75, -0.09], [1.75, 0.2154327719506772], [1.75, 0.6108655439013544], [1.75, 0.8054327719506772], [1.75, 1.0], [2.0, -0.09], [2.0, 0.20242640687119287], [2.0, 0.5848528137423857], [2.0, 0.7924264068711928], [2.0, 1.0], [2.25, -0.09], [2.25, 0.18296100594190542], [2.25, 0.5459220118838108], [2.25, 0.7729610059419054], [2.25, 1.0], [2.5, 0.0], [2.5, 0.25], [2.5, 0.5], [2.5, 0.75], [2.5, 1.0], [2.75, 0.0], [2.75, 0.25], [2.75, 0.5], [2.75, 0.75], [2.75, 1.0], [3.0, 0.0], [3.0, 0.25], [3.0, 0.5], [3.0, 0.75], [3.0, 1.0], [0.75, 0.45592201188381076], [1.0, 0.49485281374238577], [1.25, 0.5208655439013544], [1.5, 0.53], [1.75, 0.5208655439013544], [2.0, 0.49485281374238577], [2.25, 0.4559220118838109]], "tris": [[0, 1, 2], [0, 2, 3], [3, 2, 4], [3, 4, 5], [5, 4, 6], [5, 6, 7], [7, 6, 8], [7, 8, 9], [1, 10, 11], [1, 11, 2], [2, 11, 12], [2, 12, 4], [4, 12, 13], [4, 13, 6], [6, 13, 14], [6, 14, 8], [10, 15, 16], [10, 16, 11], [11, 16, 17], [11, 17, 12], [12, 17, 18], [12, 18, 13], [13, 18, 19], [13, 19, 14], [15, 20, 21], [15, 21, 16], [16, 21, 22], [16, 22, 17], [17, 22, 23], [17, 23, 18], [18, 23, 24], [18, 24, 19], [20, 25, 26], [20, 26, 21], [21, 26, 27], [21, 27, 22], [22, 27, 28], [22, 28, 23], [23, 28, 29], [23, 29, 24], [25, 30, 31], [25, 31, 26], [26, 31, 32], [26, 32, 27], [27, 32, 33], [27, 33, 28], [28, 33, 34], [28, 34, 29], [30, 35, 36], [30, 36, 31], [31, 36, 37], [31, 37, 32], [32, 37, 38], [32, 38, 33], [33, 38, 39], [33, 39, 34], [35, 40, 41], [35, 41, 36], [36, 41, 42], [36, 42, 37], [37, 42, 43], [37, 43, 38], [38, 43, 44], [38, 44, 39], [40, 45, 46], [40, 46, 41], [41, 46, 47], [41, 47, 42], [42, 47, 48], [42, 48, 43], [43, 48, 49], [43, 49, 44], [45, 50, 51], [45, 51, 46], [46, 51, 52], [46, 52, 47], [47, 52, 53], [47, 53, 48], [48, 53, 54], [48, 54, 49], [50, 55, 56], [50, 56, 51], [51, 56, 57], [51, 57, 52], [52, 57, 58], [52, 58, 53], [53, 58, 59], [53, 59, 54], [55, 60, 61], [55, 61, 56], [56, 61, 62], [56, 62, 57], [57, 62, 63], [57, 63, 58], [58, 63, 64], [58, 64, 59]], "moved_tris": [[0, 1, 2], [0, 2, 3], [3, 2, 4], [3, 4, 5], [5, 4, 6], [5, 6, 7], [7, 6, 8], [7, 8, 9], [1, 10, 11], [1, 11, 2], [2, 11, 12], [2, 12, 4], [4, 12, 13], [4, 13, 6], [6, 13, 14], [6, 14, 8], [10, 15, 16], [10, 16, 11], [11, 16, 65], [11, 65, 12], [12, 17, 18], [12, 18, 13], [13, 18, 19], [13, 19, 14], [15, 20, 21], [15, 21, 16], [16, 21, 66], [16, 66, 65], [17, 22, 23], [17, 23, 18], [18, 23, 24], [18, 24, 19], [20, 25, 26], [20, 26, 21], [21, 26, 67], [21, 67, 66], [22, 27, 28], [22, 28, 23], [23, 28, 29], [23, 29, 24], [25, 30, 31], [25, 31, 26], [26, 31, 68], [26, 68, 67], [27, 32, 33], [27, 33, 28], [28, 33, 34], [28, 34, 29], [30, 35, 36], [30, 36, 31], [31, 36, 69], [31, 69, 68], [32, 37, 38], [32, 38, 33], [33, 38, 39], [33, 39, 34], [35, 40, 41], [35, 41, 36], [36, 41, 70], [36, 70, 69], [37, 42, 43], [37, 43, 38], [38, 43, 44], [38, 44, 39], [40, 45, 46], [40, 46, 41], [41, 46, 71], [41, 71, 70], [42, 47, 48], [42, 48, 43], [43, 48, 49], [43, 49, 44], [45, 50, 51], [45, 51, 46], [46, 51, 52], [46, 52, 71], [47, 52, 53], [47, 53, 48], [48, 53, 54], [48, 54, 49], [50, 55, 56], [50, 56, 51], [51, 56, 57], [51, 57, 52], [52, 57, 58], [52, 58, 53], [53, 58, 59], [53, 59, 54], [55, 60, 61], [55, 61, 56], [56, 61, 62], [56, 62, 57], [57, 62, 63], [57, 63, 58], [58, 63, 64], [58, 64, 59]], "side": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], "chain": [12, 17, 22, 27, 32, 37, 42, 47, 52], "tips": [12, 52], "interior": [17, 22, 27, 32, 37, 42, 47], "replicas": {"17": 65, "22": 66, "27": 67, "32": 68, "37": 69, "42": 70, "47": 71}, "delta": 0.09} \ No newline at end of file diff --git a/articles/putting-a-fault-in-a-mesh/figures/fault-anatomy.png b/articles/putting-a-fault-in-a-mesh/figures/fault-anatomy.png new file mode 100644 index 0000000..c24ad25 Binary files /dev/null and b/articles/putting-a-fault-in-a-mesh/figures/fault-anatomy.png differ diff --git a/articles/putting-a-fault-in-a-mesh/figures/s_fault_geometry.png b/articles/putting-a-fault-in-a-mesh/figures/s_fault_geometry.png new file mode 100644 index 0000000..0fa79c3 Binary files /dev/null and b/articles/putting-a-fault-in-a-mesh/figures/s_fault_geometry.png differ diff --git a/articles/putting-a-fault-in-a-mesh/figures/sf_note_stress_slip.png b/articles/putting-a-fault-in-a-mesh/figures/sf_note_stress_slip.png new file mode 100644 index 0000000..31423fe Binary files /dev/null and b/articles/putting-a-fault-in-a-mesh/figures/sf_note_stress_slip.png differ diff --git a/articles/putting-a-fault-in-a-mesh/metadata.yml b/articles/putting-a-fault-in-a-mesh/metadata.yml new file mode 100644 index 0000000..b93ef7f --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/metadata.yml @@ -0,0 +1,33 @@ +# Validated in CI against schemas/article-metadata.schema.json. +# `pixi run validate` checks this and the cross-file invariants a schema cannot +# express -- that the article file is named .md, that canonical_path +# matches the slug, and that no legacy DOI is ever paired with a new registrant. +id: UWTN 2026-017 +slug: putting-a-fault-in-a-mesh +title: "Faults: to mesh or not to mesh?" +article_type: technical-note +status: published +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliation: Australian National University + - name: Thyagarajulu Gollapalli + orcid: 0000-0001-9394-4104 + affiliation: Australian National University +publication_date: 2026-09-18 +version: 1.0.0 +# The deposit writes archive_doi and repository_record_id when the note is +# published; leave them out until then. `doi` and `doi_registrant` were here +# once and are not fields the schema knows -- every note made from this +# template failed `pixi run validate` on all three of them. +license: CC-BY-4.0 +canonical_path: /putting-a-fault-in-a-mesh/ +legacy_paths: [] +# Facets, from vocabulary.yml. Both keys must be present even when empty: a +# note with no subject is normal -- many are purely about method. +subjects: +methods: +ghost_tags: + - Underworld Code +figures: 3 +source: native diff --git a/articles/putting-a-fault-in-a-mesh/putting-a-fault-in-a-mesh.md b/articles/putting-a-fault-in-a-mesh/putting-a-fault-in-a-mesh.md new file mode 100644 index 0000000..07567b6 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/putting-a-fault-in-a-mesh.md @@ -0,0 +1,154 @@ +--- +title: "Faults: to mesh or not to mesh?" +description: >- + A fault is a discontinuity and a mesh represents continuity, so the first + question is whether the fault has to go into the mesh at all. Shear bands and + a painted director cost no mesh work; cutting and conforming cost a great + deal. What each of the four representations asks of the mesh and of the + constitutive model, where they agree, and why they part company at a junction. +date: 2026-09-18 +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliations: + - Australian National University + - name: Thyagarajulu Gollapalli + orcid: 0000-0001-9394-4104 + affiliations: + - Australian National University +license: CC-BY-4.0 +bibliography: + - references.bib +keywords: + - Underworld Code + - Tricks of the Trade + - meshing +exports: + - format: typst + logo: ../../static/uwtn-logo.png + series: "Underworld Technical Notes" + origin_url: https://www.underworldcode.org/putting-a-fault-in-a-mesh/ + template: ../../templates/pdf + output: putting-a-fault-in-a-mesh.pdf + article_id: UWTN 2026-017 + article_version: 1.0.0 + software_version: underworld3 0.0.0 +--- +Faults can dominate the dynamic behaviour of Earth systems at scales from the entire planet to a few 10s of metres. They are an extreme example of localisation: feedback between forcing and response that results in self-reinforcing weakening that does not have a brake at the largest scale. Faults are *extreme* in the sense that their natural thickness is orders of magnitude below that of a typical tectonic simulation. + +At the tectonic scale, a fault is an infinitesimally thin surface across which the rock moves discontinuously by overcoming a frictional resistance. A finite element mesh (the mesh we use in Underworld) is a mechanism for representing continuous fields and there is not a native mechanism that perfectly represents a fault. + +There are lots of different potential solutions to this difficulty. They include: 1) ignoring it by using a continuum model of the fault, 2) adding extra interpolation functions that represent discontinuities, 3) splitting the mesh along the line of the fault and dealing with it as a surface, 4) refining the mesh enough that the fault width is invisibly small at the model length scale. + + +## Faults in numerical models + +Localisation is something a rheological law produces all by itself: if we give the material a yield stress or a strain-rate-weakening viscosity, shear bands appear where the stress requires them, with an orientation determined by the stress, and at whatever width the physics and the mesh between them allow. Do we really need to do anything more than this to represent faults ? The answer is yes and there are two main reasons. + +First, a fault is not in the same category as a shear band or a damage-zone. At the lithospheric scale a fault is persistent through changes in the tectonic loading. Faults localise far more sharply than any band a lithosphere-scale mesh resolves, down to a gouge zone of metres or even less. They are self-reinforcing: once a fault has accumulated slip it has juxtaposed distinct rock units; the weakness becomes structural and persistent even when the load is absent. Faults have history. + +Second, conceptually, at the tectonic scale fault is an infinitesimally thin surface across which the rock moves discontinuously by overcoming a frictional resistance. A finite element mesh (the mesh we use in Underworld) is a mechanism for representing continuous fields and there is not a native mechanism that perfectly represents a fault. Faults are sub-grid objects with their own constitutive properties. + +That is why it is common to impose plate boundaries in a mantle model as prior knowledge [@Davies_1988], and why those plate boundaries often have additional evolution rules. + +It is also why crustal models of stress build-up and release always try to include known faults. They are very fine-scale structures, they reflect complex geological history and they are not simply emergent from the imposed loading. + + +## Describe the fault not the implementation + +The fault geometry and its numerical representation ought to be decoupled as far as possible. +It is the thesis of both @Zhong_1995 and @Sandiford_2019: what the fault surface is, and what the solver is asked to do with it, are distinct choices. The first is a statement about the Earth: this is the fault surface, sampled as a polyline in two dimensions or a triangulated sheet in three, curved as it likes. The second is a modelling choice: this surface is to become a slippery interface, or a weak band of a stated width, or a direction of easy shear painted into the rheology. The same surface supports all of them and is agnostic to the implementation. + + +```{figure} figures/s_fault_geometry.png +:width: 60% +:name: fig-s-fault-geometry +:alt: A square domain with a red fault trace running from lower left to upper right, gently S-bent in the middle, drawn inside a pale blue ribbon labelled w = 0.03 and annotated main (tanh S); the trace continues as a red dashed line out of each corner, and the whole region above and left of it is shaded beige and labelled STRONG (eta x contrast) against a white region below and right. Beyond the main trace's upper tip a short collinear segment carries on in its own ribbon, labelled main segment (stepover) gap = 0.010. Just past the bend, three short parallel strands in pale green ribbons climb away into the beige region at a shallow angle to the main trace, labelled splay (kissing) + en-echelon zone; the lowest of the three almost touches the main trace, labelled Y gap = 0.010, and the other two are stepped up and to the left of it. To the right of the bend a straight dark red trace in its own blue ribbon, labelled branch (through-line), runs parallel to the strike and stops short of the main trace, marked gap = 0.071. Near the lower tip a short segment parallel to the main trace sits just below it in the white region, labelled lower stepover (offset 0.035). Every trace ends in a black dot. Two black arrows, one in each half of the domain, point up-right and down-left. + +A synthetic fault network that we use to validate the different fault algorithms. +``` + +[](#fig-s-fault-geometry) shows a synthetic fault network for our 2D and 3D experiments. The network is described by a number of independently meshed vertical segments in 3D and their surface traces for 2D models. It has a number of characteristics designed to test the fault implementation: multiple faults in a single domain; curved faults; junctions; steps and segment breaks for a single fault. There is a material property jump across the main fault in this network. The system is driven by boundary shear. + + +## Three implementations + +There is no perfect choice in how we model a fault in a geodynamic context and that means keeping several possible choices on hand to see which one works best for a specific problem. +Underworld3 offers three: a weak zone, a weak zone with a direction, and a cut. Each of them can be built on the unstructured mesh, or on a mesh modified to conform to the fault. + +| | Mesh untouched | Mesh conforms to the fault | +|---|---|---| +| Weak zone | painted band — the mesh sets the width | ribbon — the width is prescribed | +| TI weak zone ($\eta_1$) | painted, director from the distance gradient | TI ribbon — the width is prescribed | +| Cut | XFEM: enrichment functions carry the jump through the elements (not in Underworld) | split nodes with additional degrees of freedom | + +The **non-conforming, transversely isotropic (TI)** rheology requires no changes to the mesh to represent the mechanics of the fault. It simply creates a near-fault band of material that has a lower frictional strength parallel to the fault: the material there resists shear on the fault plane with a second, smaller viscosity $\eta_1$, and everything else with the background viscosity $\eta_0$. There is a zone of weakness that depends on the perpendicular distance to the *fault object*, and the internal orientation is given by the perpendicular vector to the fault surface. The band is whichever cells fall within $w/2$ of the fault, so its width is set by the mesh as much as by $w$. The resulting zone of weakness crosses element boundaries and can produce anomalous stress concentrations along the fault. These are mitigated at the large scale by refining the triangulation and smoothing the fault's influence function but they never disappear [@Yang_2021]. This is the representation of @Sharples_2015 and of @Sandiford_2019. + +Material property jumps can be incorporated into finite element representations when they lie along element boundaries. It is therefore possible to represent a fault as a rheologically distinct volume if we are prepared to remesh. The **weak ribbon** is a band of width $w$ meshed along the fault trace — its vertices are the fault's own points offset by $w/2$ either side — with a contrasting rheology in that region (a weak zone, or a zone with its own plasticity coefficients). For this we do need the ability to remesh so it is more difficult to implement in models where the fault system evolves. In this model, the fault, as a volume under normal stress, can deform internally and violate the frictional surface *approximation*. This is primarily an issue when $w$ is significantly larger than the fault's true physical width. + +The latter problem can be alleviated by combining the first two approaches. A **TI ribbon** is the same band but using the transversely isotropic frictional model of @Moresi_2006 within the remeshed band: the same rheology as the non-conforming case, with the remesh taking control of the width. This has the same desirable properties from the finite element solver's point of view as the weak band model, but it can transmit normal stresses across the fault without internal flow. + + +The **split** algorithm is the one that embraces the notion of the fault as an embedded surface. It produces a cut through the mesh to create a new internal surface boundary. In general this is also a remeshing step: the mesh is cut along the fault line, dividing the elements it crosses, and degrees of freedom are added along the cut; @Zhong_1995 put slippery nodes into a convection model this way, on a hexahedral grid nudged towards the fault. We make the mesh conform to the fault beforehand, so that what remains is the duplication alone. Faults are surfaces that conform to element boundaries but they are implemented as pairs of surfaces to represent the two sides of the fault. The constitutive model lies in the interaction of these two coincident surfaces. This is a good choice of model when the physical scale completely precludes resolving $w$, but there are some limitations: because the mesh is cut into sliding surfaces, there are incompatible constraints when two faults meet or cross. + + +```{figure} figures/fault-anatomy.png +:width: 72% +:name: fig-fault-anatomy +:alt: Three panels, each the same rectangular triangulation twelve cells across and four deep, pale grey, stacked vertically. (a) The grid is flat and a red arch is drawn across the middle two thirds without dots, cutting through the triangles; the tinted green cells form a ragged band around it, two rows deep on the flanks and three at the crest, with a saw-toothed outline, each with a short dark-green stroke tilted perpendicular to the arch. (b) The grid is gently bowed upward in its middle rows so that the same red arch, now with nine dots, runs along mesh edges; the cells immediately above and below it are tinted green, each with a stroke perpendicular to the local trace; a bracket at the right marks the two-cell band height as w. (c) The bowed grid again, with the row of cells above the arch tinted pale blue and the row below pale pink; the pink block has dropped, so the red line has opened into two — a solid upper arch with seven filled dots labelled v-plus (original) and Gamma-plus, and a dashed lower arch with open circles labelled v-minus (replica) and Gamma-minus — still meeting at a black ringed vertex labelled tip at each end, with the cells at the two ends sheared where the block has dropped. + +One fault, three strategies, on one mesh. (a) The non-conforming paint: the trace across the flat grid, and the band is whichever cells fall within $w/2$ of it. (b) The ribbon: the mesh is bent so that the same trace runs along element edges, and the band is the cells either side of it, each carrying a director; nothing is duplicated and every field is continuous. (c) The split, on the bent mesh: each interior vertex of the chain is duplicated and the cells on the Minus side are rewired to the replica; the tips are not duplicated. The copies are coincident — the lower block is pulled away only so that they can be seen. +``` + +[](#fig-fault-anatomy) compares the three strategies on one mesh. Panel (a) shows the fault running through the mesh with elements identified as fault / not-fault depending on their centroid distance to the fault. Panel (b) shows the ribbon near the fault constructed so that the fault volume itself is defined by element boundaries. Panel (c) is what the split does. The trace has first been made a chain of element edges, so that the two cells at every facet share it and every field is continuous across it, as anywhere else in the mesh. The split then duplicates each interior vertex of the chain. The original stays with the cells on one side, which we label Plus; the cells on the other side, Minus, are rewired to a replica at the same position. In this simple example, no elements are added or divided, but, across the fault, the two sides no longer share a degree of freedom, so the velocity is free to jump. + +The fault is just the pair of surfaces and the condition we impose between each vertex and its replica. The simplest condition is no-opening: the two normal velocities of a pair are equal and the tangential velocities are free, which is a frictionless slippery interface. The slip rate is read off the pair as the tangential jump. A friction law is a relation between that jump and the traction the pair carries, and it lives on the pair as well. + +The two end vertices of the chain are not duplicated. Slip therefore goes to zero at the tips, which is the crack condition, and the front of the fault needs no treatment of its own. This is also a limitation: a vertex that belonged to two chains would be a tip of each, pinned on both, so two cuts cannot meet. + +## Difficulties with branching faults + +```{figure} figures/sf_note_stress_slip.png +:width: 64% +:name: fig-stress-and-slip +:alt: A three-by-two grid of panels; the left column is the cut, the right the band, both at w = 0.005. The top row shows the whole square domain coloured by log10 of the second stress invariant on a black-purple-orange-white scale from -1.18 to 0.99: a flat mid-orange background, a dark low-stress lobe flanking each trace, and a bright concentration at every tip. Each trace is a tube coloured by its signed slip rate on a blue-grey-green scale, blue sinistral down to -0.029 and green dextral up to 0.54: the S-bent main trace is dark green, the branch and the lower stepover segment mid green, the three en-echelon strands pale, and the splay nearest the main carries a short blue reach where it meets it. The two columns look the same at this scale. The middle row zooms the Y junction with the mesh in faint white: in the cut the splay's tube stops a short distance from the main trace and the stress there is smooth; in the band it runs into the main trace and the two merge, with a bright knot at the join and dark lobes either side. The bottom row zooms the upper stepover, where the main trace ends and a collinear segment carries on one element beyond it. In the cut, a bright four-lobed stress concentration sits on the one-element bridge between the two tips, the continuation segment's tube is noticeably paler than the main trace's, and a dark lobe lies along its far side. In the band, one green tube runs straight through with no break and no change of shade, a narrow dark band runs along it, and there is no concentration at the old tip at all — only a single dark speck on the line where the two tips were. + +The same network as a cut (left) and as a band (right), both at $w = 0.005$ with $\eta_1/w = 0.1$. Background: log10 of the second stress invariant on one scale. Traces: each representation's own slip rate. Below, the Y junction, and the collinear stepover closed to one element — the least gap a cut can leave. +``` + +Consider one un-branched strand by itself. If we look at the slip rate across the fault for the transversely isotropic ribbon and compare it to the split mesh, we find that the two representations converge as the fault ribbon shrinks in width (in a background mesh of fixed resolution), provided the ratio $\eta_1/w$ remains fixed. The band's mechanical strength is the ratio $\eta_1/w$, not $\eta_1$: halving the width at fixed viscosity effectively doubles the interface strength (see table). + +| Main strand alone, peak slip rate | $w = 0.03$ | $w = 0.01$ | $w = 0.005$ | +|---|---|---|---| +| split | 0.5149 | 0.5154 | 0.5152 | +| TI band, $\eta_1/w = 0.1$ | — | 0.5232 | 0.5133 | +| TI band, $\eta_1 = 10^{-3}$ fixed | — | 0.5232 | 0.5000 | + + +Branching of the fault breaks this convergence because a junction is the one place where the two representations describe genuinely different objects. A band can fork: two weak zones meet and merge into one continuous weak region. Two cuts cannot: if they touch they would share a node that would carry two incompatible sets of constraints. The branch in the model is therefore represented slightly differently. The rheological bands merge smoothly in the mesh, but the cut branch stops short of the main fault. + +The band's junction is not free of choices either: in the cells the two bands share there can be only one director, and it is the orientation of whichever band was painted last. + +Where the branches touch, slip is significantly higher, with a noticeable reversal of polarity in the slip orientation along the branch. The effect is present but more muted for the cut-in branch that does not quite reach the main fault. + +The bottom row of [](#fig-stress-and-slip) is the same difficulty in its simplest form: two collinear segments of one fault, butted end to end with a single element between them. The cut has no choice about that element — two cuts that shared a vertex would pin it, so one element of intact rock is the least gap it can leave — and that ligament welds the two segments together. The segment beyond the gap carries two thirds of the slip a continuous fault carries there, the weld draws the main segment down by a few per cent as well, and refinement halves the weld's length without ever removing it. The band has no such constraint: the two weak zones overlap and become one, the slip runs through unbroken, and the stress concentration at the buried tips disappears with it. Where a fault is segmented on the scale of the mesh rather than the scale of the model, this is the case for the band. + + +## Choices + +The split-node approach is by far the most efficient representation of a discontinuous, frictional fault that works well when the fault-width is completely unreachable with meshing. It does require cutting into the mesh each time the mesh is adapted or the fault is moved. The surface conditions are well-behaved when the solver sees them. + +The transversely isotropic, meshed ribbon does a good job of faults that have branching structures or multiple, overlapping segments (or even segments that butt against each other). It requires mesh adaptation to keep the solver happy, and there is some tuning required to ensure the implementation converges to the split-node formulation. + +The non-conforming transversely isotropic fault representation trades fault fidelity and solver efficiency against simplicity. This is the choice for cases where remeshing or mesh adaptation is difficult, and fluctuations in the near-fault stress-field can be tolerated. + + + +
Comments
Discussion of these notes happens in GitHub Discussions, so it stays with the source and is searchable alongside it.
Read the discussionStart one
diff --git a/articles/putting-a-fault-in-a-mesh/references.bib b/articles/putting-a-fault-in-a-mesh/references.bib new file mode 100644 index 0000000..71327e5 --- /dev/null +++ b/articles/putting-a-fault-in-a-mesh/references.bib @@ -0,0 +1,96 @@ +% Resolved from doi.org and pinned here on purpose. Letting MyST fetch these at +% build time made the build depend on doi.org answering: it worked locally and +% failed on CI with "Citation data from doi.org was not available or malformed" +% (Zhong & Gurnis and Sandiford & Moresi on 2026-09-11, Yang et al. on +% 2026-09-18), which would have published the note with a broken citation and +% no reference entry. A deposited PDF cannot be repaired after the fact, so the +% bibliography is part of the source. + +@article{Davies_1988, + title = {Role of the lithosphere in mantle convection}, + volume = {93}, + ISSN = {0148-0227}, + url = {http://dx.doi.org/10.1029/JB093iB09p10451}, + DOI = {10.1029/JB093iB09p10451}, + number = {B9}, + journal = {Journal of Geophysical Research: Solid Earth}, + publisher = {American Geophysical Union (AGU)}, + author = {Davies, Geoffrey F.}, + year = {1988}, + month = sep, + pages = {10451--10466} +} + +@article{Zhong_1995, + title = {Mantle Convection with Plates and Mobile, Faulted Plate Margins}, + volume = {267}, + ISSN = {1095-9203}, + url = {http://dx.doi.org/10.1126/science.267.5199.838}, + DOI = {10.1126/science.267.5199.838}, + number = {5199}, + journal = {Science}, + publisher = {American Association for the Advancement of Science (AAAS)}, + author = {Zhong, Shijie and Gurnis, Michael}, + year = {1995}, + month = feb, + pages = {838--843} +} + +@article{Moresi_2006, + title = {Anisotropic viscous models of large-deformation Mohr–Coulomb failure}, + volume = {86}, + ISSN = {1478-6443}, + url = {http://dx.doi.org/10.1080/14786430500255419}, + DOI = {10.1080/14786430500255419}, + number = {21-22}, + journal = {Philosophical Magazine}, + publisher = {Informa UK Limited}, + author = {Moresi, L. and Mühlhaus, H.-B.}, + year = {2006}, + month = jul, + pages = {3287--3305} +} + +@article{Sharples_2015, + title = {Styles of rifting and fault spacing in numerical models of crustal extension}, + volume = {120}, + ISSN = {2169-9356}, + url = {http://dx.doi.org/10.1002/2014JB011813}, + DOI = {10.1002/2014JB011813}, + number = {6}, + journal = {Journal of Geophysical Research: Solid Earth}, + publisher = {American Geophysical Union (AGU)}, + author = {Sharples, W. and Moresi, L.-N. and Jadamec, M. A. and Revote, J.}, + year = {2015}, + month = jun, + pages = {4379--4404} +} + +@article{Sandiford_2019, + title = {Improving subduction interface implementation in dynamic numerical models}, + volume = {10}, + ISSN = {1869-9529}, + url = {http://dx.doi.org/10.5194/se-10-969-2019}, + DOI = {10.5194/se-10-969-2019}, + number = {3}, + journal = {Solid Earth}, + publisher = {Copernicus GmbH}, + author = {Sandiford, Dan and Moresi, Louis}, + year = {2019}, + month = jun, + pages = {969--985} +} + +@article{Yang_2021, + title = {Stress recovery for the particle-in-cell finite element method}, + volume = {311}, + ISSN = {0031-9201}, + url = {http://dx.doi.org/10.1016/j.pepi.2020.106637}, + DOI = {10.1016/j.pepi.2020.106637}, + journal = {Physics of the Earth and Planetary Interiors}, + publisher = {Elsevier BV}, + author = {Yang, Haibin and Moresi, Louis N. and Mansour, John}, + year = {2021}, + month = feb, + pages = {106637} +} diff --git a/authors.yml b/authors.yml index ec7e399..11929b7 100644 --- a/authors.yml +++ b/authors.yml @@ -50,7 +50,8 @@ authors: gthyagi: name: Thyagarajulu Gollapalli orcid: 0000-0001-9394-4104 - affiliation: Monash University + # At ANU since the start of 2026; the JOSS paper predates the move. + affiliation: Australian National University knepley: name: Matt Knepley orcid: 0000-0002-2292-0735 diff --git a/classification.yml b/classification.yml index 4925dd5..60749b4 100644 --- a/classification.yml +++ b/classification.yml @@ -318,6 +318,11 @@ introducing-the-technical-notes: subjects: [] methods: [] +putting-a-fault-in-a-mesh: + article_type: technical-note + subjects: [earthquakes-faults, tectonics-lithosphere] + methods: [meshing, rheology] + setting-up-full-multigrid: article_type: technical-note subjects: [mantle-convection] diff --git a/scripts/build_pdf.py b/scripts/build_pdf.py index b7e808c..b794b8b 100644 --- a/scripts/build_pdf.py +++ b/scripts/build_pdf.py @@ -31,6 +31,44 @@ def run(*command): return subprocess.call(list(command), cwd=ROOT) +def warm_package_cache(): + """Fetch every Typst package the PDF template imports, once, serially. + + Typst downloads a package into a per-user cache the first time any + compile imports it, and MyST compiles the archival PDFs concurrently. + On a fresh CI runner that cache is empty, so every compile that starts + together downloads the same package into the same directory, and the + losers die with "failed to decompress package" or "package not found" + having watched somebody else's download reach 100%. Two contenders got + away with it for months; when the template itself began importing + @preview/tablex (the house table style), every article became a + contender and four PDFs died on the first cold runner (PR #54). + + The template's own imports are ours to know, so they are read from the + template files rather than listed here, and compiled from a stub before + MyST starts. MyST's generated preamble may import more (subpar, for + subfigures); the serial first-target build below still covers those. + """ + import re + import tempfile + pattern = re.compile(r'@preview/[A-Za-z0-9_-]+:[0-9]+\.[0-9]+\.[0-9]+') + packages = set() + for path in sorted((ROOT / "templates" / "pdf").glob("*.typ")): + packages.update(pattern.findall(path.read_text(encoding="utf-8"))) + if not packages: + return + with tempfile.TemporaryDirectory() as tmp: + stub = pathlib.Path(tmp) / "warm.typ" + stub.write_text("".join('#import "%s"\n' % p for p in sorted(packages)) + + "warm\n", encoding="utf-8") + status = subprocess.call( + ["typst", "compile", str(stub), str(pathlib.Path(tmp) / "warm.pdf")], + cwd=ROOT) + print("typst package cache warmed: %s" % ", ".join(sorted(packages)) + if status == 0 else + "typst package cache NOT warmed (exit %d); the build may race" % status) + + def main(): if run(sys.executable, "scripts/sync_archival.py") != 0: sys.exit("could not sync the archival metadata") @@ -76,6 +114,7 @@ def main(): # # Seen on PR #17, the first build to produce two PDFs at once. It is luck # rather than design that the 43-PDF production build has not hit it. + warm_package_cache() status = 0 if len(targets) > 1: status = run("myst", "build", "--typst", targets[0]) diff --git a/static/uwtn.css b/static/uwtn.css index 710aae6..5ba2354 100644 --- a/static/uwtn.css +++ b/static/uwtn.css @@ -12,6 +12,8 @@ --uwtn-muted: #5d6b7a; --uwtn-rule: #dfe4ea; --uwtn-accent: #1a4f80; /* the PDF's theme colour, blue.darken(30%) */ + --uwtn-table-head: #1a4f80; /* table header fill: the theme blue, both modes */ + --uwtn-table-head-ink: #dbe9f7; --uwtn-paper: #ffffff; --uwtn-serif: "Iowan Old Style", "Palatino Linotype", Palatino, Georgia, "Times New Roman", serif; @@ -25,6 +27,8 @@ --uwtn-rule: #2b3541; --uwtn-accent: #7cb3e0; --uwtn-paper: #131a21; + --uwtn-table-head: #16456f; + --uwtn-table-head-ink: #d3e4f5; } } @@ -213,6 +217,39 @@ article.content h3 { letter-spacing: -0.005em; } +/* ---- tables ------------------------------------------------------------ */ + +/* The house table, matching the PDF template's tableStyle: centred, sans, + horizontal hairlines only, the header row filled in the theme blue with + light text. MyST emits a bare inside a wrapper div, so the rules + key on the elements. */ +article.content table { + margin: 1.4rem auto; + width: auto; + border-collapse: collapse; + font-family: var(--uwtn-sans); + font-size: 0.9rem; + line-height: 1.4; +} + +article.content th { + background: var(--uwtn-table-head); + color: var(--uwtn-table-head-ink); + font-weight: 600; + text-align: left; + padding: 0.45rem 0.85rem; + border: none; +} + +article.content td { + padding: 0.4rem 0.85rem; + border: none; + border-bottom: 1px solid var(--uwtn-rule); +} + +article.content th code, +article.content td code { font-size: 0.85em; } + /* Captions read as apparatus, matching the PDF: sans, smaller, lighter. */ article.content figcaption, article.content figcaption p { diff --git a/templates/pdf/template.typ b/templates/pdf/template.typ index 06f993e..c382753 100644 --- a/templates/pdf/template.typ +++ b/templates/pdf/template.typ @@ -133,4 +133,29 @@ [-IMPORTS-] +// The house table style. MyST renders every markdown table as +// `tablex(columns: n, header-rows: 1, ..tableStyle, ..columnStyle, ...)` +// with both dictionaries EMPTY in its generated myst-imports.typ, so the +// template sets the style by shadowing them here, after that import and +// before the content. Horizontal hairlines only, no verticals; the header +// row filled in the theme blue with light text; the table centred. The +// site's CSS (static/uwtn.css, "tables") carries the same design to HTML. +#import "@preview/tablex:0.0.9": tablex as tablex-base +#let uwtn-theme = blue.darken(30%) +#let uwtn-table-sans = ("Helvetica Neue", "Helvetica", "Arial") +#let tableStyle = ( + auto-vlines: false, + auto-hlines: true, + map-hlines: h => (..h, stroke: 0.3pt + luma(165)), + fill: (col, row) => if row == 0 { uwtn-theme } else { none }, + map-cells: c => if c.y == 0 { + (..c, content: text(font: uwtn-table-sans, size: 8pt, weight: "semibold", + fill: blue.lighten(88%), c.content)) + } else { + (..c, content: text(size: 9pt, c.content)) + }, + inset: (x: 8pt, y: 4.5pt), +) +#let tablex(..args) = align(center, block(above: 10pt, below: 12pt, tablex-base(..args))) + [-CONTENT-]