Aislo/resources/tester/test_532_2_stab_section_model.py

# -*- coding: utf-8 -*-
"""532-2 E03 2D 지층 모델 빌더 — 폴리곤 닫힘 · 지층 · 성토체 · 절성토 면적 = B06 · 비탈 목록 · 측점 해석 꼴.

- 순수 함수(`build_section_model`)는 만든 단면으로 ·
  B06 대조는 검증_sub4_E 복사본(읽기만 · 없으면 건너뜀) ·
  E03 분석 길 그대로(관 · 벽 측점 = `walls_by_station` 그려지는 설계선 · B06 면적도 그 선).

    ./venv/Scripts/python.exe -m pytest resources/tester/test_532_2_stab_section_model.py -q
"""

from __future__ import annotations

import sys
from pathlib import Path

import pytest
from shapely.geometry import Polygon

ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))

from common_util.common_util_station_key import chainage_key  # noqa: E402
from E01_StabRegister import E01_StabRegister_Source as source  # noqa: E402
from E03_StabSection.E03_StabSection_Engine_Model import (  # noqa: E402
    build_section_model,
    polygon_closed,
    section_model_from_b06,
    side_cases,
)
from E03_StabSection.E03_StabSection_Engine_Station import analyze_station  # noqa: E402
from E03_StabSection.E03_StabSection_Walls import (  # noqa: E402
    design_with_walls,
    walls_by_station,
    walls_of,
)

SUB4_E = ROOT / "storage/1/3/1fe7901c-d26b-443d-89a2-28172ec438bc"

# 오른쪽(−) 산 · 왼쪽(+) 골 · 편절편성: 노면 −2~2 · 우 절토 1:0.5 · 좌 성토 1:1.5
GROUND = [(x / 2, 10.0 - 0.6 * (x / 2)) for x in range(-40, 41)]  # 10 − 0.6x (좌로 내려감)
DESIGN = [
    (-20.0, 22.0),
    (-6.0, 13.6),
    (-2.6, 10.0),
    (-2.0, 10.0),
    (2.0, 10.0),
    (14.0, 2.0),
    (20.0, -2.0),
]


def _trapezoid(ground, line, sign):
    """설계선 꼭짓점 · 지반 표본 합집합 위 사다리꼴(B06 `_trapezoid_areas` 와 같은 정의) — 절토(+) / 성토(−)."""
    import numpy as np

    gx, gz = zip(*ground)
    lx, lz = zip(*line)
    xs = sorted({*gx, *lx})
    xs = [x for x in xs if lx[0] <= x <= lx[-1]]
    d = np.interp(xs, gx, gz) - np.interp(xs, lx, lz)
    area = 0.0
    for x0, x1, d0, d1 in zip(xs, xs[1:], d, d[1:]):
        if d0 * d1 < 0:
            xm = x0 + (x1 - x0) * d0 / (d0 - d1)
            area += 0.5 * (xm - x0) * max(sign * d0, 0) + 0.5 * (x1 - xm) * max(sign * d1, 0)
        else:
            area += 0.5 * (x1 - x0) * (max(sign * d0, 0) + max(sign * d1, 0))
    return area


def test_build_model_polygons_layers_fill_and_areas():
    m = build_section_model(
        GROUND,
        DESIGN,
        layers=[("weathered_soil", 0.5), ("weathered_rock", None)],
        area_range=(-20.0, 20.0),
    )
    assert m["regions"] and all(polygon_closed(r) for r in m["regions"])
    ids = {r["material_id"] for r in m["regions"]}
    assert ids == {"weathered_soil", "weathered_rock", "fill"}
    # 지층 사이 겹침 없음 · 합 = 공사 뒤 지표 아래 전체
    polys = [Polygon(r["polygon"], r.get("holes")) for r in m["regions"]]
    total = sum(p.area for p in polys)
    union = polys[0]
    for p in polys[1:]:
        union = union.union(p)
    assert abs(union.area - total) < 1e-6
    a = m["areas"]
    assert abs(a["cut_m2"] - _trapezoid(GROUND, DESIGN, 1)) < 1e-6
    assert abs(a["fill_m2"] - _trapezoid(GROUND, DESIGN, -1)) < 1e-6
    assert (
        0 < a["cut_soil_m2"] < a["cut_m2"]
        and abs(a["cut_soil_m2"] + a["cut_rock_m2"] - a["cut_m2"]) < 1e-9
    )


