Aislo/E03_StabSection/E03_StabSection_Engine_Model.py

"""E03 2D 지층 모델 빌더 — 측점 횡단(지반선 · 설계선 · 소단 · 측구 · 암 경계 · 성토체 · 벽) → 지층 폴리곤.

`build_section_model` = storage 를 모르는 순수 함수(좌표 · 값만) · `section_model_from_b06` = B06 횡단 자료 꼴
(분석 345-D ②)을 읽는 얇은 어댑터. 좌표 x = B06 offset(좌 +) · z = 표고.
설계선(design_line)이 소단 · 측구 · 노면을 이미 담음 — 여기서 다시 그리지 않음. 설계선 밖은 원지반.
암 경계 = 원지반을 연직으로 |offset| 내린 선(B06 면적 가름과 같은 정의 · `Engine_Design.py` 363-377).
"""

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

from E03_StabSection.E03_StabSection_Types import (
    DEFAULT_MATERIALS,
    DEFAULT_ROAD_SURCHARGE_KPA,
)
from E03_StabSection.E03_StabSection_Types import (
    groundwater as norm_groundwater,
)
from E03_StabSection.E03_StabSection_Types import (
    surcharges as norm_surcharges,
)

#: B06 지반 라벨 → 암 경계 아래 지층(요구서 345 6장 · 분석 D 5-2). 경계 위는 풍화토.
ROCK_LAYER_OF = {"ripping_rock": "weathered_rock", "blasting_rock": "soft_rock"}
SOIL_ID = "weathered_soil"
FILL_ID = "fill"
#: 벽(C군) — 활동면이 못 지나게 강한 물체로(콘크리트 · 석축 대표값). 물성 표에 없으면 이 값을 덧붙임.
WALL_MATERIAL = {
    "id": "wall",
    "name": "벽체",
    "gamma_t": 23.0,
    "gamma_sat": 23.0,
    "c": 1000.0,
    "phi": 45.0,
}
DEFAULT_HALF_WIDTH_M = 40.0


def _polys(geom) -> list[Polygon]:
    geom = make_valid(geom)
    if isinstance(geom, Polygon):
        parts = [geom]
    elif isinstance(geom, MultiPolygon):
        parts = list(geom.geoms)
    else:
        parts = [g for g in getattr(geom, "geoms", []) if isinstance(g, Polygon)]
    return [p for p in parts if p.area > 1e-6]


def _region(material_id: str, poly: Polygon) -> dict:
    out = {"material_id": material_id, "polygon": [list(c) for c in poly.exterior.coords]}
    if poly.interiors:
        out["holes"] = [[list(c) for c in ring.coords] for ring in poly.interiors]
    return out


def final_surface(ground: list, design_line: list) -> list[list[float]]:
    """공사 뒤 지표 = 설계선 범위는 설계선 · 밖은 원지반."""
    line = sorted((float(x), float(z)) for x, z in design_line)
    if not line:
        return [list(p) for p in ground]
    x0, x1 = line[0][0], line[-1][0]
    left = [p for p in ground if p[0] < x0 - 1e-9]
    right = [p for p in ground if p[0] > x1 + 1e-9]
    return [list(map(float, p)) for p in left + line + right]


