Building/.work/wedge/draw_wedge_cases.py

"""RocSlope3 wedge-formation examples: which 2-plane sets form a block, which don't.

Frame: x = east, y = north, z = up. Azimuths from north, clockwise.
Slope: face dips 60 deg toward 180 (due south = -y). H = 8 m, toe at y=0, crest y=4.619.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

rd = np.radians
OK, BAD = "\u2713", "\u2717"   # check / cross

# ---------------------------------------------------------------- geometry ---
def normal(dip, dd):
    d, a = rd(dip), rd(dd)
    return np.array([np.sin(a) * np.sin(d), np.cos(a) * np.sin(d), np.cos(d)])

def clip_poly(poly, n, h):
    out = []
    for i in range(len(poly)):
        A, B = poly[i], poly[(i + 1) % len(poly)]
        da, db = float(n @ A) - h, float(n @ B) - h
        if da <= 0:
            out.append(A)
        if (da < 0) != (db < 0):
            t = da / (da - db)
            out.append(A + t * (B - A))
    return out

def plane_polygon(n, c, box, L=40.0):
    a = np.array([0.13, 0.52, 0.84])
    if abs(n @ a) > 0.9:
        a = np.array([0.9, 0.2, 0.35])
    u = np.cross(n, a); u /= np.linalg.norm(u)
    v = np.cross(n, u)
    pts = [c + L * (u + v), c + L * (-u + v), c + L * (-u - v), c + L * (u - v)]
    for axis, (lo, hi) in enumerate(box):
        e = np.zeros(3); e[axis] = 1.0
        pts = clip_poly(pts, e, hi)
        e = np.zeros(3); e[axis] = -1.0
        pts = clip_poly(pts, e, -lo)
    return np.array(pts)

def line_in_box(p0, d, box):
    t0, t1 = -1e9, 1e9
    for axis, (lo, hi) in enumerate(box):
        if abs(d[axis]) < 1e-12:
            continue
        ta, tb = (lo - p0[axis]) / d[axis], (hi - p0[axis]) / d[axis]
        ta, tb = min(ta, tb), max(ta, tb)
        t0, t1 = max(t0, ta), min(t1, tb)
    return p0 + t0 * d, p0 + t1 * d

def inter_point(n1, c1, n2, c2):
    A = np.vstack([n1, n2])
    p, *_ = np.linalg.lstsq(A, np.array([c1, c2]), rcond=None)
    return p

def isect_dir(n1, n2):
    L = np.cross(n1, n2)
    return L / np.linalg.norm(L)

def down_plunge(L):                # ensure z<0, return (plunge, trend)
    if L[2] > 0:
        L = -L
    plunge = np.degrees(np.arcsin(-L[2]))
    trend = np.degrees(np.arctan2(L[0], L[1])) % 360.0
    return L, plunge, trend

# slope constants
FDIP, FDD, H = 60.0, 180.0, 8.0
YC = H / np.tan(rd(FDIP))                      # crest y = 4.6188
pbox = ((-5.7, 5.7), (-1.2, 9.3), (0.15, 8.45))  # clip box for joint sheets

# ------------------------------------------------------------------ cases ---
# case 1: WORKS  J1 55/145, J2 55/215 through a toe point on the face
F1 = np.array([0.0, 1.8, np.tan(rd(FDIP)) * 1.8])
n1a, n1b = normal(55, 145), normal(55, 215)
c1a, c1b = n1a @ F1, n1b @ F1
La, pl1, tr1 = down_plunge(isect_dir(n1a, n1b))
T1 = F1 - ((8.0 - F1[2]) / -La[2]) * La        # walk up-plunge from F to the top surface
x1a = (c1a - n1a[1] * YC - n1a[2] * 8.0) / n1a[0]
C1a, C1b = np.array([x1a, YC, 8.0]), np.array([-x1a, YC, 8.0])
print("case1 isect plunge %.1f trend %.0f | F y=%.2f | T y=%.2f | crest x=+/-%.2f"
      % (pl1, tr1, F1[1], T1[1], abs(x1a)))

# case 2: FAIL steep  J 80/145 & 80/215 through crest points x=+/-1.9
n2a, n2b = normal(80, 145), normal(80, 215)
P2a, P2b = np.array([-1.9, YC, 8.0]), np.array([1.9, YC, 8.0])
c2a, c2b = n2a @ P2a, n2b @ P2b
L2, pl2, tr2 = down_plunge(isect_dir(n2a, n2b))
p2 = inter_point(n2a, c2a, n2b, c2b)
s2a, s2b = line_in_box(p2, L2, ((-6, 6), (-2.2, 11), (0.0, 8.0)))
print("case2 isect plunge %.1f trend %.0f | line %s -> %s" %
      (pl2, tr2, np.round(s2a, 2), np.round(s2b, 2)))

# case 3: FAIL dips into slope  J 55/035 & 55/325 through crest points x=+/-1.9
n3a, n3b = normal(55, 35), normal(55, 325)
P3a, P3b = np.array([-1.9, YC, 8.0]), np.array([1.9, YC, 8.0])
c3a, c3b = n3a @ P3a, n3b @ P3b
L3, pl3, tr3 = down_plunge(isect_dir(n3a, n3b))
p3 = inter_point(n3a, c3a, n3b, c3b)
s3a, s3b = line_in_box(p3, L3, ((-6, 6), (-2.2, 11), (0.0, 8.0)))
print("case3 isect plunge %.1f trend %.0f | line %s -> %s" %
      (pl3, tr3, np.round(s3a, 2), np.round(s3b, 2)))

# case 4: FAIL no physical intersection - small offset patches on case-1 planes
def patch(n, c, center, h=1.15):
    a = np.array([0.13, 0.52, 0.84])
    if abs(n @ a) > 0.9:
        a = np.array([0.9, 0.2, 0.35])
    u = np.cross(n, a); u /= np.linalg.norm(u)
    v = np.cross(n, u)
    return np.array([center + h * (u + v), center + h * (-u + v),
                     center + h * (-u - v), center + h * (u - v)])

PA = np.array([-1.8, 2.2, (c1a - n1a[0] * -1.8 - n1a[1] * 2.2) / n1a[2]])
PB = np.array([1.8, 2.2, (c1b - n1b[0] * 1.8 - n1b[1] * 2.2) / n1b[2]])
pa4, pb4 = patch(n1a, c1a, PA), patch(n1b, c1b, PB)
print("case4 patch centers", np.round(PA, 2), np.round(PB, 2))

# ---------------------------------------------------------------- drawing ---
BLUE, ORNG, RED, GRN = "#1f77b4", "#ff7f0e", "#c0392b", "#1e8449"
fig = plt.figure(figsize=(17.5, 9.8))

titles = [
    (OK + "  WORKS \u2014 wedge daylights in the face", GRN,
     "J1 55/145 + J2 55/215  \u2192  line of intersection 49/180,\nshallower than the 60/180 face \u2192 toe daylights mid-face"),
    (BAD + "  FAILS \u2014 intersection too steep", RED,
     "J1 80/145 + J2 80/215  \u2192  intersection 78/180 is steeper\nthan the face (60/180) \u2192 toe never reaches the face"),
    (BAD + "  FAILS \u2014 planes dip into the slope", RED,
     "J1 55/035 + J2 55/325  \u2192  intersection plunges backward\n(49/360, up-hill) \u2192 tapered wedge, no release"),
    (BAD + "  FAILS \u2014 planes never meet", RED,
     "orientations like case 1, but finite planes too small / offset\n\u2192 no intersection line inside the slope \u2192 no closed block"),
]

for i in range(4):
    xc = (i + 0.5) / 4
    fig.text(xc, 0.955, titles[i][0], ha="center", va="top", fontsize=14,
             fontweight="bold", color=titles[i][1])
    fig.text(xc, 0.916, titles[i][2], ha="center", va="top", fontsize=10.2,
             color="#34495e")

fig.text(0.5, 0.999, "RocSlope3: when do two planes form a wedge block?   "
         "(slope face dips 60\u00b0 toward azimuth 180\u00b0 \u2014 dip / dip-direction given as \u03b4/\u03b1)",
         ha="center", va="top", fontsize=15, fontweight="bold")

# ------------------------------------------------------------- map views ----
for i in range(4):
    ax = fig.add_subplot(2, 4, i + 1)
    ax.set_xlim(-6.7, 7.0); ax.set_ylim(-3.0, 8.6)
    ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
    for s in ax.spines.values():
        s.set_edgecolor("#bbbbbb")
    ax.add_patch(plt.Rectangle((-6.2, YC), 12.9, 8.6 - YC, fc="#e4e8ec", ec="none", zorder=0))
    ax.add_patch(plt.Rectangle((-6.2, 0), 12.9, YC, fc="#ddd0b6", ec="none", zorder=0))
    ax.axhline(YC, color="#3c3c3c", lw=2.2, zorder=2)
    ax.axhline(0, color="#6b6b6b", lw=1.3, zorder=2)
    ax.text(6.3, YC + 0.15, "crest", fontsize=8.5, color="#3c3c3c")
    ax.text(-6.55, 0.15, "toe", fontsize=8.5, color="#6b6b6b", ha="left")
    ax.text(-6.45, 6.9, "rock", fontsize=9, color="#7a8288", ha="left")
    ax.text(-6.45, -2.0, "air", fontsize=9, color="#9aa1a6", ha="left")

    def face_arrow(x=0.0):
        ax.annotate("", xy=(x, -1.15), xytext=(x, 3.7),
                    arrowprops=dict(arrowstyle="-|>", lw=2.6, color="#222222"))
        ax.text(x + 0.35, -1.85, "face 60/180", fontsize=9.5,
                fontweight="bold", color="#222222", ha="left")

    def strike_trace(x0, dd, color):
        s = np.array([np.cos(rd(dd)), -np.sin(rd(dd))])
        p0 = np.array([x0, YC])
        a, b = p0 - 4.6 * s, p0 + 4.6 * s
        ax.plot([a[0], b[0]], [a[1], b[1]], ls=":", lw=1.1, color=color,
                alpha=0.75, zorder=2)

    def joint_arrow(x0, dip, dd, color, label, lab_off):
        v = np.array([np.sin(rd(dd)), np.cos(rd(dd))])
        p0 = np.array([x0, YC]); tip = p0 + 2.6 * v
        ax.annotate("", xy=tuple(tip), xytext=tuple(p0),
                    arrowprops=dict(arrowstyle="-|>", lw=2.5, color=color))
        ax.text(*(tip + np.array(lab_off)), label, fontsize=9.3, color=color,
                fontweight="bold", ha="center")

    def isect_arrow(trend, good, label):
        v = np.array([np.sin(rd(trend)), np.cos(rd(trend))])
        p0 = np.array([0.0, YC]); tip = p0 + 3.9 * v
        col = GRN if good else RED
        ax.annotate("", xy=tuple(tip), xytext=tuple(p0),
                    arrowprops=dict(arrowstyle="-|>", lw=3.2, color=col,
                                    linestyle=(0, (5, 3))))
        lx = 1.15 if abs(v[0]) < 0.3 else (0.9 * np.sign(v[0]))
        ly = tip[1] + (-0.55 if v[1] < 0 else 0.45)
        ax.text(tip[0] + lx, ly, label, fontsize=9.3, color=col,
                fontweight="bold", ha="left" if lx > 0 else "right")

    if i == 0:
        strike_trace(-3.2, 145, BLUE); strike_trace(3.2, 215, ORNG)
        joint_arrow(-3.2, 55, 145, BLUE, "J1 55/145", (-0.35, -0.75))
        joint_arrow(3.2, 55, 215, ORNG, "J2 55/215", (0.35, -0.75))
        isect_arrow(180, True, "intersection 49/180 " + OK)
        face_arrow()
        ax.annotate("", xy=(5.55, 8.15), xytext=(5.55, 7.15),
                    arrowprops=dict(arrowstyle="-|>", lw=1.8, color="#222222"))
        ax.text(5.55, 8.32, "N", fontsize=10, ha="center", fontweight="bold")
    elif i == 1:
        strike_trace(-3.2, 145, BLUE); strike_trace(3.2, 215, ORNG)
        joint_arrow(-3.2, 80, 145, BLUE, "J1 80/145", (-0.35, -0.75))
        joint_arrow(3.2, 80, 215, ORNG, "J2 80/215", (0.35, -0.75))
        isect_arrow(180, False, "78/180 " + BAD + "  >  face 60\u00b0")
        face_arrow()
    elif i == 2:
        strike_trace(-3.2, 35, BLUE); strike_trace(3.2, 325, ORNG)
        joint_arrow(-3.2, 55, 35, BLUE, "J1 55/035", (-0.45, 0.55))
        joint_arrow(3.2, 55, 325, ORNG, "J2 55/325", (0.45, 0.55))
        isect_arrow(360, False, "49/360 " + BAD + " backward")
        face_arrow()
    else:
        sB = np.array([np.cos(rd(145)), -np.sin(rd(145))])
        pB = np.array([-1.15, 2.0]); a, b = pB - 1.55 * sB, pB + 1.55 * sB
        ax.plot([a[0], b[0]], [a[1], b[1]], lw=4.5, color=BLUE, solid_capstyle="round")
        ax.text(-3.6, 1.1, "J1", fontsize=9.5, color=BLUE, fontweight="bold")
        sO = np.array([np.cos(rd(215)), -np.sin(rd(215))])
        pO = np.array([1.15, 2.0]); a, b = pO - 1.55 * sO, pO + 1.55 * sO
        ax.plot([a[0], b[0]], [a[1], b[1]], lw=4.5, color=ORNG, solid_capstyle="round")
        ax.text(3.6, 1.1, "J2", fontsize=9.5, color=ORNG, fontweight="bold")
        ax.add_patch(Ellipse((0, 2.0), 2.1, 2.7, fill=False, lw=2.0,
                             ec=RED, ls="--", zorder=3))
        ax.text(0, 3.75, "gap " + BAD + " planes never meet", fontsize=9.6,
                color=RED, fontweight="bold", ha="center")
        face_arrow(x=-4.9)

fig.text(0.5, 0.487, "map view (looking down; arrows = dip directions)", ha="center",
         fontsize=11, style="italic", color="#5d6d7e")

# ------------------------------------------------------------- 3d views -----
rock_faces = [
    ([(-6, 0, 0), (6, 0, 0), (6, YC, 8), (-6, YC, 8)], "#d8c9ad"),   # face
    ([(-6, YC, 8), (6, YC, 8), (6, 11, 8), (-6, 11, 8)], "#eef1f3"), # top
    ([(-6, 0, 0), (6, 0, 0), (6, 11, 0), (-6, 11, 0)], "#cfd4d8"),   # bottom
    ([(-6, 11, 0), (6, 11, 0), (6, 11, 8), (-6, 11, 8)], "#dfe3e6"), # back
    ([(-6, 0, 0), (-6, YC, 8), (-6, 11, 8), (-6, 11, 0)], "#dfe3e6"),
    ([(6, 0, 0), (6, YC, 8), (6, 11, 8), (6, 11, 0)], "#d7dbdf"),
]
rock_alpha = 0.16

wire = [  # crisp prism outline drawn on top
    [(-6, 0, 0), (6, 0, 0)], [(6, 0, 0), (6, YC, 8)], [(6, YC, 8), (-6, YC, 8)],
    [(-6, YC, 8), (-6, 0, 0)], [(-6, YC, 8), (-6, 11, 8)], [(6, YC, 8), (6, 11, 8)],
    [(-6, 11, 8), (6, 11, 8)], [(-6, 0, 0), (-6, 11, 0)], [(6, 0, 0), (6, 11, 0)],
    [(-6, 11, 0), (6, 11, 0)], [(-6, 11, 0), (-6, 11, 8)], [(6, 11, 0), (6, 11, 8)],
]

def add_3d(i):
    ax = fig.add_subplot(2, 4, i + 5, projection="3d")
    ax.set_xlim(-6, 6); ax.set_ylim(-2.2, 11); ax.set_zlim(0, 8.6)
    ax.set_box_aspect((12, 13.2, 8.6))
    ax.view_init(elev=16, azim=-64)
    ax.set_axis_off()
    return ax

def draw_polys(ax, polys):     # one collection so polygons are depth-sorted together
    coll = Poly3DCollection([p for p, _, _, _ in polys],
                            facecolors=[c for _, c, _, _ in polys],
                            edgecolors=[e for _, _, e, _ in polys],
                            linewidths=[w for _, _, _, w in polys])
    ax.add_collection3d(coll)

def rock_polys():
    return [(np.array(pts), matplotlib.colors.to_rgba(col, rock_alpha),
             matplotlib.colors.to_rgba("#8a9095", 0.5), 0.5)
            for pts, col in rock_faces]

import matplotlib.colors

def line3(ax, a, b, color, style="--", lw=2.6):
    ax.plot([a[0], b[0]], [a[1], b[1]], [a[2], b[2]], ls=style, lw=lw, color=color)

def wire3(ax):
    for a, b in wire:
        ax.plot([a[0], b[0]], [a[1], b[1]], [a[2], b[2]], color="#64707a", lw=1.2)

def cap(ax, s, col):
    ax.text2D(0.04, 0.05, s, transform=ax.transAxes, fontsize=9.3, color=col,
              fontweight="bold", va="bottom", ha="left",
              bbox=dict(fc="white", ec="none", alpha=0.85, pad=1.5))

# case 1
ax = add_3d(0)
polys = rock_polys()
polys.append((plane_polygon(n1a, c1a, pbox), matplotlib.colors.to_rgba(BLUE, 0.40),
              matplotlib.colors.to_rgba(BLUE, 0.9), 1.4))
polys.append((plane_polygon(n1b, c1b, pbox), matplotlib.colors.to_rgba(ORNG, 0.40),
              matplotlib.colors.to_rgba(ORNG, 0.9), 1.4))
tet = [[F1, T1, C1a], [F1, T1, C1b], [F1, C1a, C1b], [T1, C1a, C1b]]
for t in tet:
    polys.append((np.array(t), matplotlib.colors.to_rgba("#e74c3c", 0.62),
                  matplotlib.colors.to_rgba("#7b241c", 0.95), 1.6))
draw_polys(ax, polys)
wire3(ax)
ax.scatter(*F1, s=150, marker="*", color=GRN, depthshade=False, zorder=10)
cap(ax, "J1 55/145 (blue) + J2 55/215 (orange)\nred = wedge block \u00b7 star = daylighting toe " + OK, GRN)

# case 2
ax = add_3d(1)
polys = rock_polys()
polys.append((plane_polygon(n2a, c2a, pbox), matplotlib.colors.to_rgba(BLUE, 0.40),
              matplotlib.colors.to_rgba(BLUE, 0.9), 1.4))
polys.append((plane_polygon(n2b, c2b, pbox), matplotlib.colors.to_rgba(ORNG, 0.40),
              matplotlib.colors.to_rgba(ORNG, 0.9), 1.4))
draw_polys(ax, polys)
wire3(ax)
line3(ax, s2a, s2b, RED)
ax.scatter(*s2b, s=150, marker="X", color=RED, depthshade=False, zorder=10)
cap(ax, "steep joints 80/145 + 80/215\ndashed = intersection line: buried, hits the bottom " + BAD, RED)

# case 3
ax = add_3d(2)
polys = rock_polys()
polys.append((plane_polygon(n3a, c3a, pbox), matplotlib.colors.to_rgba(BLUE, 0.40),
              matplotlib.colors.to_rgba(BLUE, 0.9), 1.4))
polys.append((plane_polygon(n3b, c3b, pbox), matplotlib.colors.to_rgba(ORNG, 0.40),
              matplotlib.colors.to_rgba(ORNG, 0.9), 1.4))
draw_polys(ax, polys)
wire3(ax)
line3(ax, s3a, s3b, RED)
ax.scatter(*s3b, s=150, marker="X", color=RED, depthshade=False, zorder=10)
cap(ax, "joints dip INTO the slope (55/035 + 55/325)\ndashed = intersection plunges backward, up-hill " + BAD, RED)

# case 4
ax = add_3d(3)
polys = rock_polys()
polys.append((pa4, matplotlib.colors.to_rgba(BLUE, 0.55),
              matplotlib.colors.to_rgba(BLUE, 0.95), 2.2))
polys.append((pb4, matplotlib.colors.to_rgba(ORNG, 0.55),
              matplotlib.colors.to_rgba(ORNG, 0.95), 2.2))
draw_polys(ax, polys)
wire3(ax)
w1, w2 = line_in_box(inter_point(n1a, c1a, n1b, c1b), La,
                     ((-6, 6), (-1.2, 9.3), (0.2, 8.4)))
line3(ax, w1, w2, RED, style=(0, (2, 3)), lw=1.8)
ax.scatter(0, 2.2, 4.6, s=150, marker="X", color=RED, depthshade=False, zorder=10)
cap(ax, "small / offset planes pass by each other\nthin dashed = where the intersection line would be " + BAD, RED)

fig.text(0.5, 0.028, "A valid block needs BOTH:  (1) the two planes physically intersect inside the slope volume, and  "
         "(2) the wedge daylights into the free face \u2014\n"
         "i.e. the line of intersection plunges out of the slope (trend \u2248 face dip direction) at an angle shallower than the face dip. "
         "In RocSlope3, enlarge joint persistence and reposition if planes miss.",
         ha="center", fontsize=10.8, color="#34495e")

fig.subplots_adjust(left=0.02, right=0.98, top=0.878, bottom=0.075,
                    wspace=0.06, hspace=0.10)
import pathlib
out = pathlib.Path(__file__).with_name("rocslope3_wedge_cases.png")
fig.savefig(out, dpi=165, facecolor="white")
print("saved", out)