Aislo/E03_StabSection/E03_StabSection_Engine_Support.py

"""E03 보강재(네일 · 앵커) — 활동면을 지나는 보강력을 절편 바닥 점하중으로(Slide2 Support 와 같은 꼴).

입력 꼴(E03 보강 검토 화면과 한 벌) supports = [{
  "type": "nail" | "anchor", "x", "z" (머리 · 모형 좌표 m), "angle_deg" (수평 아래로 +), "length_m",
  "spacing_m" (종방향 간격 · 단위 폭 힘 = 힘 / 간격), "tensile_kN" (강재 인장 한도),
  "bond_kPa" (그라우트 · 지반 부착 · 네일 = 전 길이 · 앵커 = 끝 정착장만), "bond_length_m" (앵커 정착장),
  "hole_dia_mm", "prestress_kN" (앵커만 · 0 이면 수동), "plate_kN"? (머리 판 · 없으면 tensile_kN),
  "side"? ("left" | "right" · 고르는 것은 `Engine_Station`)}]
힘 크기 = 지점마다 Slide2 Force Diagram 의 최소 — 인장(tensile) · 인발(pullout = 부착 × 활동면 밖 부착 길이) ·
  벗겨짐(stripping = 판 + 부착 × 활동면 안 부착 길이) · 부착 = π·D·bond_kPa (kN/m) · 자유장은 부착 없음.
  인발 안전율은 입력 부착값에 넣어 둠(Slide2 와 같음 — 엔진은 나누지 않음).
방향 = 보강재와 나란히(Slide2 Parallel to Reinforcement) · 비탈 속으로(활동 방향 반대 = −x) 아래로 angle.
능동 · 수동 = Slide2 Active/Passive Force Application:
  능동(프리스트레스 앵커) = 힘 그대로 — 활동력을 줄임 F = R / (D − T) · 힘 = min(프리스트레스, 한도)
  수동(네일 · 프리스트레스 0 앵커) = 힘 / F — 저항을 늘림 F = (R + T) / D
머리가 활동 덩이 밖이면(활동면이 보강재를 안 지나거나 덩이 밖에서 드나듦) 0 = not_crossing.
절편 칸(꼴 (S, n) · 활동 방향 +x 틀): sup_va · sup_ha = 능동 연직(아래 +) · 수평(+x) · sup_vp · sup_hp = 수동(× 1/F 전)
  · sup_xva · sup_zha · sup_xvp · sup_zhp = 모멘트용 x·연직 · z·수평. 보강재 칸(꼴 (S, K)) sup_T · sup_mode · sup_x · sup_z.
"""

import numpy as np

MODES = ("tensile", "pullout", "stripping", "not_crossing", "prestress")
_NOT, _PRE = 3, 4
_FIELDS = ("va", "ha", "vp", "hp", "xva", "zha", "xvp", "zhp")


def prepare(raw: list[dict] | None, direction: int) -> dict | None:
    """입력 꼴 → 배열(활동 방향 +x 틀 · x 에 direction 곱함) · 없으면 None."""
    if not raw:
        return None
    s = 1.0 if direction >= 0 else -1.0
    cols = {k: [] for k in ("hx", "hz", "dx", "dz", "L", "free", "q", "tens", "plate", "sp", "pre")}
    active = []
    for r in raw:
        kind = r.get("type", "nail")
        if kind not in ("nail", "anchor"):
            raise ValueError(f"보강재 종류 모름: {kind}")
        L, sp = float(r["length_m"]), float(r["spacing_m"])
        tens, bond = float(r["tensile_kN"]), float(r["bond_kPa"])
        dia = float(r["hole_dia_mm"]) / 1000.0
        lb = float(r.get("bond_length_m") or 0.0) if kind == "anchor" else L
        pre = float(r.get("prestress_kN") or 0.0) if kind == "anchor" else 0.0
        plate = float(r["plate_kN"]) if r.get("plate_kN") is not None else tens
        if L <= 0 or sp <= 0 or dia <= 0 or tens < 0 or bond < 0 or plate < 0 or pre < 0:
            raise ValueError(f"보강재 값 범위 밖: {r}")
        if not 0 < lb <= L:
            raise ValueError(f"앵커 정착장은 0 초과 · 길이 이하: {r}")
        a = np.radians(float(r["angle_deg"]))
        for k, v in (
            ("hx", s * float(r["x"])),
            ("hz", float(r["z"])),
            ("dx", -np.cos(a)),
            ("dz", -np.sin(a)),
            ("L", L),
            ("free", L - lb),
            ("q", np.pi * dia * bond),
            ("tens", tens),
            ("plate", plate),
            ("sp", sp),
            ("pre", pre),
        ):
            cols[k].append(v)
        active.append(pre > 0)
    out = {k: np.array(v, float) for k, v in cols.items()}
    out["active"] = np.array(active, bool)
    return out


