From b97ddf4783344ddbe0f0a3dd55c8582a5fef4e73 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Sun, 25 Oct 2026 20:13:33 +0900 Subject: [PATCH] =?UTF-8?q?feat(533-1=20=C2=B7=20533-2):=20E04=203D=20?= =?UTF-8?q?=EA=B8=B0=EB=91=A5=20=EC=86=94=EB=B2=84=20=C2=B7=20=EA=B5=AC=20?= =?UTF-8?q?=C2=B7=20=ED=83=80=EC=9B=90=EC=B2=B4=20=ED=83=90=EC=83=89=20?= =?UTF-8?q?=C2=B7=20=ED=95=A9=EC=84=B1=20DEM=20=C2=B7=20=EC=A7=80=EC=B8=B5?= =?UTF-8?q?=20=EA=B2=A9=EC=9E=90?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Column: Hovland · Hungr 3D Bishop · Hungr 3D Janbu 간편 + 보정(f0) · 활동 방향 자동(지표 낙하) · r_u · 지하수위면 · kh · 노면 하중 · 못 뚫는 층에서 잘림 · 부분 기둥 - Search: Scoops3D 방식 구 격자(노드 = 지표 + h · 반지름 증분 · 부피/면적 제약 · 거친→고운) + 타원체(축 방향 ratio 배) · 덩이 나눔 · 셀별 최소 F · 임계 체적 래스터 - 잘린 면 ≥ 10 % → Janbu 보정 전환 · Bishop 이 낮으면 Bishop - Grid: 측점 사이 선직면 보간 δ(설계 − 원지형) · 띠 밖 원지형 · 흙 두께(암 경계 · 토사 3 m) · 리핑암 뚫림 · 발파암 못 뚫음 · 노면 하중 · 1 m 정사각 재표본 - 공용 꼴 = E03_StabSection_Types import(칸 이름 그대로) · 새 의존 없음(NumPy · SciPy) - 시험: Gens 구 1.402 · Zhang 1988 경우 1 · 2(Hungr 2.167 · 1.620) · Scoops3D 제방 1.32 · D&G 1a 1.04 · 합성 DEM 단면 · 띠 밖 · 지층 Co-Authored-By: Claude Opus 5.5 --- E04_Stab3D/E04_Stab3D_Engine_Column.py | 377 +++++++++++++++++ E04_Stab3D/E04_Stab3D_Engine_Grid.py | 273 ++++++++++++ E04_Stab3D/E04_Stab3D_Engine_Search.py | 420 +++++++++++++++++++ resources/tester/test_533_1_stab3d_grid.py | 105 +++++ resources/tester/test_533_2_stab3d_column.py | 105 +++++ resources/tester/test_533_2_stab3d_search.py | 60 +++ 6 files changed, 1340 insertions(+) create mode 100644 E04_Stab3D/E04_Stab3D_Engine_Column.py create mode 100644 E04_Stab3D/E04_Stab3D_Engine_Grid.py create mode 100644 E04_Stab3D/E04_Stab3D_Engine_Search.py create mode 100644 resources/tester/test_533_1_stab3d_grid.py create mode 100644 resources/tester/test_533_2_stab3d_column.py create mode 100644 resources/tester/test_533_2_stab3d_search.py diff --git a/E04_Stab3D/E04_Stab3D_Engine_Column.py b/E04_Stab3D/E04_Stab3D_Engine_Column.py new file mode 100644 index 000000000..2b901a17d --- /dev/null +++ b/E04_Stab3D/E04_Stab3D_Engine_Column.py @@ -0,0 +1,377 @@ +"""E04 3D 기둥 솔버 — Hovland · Hungr 3D Bishop 간편 · Hungr 3D Janbu 간편 + 보정. + +근거(수식만 · 코드는 자체): +- Reid 외(2015) USGS TM 14-A1 Scoops3D 2.1~2.5 — 기둥 기하 · 모멘트 평형 · 간극수압 · 특수 경우. +- Hungr(1987) Géotechnique 37(1) · Hungr·Salgado·Byrne(1989) CGJ 26(4) — 3D Bishop · 3D Janbu. +- Hovland(1977) ASCE GT9 — 기둥 사이 힘 0 · N = W cosψ. +- Janbu(1973) 보정 계수 f0 = 1 + b1 (d/L − 1.4 (d/L)²). +- 지식DB `original/안정성평가/05_3D_사면안정해석_논문_기술정보.md` 1장. + +꼴: 기둥은 납작한 배열 + `group`(활동체 번호) — 여러 활동체를 한 번에 셈(np.bincount). +좌표: 활동 방향 d(수평 단위) · 회전축 = 중심을 지나는 수평선 ⊥ d. +모멘트(내리미는 쪽 +) — 자중 −p·W · 지진 e·kh·W · 수직력 N·m_n · 저항 전단 −S·arm_s. + p = (P − C)·d · e = 축 높이 − 기둥 무게중심 높이 · g = 바닥면의 d 방향 기울기 + arm_s = (p·g − z_b')/√(1+g²) · m_n = cosψ (p + z_b'·g) + 구 · 원형 단면 타원체는 수직력이 축을 지나 m_n = 0 → Scoops3D 식 2.22 그대로. +Bishop 수직력(식 2.31): N = (W+Q − (cA − uA tanφ) sinα / F) / mα · mα = cosψ + sinα tanφ / F + → 저항 S_res = [c A_h + (W+Q − u A_h) tanφ] / mα (식 2.33 · A_h = 수평 면적). +Janbu 힘 평형(d 방향): F0 = Σ S_res / cosα ÷ Σ((W+Q) tanα + kh W) · 보정 F = f0 F0. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field + +import numpy as np +from scipy import ndimage + +from E03_StabSection.E03_StabSection_Types import GAMMA_W, Groundwater, Material, groundwater + +MA_MIN = 0.2 # mα 하한 경고(Duncan & Wright) — 밑이면 결과에 표시 +F_CAP = 100.0 # Scoops3D 2.5 — 너무 큰 F 는 100 +TOL = 1e-4 +MAX_ITER = 50 + + +@dataclass +class Strata: + """지층 격자 — 셀 중심 축 x(nx) · y(ny) · 지표 z_top(ny, nx) · 층 바닥(K−1, ny, nx). + + 층 k 는 위(층 k−1 바닥 · 0 층은 지표)부터 bottoms[k] 까지 · 맨 아래 층은 끝없음. + hard[k] = 활동면이 못 뚫는 층(암반) — 그 윗면에서 잘림(Search 머리 주석 규칙). + z_water = 지하수위면(mode "line" 의 3D 꼴) · q = 상재하중(kPa) 격자. + """ + + x: np.ndarray + y: np.ndarray + z_top: np.ndarray + bottoms: np.ndarray + materials: list[Material] + hard: list[bool] = field(default_factory=list) + z_water: np.ndarray | None = None + q: np.ndarray | None = None + + def __post_init__(self) -> None: + self.x = np.asarray(self.x, float) + self.y = np.asarray(self.y, float) + self.z_top = np.asarray(self.z_top, float) + self.bottoms = np.asarray(self.bottoms, float).reshape(-1, *self.z_top.shape) + if len(self.materials) != len(self.bottoms) + 1: + raise ValueError("물성 수 = 층 바닥 수 + 1 이어야 함") + self.hard = list(self.hard) or [False] * len(self.materials) + self.dx = float(self.x[1] - self.x[0]) + if abs(float(self.y[1] - self.y[0]) - self.dx) > 1e-4 * self.dx: + raise ValueError("격자 간격은 x · y 같아야 함") + # 층 바닥은 지표 아래 · 위에서 아래로 내려감 + tops = np.concatenate([self.z_top[None], self.bottoms]) + self.bottoms = np.minimum.accumulate(tops, axis=0)[1:] + self.gamma_t = np.array([m["gamma_t"] for m in self.materials]) + self.gamma_sat = np.array([m["gamma_sat"] for m in self.materials]) + self.c = np.array([m["c"] for m in self.materials]) + self.tan_phi = np.tan(np.radians([m["phi"] for m in self.materials])) + # 못 뚫는 층 윗면(없으면 −inf) · 그 기울기 + tops = np.concatenate([self.z_top[None], self.bottoms]) + first = next((k for k, h in enumerate(self.hard) if h), None) + hard_top = np.full(self.z_top.shape, -np.inf) if first is None else tops[first].copy() + self.hard_top = hard_top + finite = np.where(np.isfinite(hard_top), hard_top, 0.0) + self.hard_gy, self.hard_gx = np.gradient(finite, self.dx) + self.gy, self.gx = np.gradient(self.z_top, self.dx) + # 가장자리 · 빈칸 이웃 — 여기 닿은 활동체는 질량이 잘려 버림(Search) + edge = ndimage.binary_dilation(~np.isfinite(self.z_top)) + edge[[0, -1], :] = True + edge[:, [0, -1]] = True + self.edge = edge + self._z_pad = np.pad(self.z_top, 1, mode="edge") + self._shifted: dict = {} + + def ground_at(self, iy: np.ndarray, ix: np.ndarray, fy: float, fx: float) -> np.ndarray: + """셀 (iy, ix) 중심에서 (fy, fx) 칸만큼 옮긴 점의 지표(쌍선형 · |f| ≤ 0.5).""" + if fx == 0 and fy == 0: + return self.z_top[iy, ix] + key = (fy, fx) + if key not in self._shifted: # 옮긴 점 지표 격자를 한 번 만들어 둠 + sx, sy = int(np.sign(fx)), int(np.sign(fy)) + ax, ay = abs(fx), abs(fy) + zp = self._z_pad + ny, nx = self.z_top.shape + z00 = zp[1 : ny + 1, 1 : nx + 1] + z01 = zp[1 : ny + 1, 1 + sx : nx + 1 + sx] + z10 = zp[1 + sy : ny + 1 + sy, 1 : nx + 1] + z11 = zp[1 + sy : ny + 1 + sy, 1 + sx : nx + 1 + sx] + self._shifted[key] = (1 - ay) * ((1 - ax) * z00 + ax * z01) + ay * ( + (1 - ax) * z10 + ax * z11 + ) + return self._shifted[key][iy, ix] + + +@dataclass +class Columns: + """활동체 묶음의 기둥 — 납작한 배열(길이 n) · group 은 0..G−1.""" + + iy: np.ndarray + ix: np.ndarray + group: np.ndarray + z_slip: np.ndarray # 시험 활동면 높이(잘리기 전) + gx: np.ndarray # 활동면 기울기 ∂z/∂x + gy: np.ndarray # ∂z/∂y + center: np.ndarray # (G, 3) 회전축이 지나는 점 + direction: np.ndarray | None = None # (G, 2) 활동 방향 · None 이면 지표 낙하 방향(Scoops3D) + # 부분 기둥(셀 안 표본) — 대표 위치 · 걸친 비율 · 걸친 부분 평균 깊이 + # None = 셀 중심 · 한 칸 · 지표 − 바닥 + px: np.ndarray | None = None + py: np.ndarray | None = None + frac: np.ndarray | None = None + depth: np.ndarray | None = None + + +def _gsum(group: np.ndarray, values: np.ndarray, n: int) -> np.ndarray: + return np.bincount(group, weights=values, minlength=n) + + +def fall_direction(strata: Strata, cols: Columns, n: int) -> np.ndarray: + """활동체 지표의 평균 낙하 방향(Scoops3D 2.1) — (n, 2) 단위 벡터.""" + sx = _gsum(cols.group, -strata.gx[cols.iy, cols.ix], n) + sy = _gsum(cols.group, -strata.gy[cols.iy, cols.ix], n) + norm = np.hypot(sx, sy) + norm[norm == 0] = 1.0 + return np.column_stack([sx / norm, sy / norm]) + + +def _column_state(strata: Strata, cols: Columns, gw: Groundwater, kh: float) -> dict: + """기둥마다 자중 · 바닥 물성 · 간극수압 · 기하 — 활동면을 못 뚫는 층에서 자름.""" + iy, ix, grp = cols.iy, cols.ix, cols.group + hard_top = strata.hard_top[iy, ix] + truncated = cols.z_slip < hard_top + z_b = np.where(truncated, hard_top, cols.z_slip) + gx = np.where(truncated, strata.hard_gx[iy, ix], cols.gx) + gy = np.where(truncated, strata.hard_gy[iy, ix], cols.gy) + z_top = strata.z_top[iy, ix] if cols.depth is None else z_b + cols.depth + depth = np.maximum(z_top - z_b, 0.0) + x = strata.x[ix] if cols.px is None else cols.px + y = strata.y[iy] if cols.py is None else cols.py + + # 층별 두께 → 자중(지하수위 아래는 γsat) · 부분 기둥은 걸친 비율만큼 면적 + a_h = strata.dx**2 * (np.ones(len(iy)) if cols.frac is None else cols.frac) + k_n = len(strata.materials) + tops = np.concatenate([z_top[None], strata.bottoms[:, iy, ix]]) + bots = np.concatenate([strata.bottoms[:, iy, ix], np.full((1, len(iy)), -np.inf)]) + upper = np.minimum(tops, z_top) + lower = np.maximum(bots, z_b) + thick = np.clip(upper - lower, 0.0, None) + if gw.get("mode") == "line" and strata.z_water is not None: + zw = strata.z_water[iy, ix] + wet = np.clip(np.minimum(upper, zw) - lower, 0.0, None) + weight = a_h * (strata.gamma_t[:, None] * (thick - wet) + strata.gamma_sat[:, None] * wet) + else: + weight = a_h * strata.gamma_t[:, None] * thick + w = weight.sum(0) + + # 바닥이 놓인 층 — 경계 위면 위층(잘린 면 = 암반 위 흙의 강도) + layer = (strata.bottoms[:, iy, ix] >= z_b + 1e-9).sum(0) if k_n > 1 else np.zeros(len(iy), int) + c = strata.c[layer] + tan_phi = strata.tan_phi[layer] + + q = strata.q[iy, ix] * a_h if strata.q is not None else np.zeros_like(w) + mode = gw.get("mode", "none") + if mode == "ru": + u_ah = gw["ru"] * w # Scoops3D 식 2.37 — u A_h = ru W + elif mode == "line" and strata.z_water is not None: + u_ah = GAMMA_W * np.maximum(strata.z_water[iy, ix] - z_b, 0.0) * a_h + else: + u_ah = np.zeros_like(w) + + n_g = len(cols.center) + d = cols.direction if cols.direction is not None else fall_direction(strata, cols, n_g) + dxg, dyg = d[grp, 0], d[grp, 1] + xc, yc, zc = cols.center[grp, 0], cols.center[grp, 1], cols.center[grp, 2] + p = (x - xc) * dxg + (y - yc) * dyg + g = gx * dxg + gy * dyg + sec_g = np.sqrt(1.0 + g * g) + cos_psi = 1.0 / np.sqrt(1.0 + gx * gx + gy * gy) + zr = z_b - zc + return { + "w": w, + "q": q, + "c": c, + "tan_phi": tan_phi, + "u_ah": u_ah, + "a_h": a_h, + "cos_psi": cos_psi, + "sin_a": -g / sec_g, + "cos_a": 1.0 / sec_g, + "tan_a": -g, + "p": p, + "arm_s": (p * g - zr) / sec_g, + "m_n": cos_psi * (p + zr * g), + "e": zc - 0.5 * (z_top + z_b), + "kh": kh, + "depth": depth, + "z_b": z_b, + "z_top": z_top, + "truncated": truncated, + "direction": d, + "group": grp, + "n": n_g, + "layer": layer, + } + + +def _s_res(st: dict, f: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Bishop 수직력에서 바닥 저항력 S_res(F) · mα (Scoops3D 식 2.33 · 2.43).""" + fg = np.maximum(f[st["group"]], 1e-3) + m_a = st["cos_psi"] + st["sin_a"] * st["tan_phi"] / fg + fric = np.maximum(st["w"] + st["q"] - st["u_ah"], 0.0) * st["tan_phi"] + return (st["c"] * st["a_h"] + fric) / m_a, m_a + + +def _hovland(st: dict) -> np.ndarray: + """N = (W+Q) cosψ · 기둥 사이 힘 0 — 반복 없음(Hovland 1977 · 지식DB 05 식 ①).""" + grp, n = st["group"], st["n"] + a = st["a_h"] / st["cos_psi"] + wq = st["w"] + st["q"] + nn = wq * st["cos_psi"] + fric = np.maximum(nn - st["u_ah"] / st["cos_psi"], 0.0) * st["tan_phi"] + res = _gsum(grp, st["arm_s"] * (st["c"] * a + fric), n) + drive = _gsum(grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + nn * st["m_n"], n) + return _finish(res, drive) + + +def _finish(res: np.ndarray, drive: np.ndarray) -> np.ndarray: + out = np.full(res.shape, F_CAP) + ok = drive > 0 + out[ok] = np.minimum(np.maximum(res[ok], 0.0) / drive[ok], F_CAP) + return out + + +def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + grp, n = st["group"], st["n"] + wq = st["w"] + st["q"] + # 거꾸로 선 바닥(α<0)이 있으면 처음 값을 올림 — Scoops3D 식 2.34 · 2.35 + neg = np.where(st["sin_a"] < 0, -st["sin_a"] * st["tan_phi"] / st["cos_psi"], 0.0) + lift = np.zeros(n) + np.maximum.at(lift, grp, neg) + f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0) + lift + conv = np.zeros(n, bool) + for _ in range(MAX_ITER): + s_res, _ = _s_res(st, f) + s = s_res / np.maximum(f[grp], 1e-3) + nn = (wq - s * st["sin_a"]) / st["cos_psi"] # 연직 평형에서 N + res = _gsum(grp, st["arm_s"] * s_res, n) + drive = _gsum(grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + nn * st["m_n"], n) + f_new = _finish(res, drive) + done = np.abs(f_new - f) < TOL + f = np.where(conv, f, f_new) + conv |= done + if conv.all(): + break + _, m_a = _s_res(st, f) + ma_min = _gmin(grp, m_a, n) + # 못 쓰는 면 — 수렴 안 함 · mα ≤ 0(바닥 수직력이 음) · F ≤ 0 → nan(후보에서 빠짐) + f[~conv | (ma_min <= 0) | (f <= 0)] = np.nan + return f, ma_min + + +def _gmin(grp: np.ndarray, values: np.ndarray, n: int) -> np.ndarray: + out = np.full(n, np.inf) + np.minimum.at(out, grp, np.nan_to_num(values, nan=-np.inf)) + return out + + +def _janbu(st: dict, f0: np.ndarray) -> np.ndarray: + grp, n = st["group"], st["n"] + drive = _gsum(grp, (st["w"] + st["q"]) * st["tan_a"] + st["kh"] * st["w"], n) + f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0) + conv = np.zeros(n, bool) + for _ in range(MAX_ITER): + s_res, _ = _s_res(st, f) + f_new = _finish(_gsum(grp, s_res / st["cos_a"], n), drive) + done = np.abs(f_new - f) < TOL + f = np.where(conv, f, f_new) + conv |= done + if conv.all(): + break + f[~conv | (f <= 0)] = np.nan + return f + + +def janbu_f0(st: dict) -> np.ndarray: + """Janbu(1973) 보정 계수 — d/L 은 활동 방향 단면(현 = 윗끝 · 아랫끝 지표를 잇는 선)에서. + + b1 = 0.69(c=0) · 0.31(φ=0) · 0.50(c−φ). L = 활동 방향 길이 · d = 바닥에서 현까지 가장 먼 거리. + """ + grp, n = st["group"], st["n"] + s = st["p"] # 축 기준 d 방향 위치(같은 묶음 안 비교만 씀) + lo = np.full(n, np.inf) + hi = np.full(n, -np.inf) + np.minimum.at(lo, grp, s) + np.maximum.at(hi, grp, s) + # 양 끝 지표 높이 — 끝 기둥들 가운데 가장 얕은 바닥 위 지표 + at_lo = s <= lo[grp] + 1e-9 + at_hi = s >= hi[grp] - 1e-9 + z_lo = np.full(n, -np.inf) + z_hi = np.full(n, -np.inf) + np.maximum.at(z_lo, grp[at_lo], st["z_top"][at_lo]) + np.maximum.at(z_hi, grp[at_hi], st["z_top"][at_hi]) + length = np.maximum(hi - lo, 1e-9) + # 현: (lo, z_lo) → (hi, z_hi) · 바닥점과의 수직 거리 + dz = z_hi - z_lo + norm = np.hypot(length, dz) + dist = np.abs((s - lo[grp]) * dz[grp] - (st["z_b"] - z_lo[grp]) * length[grp]) / norm[grp] + depth = np.zeros(n) + np.maximum.at(depth, grp, dist) + ratio = depth / norm + has_c = _gsum(grp, (st["c"] > 0).astype(float), n) > 0 + has_phi = _gsum(grp, (st["tan_phi"] > 0).astype(float), n) > 0 + b1 = np.where(has_c & has_phi, 0.50, np.where(has_phi, 0.69, 0.31)) + return 1.0 + b1 * (ratio - 1.4 * ratio**2) + + +def solve( + strata: Strata, + cols: Columns, + groundwater_raw: dict | None = None, + seismic_kh: float = 0.0, +) -> dict: + """활동체마다 안전율 — hovland · bishop · janbu(간편) · janbu_corrected · 부피 · 잘림 비율. + + group 은 0..G−1 · 활동면이 지표 아래인 기둥만 셈(밖은 버림). + """ + raw = groundwater_raw or {} + # 3D 의 "line" = 지하수위면 격자(Strata.z_water) — 2D 선(line 칸)은 없어도 됨 + gw = {"mode": "line"} if raw.get("mode") == "line" and "line" not in raw else groundwater(raw) + if gw["mode"] == "line" and strata.z_water is None: + raise ValueError("지하수 mode line — Strata.z_water(지하수위면 격자)가 없음") + st = _column_state(strata, cols, gw, float(seismic_kh)) + keep = st["depth"] > 0 + n = st["n"] + st = { + k: (v[keep] if isinstance(v, np.ndarray) and v.shape == keep.shape else v) + for k, v in st.items() + } + grp = st["group"] + count = np.bincount(grp, minlength=n) + f_hov = _hovland(st) + f_bis, ma_min = _bishop(st, f_hov) + f_jan = _janbu(st, f_bis) + f0 = janbu_f0(st) + area = _gsum(grp, st["a_h"], n) + trunc_area = _gsum(grp, st["truncated"] * st["a_h"], n) + empty = count == 0 + out = { + "hovland": f_hov, + "bishop": f_bis, + "janbu": f_jan, + "janbu_corrected": f_jan * f0, + "janbu_f0": f0, + "volume": _gsum(grp, st["depth"] * st["a_h"], n), + "weight": _gsum(grp, st["w"], n), + "area": area, + "truncated_ratio": np.divide(trunc_area, area, out=np.zeros(n), where=area > 0), + "ma_min": ma_min, + "columns": count, + "direction": st["direction"], + } + for key in ("hovland", "bishop", "janbu", "janbu_corrected"): + out[key] = np.where(empty, np.nan, out[key]) + return out diff --git a/E04_Stab3D/E04_Stab3D_Engine_Grid.py b/E04_Stab3D/E04_Stab3D_Engine_Grid.py new file mode 100644 index 000000000..6c5217624 --- /dev/null +++ b/E04_Stab3D/E04_Stab3D_Engine_Grid.py @@ -0,0 +1,273 @@ +"""E04 3D 입력 격자 — 합성 DEM · 지층 경계 격자 · 노면 하중 격자. + +합성 DEM(원지형 + B06 깎기 · 쌓기): +- 측점 i · i+1 횡단선을 t 로 잇는 선직면에 셀을 앉힘 + — P = C(t) + o·L(t) 를 t 의 2차식으로 풂 · 측점 위(t=0 · 1)는 횡단선 그대로. +- δ(o) = 설계선 z − 원지형 z(같은 offset) · 측점 사이는 t 선형 · 합성 = 원지형 + δ. + → 측점 위 합성 단면 = B06 설계선 · δ=0 인 곳과 띠 밖은 원지형 그대로. +- 띠 = 두 측점 설계선 offset 범위의 합(B06 ±20 m · 사면이 안 닫히면 더 넓음) + · 한쪽 설계선 밖은 δ=0 으로 이음 · 두 구간에 걸린 셀은 |o| 작은 쪽. + +지층 경계(정함 · 보고): +- 흙 두께 T = 암 측점 |암 경계 offset|(B06 평행선) · 토사 측점 = 기본 3 m + (국내 산지 풍화토 · 붕적토 1~3 m 의 큰 쪽). +- 띠 안 = 두 측점 T 를 t 선형 · 띠 밖 · 노선 끝 밖 = 가장 가까운 측점 값. +- 암 윗면 = 원지형 − T(깎기로 드러나면 지표) · 암 종류 = 가까운 측점 라벨 + (토사 측점 아래 = 리핑암). +- 층(위→아래) = 성토 · 풍화토 · 리핑암(풍화암 · 뚫림) · 발파암(연암 · 못 뚫음). +노면 하중 = 두 측점 노면 끝(road_edges) 사이(t 보간) 띠 안 셀에 q(kPa). +""" + +from __future__ import annotations + +import numpy as np +from scipy.interpolate import RegularGridInterpolator +from scipy.spatial import cKDTree + +from E03_StabSection.E03_StabSection_Types import DEFAULT_MATERIALS, DEFAULT_ROAD_SURCHARGE_KPA +from E04_Stab3D.E04_Stab3D_Engine_Column import Strata + +DEFAULT_SOIL_THICKNESS_M = 3.0 +ROCK_KINDS = ("ripping_rock", "blasting_rock") +#: 층 → 물성 id(E03 기본표) · 못 뚫는 층 +LAYER_MATERIALS = ("fill", "weathered_soil", "weathered_rock", "soft_rock") +LAYER_HARD = (False, False, False, True) + + +def station_from_b06(cross: dict, design: dict) -> dict: + """B06 횡단 파일(`cross_NNNNNm.json`) + 설계(`cross_sections.data.design`) → 측점 한 줄.""" + ground = [[s["offset_m"], s["z"]] for s in cross["samples"] if s.get("valid", True)] + edges = design.get("road_edges") or {} + rb = design.get("rock_boundary_offset_m") + return { + "chainage_m": float(cross["chainage_m"]), + "center": [float(cross["center_x"]), float(cross["center_y"])], + "left": [float(v) for v in cross["frame"]["left_xy"]], + "ground": ground, + "design": [[p["offset_m"], p["elevation_m"]] for p in design["design_line"]], + "ground_type": design.get("ground_type", "soil"), + "rock_depth_m": abs(float(rb)) + if rb is not None and design.get("ground_type") in ROCK_KINDS + else None, + "road": sorted(float(e["offset_m"]) for e in edges.values()) if len(edges) == 2 else None, + } + + +def _delta( + st: dict, center: np.ndarray, left: np.ndarray, ground_at +) -> tuple[np.ndarray, np.ndarray]: + """설계선 offset 위 δ = 설계 − 원지형(쌍선형) — B06 이 반폭 밖으로 늘린 설계선도 받음.""" + d = np.asarray(st["design"], float) + d = d[np.argsort(d[:, 0], kind="stable")] + pts = center[None] + d[:, :1] * left[None] + zg = ground_at(pts) + delta = d[:, 1] - zg + ok = np.isfinite(delta) + return d[ok, 0], delta[ok] + + +def _segment_map(xx, yy, c0, c1, l0, l1): + """선직면 P = C0 + tD + o(L0 + tM) 에서 (t, o[m]) — t∈[0,1] 아닌 곳은 nan.""" + d = c1 - c0 + m = l1 - l0 + qx, qy = xx - c0[0], yy - c0[1] + + def cross(ax, ay, bx, by): + return ax * by - ay * bx + + a = -cross(d[0], d[1], m[0], m[1]) + b = cross(qx, qy, m[0], m[1]) - cross(d[0], d[1], l0[0], l0[1]) + c = cross(qx, qy, l0[0], l0[1]) + if abs(a) < 1e-12: + t = np.where(np.abs(b) > 1e-12, -c / np.where(b == 0, 1, b), np.nan) + else: + disc = np.sqrt(np.maximum(b * b - 4 * a * c, 0.0)) + t1 = (-b + disc) / (2 * a) + t2 = (-b - disc) / (2 * a) + ok1 = (t1 >= -1e-9) & (t1 <= 1 + 1e-9) + t = np.where(ok1, t1, t2) + t = np.where((t >= -1e-9) & (t <= 1 + 1e-9), np.clip(t, 0, 1), np.nan) + lx, ly = l0[0] + t * m[0], l0[1] + t * m[1] + norm = np.hypot(lx, ly) + o = ((qx - t * d[0]) * lx + (qy - t * d[1]) * ly) / norm + return t, o + + +def build_grid( + dem: dict, + stations: list[dict], + *, + bbox: tuple[float, float, float, float] | None = None, + default_soil_thickness_m: float = DEFAULT_SOIL_THICKNESS_M, + road_q_kpa: float = DEFAULT_ROAD_SURCHARGE_KPA, + cell_m: float | None = 1.0, +) -> dict: + """원지형 DEM(x · y · z · valid_mask) + 측점 → 합성 · 지층 격자. bbox=(x0, x1, y0, y1) 로 자름. + + cell_m — 정사각 격자로 다시 뜸(쌍선형 · B04 축은 float32 라 간격이 흔들림) + · None = DEM 축 그대로. + """ + x = np.asarray(dem["x"], float) + y = np.asarray(dem["y"], float) + z = np.asarray(dem["z"], float) + valid = np.asarray(dem.get("valid_mask", np.isfinite(z)), bool) + if x[0] > x[-1]: + x, z, valid = x[::-1], z[:, ::-1], valid[:, ::-1] + if y[0] > y[-1]: + y, z, valid = y[::-1], z[::-1], valid[::-1] + z_orig = np.where(valid, z, np.nan) + if cell_m: + x0, x1, y0, y1 = bbox or (x[0], x[-1], y[0], y[-1]) + x0, y0 = max(x0, x[0]), max(y0, y[0]) + xs = x0 + np.arange(int((min(x1, x[-1]) - x0) // cell_m) + 1) * cell_m + ys = y0 + np.arange(int((min(y1, y[-1]) - y0) // cell_m) + 1) * cell_m + f = RegularGridInterpolator((y, x), z_orig, bounds_error=False, fill_value=np.nan) + gy_, gx_ = np.meshgrid(ys, xs, indexing="ij") + z_orig = f(np.column_stack([gy_.ravel(), gx_.ravel()])).reshape(gy_.shape) + x, y = xs, ys + elif bbox is not None: + sx = (x >= bbox[0]) & (x <= bbox[1]) + sy = (y >= bbox[2]) & (y <= bbox[3]) + x, y, z_orig = x[sx], y[sy], z_orig[np.ix_(sy, sx)] + z = z_orig + xx, yy = np.meshgrid(x, y) + sts = sorted(stations, key=lambda s: s["chainage_m"]) + centers = np.array([s["center"] for s in sts], float) + lefts = np.array([s["left"] for s in sts], float) + lefts /= np.hypot(lefts[:, 0], lefts[:, 1])[:, None] + interp = RegularGridInterpolator((y, x), z_orig, bounds_error=False, fill_value=np.nan) + + def ground_at(pts: np.ndarray) -> np.ndarray: + return interp(np.column_stack([pts[:, 1], pts[:, 0]])) + + profiles = [_delta(s, centers[i], lefts[i], ground_at) for i, s in enumerate(sts)] + thick = np.array( + [ + s["rock_depth_m"] if s.get("rock_depth_m") is not None else default_soil_thickness_m + for s in sts + ] + ) + kind = np.array([1 if s.get("ground_type") == "blasting_rock" else 0 for s in sts]) + + seg = np.full(z.shape, -1) + t_map = np.full(z.shape, np.nan) + o_map = np.full(z.shape, np.nan) + dz = np.zeros(z.shape) + road = np.zeros(z.shape, bool) + for i in range(len(sts) - 1): + o0, d0 = profiles[i] + o1, d1 = profiles[i + 1] + lo, hi = min(o0[0], o1[0]), max(o0[-1], o1[-1]) + reach = max(abs(lo), abs(hi)) + corners = np.array( + [centers[i] + s * reach * lefts[i] for s in (-1, 1)] + + [centers[i + 1] + s * reach * lefts[i + 1] for s in (-1, 1)] + ) + jx = np.flatnonzero((x >= corners[:, 0].min() - 1) & (x <= corners[:, 0].max() + 1)) + jy = np.flatnonzero((y >= corners[:, 1].min() - 1) & (y <= corners[:, 1].max() + 1)) + if not len(jx) or not len(jy): + continue + sub = np.ix_(jy, jx) + t, o = _segment_map(xx[sub], yy[sub], centers[i], centers[i + 1], lefts[i], lefts[i + 1]) + inside = np.isfinite(t) & (o >= lo) & (o <= hi) + better = inside & ( + ~np.isfinite(o_map[sub]) | (np.abs(o) < np.abs(np.nan_to_num(o_map[sub], nan=np.inf))) + ) + if not better.any(): + continue + tb, ob = t[better], o[better] + val = (1 - tb) * np.interp(ob, o0, d0, 0.0, 0.0) + tb * np.interp(ob, o1, d1, 0.0, 0.0) + for arr, v in ((seg, i), (t_map, tb), (o_map, ob), (dz, val)): + blk = arr[sub] + blk[better] = v + arr[sub] = blk + r0, r1 = sts[i].get("road"), sts[i + 1].get("road") + if r0 and r1: + a = (1 - tb) * r0[0] + tb * r1[0] + b = (1 - tb) * r0[1] + tb * r1[1] + blk = road[sub] + blk[better] = (ob >= a) & (ob <= b) + road[sub] = blk + + in_band = seg >= 0 + dz = np.where(np.isfinite(z_orig), dz, 0.0) + z_syn = z_orig + dz + + # 흙 두께 · 암 종류 — 띠 안 = t 보간 · 밖 = 가까운 측점 + near = cKDTree(centers).query(np.column_stack([xx.ravel(), yy.ravel()]))[1].reshape(z.shape) + t_safe = np.nan_to_num(t_map) + s_safe = np.maximum(seg, 0) + t_band = (1 - t_safe) * thick[s_safe] + t_safe * thick[np.minimum(s_safe + 1, len(sts) - 1)] + soil_t = np.where(in_band, t_band, thick[near]) + k_band = np.where(t_safe < 0.5, kind[s_safe], kind[np.minimum(s_safe + 1, len(sts) - 1)]) + rock_kind = np.where(in_band, k_band, kind[near]) + return { + "x": x, + "y": y, + "z_orig": z_orig, + "z": z_syn, + "dz": dz, + "in_band": in_band, + "segment": seg, + "t": t_map, + "offset": o_map, + "rock_top": z_orig - soil_t, + "soil_thickness": soil_t, + "rock_kind": rock_kind, # 0 리핑암 · 1 발파암 + "q": np.where(road & in_band, road_q_kpa, 0.0), + "stations": sts, + } + + +def strata_from_grid( + grid: dict, + materials: list[dict] | None = None, + *, + water_depth_m: float | None = None, + with_surcharge: bool = True, +) -> Strata: + """격자 → 기둥 솔버 지층(층 차례 = 머리 주석) · 물성은 id 로 고름(E03 꼴 · 빠지면 기본표). + + water_depth_m — 지하수위면 = 합성 지표 − 깊이(0 = 지표까지 포화 · KDS 우기 규칙) + · None = 지하수 없음. + """ + by_id = {m["id"]: m for m in DEFAULT_MATERIALS} + by_id.update({m["id"]: m for m in materials or []}) + z = grid["z"] + fill_bot = np.minimum(z, grid["z_orig"]) + soil_bot = np.minimum(grid["rock_top"], fill_bot) + blasting = grid["rock_kind"] == 1 + rip_bot = np.where(blasting, soil_bot, -1e6) + return Strata( + grid["x"], + grid["y"], + z, + np.stack([fill_bot, soil_bot, rip_bot]), + [by_id[k] for k in LAYER_MATERIALS], + hard=list(LAYER_HARD), + z_water=None if water_depth_m is None else z - water_depth_m, + q=grid["q"] if with_surcharge else None, + ) + + +def section_check(grid: dict, dem_sampler, station: dict) -> dict: + """측점 횡단선 위 합성 DEM(쌍선형) − B06 설계선 · 원지형 − B06 지반선 — 오차 기록용.""" + d = np.asarray(station["design"], float) + c = np.asarray(station["center"], float) + left = np.asarray(station["left"], float) + left = left / np.hypot(*left) + pts = c[None] + d[:, :1] * left[None] + f = RegularGridInterpolator( + (grid["y"], grid["x"]), grid["z"], bounds_error=False, fill_value=np.nan + ) + z_syn = f(np.column_stack([pts[:, 1], pts[:, 0]])) + err = z_syn - d[:, 1] + g = np.asarray(station["ground"], float) + gp = c[None] + g[:, :1] * left[None] + z0, _ = dem_sampler(gp) + gerr = z0 - g[:, 1] + return { + "max_abs": float(np.nanmax(np.abs(err))), + "rms": float(np.sqrt(np.nanmean(err**2))), + "ground_max_abs": float(np.nanmax(np.abs(gerr))), + } diff --git a/E04_Stab3D/E04_Stab3D_Engine_Search.py b/E04_Stab3D/E04_Stab3D_Engine_Search.py new file mode 100644 index 000000000..2ce1b88be --- /dev/null +++ b/E04_Stab3D/E04_Stab3D_Engine_Search.py @@ -0,0 +1,420 @@ +"""E04 3D 활동면 탐색 — Scoops3D 방식 구 격자 + 타원체 · 셀별 최소 안전율 · 임계 체적. + +근거: Reid 외(2015) USGS TM 14-A1 3.2(탐색 격자 · 반지름 증분 · 부피/면적 제약 · +거친→고운 탐색) — 방식만 · 코드 자체. +다른 점(정하고 보고): +- 노드 높이 = 노드 아래 지표 + h — 고저차 큰 산지에서 땅속 · 너무 높은 노드를 덜어 냄. +- 타원체 = 활동 단면(d−z)은 원 · 회전축 방향 반축만 ratio 배(Hungr 1989 꼴) + — 수직력이 축을 지나 Bishop 그대로. 방향 = 노드 지표 낙하 방향(5×5 평활). + 구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1). +- 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈(ndimage.label · 8 이웃). +- 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음). + +다스리는 안전율(governing): +- 못 뚫는 층(기본 발파암 = 연암)에 잘린 바닥 면적 비 ≥ 10 % → Janbu 보정으로 전환 + · 단 Bishop 이 더 낮으면 Bishop(작은 쪽). Hungr 1989 = 약층 복합면에서 Bishop 도 + 쓸 만함 · Zhang 1988 약층 경우 Janbu 보정 1.70 > Bishop 1.62 라 전환만은 덜 보수적. +- 10 % 문턱 = 기둥 몇 개 닿아 방법이 바뀌는 지도 이음매를 막음(Zhang 약층 잘림 비 33 %). +- 그 밖 → Hungr 3D Bishop · Bishop 이 수렴 못 하면 그 후보는 버림(Scoops3D 2.3.2). +""" + +from __future__ import annotations + +import math +import time +from dataclasses import dataclass, field + +import numpy as np +from scipy import ndimage + +from E04_Stab3D.E04_Stab3D_Engine_Column import Columns, Strata, solve + +TRUNCATION_JANBU = 0.10 +_STRUCT = np.zeros((3, 3, 3), bool) +_STRUCT[1] = True # 같은 후보 안에서만 8 이웃 · 후보끼리는 안 이음 +BATCH_CELLS = 3_000_000 # 한 번에 만드는 (후보 × 창) 칸 수 — 메모리 약 0.5 GB 아래 + + +@dataclass +class SearchParams: + """탐색 조건 — 길이 m · 부피 m³ · 면적 m².""" + + box: tuple[float, float, float, float] | None = None # (x0, x1, y0, y1) · None = 격자 전체 + node_mask: np.ndarray | None = ( + None # (ny, nx) — 노드 아래 셀이 False 면 노드 뺌(노선 띠만 등 · Scoops3D 파일 탐색) + ) + spacing: float = 2.0 # 가장 고운 노드 수평 간격 + h_min: float = 2.0 # 노드 높이(지표 위) 범위 · 간격 + h_max: float = 30.0 + h_step: float = 2.0 + coarse: int = 4 # 거친 격자 배수(Scoops3D multres) + dr: float = 1.0 # 반지름 증분 + max_depth: float = 15.0 # 반지름 상한 = h + max_depth + v_min: float = 10.0 + v_max: float = 10_000.0 + a_min: float = 0.0 + a_max: float = math.inf + ratios: tuple[float, ...] = (1.0, 2.0) # 1 = 구 · 그 밖 = 축 방향 반축 배수(타원체) + sub: int = 2 # 셀 안 표본 sub×sub(부분 기둥) + refine_rounds: int = 4 # 가장 고운 간격에서 더 내려가 보는 횟수 상한 + refine_below: float = 2.0 # 고운 탐색은 이 F 아래로 고친 노드만(높은 곳 지도는 거친 값) + groundwater: dict | None = None + seismic_kh: float = 0.0 + + +@dataclass +class SearchResult: + f_min: np.ndarray # 셀별 최소 안전율(다스리는 값) · 닿지 않은 셀 = nan + volume: np.ndarray # 그 임계 활동체 부피 + method: np.ndarray # 0 = Bishop · 1 = Janbu 보정 + best: dict # 전체 최소 활동체(중심 · 반지름 · 비 · 방향 · 방법별 F · 부피) + trials: int = 0 + evaluated: int = 0 + seconds: float = 0.0 + stats: dict = field(default_factory=dict) + + +def trial_columns( + strata: Strata, + centers: np.ndarray, + r: float, + ratio: float = 1.0, + directions: np.ndarray | None = None, + v_min: float = 0.0, + v_max: float = math.inf, + a_min: float = 0.0, + a_max: float = math.inf, + sub: int = 2, +) -> tuple[Columns, np.ndarray, np.ndarray]: + """후보 T 개(같은 r · ratio)의 기둥 — (Columns, 덩이의 후보 번호, 후보별 가장 큰 덩이 부피). + + 부분 기둥 — 셀마다 sub×sub 점에서 활동면 · 지표(쌍선형)를 보고 걸친 비율 · 평균 깊이로 + 기둥을 만듦(Scoops3D 2.1 부분 기둥과 같은 뜻 · 가장자리 적분 오차 줄임). sub=1 = 셀 중심. + """ + centers = np.atleast_2d(np.asarray(centers, float)) + n_t = len(centers) + dx = strata.dx + ny, nx = strata.z_top.shape + w = int(math.ceil(r * max(ratio, 1.0) / dx)) + 1 + off = np.arange(-w, w + 1) + ix0 = np.rint((centers[:, 0] - strata.x[0]) / dx).astype(int) + iy0 = np.rint((centers[:, 1] - strata.y[0]) / dx).astype(int) + ixg = ix0[:, None, None] + off[None, None, :] + iyg = iy0[:, None, None] + off[None, :, None] + inb = (ixg >= 0) & (ixg < nx) & (iyg >= 0) & (iyg < ny) + ixc, iyc, inb = np.broadcast_arrays(np.clip(ixg, 0, nx - 1), np.clip(iyg, 0, ny - 1), inb) + shape = ixc.shape + d = None if ratio == 1.0 else np.asarray(directions, float) + k2 = 1.0 / (ratio * ratio) + + def surface(t, xs, ys): + """후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작 · 창 모두).""" + dxp = xs - centers[t, 0] + dyp = ys - centers[t, 1] + if d is None: + rem = r * r - dxp * dxp - dyp * dyp + root = np.sqrt(np.maximum(rem, 1e-12)) + return rem, centers[t, 2] - root, dxp / root, dyp / root + dxd, dyd = d[t, 0], d[t, 1] + pp = dxp * dxd + dyp * dyd + qq = -dxp * dyd + dyp * dxd # 축 방향 a = (−d_y, d_x) + rem = r * r - pp * pp - k2 * qq * qq + root = np.sqrt(np.maximum(rem, 1e-12)) + gx = (pp * dxd - k2 * qq * dyd) / root + gy = (pp * dyd + k2 * qq * dxd) / root + return rem, centers[t, 2] - root, gx, gy + + # 1) 셀 중심으로 거르기 — 걸친 셀과 그 이웃만 셀 안 표본을 봄(창 대부분은 허공 · 땅속 밖) + t_grid = np.broadcast_to(np.arange(n_t)[:, None, None], shape) + xs0 = strata.x[0] + ixg * dx + ys0 = strata.y[0] + iyg * dx + rem0, zs0, _, _ = surface(t_grid, xs0, ys0) + in0 = np.flatnonzero(inb & (rem0 > 0)) + yy0, xx0 = iyc.ravel()[in0], ixc.ravel()[in0] + zb0 = np.maximum(zs0.ravel()[in0], strata.hard_top[yy0, xx0]) + act0 = np.zeros(shape, bool) + act0.ravel()[in0[zb0 < strata.z_top[yy0, xx0]]] = True + cand = act0 if sub == 1 else ndimage.binary_dilation(act0, structure=_STRUCT) & inb + ci = np.flatnonzero(cand) + big = np.zeros(n_t) + e = np.zeros(0, int) + empty = (Columns(e, e, e, np.zeros(0), np.zeros(0), np.zeros(0), np.zeros((0, 3))), e, big) + if not len(ci): + return empty + tt = t_grid.ravel()[ci] + cy = iyc.ravel()[ci] + cx = ixc.ravel()[ci] + hard = strata.hard_top[cy, cx] + + # 2) 부분 기둥 — sub×sub 점 + n_act = np.zeros(len(ci)) + acc = {k: np.zeros(len(ci)) for k in ("depth", "zs", "gx", "gy", "px", "py")} + fr = ((np.arange(sub) + 0.5) / sub - 0.5) if sub > 1 else np.zeros(1) + for fy in fr: + for fx in fr: + xs = strata.x[0] + (cx + fx) * dx + ys = strata.y[0] + (cy + fy) * dx + rem, z_s, gx, gy = surface(tt, xs, ys) + zg = strata.ground_at(cy, cx, fy, fx) + z_b = np.maximum(z_s, hard) + act = (rem > 0) & (z_b < zg) + n_act += act + acc["depth"] += np.where(act, zg - z_b, 0.0) + for key, val in (("zs", z_s), ("gx", gx), ("gy", gy), ("px", xs), ("py", ys)): + acc[key] += np.where(act, val, 0.0) + keep = n_act > 0 + ci, tt, cy, cx, n_act = ci[keep], tt[keep], cy[keep], cx[keep], n_act[keep] + acc = {k: v[keep] for k, v in acc.items()} + active = np.zeros(shape, bool) + active.ravel()[ci] = True + labels, n_lab = ndimage.label(active, structure=_STRUCT) + if n_lab == 0: + return empty + lab = labels.ravel()[ci] + a_h = dx * dx + depth_cell = acc["depth"] / (sub * sub) + vol = np.bincount(lab, weights=depth_cell, minlength=n_lab + 1) * a_h + area = np.bincount(lab, weights=n_act / (sub * sub), minlength=n_lab + 1) * a_h + edge = np.bincount(lab, weights=strata.edge[cy, cx].astype(float), minlength=n_lab + 1) + trial_of = np.zeros(n_lab + 1, int) + trial_of[lab] = tt + np.maximum.at(big, trial_of[1:], vol[1:]) + ok = (vol >= v_min) & (vol <= v_max) & (area >= a_min) & (area <= a_max) & (edge == 0) + ok[0] = False + keep_lab = np.flatnonzero(ok) + remap = np.full(n_lab + 1, -1) + remap[keep_lab] = np.arange(len(keep_lab)) + sel = remap[lab] >= 0 + na_s = n_act[sel] + mean = {k: acc[k][sel] / na_s for k in ("zs", "gx", "gy", "px", "py")} + frac = na_s / (sub * sub) + trial = trial_of[keep_lab] + cols = Columns( + iy=cy[sel], + ix=cx[sel], + group=remap[lab][sel], + z_slip=mean["zs"], + gx=mean["gx"], + gy=mean["gy"], + center=centers[trial], + direction=None if d is None else d[trial], + px=mean["px"], + py=mean["py"], + frac=frac, + depth=depth_cell[sel] / frac, + ) + return cols, trial, big + + +def governing(res: dict) -> tuple[np.ndarray, np.ndarray]: + """다스리는 안전율 · 방법(0 Bishop · 1 Janbu 보정) — 머리 주석 규칙.""" + jc = np.where(np.isfinite(res["janbu_corrected"]), res["janbu_corrected"], np.inf) + use_j = (res["truncated_ratio"] >= TRUNCATION_JANBU) & (jc < res["bishop"]) + f = np.where(use_j, jc, res["bishop"]) + f = np.where(np.isfinite(res["bishop"]), f, np.nan) + return f, use_j.astype(np.int8) + + +def evaluate_surface( + strata: Strata, + center, + r: float, + ratio: float = 1.0, + direction=None, + groundwater: dict | None = None, + seismic_kh: float = 0.0, + sub: int = 2, +) -> dict: + """활동면 하나(구 · 타원체) — 가장 큰 덩이의 방법별 안전율 · 부피.""" + dirs = None if direction is None else np.asarray(direction, float)[None] + cols, _, _ = trial_columns(strata, np.asarray(center, float)[None], r, ratio, dirs, sub=sub) + if not len(cols.group): + raise ValueError("활동면이 지표를 자르지 않음") + res = solve(strata, cols, groundwater, seismic_kh) + k = int(np.argmax(res["volume"])) + out = {key: (v[k] if isinstance(v, np.ndarray) else v) for key, v in res.items()} + out["governing"], out["method"] = (float(x[k]) for x in governing(res)) + return out + + +def _local_fall(strata: Strata) -> tuple[np.ndarray, np.ndarray]: + zs = ndimage.uniform_filter(np.nan_to_num(strata.z_top, nan=float(np.nanmean(strata.z_top))), 5) + gy, gx = np.gradient(zs, strata.dx) + norm = np.hypot(gx, gy) + norm[norm == 0] = 1.0 + return -gx / norm, -gy / norm + + +def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: + """거친 격자 → 고운 격자(Scoops3D 3.2.2) 로 셀별 최소 안전율 · 임계 체적.""" + p = params or SearchParams() + t0 = time.perf_counter() + ny, nx = strata.z_top.shape + x0, x1, y0, y1 = p.box or (strata.x[0], strata.x[-1], strata.y[0], strata.y[-1]) + nxs = int((x1 - x0) // p.spacing) + 1 + nys = int((y1 - y0) // p.spacing) + 1 + nhs = int(round((p.h_max - p.h_min) / p.h_step)) + 1 + fall_x, fall_y = _local_fall(strata) + + f_min = np.full((ny, nx), np.inf) + volume = np.full((ny, nx), np.nan) + method = np.full((ny, nx), -1, np.int8) + node_map = np.full((ny, nx), -1, np.int64) + best = {"governing": np.inf} + seen: set[int] = set() + counters = {"trials": 0, "evaluated": 0} + + def node_xyz(ids: np.ndarray) -> np.ndarray: + ix = ids % nxs + iy = (ids // nxs) % nys + ih = ids // (nxs * nys) + x = x0 + ix * p.spacing + y = y0 + iy * p.spacing + cx = np.clip(np.rint((x - strata.x[0]) / strata.dx).astype(int), 0, nx - 1) + cy = np.clip(np.rint((y - strata.y[0]) / strata.dx).astype(int), 0, ny - 1) + z = strata.z_top[cy, cx] + p.h_min + ih * p.h_step + return np.column_stack([x, y, z]), cy, cx, p.h_min + ih * p.h_step + + def run(ids: np.ndarray) -> set[int]: + """노드들의 모든 반지름 · 비 — 셀 최소를 고친 노드 번호를 돌려줌.""" + improved: set[int] = set() + xyz0, cy0, cx0, _ = node_xyz(ids) + keep = np.isfinite(xyz0[:, 2]) + if p.node_mask is not None: + keep &= p.node_mask[cy0, cx0] + ids = ids[keep] + if not len(ids): + return improved + xyz, cy, cx, h = node_xyz(ids) + dirs = np.column_stack([fall_x[cy, cx], fall_y[cy, cx]]) + r_lo = np.maximum(p.dr, np.floor(0.5 * h / p.dr) * p.dr) + r_hi = h + p.max_depth + for ratio in p.ratios: + done = np.zeros(len(ids), bool) + for r in np.arange(r_lo.min(), r_hi.max() + 1e-9, p.dr): + live = np.flatnonzero(~done & (r_lo <= r + 1e-9) & (r <= r_hi + 1e-9)) + if not len(live): + continue + w_now = (2 * math.ceil(r * max(ratio, 1) / strata.dx) + 3) ** 2 + chunk = max(1, BATCH_CELLS // max(w_now, 1)) + for s in range(0, len(live), chunk): + part = live[s : s + chunk] + d_part = dirs[part] if ratio != 1.0 else None + cols, trial, big = trial_columns( + strata, + xyz[part], + r, + ratio, + d_part, + p.v_min, + p.v_max, + p.a_min, + p.a_max, + p.sub, + ) + done[part[big > p.v_max]] = True + counters["trials"] += len(part) + if not len(trial): + continue + res = solve(strata, cols, p.groundwater, p.seismic_kh) + counters["evaluated"] += len(trial) + f, meth = governing(res) + improved |= _update( + f, + meth, + res, + cols, + trial, + ids[part], + f_min, + volume, + method, + node_map, + p.refine_below, + ) + k = ( + int(np.nanargmin(np.where(np.isfinite(f), f, np.inf))) + if np.isfinite(f).any() + else -1 + ) + if k >= 0 and f[k] < best["governing"]: + best.clear() + best.update( + { + key: (v[k].tolist() if isinstance(v, np.ndarray) else v) + for key, v in res.items() + } + ) + best.update( + governing=float(f[k]), + method=int(meth[k]), + center=xyz[part][trial[k]].tolist(), + radius=float(r), + ratio=float(ratio), + ) + return improved + + step = max(1, int(p.coarse)) + grid = np.stack( + np.meshgrid( + np.arange(0, nhs, step), np.arange(0, nys, step), np.arange(0, nxs, step), indexing="ij" + ), + -1, + ).reshape(-1, 3) + ids = (grid[:, 0] * nys + grid[:, 1]) * nxs + grid[:, 2] + seen.update(ids.tolist()) + seeds = run(ids) + rounds = 0 + while seeds and (step > 1 or rounds < p.refine_rounds): + if step == 1: + rounds += 1 + step = max(1, step // 2) + cand = set() + for nid in seeds: + ih, rem = divmod(nid, nys * nxs) + iy, ix = divmod(rem, nxs) + for dh in (-step, 0, step): + for dy in (-step, 0, step): + for dxs in (-step, 0, step): + h2, y2, x2 = ih + dh, iy + dy, ix + dxs + if 0 <= h2 < nhs and 0 <= y2 < nys and 0 <= x2 < nxs: + cid = (h2 * nys + y2) * nxs + x2 + if cid not in seen: + cand.add(cid) + if not cand: + if step == 1: + break + continue + seen.update(cand) + seeds = run(np.fromiter(cand, np.int64)) + + f_out = np.where(np.isfinite(f_min), f_min, np.nan) + return SearchResult( + f_min=f_out, + volume=volume, + method=method, + best=best, + trials=counters["trials"], + evaluated=counters["evaluated"], + seconds=time.perf_counter() - t0, + stats={"nodes": len(seen), "cells": int(np.isfinite(f_out).sum())}, + ) + + +def _update( + f, meth, res, cols, trial, node_ids, f_min, volume, method, node_map, refine_below=np.inf +) -> set[int]: + """셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김).""" + fc = f[cols.group] + ok = np.isfinite(fc) + iy, ix, grp, fc = cols.iy[ok], cols.ix[ok], cols.group[ok], fc[ok] + order = np.argsort(-fc, kind="stable") + iy, ix, grp, fc = iy[order], ix[order], grp[order], fc[order] + better = fc < f_min[iy, ix] * (1 - 1e-6) + iy, ix, grp, fc = iy[better], ix[better], grp[better], fc[better] + f_min[iy, ix] = fc + volume[iy, ix] = res["volume"][grp] + method[iy, ix] = meth[grp] + nid = node_ids[trial[grp]] + node_map[iy, ix] = nid + return set(np.unique(nid[fc < refine_below]).tolist()) diff --git a/resources/tester/test_533_1_stab3d_grid.py b/resources/tester/test_533_1_stab3d_grid.py new file mode 100644 index 000000000..8e4ced245 --- /dev/null +++ b/resources/tester/test_533_1_stab3d_grid.py @@ -0,0 +1,105 @@ +"""533-1 E04 합성 DEM · 지층 격자 — 측점 위 단면 = B06 설계선 · 띠 밖 = 원지형 · 보간 · 지층.""" + +import numpy as np +import pytest + +from E04_Stab3D.E04_Stab3D_Engine_Grid import _segment_map, build_grid, strata_from_grid + + +def _dem(dx=0.5): + x = np.arange(-30, 90 + 1e-9, dx) + y = np.arange(-40, 40 + 1e-9, dx) + xx, yy = np.meshgrid(x, y) + z = 100.0 + 0.3 * yy + 0.02 * xx # 왼쪽(+y)이 오르막 + return {"x": x, "y": y, "z": z, "valid_mask": np.ones(z.shape, bool)}, xx, yy + + +def _station(dem, cx, depth_cut, ground_type="soil", rock=None): + """x=cx · 왼쪽 = +y · 노면 ±2 m 를 depth_cut 낮춘 설계선(사면 1:1 로 지반에 닿음).""" + offs = np.arange(-20, 20.01, 0.5) + g = 100.0 + 0.3 * offs + 0.02 * cx + design = g.copy() + road = np.abs(offs) <= 2.0 + design[road] = 100.0 + 0.02 * cx - depth_cut + for k in np.flatnonzero(~road): + edge = 2.0 * np.sign(offs[k]) + design[k] = ( + min(g[k], 100.0 + 0.02 * cx - depth_cut + abs(offs[k] - edge)) + if depth_cut > 0 + else g[k] + ) + return { + "chainage_m": float(cx), + "center": [float(cx), 0.0], + "left": [0.0, 1.0], + "ground": np.column_stack([offs, g]).tolist(), + "design": np.column_stack([offs, design]).tolist(), + "ground_type": ground_type, + "rock_depth_m": rock, + "road": [-2.0, 2.0], + } + + +def _grid(): + dem, xx, yy = _dem() + sts = [ + _station(dem, 0, 1.0), + _station(dem, 20, 3.0, "blasting_rock", 0.5), + _station(dem, 40, 3.0, "ripping_rock", 1.0), + ] + return dem, build_grid(dem, sts, cell_m=0.5), xx, yy, sts + + +def test_station_section_equals_design_and_outside_unchanged(): + dem, g, xx, yy, sts = _grid() + for st in sts: + col = np.isclose(g["x"], st["center"][0]) + d = np.asarray(st["design"]) + z_line = np.interp(g["y"], d[:, 0], d[:, 1]) + inside = np.abs(g["y"]) <= 20 + assert np.abs(g["z"][inside, col][:, None].ravel() - z_line[inside]).max() < 1e-9 + out = ~g["in_band"] + assert np.array_equal(g["z"][out], g["z_orig"][out]) + assert np.array_equal(g["z"][:, g["x"] < 0], dem["z"][:, g["x"] < 0]) # 첫 측점 앞 = 원지형 + + +def test_between_stations_linear_delta(): + _, g, xx, yy, sts = _grid() + j = np.flatnonzero(np.isclose(g["x"], 10.0))[0] + i = np.flatnonzero(np.isclose(g["y"], 0.0))[0] + # 노면 가운데 δ: 측점 0 −1 m · 측점 20 −3 m → 중간 −2 m + assert g["dz"][i, j] == pytest.approx(-2.0, abs=1e-9) + assert g["q"][i, j] == pytest.approx(10.0) + + +def test_curved_segment_map_round_trip(): + c0, c1 = np.array([0.0, 0.0]), np.array([19.0, 6.0]) + l0 = np.array([0.0, 1.0]) + a = np.radians(35.0) + l1 = np.array([-np.sin(a), np.cos(a)]) + t_true = np.array([0.0, 0.3, 0.75, 1.0]) + o_true = np.array([-15.0, 4.0, 12.0, -7.0]) + lt = l0[None] + t_true[:, None] * (l1 - l0)[None] + lt /= np.hypot(lt[:, 0], lt[:, 1])[:, None] + p = c0[None] + t_true[:, None] * (c1 - c0)[None] + o_true[:, None] * lt + t, o = _segment_map(p[:, 0], p[:, 1], c0, c1, l0, l1) + assert np.allclose(t, t_true, atol=1e-9) and np.allclose(o, o_true, atol=1e-9) + + +def test_rock_boundary_and_layers(): + _, g, xx, yy, sts = _grid() + j0 = np.flatnonzero(np.isclose(g["x"], 20.0))[0] + jm = np.flatnonzero(np.isclose(g["x"], 30.0))[0] + i = np.flatnonzero(np.isclose(g["y"], 10.0))[0] + assert g["soil_thickness"][i, j0] == pytest.approx(0.5) + assert g["soil_thickness"][i, jm] == pytest.approx(0.75) # 0.5 → 1.0 가운데 + far = np.flatnonzero(np.isclose(g["y"], 35.0))[0] + assert not g["in_band"][far, j0] + assert g["soil_thickness"][far, j0] == pytest.approx(0.5) # 띠 밖 = 가까운 측점 + assert g["soil_thickness"][i, 0] == pytest.approx(3.0) # 토사 측점 = 기본 두께 + st = strata_from_grid(g) + assert g["rock_kind"][i, j0] == 1 # 발파암 → 못 뚫는 층 윗면 = 암 윗면(깎기로 드러나면 지표) + assert st.hard_top[i, j0] == pytest.approx(min(g["rock_top"][i, j0], g["z"][i, j0])) + assert np.isneginf(st.hard_top[i, jm + 20]) or st.hard_top[i, jm + 20] < -1e5 # 리핑암 = 뚫림 + ii = np.flatnonzero(np.isclose(g["y"], 0.0))[0] + assert st.bottoms[0, ii, j0] == pytest.approx(g["z"][ii, j0]) # 깎기 = 성토 두께 0 diff --git a/resources/tester/test_533_2_stab3d_column.py b/resources/tester/test_533_2_stab3d_column.py new file mode 100644 index 000000000..b85de1337 --- /dev/null +++ b/resources/tester/test_533_2_stab3d_column.py @@ -0,0 +1,105 @@ +"""533-2 E04 3D 기둥 솔버 — 문헌 기준 문제 ±3 %. + +- Gens 외(1988) · Hungr 외(1989) · Scoops3D 6.1: 2:1 평면 사면 · 구 R=1 · 중심-지표 0.5 + · c=0.1 γ=1 φ=0 → 해석해 1.402 · 부피 0.654. +- Zhang(1988) 경우 1 · 2(PLoS ONE 18(7) e0287998 그림 5 · 표 2): 1:2 사면 H 12.2 m + · 타원체 중심 (6.1, 0, 21.3) + Rx=Rz=24.4 · Ry=78.2 · c 29 kPa φ 20° γ 18.8 → Hungr 3D Bishop 2.167 · Zhang 2.122 / + 약층(c 0 · φ 10°, 발끝 아래 0.94 m — 그림 화소로 잼) → Hungr 1.620 · Zhang 1.533. +""" + +import numpy as np +import pytest + +from E04_Stab3D.E04_Stab3D_Engine_Column import Strata +from E04_Stab3D.E04_Stab3D_Engine_Search import evaluate_surface, governing + + +def _mat(gamma, c, phi, mid="m"): + return {"id": mid, "name": mid, "gamma_t": gamma, "gamma_sat": gamma, "c": c, "phi": phi} + + +def _gens(dx=0.02, angle=0.0): + """2:1 평면 사면 — 아래로 향하는 방향을 angle(rad) 만큼 돌림.""" + beta = np.arctan(0.5) + x = np.arange(-2, 2 + 1e-9, dx) + xx, yy = np.meshgrid(x, x) + down = np.array([np.cos(angle), np.sin(angle)]) + z = -(xx * down[0] + yy * down[1]) * np.tan(beta) + st = Strata(x, x, z, np.zeros((0,) + z.shape), [_mat(1.0, 0.1, 0.0)]) + c = np.r_[0.5 * np.sin(beta) * down, 0.5 * np.cos(beta)] + return st, c + + +def test_gens_sphere_analytic(): + st, c = _gens() + r = evaluate_surface(st, c, 1.0) + assert r["bishop"] == pytest.approx(1.402, rel=0.01) + assert r["hovland"] == pytest.approx(1.402, rel=0.01) + assert r["volume"] == pytest.approx(np.pi * 0.25 * 2.5 / 3, rel=0.01) # 구 모자 부피 0.6545 + + +def test_rotated_slope_same_factor(): + """Scoops3D 6.7 대칭 — 낙하 방향이 격자 축과 어긋나도 같은 F.""" + f0 = evaluate_surface(*_gens(0.02), 1.0)["bishop"] + f1 = evaluate_surface(*_gens(0.02, np.radians(33.0)), 1.0)["bishop"] + assert f1 == pytest.approx(f0, rel=0.01) + + +def _zhang(weak: bool, dx=0.5): + x = np.arange(-15, 45 + 1e-9, dx) + y = np.arange(-90, 90 + 1e-9, dx) + xx, _ = np.meshgrid(x, y) + z = np.clip(xx * 0.5, 0, 12.2) + soil = _mat(18.8, 29.0, 20.0, "soil") + if not weak: + return Strata(x, y, z, np.zeros((0,) + z.shape), [soil]) + zw = -0.94 + bottoms = np.stack([np.full(z.shape, zw + 0.3), np.full(z.shape, zw)]) + mats = [soil, _mat(18.8, 0.0, 10.0, "weak"), _mat(18.8, 1000.0, 45.0, "hard")] + return Strata(x, y, z, bottoms, mats, hard=[False, False, True]) + + +def test_zhang_case1_ellipsoid(): + r = evaluate_surface(_zhang(False), [6.1, 0.0, 21.3], 24.4, 78.2 / 24.4, [-1.0, 0.0]) + assert r["bishop"] == pytest.approx(2.167, rel=0.03) # Hungr 외 1989 + assert r["bishop"] == pytest.approx(2.122, rel=0.03) # Zhang 1988 + assert r["hovland"] < r["bishop"] # Hovland 은 보수적 + + +def test_zhang_case2_weak_layer_truncated(): + r = evaluate_surface(_zhang(True), [6.1, 0.0, 21.3], 24.4, 78.2 / 24.4, [-1.0, 0.0]) + assert r["bishop"] == pytest.approx(1.620, rel=0.03) # Hungr 외 1989 + assert r["truncated_ratio"] > 0.1 + # 잘린 면 → Janbu 보정으로 전환하되 Bishop 이 낮으면 Bishop(작은 쪽) + assert r["governing"] == pytest.approx(min(r["bishop"], r["janbu_corrected"])) + assert 1.45 < r["janbu_corrected"] < 1.80 # 문헌 폭 1.53~1.67 둘레 + + +def test_pore_pressure_and_seismic(): + st, c = _gens(0.04) + dry = evaluate_surface(st, c, 1.0) + # φ=0 이면 간극수압은 저항에 안 들어감 + assert evaluate_surface(st, c, 1.0, groundwater={"mode": "ru", "ru": 0.3})[ + "bishop" + ] == pytest.approx(dry["bishop"]) + # c−φ 흙: r_u 가 F 를 낮춤 · 깊은 지하수위는 영향 없음 · kh 는 낮춤 + st2 = Strata(st.x, st.y, st.z_top, st.bottoms, [_mat(1.0, 0.05, 25.0)], z_water=st.z_top - 5.0) + base = evaluate_surface(st2, c, 1.0)["bishop"] + assert evaluate_surface(st2, c, 1.0, groundwater={"mode": "ru", "ru": 0.3})["bishop"] < base + assert evaluate_surface(st2, c, 1.0, groundwater={"mode": "line"})["bishop"] == pytest.approx( + base + ) + st2.z_water = st2.z_top + 0.0 + assert evaluate_surface(st2, c, 1.0, groundwater={"mode": "line"})["bishop"] < base + assert evaluate_surface(st2, c, 1.0, seismic_kh=0.1)["bishop"] < base + + +def test_governing_picks_bishop_when_not_truncated(): + res = { + "bishop": np.array([1.5, 1.5]), + "janbu_corrected": np.array([1.4, 1.4]), + "truncated_ratio": np.array([0.0, 0.5]), + } + f, meth = governing(res) + assert f.tolist() == [1.5, 1.4] and meth.tolist() == [0, 1] diff --git a/resources/tester/test_533_2_stab3d_search.py b/resources/tester/test_533_2_stab3d_search.py new file mode 100644 index 000000000..1ef54b604 --- /dev/null +++ b/resources/tester/test_533_2_stab3d_search.py @@ -0,0 +1,60 @@ +"""533-2 E04 활동면 탐색 — Scoops3D 검증 표 전역 최소 ±3 %. + +Reid 외(2015) Scoops3D 표 6.8 · 6.11(CLARA-W 와 같은 구): H 10 m 제방 +- 30° · λ 0.125(γ 19 · φ 20° · c 8.64 kPa) → Bishop 1.32 · 부피 1,124 m³ +- Donald & Giam 1a 2:1 · γ 20 · φ 19.6° · c 3 kPa → 1.04 +탐색 격자는 시험 시간 때문에 문헌 임계 중심 둘레로 좁힘. +""" + +import numpy as np +import pytest + +from E04_Stab3D.E04_Stab3D_Engine_Column import Strata +from E04_Stab3D.E04_Stab3D_Engine_Search import SearchParams, search + + +def _embankment(run, gamma, phi, c, dx=1.0): + x = np.arange(-25, run + 35 + 1e-9, dx) + y = np.arange(-40, 40 + 1e-9, dx) + xx, _ = np.meshgrid(x, y) + z = np.clip(xx / run * 10.0, 0, 10.0) + mat = {"id": "m", "name": "m", "gamma_t": gamma, "gamma_sat": gamma, "c": c, "phi": phi} + return Strata(x, y, z, np.zeros((0,) + z.shape), [mat]) + + +def test_embankment_30deg_global_minimum(): + st = _embankment(10 / np.tan(np.radians(30)), 19, 20, 8.64) + p = SearchParams( + box=(-2, 10, -1, 1), + h_min=4, + h_max=24, + spacing=1.0, + coarse=2, + max_depth=20, + v_min=50, + v_max=2e4, + ratios=(1.0,), + ) + r = search(st, p) + assert r.best["governing"] == pytest.approx(1.32, rel=0.03) + assert r.best["volume"] == pytest.approx(1124, rel=0.15) + f = r.f_min + assert np.nanmin(f) == pytest.approx(r.best["governing"]) + # 셀별 최소 지도: 임계 활동체 안 셀은 전역 최소 · 부피 래스터가 같이 채워짐 + assert np.isfinite(r.volume[np.isfinite(f)]).all() + + +def test_donald_giam_1a(): + st = _embankment(20, 20, 19.6, 3.0) + p = SearchParams( + box=(-4, 6, -1, 1), + h_min=16, + h_max=30, + spacing=1.0, + coarse=2, + max_depth=20, + v_min=50, + v_max=3e4, + ratios=(1.0,), + ) + assert search(st, p).best["governing"] == pytest.approx(1.04, rel=0.03)