Faults: to mesh or not to mesh? (UWTN 2026-017)

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

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)])
+})

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")

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)}")