From 87928cb26739a3e75081ddb851f280b9615948f4 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Mon, 26 Oct 2026 03:18:35 +0900 Subject: [PATCH] auto: 2026-10-26 03:18 (ESD_LAPTOP) --- .../E03_StabSection_Engine_Describe.py | 114 +++++ E03_StabSection/E03_StabSection_Engine_GLE.py | 198 ++++++++ E03_StabSection/E03_StabSection_Engine_LEM.py | 334 ++++++++----- .../E03_StabSection_Engine_Search.py | 459 ++++++++---------- .../E03_StabSection_Engine_Station.py | 35 +- .../E03_StabSection_Engine_Surface.py | 185 +++++++ .../tester/test_532_2_stab_lem_benchmarks.py | 10 +- .../tester/test_532_2_stab_mp_outputs.py | 141 ++++++ 8 files changed, 1101 insertions(+), 375 deletions(-) create mode 100644 E03_StabSection/E03_StabSection_Engine_Describe.py create mode 100644 E03_StabSection/E03_StabSection_Engine_GLE.py create mode 100644 E03_StabSection/E03_StabSection_Engine_Surface.py create mode 100644 resources/tester/test_532_2_stab_mp_outputs.py diff --git a/E03_StabSection/E03_StabSection_Engine_Describe.py b/E03_StabSection/E03_StabSection_Engine_Describe.py new file mode 100644 index 000000000..0ca4ba4a0 --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_Describe.py @@ -0,0 +1,114 @@ +"""E03 결과 꾸미기 — 임계면 글상자 · 절편 자유물체도 · F–λ 곡선 · Safety Map 상위 면. + +좌표는 원래 방향(x = B06 offset · 좌 +)으로 되돌려 냄. 출력 칸은 덧붙이기만(E03 화면 · E05 보고서가 읽음). +""" + +import numpy as np + +from E03_StabSection import E03_StabSection_Engine_GLE as gle +from E03_StabSection import E03_StabSection_Engine_LEM as lem +from E03_StabSection.E03_StabSection_Engine_Surface import ( + circle_slices, + evaluate, + poly_slices, + three_point_center, +) + +#: Safety Map 면 수 · 면마다 꼭짓점 수(원호) +SAFETY_MAP_TOP = 50 +_MAP_POINTS = 25 + + +def _curve(sl, center, method, f=None) -> dict: + """F_m · F_f 대 λ (화면 그래프) — Spencer = 상수 f · M-P = 반정현.""" + if method == "spencer": + grid, fm, ff = lem.lambda_curve( + sl, center, np.linspace(*lem.LAMBDA_RANGE, gle.CURVE_POINTS) + ) + fm, ff = fm[:, 0], ff[:, 0] + else: + grid, fm, ff = gle.lambda_curves(sl, center, f) + fm, ff = fm[0], ff[0] + num = lambda v: [float(x) if np.isfinite(x) else None for x in v] # noqa: E731 + return {"lambda": num(grid), "Fm": num(fm), "Ff": num(ff)} + + +def describe(prep, kind, geom, method, n, min_depth, seed_F=None) -> dict: + """임계면 하나 → 안전율 · 좌표 · 절편표(자유물체도 힘) · 수렴 · (Spencer · M-P) λ · F–λ 곡선.""" + s = prep.direction + if kind == "circle": + xc, zc, R, xE, xX = geom + sl = circle_slices(prep, [xc], [zc], [R], [xE], [xX], n, min_depth) + center = np.array([[xc, zc]]) + pts = np.column_stack([sl["xe"][0], sl["ze"][0]]) + else: + vx, vz = geom + sl = poly_slices(prep, vx[None, :], vz[None, :], n, min_depth) + center = three_point_center(vx[None, :], vz[None, :]) + pts = np.column_stack([vx, vz]) + seed = None if seed_F is None else np.array([seed_F], float) + res = evaluate(sl, center, (method,), mp_seed=seed)[method] + F = float(res["F"][0]) + lam = float(res["lambda"][0]) if method in ("spencer", "morgenstern_price") else 0.0 + f = gle.interslice_function(sl["xe"]) if method == "morgenstern_price" else None + rows = lem.slice_table(sl, 0, F, lam, method, f) if np.isfinite(F) else [] + for r in rows: + r["material"] = ( + prep.materials[r["material"]]["id"] if r["material"] < len(prep.materials) else None + ) + # 행 순서 = 머리 → 끝 · E · X 의 left = 머리 쪽 경계 · right = 끝 쪽 경계(해석 틀 그대로) + r["x_left"], r["x_right"] = sorted((s * r["x_left"], s * r["x_right"])) + r["x_mid"] = s * r["x_mid"] + out = { + "F": F if np.isfinite(F) else None, + "kind": kind, + "points": [[float(s * x), float(z)] for x, z in (pts if s > 0 else pts[::-1])], + "entry": [float(s * sl["xe"][0, 0]), float(sl["ze"][0, 0])], + "exit": [float(s * sl["xe"][0, -1]), float(sl["ze"][0, -1])], + "iterations": int(res["iterations"][0]), + "converged": bool(res["converged"][0]), + "slices": rows, + } + if kind == "circle": + out["center"] = [float(s * center[0, 0]), float(center[0, 1])] + out["radius"] = float(geom[2]) + else: + out["moment_center"] = [float(s * center[0, 0]), float(center[0, 1])] + if method in ("spencer", "morgenstern_price"): + out["lambda"] = lam + out["theta_deg"] = float(np.degrees(np.arctan(lam))) + if np.isfinite(F): + out["lambda_curve"] = _curve(sl, center, method, f) + if method == "morgenstern_price": + out["interslice_function"] = "half_sine" + if method == "janbu": + out["F0"] = float(res["F0"][0]) if np.isfinite(res["F0"][0]) else None + out["f0"] = float(res["f0"][0]) + if prep.crack_depth > 0: + out["crack_water_force"] = float(sl["Pw"][0]) + return out + + +def safety_map(prep, circles, F, method, top=SAFETY_MAP_TOP) -> dict: + """Safety Map — 거친 탐색 원호 중 안전율 낮은 top 개(F · 중심 · 반지름 · 진입 · 진출 · 꼭짓점).""" + s = prep.direction + xc, zc, R, xE, xX = circles + order = np.argsort(np.where(np.isfinite(F), F, np.inf))[:top] + order = order[np.isfinite(F[order])] + t = np.linspace(0.0, 1.0, _MAP_POINTS) + out = [] + for k in order: + x = xE[k] + (xX[k] - xE[k]) * t + z = zc[k] - np.sqrt(np.clip(R[k] ** 2 - (x - xc[k]) ** 2, 0, None)) + pts = np.column_stack([s * x, z]) + out.append( + { + "F": float(F[k]), + "center": [float(s * xc[k]), float(zc[k])], + "radius": float(R[k]), + "entry": [float(s * xE[k]), float(z[0])], + "exit": [float(s * xX[k]), float(z[-1])], + "points": [[float(a), float(b)] for a, b in (pts if s > 0 else pts[::-1])], + } + ) + return {"method": method, "surfaces": out} diff --git a/E03_StabSection/E03_StabSection_Engine_GLE.py b/E03_StabSection/E03_StabSection_Engine_GLE.py new file mode 100644 index 000000000..65f4ddf56 --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_GLE.py @@ -0,0 +1,198 @@ +"""E03 엄밀해 — Morgenstern-Price(GLE) · 절편간 힘 X = λ·f(x)·E · f = 반정현(기본) 또는 상수. + +식 = LEM 과 같은 절편 평형(`E03_StabSection_Engine_LEM` 머리글) — 절편간 힘 함수만 경계마다 다름: + 절편 i 연직 평형 → N·(m_α − λ·f_{i+1}·A) = W − c'·sinα/F − X_i + λ·f_{i+1}·(E_i + B0) + 수평 평형 → E_{i+1} = E_i + B0 + N·A · X_{i+1} = λ·f_{i+1}·E_{i+1} · 머리 E_0 = 균열 수압 · X_0 = 0 + f ≡ 1 이면 Spencer 와 같은 식(시험으로 확인). F_m · F_f 는 LEM 과 같은 전체 평형식. +풀이 = λ 훑기(F_m(λ) · F_f(λ) 곡선 — 화면 그래프 자료) → 교점 → (F, λ) 2 차 뉴턴(수치 야코비안). +계산은 임계면 몇 개만(탐색 고리 밖) — 면을 λ 갯수만큼 복제해 한 번에 셈. 외부 코드 없음(NumPy 만). +""" + +import numpy as np + +from E03_StabSection import E03_StabSection_Engine_LEM as lem + +FUNCTIONS = ("half_sine", "constant") +_TINY = 1e-9 +CURVE_POINTS = 21 +#: M-P λ 범위 — 반정현 평균 2/π 라 평균 X/E 가 Spencer 범위(LEM.LAMBDA_RANGE)와 같게 π/2 배 +LAMBDA_RANGE = tuple(v * np.pi / 2 for v in lem.LAMBDA_RANGE) + + +def interslice_function(xe: np.ndarray, kind: str = "half_sine") -> np.ndarray: + """절편 경계마다 f(x) — 반정현 sin(π·(x − x_진입)/(x_진출 − x_진입)) · 상수 1.""" + if kind == "constant": + return np.ones_like(xe) + if kind != "half_sine": + raise ValueError(f"절편간 힘 함수 모름: {kind}") + x0, x1 = xe[:, :1], xe[:, -1:] + return np.sin(np.pi * (xe - x0) / np.maximum(x1 - x0, _TINY)) + + +def forces(sl: dict, F, lam, f): + """절편마다 N · 분모 · 경계 E · X (꼴 (S, n) · (S, n) · (S, n+1) · (S, n+1)) — 머리부터 차례로.""" + p = lem._pre(sl) + F = np.asarray(F, float) + lam = np.asarray(lam, float) + g = 1.0 / F[:, None] + m = sl["ca"] + g * p["sat"] + A = sl["sa"] - g * p["tca"] + B0 = p["khws"] - g * p["cuca"] + rhs = sl["W"] - g * p["cusa"] + S, n = rhs.shape + N = np.empty((S, n)) + den = np.empty((S, n)) + E = np.empty((S, n + 1)) + X = np.zeros((S, n + 1)) + E[:, 0] = sl["Pw"] + for i in range(n): + lf = lam * f[:, i + 1] + d = m[:, i] - lf * A[:, i] + den[:, i] = d + N[:, i] = (rhs[:, i] - X[:, i] + lf * (E[:, i] + B0[:, i])) / np.where( + np.abs(d) < _TINY, _TINY, d + ) + E[:, i + 1] = E[:, i] + B0[:, i] + N[:, i] * A[:, i] + X[:, i + 1] = lf * E[:, i + 1] + return N, den, E, X + + +def _fm_ff(sl, center, F, lam, f): + N, den, *_ = forces(sl, F, lam, f) + at, gn, d0, ma = lem._mom(sl, center) + p = lem._pre(sl) + fm = lem._ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1)) + ff = lem._ratio(p["rc0"] + (p["tca"] * N).sum(1), p["fd0"] + (sl["sa"] * N).sum(1), neg=False) + return fm, ff, den + + +def _repeat(sl: dict, k: int) -> dict: + """면마다 k 번 복제(λ 훑기 · 야코비안을 한 번에) — 미리 묶은 칸은 버림.""" + S = len(sl["W"]) + return { + key: (np.repeat(v, k, axis=0) if isinstance(v, np.ndarray) and v.shape[:1] == (S,) else v) + for key, v in sl.items() + if key not in ("_pre", "_mom") + } + + +def _solve_f(fun, F0, iters=60): + """F = fun(F) — 수치 미분 뉴턴 · 안 되면 반쯤 섞은 고정점(LEM `_iterate` 와 같은 안전장치).""" + F = np.asarray(F0, float).copy() + done = np.zeros(len(F), bool) + for _ in range(iters): + f1 = fun(F) + h1 = f1 - F + dlt = 1e-6 * np.maximum(np.abs(F), 0.1) + dh = ((fun(F + dlt) - (F + dlt)) - h1) / dlt + newton = F - h1 / np.where(np.abs(dh) < 1e-12, -1.0, dh) + use = np.isfinite(newton) & (dh < -1e-9) & (newton > 0.5 * F) & (newton < 2.0 * F) + Fn = np.where(use, newton, 0.5 * (F + f1)) + Fn = np.where(np.isfinite(Fn) & (Fn > 0), Fn, np.nan) + conv = np.abs(h1) < lem.TOL * np.maximum(1.0, np.abs(F)) + F = np.where(done | conv, F, Fn) + done |= conv | ~np.isfinite(F) + if done.all(): + break + return np.where(done, F, np.nan) + + +def lambda_curves(sl, center, f, grid=None, F0=None): + """λ 마다 F_m(λ) · F_f(λ) — 꼴 (G,) · (S, G) · (S, G).""" + grid = np.linspace(*LAMBDA_RANGE, CURVE_POINTS) if grid is None else np.asarray(grid) + S, G = len(sl["W"]), len(grid) + rep = _repeat(sl, G) + rc = np.repeat(center, G, axis=0) + rf = np.repeat(f, G, axis=0) + lam = np.tile(grid, S) + start = np.full(S * G, 1.5) if F0 is None else np.repeat(F0, G) + fm = _solve_f(lambda F: _fm_ff(rep, rc, F, lam, rf)[0], start) + ff = _solve_f(lambda F: _fm_ff(rep, rc, F, lam, rf)[1], start) + return grid, fm.reshape(S, G), ff.reshape(S, G) + + +def _crossing(grid, fm, ff): + """F_m − F_f 부호가 바뀌는 곳 중 λ = 0 가까운 쪽 → (F, λ) · 없으면 nan.""" + S = len(fm) + F, L = np.full(S, np.nan), np.full(S, np.nan) + best = np.full(S, np.inf) + d = fm - ff + for k in range(len(grid) - 1): + d0, d1 = d[:, k], d[:, k + 1] + cross = np.isfinite(d0) & np.isfinite(d1) & (np.sign(d0) != np.sign(d1)) + t = np.where(cross, d0 / np.where(d0 == d1, 1, d0 - d1), 0) + lc = grid[k] + t * (grid[k + 1] - grid[k]) + fc = fm[:, k] + t * (fm[:, k + 1] - fm[:, k]) + better = cross & (np.abs(lc) < best) + best = np.where(better, np.abs(lc), best) + F, L = np.where(better, fc, F), np.where(better, lc, L) + return F, L + + +def _newton2(sl, center, f, F, lam, iters=40): + """(F_m − F, F_f − F) = 0 — 세 벌 복제(F · F+h · λ+h)로 수치 야코비안 · 걸음 제한 LEM Spencer 와 같음.""" + S = len(F) + rep, rc, rf = _repeat(sl, 3), np.repeat(center, 3, axis=0), np.repeat(f, 3, axis=0) + F, lam = F.copy(), lam.copy() + done = ~np.isfinite(F) + it = np.zeros(S, int) + for _ in range(iters): + hF = 1e-6 * np.maximum(np.abs(F), 0.1) + hl = 1e-6 + Fs = np.column_stack([F, F + hF, F]).ravel() + ls = np.column_stack([lam, lam, lam + hl]).ravel() + fm, ff, _ = _fm_ff(rep, rc, np.where(np.isfinite(Fs), Fs, 1.0), ls, rf) + rm = (fm - Fs).reshape(S, 3) + rf_ = (ff - Fs).reshape(S, 3) + j11, j12 = (rm[:, 1] - rm[:, 0]) / hF, (rm[:, 2] - rm[:, 0]) / hl + j21, j22 = (rf_[:, 1] - rf_[:, 0]) / hF, (rf_[:, 2] - rf_[:, 0]) / hl + det = j11 * j22 - j12 * j21 + det = np.where(np.abs(det) < 1e-14, np.nan, det) + dF = -(j22 * rm[:, 0] - j12 * rf_[:, 0]) / det + dl = -(-j21 * rm[:, 0] + j11 * rf_[:, 0]) / det + scale = np.minimum( + 1.0, + np.minimum( + 0.5 * np.abs(F) / np.maximum(np.abs(dF), _TINY), + 0.3 / np.maximum(np.abs(dl), _TINY), + ), + ) + F = np.where(done, F, F + scale * dF) + lam = np.where(done, lam, lam + scale * dl) + it += ~done + small = (np.abs(dF) < lem.TOL * np.maximum(1, np.abs(F))) & (np.abs(dl) < lem.TOL) + done |= small | ~np.isfinite(F) | (F < lem.F_MIN) + if done.all(): + break + return F, lam, it + + +def morgenstern_price(sl: dict, center: np.ndarray, kind="half_sine", F0=None, curve=True): + """M-P (F, λ) · 수렴 · (curve) F_m · F_f–λ 곡선. F0 = 출발 참고값(Spencer · Bishop) — 엉뚱한 근 거름.""" + f = interslice_function(sl["xe"], kind) + grid, fm, ff = lambda_curves(sl, center, f, F0=F0) + F, lam = _crossing(grid, fm, ff) + F, lam, it = _newton2(sl, center, f, F, lam) + ok = np.isfinite(F) + safe = np.where(ok & (F > 0), F, 1.0) + rm, rff, den = _fm_ff(sl, center, safe, np.where(ok, lam, 0.0), f) + p = lem._pre(sl) + N = forces(sl, safe, np.where(ok, lam, 0.0), f)[0] + conv = ok & (np.abs(rm - F) < 1e-4) & (np.abs(rff - F) < 1e-4) + conv &= (lam >= LAMBDA_RANGE[0]) & (lam <= LAMBDA_RANGE[1]) + if F0 is not None: + conv &= (F > 0.5 * F0) & (F < 2.0 * F0) + malpha_ok = np.all(den >= lem.MALPHA_MIN, axis=1) + R = p["cu"] + N * sl["tanphi"] + good = sl["valid"] & conv & (F > lem.F_MIN) & malpha_ok & (R.sum(1) > 0) + out = { + "F": np.where(good, F, np.nan), + "lambda": np.where(good, lam, np.nan), + "iterations": it, + "converged": conv, + "malpha_ok": malpha_ok, + "function": kind, + } + if curve: + out["curve"] = (grid, fm, ff) + return out diff --git a/E03_StabSection/E03_StabSection_Engine_LEM.py b/E03_StabSection/E03_StabSection_Engine_LEM.py index b68ec84ba..79e755295 100644 --- a/E03_StabSection/E03_StabSection_Engine_LEM.py +++ b/E03_StabSection/E03_StabSection_Engine_LEM.py @@ -25,6 +25,8 @@ LAMBDA_RANGE = (-1.25, 0.25) #: 이보다 작은 F = 1/F 항이 터진 엉뚱한 근(실제 비탈 F 는 이 위) — 버림. F_MIN = 0.05 TOL = 1e-6 +#: Janbu 보정 f0 식의 d/L 상한(지식DB 02 §2.3) +JANBU_DL_MAX = 0.35 MAX_ITER = 100 @@ -116,6 +118,12 @@ class Prepared: # 빈 칸 = 아주 낮은 두께 0 구간(nan 처리 없이 바로 셈) for name in ("bot_a", "bot_b", "top_a", "top_b"): setattr(self, name, np.nan_to_num(getattr(self, name), nan=_PAD_Z)) + # 절편 계산이 한 번에 모아 읽는 표 (띠, 6, 구간) — 바닥 · 바닥 기울기 · 윗면 · 윗면 기울기 · γt · γsat − γt + gt, gs = self.gamma_t[self.mat], self.gamma_sat[self.mat] + self.table = np.stack( + [self.bot_a, self.bot_b - self.bot_a, self.top_a, self.top_b - self.top_a, gt, gs - gt], + axis=1, + ) # ---- 조회 ------------------------------------------------------------- def surface(self, x: np.ndarray) -> np.ndarray: @@ -147,25 +155,36 @@ def slices(prep: Prepared, xe: np.ndarray, ze: np.ndarray) -> dict: xm = 0.5 * (xe[:, :-1] + xe[:, 1:]) zb = 0.5 * (ze[:, :-1] + ze[:, 1:]) alpha = np.arctan2(ze[:, :-1] - ze[:, 1:], np.maximum(b, _TINY)) - bot, top, mat = prep.columns(xm) + j = np.clip(np.searchsorted(prep.X, xm, side="right") - 1, 0, len(prep.xa) - 1) + t = ((xm - prep.xa[j]) / (prep.xb[j] - prep.xa[j]))[..., None] + G = prep.table[j] + mat = prep.mat[j] + bot = G[..., 0, :] + t * G[..., 1, :] + top = G[..., 2, :] + t * G[..., 3, :] + gt, dgs = G[..., 4, :], G[..., 5, :] lo = np.maximum(bot, zb[..., None]) thick = np.clip(top - lo, 0.0, None) zt = prep.surface(xm) if prep.gw_mode == "line": zw = np.minimum(prep.water(xm), zt) wet = np.clip(np.minimum(top, zw[..., None]) - lo, 0.0, None) - dry = thick - wet u = GAMMA_W * np.clip(zw - zb, 0.0, None) - wet_mid = lo + 0.5 * wet - dry_mid = lo + wet + 0.5 * dry - g_dry, g_wet = prep.gamma_t[mat] * dry, prep.gamma_sat[mat] * wet + weight_per_b = np.sum(gt * thick + dgs * wet, axis=-1) else: - g_dry, g_wet = prep.gamma_t[mat] * thick, 0.0 * thick - dry_mid = wet_mid = lo + 0.5 * thick + wet = None u = np.zeros_like(zb) - weight_per_b = np.sum(g_dry + g_wet, axis=-1) - moment_z = np.sum(g_dry * dry_mid + g_wet * wet_mid, axis=-1) - zg = np.where(weight_per_b > 0, moment_z / np.maximum(weight_per_b, _TINY), zb) + weight_per_b = np.sum(gt * thick, axis=-1) + if prep.kh > 0: # 무게중심 높이는 지진 관성력 팔에만 쓰임 + if wet is None: + moment_z = np.sum(gt * thick * (lo + 0.5 * thick), axis=-1) + else: + dry = thick - wet + moment_z = np.sum( + gt * dry * (lo + wet + 0.5 * dry) + (gt + dgs) * wet * (lo + 0.5 * wet), axis=-1 + ) + zg = np.where(weight_per_b > 0, moment_z / np.maximum(weight_per_b, _TINY), zb) + else: + zg = zb Ws = weight_per_b * b if prep.gw_mode == "ru": u = prep.ru * weight_per_b @@ -225,77 +244,130 @@ def add_crack_water(sl: dict, depth: np.ndarray, fill: float) -> None: # ---- 방법 ------------------------------------------------------------------- -def _parts(sl: dict, F, lam): - """주어진 F · λ 의 절편 수직력 N · 저항 R · 분모(m_α − λA). - - 유효 수직력 음(인장) 절편도 그대로 둠(Slide2 와 같음 — 0 으로 자르면 ACADS · Arai & Tagyo 대조가 3~6 % 벗어남). - 대신 저항 합이 0 이하인 면은 `_finish` 가 버림. - """ - F = np.asarray(F, float)[:, None] - lam = np.asarray(lam, float)[:, None] - W = sl["W"].copy() - W[:, 0] += lam[:, 0] * sl["Pw"] # 머리 경계 수압은 전단 없음 → X_0 = 0 - t = sl["tanphi"] - N, den, A, B0 = _normal(sl, F, lam, W, t) - R = sl["c"] * sl["l"] + (N - sl["u"] * sl["l"]) * t - return N, R, den, A, B0 +# 한 벌 식을 F 와 무관한 칸(_pre) · 모멘트 중심 칸(_mom) 으로 미리 묶어 두고, F · λ 마다는 N 만 새로 셈. +# g = 1/F 로 놓으면 N = (P0 − g·P1)/(Q0 + g·Q1) — dN/dg · dN/dλ 가 닫힌 꼴이라 뉴턴 야코비안을 해석적으로 씀. +def _pre(sl: dict) -> dict: + p = sl.get("_pre") + if p is None: + t, sa, ca = sl["tanphi"], sl["sa"], sl["ca"] + cu = sl["c"] * sl["l"] - sl["u"] * sl["l"] * t + khws = sl["kh"] * sl["Ws"] + lw = khws.copy() + lw[:, 0] += sl["Pw"] # 머리 경계 수압은 전단 없음 → X_0 = 0 + p = { + "cu": cu, + "cusa": cu * sa, + "cuca": cu * ca, + "sat": sa * t, + "tca": t * ca, + "khws": khws, + "lw": lw, + "fd0": khws.sum(1) + sl["Pw"], + "rc0": (cu * ca).sum(1), + } + sl["_pre"] = p + return p -def _normal(sl, F, lam, W, t): - sa, ca = sl["sa"], sl["ca"] - cu = sl["c"] * sl["l"] - sl["u"] * sl["l"] * t - m = ca + sa * t / F - A = sa - t * ca / F - B0 = sl["kh"] * sl["Ws"] - cu * ca / F - den = m - lam * A - N = (W + lam * B0 - cu * sa / F) / np.where(np.abs(den) < _TINY, _TINY, den) - return N, den, A, B0 - - -def _moment_terms(sl: dict, center: np.ndarray, N: np.ndarray): +def _mom(sl: dict, center: np.ndarray): + """모멘트 팔(a) · 분모 N 계수(gN) · 상수 칸 — 같은 중심 배열이면 다시 쓰고 아니면 새로 셈.""" + cached = sl.get("_mom") + if cached is not None and cached[0] is center: + return cached[1] + p = _pre(sl) xo, zo = center[:, 0:1], center[:, 1:2] - a = (sl["xm"] - xo) * sl["sa"] + (sl["zb"] - zo) * sl["ca"] - D = ( - (sl["xw"] - xo) * sl["W"] - + (sl["zg"] - zo) * sl["kh"] * sl["Ws"] - - (sl["xm"] - xo) * N * sl["ca"] - + (sl["zb"] - zo) * N * sl["sa"] - ).sum(1) + (sl["zp"] - center[:, 1]) * sl["Pw"] - return a, D + dx, dz = sl["xm"] - xo, sl["zb"] - zo + a = dx * sl["sa"] + dz * sl["ca"] + gn = dz * sl["sa"] - dx * sl["ca"] + d0 = ((sl["xw"] - xo) * sl["W"] + (sl["zg"] - zo) * p["khws"]).sum(1) + d0 = d0 + (sl["zp"] - center[:, 1]) * sl["Pw"] + out = (a * sl["tanphi"], gn, d0, (a * p["cu"]).sum(1)) + sl["_mom"] = (center, out) + return out -def _force_den(sl: dict, N: np.ndarray): - return (N * sl["sa"] + sl["kh"] * sl["Ws"]).sum(1) + sl["Pw"] +def _normal(sl, F, lam, grad=0): + """절편 수직력 N · 분모(m_α − λA) · grad ≥ 1 dN/dg · grad 2 dN/dλ — F · λ 꼴 (S,).""" + p = _pre(sl) + g = 1.0 / np.asarray(F, float)[:, None] + L = np.asarray(lam, float)[:, None] + P1 = p["cusa"] + L * p["cuca"] + Q1 = p["sat"] + L * p["tca"] + den = sl["ca"] - L * sl["sa"] + g * Q1 + safe = np.where(np.abs(den) < _TINY, _TINY, den) + N = (sl["W"] + L * p["lw"] - g * P1) / safe + if not grad: + return N, den + dg = -(P1 + N * Q1) / safe + if grad == 1: + return N, den, dg, None + dl = (p["lw"] - g * p["cuca"] + N * (sl["sa"] - g * p["tca"])) / safe + return N, den, dg, dl -def _fm(sl, center, F, lam): - N, R, den, *_ = _parts(sl, F, lam) - a, D = _moment_terms(sl, center, N) - return (a * R).sum(1) / np.where(np.abs(D) < _TINY, -_TINY, D), den +def _parts(sl: dict, F, lam): + """주어진 F · λ 의 N · 저항 R · 분모 · A · B0(절편표용). + + 유효 수직력 음(인장) 절편도 그대로 둠(Slide2 와 같음 — 0 으로 자르면 ACADS · Arai & Tagyo + 대조가 3~6 % 벗어남). 대신 저항 합이 0 이하인 면은 `_finish` 가 버림. + """ + p = _pre(sl) + F = np.asarray(F, float) + N, den = _normal(sl, F, lam) + g = 1.0 / F[:, None] + A = sl["sa"] - p["tca"] * g + B0 = p["khws"] - p["cuca"] * g + return N, p["cu"] + N * sl["tanphi"], den, A, B0 -def _ff(sl, F, lam): - N, R, den, *_ = _parts(sl, F, lam) - return (R * sl["ca"]).sum(1) / np.where( - np.abs(_force_den(sl, N)) < _TINY, _TINY, _force_den(sl, N) - ), den +def _ratio(num, den, dnum=None, dden=None, neg=True): + """f = num/den (작은 분모는 ∓TINY) · 도함수 있으면 (f, df).""" + d = np.where(np.abs(den) < _TINY, -_TINY if neg else _TINY, den) + f = num / d + return f if dnum is None else (f, (dnum - f * dden) / d) + + +def _balance(cf, cd, c0, d0, nd, F, neg): + """f = (c0 + Σcf·N)/(d0 + Σcd·N) 와 df/dF · df/dλ (nd = _normal 결과 · 없는 도함수는 None).""" + N, _den, dg, dl = nd + num, den = c0 + (cf * N).sum(1), d0 + (cd * N).sum(1) + if dg is None: + return _ratio(num, den, neg=neg), None, None + f, fg = _ratio(num, den, (cf * dg).sum(1), (cd * dg).sum(1), neg) + fl = None if dl is None else _ratio(num, den, (cf * dl).sum(1), (cd * dl).sum(1), neg)[1] + return f, -fg / np.asarray(F, float) ** 2, fl + + +def _fm(sl, center, F, lam, grad=0, nd=None): + """모멘트 평형 F_m = Σ a·R / D → (F_m, dF_m/dF, dF_m/dλ, 분모).""" + at, gn, d0, ma = _mom(sl, center) + nd = nd or _normal(sl, F, lam, grad) + if len(nd) == 2: + nd = (*nd, None, None) + return (*_balance(at, gn, ma, d0, nd, F, True), nd[1]) + + +def _ff(sl, F, lam, grad=0, nd=None): + """힘 평형 F_f = Σ R·cosα / (Σ(N·sinα + kh·Ws) + Pw) → (F_f, dF_f/dF, dF_f/dλ, 분모).""" + p = _pre(sl) + nd = nd or _normal(sl, F, lam, grad) + if len(nd) == 2: + nd = (*nd, None, None) + return (*_balance(p["tca"], sl["sa"], p["rc0"], p["fd0"], nd, F, False), nd[1]) def _iterate(fun, F0): - """F = fun(F) 풀이(묶음 전체) — 수치 미분 뉴턴(h = fun(F) − F) · 안 되면 반쯤 섞은 고정점 걸음. + """F = fun(F) 풀이(묶음 전체) — 해석 도함수 뉴턴 · 안 되면 반쯤 섞은 고정점 걸음. - 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남. + fun(F) → (f, df/dF, …). 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남. """ F = np.asarray(F0, float).copy() done = np.zeros(len(F), bool) iters = np.zeros(len(F), int) - den = None for _ in range(MAX_ITER): - f1, den = fun(F) + f1, df = fun(F)[:2] h1 = f1 - F - delta = 1e-6 * np.maximum(np.abs(F), 0.1) - f2, _ = fun(F + delta) - dh = ((f2 - (F + delta)) - h1) / delta + dh = df - 1.0 newton = F - h1 / np.where(np.abs(dh) < 1e-12, -1.0, dh) relaxed = 0.5 * (F + f1) use = np.isfinite(newton) & (dh < -1e-9) & (newton > 0.5 * F) & (newton < 2.0 * F) @@ -307,17 +379,20 @@ def _iterate(fun, F0): done |= conv | ~np.isfinite(F) if done.all(): break - f1, den = fun(np.where(np.isfinite(F), F, 1.0)) - return F, iters, done & np.isfinite(F), den + return F, iters, done & np.isfinite(F) def ordinary(sl: dict, center: np.ndarray) -> dict: """Fellenius — 절편간 힘 무시 · 바닥 수직 분해 N = W·cosα − kh·Ws·sinα · 모멘트 평형.""" - N = sl["W"] * sl["ca"] - sl["kh"] * sl["Ws"] * sl["sa"] - R = sl["c"] * sl["l"] + (N - sl["u"] * sl["l"]) * sl["tanphi"] - a, D = _moment_terms(sl, center, N) - F = (a * R).sum(1) / np.where(np.abs(D) < _TINY, -_TINY, D) - ok = sl["valid"] & np.isfinite(F) & (F > F_MIN) & (R.sum(1) > 0) + N = sl["W"] * sl["ca"] - _pre(sl)["khws"] * sl["sa"] + at, gn, d0, ma = _mom(sl, center) + F = _ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1)) + ok = ( + sl["valid"] + & np.isfinite(F) + & (F > F_MIN) + & ((_pre(sl)["cu"] + N * sl["tanphi"]).sum(1) > 0) + ) # 같은 꼴 거르기 — 진출부 급경사(m_α < 0.2) 면은 Bishop 처럼 버림 m = sl["ca"] + sl["sa"] * sl["tanphi"] / np.where(ok, F, 1.0)[:, None] ok &= np.all(m >= MALPHA_MIN, axis=1) @@ -327,7 +402,7 @@ def ordinary(sl: dict, center: np.ndarray) -> dict: def bishop(sl: dict, center: np.ndarray, F0=None) -> dict: F0 = np.full(len(sl["W"]), 1.5) if F0 is None else F0 lam = np.zeros(len(F0)) - F, it, conv, _ = _iterate(lambda f: _fm(sl, center, f, lam), F0) + F, it, conv = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0) return _finish(sl, F, it, conv, lam) @@ -335,7 +410,7 @@ def janbu(sl: dict, F0=None) -> dict: """Janbu 간편(F_0) + 보정 f_0 — 지식DB 02 §2.3 · b1 = 0.31(c=0) · 0.69(φ=0) · 0.50.""" F0 = np.full(len(sl["W"]), 1.5) if F0 is None else F0 lam = np.zeros(len(F0)) - F, it, conv, _ = _iterate(lambda f: _ff(sl, f, lam), F0) + F, it, conv = _iterate(lambda f: _ff(sl, f, lam, 1), F0) out = _finish(sl, F, it, conv, lam) out["F0"] = out["F"] out["f0"] = janbu_f0(sl) @@ -350,7 +425,8 @@ def janbu_f0(sl: dict) -> np.ndarray: dist = np.abs(dz * (xe - xe[:, :1]) - dx * (ze - ze[:, :1]))[:, :] / np.maximum( L[:, None], _TINY ) - ratio = dist.max(1) / np.maximum(L, _TINY) + # d/L ≤ 0.35 에서만 식이 맞음(지식DB 02 §2.3) — 넘으면 0.35 값(f0 곡선 꼭대기 근처) + ratio = np.minimum(dist.max(1) / np.maximum(L, _TINY), JANBU_DL_MAX) c_all0 = np.all(sl["c"] <= 0, axis=1) phi_all0 = np.all(sl["tanphi"] <= 0, axis=1) b1 = np.where(c_all0, 0.31, np.where(phi_all0, 0.69, 0.50)) @@ -358,7 +434,7 @@ def janbu_f0(sl: dict) -> np.ndarray: def spencer(sl: dict, center: np.ndarray, F0=None, scan: bool = True) -> dict: - """Spencer — 모멘트 · 힘 두 평형을 함께 맞추는 (F, λ) · 뉴턴(수치 야코비안) · 안 맞으면 λ 훑기. + """Spencer — 모멘트 · 힘 두 평형을 함께 맞추는 (F, λ) · 뉴턴(해석 야코비안) · 안 맞으면 λ 훑기. scan=False — 탐색 고리 안(시험 면 대부분이 버려짐)에서는 훑기를 건너뜀(시간 대부분이 여기). """ @@ -381,27 +457,27 @@ def spencer(sl: dict, center: np.ndarray, F0=None, scan: bool = True) -> dict: return out -def _residual(sl, center, F, lam): - N, R, den, *_ = _parts(sl, F, lam) - a, D = _moment_terms(sl, center, N) - fd = _force_den(sl, N) - rm = ((a * R).sum(1) - F * D) / np.maximum(np.abs(D), _TINY) - rf = ((R * sl["ca"]).sum(1) - F * fd) / np.maximum(np.abs(fd), _TINY) - return rm, rf - - def _spencer_newton(sl, center, F, lam, iters=40): + """(F_m − F, F_f − F) = 0 뉴턴 — 걸음 크기 제한(F 절반 · λ 0.3) · 수렴 = 두 잔차 < 1e-4 · λ 범위 안. + + 끝난 면은 빼고(남은 면이 절반 아래면 묶음을 줄임) · F < F_MIN 으로 빠진 면은 버림(어차피 `_finish` 가 버림). + """ F, lam = F.copy(), lam.copy() it = np.zeros(len(F), int) done = np.zeros(len(F), bool) + idx, w, wc = np.arange(len(F)), sl, center for _ in range(iters): - rm, rf = _residual(sl, center, F, lam) - hF = 1e-6 * np.maximum(F, 0.1) - hl = 1e-6 - rmF, rfF = _residual(sl, center, F + hF, lam) - rml, rfl = _residual(sl, center, F, lam + hl) - j11, j12 = (rmF - rm) / hF, (rml - rm) / hl - j21, j22 = (rfF - rf) / hF, (rfl - rf) / hl + act = ~done[idx] + if not act.any(): + break + if act.sum() < 0.5 * len(idx): + idx, w, wc = idx[act], _subset(w, act), wc[act] + f, lm, dn = F[idx], lam[idx], done[idx] + nd = _normal(w, f, lm, 2) + fm, fmF, fml, _ = _fm(w, wc, f, lm, nd=nd) + ff, ffF, ffl, _ = _ff(w, f, lm, nd=nd) + rm, rf = fm - f, ff - f + j11, j12, j21, j22 = fmF - 1.0, fml, ffF - 1.0, ffl det = j11 * j22 - j12 * j21 det = np.where(np.abs(det) < 1e-14, np.nan, det) dF = -(j22 * rm - j12 * rf) / det @@ -409,34 +485,41 @@ def _spencer_newton(sl, center, F, lam, iters=40): scale = np.minimum( 1.0, np.minimum( - 0.5 * F / np.maximum(np.abs(dF), _TINY), 0.3 / np.maximum(np.abs(dl), _TINY) + 0.5 * f / np.maximum(np.abs(dF), _TINY), 0.3 / np.maximum(np.abs(dl), _TINY) ), ) - F = np.where(done, F, F + scale * dF) - lam = np.where(done, lam, lam + scale * dl) - it += ~done - small = (np.abs(dF) < TOL * np.maximum(1, np.abs(F))) & (np.abs(dl) < TOL) - done |= small | ~np.isfinite(F) - if done.all(): - break - rm, rf = _residual(sl, center, F, lam) + f = np.where(dn, f, f + scale * dF) + F[idx], lam[idx] = f, np.where(dn, lm, lm + scale * dl) + it[idx] += ~dn + small = (np.abs(dF) < TOL * np.maximum(1, np.abs(f))) & (np.abs(dl) < TOL) + done[idx] = dn | small | ~np.isfinite(f) | (f < F_MIN) + nd = _normal(sl, F, lam) + rm = _fm(sl, center, F, lam, nd=nd)[0] - F + rf = _ff(sl, F, lam, nd=nd)[0] - F conv = done & np.isfinite(F) & (F > 0) & (np.abs(rm) < 1e-4) & (np.abs(rf) < 1e-4) conv &= (lam >= LAMBDA_RANGE[0]) & (lam <= LAMBDA_RANGE[1]) return F, lam, it, conv -def _spencer_scan(sl, center): - """λ 를 훑어 F_m(λ) − F_f(λ) 부호가 바뀌는 곳(λ=0 가까운 쪽)을 출발점으로.""" - grid = np.linspace(LAMBDA_RANGE[0], LAMBDA_RANGE[1], 16) +def lambda_curve(sl, center, grid=None): + """λ 마다 모멘트 · 힘 평형 F_m(λ) · F_f(λ) — 꼴 (len(grid), S) 둘(화면 F–λ 그래프 · 훑기).""" + grid = np.linspace(LAMBDA_RANGE[0], LAMBDA_RANGE[1], 16) if grid is None else grid S = len(sl["W"]) - diff = np.full((len(grid), S), np.nan) fm_all = np.full((len(grid), S), np.nan) + ff_all = np.full((len(grid), S), np.nan) F0 = np.full(S, 1.5) for k, g in enumerate(grid): lam = np.full(S, g) - fm, *_ = _iterate(lambda f: _fm(sl, center, f, lam), F0) - ff, *_ = _iterate(lambda f: _ff(sl, f, lam), F0) - diff[k], fm_all[k] = fm - ff, fm + fm_all[k] = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0)[0] + ff_all[k] = _iterate(lambda f: _ff(sl, f, lam, 1), F0)[0] + return grid, fm_all, ff_all + + +def _spencer_scan(sl, center): + """λ 를 훑어 F_m(λ) − F_f(λ) 부호가 바뀌는 곳(λ=0 가까운 쪽)을 출발점으로.""" + grid, fm_all, ff_all = lambda_curve(sl, center) + S = len(sl["W"]) + diff = fm_all - ff_all best_F, best_l = np.full(S, 1.5), np.zeros(S) best_gap = np.full(S, np.inf) for k in range(len(grid) - 1): @@ -453,15 +536,18 @@ def _spencer_scan(sl, center): def _subset(sl: dict, mask: np.ndarray) -> dict: - out = { + """면 묶음 일부 — 미리 묶은 칸(_pre · _mom)은 버림(부분 꼴이 달라 다시 셈).""" + return { k: (v[mask] if isinstance(v, np.ndarray) and v.shape[:1] == mask.shape else v) for k, v in sl.items() + if k not in ("_pre", "_mom") } - return out def _finish(sl, F, it, conv, lam): - _, R, den, *_ = _parts(sl, np.where(np.isfinite(F) & (F > 0), F, 1.0), lam) + Fs = np.where(np.isfinite(F) & (F > 0), F, 1.0) + N, den = _normal(sl, Fs, lam) + R = _pre(sl)["cu"] + N * sl["tanphi"] malpha_ok = np.all(den >= MALPHA_MIN, axis=1) ok = sl["valid"] & conv & np.isfinite(F) & (F > F_MIN) & malpha_ok & (R.sum(1) > 0) return { @@ -472,19 +558,33 @@ def _finish(sl, F, it, conv, lam): } -def slice_table(sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bishop") -> list[dict]: - """한 면(k 번째)의 절편표 — 폭 · 높이 · 경사 · 무게 · 수압 · 강도 · N' · 동원 전단 · (Spencer) E · X.""" +def slice_table( + sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bishop", f=None +) -> list[dict]: + """한 면(k 번째)의 절편표 + 자유물체도 힘 — 행 = 머리 → 끝 · left = 머리 쪽 경계. + + 폭 · 높이 · 경사 · 무게 · 수압 · 강도 · N(전) · N' · 동원 전단 · 저항 · 지진력 · 경계 E · X. + M-P 는 f(경계마다 절편간 힘 함수)를 받아 차례 계산(`Engine_GLE.forces`). + """ one = _subset(sl, np.arange(len(sl["W"])) == k) if method == "ordinary": N = one["W"] * one["ca"] - one["kh"] * one["Ws"] * one["sa"] E = np.zeros((1, N.shape[1] + 1)) + X = E + elif method == "morgenstern_price": + from E03_StabSection.E03_StabSection_Engine_GLE import forces + + N, _den, E, X = forces(one, np.array([F]), np.array([lam]), f[k : k + 1]) else: N, R, _den, A, B0 = _parts(one, np.array([F]), np.array([lam])) dE = B0 + N * A E = np.concatenate([one["Pw"][:, None], one["Pw"][:, None] + np.cumsum(dE, axis=1)], axis=1) + X = lam * E + X[:, 0] = 0.0 R = one["c"] * one["l"] + (N - one["u"] * one["l"]) * one["tanphi"] rows = [] for i in range(N.shape[1]): + ul = float(one["u"][0, i] * one["l"][0, i]) rows.append( { "x_left": float(one["xe"][0, i]), @@ -499,10 +599,20 @@ def slice_table(sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bis "material": int(one["mat"][0, i]), "c": float(one["c"][0, i]), "phi_deg": float(np.degrees(np.arctan(one["tanphi"][0, i]))), - "normal_eff": float(N[0, i] - one["u"][0, i] * one["l"][0, i]), + "normal_eff": float(N[0, i]) - ul, "shear_mobilized": float(R[0, i] / F), "E_left": float(E[0, i]), - "X_left": float(lam * E[0, i]), + "X_left": float(X[0, i]), + # 덧(2차) — 자유물체도 + "x_mid": float(one["xm"][0, i]), + "z_base": float(one["zb"][0, i]), + "z_top": float(one["zt"][0, i]), + "normal_total": float(N[0, i]), + "pore_force": ul, + "shear_strength": float(R[0, i]), + "seismic_force": float(one["kh"] * one["Ws"][0, i]), + "E_right": float(E[0, i + 1]), + "X_right": float(X[0, i + 1]), } ) return rows diff --git a/E03_StabSection/E03_StabSection_Engine_Search.py b/E03_StabSection/E03_StabSection_Engine_Search.py index 3ce607bb3..014af338c 100644 --- a/E03_StabSection/E03_StabSection_Engine_Search.py +++ b/E03_StabSection/E03_StabSection_Engine_Search.py @@ -2,8 +2,8 @@ 활동 방향 +x 로 맞춘 `Prepared` 위에서 찾고 좌표를 원래 방향으로 되돌려 냄. 다듬기 · 비원호 = 차분진화(best/1/bin · 무리 전체를 한 번에 평가 · seed 고정 → 같은 입력 같은 답) -+ 비원호 마지막 국부 다듬기 = SciPy Nelder-Mead. SciPy differential_evolution 은 개체마다 파이썬 돌연변이 고리라 -한 비탈 시간의 30 % 를 먹어 직접 묶음 꼴로 씀. ++ 비원호 마지막 국부 다듬기 = 묶음 Nelder-Mead(후보 넷을 한 번에). SciPy 최적화는 점마다 파이썬 고리라 느려 +직접 묶음 꼴로 씀. 방법마다 찾은 면을 모든 방법으로 다시 셈(교차) → 한 방법 탐색이 놓친 골짜기를 메움. 상용 SW 대응: Grid & Radius(SLOPE/W · Slide) · Entry & Exit(SLOPE/W) · Auto Refine(Slide) · Block Search + Optimize(Slide · SLOPE/W 꼭짓점 최적화). """ @@ -11,181 +11,38 @@ Block Search + Optimize(Slide · SLOPE/W 꼭짓점 최적화). import time import numpy as np -from scipy.optimize import minimize from E03_StabSection import E03_StabSection_Engine_LEM as lem +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, +) -METHODS = ("ordinary", "bishop", "janbu", "spencer") _BIG = 1e3 -_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) -> dict: - 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" in methods or "spencer" in 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) - 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 +# 탐색 예산(실험 손잡이) +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 # ---- 원호 --------------------------------------------------------------------- -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 - - def entry_exit_circles(prep, limits, n_points=25, thetas_deg=(10, 20, 35, 50, 65, 80, 100, 120)): xl, xr = limits pts = np.linspace(xl, xr, n_points) @@ -237,9 +94,9 @@ def _circle_eval(prep, xc, zc, R, xE, xX, n, min_depth, methods, chunk=4000): return res -def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=60): +def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=None): """원호 다듬기 — (진입 x, 진출 x, 중심각) 차분진화 · 첫 무리 = 거친 탐색의 좋은 면.""" - fun = _objective(method) + fun = objective(method) lo = np.array([limits[0], limits[0], np.radians(3)]) hi = np.array([limits[1], limits[1], np.radians(170)]) @@ -252,9 +109,35 @@ def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=60) F[ok] = np.where(np.isfinite(f), f, _BIG) return F - x, fmin, _pop, _f = _de(obj, lo, hi, _population(seeds, lo, hi, 30, seed), maxiter, seed) + 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]) - return float(fmin), (float(xc[0]), float(zc[0]), float(R[0]), float(x[0]), float(x[1])) + 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], limits[0], 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): @@ -310,6 +193,66 @@ def _local_search(obj, starts, lo, hi, step, maxiter, seed, lam=24, starts_keep= 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) @@ -352,10 +295,10 @@ def layer_bottom_depth(prep, x, material): def optimize_polyline( - prep, limits, method, seeds, n, min_depth, m=4, seed=0, maxiter=60, polish=True + prep, limits, method, seeds, n, min_depth, m=4, seed=0, maxiter=None, polish=True ): """비원호 꼭짓점 최적화 — 차분진화 뒤 가장 좋은 면을 Nelder-Mead 로 국부 다듬기(상용 SAO · Optimize 대응).""" - fun = _objective(method) + fun = objective(method) xl, xr = limits depth_max = max(prep.surface(np.linspace(xl, xr, 50)).max() - prep.bottom_z, 1.0) lo = np.array([xl, xl] + [0.02] * m + [0.0] * m) @@ -373,16 +316,15 @@ def optimize_polyline( 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, seed) - if polish and fmin < _BIG: - r = minimize( - lambda q: float(obj(np.clip(q, lo, hi))[0]), - x, - method="Nelder-Mead", - options={"maxfev": 200, "xatol": 1e-4, "fatol": 1e-6}, - ) - if r.fun < fmin: - x, fmin = np.clip(r.x, lo, hi), float(r.fun) + 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]] @@ -452,53 +394,6 @@ def block_seeds(prep, limits, m, n_points=12, angles=((55, 35), (65, 25))): return out -# ---- 결과 꾸미기 ---------------------------------------------------------------- -def describe(prep, kind, geom, method, n, min_depth) -> dict: - s = prep.direction - if kind == "circle": - xc, zc, R, xE, xX = geom - sl = circle_slices(prep, [xc], [zc], [R], [xE], [xX], n, min_depth) - center = np.array([[xc, zc]]) - xs = sl["xe"][0] - pts = np.column_stack([xs, sl["ze"][0]]) - else: - vx, vz = geom - sl = poly_slices(prep, vx[None, :], vz[None, :], n, min_depth) - center = three_point_center(vx[None, :], vz[None, :]) - pts = np.column_stack([vx, vz]) - res = evaluate(sl, center, (method,))[method] - F = float(res["F"][0]) - lam = float(res.get("lambda", [0.0])[0]) if method == "spencer" else 0.0 - rows = lem.slice_table(sl, 0, F, lam, method) if np.isfinite(F) else [] - for r in rows: - r["material"] = ( - prep.materials[r["material"]]["id"] if r["material"] < len(prep.materials) else None - ) - r["x_left"], r["x_right"] = sorted((s * r["x_left"], s * r["x_right"])) - out = { - "F": F if np.isfinite(F) else None, - "kind": kind, - "points": [[float(s * x), float(z)] for x, z in (pts if s > 0 else pts[::-1])], - "entry": [float(s * sl["xe"][0, 0]), float(sl["ze"][0, 0])], - "exit": [float(s * sl["xe"][0, -1]), float(sl["ze"][0, -1])], - "iterations": int(res["iterations"][0]), - "converged": bool(res["converged"][0]), - "slices": rows, - } - if kind == "circle": - out["center"] = [float(s * center[0, 0]), float(center[0, 1])] - out["radius"] = float(geom[2]) - if method == "spencer": - out["lambda"] = lam - out["theta_deg"] = float(np.degrees(np.arctan(lam))) - if method == "janbu": - out["F0"] = float(res["F0"][0]) if np.isfinite(res["F0"][0]) else None - out["f0"] = float(res["f0"][0]) - if prep.crack_depth > 0: - out["crack_water_force"] = float(sl["Pw"][0]) - return out - - # ---- 한 비탈 해석 ---------------------------------------------------------------- def analyze( model: dict, @@ -526,14 +421,15 @@ def analyze( lim = tuple(sorted((s * float(limits[0]), s * float(limits[1])))) results: dict = {m: {} for m in methods} counts = {} - contour = None + 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 != "spencer") or ("bishop",) + 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]]))) @@ -553,17 +449,44 @@ def analyze( for row in cmin.reshape(len(gz), len(gx)) ], } - for m in methods: + 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)]) - fmin, geom = refine_circle(prep, lim, m, seeds, n_slices, min_depth, seed) + 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: - results[m]["circular"] = describe(prep, "circle", geom, m, n_slices, min_depth) 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) if noncircular: mv = polyline_vertices seeds = [] @@ -576,35 +499,73 @@ def analyze( 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) + # 힘 평형 방법(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 + prep, lim, m, seeds, n_slices, min_depth, mv, seed, iters ) seeds = [*top, *seeds] if fmin < _BIG: - results[m]["noncircular"] = describe(prep, "polyline", geom, m, n_slices, min_depth) found.append(geom) - for m in [k for k in ("ordinary", "bishop") if k in methods]: - cands = [describe(prep, "polyline", g, m, n_slices, min_depth) for g in found] - cands = [c for c in cands if c["F"] is not None] - if cands: - results[m]["noncircular"] = min(cands, key=lambda c: c["F"]) + 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: - cands = [v for v in results[m].values() if v.get("F") is not None] + # 비원호 면 위 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 return { "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), } +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)) @@ -630,5 +591,5 @@ def evaluate_specified( 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) + out[m] = _describe(prep, kind, geom, m, n_slices, 0.0) return out diff --git a/E03_StabSection/E03_StabSection_Engine_Station.py b/E03_StabSection/E03_StabSection_Engine_Station.py index e6f9438dd..dbfdae283 100644 --- a/E03_StabSection/E03_StabSection_Engine_Station.py +++ b/E03_StabSection/E03_StabSection_Engine_Station.py @@ -25,10 +25,16 @@ def analyze_station( methods=METHODS, noncircular: bool = True, min_depth: float = MIN_DEPTH_M, + case_workers: int = 1, **model_kw, ) -> dict: + """측점 하나 — case_workers > 1 이면 비탈 × 조건 경우를 프로세스로 나눠 돎(측점 하나만 다시 셀 때). + + 서버 프로세스 안에서는 1(차례로) 그대로 둠 — 자식이 main.py 를 다시 읽음(`Worker` 머리글). + """ t0 = time.perf_counter() out = {"chainage_m": cross.get("chainage_m"), "label": cross.get("label"), "cases": []} + jobs = [] for cond in conditions: model = section_model_from_b06( cross, @@ -41,17 +47,24 @@ def analyze_station( **model_kw, ) out.setdefault("areas", model["areas"]) - for side in model["sides"]: - res = analyze( - model, - direction=side["direction"], - limits=side["limits"], - methods=methods, - noncircular=noncircular, - min_depth=min_depth, - ) - mins = {m: (res["methods"][m]["min"] or {}).get("F") for m in methods} - out["cases"].append({"condition": cond, **side, "min_F": mins, "result": res}) + jobs += [(cond, model, side) for side in model["sides"]] + kw = {"methods": methods, "noncircular": noncircular, "min_depth": min_depth} + if case_workers > 1 and len(jobs) > 1: + from concurrent.futures import ProcessPoolExecutor + + with ProcessPoolExecutor(max_workers=min(case_workers, len(jobs))) as pool: + futures = [ + pool.submit(analyze, m, direction=sd["direction"], limits=sd["limits"], **kw) + for _, m, sd in jobs + ] + results = [f.result() for f in futures] + else: + results = [ + analyze(m, direction=sd["direction"], limits=sd["limits"], **kw) for _, m, sd in jobs + ] + for (cond, _m, side), res in zip(jobs, results): + mins = {m: (res["methods"][m]["min"] or {}).get("F") for m in methods} + out["cases"].append({"condition": cond, **side, "min_F": mins, "result": res}) out["elapsed_s"] = round(time.perf_counter() - t0, 3) return out diff --git a/E03_StabSection/E03_StabSection_Engine_Surface.py b/E03_StabSection/E03_StabSection_Engine_Surface.py new file mode 100644 index 000000000..ff78a96d1 --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_Surface.py @@ -0,0 +1,185 @@ +"""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 diff --git a/resources/tester/test_532_2_stab_lem_benchmarks.py b/resources/tester/test_532_2_stab_lem_benchmarks.py index bc5da1f29..676d0a48f 100644 --- a/resources/tester/test_532_2_stab_lem_benchmarks.py +++ b/resources/tester/test_532_2_stab_lem_benchmarks.py @@ -180,9 +180,8 @@ def test_acads_3a_weak_seam_noncircular_search(): near(nonc(r, "spencer"), 1.258) -def test_acads_4_inclined_seam_water_surcharge(): +def _acads4_model(): # S2VM #9 (ACADS 4) — 기운 약층 · 지하수위선 · 등분포 20 kPa + 선형 20→40 kPa - # 최적화 비원호 Spencer 0.707 · Janbu 보정 0.699 (약층 좌표 = 저해상 그림 판독 — 출처 ±5 % 권함) ext = Polygon([[20, 15], [84, 15], [84, 40], [67.5, 40], [43, 27.75], [20, 27.75]]) seam = Polygon([[20, 18.26], [84, 36.2], [84, 36.8], [20, 18.88]]) regions = [{"material_id": "seam", "polygon": list(ext.intersection(seam).exterior.coords)}] + [ @@ -208,7 +207,12 @@ def test_acads_4_inclined_seam_water_surcharge(): ], }, } - r = analyze(m, direction=-1, methods=("janbu", "spencer")) + return m + + +def test_acads_4_inclined_seam_water_surcharge(): + # S2VM #9 — 최적화 비원호 Spencer 0.707 · Janbu 보정 0.699 (약층 좌표 = 저해상 그림 판독 — 출처 ±5 % 권함) + r = analyze(_acads4_model(), direction=-1, methods=("janbu", "spencer")) near(nonc(r, "spencer"), 0.707) near(nonc(r, "janbu"), 0.699) diff --git a/resources/tester/test_532_2_stab_mp_outputs.py b/resources/tester/test_532_2_stab_mp_outputs.py new file mode 100644 index 000000000..71d592357 --- /dev/null +++ b/resources/tester/test_532_2_stab_mp_outputs.py @@ -0,0 +1,141 @@ +# -*- coding: utf-8 -*- +"""532-2 2차 E03 2D — Morgenstern-Price(반정현) 대조 · 결과 자료(자유물체도 · F–λ 곡선 · Safety Map). + +기대값 = Rocscience Slide2 Slope Stability Verification Manual(지식DB 참고문헌) 원문 GLE · M-P 값. +- #1 ACADS 1a GLE 0.986 · #3 ACADS 1c GLE 1.374 · #8 ACADS 3b 지정 비원호 GLE 1.262 · + #9 ACADS 4 블록 + 최적화 GLE 0.683 · #21 F&K 지정 원 M-P 건기 2.075 · r_u 1.760(Slide 열). +- 뺀 것: #7 · #15 — 원문 GLE(1.275 · 0.437)가 같은 표 Spencer(1.246 · 0.409)보다 높음(Slide 탐색이 GLE + 임계면을 놓친 값 · 같은 문제 #8 지정면 GLE 는 Spencer 보다 낮음) → 대조 기준으로 못 씀. + + ./venv/Scripts/python.exe -m pytest resources/tester/test_532_2_stab_mp_outputs.py -q +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import numpy as np +import pytest + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from E03_StabSection import E03_StabSection_Engine_GLE as gle # noqa: E402 +from E03_StabSection import E03_StabSection_Engine_LEM as lem # noqa: E402 +from E03_StabSection.E03_StabSection_Engine_Search import ( # noqa: E402 + analyze, + evaluate_specified, +) +from E03_StabSection.E03_StabSection_Engine_Surface import ( # noqa: E402 + poly_slices, + three_point_center, +) +from test_532_2_stab_lem_benchmarks import ( # noqa: E402 + ACADS1C_MATS, + ACADS1C_REGIONS, + ACADS_EXT, + FK_EXT, + _acads3_model, + _acads4_model, + circ, + mats, + near, + nonc, +) + +MP = "morgenstern_price" +SEAM_3B = [[41.85, 27.75], [44.0, 26.5], [63.5, 27.0], [73.31, 40.0]] + + +def test_mp_constant_function_equals_spencer(): + # f(x) ≡ 1 이면 M-P 식 = Spencer 식 — 같은 F · λ + prep = lem.Prepared(_acads3_model(), -1) + pts = sorted((-x, z) for x, z in SEAM_3B) + vx, vz = np.array([[p[0] for p in pts]]), np.array([[p[1] for p in pts]]) + sl = poly_slices(prep, vx, vz, 50) + center = three_point_center(vx, vz) + sp = lem.spencer(sl, center, np.array([1.3])) + mp = gle.morgenstern_price(sl, center, "constant", F0=np.array([1.3])) + assert abs(mp["F"][0] - sp["F"][0]) < 1e-6 + assert abs(mp["lambda"][0] - sp["lambda"][0]) < 1e-5 + + +def test_mp_specified_surfaces_slide2(): + # #8 지정 비원호 GLE 1.262 · #21 지정 원 M-P 건기 2.075 · r_u 0.25 1.760 + r = evaluate_specified(_acads3_model(), polyline=SEAM_3B, direction=-1, methods=(MP,)) + near(r[MP]["F"], 1.262, 0.01) + m = {"materials": mats(("s", 120, 600, 20)), "regions": [{"material_id": "s", "polygon": FK_EXT}]} + near(evaluate_specified(m, circle=[120, 90, 80], direction=1, methods=(MP,))[MP]["F"], 2.075, 0.01) + m["groundwater"] = {"mode": "ru", "ru": 0.25} + near(evaluate_specified(m, circle=[120, 90, 80], direction=1, methods=(MP,))[MP]["F"], 1.760, 0.01) + + +def test_mp_search_slide2(): + # #1 원호 GLE 0.986 · #3 원호 GLE 1.374 + m = {"materials": mats(("s", 20, 3, 19.6)), "regions": [{"material_id": "s", "polygon": ACADS_EXT}]} + near(circ(analyze(m, direction=-1, noncircular=False), MP), 0.986, 0.01) + r = analyze({"materials": ACADS1C_MATS, "regions": ACADS1C_REGIONS}, direction=-1, noncircular=False) + near(circ(r, MP), 1.374, 0.01) + + +def test_mp_acads4_noncircular(): + # #9 블록 + 최적화 GLE 0.683 (약층 좌표 = 그림 판독 → ±2 %) + r = analyze(_acads4_model(), direction=-1, methods=("janbu", "spencer", MP)) + near(nonc(r, MP), 0.683, 0.02) + + +def _closure(rows, lam, f_last=1.0): + weight = sum(r["weight"] for r in rows) + for a, b in zip(rows, rows[1:]): + assert abs(a["E_right"] - b["E_left"]) < 1e-9 and abs(a["X_right"] - b["X_left"]) < 1e-9 + assert abs(rows[-1]["E_right"]) < 1e-3 * weight # 끝 절편에서 힘 평형이 닫힘 + assert abs(rows[0]["X_left"]) < 1e-12 + for r in rows: + # 연직 평형 — N·cosα + S·sinα = W + (X_right − X_left) (지진 · 균열 없음) + a = np.radians(r["alpha_deg"]) + lhs = r["normal_total"] * np.cos(a) + r["shear_mobilized"] * np.sin(a) + assert abs(lhs - (r["weight"] + r["X_right"] - r["X_left"])) < 1e-6 * weight + + +def test_free_body_and_lambda_curve_outputs(): + m = {"materials": mats(("s", 20, 3, 19.6)), "regions": [{"material_id": "s", "polygon": ACADS_EXT}]} + r = evaluate_specified(m, circle=[29.104, 55.527, 30.527], direction=-1, methods=("spencer", MP)) + for key in ("spencer", MP): + out = r[key] + rows = out["slices"] + assert {"normal_total", "pore_force", "shear_strength", "E_right", "X_right", "x_mid"} <= rows[0].keys() + _closure(rows, out["lambda"]) + curve = out["lambda_curve"] + lam, fm, ff = (np.array(curve[k], float) for k in ("lambda", "Fm", "Ff")) + assert len(lam) == gle.CURVE_POINTS and np.isfinite(fm).sum() > 10 + # 곡선 교점 = 해 (F · λ) + d = fm - ff + k = np.where(np.isfinite(d[:-1]) & np.isfinite(d[1:]) & (np.sign(d[:-1]) != np.sign(d[1:])))[0] + assert len(k) + j = k[np.argmin(np.abs(lam[k] - out["lambda"]))] + t = d[j] / (d[j] - d[j + 1]) + assert abs(lam[j] + t * (lam[j + 1] - lam[j]) - out["lambda"]) < 0.02 + assert abs(fm[j] + t * (fm[j + 1] - fm[j]) - out["F"]) < 0.01 * out["F"] + assert abs(r[MP]["slices"][-1]["X_right"]) < 1e-9 # 반정현 f = 0 (끝) + + +def test_safety_map_and_reference_only(): + m = {"materials": mats(("s", 17.64, 9.8, 10)), "regions": [{"material_id": "s", "polygon": [[0, 0], [25, 0], [25, 10], [15, 10], [5, 5], [0, 5]]}]} + r = analyze(m, direction=-1) + smap = r["safety_map"]["surfaces"] + assert 0 < len(smap) <= 50 + F = [s["F"] for s in smap] + assert F == sorted(F) and F[0] >= 0.97 * circ(r, "bishop") + assert all(len(s["points"]) == 25 and s["points"][0][0] < s["points"][-1][0] for s in smap) + # 비원호 Fellenius · Bishop = 참고 칸 · 최소 F 는 원호에서만 + for k in ("ordinary", "bishop"): + res = r["methods"][k] + assert res["noncircular"]["reference_only"] is True + assert res["min"]["kind"] == "circle" + assert r["methods"]["spencer"]["min"]["F"] <= circ(r, "spencer") + 1e-12 + + +if __name__ == "__main__": + sys.exit(pytest.main([__file__, "-q"]))