eomsangdonandClaude Opus 5.5
f341e0465f
- LEM 절편: 바닥 물성 = 벽 인 절편이 하나라도 있으면 그 면 무효 · 벽 밑으로 도는 면만 후보 - 시험: 벽 한가운데 지나는 꺾은선 = 값 없음 · 벽 바닥 아래 꺾은선 = 값 있음 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
623 lines
27 KiB
Python
623 lines
27 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
|
||
#: Janbu 보정 f0 식의 d/L 상한(지식DB 02 §2.3)
|
||
JANBU_DL_MAX = 0.35
|
||
MAX_ITER = 100
|
||
WALL_ID = "wall" # Engine_Model.WALL_MATERIAL id
|
||
|
||
|
||
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])
|
||
# 벽(옹벽 · 석축) — 바닥이 벽 안을 지나는 활동면은 버림(벽 밑 · 뒤로 돌아가는 면만 후보)
|
||
self.blocked = np.array([m["id"] == WALL_ID for m in mats] + [False])
|
||
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))
|
||
# 절편 계산이 한 번에 모아 읽는 표 (띠, 6, 구간) — 바닥 · 바닥 기울기 · 윗면 · 윗면 기울기 · γt · γsat − γt
|
||
gt, gs = self.gamma_t[self.mat], self.gamma_sat[self.mat]
|
||
self.table = np.stack(
|
||
[self.bot_a, self.bot_b - self.bot_a, self.top_a, self.top_b - self.top_a, gt, gs - gt],
|
||
axis=1,
|
||
)
|
||
|
||
# ---- 조회 -------------------------------------------------------------
|
||
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))
|
||
j = np.clip(np.searchsorted(prep.X, xm, side="right") - 1, 0, len(prep.xa) - 1)
|
||
t = ((xm - prep.xa[j]) / (prep.xb[j] - prep.xa[j]))[..., None]
|
||
G = prep.table[j]
|
||
mat = prep.mat[j]
|
||
bot = G[..., 0, :] + t * G[..., 1, :]
|
||
top = G[..., 2, :] + t * G[..., 3, :]
|
||
gt, dgs = G[..., 4, :], G[..., 5, :]
|
||
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)
|
||
u = GAMMA_W * np.clip(zw - zb, 0.0, None)
|
||
weight_per_b = np.sum(gt * thick + dgs * wet, axis=-1)
|
||
else:
|
||
wet = None
|
||
u = np.zeros_like(zb)
|
||
weight_per_b = np.sum(gt * thick, axis=-1)
|
||
if prep.kh > 0: # 무게중심 높이는 지진 관성력 팔에만 쓰임
|
||
if wet is None:
|
||
moment_z = np.sum(gt * thick * (lo + 0.5 * thick), axis=-1)
|
||
else:
|
||
dry = thick - wet
|
||
moment_z = np.sum(
|
||
gt * dry * (lo + wet + 0.5 * dry) + (gt + dgs) * wet * (lo + 0.5 * wet), axis=-1
|
||
)
|
||
zg = np.where(weight_per_b > 0, moment_z / np.maximum(weight_per_b, _TINY), zb)
|
||
else:
|
||
zg = 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)
|
||
valid &= ~np.any(prep.blocked[base_mat] & found, 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
|
||
|
||
|
||
# ---- 방법 -------------------------------------------------------------------
|
||
# 한 벌 식을 F 와 무관한 칸(_pre) · 모멘트 중심 칸(_mom) 으로 미리 묶어 두고, F · λ 마다는 N 만 새로 셈.
|
||
# g = 1/F 로 놓으면 N = (P0 − g·P1)/(Q0 + g·Q1) — dN/dg · dN/dλ 가 닫힌 꼴이라 뉴턴 야코비안을 해석적으로 씀.
|
||
def _pre(sl: dict) -> dict:
|
||
p = sl.get("_pre")
|
||
if p is None:
|
||
t, sa, ca = sl["tanphi"], sl["sa"], sl["ca"]
|
||
cu = sl["c"] * sl["l"] - sl["u"] * sl["l"] * t
|
||
khws = sl["kh"] * sl["Ws"]
|
||
lw = khws.copy()
|
||
lw[:, 0] += sl["Pw"] # 머리 경계 수압은 전단 없음 → X_0 = 0
|
||
p = {
|
||
"cu": cu,
|
||
"cusa": cu * sa,
|
||
"cuca": cu * ca,
|
||
"sat": sa * t,
|
||
"tca": t * ca,
|
||
"khws": khws,
|
||
"lw": lw,
|
||
"fd0": khws.sum(1) + sl["Pw"],
|
||
"rc0": (cu * ca).sum(1),
|
||
}
|
||
sl["_pre"] = p
|
||
return p
|
||
|
||
|
||
def _mom(sl: dict, center: np.ndarray):
|
||
"""모멘트 팔(a) · 분모 N 계수(gN) · 상수 칸 — 같은 중심 배열이면 다시 쓰고 아니면 새로 셈."""
|
||
cached = sl.get("_mom")
|
||
if cached is not None and cached[0] is center:
|
||
return cached[1]
|
||
p = _pre(sl)
|
||
xo, zo = center[:, 0:1], center[:, 1:2]
|
||
dx, dz = sl["xm"] - xo, sl["zb"] - zo
|
||
a = dx * sl["sa"] + dz * sl["ca"]
|
||
gn = dz * sl["sa"] - dx * sl["ca"]
|
||
d0 = ((sl["xw"] - xo) * sl["W"] + (sl["zg"] - zo) * p["khws"]).sum(1)
|
||
d0 = d0 + (sl["zp"] - center[:, 1]) * sl["Pw"]
|
||
out = (a * sl["tanphi"], gn, d0, (a * p["cu"]).sum(1))
|
||
sl["_mom"] = (center, out)
|
||
return out
|
||
|
||
|
||
def _normal(sl, F, lam, grad=0):
|
||
"""절편 수직력 N · 분모(m_α − λA) · grad ≥ 1 dN/dg · grad 2 dN/dλ — F · λ 꼴 (S,)."""
|
||
p = _pre(sl)
|
||
g = 1.0 / np.asarray(F, float)[:, None]
|
||
L = np.asarray(lam, float)[:, None]
|
||
P1 = p["cusa"] + L * p["cuca"]
|
||
Q1 = p["sat"] + L * p["tca"]
|
||
den = sl["ca"] - L * sl["sa"] + g * Q1
|
||
safe = np.where(np.abs(den) < _TINY, _TINY, den)
|
||
N = (sl["W"] + L * p["lw"] - g * P1) / safe
|
||
if not grad:
|
||
return N, den
|
||
dg = -(P1 + N * Q1) / safe
|
||
if grad == 1:
|
||
return N, den, dg, None
|
||
dl = (p["lw"] - g * p["cuca"] + N * (sl["sa"] - g * p["tca"])) / safe
|
||
return N, den, dg, dl
|
||
|
||
|
||
def _parts(sl: dict, F, lam):
|
||
"""주어진 F · λ 의 N · 저항 R · 분모 · A · B0(절편표용).
|
||
|
||
유효 수직력 음(인장) 절편도 그대로 둠(Slide2 와 같음 — 0 으로 자르면 ACADS · Arai & Tagyo
|
||
대조가 3~6 % 벗어남). 대신 저항 합이 0 이하인 면은 `_finish` 가 버림.
|
||
"""
|
||
p = _pre(sl)
|
||
F = np.asarray(F, float)
|
||
N, den = _normal(sl, F, lam)
|
||
g = 1.0 / F[:, None]
|
||
A = sl["sa"] - p["tca"] * g
|
||
B0 = p["khws"] - p["cuca"] * g
|
||
return N, p["cu"] + N * sl["tanphi"], den, A, B0
|
||
|
||
|
||
def _ratio(num, den, dnum=None, dden=None, neg=True):
|
||
"""f = num/den (작은 분모는 ∓TINY) · 도함수 있으면 (f, df)."""
|
||
d = np.where(np.abs(den) < _TINY, -_TINY if neg else _TINY, den)
|
||
f = num / d
|
||
return f if dnum is None else (f, (dnum - f * dden) / d)
|
||
|
||
|
||
def _balance(cf, cd, c0, d0, nd, F, neg):
|
||
"""f = (c0 + Σcf·N)/(d0 + Σcd·N) 와 df/dF · df/dλ (nd = _normal 결과 · 없는 도함수는 None)."""
|
||
N, _den, dg, dl = nd
|
||
num, den = c0 + (cf * N).sum(1), d0 + (cd * N).sum(1)
|
||
if dg is None:
|
||
return _ratio(num, den, neg=neg), None, None
|
||
f, fg = _ratio(num, den, (cf * dg).sum(1), (cd * dg).sum(1), neg)
|
||
fl = None if dl is None else _ratio(num, den, (cf * dl).sum(1), (cd * dl).sum(1), neg)[1]
|
||
return f, -fg / np.asarray(F, float) ** 2, fl
|
||
|
||
|
||
def _fm(sl, center, F, lam, grad=0, nd=None):
|
||
"""모멘트 평형 F_m = Σ a·R / D → (F_m, dF_m/dF, dF_m/dλ, 분모)."""
|
||
at, gn, d0, ma = _mom(sl, center)
|
||
nd = nd or _normal(sl, F, lam, grad)
|
||
if len(nd) == 2:
|
||
nd = (*nd, None, None)
|
||
return (*_balance(at, gn, ma, d0, nd, F, True), nd[1])
|
||
|
||
|
||
def _ff(sl, F, lam, grad=0, nd=None):
|
||
"""힘 평형 F_f = Σ R·cosα / (Σ(N·sinα + kh·Ws) + Pw) → (F_f, dF_f/dF, dF_f/dλ, 분모)."""
|
||
p = _pre(sl)
|
||
nd = nd or _normal(sl, F, lam, grad)
|
||
if len(nd) == 2:
|
||
nd = (*nd, None, None)
|
||
return (*_balance(p["tca"], sl["sa"], p["rc0"], p["fd0"], nd, F, False), nd[1])
|
||
|
||
|
||
def _iterate(fun, F0):
|
||
"""F = fun(F) 풀이(묶음 전체) — 해석 도함수 뉴턴 · 안 되면 반쯤 섞은 고정점 걸음.
|
||
|
||
fun(F) → (f, df/dF, …). 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남.
|
||
"""
|
||
F = np.asarray(F0, float).copy()
|
||
done = np.zeros(len(F), bool)
|
||
iters = np.zeros(len(F), int)
|
||
for _ in range(MAX_ITER):
|
||
f1, df = fun(F)[:2]
|
||
h1 = f1 - F
|
||
dh = df - 1.0
|
||
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
|
||
return F, iters, done & np.isfinite(F)
|
||
|
||
|
||
def ordinary(sl: dict, center: np.ndarray) -> dict:
|
||
"""Fellenius — 절편간 힘 무시 · 바닥 수직 분해 N = W·cosα − kh·Ws·sinα · 모멘트 평형."""
|
||
N = sl["W"] * sl["ca"] - _pre(sl)["khws"] * sl["sa"]
|
||
at, gn, d0, ma = _mom(sl, center)
|
||
F = _ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1))
|
||
ok = (
|
||
sl["valid"]
|
||
& np.isfinite(F)
|
||
& (F > F_MIN)
|
||
& ((_pre(sl)["cu"] + N * sl["tanphi"]).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, 1), 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, 1), 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
|
||
)
|
||
# d/L ≤ 0.35 에서만 식이 맞음(지식DB 02 §2.3) — 넘으면 0.35 값(f0 곡선 꼭대기 근처)
|
||
ratio = np.minimum(dist.max(1) / np.maximum(L, _TINY), JANBU_DL_MAX)
|
||
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 _spencer_newton(sl, center, F, lam, iters=40):
|
||
"""(F_m − F, F_f − F) = 0 뉴턴 — 걸음 크기 제한(F 절반 · λ 0.3) · 수렴 = 두 잔차 < 1e-4 · λ 범위 안.
|
||
|
||
끝난 면은 빼고(남은 면이 절반 아래면 묶음을 줄임) · F < F_MIN 으로 빠진 면은 버림(어차피 `_finish` 가 버림).
|
||
"""
|
||
F, lam = F.copy(), lam.copy()
|
||
it = np.zeros(len(F), int)
|
||
done = np.zeros(len(F), bool)
|
||
idx, w, wc = np.arange(len(F)), sl, center
|
||
for _ in range(iters):
|
||
act = ~done[idx]
|
||
if not act.any():
|
||
break
|
||
if act.sum() < 0.5 * len(idx):
|
||
idx, w, wc = idx[act], _subset(w, act), wc[act]
|
||
f, lm, dn = F[idx], lam[idx], done[idx]
|
||
nd = _normal(w, f, lm, 2)
|
||
fm, fmF, fml, _ = _fm(w, wc, f, lm, nd=nd)
|
||
ff, ffF, ffl, _ = _ff(w, f, lm, nd=nd)
|
||
rm, rf = fm - f, ff - f
|
||
j11, j12, j21, j22 = fmF - 1.0, fml, ffF - 1.0, ffl
|
||
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(dn, f, f + scale * dF)
|
||
F[idx], lam[idx] = f, np.where(dn, lm, lm + scale * dl)
|
||
it[idx] += ~dn
|
||
small = (np.abs(dF) < TOL * np.maximum(1, np.abs(f))) & (np.abs(dl) < TOL)
|
||
done[idx] = dn | small | ~np.isfinite(f) | (f < F_MIN)
|
||
nd = _normal(sl, F, lam)
|
||
rm = _fm(sl, center, F, lam, nd=nd)[0] - F
|
||
rf = _ff(sl, F, lam, nd=nd)[0] - F
|
||
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 lambda_curve(sl, center, grid=None):
|
||
"""λ 마다 모멘트 · 힘 평형 F_m(λ) · F_f(λ) — 꼴 (len(grid), S) 둘(화면 F–λ 그래프 · 훑기)."""
|
||
grid = np.linspace(LAMBDA_RANGE[0], LAMBDA_RANGE[1], 16) if grid is None else grid
|
||
S = len(sl["W"])
|
||
fm_all = np.full((len(grid), S), np.nan)
|
||
ff_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_all[k] = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0)[0]
|
||
ff_all[k] = _iterate(lambda f: _ff(sl, f, lam, 1), F0)[0]
|
||
return grid, fm_all, ff_all
|
||
|
||
|
||
def _spencer_scan(sl, center):
|
||
"""λ 를 훑어 F_m(λ) − F_f(λ) 부호가 바뀌는 곳(λ=0 가까운 쪽)을 출발점으로."""
|
||
grid, fm_all, ff_all = lambda_curve(sl, center)
|
||
S = len(sl["W"])
|
||
diff = fm_all - ff_all
|
||
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:
|
||
"""면 묶음 일부 — 미리 묶은 칸(_pre · _mom)은 버림(부분 꼴이 달라 다시 셈)."""
|
||
return {
|
||
k: (v[mask] if isinstance(v, np.ndarray) and v.shape[:1] == mask.shape else v)
|
||
for k, v in sl.items()
|
||
if k not in ("_pre", "_mom")
|
||
}
|
||
|
||
|
||
def _finish(sl, F, it, conv, lam):
|
||
Fs = np.where(np.isfinite(F) & (F > 0), F, 1.0)
|
||
N, den = _normal(sl, Fs, lam)
|
||
R = _pre(sl)["cu"] + N * sl["tanphi"]
|
||
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", f=None
|
||
) -> list[dict]:
|
||
"""한 면(k 번째)의 절편표 + 자유물체도 힘 — 행 = 머리 → 끝 · left = 머리 쪽 경계.
|
||
|
||
폭 · 높이 · 경사 · 무게 · 수압 · 강도 · N(전) · N' · 동원 전단 · 저항 · 지진력 · 경계 E · X.
|
||
M-P 는 f(경계마다 절편간 힘 함수)를 받아 차례 계산(`Engine_GLE.forces`).
|
||
"""
|
||
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))
|
||
X = E
|
||
elif method == "morgenstern_price":
|
||
from E03_StabSection.E03_StabSection_Engine_GLE import forces
|
||
|
||
N, _den, E, X = forces(one, np.array([F]), np.array([lam]), f[k : k + 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)
|
||
X = lam * E
|
||
X[:, 0] = 0.0
|
||
R = one["c"] * one["l"] + (N - one["u"] * one["l"]) * one["tanphi"]
|
||
rows = []
|
||
for i in range(N.shape[1]):
|
||
ul = float(one["u"][0, i] * one["l"][0, i])
|
||
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]) - ul,
|
||
"shear_mobilized": float(R[0, i] / F),
|
||
"E_left": float(E[0, i]),
|
||
"X_left": float(X[0, i]),
|
||
# 덧(2차) — 자유물체도
|
||
"x_mid": float(one["xm"][0, i]),
|
||
"z_base": float(one["zb"][0, i]),
|
||
"z_top": float(one["zt"][0, i]),
|
||
"normal_total": float(N[0, i]),
|
||
"pore_force": ul,
|
||
"shear_strength": float(R[0, i]),
|
||
"seismic_force": float(one["kh"] * one["Ws"][0, i]),
|
||
"E_right": float(E[0, i + 1]),
|
||
"X_right": float(X[0, i + 1]),
|
||
}
|
||
)
|
||
return rows
|