Files
Aislo/E03_StabSection/E03_StabSection_Engine_Describe.py
T

115 lines
4.8 KiB
Python

"""E03 결과 꾸미기 — 임계면 글상자 · 절편 자유물체도 · F–λ 곡선 · Safety Map 상위 면.
좌표는 원래 방향(x = B06 offset · 좌 +)으로 되돌려 냄. 출력 칸은 덧붙이기만(E03 화면 · E05 보고서가 읽음).
"""
import numpy as np
from E03_StabSection import E03_StabSection_Engine_GLE as gle
from E03_StabSection import E03_StabSection_Engine_LEM as lem
from E03_StabSection.E03_StabSection_Engine_Surface import (
circle_slices,
evaluate,
poly_slices,
three_point_center,
)
#: Safety Map 면 수 · 면마다 꼭짓점 수(원호)
SAFETY_MAP_TOP = 50
_MAP_POINTS = 25
def _curve(sl, center, method, f=None) -> dict:
"""F_m · F_f 대 λ (화면 그래프) — Spencer = 상수 f · M-P = 반정현."""
if method == "spencer":
grid, fm, ff = lem.lambda_curve(
sl, center, np.linspace(*lem.LAMBDA_RANGE, gle.CURVE_POINTS)
)
fm, ff = fm[:, 0], ff[:, 0]
else:
grid, fm, ff = gle.lambda_curves(sl, center, f)
fm, ff = fm[0], ff[0]
num = lambda v: [float(x) if np.isfinite(x) else None for x in v] # noqa: E731
return {"lambda": num(grid), "Fm": num(fm), "Ff": num(ff)}
def describe(prep, kind, geom, method, n, min_depth, seed_F=None) -> dict:
"""임계면 하나 → 안전율 · 좌표 · 절편표(자유물체도 힘) · 수렴 · (Spencer · M-P) λ · F–λ 곡선."""
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]])
pts = np.column_stack([sl["xe"][0], 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])
seed = None if seed_F is None else np.array([seed_F], float)
res = evaluate(sl, center, (method,), mp_seed=seed)[method]
F = float(res["F"][0])
lam = float(res["lambda"][0]) if method in ("spencer", "morgenstern_price") else 0.0
f = gle.interslice_function(sl["xe"]) if method == "morgenstern_price" else None
rows = lem.slice_table(sl, 0, F, lam, method, f) if np.isfinite(F) else []
for r in rows:
r["material"] = (
prep.materials[r["material"]]["id"] if r["material"] < len(prep.materials) else None
)
# 행 순서 = 머리 → 끝 · E · X 의 left = 머리 쪽 경계 · right = 끝 쪽 경계(해석 틀 그대로)
r["x_left"], r["x_right"] = sorted((s * r["x_left"], s * r["x_right"]))
r["x_mid"] = s * r["x_mid"]
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])
else:
out["moment_center"] = [float(s * center[0, 0]), float(center[0, 1])]
if method in ("spencer", "morgenstern_price"):
out["lambda"] = lam
out["theta_deg"] = float(np.degrees(np.arctan(lam)))
if np.isfinite(F):
out["lambda_curve"] = _curve(sl, center, method, f)
if method == "morgenstern_price":
out["interslice_function"] = "half_sine"
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 safety_map(prep, circles, F, method, top=SAFETY_MAP_TOP) -> dict:
"""Safety Map — 거친 탐색 원호 중 안전율 낮은 top 개(F · 중심 · 반지름 · 진입 · 진출 · 꼭짓점)."""
s = prep.direction
xc, zc, R, xE, xX = circles
order = np.argsort(np.where(np.isfinite(F), F, np.inf))[:top]
order = order[np.isfinite(F[order])]
t = np.linspace(0.0, 1.0, _MAP_POINTS)
out = []
for k in order:
x = xE[k] + (xX[k] - xE[k]) * t
z = zc[k] - np.sqrt(np.clip(R[k] ** 2 - (x - xc[k]) ** 2, 0, None))
pts = np.column_stack([s * x, z])
out.append(
{
"F": float(F[k]),
"center": [float(s * xc[k]), float(zc[k])],
"radius": float(R[k]),
"entry": [float(s * xE[k]), float(z[0])],
"exit": [float(s * xX[k]), float(z[-1])],
"points": [[float(a), float(b)] for a, b in (pts if s > 0 else pts[::-1])],
}
)
return {"method": method, "surfaces": out}