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