Aislo/E03_StabSection/E03_StabSection_WallCheck.py

"""E03 벽(옹벽 · 석축) 외적 안정 — 활동 · 전도 · 지지력 (요구 345 · 351 · PLAN 536-3).

storage 를 모르는 순수 함수. 벽 폴리곤(2D 모형 `walls` · x = B06 offset · z = 표고) + 배면 · 기초 물성
→ 하중 · 안전율 · 판정. 기준 = KDS 11 80 05(지식DB `구조물_안정계산.md` §2 · `original/안정성평가/04`).

- 토압 = Coulomb(벽 배면 = 뒷굽 바닥 ~ 배면 꼭대기 가상면 · 벽마찰 δ) 또는 Rankine(뒷굽 연직 가상면 ·
  δ = β) · 지진 시 = Mononobe-Okabe(정적 H/3 + 증분 0.6H · Seed-Whitman) + 벽 관성 kh·W.
- 가상면과 벽 사이 흙 = 벽 무게에 더함 · 배면 점착 = 0(보수 쪽).
- 상재하중 q = 배면 등분포(노면 하중) · 수압 = 배면 수위 hw(바닥부터) 정수압 + 바닥 삼각 양압력 ·
  수위 아래 토압 = 수중 단위중량.
- 활동 = (V·μ + c_b·B) / ΣH(수동토압 빼고) · 전도 = 앞굽 기준 저항 / 전도 모멘트 ·
  지지력 = Meyerhof 극한(유효폭 B' · 경사 계수) / 최대 지반반력(사다리꼴 · 삼각).
"""

from __future__ import annotations

import math
from typing import Any

import numpy as np
from shapely.geometry import LineString, Polygon, box
from shapely.ops import unary_union

GAMMA_W = 9.81
#: KDS 11 80 05 기준안전율 — 평상시 · 지진 시(활동은 수동토압 뺀 값)
REQUIRED = {
    "normal": {"sliding": 1.5, "overturning": 2.0, "bearing": 3.0},
    "seismic": {"sliding": 1.2, "overturning": 1.5, "bearing": 2.0},
}
#: 벽 단위중량(kN/m³) — 콘크리트 · 찰쌓기 23 · 메쌓기(공극) 20
WALL_GAMMA = {"default": 23.0, "dry_masonry": 20.0}
_N = 400  # 토압 적분 칸


def ka_coulomb(phi: float, delta: float, theta: float, beta: float, psi: float = 0.0) -> float:
    """Coulomb 주동토압계수(ψ > 0 = Mononobe-Okabe KAE · kv = 0). 각 = 도.

    θ = 배면이 연직에서 앞쪽(벽 쪽)으로 누운 각(쐐기가 커지는 쪽 +) · β = 배면 흙 경사 · δ = 벽마찰."""
    f, d, t, b, p = (math.radians(v) for v in (phi, delta, theta, beta, psi))
    under = math.sin(f + d) * math.sin(f - b - p) / (math.cos(d + t + p) * math.cos(b - t))
    root = math.sqrt(max(under, 0.0))
    return math.cos(f - t - p) ** 2 / (
        math.cos(p) * math.cos(t) ** 2 * math.cos(d + t + p) * (1.0 + root) ** 2
    )


def ka_rankine(phi: float, beta: float = 0.0) -> float:
    """Rankine 주동토압계수(경사 배면 β · 합력은 β 로 기욺)."""
    f, b = math.radians(phi), math.radians(beta)
    root = math.sqrt(max(math.cos(b) ** 2 - math.cos(f) ** 2, 0.0))
    return math.cos(b) * (math.cos(b) - root) / (math.cos(b) + root)


def bearing_factors(phi: float) -> tuple[float, float, float]:
    """Meyerhof 지지력 계수(Nc · Nq · Nγ) — φ 도."""
    f = math.radians(phi)
    if f < 1e-9:
        return 5.14, 1.0, 0.0
    nq = math.exp(math.pi * math.tan(f)) * math.tan(math.pi / 4 + f / 2) ** 2
    return (nq - 1.0) / math.tan(f), nq, (nq - 1.0) * math.tan(1.4 * f)