def test_side_cases_roles_directions_limits():
    design = {
        "section_mode": "right_cut",
        "road_edges": {"left": {"offset_m": 2.0}, "right": {"offset_m": -2.0}},
    }
    left, right = side_cases(design, GROUND, DESIGN)
    assert (left["role"], left["direction"]) == ("fill", 1)  # 좌 성토 → 바깥(+x)으로
    assert (right["role"], right["direction"]) == ("cut", 1)  # 우 절토 → 노면(+x)으로
    assert left["limits"][0] == -2.0 and right["limits"][1] == 2.0
    assert left["height_m"] > 0 and right["height_m"] > 0


def test_analyze_station_shape():
    cross = {
        "chainage_m": 0.0,
        "label": "STA.0+000",
        "samples": [{"offset_m": x, "elevation_m": z} for x, z in GROUND],
    }
    design = {
        "section_mode": "right_cut",
        "ground_type": "ripping_rock",
        "rock_boundary_offset_m": -0.5,
        "design_line": [{"offset_m": x, "elevation_m": z} for x, z in DESIGN],
        "road_edges": {
            "left": {"offset_m": 2.0, "elevation_m": 10.0},
            "right": {"offset_m": -2.0, "elevation_m": 10.0},
        },
    }
    res = analyze_station(
        cross, design, conditions=("dry",), methods=("bishop",), noncircular=False
    )
    assert len(res["cases"]) == 2
    for case in res["cases"]:
        F = case["min_F"]["bishop"]
        assert F is not None and 0.3 < F < 20
        best = case["result"]["methods"]["bishop"]["circular"]
        assert best["slices"] and best["center"] and best["radius"] > 0
        lo, hi = case["limits"]
        assert (
            lo - 1e-6 <= min(best["entry"][0], best["exit"][0])
            and max(best["entry"][0], best["exit"][0]) <= hi + 1e-6
        )
    assert res["cases"][0]["result"]["contour"]["F"]


@pytest.mark.skipif(not source.source_dir(SUB4_E).is_dir(), reason="검증_sub4_E 복사본 없음")
def test_b06_areas_match_all_stations():
    """검증_sub4_E 전 측점 — 모델 절토 · 성토 · 토사 · 암 면적 = B06 저장값(반올림 4 자리 안).

    관 · 벽 측점은 E03 분석 길 그대로(`design_with_walls` · `walls_of`) — B06 [저장] 면적이
    구조물 트림 뒤 「그려지는 설계선」 기준이라 원 설계선으로 재면 어긋남."""
    root = source.source_dir(SUB4_E)
    found = walls_by_station(root)
    stations = source.stations(root)
    assert stations and found  # 관 측점이 있어야 트림 길도 잼
    for station in stations:
        cross = source.cross_section(root, station["chainage_m"])
        design = cross["design"]
        name = cross.get("label")
        entry = found.get(chainage_key(station["chainage_m"]))
        m = section_model_from_b06(cross, design_with_walls(design, entry), walls=walls_of(entry))
        assert all(polygon_closed(r) for r in m["regions"]), name
        a = m["areas"]
        assert abs(a["cut_m2"] - design["cut_area_m2"]) < 2e-3, name
        assert abs(a["fill_m2"] - design["fill_area_m2"]) < 2e-3, name
        assert abs(a["cut_soil_m2"] - design["cut_soil_area_m2"]) < 2e-3, name
        assert abs(a["cut_rock_m2"] - design["cut_rock_area_m2"]) < 2e-3, name


if __name__ == "__main__":
    sys.exit(pytest.main([__file__, "-q"]))