Aislo/resources/tester/test_533_1_stab3d_grid.py

"""533-1 E04 합성 DEM · 지층 격자 — 측점 위 단면 = B06 설계선 · 띠 밖 = 원지형 · 보간 · 지층."""

import numpy as np
import pytest

from E04_Stab3D.E04_Stab3D_Engine_Grid import _segment_map, build_grid, strata_from_grid


def _dem(dx=0.5):
    x = np.arange(-30, 90 + 1e-9, dx)
    y = np.arange(-40, 40 + 1e-9, dx)
    xx, yy = np.meshgrid(x, y)
    z = 100.0 + 0.3 * yy + 0.02 * xx  # 왼쪽(+y)이 오르막
    return {"x": x, "y": y, "z": z, "valid_mask": np.ones(z.shape, bool)}, xx, yy


def _station(dem, cx, depth_cut, ground_type="soil", rock=None):
    """x=cx · 왼쪽 = +y · 노면 ±2 m 를 depth_cut 낮춘 설계선(사면 1:1 로 지반에 닿음)."""
    offs = np.arange(-20, 20.01, 0.5)
    g = 100.0 + 0.3 * offs + 0.02 * cx
    design = g.copy()
    road = np.abs(offs) <= 2.0
    design[road] = 100.0 + 0.02 * cx - depth_cut
    for k in np.flatnonzero(~road):
        edge = 2.0 * np.sign(offs[k])
        design[k] = (
            min(g[k], 100.0 + 0.02 * cx - depth_cut + abs(offs[k] - edge))
            if depth_cut > 0
            else g[k]
        )
    return {
        "chainage_m": float(cx),
        "center": [float(cx), 0.0],
        "left": [0.0, 1.0],
        "ground": np.column_stack([offs, g]).tolist(),
        "design": np.column_stack([offs, design]).tolist(),
        "ground_type": ground_type,
        "rock_depth_m": rock,
        "road": [-2.0, 2.0],
    }


def _grid():
    dem, xx, yy = _dem()
    sts = [
        _station(dem, 0, 1.0),
        _station(dem, 20, 3.0, "blasting_rock", 0.5),
        _station(dem, 40, 3.0, "ripping_rock", 1.0),
    ]
    return dem, build_grid(dem, sts, cell_m=0.5), xx, yy, sts


def test_station_section_equals_design_and_outside_unchanged():
    dem, g, xx, yy, sts = _grid()
    for st in sts:
        col = np.isclose(g["x"], st["center"][0])
        d = np.asarray(st["design"])
        z_line = np.interp(g["y"], d[:, 0], d[:, 1])
        inside = np.abs(g["y"]) <= 20
        assert np.abs(g["z"][inside, col][:, None].ravel() - z_line[inside]).max() < 1e-9
    out = ~g["in_band"]
    assert np.array_equal(g["z"][out], g["z_orig"][out])
    assert np.array_equal(g["z"][:, g["x"] < 0], dem["z"][:, g["x"] < 0])  # 첫 측점 앞 = 원지형


def test_between_stations_linear_delta():
    _, g, xx, yy, sts = _grid()
    j = np.flatnonzero(np.isclose(g["x"], 10.0))[0]
    i = np.flatnonzero(np.isclose(g["y"], 0.0))[0]
    # 노면 가운데 δ: 측점 0 −1 m · 측점 20 −3 m → 중간 −2 m
    assert g["dz"][i, j] == pytest.approx(-2.0, abs=1e-9)
    assert g["q"][i, j] == pytest.approx(10.0)


def test_curved_segment_map_round_trip():
    c0, c1 = np.array([0.0, 0.0]), np.array([19.0, 6.0])
    l0 = np.array([0.0, 1.0])
    a = np.radians(35.0)
    l1 = np.array([-np.sin(a), np.cos(a)])
    t_true = np.array([0.0, 0.3, 0.75, 1.0])
    o_true = np.array([-15.0, 4.0, 12.0, -7.0])
    lt = l0[None] + t_true[:, None] * (l1 - l0)[None]
    lt /= np.hypot(lt[:, 0], lt[:, 1])[:, None]
    p = c0[None] + t_true[:, None] * (c1 - c0)[None] + o_true[:, None] * lt
    t, o = _segment_map(p[:, 0], p[:, 1], c0, c1, l0, l1)
    assert np.allclose(t, t_true, atol=1e-9) and np.allclose(o, o_true, atol=1e-9)


def test_rock_boundary_and_layers():
    _, g, xx, yy, sts = _grid()
    j0 = np.flatnonzero(np.isclose(g["x"], 20.0))[0]
    jm = np.flatnonzero(np.isclose(g["x"], 30.0))[0]
    i = np.flatnonzero(np.isclose(g["y"], 10.0))[0]
    assert g["soil_thickness"][i, j0] == pytest.approx(0.5)
    assert g["soil_thickness"][i, jm] == pytest.approx(0.75)  # 0.5 → 1.0 가운데
    far = np.flatnonzero(np.isclose(g["y"], 35.0))[0]
    assert not g["in_band"][far, j0]
    assert g["soil_thickness"][far, j0] == pytest.approx(0.5)  # 띠 밖 = 가까운 측점
    assert g["soil_thickness"][i, 0] == pytest.approx(3.0)  # 토사 측점 = 기본 두께
    st = strata_from_grid(g)
    assert g["rock_kind"][i, j0] == 1  # 발파암 → 못 뚫는 층 윗면 = 암 윗면(깎기로 드러나면 지표)
    assert st.hard_top[i, j0] == pytest.approx(min(g["rock_top"][i, j0], g["z"][i, j0]))
    assert np.isneginf(st.hard_top[i, jm + 20]) or st.hard_top[i, jm + 20] < -1e5  # 리핑암 = 뚫림
    ii = np.flatnonzero(np.isclose(g["y"], 0.0))[0]
    assert st.bottoms[0, ii, j0] == pytest.approx(g["z"][ii, j0])  # 깎기 = 성토 두께 0


def test_unclosed_slope_extended_to_ground():
    """B06 사면이 반폭 안에서 안 닫힘(끝 δ ≠ 0) → 끝 기울기로 원지형까지 늘려 띠 끝 절벽 없앰."""
    dem, xx, yy = _dem()
    sts = [_station(dem, cx, 1.0) for cx in (0, 20)]
    for st in sts:  # 오른쪽(−y) 끝 5 m 를 1:1 성토로 바꿔 −20 m 에서 2 m 떠 있게 함
        d = np.asarray(st["design"])
        g = np.asarray(st["ground"])
        tail = d[:, 0] <= -15
        d[tail, 1] = g[tail, 1] + 2.0 - (d[tail, 0] + 20) * 0.4
        st["design"] = d.tolist()
    g = build_grid(dem, sts, cell_m=0.5)
    j = np.flatnonzero(np.isclose(g["x"], 10.0))[0]
    col = g["dz"][:, j]
    # 늘린 사면이 원지형과 만나는 곳 너머는 0 · 이웃 칸 사이 단차 1 m 아래
    assert np.abs(np.diff(col[np.isfinite(col)])).max() < 1.0
    assert g["in_band"][np.flatnonzero(np.isclose(g["y"], -21.0))[0], j]