def build_section_model(
    ground: list,
    design_line: list,
    *,
    layers: list,
    materials: list[dict] | None = None,
    fill_material: str = FILL_ID,
    walls: list[dict] | None = None,
    groundwater: dict | None = None,
    seismic_kh: float = 0.0,
    surcharge: list[dict] | None = None,
    bottom_z: float | None = None,
    area_range: tuple[float, float] | None = None,
) -> dict:
    """원지반 · 설계선 · 지층(위부터 [(material_id, 원지반 아래 바닥 깊이 m | None=모형 바닥)]) → 해석 모형.

    돌려주는 값 = LEM 모형(materials · regions · surface · groundwater · seismic_kh · surcharge)
    + `areas`(area_range 안 절토 · 성토 · 절토 토사 · 절토 암 ㎡ — B06 면적 대조용).
    """
    ground = sorted((float(x), float(z)) for x, z in ground)
    surface = final_surface(ground, design_line)
    xs = [p[0] for p in ground]
    zs = [p[1] for p in ground] + [p[1] for p in surface]
    relief = max(zs) - min(zs)
    zb = bottom_z if bottom_z is not None else min(zs) - max(5.0, 0.5 * relief)
    x0, x1 = xs[0], xs[-1]
    bottom = [(x1, zb), (x0, zb)]
    below_final = make_valid(Polygon(surface + bottom))
    ground_body = make_valid(Polygon(ground + bottom))
    natural = below_final.intersection(ground_body)
    fill = below_final.difference(ground_body)
    cut = ground_body.difference(below_final)

    def band(d_top: float, d_bot: float | None) -> Polygon:
        top = [(x, z - d_top) for x, z in ground]
        if d_bot is None:
            low = [(x1, zb - 1.0), (x0, zb - 1.0)]
        else:
            low = [(x, z - d_bot) for x, z in reversed(ground)]
        return make_valid(Polygon(top + low))

    wall_geom = unary_union([Polygon(w["polygon"]) for w in walls or []]) if walls else None
    regions: list[dict] = []
    bands: list = []
    depth = 0.0
    for material_id, d_bot in layers:
        b = band(depth, None if d_bot is None else float(d_bot))
        bands.append(b)
        part = natural.intersection(b)
        if wall_geom is not None:
            part = part.difference(wall_geom)
        regions += [_region(material_id, p) for p in _polys(part)]
        if d_bot is None:
            break
        depth = float(d_bot)
    fill_part = fill.difference(wall_geom) if wall_geom is not None else fill
    regions += [_region(fill_material, p) for p in _polys(fill_part)]
    mats = [dict(m) for m in (materials or DEFAULT_MATERIALS)]
    if wall_geom is not None:
        for w in walls:
            regions += [
                _region(w.get("material_id", "wall"), p) for p in _polys(Polygon(w["polygon"]))
            ]
        if not any(m["id"] == "wall" for m in mats):
            mats.append(dict(WALL_MATERIAL))
    a0, a1 = area_range or (x0, x1)
    clip = box(a0, zb - 1.0, a1, max(zs) + 100.0)
    cut_in = cut.intersection(clip)
    cut_soil = cut_in.intersection(bands[0]).area if len(bands) > 1 else cut_in.area
    return {
        "materials": mats,
        "regions": regions,
        "surface": surface,
        "groundwater": norm_groundwater(groundwater),
        "seismic_kh": float(seismic_kh or 0.0),
        "surcharge": norm_surcharges(surcharge),
        "bottom_z": zb,
        "areas": {
            "cut_m2": cut_in.area,
            "fill_m2": fill.intersection(clip).area,
            "cut_soil_m2": cut_soil,
            "cut_rock_m2": cut_in.area - cut_soil,
        },
    }


# ---- B06 어댑터 -----------------------------------------------------------------
def ground_from_cross(cross: dict) -> list[tuple[float, float]]:
    return sorted(
        (float(s["offset_m"]), float(s["elevation_m"]))
        for s in cross.get("samples") or []
        if s.get("valid") is not False and s.get("elevation_m") is not None
    )


def widen_ground(cross: dict, sampler, half_width: float = DEFAULT_HALF_WIDTH_M, step: float = 0.5):
    """해석 폭 — 횡단 표본 바깥을 지형에서 다시 떠 ±half_width 까지 넓힘(지형 밖이 나오면 거기서 멈춤)."""
    ground = ground_from_cross(cross)
    frame = cross["frame"]
    origin = np.array([float(frame["origin"]["x"]), float(frame["origin"]["y"])])
    axis = np.array([float(v) for v in frame["left_xy"]])
    added = []
    for sign, edge in ((-1.0, ground[0][0]), (1.0, ground[-1][0])):
        n = int(max(0.0, half_width - abs(edge)) / step)
        if n <= 0:
            continue
        offs = edge + sign * step * np.arange(1, n + 1)
        z, valid = sampler.sample_xy(origin[None, :] + axis[None, :] * offs[:, None])
        stop = int(np.argmin(valid)) if not valid.all() else len(valid)
        added += [(float(o), float(e)) for o, e in zip(offs[:stop], z[:stop])]
    return sorted(ground + added)