def _pressure(
    ka: float, H: float, gamma: float, gamma_sat: float, q: float, hw: float
) -> tuple[float, float]:
    """연직 높이 H 가상면 토압 합력 · 바닥에서 작용 높이 — σ' = 수위 위 γ · 아래 γ' 누적 + q."""
    z = np.linspace(0.0, H, _N + 1)  # 꼭대기에서 아래로
    dry = np.minimum(z, H - hw)
    wet = np.maximum(z - (H - hw), 0.0)
    p = ka * (gamma * dry + (gamma_sat - GAMMA_W) * wet + q)
    total = float(np.trapz(p, z))
    arm = float(np.trapz(p * (H - z), z)) / total if total > 0 else H / 3
    return total, arm


def check_wall(
    polygon: list,
    back: int,
    *,
    backfill: dict[str, float],
    foundation: dict[str, float],
    method: str = "coulomb",
    delta_ratio: float = 2.0 / 3.0,
    beta_deg: float = 0.0,
    surcharge_kpa: float = 10.0,
    water_ratio: float = 0.0,
    kh: float = 0.0,
    wall_gamma: float = WALL_GAMMA["default"],
    base_friction: float | None = None,
    base_adhesion: float = 0.0,
    embed_m: float = 0.0,
) -> dict[str, Any]:
    """벽 하나 외적 안정. back = 배면(흙) 쪽 x 방향(+1 · −1). 물성 = {gamma_t, gamma_sat, c, phi}.

    kh > 0 이면 지진 시(기준 `REQUIRED["seismic"]`) · 아니면 평상시. 돌려줌 = 하중 · 안전율 · 판정."""
    wall = Polygon(polygon)
    xs = np.array([p[0] for p in polygon], dtype=float)
    zs = np.array([p[1] for p in polygon], dtype=float)
    zb, zt = float(zs.min()), float(zs.max())
    H = zt - zb
    base = xs[zs <= zb + 1e-3]
    heel = float(base.max() if back > 0 else base.min())
    toe = float(base.min() if back > 0 else base.max())
    B = abs(heel - toe)
    phi, gamma = float(backfill["phi"]), float(backfill["gamma_t"])
    gamma_sat = float(backfill.get("gamma_sat") or gamma)
    hw = min(max(water_ratio, 0.0), 1.0) * H

    if method == "rankine":
        x_top, theta, delta = heel, 0.0, beta_deg
    else:
        top = xs[zs >= zt - 1e-3]
        x_top = float(top.max() if back > 0 else top.min())  # 배면 꼭대기 모서리
        theta = math.degrees(math.atan2(-(x_top - heel) * back, H))
        delta = delta_ratio * phi
    # 가상면과 벽 사이 흙 — 벽 쪽 반평면 ∩ 벽 상자 − 벽 · 가상면에 붙은 조각만(앞면 앞 흙 빼고)
    reach = 10.0 * (H + B + 1.0)
    side = [(heel, zb), (x_top, zt), (x_top - back * reach, zt), (heel - back * reach, zb)]
    plane = LineString([(heel, zb), (x_top, zt)])
    rest = Polygon(side).intersection(box(*wall.bounds)).difference(wall)
    parts = [g for g in getattr(rest, "geoms", [rest]) if isinstance(g, Polygon)]
    soil = unary_union([g for g in parts if g.area > 1e-9 and g.distance(plane) < 1e-6])
    soil_w = soil.area * gamma
    soil_x = soil.centroid.x if soil_w > 0 else heel
    soil_z = soil.centroid.y if soil_w > 0 else zb

    ka = (
        ka_rankine(phi, beta_deg)
        if method == "rankine"
        else ka_coulomb(phi, delta, theta, beta_deg)
    )
    q_eff = surcharge_kpa * (
        1.0
        if method == "rankine"
        else math.cos(math.radians(theta)) / math.cos(math.radians(theta - beta_deg))
    )
    Pa, ya = _pressure(ka, H, gamma, gamma_sat, q_eff, hw)
    loads = [("earth", Pa, ya)]
    kae = None
    if kh > 0:
        psi = math.degrees(math.atan(kh))
        kae = ka_coulomb(phi, delta, theta, beta_deg, psi)
        Pae, _ = _pressure(kae, H, gamma, gamma_sat, q_eff, hw)
        loads.append(("seismic_earth", max(Pae - Pa, 0.0), 0.6 * H))
    incline = math.radians(delta + theta)

    def plane_x(y: float) -> float:  # 가상면 위 높이 y 자리 x
        return heel + (x_top - heel) * (y / H if H > 0 else 0.0)

    def arm(x: float) -> float:  # 앞굽에서 배면 쪽 수평 팔 길이
        return (x - toe) * back

    W = wall.area * wall_gamma
    xw, zw = wall.centroid.x, wall.centroid.y
    Ph = sum(P for _, P, _ in loads) * math.cos(incline)
    Pv = sum(P for _, P, _ in loads) * math.sin(incline)
    m_res = W * arm(xw) + soil_w * arm(soil_x)
    m_over = 0.0
    for _, P, y in loads:
        m_res += P * math.sin(incline) * arm(plane_x(y))
        m_over += P * math.cos(incline) * y
    Pw = 0.5 * GAMMA_W * hw**2
    U = 0.5 * GAMMA_W * hw * B
    m_over += Pw * hw / 3.0 + U * (2.0 * B / 3.0)
    inertia = kh * (W + soil_w)
    m_over += kh * W * (zw - zb) + kh * soil_w * (soil_z - zb)
    H_total = Ph + Pw + inertia
    V = W + soil_w + Pv - U

    phi_f = float(foundation["phi"])
    mu = base_friction if base_friction is not None else math.tan(math.radians(2.0 / 3.0 * phi_f))
    sliding = (V * mu + base_adhesion * B) / H_total if H_total > 0 else math.inf
    overturning = m_res / m_over if m_over > 0 else math.inf
    xbar = (m_res - m_over) / V if V > 0 else -math.inf
    e = B / 2.0 - xbar
    if xbar <= 0:
        q_max, b_eff = math.inf, 0.0
    elif abs(e) <= B / 6.0:
        q_max, b_eff = V / B * (1.0 + 6.0 * abs(e) / B), B - 2.0 * abs(e)
    else:
        q_max, b_eff = 2.0 * V / (3.0 * min(xbar, B - xbar)), B - 2.0 * abs(e)
    nc, nq, ng = bearing_factors(phi_f)
    alpha = math.degrees(math.atan2(H_total, V)) if V > 0 else 90.0
    i_cq = (1.0 - alpha / 90.0) ** 2
    i_g = (1.0 - alpha / phi_f) ** 2 if phi_f > 0 and alpha < phi_f else 0.0
    gf = float(foundation["gamma_t"])
    q_ult = (
        float(foundation.get("c") or 0.0) * nc * i_cq
        + gf * max(embed_m, 0.0) * nq * i_cq
        + 0.5 * gf * max(b_eff, 0.0) * ng * i_g
    )
    bearing = q_ult / q_max if q_max > 0 and math.isfinite(q_max) else 0.0
    case = "seismic" if kh > 0 else "normal"
    fs = {"sliding": sliding, "overturning": overturning, "bearing": bearing}
    need = REQUIRED[case]
    return {
        "case": case,
        "method": method,
        "geometry": {"H": H, "B": B, "toe_x": toe, "heel_x": heel, "theta_deg": theta},
        "coefficients": {"ka": ka, "kae": kae, "delta_deg": delta, "beta_deg": beta_deg},
        "loads": {
            "W": W,
            "W_soil": soil_w,
            "Pa": Pa,
            "Pa_arm": ya,
            "Ph": Ph,
            "Pv": Pv,
            "Pw": Pw,
            "U": U,
            "inertia": inertia,
            "V": V,
            "H": H_total,
            "M_resist": m_res,
            "M_overturn": m_over,
        },
        "base": {"e": e, "q_max": q_max, "q_ult": q_ult, "B_eff": b_eff, "mu": mu, "alpha": alpha},
        "fs": fs,
        "required": dict(need),
        "ok": {k: bool(fs[k] >= need[k]) for k in fs},
    }