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