def tips(sup: dict) -> np.ndarray:
    """보강재 끝 x (활동 방향 +x 틀) — 탐색 한계를 끝 뒤까지 넓힘."""
    return sup["hx"] + sup["L"] * sup["dx"]


def apply(sl: dict, sup: dict) -> None:
    """절편 묶음(sl)에 보강력 칸을 덧붙임 — 면마다 보강재마다 첫 교점 · 힘 · 파괴 꼴."""
    xe, ze = sl["xe"], sl["ze"]
    S, n = sl["b"].shape
    hx, hz, dx, dz = (sup[k][None, :, None] for k in ("hx", "hz", "dx", "dz"))
    ax, az = xe[:, None, :-1], ze[:, None, :-1]
    ex, ez = np.diff(xe, axis=1)[:, None, :], np.diff(ze, axis=1)[:, None, :]
    den = dx * ez - dz * ex
    safe = np.where(np.abs(den) < 1e-12, np.nan, den)
    wx, wz = ax - hx, az - hz
    t = (wx * ez - wz * ex) / safe  # 보강재 위 머리에서 거리
    u = (wx * dz - wz * dx) / safe  # 활동면 마디 위 비율
    hit = (u >= 0) & (u < 1) & (t > 0) & (t <= sup["L"][None, :, None])
    tt = np.where(hit, t, np.inf)
    j = np.argmin(tt, axis=-1)  # (S, K) 첫 교점 절편
    tc = np.take_along_axis(tt, j[..., None], -1)[..., 0]
    # 머리가 활동 덩이 안(진입 ~ 진출 사이 · 활동면 위)일 때만 — 첫 교점 = 덩이에서 나가는 점
    k_head = np.clip((xe[:, None, :] <= sup["hx"][None, :, None]).sum(-1) - 1, 0, n - 1)
    x0 = np.take_along_axis(xe, k_head, 1)
    x1 = np.take_along_axis(xe, k_head + 1, 1)
    z0 = np.take_along_axis(ze, k_head, 1)
    z1 = np.take_along_axis(ze, k_head + 1, 1)
    zslip = z0 + (sup["hx"][None, :] - x0) * (z1 - z0) / np.maximum(x1 - x0, 1e-12)
    inside = (sup["hx"] >= xe[:, :1]) & (sup["hx"] <= xe[:, -1:]) & (sup["hz"] > zslip)
    cross = inside & np.isfinite(tc) & sl["valid"][:, None]
    tc = np.where(cross, tc, 0.0)
    lb = sup["L"] - sup["free"]
    bond_in = np.clip(tc - sup["free"], 0.0, lb)
    caps = np.stack(
        [
            np.broadcast_to(sup["tens"], tc.shape),
            sup["q"] * (lb - bond_in),
            sup["plate"] + sup["q"] * bond_in,
        ],
        axis=-1,
    )
    mode = np.argmin(caps, axis=-1)
    T = np.min(caps, axis=-1)
    pre_rules = sup["active"] & (sup["pre"] < T)
    T = np.where(pre_rules, sup["pre"], T)
    mode = np.where(pre_rules, _PRE, mode)
    T = np.where(cross, T, 0.0) / sup["sp"]
    mode = np.where(cross, mode, _NOT)
    xs = sup["hx"] + tc * sup["dx"]
    zs = sup["hz"] + tc * sup["dz"]
    acc = {f: np.zeros((S, n)) for f in _FIELDS}
    rows = np.arange(S)
    for k in range(len(sup["L"])):
        v, h = -T[:, k] * sup["dz"][k], T[:, k] * sup["dx"][k]  # 아래 + · +x(비탈 속 = −)
        tag = "a" if sup["active"][k] else "p"
        for f, val in ((f"v{tag}", v), (f"h{tag}", h), (f"xv{tag}", xs[:, k] * v)):
            np.add.at(acc[f], (rows, j[:, k]), val)
        np.add.at(acc[f"zh{tag}"], (rows, j[:, k]), zs[:, k] * h)
    for f in _FIELDS:
        sl[f"sup_{f}"] = acc[f]
    sl["sup_T"], sl["sup_mode"], sl["sup_x"], sl["sup_z"] = T, mode, xs, zs


def forces(sl: dict, k: int, F: float, active: np.ndarray, direction: int) -> list[dict]:
    """면 k 의 보강재별 결과 칸(support_forces) — 한도(단위 폭) · 실제 넣은 힘(수동 = / F) · 파괴 꼴 · 교점."""
    out = []
    for i, (T, m) in enumerate(zip(sl["sup_T"][k], sl["sup_mode"][k])):
        crossing = int(m) != _NOT
        out.append(
            {
                "index": i,
                "mobilized_kN_per_m": float(T),
                "applied_kN_per_m": float(T if active[i] else T / F),
                "mode": MODES[int(m)],
                "force_application": "active" if active[i] else "passive",
                "point": (
                    [float(direction * sl["sup_x"][k, i]), float(sl["sup_z"][k, i])]
                    if crossing
                    else None
                ),
            }
        )
    return out