From bfd87f3e3eec64595adafad7fe4261348179aaef Mon Sep 17 00:00:00 2001 From: umsangdon Date: Sun, 25 Oct 2026 20:07:14 +0900 Subject: [PATCH] =?UTF-8?q?feat(E03):=202D=20=EC=A7=80=EC=B8=B5=20?= =?UTF-8?q?=EB=AA=A8=EB=8D=B8=20=EB=B9=8C=EB=8D=94=20=C2=B7=20=ED=95=9C?= =?UTF-8?q?=EA=B3=84=ED=8F=89=ED=98=95=20=EC=86=94=EB=B2=84(Fellenius=20?= =?UTF-8?q?=C2=B7=20Bishop=20=C2=B7=20Janbu=20=EB=B3=B4=EC=A0=95=20=C2=B7?= =?UTF-8?q?=20Spencer)=20=C2=B7=20=EC=9B=90=ED=98=B8=20=C2=B7=20=EB=B9=84?= =?UTF-8?q?=EC=9B=90=ED=98=B8=20=ED=83=90=EC=83=89=20(532-2)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 모델: B06 횡단(지반선 · 설계선 · 암 경계 · 성토체 · 벽) → 지층 폴리곤 · 지하수 · 노면 하중 · 비탈별 방향 · 탐색 한계 - 솔버: 한 벌 절편 식(λ = 0 Bishop · Janbu / Spencer 모멘트 · 힘 동시) · 간극수압 none · r_u · 지하수위선 · kh · 상재하중 · 인장균열 - 탐색: Entry & Exit · Grid & Radius · 차분진화 다듬기 · 블록 씨앗 + 꼭짓점 국부 최적화 + Nelder-Mead - 출력: 방법별 최소 안전율 · 임계면 좌표 · 절편표 · 등안전율 격자 · 수렴 정보 - 시험: Slide2 검증 예제 12 문제(ACADS 1a~4 · F&K · Loukidis · Arai & Tagyo · Yamagami & Ueta) ±3 % · 면적 = B06 Co-Authored-By: Claude Opus 5.5 --- E03_StabSection/E03_StabSection_Engine_LEM.py | 508 ++++++++++++++ .../E03_StabSection_Engine_Model.py | 293 ++++++++ .../E03_StabSection_Engine_Search.py | 634 ++++++++++++++++++ .../E03_StabSection_Engine_Station.py | 75 +++ .../tester/test_532_2_stab_lem_benchmarks.py | 348 ++++++++++ .../tester/test_532_2_stab_section_model.py | 155 +++++ 6 files changed, 2013 insertions(+) create mode 100644 E03_StabSection/E03_StabSection_Engine_LEM.py create mode 100644 E03_StabSection/E03_StabSection_Engine_Model.py create mode 100644 E03_StabSection/E03_StabSection_Engine_Search.py create mode 100644 E03_StabSection/E03_StabSection_Engine_Station.py create mode 100644 resources/tester/test_532_2_stab_lem_benchmarks.py create mode 100644 resources/tester/test_532_2_stab_section_model.py diff --git a/E03_StabSection/E03_StabSection_Engine_LEM.py b/E03_StabSection/E03_StabSection_Engine_LEM.py new file mode 100644 index 000000000..b68ec84ba --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_LEM.py @@ -0,0 +1,508 @@ +"""E03 2D 한계평형해석(LEM) — 절편 · 네 방법(Fellenius · Bishop 간편 · Janbu 간편+보정 · Spencer). + +좌표 x 오른쪽 · z 위 · 활동 방향 +x(왼쪽 위 → 오른쪽 아래). 반대쪽 비탈은 `Prepared(direction=-1)` +이 모형을 좌우로 뒤집어 맞춤 — 결과 좌표를 되돌리는 것은 부르는 쪽 몫(`unflip`). +식 = 지식DB `original/안정성평가/02` 절편 힘 · 모멘트 평형. 네 방법 모두 한 벌 식에서 나옴: + 절편 i 연직 평형 → N = [W + λ·B0 − (c·l − u·l·tanφ)·sinα/F] / (m_α − λ·A) + m_α = cosα + sinα·tanφ/F · A = sinα − tanφ·cosα/F · B0 = kh·Ws − (c·l − u·l·tanφ)·cosα/F + λ = 절편간 힘 기울기(X = λ·E) — Bishop · Janbu 는 λ = 0 · Spencer 는 모멘트 · 힘 평형을 함께 맞추는 λ. + 모멘트 F_m = Σ a·R / D (점 O 둘레) · 힘 F_f = Σ R·cosα / (Σ(N·sinα + kh·Ws) + Pw) · R = c·l + (N − u·l)·tanφ +외부 코드 없음(NumPy 만) — 상용 · 오픈소스 결과는 대조 시험에서만 씀. +""" + +import numpy as np + +from E03_StabSection.E03_StabSection_Types import GAMMA_W, groundwater, material, surcharges + +_TINY = 1e-9 +_PAD_Z = -1e9 +#: m_α 가 이 값보다 작은 절편이 있으면 그 활동면을 버림(상용 SW 관례 0.2 — 진출부 급경사 발산 막음). +MALPHA_MIN = 0.2 +#: Spencer 절편간 힘 기울기 λ = X/E 범위 — 여기 부호(X = 윗 절편이 아랫 절편에 주는 연직력, 위 +)로 +#: 윗 덩이가 아랫 덩이를 내리누르면 λ < 0. 아래 한도 −1.25 = SLOPE/W GLE λ 범위 · 위는 작은 + 만 허용 +#: (아랫 덩이를 들어 올리는 해 = 비현실 근 — Yamagami & Ueta 대조에서 λ ≈ +1 근이 F 를 9 % 낮췄음). +LAMBDA_RANGE = (-1.25, 0.25) +#: 이보다 작은 F = 1/F 항이 터진 엉뚱한 근(실제 비탈 F 는 이 위) — 버림. +F_MIN = 0.05 +TOL = 1e-6 +MAX_ITER = 100 + + +class Prepared: + """모형 → 절편 계산용 배열. 지층 = 꼭짓점 x 사이 띠마다 연직 구간(바닥 · 윗면 · 물성) — 띠 안은 선형이라 정확.""" + + def __init__(self, model: dict, direction: int = 1): + self.direction = 1 if direction >= 0 else -1 + s = float(self.direction) + mats = [material(m) for m in model["materials"]] + index = {m["id"]: i for i, m in enumerate(mats)} + self.materials = mats + pad = len(mats) # 빈 칸 번호 + self.gamma_t = np.array([m["gamma_t"] for m in mats] + [0.0]) + self.gamma_sat = np.array([m["gamma_sat"] for m in mats] + [0.0]) + self.c = np.array([m["c"] for m in mats] + [np.nan]) + self.tanphi = np.array([np.tan(np.radians(m["phi"])) for m in mats] + [np.nan]) + polys = [] + for region in model["regions"]: + rings = [region["polygon"], *region.get("holes", [])] + edges = [] + for ring in rings: + pts = np.asarray(ring, dtype=float) * [s, 1.0] + if np.allclose(pts[0], pts[-1]): + pts = pts[:-1] + edges.append(np.column_stack([pts, np.roll(pts, -1, axis=0)])) + polys.append((np.concatenate(edges), index[region["material_id"]])) + self._build_columns(polys, pad) + gw = groundwater(model.get("groundwater")) + self.gw_mode = gw["mode"] + self.ru = float(gw.get("ru", 0.0)) + if self.gw_mode == "line": + line = sorted([s * x, z] for x, z in gw["line"]) + self.water_x = np.array([p[0] for p in line]) + self.water_z = np.array([p[1] for p in line]) + self.kh = float(model.get("seismic_kh") or 0.0) + loads = surcharges(model.get("surcharge")) + self.loads = np.array( + [sorted((s * q["x0"], s * q["x1"])) + [q["q"]] for q in loads], dtype=float + ).reshape(-1, 3) + crack = model.get("tension_crack") or {} + self.crack_depth = float(crack.get("depth") or 0.0) + self.crack_water = float(crack.get("water") or 0.0) + + def _build_columns(self, polys: list, pad: int) -> None: + xs = np.unique(np.concatenate([e[:, 0] for e, _ in polys])) + self.X = xs + gap = np.diff(xs) + self.xa = xs[:-1] + 1e-7 * gap + self.xb = xs[1:] - 1e-7 * gap + q = np.concatenate([self.xa, self.xb]) + bots, tops, mats = [], [], [] + for edges, mat in polys: + x1, z1, x2, z2 = edges.T + lo, hi = np.minimum(x1, x2), np.maximum(x1, x2) + cross = (lo[None, :] < q[:, None]) & (q[:, None] < hi[None, :]) + with np.errstate(invalid="ignore", divide="ignore"): + z = z1 + (q[:, None] - x1) * (z2 - z1) / (x2 - x1) + z = np.sort(np.where(cross, z, np.nan), axis=1) + count = int(np.max(np.sum(cross, axis=1))) if len(q) else 0 + if count < 2: + continue + bots.append(z[:, 0 : count - 1 : 2]) + tops.append(z[:, 1:count:2]) + mats.append(np.full(bots[-1].shape, mat)) + bot = np.concatenate(bots, axis=1) + top = np.concatenate(tops, axis=1) + mat = np.concatenate(mats, axis=1) + mat = np.where(np.isnan(bot), pad, mat) + half = len(self.xa) + order = np.argsort(np.where(np.isnan(bot[:half]), np.inf, bot[:half]), axis=1) + self.bot_a, self.bot_b = (np.take_along_axis(bot[k], order, 1) for k in _halves(half)) + self.top_a, self.top_b = (np.take_along_axis(top[k], order, 1) for k in _halves(half)) + self.mat = np.take_along_axis(mat[:half], order, 1) + # 지표면 꺾은선(정확) — 띠 양 끝으로 늘린 맨 윗면 · 연직 면은 같은 x 두 점 + ta = np.nanmax(self.top_a, axis=1) + tb = np.nanmax(self.top_b, axis=1) + slope = (tb - ta) / (self.xb - self.xa) + left = ta + (xs[:-1] - self.xa) * slope + right = tb + (xs[1:] - self.xb) * slope + sx = np.column_stack([xs[:-1], xs[1:]]).ravel() + sz = np.column_stack([left, right]).ravel() + keep = np.ones(len(sx), bool) + keep[1:] = (np.abs(np.diff(sx)) > 1e-9) | (np.abs(np.diff(sz)) > 1e-6) + self.surface_x, self.surface_z = sx[keep], sz[keep] + ba = np.nanmin(self.bot_a, axis=1) + bb = np.nanmin(self.bot_b, axis=1) + self.bottom_z = float(min(np.min(ba), np.min(bb))) + # 빈 칸 = 아주 낮은 두께 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)) + + # ---- 조회 ------------------------------------------------------------- + def surface(self, x: np.ndarray) -> np.ndarray: + """지표면 높이(연직 면은 낮은 쪽 · 모형 밖은 nan).""" + x = np.asarray(x, dtype=float) + z = np.interp(x, self.surface_x, self.surface_z) + return np.where((x < self.X[0]) | (x > self.X[-1]), np.nan, z) + + def columns(self, x: np.ndarray): + j = np.clip(np.searchsorted(self.X, x, side="right") - 1, 0, len(self.xa) - 1) + t = ((x - self.xa[j]) / (self.xb[j] - self.xa[j]))[..., None] + bot = self.bot_a[j] + t * (self.bot_b[j] - self.bot_a[j]) + top = self.top_a[j] + t * (self.top_b[j] - self.top_a[j]) + return bot, top, self.mat[j] + + def water(self, x: np.ndarray) -> np.ndarray: + return np.interp(x, self.water_x, self.water_z) + + +def _halves(half: int): + return (slice(0, half), slice(half, None)) + + +def slices(prep: Prepared, xe: np.ndarray, ze: np.ndarray) -> dict: + """활동면(절편 경계 x · 바닥 z, 꼴 (S, n+1)) → 절편 값 (S, n) · 잘못된 면은 valid False.""" + xe = np.atleast_2d(np.asarray(xe, float)) + ze = np.atleast_2d(np.asarray(ze, float)) + b = np.diff(xe, axis=1) + 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) + 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 + else: + g_dry, g_wet = prep.gamma_t[mat] * thick, 0.0 * thick + dry_mid = wet_mid = lo + 0.5 * thick + 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) + Ws = weight_per_b * b + if prep.gw_mode == "ru": + u = prep.ru * weight_per_b + # 바닥 물성 = 바닥 바로 위 구간(약층 바닥을 따라가는 면은 그 약층) + probe = zb[..., None] + 1e-6 + inside = (bot <= probe) & (probe < top) & (top - bot > 1e-9) + base_k = np.argmax(inside, axis=-1) + base_mat = np.take_along_axis(mat, base_k[..., None], -1)[..., 0] + found = np.any(inside, axis=-1) + c = prep.c[base_mat] + tanphi = prep.tanphi[base_mat] + # 상재하중 — 절편 폭과 겹친 길이만큼 정확히 · 작용점 = 겹친 구간 가운데 + Q = np.zeros_like(zb) + Mq = np.zeros_like(zb) + for x0, x1, q in prep.loads: + a = np.clip(xe[:, :-1], x0, x1) + bb = np.clip(xe[:, 1:], x0, x1) + ov = bb - a + Q += q * ov + Mq += q * ov * 0.5 * (a + bb) + W = Ws + Q + xw = np.where(W > 0, (Ws * xm + Mq) / np.maximum(W, _TINY), xm) + valid = np.all(found & (zb >= prep.bottom_z - 1e-6) & np.isfinite(zt), axis=1) + valid &= np.all(b > 0, axis=1) + return { + "xe": xe, + "ze": ze, + "b": b, + "xm": xm, + "zb": zb, + "zt": zt, + "alpha": alpha, + "sa": np.sin(alpha), + "ca": np.cos(alpha), + "l": b / np.cos(alpha), + "W": W, + "Ws": Ws, + "Q": Q, + "xw": xw, + "zg": zg, + "u": u, + "c": np.where(found, c, 0.0), + "tanphi": np.where(found, tanphi, 0.0), + "mat": base_mat, + "valid": valid, + "kh": prep.kh, + "Pw": np.zeros(len(xe)), + "zp": ze[:, 0].copy(), + } + + +def add_crack_water(sl: dict, depth: np.ndarray, fill: float) -> None: + """인장균열 물 — 머리 쪽 연직 면에 수평 수압 Pw = ½·γw·(f·z_c)² · 바닥에서 f·z_c/3.""" + hw = fill * np.asarray(depth, float) + sl["Pw"] = 0.5 * GAMMA_W * hw**2 + sl["zp"] = sl["ze"][:, 0] + hw / 3.0 + + +# ---- 방법 ------------------------------------------------------------------- +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 + + +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): + 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 + + +def _force_den(sl: dict, N: np.ndarray): + return (N * sl["sa"] + sl["kh"] * sl["Ws"]).sum(1) + sl["Pw"] + + +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 _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 _iterate(fun, F0): + """F = fun(F) 풀이(묶음 전체) — 수치 미분 뉴턴(h = fun(F) − F) · 안 되면 반쯤 섞은 고정점 걸음. + + 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −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) + h1 = f1 - F + delta = 1e-6 * np.maximum(np.abs(F), 0.1) + f2, _ = fun(F + delta) + dh = ((f2 - (F + delta)) - h1) / delta + 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) + Fn = np.where(use, newton, relaxed) + Fn = np.where(np.isfinite(Fn) & (Fn > 0), Fn, np.nan) + conv = np.abs(h1) < TOL * np.maximum(1.0, np.abs(F)) + F = np.where(done | conv, F, Fn) + iters += ~done + 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 + + +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) + # 같은 꼴 거르기 — 진출부 급경사(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) + return {"F": np.where(ok, F, np.nan), "iterations": np.ones(len(F), int), "converged": ok} + + +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) + return _finish(sl, F, it, conv, lam) + + +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) + out = _finish(sl, F, it, conv, lam) + out["F0"] = out["F"] + out["f0"] = janbu_f0(sl) + out["F"] = out["F0"] * out["f0"] + return out + + +def janbu_f0(sl: dict) -> np.ndarray: + xe, ze = sl["xe"], sl["ze"] + dx, dz = xe[:, -1:] - xe[:, :1], ze[:, -1:] - ze[:, :1] + L = np.hypot(dx, dz)[:, 0] + dist = np.abs(dz * (xe - xe[:, :1]) - dx * (ze - ze[:, :1]))[:, :] / np.maximum( + L[:, None], _TINY + ) + ratio = dist.max(1) / np.maximum(L, _TINY) + 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)) + return 1.0 + b1 * (ratio - 1.4 * ratio**2) + + +def spencer(sl: dict, center: np.ndarray, F0=None, scan: bool = True) -> dict: + """Spencer — 모멘트 · 힘 두 평형을 함께 맞추는 (F, λ) · 뉴턴(수치 야코비안) · 안 맞으면 λ 훑기. + + scan=False — 탐색 고리 안(시험 면 대부분이 버려짐)에서는 훑기를 건너뜀(시간 대부분이 여기). + """ + S = len(sl["W"]) + F = np.full(S, 1.5) if F0 is None else np.asarray(F0, float).copy() + start = F.copy() + lam = np.zeros(S) + # 엉뚱한 근(F → 0 쪽) 거르기 — 엄밀해는 출발값(Bishop)의 절반~두 배 안에 있음 + near = lambda f, ref: (f > 0.5 * ref) & (f < 2.0 * ref) # noqa: E731 + F, lam, it, conv = _spencer_newton(sl, center, F, lam) + conv &= near(F, start) + bad = ~conv & sl["valid"] + if scan and bad.any(): + sub = _subset(sl, bad) + Fs, ls = _spencer_scan(sub, center[bad]) + F2, l2, it2, c2 = _spencer_newton(sub, center[bad], Fs, ls) + F[bad], lam[bad], it[bad], conv[bad] = F2, l2, it[bad] + it2, c2 & near(F2, start[bad]) + out = _finish(sl, F, it, conv, lam) + out["lambda"] = np.where(np.isfinite(out["F"]), lam, np.nan) + 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, lam = F.copy(), lam.copy() + it = np.zeros(len(F), int) + done = np.zeros(len(F), bool) + 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 + det = j11 * j22 - j12 * j21 + det = np.where(np.abs(det) < 1e-14, np.nan, det) + dF = -(j22 * rm - j12 * rf) / det + dl = -(-j21 * rm + j11 * rf) / det + scale = np.minimum( + 1.0, + np.minimum( + 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) + 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) + S = len(sl["W"]) + diff = np.full((len(grid), S), np.nan) + fm_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 + 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): + d0, d1 = diff[k], diff[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) + lam_c = grid[k] + t * (grid[k + 1] - grid[k]) + F_c = fm_all[k] + t * (fm_all[k + 1] - fm_all[k]) + better = cross & (np.abs(lam_c) < best_gap) + best_gap = np.where(better, np.abs(lam_c), best_gap) + best_F = np.where(better, F_c, best_F) + best_l = np.where(better, lam_c, best_l) + return best_F, best_l + + +def _subset(sl: dict, mask: np.ndarray) -> dict: + out = { + k: (v[mask] if isinstance(v, np.ndarray) and v.shape[:1] == mask.shape else v) + for k, v in sl.items() + } + return out + + +def _finish(sl, F, it, conv, lam): + _, R, den, *_ = _parts(sl, np.where(np.isfinite(F) & (F > 0), F, 1.0), lam) + 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 { + "F": np.where(ok, F, np.nan), + "iterations": it, + "converged": conv, + "malpha_ok": malpha_ok, + } + + +def slice_table(sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bishop") -> list[dict]: + """한 면(k 번째)의 절편표 — 폭 · 높이 · 경사 · 무게 · 수압 · 강도 · N' · 동원 전단 · (Spencer) E · X.""" + 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)) + 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) + R = one["c"] * one["l"] + (N - one["u"] * one["l"]) * one["tanphi"] + rows = [] + for i in range(N.shape[1]): + rows.append( + { + "x_left": float(one["xe"][0, i]), + "x_right": float(one["xe"][0, i + 1]), + "width": float(one["b"][0, i]), + "height": float(one["zt"][0, i] - one["zb"][0, i]), + "alpha_deg": float(np.degrees(one["alpha"][0, i])), + "weight": float(one["W"][0, i]), + "surcharge": float(one["Q"][0, i]), + "pore_pressure": float(one["u"][0, i]), + "base_length": float(one["l"][0, i]), + "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]), + "shear_mobilized": float(R[0, i] / F), + "E_left": float(E[0, i]), + "X_left": float(lam * E[0, i]), + } + ) + return rows diff --git a/E03_StabSection/E03_StabSection_Engine_Model.py b/E03_StabSection/E03_StabSection_Engine_Model.py new file mode 100644 index 000000000..c706dca2c --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_Model.py @@ -0,0 +1,293 @@ +"""E03 2D 지층 모델 빌더 — 측점 횡단(지반선 · 설계선 · 소단 · 측구 · 암 경계 · 성토체 · 벽) → 지층 폴리곤. + +`build_section_model` = storage 를 모르는 순수 함수(좌표 · 값만) · `section_model_from_b06` = B06 횡단 자료 꼴 +(분석 345-D ②)을 읽는 얇은 어댑터. 좌표 x = B06 offset(좌 +) · z = 표고. +설계선(design_line)이 소단 · 측구 · 노면을 이미 담음 — 여기서 다시 그리지 않음. 설계선 밖은 원지반. +암 경계 = 원지반을 연직으로 |offset| 내린 선(B06 면적 가름과 같은 정의 · `Engine_Design.py` 363-377). +""" + +import numpy as np +from shapely.geometry import LineString, MultiPolygon, Polygon, box +from shapely.ops import unary_union +from shapely.validation import make_valid + +from E03_StabSection.E03_StabSection_Types import ( + DEFAULT_MATERIALS, + DEFAULT_ROAD_SURCHARGE_KPA, +) +from E03_StabSection.E03_StabSection_Types import ( + groundwater as norm_groundwater, +) +from E03_StabSection.E03_StabSection_Types import ( + surcharges as norm_surcharges, +) + +#: B06 지반 라벨 → 암 경계 아래 지층(요구서 345 6장 · 분석 D 5-2). 경계 위는 풍화토. +ROCK_LAYER_OF = {"ripping_rock": "weathered_rock", "blasting_rock": "soft_rock"} +SOIL_ID = "weathered_soil" +FILL_ID = "fill" +#: 벽(C군) — 활동면이 못 지나게 강한 물체로(콘크리트 · 석축 대표값). 물성 표에 없으면 이 값을 덧붙임. +WALL_MATERIAL = { + "id": "wall", + "name": "벽체", + "gamma_t": 23.0, + "gamma_sat": 23.0, + "c": 1000.0, + "phi": 45.0, +} +DEFAULT_HALF_WIDTH_M = 40.0 + + +def _polys(geom) -> list[Polygon]: + geom = make_valid(geom) + if isinstance(geom, Polygon): + parts = [geom] + elif isinstance(geom, MultiPolygon): + parts = list(geom.geoms) + else: + parts = [g for g in getattr(geom, "geoms", []) if isinstance(g, Polygon)] + return [p for p in parts if p.area > 1e-6] + + +def _region(material_id: str, poly: Polygon) -> dict: + out = {"material_id": material_id, "polygon": [list(c) for c in poly.exterior.coords]} + if poly.interiors: + out["holes"] = [[list(c) for c in ring.coords] for ring in poly.interiors] + return out + + +def final_surface(ground: list, design_line: list) -> list[list[float]]: + """공사 뒤 지표 = 설계선 범위는 설계선 · 밖은 원지반.""" + line = sorted((float(x), float(z)) for x, z in design_line) + if not line: + return [list(p) for p in ground] + x0, x1 = line[0][0], line[-1][0] + left = [p for p in ground if p[0] < x0 - 1e-9] + right = [p for p in ground if p[0] > x1 + 1e-9] + return [list(map(float, p)) for p in left + line + right] + + +def build_section_model( + ground: list, + design_line: list, + *, + layers: list, + materials: list[dict] | None = None, + fill_material: str = FILL_ID, + walls: list[dict] | None = None, + groundwater: dict | None = None, + seismic_kh: float = 0.0, + surcharge: list[dict] | None = None, + bottom_z: float | None = None, + area_range: tuple[float, float] | None = None, +) -> dict: + """원지반 · 설계선 · 지층(위부터 [(material_id, 원지반 아래 바닥 깊이 m | None=모형 바닥)]) → 해석 모형. + + 돌려주는 값 = LEM 모형(materials · regions · surface · groundwater · seismic_kh · surcharge) + + `areas`(area_range 안 절토 · 성토 · 절토 토사 · 절토 암 ㎡ — B06 면적 대조용). + """ + ground = sorted((float(x), float(z)) for x, z in ground) + surface = final_surface(ground, design_line) + xs = [p[0] for p in ground] + zs = [p[1] for p in ground] + [p[1] for p in surface] + relief = max(zs) - min(zs) + zb = bottom_z if bottom_z is not None else min(zs) - max(5.0, 0.5 * relief) + x0, x1 = xs[0], xs[-1] + bottom = [(x1, zb), (x0, zb)] + below_final = make_valid(Polygon(surface + bottom)) + ground_body = make_valid(Polygon(ground + bottom)) + natural = below_final.intersection(ground_body) + fill = below_final.difference(ground_body) + cut = ground_body.difference(below_final) + + def band(d_top: float, d_bot: float | None) -> Polygon: + top = [(x, z - d_top) for x, z in ground] + if d_bot is None: + low = [(x1, zb - 1.0), (x0, zb - 1.0)] + else: + low = [(x, z - d_bot) for x, z in reversed(ground)] + return make_valid(Polygon(top + low)) + + wall_geom = unary_union([Polygon(w["polygon"]) for w in walls or []]) if walls else None + regions: list[dict] = [] + bands: list = [] + depth = 0.0 + for material_id, d_bot in layers: + b = band(depth, None if d_bot is None else float(d_bot)) + bands.append(b) + part = natural.intersection(b) + if wall_geom is not None: + part = part.difference(wall_geom) + regions += [_region(material_id, p) for p in _polys(part)] + if d_bot is None: + break + depth = float(d_bot) + fill_part = fill.difference(wall_geom) if wall_geom is not None else fill + regions += [_region(fill_material, p) for p in _polys(fill_part)] + mats = [dict(m) for m in (materials or DEFAULT_MATERIALS)] + if wall_geom is not None: + for w in walls: + regions += [ + _region(w.get("material_id", "wall"), p) for p in _polys(Polygon(w["polygon"])) + ] + if not any(m["id"] == "wall" for m in mats): + mats.append(dict(WALL_MATERIAL)) + a0, a1 = area_range or (x0, x1) + clip = box(a0, zb - 1.0, a1, max(zs) + 100.0) + cut_in = cut.intersection(clip) + cut_soil = cut_in.intersection(bands[0]).area if len(bands) > 1 else cut_in.area + return { + "materials": mats, + "regions": regions, + "surface": surface, + "groundwater": norm_groundwater(groundwater), + "seismic_kh": float(seismic_kh or 0.0), + "surcharge": norm_surcharges(surcharge), + "bottom_z": zb, + "areas": { + "cut_m2": cut_in.area, + "fill_m2": fill.intersection(clip).area, + "cut_soil_m2": cut_soil, + "cut_rock_m2": cut_in.area - cut_soil, + }, + } + + +# ---- B06 어댑터 ----------------------------------------------------------------- +def ground_from_cross(cross: dict) -> list[tuple[float, float]]: + return sorted( + (float(s["offset_m"]), float(s["elevation_m"])) + for s in cross.get("samples") or [] + if s.get("valid") is not False and s.get("elevation_m") is not None + ) + + +def widen_ground(cross: dict, sampler, half_width: float = DEFAULT_HALF_WIDTH_M, step: float = 0.5): + """해석 폭 — 횡단 표본 바깥을 지형에서 다시 떠 ±half_width 까지 넓힘(지형 밖이 나오면 거기서 멈춤).""" + ground = ground_from_cross(cross) + frame = cross["frame"] + origin = np.array([float(frame["origin"]["x"]), float(frame["origin"]["y"])]) + axis = np.array([float(v) for v in frame["left_xy"]]) + added = [] + for sign, edge in ((-1.0, ground[0][0]), (1.0, ground[-1][0])): + n = int(max(0.0, half_width - abs(edge)) / step) + if n <= 0: + continue + offs = edge + sign * step * np.arange(1, n + 1) + z, valid = sampler.sample_xy(origin[None, :] + axis[None, :] * offs[:, None]) + stop = int(np.argmin(valid)) if not valid.all() else len(valid) + added += [(float(o), float(e)) for o, e in zip(offs[:stop], z[:stop])] + return sorted(ground + added) + + +def side_cases(design: dict, ground: list, line: list) -> list[dict]: + """좌 · 우 비탈마다 역할(절토 · 성토) · 활동 방향 · 탐색 한계. 좌 = +offset. + + 탐색 한계 = 노면 반대쪽 끝 ~ 비탈 끝(설계선이 원지반과 만나는 점) 바깥 max(5 m, 2H) — 설계 비탈을 + 지나는 면만 봄(바깥 자연 비탈 전체 붕괴는 E 광역 · 2차 몫). H = 노면 끝 ~ 비탈 끝 높이차. + """ + mode = design.get("section_mode") or "right_cut" + roles = { + "left_cut": ("cut", "fill"), + "right_cut": ("fill", "cut"), + "both_cut": ("cut", "cut"), + "both_fill": ("fill", "fill"), + }[mode] + edges = design.get("road_edges") or {} + left_edge = float((edges.get("left") or {}).get("offset_m", 2.0)) + right_edge = float((edges.get("right") or {}).get("offset_m", -2.0)) + gx = np.array([p[0] for p in ground]) + gz = np.array([p[1] for p in ground]) + lx = np.array([p[0] for p in line]) + lz = np.array([p[1] for p in line]) + off = np.abs(lz - np.interp(lx, gx, gz)) > 0.005 + out = [] + for side, role in zip(("left", "right"), roles): + inward = role == "cut" # 절토 = 산 쪽에서 노면으로 · 성토 = 노면에서 바깥으로 + sign = 1.0 if side == "left" else -1.0 + edge = left_edge if side == "left" else right_edge + # 비탈 끝 = 원지반에서 떨어진 가장 바깥 꼭짓점의 바로 바깥 꼭짓점(다시 원지반에 붙는 점) + idx = [k for k, (x, o) in enumerate(zip(lx, off)) if o and sign * x >= sign * edge - 1e-9] + if idx: + k = max(idx) if sign > 0 else min(idx) + k = min(max(k + int(sign), 0), len(lx) - 1) + day = float(lx[k]) + else: + day = edge + H = abs(float(np.interp(day, lx, lz)) - float(np.interp(edge, lx, lz))) + reach = day + sign * max(5.0, 2.0 * H) + reach = min(max(reach, gx[0]), gx[-1]) + inner = right_edge if side == "left" else left_edge + limits = tuple(sorted((inner, reach))) + direction = int(-sign) if inward else int(sign) + out.append( + {"side": side, "role": role, "direction": direction, "limits": limits, "height_m": H} + ) + return out + + +def section_model_from_b06( + cross: dict, + design: dict, + *, + materials: list[dict] | None = None, + layers: list | None = None, + ground: list | None = None, + wet: bool = False, + seismic_kh: float = 0.0, + road_surcharge_kpa: float = DEFAULT_ROAD_SURCHARGE_KPA, + walls: list[dict] | None = None, +) -> dict: + """B06 횡단 파일(samples · frame) + 설계(design) → 해석 모형 + 비탈 목록(`sides`). + + ground = 넓혀 다시 뜬 지반선(`widen_ground`) · 없으면 횡단 표본 그대로. + layers 기본 = 토사 측점 풍화토 한 층 · 암 측점 풍화토(암 경계까지) + 풍화암 / 연암. + wet = 우기 — 지하수위 = 공사 뒤 지표(KDS 11 70 05 표 4.3-3 · 토사 · 풍화암 지표 포화). + 노면 하중 = 노면 끝 사이 등분포. 벽 도형(walls)은 B06 화면 산식(TS)이 그리므로 부르는 쪽이 넘김. + """ + ground = ground or ground_from_cross(cross) + line = [ + (float(p["offset_m"]), float(p["elevation_m"])) for p in design.get("design_line") or [] + ] + if layers is None: + rock = ROCK_LAYER_OF.get(design.get("ground_type")) + offset = design.get("rock_boundary_offset_m") + if rock and offset is not None: + layers = [(SOIL_ID, abs(float(offset))), (rock, None)] + elif rock: + layers = [(rock, None)] + else: + layers = [(SOIL_ID, None)] + edges = design.get("road_edges") or {} + loads = [] + if road_surcharge_kpa and edges.get("left") and edges.get("right"): + loads = [ + { + "x0": edges["right"]["offset_m"], + "x1": edges["left"]["offset_m"], + "q": road_surcharge_kpa, + } + ] + surface = final_surface(ground, line) + gw = {"mode": "line", "line": surface} if wet else {"mode": "none"} + model = build_section_model( + ground, + line, + layers=layers, + materials=materials, + walls=walls, + groundwater=gw, + seismic_kh=seismic_kh, + surcharge=loads, + area_range=(line[0][0], line[-1][0]) if line else None, + ) + model["sides"] = side_cases(design, ground, line or ground) + model["station"] = {"chainage_m": cross.get("chainage_m"), "label": cross.get("label")} + return model + + +def polygon_closed(region: dict) -> bool: + """폴리곤 닫힘 · 자기 교차 없음 · 넓이 > 0 (시험용).""" + ring = region["polygon"] + poly = Polygon(ring, region.get("holes")) + return bool(poly.is_valid and poly.area > 0 and LineString(ring).is_closed) diff --git a/E03_StabSection/E03_StabSection_Engine_Search.py b/E03_StabSection/E03_StabSection_Engine_Search.py new file mode 100644 index 000000000..3ce607bb3 --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_Search.py @@ -0,0 +1,634 @@ +"""E03 임계 활동면 탐색 — 원호(Grid & Radius · Entry & Exit · 다듬기) · 비원호(블록 + 꼭짓점 최적화). + +활동 방향 +x 로 맞춘 `Prepared` 위에서 찾고 좌표를 원래 방향으로 되돌려 냄. +다듬기 · 비원호 = 차분진화(best/1/bin · 무리 전체를 한 번에 평가 · seed 고정 → 같은 입력 같은 답) ++ 비원호 마지막 국부 다듬기 = SciPy Nelder-Mead. SciPy differential_evolution 은 개체마다 파이썬 돌연변이 고리라 +한 비탈 시간의 30 % 를 먹어 직접 묶음 꼴로 씀. +상용 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 scipy.optimize import minimize + +from E03_StabSection import E03_StabSection_Engine_LEM as lem + +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 + + +# ---- 원호 --------------------------------------------------------------------- +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) + 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) + return xc[ok], zc[ok], R[ok], xE[ok], xX[ok] + + +def grid_centers(prep, limits, nx=15, nz=15): + xl, xr = limits + 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) + 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=60): + """원호 다듬기 — (진입 x, 진출 x, 중심각) 차분진화 · 첫 무리 = 거친 탐색의 좋은 면.""" + fun = _objective(method) + lo = np.array([limits[0], limits[0], 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 + + x, fmin, _pop, _f = _de(obj, lo, hi, _population(seeds, lo, hi, 30, seed), maxiter, seed) + 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])) + + +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 _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=60, polish=True +): + """비원호 꼭짓점 최적화 — 차분진화 뒤 가장 좋은 면을 Nelder-Mead 로 국부 다듬기(상용 SAO · Optimize 대응).""" + 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) + 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, 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) + 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 + 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) + 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 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, + *, + 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])))) + results: dict = {m: {} for m in methods} + counts = {} + contour = 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",) + 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)) + ], + } + for m in methods: + 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 fmin < _BIG: + results[m]["circular"] = describe(prep, "circle", geom, m, n_slices, min_depth) + best_circles[m] = geom + 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) + found = [] + for m in [k for k in ("janbu", "spencer") if k in methods]: + fmin, _x, geom, top = optimize_polyline( + prep, lim, m, seeds, n_slices, min_depth, mv, seed + ) + 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: + cands = [v for v in results[m].values() if v.get("F") is not None] + 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, + "counts": counts, + "elapsed_s": round(time.perf_counter() - t0, 3), + } + + +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 diff --git a/E03_StabSection/E03_StabSection_Engine_Station.py b/E03_StabSection/E03_StabSection_Engine_Station.py new file mode 100644 index 000000000..e6f9438dd --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_Station.py @@ -0,0 +1,75 @@ +"""E03 측점 하나 해석 — 좌 · 우 비탈 × 조건(건기 · 우기 · 지진) × 방법. + +조건 = KDS 11 70 05 표 4.3-3: 건기 u = 0 · 우기 지하수위 = 지표 · 지진 = 건기 수위 + kh. +판정(기준안전율 대조)은 화면 · 보고서 몫 — 여기는 값만 냄. +""" + +import time + +from E03_StabSection.E03_StabSection_Engine_Model import section_model_from_b06 +from E03_StabSection.E03_StabSection_Engine_Search import METHODS, analyze + +#: 활동면 최소 깊이(m) — 지표를 얇게 긁는 면(표토 박리)은 사면 붕괴로 보지 않음. +MIN_DEPTH_M = 0.3 + + +def analyze_station( + cross: dict, + design: dict, + *, + materials=None, + layers=None, + ground=None, + conditions=("dry", "wet"), + seismic_kh: float = 0.0, + methods=METHODS, + noncircular: bool = True, + min_depth: float = MIN_DEPTH_M, + **model_kw, +) -> dict: + t0 = time.perf_counter() + out = {"chainage_m": cross.get("chainage_m"), "label": cross.get("label"), "cases": []} + for cond in conditions: + model = section_model_from_b06( + cross, + design, + materials=materials, + layers=layers, + ground=ground, + wet=cond == "wet", + seismic_kh=seismic_kh if cond == "seismic" else 0.0, + **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}) + out["elapsed_s"] = round(time.perf_counter() - t0, 3) + return out + + +def analyze_stations(items: list[dict], workers: int | None = None, **kw) -> list[dict]: + """전 측점 — 측점마다 따로 도니 프로세스로 나눠 돎. items = [{"cross", "design", "ground"?}, …] · 순서 그대로.""" + import os + from concurrent.futures import ProcessPoolExecutor + + workers = workers or max(1, (os.cpu_count() or 2) - 2) + if workers == 1 or len(items) < 2: + return [ + analyze_station(it["cross"], it["design"], ground=it.get("ground"), **kw) + for it in items + ] + with ProcessPoolExecutor(max_workers=workers) as pool: + futures = [ + pool.submit(analyze_station, it["cross"], it["design"], ground=it.get("ground"), **kw) + for it in items + ] + return [f.result() for f in futures] diff --git a/resources/tester/test_532_2_stab_lem_benchmarks.py b/resources/tester/test_532_2_stab_lem_benchmarks.py new file mode 100644 index 000000000..bc5da1f29 --- /dev/null +++ b/resources/tester/test_532_2_stab_lem_benchmarks.py @@ -0,0 +1,348 @@ +# -*- coding: utf-8 -*- +"""532-2 E03 2D 한계평형 솔버 — 공개 검증 예제 대조(±3 %). + +출처(기대값 = 원문 그대로): +- S2VM = Rocscience Slide2 Slope Stability Verification Manual (2022) — 지식DB + `original/안정성평가/참고문헌/Rocscience_Slide2_SlopeStabilityVerification.pdf` · 문제 번호 = 그 책 번호. +- ACADS = Giam & Donald (1989) ACADS 사면 시험문제(S2VM #1~#9 에 실림). +- FK = Fredlund & Krahn (1977) Can. Geotech. J. 14(3) (S2VM #21). +- AT = Arai & Tagyo (1985) · YU = Yamagami & Ueta (1988) · LOUK = Loukidis et al. (2003) (S2VM #14~#17 · #62). +- 지정 활동면(원 · 꺾은선)은 방법 대조 · 탐색 문제는 상용 SW 탐색 결과 대조. +- Ordinary(Fellenius) 간극수압 = F&K 식 N = W·cosα − u·l (지식DB 02 §2.1) — Slide 는 (W − u·b)·cosα 라 r_u 문제 값이 다름. + + ./venv/Scripts/python.exe -m pytest resources/tester/test_532_2_stab_lem_benchmarks.py -q +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import pytest +from shapely.geometry import Polygon, box + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) + +from E03_StabSection.E03_StabSection_Engine_Search import analyze, evaluate_specified # noqa: E402 + +TOL = 0.03 + + +def mats(*rows): + return [{"id": i, "name": i, "gamma_t": g, "c": c, "phi": p} for i, g, c, p in rows] + + +def near(got, expected, tol=TOL): + assert got is not None, "안전율 없음" + assert abs(got / expected - 1) <= tol, ( + f"{got:.4f} vs {expected} ({(got / expected - 1) * 100:+.2f} %)" + ) + + +def circ(res, method): + return (res["methods"][method].get("circular") or {}).get("F") + + +def nonc(res, method): + return (res["methods"][method].get("noncircular") or {}).get("F") + + +ACADS_EXT = [[20, 20], [70, 20], [70, 35], [50, 35], [30, 25], [20, 25]] +ACADS1C_REGIONS = [ + { + "material_id": "soil1", + "polygon": [[30, 25], [50, 35], [70, 35], [70, 31], [54, 31], [50, 29], [40, 27]], + }, + { + "material_id": "soil2", + "polygon": [[40, 27], [50, 29], [54, 31], [70, 31], [70, 24], [52, 24]], + }, + { + "material_id": "soil3", + "polygon": [[20, 20], [70, 20], [70, 24], [52, 24], [40, 27], [30, 25], [20, 25]], + }, +] +ACADS1C_MATS = mats(("soil1", 19.5, 0, 38), ("soil2", 19.5, 5.3, 23), ("soil3", 19.5, 7.2, 20)) + + +def test_acads_1a_homogeneous_circular_search(): + # S2VM #1 (ACADS 1a) — Bishop 0.987 · Spencer 0.986 · Janbu 보정 0.990 · 심판값 1.00 + m = { + "materials": mats(("s", 20, 3, 19.6)), + "regions": [{"material_id": "s", "polygon": ACADS_EXT}], + } + r = analyze(m, direction=-1, noncircular=False) + near(circ(r, "bishop"), 0.987) + near(circ(r, "spencer"), 0.986) + near(circ(r, "janbu"), 0.990) + + +def test_acads_1b_tension_crack_full_of_water(): + # S2VM #2 (ACADS 1b) — 균열 깊이 4 m · 물 가득 — Bishop 1.596 · Spencer 1.592 · Janbu 보정 1.489 + m = { + "materials": mats(("s", 20, 32, 10)), + "regions": [{"material_id": "s", "polygon": ACADS_EXT}], + "tension_crack": {"depth": 4.0, "water": 1.0}, + } + r = analyze(m, direction=-1, noncircular=False) + near(circ(r, "bishop"), 1.596) + near(circ(r, "spencer"), 1.592) + near(circ(r, "janbu"), 1.489) + + +def test_acads_1c_three_layers(): + # S2VM #3 (ACADS 1c) — Bishop 1.405 · Spencer 1.375 · Janbu 보정 1.357 + r = analyze( + {"materials": ACADS1C_MATS, "regions": ACADS1C_REGIONS}, direction=-1, noncircular=False + ) + near(circ(r, "bishop"), 1.405) + near(circ(r, "spencer"), 1.375) + near(circ(r, "janbu"), 1.357) + + +def test_acads_1d_three_layers_seismic(): + # S2VM #4 (ACADS 1d) — kh 0.15 — Bishop 1.016 · Spencer 0.991 · Janbu 보정 0.965 + m = {"materials": ACADS1C_MATS, "regions": ACADS1C_REGIONS, "seismic_kh": 0.15} + r = analyze(m, direction=-1, noncircular=False) + near(circ(r, "bishop"), 1.016) + near(circ(r, "spencer"), 0.991) + near(circ(r, "janbu"), 0.965) + + +def test_acads_2b_embankment_dam_specified_circle(): + # S2VM #6 (ACADS 2b Talbingo 댐 · 성토) — 지정 원 (100.3, 291.0, R 278.8) + # Bishop 2.208 · Spencer 2.292 · Janbu 보정 2.073 · 심판값 2.29 + ext = [ + [0, 0], [315.5, 162], [319.5, 162], [321.6, 162], [327.6, 162], [386.9, 130.6], [394.1, 130.6], + [453.4, 97.9], [460.6, 97.9], [515, 65.3], [521.1, 65.3], [577.9, 31.4], [585.1, 31.4], [648, 0], + ] # fmt: skip + core = [ + [200.7, 0], + [311.9, 130.6], + [319.5, 162], + [321.6, 162], + [328.8, 146.1], + [331.3, 130.6], + [307.1, 0], + ] + rest = Polygon(ext).difference(Polygon(core)) + regions = [{"material_id": "core", "polygon": core}] + [ + {"material_id": "rockfill", "polygon": list(p.exterior.coords)} + for p in getattr(rest, "geoms", [rest]) + ] + m = {"materials": mats(("rockfill", 20.4, 0, 45), ("core", 18.1, 85, 23)), "regions": regions} + r = evaluate_specified(m, circle=[100.3, 291.0, 278.8], direction=-1) + near(r["bishop"]["F"], 2.208) + near(r["spencer"]["F"], 2.292) + near(r["janbu"]["F"], 2.073) + + +def _acads3_model(): + ext = Polygon([[20, 20], [84, 20], [84, 40], [67.5, 40], [43, 27.75], [20, 27.75]]) + regions = [ + { + "material_id": "soil", + "polygon": list(ext.intersection(box(0, 27, 100, 100)).exterior.coords), + }, + { + "material_id": "seam", + "polygon": list(ext.intersection(box(0, 26.5, 100, 27)).exterior.coords), + }, + { + "material_id": "soil", + "polygon": list(ext.intersection(box(0, 0, 100, 26.5)).exterior.coords), + }, + ] + return { + "materials": mats(("soil", 18.84, 28.5, 20), ("seam", 18.84, 0, 10)), + "regions": regions, + "groundwater": {"mode": "line", "line": [[20, 26.5], [84, 26.5]]}, + } + + +def test_acads_3b_weak_seam_specified_polyline(): + # S2VM #8 (ACADS 3b) — 약층 + 지하수위 · 지정 비원호 — Spencer 1.277 · Janbu 보정 1.294 + r = evaluate_specified( + _acads3_model(), + polyline=[[41.85, 27.75], [44.0, 26.5], [63.5, 27.0], [73.31, 40.0]], + direction=-1, + methods=("janbu", "spencer"), + ) + near(r["spencer"]["F"], 1.277) + near(r["janbu"]["F"], 1.294) + + +def test_acads_3a_weak_seam_noncircular_search(): + # S2VM #7 (ACADS 3a) — 블록 탐색 Janbu 보정 1.275 · Spencer 1.258(심판값 1.24~1.27) + r = analyze(_acads3_model(), direction=-1, methods=("janbu", "spencer")) + near(nonc(r, "janbu"), 1.275) + near(nonc(r, "spencer"), 1.258) + + +def test_acads_4_inclined_seam_water_surcharge(): + # 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)}] + [ + {"material_id": "soil", "polygon": list(p.exterior.coords)} + for p in ext.difference(seam).geoms + ] + linear = [{"x0": 70 + i, "x1": 71 + i, "q": 20 + 20 * (i + 0.5) / 10} for i in range(10)] + m = { + "materials": mats(("soil", 18.84, 28.5, 20), ("seam", 18.84, 0, 10)), + "regions": regions, + "surcharge": [{"x0": 23, "x1": 43, "q": 20}] + linear, + "groundwater": { + "mode": "line", + "line": [ + [20, 27.75], + [43, 27.75], + [49, 29.8], + [60, 34], + [66, 35.8], + [74, 37.6], + [80, 38.4], + [84, 38.4], + ], + }, + } + r = analyze(m, direction=-1, methods=("janbu", "spencer")) + near(nonc(r, "spencer"), 0.707) + near(nonc(r, "janbu"), 0.699) + + +FK_EXT = [[0, 0], [180, 0], [180, 20], [140, 20], [60, 60], [0, 60]] + + +def test_fredlund_krahn_specified_circle_dry_and_ru(): + # S2VM #21 (F&K 1977, ft · lbf) — 지정 원 (120, 90, R 80) + # 건기 Ordinary 1.928 · Bishop 2.080 · Spencer 2.073 / r_u 0.25 Ordinary 1.607 · Bishop 1.766 · Spencer 1.761 + m = { + "materials": mats(("s", 120, 600, 20)), + "regions": [{"material_id": "s", "polygon": FK_EXT}], + } + r = evaluate_specified( + m, circle=[120, 90, 80], direction=1, methods=("ordinary", "bishop", "spencer") + ) + near(r["ordinary"]["F"], 1.928) + near(r["bishop"]["F"], 2.080) + near(r["spencer"]["F"], 2.073) + m["groundwater"] = {"mode": "ru", "ru": 0.25} + r = evaluate_specified( + m, circle=[120, 90, 80], direction=1, methods=("ordinary", "bishop", "spencer") + ) + near(r["ordinary"]["F"], 1.607) + near(r["bishop"]["F"], 1.766) + near(r["spencer"]["F"], 1.761) + + +def test_loukidis_seismic_with_ru(): + # S2VM #62 — 건기 kh 0.432: Bishop 0.991 · Spencer 1.001 / r_u 0.5 + kh 0.132: Bishop 0.987 · Spencer 1.001 + ext = [[-50, -25], [150, -25], [150, 25], [75, 25], [0, 0], [-50, 0]] + m = { + "materials": mats(("c", 20, 25, 30)), + "regions": [{"material_id": "c", "polygon": ext}], + "seismic_kh": 0.432, + } + r = analyze(m, direction=-1, noncircular=False, methods=("bishop", "spencer")) + near(circ(r, "bishop"), 0.991) + near(circ(r, "spencer"), 1.001) + m.update(seismic_kh=0.132, groundwater={"mode": "ru", "ru": 0.5}) + r = analyze(m, direction=-1, noncircular=False, methods=("bishop", "spencer")) + near(circ(r, "bishop"), 0.987) + near(circ(r, "spencer"), 1.001) + + +AT_EXT = [[0, 0], [66, 0], [66, 35], [48, 35], [18, 15], [0, 15]] + + +def test_arai_tagyo_3_water_table(): + # S2VM #16 — 지하수위선 · 원호 Bishop 1.118 · Spencer 1.118 · Janbu 보정 1.131 / + # 비원호 Janbu 보정 1.050 · Spencer 1.094 + m = { + "materials": mats(("s", 18.82, 41.65, 15)), + "regions": [{"material_id": "s", "polygon": AT_EXT}], + "groundwater": {"mode": "line", "line": [[0, 15], [18, 15], [30, 23], [48, 29], [66, 32]]}, + } + r = analyze(m, direction=-1) + near(circ(r, "bishop"), 1.118) + near(circ(r, "spencer"), 1.118) + near(circ(r, "janbu"), 1.131) + near(nonc(r, "janbu"), 1.050) + near(nonc(r, "spencer"), 1.094) + + +def test_arai_tagyo_2_weak_middle_layer(): + # S2VM #15 — 약한 가운데 층 · Bishop 0.420 · Spencer 0.409 · Janbu 보정 0.423 + regions = [ + {"material_id": "upper", "polygon": [[24, 19], [48, 35], [72, 35]]}, + { + "material_id": "middle", + "polygon": [[0, 3], [96, 35], [72, 35], [24, 19], [18, 15], [0, 15]], + }, + {"material_id": "lower", "polygon": [[0, 3], [96, 3], [96, 35]]}, + ] + m = { + "materials": mats( + ("upper", 18.82, 29.4, 12), ("middle", 18.82, 9.8, 5), ("lower", 18.82, 294, 40) + ), + "regions": regions, + } + r = analyze(m, direction=-1, noncircular=False) + near(circ(r, "bishop"), 0.420) + near(circ(r, "spencer"), 0.409) + near(circ(r, "janbu"), 0.423) + + +def test_yamagami_ueta_ordinary_vs_bishop(): + # S2VM #17 — Bishop 1.344 · Ordinary 1.278 / 비원호 Spencer 1.325 + ext = [[0, 0], [25, 0], [25, 10], [15, 10], [5, 5], [0, 5]] + m = { + "materials": mats(("s", 17.64, 9.8, 10)), + "regions": [{"material_id": "s", "polygon": ext}], + } + r = analyze(m, direction=-1, methods=("ordinary", "bishop", "spencer")) + near(circ(r, "bishop"), 1.344) + near(circ(r, "ordinary"), 1.278) + near(nonc(r, "spencer"), 1.325) + + +def test_mirror_direction_same_factor(): + # 좌우 뒤집은 같은 비탈 = 같은 안전율 · 임계면 좌표도 거울 + m = { + "materials": mats(("s", 20, 3, 19.6)), + "regions": [{"material_id": "s", "polygon": ACADS_EXT}], + } + mirrored = {**m, "regions": [{"material_id": "s", "polygon": [[-x, z] for x, z in ACADS_EXT]}]} + a = evaluate_specified(m, circle=[29.104, 55.527, 30.527], direction=-1) + b = evaluate_specified(mirrored, circle=[-29.104, 55.527, 30.527], direction=1) + for k in ("ordinary", "bishop", "janbu", "spencer"): + assert abs(a[k]["F"] - b[k]["F"]) < 1e-9 + assert abs(a["bishop"]["entry"][0] + b["bishop"]["entry"][0]) < 1e-9 + near(a["bishop"]["F"], 0.987) # S2VM #1 Bishop 임계 원 그대로 + + +def test_slice_table_and_convergence_info(): + 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",))[ + "spencer" + ] + assert r["converged"] and r["iterations"] > 0 and -1.25 <= r["lambda"] <= 0.25 + rows = r["slices"] + assert len(rows) == 50 and all(row["width"] > 0 for row in rows) + # 절편간 힘 E — 머리(0) 에서 시작해 발끝에서 다시 0 근처(힘 평형이 닫힘) + weights = sum(row["weight"] for row in rows) + assert abs(rows[0]["E_left"]) < 1e-9 + assert {"alpha_deg", "pore_pressure", "normal_eff", "shear_mobilized", "material"} <= rows[ + 0 + ].keys() + assert weights > 0 + + +if __name__ == "__main__": + sys.exit(pytest.main([__file__, "-q"])) diff --git a/resources/tester/test_532_2_stab_section_model.py b/resources/tester/test_532_2_stab_section_model.py new file mode 100644 index 000000000..909ed9544 --- /dev/null +++ b/resources/tester/test_532_2_stab_section_model.py @@ -0,0 +1,155 @@ +# -*- coding: utf-8 -*- +"""532-2 E03 2D 지층 모델 빌더 — 폴리곤 닫힘 · 지층 · 성토체 · 절성토 면적 = B06 · 비탈 목록 · 측점 해석 꼴. + +- 순수 함수(`build_section_model`)는 만든 단면으로 · B06 대조는 검증_sub4 초기값 스냅샷(읽기만 · 없으면 건너뜀). + + ./venv/Scripts/python.exe -m pytest resources/tester/test_532_2_stab_section_model.py -q +""" + +from __future__ import annotations + +import json +import sys +from pathlib import Path + +import pytest +from shapely.geometry import Polygon + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) + +from E03_StabSection.E03_StabSection_Engine_Model import ( # noqa: E402 + build_section_model, + polygon_closed, + section_model_from_b06, + side_cases, +) +from E03_StabSection.E03_StabSection_Engine_Station import analyze_station # noqa: E402 + +SUB4 = ROOT / "storage/1/3/c725d747-50ae-46b0-87ce-4d4e56a53cf3/initial_snapshot" + +# 오른쪽(−) 산 · 왼쪽(+) 골 · 편절편성: 노면 −2~2 · 우 절토 1:0.5 · 좌 성토 1:1.5 +GROUND = [(x / 2, 10.0 - 0.6 * (x / 2)) for x in range(-40, 41)] # 10 − 0.6x (좌로 내려감) +DESIGN = [ + (-20.0, 22.0), + (-6.0, 13.6), + (-2.6, 10.0), + (-2.0, 10.0), + (2.0, 10.0), + (14.0, 2.0), + (20.0, -2.0), +] + + +def _trapezoid(ground, line, sign): + """설계선 꼭짓점 · 지반 표본 합집합 위 사다리꼴(B06 `_trapezoid_areas` 와 같은 정의) — 절토(+) / 성토(−).""" + import numpy as np + + gx, gz = zip(*ground) + lx, lz = zip(*line) + xs = sorted({*gx, *lx}) + xs = [x for x in xs if lx[0] <= x <= lx[-1]] + d = np.interp(xs, gx, gz) - np.interp(xs, lx, lz) + area = 0.0 + for x0, x1, d0, d1 in zip(xs, xs[1:], d, d[1:]): + if d0 * d1 < 0: + xm = x0 + (x1 - x0) * d0 / (d0 - d1) + area += 0.5 * (xm - x0) * max(sign * d0, 0) + 0.5 * (x1 - xm) * max(sign * d1, 0) + else: + area += 0.5 * (x1 - x0) * (max(sign * d0, 0) + max(sign * d1, 0)) + return area + + +def test_build_model_polygons_layers_fill_and_areas(): + m = build_section_model( + GROUND, + DESIGN, + layers=[("weathered_soil", 0.5), ("weathered_rock", None)], + area_range=(-20.0, 20.0), + ) + assert m["regions"] and all(polygon_closed(r) for r in m["regions"]) + ids = {r["material_id"] for r in m["regions"]} + assert ids == {"weathered_soil", "weathered_rock", "fill"} + # 지층 사이 겹침 없음 · 합 = 공사 뒤 지표 아래 전체 + polys = [Polygon(r["polygon"], r.get("holes")) for r in m["regions"]] + total = sum(p.area for p in polys) + union = polys[0] + for p in polys[1:]: + union = union.union(p) + assert abs(union.area - total) < 1e-6 + a = m["areas"] + assert abs(a["cut_m2"] - _trapezoid(GROUND, DESIGN, 1)) < 1e-6 + assert abs(a["fill_m2"] - _trapezoid(GROUND, DESIGN, -1)) < 1e-6 + assert ( + 0 < a["cut_soil_m2"] < a["cut_m2"] + and abs(a["cut_soil_m2"] + a["cut_rock_m2"] - a["cut_m2"]) < 1e-9 + ) + + +def test_side_cases_roles_directions_limits(): + design = { + "section_mode": "right_cut", + "road_edges": {"left": {"offset_m": 2.0}, "right": {"offset_m": -2.0}}, + } + left, right = side_cases(design, GROUND, DESIGN) + assert (left["role"], left["direction"]) == ("fill", 1) # 좌 성토 → 바깥(+x)으로 + assert (right["role"], right["direction"]) == ("cut", 1) # 우 절토 → 노면(+x)으로 + assert left["limits"][0] == -2.0 and right["limits"][1] == 2.0 + assert left["height_m"] > 0 and right["height_m"] > 0 + + +def test_analyze_station_shape(): + cross = { + "chainage_m": 0.0, + "label": "STA.0+000", + "samples": [{"offset_m": x, "elevation_m": z} for x, z in GROUND], + } + design = { + "section_mode": "right_cut", + "ground_type": "ripping_rock", + "rock_boundary_offset_m": -0.5, + "design_line": [{"offset_m": x, "elevation_m": z} for x, z in DESIGN], + "road_edges": { + "left": {"offset_m": 2.0, "elevation_m": 10.0}, + "right": {"offset_m": -2.0, "elevation_m": 10.0}, + }, + } + res = analyze_station( + cross, design, conditions=("dry",), methods=("bishop",), noncircular=False + ) + assert len(res["cases"]) == 2 + for case in res["cases"]: + F = case["min_F"]["bishop"] + assert F is not None and 0.3 < F < 20 + best = case["result"]["methods"]["bishop"]["circular"] + assert best["slices"] and best["center"] and best["radius"] > 0 + lo, hi = case["limits"] + assert ( + lo - 1e-6 <= min(best["entry"][0], best["exit"][0]) + and max(best["entry"][0], best["exit"][0]) <= hi + 1e-6 + ) + assert res["cases"][0]["result"]["contour"]["F"] + + +@pytest.mark.skipif(not (SUB4 / "db.json").exists(), reason="검증_sub4 스냅샷 없음") +def test_b06_areas_match_all_stations(): + """검증_sub4 전 측점 — 모델 절토 · 성토 · 토사 · 암 면적 = B06 저장값(반올림 4 자리 안).""" + rows = json.loads((SUB4 / "db.json").read_text(encoding="utf-8"))["cross_sections"] + assert rows + for row in rows: + design = json.loads(row["data"])["design"] + name = Path(row["cross_section_file_path"]).name + cross = json.loads( + (SUB4 / "B06_Section__cross_sections" / name).read_text(encoding="utf-8") + ) + m = section_model_from_b06(cross, design) + assert all(polygon_closed(r) for r in m["regions"]), name + a = m["areas"] + assert abs(a["cut_m2"] - design["cut_area_m2"]) < 2e-3, name + assert abs(a["fill_m2"] - design["fill_area_m2"]) < 2e-3, name + assert abs(a["cut_soil_m2"] - design["cut_soil_area_m2"]) < 2e-3, name + assert abs(a["cut_rock_m2"] - design["cut_rock_area_m2"]) < 2e-3, name + + +if __name__ == "__main__": + sys.exit(pytest.main([__file__, "-q"]))