186 lines
7.8 KiB
Python
186 lines
7.8 KiB
Python
"""E03 활동면 → 절편 · 방법 묶음 평가 — 원호(중심 · 반지름) · 꺾은선(꼭짓점) · 인장균열 · 모멘트 중심.
|
|
|
|
활동 방향 +x 로 맞춘 `Prepared` 위의 좌표. 탐색(`Engine_Search`) · 결과 꾸미기(`Engine_Describe`)가 같이 씀.
|
|
"""
|
|
|
|
import numpy as np
|
|
|
|
from E03_StabSection import E03_StabSection_Engine_GLE as gle
|
|
from E03_StabSection import E03_StabSection_Engine_LEM as lem
|
|
|
|
METHODS = ("ordinary", "bishop", "janbu", "spencer", "morgenstern_price")
|
|
_SURFACE_TOL = 1e-3 # 활동면이 지표면 위로 나온 높이 허용(m)
|
|
|
|
|
|
# ---- 면 → 절편 ---------------------------------------------------------------
|
|
def _crack_start(prep, xE, xX, base_fn):
|
|
"""인장균열 — 머리 쪽에서 활동면 깊이가 균열 깊이에 닿는 x (없으면 원래 진입점)."""
|
|
if prep.crack_depth <= 0:
|
|
return xE, np.zeros(len(xE))
|
|
t = np.linspace(0.0, 1.0, 201)
|
|
xs = xE[:, None] + (xX - xE)[:, None] * t
|
|
depth = prep.surface(xs) - base_fn(xs)
|
|
reach = depth >= prep.crack_depth
|
|
k = np.where(reach.any(1), np.argmax(reach, 1), 0)
|
|
k0 = np.maximum(k - 1, 0)
|
|
rows = np.arange(len(xE))
|
|
d0, d1 = depth[rows, k0], depth[rows, k]
|
|
w = np.where(d1 > d0, (prep.crack_depth - d0) / np.maximum(d1 - d0, 1e-12), 0.0)
|
|
x_new = np.where(k > 0, xs[rows, k0] + np.clip(w, 0, 1) * (xs[rows, k] - xs[rows, k0]), xE)
|
|
return x_new, np.where(k > 0, prep.crack_depth, 0.0)
|
|
|
|
|
|
def _finish_slices(prep, xe, ze, crack, min_depth):
|
|
sl = lem.slices(prep, xe, ze)
|
|
if prep.crack_depth > 0:
|
|
lem.add_crack_water(sl, crack, prep.crack_water)
|
|
height = sl["zt"] - sl["zb"]
|
|
sl["valid"] &= np.all(height > -_SURFACE_TOL, axis=1)
|
|
sl["valid"] &= np.nanmax(height, axis=1) >= max(min_depth, 1e-3)
|
|
return sl
|
|
|
|
|
|
def circle_slices(prep, xc, zc, R, xE, xX, n, min_depth=0.0):
|
|
xc, zc, R = (np.asarray(v, float) for v in (xc, zc, R))
|
|
|
|
def base(x):
|
|
return zc[:, None] - np.sqrt(np.clip(R[:, None] ** 2 - (x - xc[:, None]) ** 2, 0, None))
|
|
|
|
xE, crack = _crack_start(prep, np.asarray(xE, float), np.asarray(xX, float), base)
|
|
t = np.linspace(0.0, 1.0, n + 1)
|
|
xe = xE[:, None] + (xX - xE)[:, None] * t
|
|
sl = _finish_slices(prep, xe, base(xe), crack, min_depth)
|
|
sl["valid"] &= (xE >= xc - R) & (xX <= xc + R)
|
|
return sl
|
|
|
|
|
|
def poly_base(vx, vz, x):
|
|
"""꼭짓점(S, V) 꺾은선을 x(S, m) 에서 — 행마다 선형 보간(양 끝 밖은 끝 마디 연장)."""
|
|
z = None
|
|
for k in range(vx.shape[1] - 1):
|
|
x0, x1 = vx[:, k : k + 1], vx[:, k + 1 : k + 2]
|
|
t = (x - x0) / np.where(x1 - x0 > 1e-12, x1 - x0, 1e-12)
|
|
seg = vz[:, k : k + 1] + t * (vz[:, k + 1 : k + 2] - vz[:, k : k + 1])
|
|
z = seg if z is None else np.where(x >= x0, seg, z)
|
|
return z
|
|
|
|
|
|
def poly_slices(prep, vx, vz, n, min_depth=0.0):
|
|
vx, vz = np.asarray(vx, float), np.asarray(vz, float)
|
|
|
|
def base(x):
|
|
return poly_base(vx, vz, x)
|
|
|
|
xE, crack = _crack_start(prep, vx[:, 0], vx[:, -1], base)
|
|
# 절편 경계 = 고른 간격 + 안쪽 꼭짓점(꺾임이 절편 안에 묻히지 않게 — SLOPE/W 와 같은 생각)
|
|
t = np.linspace(0.0, 1.0, n + 1 - (vx.shape[1] - 2))
|
|
xe = xE[:, None] + (vx[:, -1] - xE)[:, None] * t
|
|
inner = np.clip(vx[:, 1:-1], xE[:, None], vx[:, -1:])
|
|
xe = np.sort(np.concatenate([xe, inner], axis=1), axis=1)
|
|
sl = _finish_slices(prep, xe, base(xe), crack, min_depth)
|
|
sl["valid"] &= np.all(np.diff(vx, axis=1) > 1e-6, axis=1)
|
|
return sl
|
|
|
|
|
|
def three_point_center(vx, vz):
|
|
"""비원호 모멘트 중심 = 진입 · 가장 깊은 꼭짓점 · 진출을 지나는 원의 중심(곧은 면은 위로 멀리)."""
|
|
p1 = np.column_stack([vx[:, 0], vz[:, 0]])
|
|
p3 = np.column_stack([vx[:, -1], vz[:, -1]])
|
|
chord_z = vz[:, :1] + (vx - vx[:, :1]) * (
|
|
(vz[:, -1:] - vz[:, :1]) / np.maximum(vx[:, -1:] - vx[:, :1], 1e-9)
|
|
)
|
|
k = np.argmax(chord_z - vz, axis=1)
|
|
p2 = np.column_stack([vx[np.arange(len(vx)), k], vz[np.arange(len(vx)), k]])
|
|
ax, ay = p1.T
|
|
bx, by = p2.T
|
|
cx, cy = p3.T
|
|
d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
|
|
ux = (
|
|
(ax**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)
|
|
) / np.where(np.abs(d) < 1e-9, 1, d)
|
|
uy = (
|
|
(ax**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax - cx) + (cx**2 + cy**2) * (bx - ax)
|
|
) / np.where(np.abs(d) < 1e-9, 1, d)
|
|
far = np.abs(d) < 1e-9
|
|
span = np.hypot(cx - ax, cy - ay)
|
|
ux = np.where(far, 0.5 * (ax + cx), ux)
|
|
uy = np.where(far, np.maximum(ay, cy) + 10 * span, uy)
|
|
return np.column_stack([ux, uy])
|
|
|
|
|
|
# ---- 방법 묶음 평가 ------------------------------------------------------------
|
|
def evaluate(sl, center, methods=METHODS, scan=True, mp_seed=None) -> dict:
|
|
"""방법별 F 묶음 — 출발값 차례 Fellenius → Bishop → (Janbu · Spencer · M-P).
|
|
|
|
M-P = 임계면 몇 개에만(차례 계산이라 느림 · 탐색 고리에서 쓰지 않음) · 출발 참고값 = mp_seed 또는 Bishop.
|
|
"""
|
|
out = {}
|
|
base = lem.ordinary(sl, center)
|
|
seed = np.where(np.isfinite(base["F"]), base["F"], 1.5)
|
|
if "ordinary" in methods:
|
|
out["ordinary"] = base
|
|
if {"bishop", "spencer", "morgenstern_price"} & set(methods):
|
|
out["bishop"] = lem.bishop(sl, center, seed)
|
|
seed = np.where(np.isfinite(out["bishop"]["F"]), out["bishop"]["F"], seed)
|
|
if "janbu" in methods:
|
|
out["janbu"] = lem.janbu(sl, seed)
|
|
if "spencer" in methods:
|
|
out["spencer"] = lem.spencer(sl, center, seed, scan)
|
|
if "morgenstern_price" in methods:
|
|
ref = seed if mp_seed is None else np.where(np.isfinite(mp_seed), mp_seed, seed)
|
|
out["morgenstern_price"] = gle.morgenstern_price(sl, center, F0=ref, curve=False)
|
|
return {m: out[m] for m in methods}
|
|
|
|
|
|
def objective(method):
|
|
def fun(sl, center):
|
|
return evaluate(sl, center, (method,), scan=False)[method]["F"]
|
|
|
|
return fun
|
|
|
|
|
|
# ---- 원호 기하 ---------------------------------------------------------------
|
|
def circle_from_chord(prep, xE, xX, theta):
|
|
"""진입 · 진출 x 와 중심각 θ → 중심 · 반지름(중심은 현 위쪽 · 두 끝이 중심보다 아래여야)."""
|
|
zE, zX = prep.surface(xE), prep.surface(xX)
|
|
dx, dz = xX - xE, zX - zE
|
|
c = np.hypot(dx, dz)
|
|
R = c / (2 * np.sin(theta / 2))
|
|
d = R * np.cos(theta / 2)
|
|
nx, nz = -dz / np.maximum(c, 1e-12), dx / np.maximum(c, 1e-12)
|
|
flip = nz < 0
|
|
nx, nz = np.where(flip, -nx, nx), np.where(flip, -nz, nz)
|
|
xc = 0.5 * (xE + xX) + d * nx
|
|
zc = 0.5 * (zE + zX) + d * nz
|
|
ok = (zE <= zc) & (zX <= zc) & (dx > 0) & np.isfinite(zE) & np.isfinite(zX)
|
|
return xc, zc, R, ok
|
|
|
|
|
|
def circle_cross_surface(prep, xc, zc, R):
|
|
"""원 아래 반쪽과 지표면의 교점 중 지표면 아래로 들어간 가장 넓은 구간 → (진입, 진출)."""
|
|
sx, sz = prep.surface_x, prep.surface_z
|
|
ax, az, bx, bz = sx[:-1], sz[:-1], sx[1:], sz[1:]
|
|
dx, dz = bx - ax, bz - az
|
|
fx, fz = ax[None, :] - xc[:, None], az[None, :] - zc[:, None]
|
|
a = dx**2 + dz**2
|
|
b = 2 * (fx * dx + fz * dz)
|
|
c = fx**2 + fz**2 - R[:, None] ** 2
|
|
disc = b**2 - 4 * a * c
|
|
root = np.sqrt(np.clip(disc, 0, None))
|
|
roots = []
|
|
for sign in (-1, 1):
|
|
t = (-b + sign * root) / (2 * np.maximum(a, 1e-12))
|
|
x = ax + t * dx
|
|
z = az + t * dz
|
|
keep = (disc >= 0) & (t >= 0) & (t <= 1) & (z <= zc[:, None] + 1e-9)
|
|
roots.append(np.where(keep, x, np.nan))
|
|
xs = np.sort(np.concatenate(roots, axis=1), axis=1)
|
|
mids = 0.5 * (xs[:, :-1] + xs[:, 1:])
|
|
arc = zc[:, None] - np.sqrt(np.clip(R[:, None] ** 2 - (mids - xc[:, None]) ** 2, 0, None))
|
|
below = np.isfinite(mids) & (prep.surface(np.nan_to_num(mids)) > arc)
|
|
width = np.where(below, xs[:, 1:] - xs[:, :-1], -1.0)
|
|
k = np.argmax(width, axis=1)
|
|
rows = np.arange(len(xc))
|
|
ok = width[rows, k] > 0
|
|
return xs[rows, k], xs[rows, k + 1], ok
|