eomsangdonandClaude Opus 5.5
bfd87f3e3e
- 모델: 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>
509 lines
22 KiB
Python
509 lines
22 KiB
Python
"""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
|