eomsangdon/Aislo
1
0
Fork 0
Code Issues Pull Requests 2 Actions Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/sub_laptop_5' into main_laptop_1

Browse Source
This commit is contained in:
eomsangdon
2026-10-25 20:15:38 +09:00
parent 770ef96f7f b97ddf4783
commit 1d7b1bc793
6 changed files with 1340 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
E04_Stab3D/E04_Stab3D_Engine_Column.py
+377
View File
@@ -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
E04_Stab3D/E04_Stab3D_Engine_Grid.py
+273
View File
@@ -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))),
}
E04_Stab3D/E04_Stab3D_Engine_Search.py
+420
View File
@@ -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())
resources/tester/test_533_1_stab3d_grid.py
+105
View File
@@ -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
resources/tester/test_533_2_stab3d_column.py
+105
View File
@@ -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]
resources/tester/test_533_2_stab3d_search.py
+60
View File
@@ -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)