def side_cases(design: dict, ground: list, line: list) -> list[dict]:
    """좌 · 우 비탈마다 역할(절토 · 성토) · 활동 방향 · 탐색 한계. 좌 = +offset.

    탐색 한계 = 노면 반대쪽 끝 ~ 비탈 끝(설계선이 원지반과 만나는 점) 바깥 max(5 m, 2H) — 설계 비탈을
    지나는 면만 봄(바깥 자연 비탈 전체 붕괴는 E 광역 · 2차 몫). H = 노면 끝 ~ 비탈 끝 높이차.
    """
    mode = design.get("section_mode") or "right_cut"
    roles = {
        "left_cut": ("cut", "fill"),
        "right_cut": ("fill", "cut"),
        "both_cut": ("cut", "cut"),
        "both_fill": ("fill", "fill"),
    }[mode]
    edges = design.get("road_edges") or {}
    left_edge = float((edges.get("left") or {}).get("offset_m", 2.0))
    right_edge = float((edges.get("right") or {}).get("offset_m", -2.0))
    gx = np.array([p[0] for p in ground])
    gz = np.array([p[1] for p in ground])
    lx = np.array([p[0] for p in line])
    lz = np.array([p[1] for p in line])
    off = np.abs(lz - np.interp(lx, gx, gz)) > 0.005
    out = []
    for side, role in zip(("left", "right"), roles):
        inward = role == "cut"  # 절토 = 산 쪽에서 노면으로 · 성토 = 노면에서 바깥으로
        sign = 1.0 if side == "left" else -1.0
        edge = left_edge if side == "left" else right_edge
        # 비탈 끝 = 원지반에서 떨어진 가장 바깥 꼭짓점의 바로 바깥 꼭짓점(다시 원지반에 붙는 점)
        idx = [k for k, (x, o) in enumerate(zip(lx, off)) if o and sign * x >= sign * edge - 1e-9]
        if idx:
            k = max(idx) if sign > 0 else min(idx)
            k = min(max(k + int(sign), 0), len(lx) - 1)
            day = float(lx[k])
        else:
            day = edge
        H = abs(float(np.interp(day, lx, lz)) - float(np.interp(edge, lx, lz)))
        reach = day + sign * max(5.0, 2.0 * H)
        reach = min(max(reach, gx[0]), gx[-1])
        inner = right_edge if side == "left" else left_edge
        limits = tuple(sorted((inner, reach)))
        direction = int(-sign) if inward else int(sign)
        out.append(
            {"side": side, "role": role, "direction": direction, "limits": limits, "height_m": H}
        )
    return out


def section_model_from_b06(
    cross: dict,
    design: dict,
    *,
    materials: list[dict] | None = None,
    layers: list | None = None,
    ground: list | None = None,
    wet: bool = False,
    seismic_kh: float = 0.0,
    road_surcharge_kpa: float = DEFAULT_ROAD_SURCHARGE_KPA,
    walls: list[dict] | None = None,
) -> dict:
    """B06 횡단 파일(samples · frame) + 설계(design) → 해석 모형 + 비탈 목록(`sides`).

    ground = 넓혀 다시 뜬 지반선(`widen_ground`) · 없으면 횡단 표본 그대로.
    layers 기본 = 토사 측점 풍화토 한 층 · 암 측점 풍화토(암 경계까지) + 풍화암 / 연암.
    wet = 우기 — 지하수위 = 공사 뒤 지표(KDS 11 70 05 표 4.3-3 · 토사 · 풍화암 지표 포화).
    노면 하중 = 노면 끝 사이 등분포. 벽 도형(walls)은 B06 화면 산식(TS)이 그리므로 부르는 쪽이 넘김.
    """
    ground = ground or ground_from_cross(cross)
    line = [
        (float(p["offset_m"]), float(p["elevation_m"])) for p in design.get("design_line") or []
    ]
    if layers is None:
        rock = ROCK_LAYER_OF.get(design.get("ground_type"))
        offset = design.get("rock_boundary_offset_m")
        if rock and offset is not None:
            layers = [(SOIL_ID, abs(float(offset))), (rock, None)]
        elif rock:
            layers = [(rock, None)]
        else:
            layers = [(SOIL_ID, None)]
    edges = design.get("road_edges") or {}
    loads = []
    if road_surcharge_kpa and edges.get("left") and edges.get("right"):
        loads = [
            {
                "x0": edges["right"]["offset_m"],
                "x1": edges["left"]["offset_m"],
                "q": road_surcharge_kpa,
            }
        ]
    surface = final_surface(ground, line)
    gw = {"mode": "line", "line": surface} if wet else {"mode": "none"}
    model = build_section_model(
        ground,
        line,
        layers=layers,
        materials=materials,
        walls=walls,
        groundwater=gw,
        seismic_kh=seismic_kh,
        surcharge=loads,
        area_range=(line[0][0], line[-1][0]) if line else None,
    )
    model["sides"] = side_cases(design, ground, line or ground)
    model["station"] = {"chainage_m": cross.get("chainage_m"), "label": cross.get("label")}
    return model


def polygon_closed(region: dict) -> bool:
    """폴리곤 닫힘 · 자기 교차 없음 · 넓이 > 0 (시험용)."""
    ring = region["polygon"]
    poly = Polygon(ring, region.get("holes"))
    return bool(poly.is_valid and poly.area > 0 and LineString(ring).is_closed)