Aislo/resources/tester/test_532_2_stab_mp_outputs.py

# -*- coding: utf-8 -*-
"""532-2 2차 E03 2D — Morgenstern-Price(반정현) 대조 · 결과 자료(자유물체도 · F–λ 곡선 · Safety Map).

기대값 = Rocscience Slide2 Slope Stability Verification Manual(지식DB 참고문헌) 원문 GLE · M-P 값.
- #1 ACADS 1a GLE 0.986 · #3 ACADS 1c GLE 1.374 · #8 ACADS 3b 지정 비원호 GLE 1.262 ·
  #9 ACADS 4 블록 + 최적화 GLE 0.683 · #21 F&K 지정 원 M-P 건기 2.075 · r_u 1.760(Slide 열).
- 뺀 것: #7 · #15 — 원문 GLE(1.275 · 0.437)가 같은 표 Spencer(1.246 · 0.409)보다 높음(Slide 탐색이 GLE
  임계면을 놓친 값 · 같은 문제 #8 지정면 GLE 는 Spencer 보다 낮음) → 대조 기준으로 못 씀.

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

from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pytest

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

from E03_StabSection import E03_StabSection_Engine_GLE as gle  # noqa: E402
from E03_StabSection import E03_StabSection_Engine_LEM as lem  # noqa: E402
from E03_StabSection.E03_StabSection_Engine_Search import (  # noqa: E402
    analyze,
    evaluate_specified,
)
from E03_StabSection.E03_StabSection_Engine_Surface import (  # noqa: E402
    poly_slices,
    three_point_center,
)
from test_532_2_stab_lem_benchmarks import (  # noqa: E402
    ACADS1C_MATS,
    ACADS1C_REGIONS,
    ACADS_EXT,
    FK_EXT,
    _acads3_model,
    _acads4_model,
    circ,
    mats,
    near,
    nonc,
)

MP = "morgenstern_price"
SEAM_3B = [[41.85, 27.75], [44.0, 26.5], [63.5, 27.0], [73.31, 40.0]]


def test_mp_constant_function_equals_spencer():
    # f(x) ≡ 1 이면 M-P 식 = Spencer 식 — 같은 F · λ
    prep = lem.Prepared(_acads3_model(), -1)
    pts = sorted((-x, z) for x, z in SEAM_3B)
    vx, vz = np.array([[p[0] for p in pts]]), np.array([[p[1] for p in pts]])
    sl = poly_slices(prep, vx, vz, 50)
    center = three_point_center(vx, vz)
    sp = lem.spencer(sl, center, np.array([1.3]))
    mp = gle.morgenstern_price(sl, center, "constant", F0=np.array([1.3]))
    assert abs(mp["F"][0] - sp["F"][0]) < 1e-6
    assert abs(mp["lambda"][0] - sp["lambda"][0]) < 1e-5


def test_mp_specified_surfaces_slide2():
    # #8 지정 비원호 GLE 1.262 · #21 지정 원 M-P 건기 2.075 · r_u 0.25 1.760
    r = evaluate_specified(_acads3_model(), polyline=SEAM_3B, direction=-1, methods=(MP,))
    near(r[MP]["F"], 1.262, 0.01)
    m = {"materials": mats(("s", 120, 600, 20)), "regions": [{"material_id": "s", "polygon": FK_EXT}]}
    near(evaluate_specified(m, circle=[120, 90, 80], direction=1, methods=(MP,))[MP]["F"], 2.075, 0.01)
    m["groundwater"] = {"mode": "ru", "ru": 0.25}
    near(evaluate_specified(m, circle=[120, 90, 80], direction=1, methods=(MP,))[MP]["F"], 1.760, 0.01)


def test_mp_search_slide2():
    # #1 원호 GLE 0.986 · #3 원호 GLE 1.374
    m = {"materials": mats(("s", 20, 3, 19.6)), "regions": [{"material_id": "s", "polygon": ACADS_EXT}]}
    near(circ(analyze(m, direction=-1, noncircular=False), MP), 0.986, 0.01)
    r = analyze({"materials": ACADS1C_MATS, "regions": ACADS1C_REGIONS}, direction=-1, noncircular=False)
    near(circ(r, MP), 1.374, 0.01)


def test_mp_acads4_noncircular():
    # #9 블록 + 최적화 GLE 0.683 (약층 좌표 = 그림 판독 · 비원호 다출발이 더 낮은 면을 찾음 → ±5 %)
    r = analyze(_acads4_model(), direction=-1, methods=("janbu", "spencer", MP))
    near(nonc(r, MP), 0.683, 0.05)


def _closure(rows, lam, f_last=1.0):
    weight = sum(r["weight"] for r in rows)
    for a, b in zip(rows, rows[1:]):
        assert abs(a["E_right"] - b["E_left"]) < 1e-9 and abs(a["X_right"] - b["X_left"]) < 1e-9
    assert abs(rows[-1]["E_right"]) < 1e-3 * weight  # 끝 절편에서 힘 평형이 닫힘
    assert abs(rows[0]["X_left"]) < 1e-12
    for r in rows:
        # 연직 평형 — N·cosα + S·sinα = W + (X_right − X_left) (지진 · 균열 없음)
        a = np.radians(r["alpha_deg"])
        lhs = r["normal_total"] * np.cos(a) + r["shear_mobilized"] * np.sin(a)
        assert abs(lhs - (r["weight"] + r["X_right"] - r["X_left"])) < 1e-6 * weight


def test_free_body_and_lambda_curve_outputs():
    m = {"materials": mats(("s", 20, 3, 19.6)), "regions": [{"material_id": "s", "polygon": ACADS_EXT}]}
    r = evaluate_specified(m, circle=[29.104, 55.527, 30.527], direction=-1, methods=("spencer", MP))
    for key in ("spencer", MP):
        out = r[key]
        rows = out["slices"]
        assert {"normal_total", "pore_force", "shear_strength", "E_right", "X_right", "x_mid"} <= rows[0].keys()
        _closure(rows, out["lambda"])
        curve = out["lambda_curve"]
        lam, fm, ff = (np.array(curve[k], float) for k in ("lambda", "Fm", "Ff"))
        assert len(lam) == gle.CURVE_POINTS and np.isfinite(fm).sum() > 10
        # 곡선 교점 = 해 (F · λ)
        d = fm - ff
        k = np.where(np.isfinite(d[:-1]) & np.isfinite(d[1:]) & (np.sign(d[:-1]) != np.sign(d[1:])))[0]
        assert len(k)
        j = k[np.argmin(np.abs(lam[k] - out["lambda"]))]
        t = d[j] / (d[j] - d[j + 1])
        assert abs(lam[j] + t * (lam[j + 1] - lam[j]) - out["lambda"]) < 0.02
        assert abs(fm[j] + t * (fm[j + 1] - fm[j]) - out["F"]) < 0.01 * out["F"]
    assert abs(r[MP]["slices"][-1]["X_right"]) < 1e-9  # 반정현 f = 0 (끝)


def test_safety_map_and_reference_only():
    m = {"materials": mats(("s", 17.64, 9.8, 10)), "regions": [{"material_id": "s", "polygon": [[0, 0], [25, 0], [25, 10], [15, 10], [5, 5], [0, 5]]}]}
    r = analyze(m, direction=-1)
    smap = r["safety_map"]["surfaces"]
    assert 0 < len(smap) <= 50
    F = [s["F"] for s in smap]
    assert F == sorted(F) and F[0] >= 0.97 * circ(r, "bishop")
    assert all(len(s["points"]) == 25 and s["points"][0][0] < s["points"][-1][0] for s in smap)
    # 비원호 Fellenius · Bishop = 참고 칸 · 최소 F 는 원호에서만
    for k in ("ordinary", "bishop"):
        res = r["methods"][k]
        assert res["noncircular"]["reference_only"] is True
        assert res["min"]["kind"] == "circle"
    assert r["methods"]["spencer"]["min"]["F"] <= circ(r, "spencer") + 1e-12


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