eomsangdonandClaude Opus 5.5
f25dac9000
- 새 Engine_Support: 첫 교점 · min(인장 · 인발 · 벗겨짐) / 간격 · 앵커 정착장만 부착 · 프리스트레스 앵커 = 능동(min(프리스트레스, 한도)) · 네일 = 수동(힘 / F) - LEM · GLE: 능동은 연직 · 수평 하중 칸 · 수동은 c' 칸(1/F) · 전체 평형 수동 항은 분자로(같은 근) · Fellenius 수동은 고정점 - Search: 보강 있으면 진입 한계만 보강재 끝 뒤까지(진출 = 원래 한계) · 판정 임계면 support_forces - Station: analyze_station(supports=None) · side 칸 · 없으면 머리 x 쪽 비탈 · supports 없으면 결과 불변 - 시험 test_536_3_stab_supports 10(쐐기 닫힌 꼴 · 힘 꼴 손계산 · Slide2 #47 · #85 · #86 · 탐색 넓힘 · 불변 · 측점 꼴) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
642 lines
28 KiB
Python
642 lines
28 KiB
Python
"""E03 임계 활동면 탐색 — 원호(Grid & Radius · Entry & Exit · 다듬기) · 비원호(블록 + 꼭짓점 최적화).
|
||
|
||
활동 방향 +x 로 맞춘 `Prepared` 위에서 찾고 좌표를 원래 방향으로 되돌려 냄.
|
||
다듬기 · 비원호 = 차분진화(best/1/bin · 무리 전체를 한 번에 평가 · seed 고정 → 같은 입력 같은 답)
|
||
+ 비원호 마지막 국부 다듬기 = 묶음 Nelder-Mead(후보 넷을 한 번에). SciPy 최적화는 점마다 파이썬 고리라 느려
|
||
직접 묶음 꼴로 씀. 방법마다 찾은 면을 모든 방법으로 다시 셈(교차) → 한 방법 탐색이 놓친 골짜기를 메움.
|
||
상용 SW 대응: Grid & Radius(SLOPE/W · Slide) · Entry & Exit(SLOPE/W) · Auto Refine(Slide) ·
|
||
Block Search + Optimize(Slide · SLOPE/W 꼭짓점 최적화).
|
||
"""
|
||
|
||
import time
|
||
|
||
import numpy as np
|
||
|
||
from E03_StabSection import E03_StabSection_Engine_LEM as lem
|
||
from E03_StabSection import E03_StabSection_Engine_Support as support
|
||
from E03_StabSection.E03_StabSection_Engine_Describe import describe, safety_map
|
||
from E03_StabSection.E03_StabSection_Engine_Surface import ( # noqa: F401 — 밖에서 Search 로 부르던 이름
|
||
METHODS,
|
||
circle_cross_surface,
|
||
circle_from_chord,
|
||
circle_slices,
|
||
evaluate,
|
||
objective,
|
||
poly_base,
|
||
poly_slices,
|
||
three_point_center,
|
||
)
|
||
|
||
_BIG = 1e3
|
||
# 탐색 예산(실험 손잡이)
|
||
LS_ITER = 30
|
||
LS_ITER_SPENCER = 60 # Spencer 비원호는 풀리는 면이 조각나 더 걸음
|
||
POLISH = "nm"
|
||
NM_ITER = 120
|
||
PATTERN_ITER = 30
|
||
DE_MAXITER = 60
|
||
DE_PATIENCE = 12
|
||
SPENCER_SEED_POOL = 300
|
||
#: 교차 재평가에 넣는 방법별 마지막 무리 상위 원 수
|
||
CROSS_TOP = 10
|
||
#: Spencer 원호 씨앗 — Bishop 임계 둘레 구름 크기
|
||
CLOUD_SIZE = 200
|
||
#: 엄밀해 「해 없음」 허용 — Bishop 임계 원에서 해 없어도 찾은 최소가 Bishop × (1 + 이 값) 안이면 값을 살림
|
||
FAR_TOL = 0.02
|
||
|
||
|
||
# ---- 원호 ---------------------------------------------------------------------
|
||
def _exit_lo(limits) -> float:
|
||
"""진출 아래 한계 — limits 셋째 칸(보강으로 진입만 넓힌 때 원래 한계) · 없으면 limits[0]."""
|
||
return limits[2] if len(limits) > 2 else limits[0]
|
||
|
||
|
||
def entry_exit_circles(prep, limits, n_points=25, thetas_deg=(10, 20, 35, 50, 65, 80, 100, 120)):
|
||
xl, xr = limits[:2]
|
||
pts = np.linspace(xl, xr, n_points)
|
||
i, j = np.triu_indices(n_points, k=1)
|
||
xE, xX = pts[i], pts[j]
|
||
xE = np.repeat(xE, len(thetas_deg))
|
||
xX = np.repeat(xX, len(thetas_deg))
|
||
th = np.tile(np.radians(thetas_deg), len(i))
|
||
xc, zc, R, ok = circle_from_chord(prep, xE, xX, th)
|
||
ok &= xX >= _exit_lo(limits) - 1e-6
|
||
return xc[ok], zc[ok], R[ok], xE[ok], xX[ok]
|
||
|
||
|
||
def grid_centers(prep, limits, nx=15, nz=15):
|
||
xl, xr = limits[:2]
|
||
z_in = prep.surface(np.linspace(xl, xr, 50))
|
||
H = max(np.nanmax(z_in) - np.nanmin(z_in), 1.0)
|
||
L = max(xr - xl, 1.0)
|
||
gx = np.linspace(xl, xr + 0.25 * L, nx)
|
||
gz = np.linspace(np.nanmax(z_in) + 0.1 * H, np.nanmax(z_in) + max(1.5 * H, 0.6 * L), nz)
|
||
return gx, gz
|
||
|
||
|
||
def grid_radius_circles(prep, limits, gx, gz, n_radii=10):
|
||
cx, cz = np.meshgrid(gx, gz)
|
||
cx, cz = cx.ravel(), cz.ravel()
|
||
sx = np.linspace(limits[0], limits[1], 80)
|
||
sz = prep.surface(sx)
|
||
dist = np.sqrt((sx[None, :] - cx[:, None]) ** 2 + (sz[None, :] - cz[:, None]) ** 2)
|
||
r0 = np.nanmin(dist, axis=1) + 0.05
|
||
r1 = np.maximum(cz - prep.bottom_z, r0 + 0.1)
|
||
frac = np.linspace(0.0, 1.0, n_radii + 1)[1:]
|
||
R = (r0[:, None] + (r1 - r0)[:, None] * frac).ravel()
|
||
owner = np.repeat(np.arange(len(cx)), n_radii)
|
||
xc, zc = cx[owner], cz[owner]
|
||
xE, xX, ok = circle_cross_surface(prep, xc, zc, R)
|
||
ok &= (xE >= limits[0] - 1e-6) & (xX <= limits[1] + 1e-6) & (xX >= _exit_lo(limits) - 1e-6)
|
||
return xc[ok], zc[ok], R[ok], xE[ok], xX[ok], owner[ok]
|
||
|
||
|
||
def _circle_eval(prep, xc, zc, R, xE, xX, n, min_depth, methods, chunk=4000):
|
||
res = {m: np.full(len(xc), np.nan) for m in methods}
|
||
for s in range(0, len(xc), chunk):
|
||
e = slice(s, s + chunk)
|
||
sl = circle_slices(prep, xc[e], zc[e], R[e], xE[e], xX[e], n, min_depth)
|
||
center = np.column_stack([xc[e], zc[e]])
|
||
out = evaluate(sl, center, methods)
|
||
for m in methods:
|
||
res[m][e] = out[m]["F"]
|
||
return res
|
||
|
||
|
||
def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=None):
|
||
"""원호 다듬기 — (진입 x, 진출 x, 중심각) 차분진화 · 첫 무리 = 거친 탐색의 좋은 면."""
|
||
fun = objective(method)
|
||
lo = np.array([limits[0], _exit_lo(limits), np.radians(3)])
|
||
hi = np.array([limits[1], limits[1], np.radians(170)])
|
||
|
||
def obj(p):
|
||
xc, zc, R, ok = circle_from_chord(prep, p[:, 0], p[:, 1], p[:, 2])
|
||
F = np.full(len(p), _BIG)
|
||
if ok.any():
|
||
sl = circle_slices(prep, xc[ok], zc[ok], R[ok], p[ok, 0], p[ok, 1], n, min_depth)
|
||
f = fun(sl, np.column_stack([xc[ok], zc[ok]]))
|
||
F[ok] = np.where(np.isfinite(f), f, _BIG)
|
||
return F
|
||
|
||
init = _population(seeds, lo, hi, 30, seed)
|
||
x, fmin, pop, f = _de(obj, lo, hi, init, maxiter or DE_MAXITER, seed, patience=DE_PATIENCE)
|
||
top = pop[np.argsort(f)[:CROSS_TOP]]
|
||
top = top[obj(top) < _BIG]
|
||
xc, zc, R, _ = circle_from_chord(prep, top[:, 0], top[:, 1], top[:, 2])
|
||
near = [
|
||
(float(a), float(b), float(c), float(d), float(e))
|
||
for a, b, c, d, e in zip(xc, zc, R, top[:, 0], top[:, 1])
|
||
]
|
||
xc, zc, R, _ = circle_from_chord(prep, x[0:1], x[1:2], x[2:3])
|
||
best = (float(xc[0]), float(zc[0]), float(R[0]), float(x[0]), float(x[1]))
|
||
return float(fmin), best, near, top
|
||
|
||
|
||
def _cloud(prep, limits, method, p0, n, min_depth, seed, k=CLOUD_SIZE, keep=10):
|
||
"""원 (진입 x, 진출 x, 중심각) 둘레 구름 k 개를 method 로 한 번에 셈 → 풀리는 좋은 것 keep 개."""
|
||
rng = np.random.default_rng(seed)
|
||
span = limits[1] - limits[0]
|
||
lo = np.array([limits[0], _exit_lo(limits), np.radians(3)])
|
||
hi = np.array([limits[1], limits[1], np.radians(170)])
|
||
p = np.clip(p0 + rng.normal(size=(k, 3)) * [0.03 * span, 0.03 * span, np.radians(4)], lo, hi)
|
||
xc, zc, R, ok = circle_from_chord(prep, p[:, 0], p[:, 1], p[:, 2])
|
||
F = np.full(k, np.inf)
|
||
if ok.any():
|
||
sl = circle_slices(prep, xc[ok], zc[ok], R[ok], p[ok, 0], p[ok, 1], n, min_depth)
|
||
f = evaluate(sl, np.column_stack([xc[ok], zc[ok]]), (method,), scan=False)[method]["F"]
|
||
F[ok] = np.where(np.isfinite(f), f, np.inf)
|
||
order = np.argsort(F)[:keep]
|
||
return p[order[np.isfinite(F[order])]].reshape(-1, 3)
|
||
|
||
|
||
def _de(obj, lo, hi, init, maxiter, seed, cr=0.7, tol=1e-6, patience=25):
|
||
"""묶음 차분진화(best/1/bin · 축척 0.5~1.0 흔들기 = SciPy 기본 전략) — 최소 x · 값 · 마지막 무리 · 값."""
|
||
rng = np.random.default_rng(seed)
|
||
pop = np.clip(np.asarray(init, float), lo, hi)
|
||
f = obj(pop)
|
||
P, D = pop.shape
|
||
best_hist = f.min()
|
||
still = 0
|
||
for _ in range(maxiter):
|
||
b = pop[np.argmin(f)]
|
||
r1 = rng.integers(0, P, P)
|
||
r2 = (r1 + rng.integers(1, P, P)) % P
|
||
trial = b + rng.uniform(0.5, 1.0) * (pop[r1] - pop[r2])
|
||
cross = rng.random((P, D)) < cr
|
||
cross[np.arange(P), rng.integers(0, D, P)] = True
|
||
trial = np.clip(np.where(cross, trial, pop), lo, hi)
|
||
ft = obj(trial)
|
||
better = ft <= f
|
||
pop[better], f[better] = trial[better], ft[better]
|
||
still = still + 1 if best_hist - f.min() < tol else 0
|
||
best_hist = min(best_hist, f.min())
|
||
if still >= patience:
|
||
break
|
||
k = int(np.argmin(f))
|
||
return pop[k].copy(), float(f[k]), pop, f
|
||
|
||
|
||
def _local_search(obj, starts, lo, hi, step, maxiter, seed, lam=24, starts_keep=8):
|
||
"""여러 출발점 (1+λ) 국부 탐색 — 출발점마다 λ 개를 흔들어 나아지면 옮기고 걸음 ×1.5 · 아니면 ×0.6.
|
||
|
||
상용 SW 의 꼭짓점 이동 최적화(SLOPE/W Optimization · Slide Surface Altering)와 같은 생각.
|
||
모든 출발점 · 후보를 한 번에 평가. 돌려주는 값 = 최소 x · 값 · 출발점별 마지막 x · 값.
|
||
"""
|
||
rng = np.random.default_rng(seed)
|
||
x = np.clip(starts[:starts_keep], lo, hi)
|
||
f = obj(x)
|
||
sig = np.tile(step, (len(x), 1))
|
||
for _ in range(maxiter):
|
||
cand = x[:, None, :] + rng.normal(size=(len(x), lam, x.shape[1])) * sig[:, None, :]
|
||
cand = np.clip(cand, lo, hi).reshape(-1, x.shape[1])
|
||
fc = obj(cand).reshape(len(x), lam)
|
||
k = np.argmin(fc, axis=1)
|
||
best = fc[np.arange(len(x)), k]
|
||
up = best < f - 1e-9
|
||
x[up] = cand.reshape(len(x), lam, -1)[np.arange(len(x)), k][up]
|
||
f[up] = best[up]
|
||
sig *= np.where(up, 1.5, 0.6)[:, None]
|
||
if np.all(sig < step * 1e-3):
|
||
break
|
||
i = int(np.argmin(f))
|
||
return x[i].copy(), float(f[i]), x, f
|
||
|
||
|
||
def _nelder_mead(obj, x0, lo, hi, maxiter, xatol=1e-4, fatol=1e-6):
|
||
"""묶음 Nelder-Mead — 반사 · 확장 · 바깥 · 안쪽 수축 넷을 한 번에 셈(판단은 표준 NM 그대로 · 줄어듦도 한 번에).
|
||
|
||
SciPy NM 은 점마다 따로 불러 한 면 계산의 파이썬 부담이 큼 — 부르는 수를 줄임. 첫 꼭짓점 = SciPy 와 같은 5 % 흔들기.
|
||
"""
|
||
D = len(x0)
|
||
sim = np.tile(np.asarray(x0, float), (D + 1, 1))
|
||
for k in range(D):
|
||
sim[k + 1, k] = sim[k + 1, k] * 1.05 if sim[k + 1, k] != 0 else 0.00025
|
||
sim = np.clip(sim, lo, hi)
|
||
fs = obj(sim)
|
||
for _ in range(maxiter):
|
||
order = np.argsort(fs, kind="stable")
|
||
sim, fs = sim[order], fs[order]
|
||
if np.max(np.abs(sim[1:] - sim[0])) <= xatol and np.max(np.abs(fs[1:] - fs[0])) <= fatol:
|
||
break
|
||
xbar = sim[:-1].mean(0)
|
||
xw = sim[-1]
|
||
cand = np.clip(
|
||
np.vstack(
|
||
[2 * xbar - xw, 3 * xbar - 2 * xw, 1.5 * xbar - 0.5 * xw, 0.5 * xbar + 0.5 * xw]
|
||
),
|
||
lo,
|
||
hi,
|
||
)
|
||
fr, fe, fc, fcc = obj(cand)
|
||
if fr < fs[0]:
|
||
sim[-1], fs[-1] = (cand[1], fe) if fe < fr else (cand[0], fr)
|
||
elif fr < fs[-2]:
|
||
sim[-1], fs[-1] = cand[0], fr
|
||
elif fr < fs[-1] and fc <= fr:
|
||
sim[-1], fs[-1] = cand[2], fc
|
||
elif fr >= fs[-1] and fcc < fs[-1]:
|
||
sim[-1], fs[-1] = cand[3], fcc
|
||
else:
|
||
sim[1:] = np.clip(sim[0] + 0.5 * (sim[1:] - sim[0]), lo, hi)
|
||
fs[1:] = obj(sim[1:])
|
||
k = int(np.argmin(fs))
|
||
return sim[k].copy(), float(fs[k])
|
||
|
||
|
||
def _pattern(obj, x, fx, lo, hi, step, iters):
|
||
"""묶음 좌표 패턴 탐색 — 축마다 ±걸음 2D 점을 한 번에 셈 · 나아지면 옮기고 그 쪽으로 한 번 더 · 아니면 걸음 ½."""
|
||
h = step.copy()
|
||
for _ in range(iters):
|
||
cand = np.clip(np.vstack([x + np.diag(h), x - np.diag(h)]), lo, hi)
|
||
f = obj(cand)
|
||
k = int(np.argmin(f))
|
||
if f[k] < fx - 1e-9:
|
||
move = cand[k] - x
|
||
x2 = np.clip(cand[k] + move, lo, hi)
|
||
f2 = float(obj(x2[None, :])[0])
|
||
x, fx = (x2, f2) if f2 < f[k] else (cand[k], float(f[k]))
|
||
else:
|
||
h *= 0.5
|
||
if np.all(h < step * 1e-3):
|
||
break
|
||
return x, fx
|
||
|
||
|
||
def _population(seeds, lo, hi, size, seed, spread=0.0):
|
||
"""첫 무리 = 씨앗 + (spread > 0 이면 씨앗 둘레 흔든 것 80 % · 고른 무작위 20 %)."""
|
||
rng = np.random.default_rng(seed)
|
||
seeds = np.clip(np.asarray(seeds, float).reshape(-1, len(lo)), lo, hi)
|
||
extra = max(size - len(seeds), 5)
|
||
rand = lo + rng.random((extra, len(lo))) * (hi - lo)
|
||
if spread > 0 and len(seeds):
|
||
k = int(0.8 * extra)
|
||
base = seeds[rng.integers(0, len(seeds), k)]
|
||
rand[:k] = np.clip(base + rng.normal(0, spread, (k, len(lo))) * (hi - lo), lo, hi)
|
||
return np.vstack([seeds, rand])
|
||
|
||
|
||
# ---- 비원호(블록 + 꼭짓점 최적화) -------------------------------------------------
|
||
def _decode(prep, p, m):
|
||
"""변수 [진입 x, 진출 x, 위치 s_1..s_m, 깊이 d_1..d_m] → 꼭짓점(S, m+2)."""
|
||
xE, xX = p[:, 0], p[:, 1]
|
||
s = np.sort(p[:, 2 : 2 + m], axis=1)
|
||
d = p[:, 2 + m : 2 + 2 * m]
|
||
ix = xE[:, None] + s * (xX - xE)[:, None]
|
||
vx = np.column_stack([xE, ix, xX])
|
||
vz = np.column_stack([prep.surface(xE), prep.surface(ix) - d, prep.surface(xX)])
|
||
return vx, vz
|
||
|
||
|
||
def _convex(vx, vz, tol_deg=2.0):
|
||
ang = np.degrees(np.arctan2(vz[:, :-1] - vz[:, 1:], np.maximum(np.diff(vx, axis=1), 1e-9)))
|
||
return np.all(np.diff(ang, axis=1) <= tol_deg, axis=1)
|
||
|
||
|
||
def layer_bottom_depth(prep, x, material):
|
||
"""x 마다 그 지층(번호) 바닥 바로 위까지 깊이 — 없으면 nan(블록 씨앗: 약층 바닥을 따라 미끄러짐)."""
|
||
x = np.atleast_1d(x)
|
||
bot, top, mat = prep.columns(x[None, :])
|
||
zt = prep.surface(x)[None, :]
|
||
z = np.where((mat == material) & (top > bot), bot, np.nan)
|
||
with np.errstate(all="ignore"):
|
||
zb = np.nanmax(z, axis=-1)
|
||
return (zt - zb - 0.02)[0]
|
||
|
||
|
||
def optimize_polyline(
|
||
prep, limits, method, seeds, n, min_depth, m=4, seed=0, maxiter=None, polish=True
|
||
):
|
||
"""비원호 꼭짓점 최적화 — 차분진화 뒤 가장 좋은 면을 Nelder-Mead 로 국부 다듬기(상용 SAO · Optimize 대응)."""
|
||
fun = objective(method)
|
||
xl, xr = limits[:2]
|
||
depth_max = max(prep.surface(np.linspace(xl, xr, 50)).max() - prep.bottom_z, 1.0)
|
||
lo = np.array([xl, _exit_lo(limits)] + [0.02] * m + [0.0] * m)
|
||
hi = np.array([xr, xr] + [0.98] * m + [depth_max] * m)
|
||
|
||
def obj(p):
|
||
p = np.atleast_2d(p)
|
||
vx, vz = _decode(prep, p, m)
|
||
ok = (p[:, 1] - p[:, 0] > 0.2) & _convex(vx, vz)
|
||
F = np.full(len(p), _BIG)
|
||
if ok.any():
|
||
sl = poly_slices(prep, vx[ok], vz[ok], n, min_depth)
|
||
f = fun(sl, three_point_center(vx[ok], vz[ok]))
|
||
F[ok] = np.where(np.isfinite(f), f, _BIG)
|
||
return F
|
||
|
||
step = np.array([0.5, 0.5] + [0.03] * m + [0.3] * m)
|
||
x, fmin, pop, f = _local_search(
|
||
obj, np.asarray(seeds, float), lo, hi, step, maxiter or LS_ITER, seed
|
||
)
|
||
if polish and fmin < _BIG and POLISH == "pattern":
|
||
x, fmin = _pattern(obj, x, fmin, lo, hi, step * 0.1, PATTERN_ITER)
|
||
elif polish and fmin < _BIG and POLISH == "nm":
|
||
xn, fn = _nelder_mead(obj, x, lo, hi, NM_ITER)
|
||
if fn < fmin:
|
||
x, fmin = xn, fn
|
||
order = np.argsort(f)
|
||
vx, vz = _decode(prep, x[None, :], m)
|
||
return fmin, x, (vx[0], vz[0]), [x, *pop[order]]
|
||
|
||
|
||
def _best_seeds(prep, seeds, m, n, min_depth, keep):
|
||
"""씨앗을 Janbu 간편으로 한 번 셈 → 좋은 것 keep 개(최적화 첫 무리 · 나쁜 씨앗이 무리를 흐리지 않게)."""
|
||
p = np.asarray(seeds, float)
|
||
vx, vz = _decode(prep, p, m)
|
||
ok = (p[:, 1] - p[:, 0] > 0.2) & _convex(vx, vz)
|
||
F = np.full(len(p), np.inf)
|
||
if ok.any():
|
||
sl = poly_slices(prep, vx[ok], vz[ok], n, min_depth)
|
||
f = lem.janbu(sl)["F"]
|
||
F[ok] = np.where(np.isfinite(f), f, np.inf)
|
||
order = np.argsort(F)[:keep]
|
||
return [seeds[i] for i in order if np.isfinite(F[i])] or seeds[:keep]
|
||
|
||
|
||
def circle_to_params(prep, circle, m):
|
||
xc, zc, R, xE, xX = circle
|
||
s = (np.arange(1, m + 1) / (m + 1)).tolist()
|
||
ix = xE + np.array(s) * (xX - xE)
|
||
arc = zc - np.sqrt(np.clip(R**2 - (ix - xc) ** 2, 0, None))
|
||
d = np.clip(prep.surface(ix) - arc, 0, None)
|
||
return [xE, xX, *s, *d.tolist()]
|
||
|
||
|
||
def _ray_to_surface(prep, x0, z0, angle_deg, sign):
|
||
"""(x0, z0) 에서 수평과 angle 로 위로(sign = −1 왼쪽 · +1 오른쪽) 쏜 반직선이 지표와 만나는 x."""
|
||
t = np.linspace(0.0, 1.0, 400)[None, :] * 200.0
|
||
dx = sign * np.cos(np.radians(angle_deg))
|
||
dz = np.sin(np.radians(angle_deg))
|
||
xs = x0[:, None] + dx * t
|
||
above = (z0[:, None] + dz * t) >= prep.surface(xs)
|
||
k = np.argmax(above, axis=1)
|
||
hit = above[np.arange(len(x0)), k]
|
||
return np.where(hit & (k > 0), xs[np.arange(len(x0)), k], np.nan)
|
||
|
||
|
||
def block_seeds(prep, limits, m, n_points=12, angles=((55, 35), (65, 25))):
|
||
"""블록 씨앗(Slide Block Search 대응) — 지층마다 바닥선 위 두 점을 잇는 토막 + 머리 쪽 능동 쐐기 ·
|
||
끝 쪽 수동 쐐기를 지표까지 올림. 안쪽 꼭짓점 = 토막 두 끝 + 사이 (m − 2) 점."""
|
||
xl, xr = limits[:2]
|
||
xs = np.linspace(xl, xr, n_points)
|
||
i, j = np.triu_indices(n_points, k=1)
|
||
xa, xb = xs[i], xs[j]
|
||
out = []
|
||
for mat in range(len(prep.materials)):
|
||
da, db = layer_bottom_depth(prep, xa, mat), layer_bottom_depth(prep, xb, mat)
|
||
ok = np.isfinite(da) & np.isfinite(db)
|
||
if not ok.any():
|
||
continue
|
||
a, b = xa[ok], xb[ok]
|
||
za, zb = prep.surface(a) - da[ok], prep.surface(b) - db[ok]
|
||
for act, pas in angles:
|
||
xE = _ray_to_surface(prep, a, za, act, -1.0)
|
||
xX = _ray_to_surface(prep, b, zb, pas, 1.0)
|
||
good = np.isfinite(xE) & np.isfinite(xX) & (xE >= xl) & (xX <= xr)
|
||
good &= xX >= _exit_lo(limits)
|
||
for e, x0, x1, xx in zip(xE[good], a[good], b[good], xX[good]):
|
||
ix = np.linspace(x0, x1, m)
|
||
d = layer_bottom_depth(prep, ix, mat)
|
||
if np.isnan(d).any():
|
||
continue
|
||
span = xx - e
|
||
out.append([e, xx, *((ix - e) / span).tolist(), *np.clip(d, 0.0, None).tolist()])
|
||
return out
|
||
|
||
|
||
# ---- 한 비탈 해석 ----------------------------------------------------------------
|
||
def analyze(
|
||
model: dict,
|
||
*,
|
||
direction: int = 1,
|
||
limits=None,
|
||
methods=METHODS,
|
||
circular=True,
|
||
noncircular=True,
|
||
n_slices=50,
|
||
coarse_slices=30,
|
||
min_depth=0.0,
|
||
grid=(15, 15, 10),
|
||
entry_exit_points=25,
|
||
polyline_vertices=4,
|
||
seed=0,
|
||
) -> dict:
|
||
"""한 비탈(활동 방향 하나)의 방법별 최소 안전율 · 임계면 · 절편표 · 등안전율 격자 · 수렴 정보."""
|
||
t0 = time.perf_counter()
|
||
prep = lem.Prepared(model, direction)
|
||
s = prep.direction
|
||
if limits is None:
|
||
lim = (float(prep.surface_x[0]), float(prep.surface_x[-1]))
|
||
else:
|
||
lim = tuple(sorted((s * float(limits[0]), s * float(limits[1]))))
|
||
if prep.supports is not None:
|
||
# 보강 — 임계면이 보강재 끝 뒤(더 깊이 · 더 멀리)로 옮겨 가는 것도 찾게 진입 한계만 끝 뒤까지 넓힘
|
||
# (진출 = 원래 한계 셋째 칸 — 설계 비탈 밖 자연 비탈만 무너지는 면은 보강 전처럼 안 봄)
|
||
reach = float(np.min(support.tips(prep.supports))) - max(
|
||
2.0, 0.25 * prep.supports["L"].max()
|
||
)
|
||
lim = (max(min(lim[0], reach), float(prep.surface_x[0])), lim[1], lim[0])
|
||
results: dict = {m: {} for m in methods}
|
||
counts = {}
|
||
contour = smap = None
|
||
best_circles: dict = {}
|
||
if circular:
|
||
a = entry_exit_circles(prep, lim, entry_exit_points)
|
||
gx, gz = grid_centers(prep, lim, grid[0], grid[1])
|
||
b = grid_radius_circles(prep, lim, gx, gz, grid[2])
|
||
coarse = tuple(np.concatenate([a[k], b[k]]) for k in range(5))
|
||
cheap = tuple(m for m in methods if m not in ("spencer", "morgenstern_price"))
|
||
cheap = cheap or ("bishop",)
|
||
F = _circle_eval(prep, *coarse, coarse_slices, min_depth, cheap)
|
||
counts["circles"] = int(len(coarse[0]))
|
||
counts["circles_valid"] = int(np.sum(np.isfinite(F[cheap[0]])))
|
||
# 등안전율 격자 — 격자 중심마다 최소(Bishop 우선)
|
||
key = "bishop" if "bishop" in F else cheap[0]
|
||
fb = F[key][len(a[0]) :]
|
||
cmin = np.full(len(gx) * len(gz), np.nan)
|
||
for o, f in zip(b[5], fb):
|
||
if np.isfinite(f) and not (f >= cmin[o]):
|
||
cmin[o] = f
|
||
contour = {
|
||
"method": key,
|
||
"x": [float(s * v) for v in gx],
|
||
"z": [float(v) for v in gz],
|
||
"F": [
|
||
[None if not np.isfinite(v) else float(v) for v in row]
|
||
for row in cmin.reshape(len(gz), len(gx))
|
||
],
|
||
}
|
||
smap = safety_map(prep, coarse, F[key], key)
|
||
if "spencer" in methods:
|
||
# Spencer 씨앗 — 거친 탐색 상위 원에서 Spencer 를 따로 셈(풀리는 원이 조각나 있어 Bishop 순위만으론 놓침)
|
||
top = np.argsort(np.where(np.isfinite(F[key]), F[key], np.inf))[:SPENCER_SEED_POOL]
|
||
top = top[np.isfinite(F[key][top])]
|
||
F["spencer"] = np.full(len(coarse[0]), np.nan)
|
||
if len(top):
|
||
c5 = [v[top] for v in coarse[:5]]
|
||
sl = circle_slices(prep, *c5, coarse_slices, min_depth)
|
||
center = np.column_stack([c5[0], c5[1]])
|
||
F["spencer"][top] = evaluate(sl, center, ("spencer",), scan=False)["spencer"]["F"]
|
||
searched = [m for m in methods if m != "morgenstern_price"]
|
||
pool, tops = [], {}
|
||
for m in searched:
|
||
score = F[m] if m in F else F[key]
|
||
order = np.argsort(np.where(np.isfinite(score), score, np.inf))[:30]
|
||
order = order[np.isfinite(score[order])]
|
||
if not len(order):
|
||
continue
|
||
seeds = np.column_stack([coarse[3][order], coarse[4][order], _theta(coarse, order)])
|
||
if m == "spencer" and "bishop" in tops:
|
||
# Spencer 가 풀리는 원은 Bishop 임계 둘레에 좁게 · 조각나 있음 — 그 둘레 구름에서 풀리는 원을 씨앗에 보탬
|
||
cloud = _cloud(prep, lim, m, tops["bishop"][0], n_slices, min_depth, seed)
|
||
seeds = np.vstack([seeds[:20], cloud, tops["bishop"]])
|
||
fmin, geom, near, top = refine_circle(prep, lim, m, seeds, n_slices, min_depth, seed)
|
||
if fmin < _BIG:
|
||
best_circles[m] = geom
|
||
pool += [geom, *near]
|
||
tops[m] = top
|
||
# 방법마다 찾은 원(+ 마지막 무리 상위)을 모든 방법으로 다시 셈 → 방법별 가장 낮은 면
|
||
# (한 방법 탐색이 놓친 골짜기 메움 · Spencer 는 풀리는 원이 조각나 있어 특히 필요)
|
||
for m in methods:
|
||
# M-P(차례 계산 · 느림)는 다른 방법이 고른 원에서만
|
||
geoms = list(best_circles.values()) if m == "morgenstern_price" else pool
|
||
g = _cross(prep, "circle", geoms, m, n_slices, min_depth)
|
||
if g is not None:
|
||
best_circles[m] = g
|
||
results[m]["circular"] = _describe(prep, "circle", g, m, n_slices, min_depth)
|
||
# 엄밀해(Spencer · M-P)가 Bishop 임계 원에서 해 없음 · 찾은 최소가 Bishop + FAR_TOL 넘게 높음 = 임계
|
||
# 둘레에서 풀리는 원이 없어 먼 원을 고른 것 → 그 방법 원호 = 「해 없음」(Bishop 임계 원 · F None · 까닭 ·
|
||
# F–λ 곡선) · 허용 안이면 찾은 값을 살림(status converged_far)
|
||
ref = best_circles.get("bishop") if "bishop" in results else None
|
||
fb = (results.get("bishop", {}).get("circular") or {}).get("F")
|
||
for m in [k for k in ("spencer", "morgenstern_price") if k in methods and ref and fb]:
|
||
cur = results[m].get("circular")
|
||
if cur and cur["F"] is not None and cur["F"] <= fb:
|
||
continue
|
||
at = _describe(prep, "circle", ref, m, n_slices, min_depth)
|
||
if at["F"] is not None:
|
||
continue
|
||
if cur and cur["F"] is not None and cur["F"] <= fb * (1 + FAR_TOL):
|
||
cur["status"] = "converged_far"
|
||
cur["reason"] = (
|
||
f"Bishop 임계 원에서 해 없음 · 다른 원 최소가 Bishop 의 +{FAR_TOL:.0%} 안이라 살림"
|
||
)
|
||
else:
|
||
results[m]["circular"] = at
|
||
if noncircular:
|
||
mv = polyline_vertices
|
||
seeds = []
|
||
circle = (
|
||
best_circles.get("janbu")
|
||
or best_circles.get("bishop")
|
||
or next(iter(best_circles.values()), None)
|
||
)
|
||
if circle:
|
||
seeds.append(circle_to_params(prep, circle, mv))
|
||
seeds += block_seeds(prep, lim, mv)
|
||
seeds = _best_seeds(prep, seeds, mv, n_slices, min_depth, 40)
|
||
# 힘 평형 방법(Janbu) 먼저 · Spencer 는 그 마지막 무리에서 이어 찾음 · 모멘트 방법(Fellenius ·
|
||
# Bishop)은 비원호 중심이 임의라 따로 최적화하지 않고 찾은 면에서 값만 셈(참고 칸 · 상용 SW 도
|
||
# 비원호 최적화 = Janbu · Spencer · GLE) · M-P 는 찾은 면에서 값만
|
||
found = []
|
||
for m in [k for k in ("janbu", "spencer") if k in methods]:
|
||
iters = LS_ITER_SPENCER if m == "spencer" else LS_ITER
|
||
fmin, _x, geom, top = optimize_polyline(
|
||
prep, lim, m, seeds, n_slices, min_depth, mv, seed, iters
|
||
)
|
||
seeds = [*top, *seeds]
|
||
if fmin < _BIG:
|
||
found.append(geom)
|
||
for m in methods:
|
||
g = _cross(prep, "polyline", found, m, n_slices, min_depth)
|
||
if g is not None:
|
||
results[m]["noncircular"] = _describe(prep, "polyline", g, m, n_slices, min_depth)
|
||
for m in methods:
|
||
# 비원호 면 위 Fellenius · Bishop = 모멘트 중심 임의 → 참고 칸(최소 F 후보에서 뺌 · Slide2 와 같음)
|
||
if m in ("ordinary", "bishop") and results[m].get("noncircular"):
|
||
results[m]["noncircular"]["reference_only"] = True
|
||
cands = [
|
||
v
|
||
for k, v in results[m].items()
|
||
if v.get("F") is not None and not v.get("reference_only")
|
||
]
|
||
results[m]["min"] = min(cands, key=lambda v: v["F"]) if cands else None
|
||
out = {
|
||
"direction": s,
|
||
"limits": sorted([s * lim[0], s * lim[1]]),
|
||
"methods": results,
|
||
"contour": contour,
|
||
"safety_map": smap,
|
||
"counts": counts,
|
||
"elapsed_s": round(time.perf_counter() - t0, 3),
|
||
}
|
||
if prep.supports is not None:
|
||
# 보강력 결과 = 판정 임계면(방법별 최소 중 가장 낮은 면) 기준
|
||
mins = [(m, v["min"]) for m, v in results.items() if v.get("min")]
|
||
if mins:
|
||
m, gov = min(mins, key=lambda kv: kv[1]["F"])
|
||
out["support_forces"] = gov.get("support_forces", [])
|
||
out["support_forces_method"] = m
|
||
else:
|
||
out["support_forces"], out["support_forces_method"] = [], None
|
||
return out
|
||
|
||
|
||
def _stack(prep, kind, geoms, n, min_depth):
|
||
if kind == "circle":
|
||
xc, zc, R, xE, xX = (np.array([g[k] for g in geoms], float) for k in range(5))
|
||
return circle_slices(prep, xc, zc, R, xE, xX, n, min_depth), np.column_stack([xc, zc])
|
||
vx = np.array([g[0] for g in geoms], float)
|
||
vz = np.array([g[1] for g in geoms], float)
|
||
return poly_slices(prep, vx, vz, n, min_depth), three_point_center(vx, vz)
|
||
|
||
|
||
def _cross(prep, kind, geoms, method, n, min_depth):
|
||
"""찾은 면들 위 method 안전율을 한 번에 → 가장 낮은 면(없으면 None)."""
|
||
if not geoms:
|
||
return None
|
||
sl, center = _stack(prep, kind, geoms, n, min_depth)
|
||
F = evaluate(sl, center, (method,))[method]["F"]
|
||
if not np.isfinite(F).any():
|
||
return None
|
||
return geoms[int(np.nanargmin(F))]
|
||
|
||
|
||
def _describe(prep, kind, geom, method, n, min_depth):
|
||
"""M-P 는 같은 면 Spencer 값을 출발 참고값으로(엉뚱한 근 거름)."""
|
||
seed = None
|
||
if method == "morgenstern_price":
|
||
sl, center = _stack(prep, kind, [geom], n, min_depth)
|
||
sp = evaluate(sl, center, ("spencer",))["spencer"]["F"][0]
|
||
seed = float(sp) if np.isfinite(sp) else None
|
||
return describe(prep, kind, geom, method, n, min_depth, seed)
|
||
|
||
|
||
def _theta(coarse, order):
|
||
xc, zc, R, xE, xX = (v[order] for v in coarse)
|
||
zE = zc - np.sqrt(np.clip(R**2 - (xE - xc) ** 2, 0, None))
|
||
zX = zc - np.sqrt(np.clip(R**2 - (xX - xc) ** 2, 0, None))
|
||
chord = np.hypot(xX - xE, zX - zE)
|
||
return 2 * np.arcsin(np.clip(chord / (2 * R), 0, 1))
|
||
|
||
|
||
def evaluate_specified(
|
||
model: dict, *, circle=None, polyline=None, direction=1, methods=METHODS, n_slices=50
|
||
):
|
||
"""지정 활동면 하나(원 [xc, zc, R] 또는 꺾은선 [[x, z], …])의 방법별 안전율 · 절편표 — 대조 · 검산용."""
|
||
prep = lem.Prepared(model, direction)
|
||
s = prep.direction
|
||
out = {}
|
||
if circle is not None:
|
||
xc, zc, R = s * float(circle[0]), float(circle[1]), float(circle[2])
|
||
xE, xX, ok = circle_cross_surface(prep, np.array([xc]), np.array([zc]), np.array([R]))
|
||
if not ok[0]:
|
||
raise ValueError("지정 원이 지표면을 지나지 않음")
|
||
geom, kind = (xc, zc, R, float(xE[0]), float(xX[0])), "circle"
|
||
else:
|
||
pts = sorted((s * float(x), float(z)) for x, z in polyline)
|
||
geom, kind = (np.array([p[0] for p in pts]), np.array([p[1] for p in pts])), "polyline"
|
||
for m in methods:
|
||
out[m] = _describe(prep, kind, geom, m, n_slices, 0.0)
|
||
return out
|