Merge remote-tracking branch 'origin/sub_laptop_5' into main_laptop_1

2026-10-25 20:45:41 +09:00
parent 8f8d1d0639 69c683b509
commit 630e62455b
4 changed files with 101 additions and 15 deletions
@@ -28,6 +28,8 @@ from E03_StabSection.E03_StabSection_Types import DEFAULT_MATERIALS, DEFAULT_ROA
 from E04_Stab3D.E04_Stab3D_Engine_Column import Strata

 DEFAULT_SOIL_THICKNESS_M = 3.0
+CLOSE_STEP_M = 0.5  # 안 닫힌 사면 늘리기 간격 · 상한
+CLOSE_MAX_M = 60.0
 ROCK_KINDS = ("ripping_rock", "blasting_rock")
 #: 층 → 물성 id(E03 기본표) · 못 뚫는 층
 LAYER_MATERIALS = ("fill", "weathered_soil", "weathered_rock", "soft_rock")
@@ -56,14 +58,42 @@ def station_from_b06(cross: dict, design: dict) -> dict:
 def _delta(
    st: dict, center: np.ndarray, left: np.ndarray, ground_at
 ) -> tuple[np.ndarray, np.ndarray]:
-    """설계선 offset 위 δ = 설계 − 원지형(쌍선형) — B06 이 반폭 밖으로 늘린 설계선도 받음."""
+    """설계선 offset 위 δ = 설계 − 원지형(쌍선형) · 양 끝 밖 δ(0 · 자료 끝이면 끝 값).
+
+    B06 이 반폭 밖으로 늘린 설계선도 받음.
+
+    끝에서 사면이 안 닫혔으면(B06 slope_unclosed · |δ| > 5 cm) 끝 기울기로 원지형과 만날 때까지
+    늘림(CLOSE_MAX_M 까지) — 안 늘리면 띠 끝에 절벽이 생겨 가짜 임계 활동면이 됨.
+    """
    d = np.asarray(st["design"], float)
    d = d[np.argsort(d[:, 0], kind="stable")]
    pts = center[None] + d[:, :1] * left[None]
-    zg = ground_at(pts)
-    delta = d[:, 1] - zg
+    delta = d[:, 1] - ground_at(pts)
    ok = np.isfinite(delta)
-    return d[ok, 0], delta[ok]
+    o, z, delta = d[ok, 0], d[ok, 1], delta[ok]
+    parts_o, parts_d = [o], [delta]
+    hold = {-1: False, 1: False}
+    for side in (-1, 1):
+        k0, k1 = (0, 1) if side < 0 else (-1, -2)
+        if len(o) < 2 or abs(delta[k0]) <= 0.05 or o[k0] == o[k1]:
+            continue
+        slope = (z[k0] - z[k1]) / (o[k0] - o[k1])
+        ext = o[k0] + side * np.arange(CLOSE_STEP_M, CLOSE_MAX_M + 1e-9, CLOSE_STEP_M)
+        de = z[k0] + slope * (ext - o[k0]) - ground_at(center[None] + ext[:, None] * left[None])
+        cross = np.flatnonzero(~np.isfinite(de) | (np.sign(de) != np.sign(delta[k0])))
+        n = cross[0] if len(cross) else len(ext)
+        if n < len(ext) and np.isfinite(de[n]):
+            ext, de = ext[: n + 1], np.append(de[:n], 0.0)  # 원지형과 만남 → 0
+        else:
+            ext, de = ext[:n], de[:n]  # 자료 끝 · 상한 — 끝 값을 이어 감(build_grid hold)
+            hold[side] = True
+        parts_o.append(ext)
+        parts_d.append(de)
+    o = np.concatenate(parts_o)
+    delta = np.concatenate(parts_d)
+    order = np.argsort(o, kind="stable")
+    o, delta = o[order], delta[order]
+    return o, delta, (delta[0] if hold[-1] else 0.0), (delta[-1] if hold[1] else 0.0)


 def _segment_map(xx, yy, c0, c1, l0, l1):
@@ -116,14 +146,15 @@ def build_grid(
    if y[0] > y[-1]:
        y, z, valid = y[::-1], z[::-1], valid[::-1]
    z_orig = np.where(valid, z, np.nan)
+    # δ 는 원래 DEM 에서 뜸(B06 지반선과 같은 면) — 다시 뜬 격자로 뜨면 cm 잡음이 끝 판정을 흔듦
+    interp = RegularGridInterpolator((y, x), z_orig, bounds_error=False, fill_value=np.nan)
    if cell_m:
        x0, x1, y0, y1 = bbox or (x[0], x[-1], y[0], y[-1])
        x0, y0 = max(x0, x[0]), max(y0, y[0])
        xs = x0 + np.arange(int((min(x1, x[-1]) - x0) // cell_m) + 1) * cell_m
        ys = y0 + np.arange(int((min(y1, y[-1]) - y0) // cell_m) + 1) * cell_m
-        f = RegularGridInterpolator((y, x), z_orig, bounds_error=False, fill_value=np.nan)
        gy_, gx_ = np.meshgrid(ys, xs, indexing="ij")
-        z_orig = f(np.column_stack([gy_.ravel(), gx_.ravel()])).reshape(gy_.shape)
+        z_orig = interp(np.column_stack([gy_.ravel(), gx_.ravel()])).reshape(gy_.shape)
        x, y = xs, ys
    elif bbox is not None:
        sx = (x >= bbox[0]) & (x <= bbox[1])
@@ -135,7 +166,6 @@ def build_grid(
    centers = np.array([s["center"] for s in sts], float)
    lefts = np.array([s["left"] for s in sts], float)
    lefts /= np.hypot(lefts[:, 0], lefts[:, 1])[:, None]
-    interp = RegularGridInterpolator((y, x), z_orig, bounds_error=False, fill_value=np.nan)

    def ground_at(pts: np.ndarray) -> np.ndarray:
        return interp(np.column_stack([pts[:, 1], pts[:, 0]]))
@@ -155,9 +185,12 @@ def build_grid(
    dz = np.zeros(z.shape)
    road = np.zeros(z.shape, bool)
    for i in range(len(sts) - 1):
-        o0, d0 = profiles[i]
-        o1, d1 = profiles[i + 1]
+        o0, d0, l0_, r0_ = profiles[i]
+        o1, d1, l1_, r1_ = profiles[i + 1]
        lo, hi = min(o0[0], o1[0]), max(o0[-1], o1[-1])
+        # 자료 끝에서 멈춘 쪽은 띠를 더 넓혀 끝 값을 이어 감(빈칸 셀은 합성도 빈칸)
+        lo -= CLOSE_MAX_M if (l0_ or l1_) else 0.0
+        hi += CLOSE_MAX_M if (r0_ or r1_) else 0.0
        reach = max(abs(lo), abs(hi))
        corners = np.array(
            [centers[i] + s * reach * lefts[i] for s in (-1, 1)]
@@ -176,7 +209,7 @@ def build_grid(
        if not better.any():
            continue
        tb, ob = t[better], o[better]
-        val = (1 - tb) * np.interp(ob, o0, d0, 0.0, 0.0) + tb * np.interp(ob, o1, d1, 0.0, 0.0)
+        val = (1 - tb) * np.interp(ob, o0, d0, l0_, r0_) + tb * np.interp(ob, o1, d1, l1_, r1_)
        for arr, v in ((seg, i), (t_map, tb), (o_map, ob), (dz, val)):
            blk = arr[sub]
            blk[better] = v
@@ -54,10 +54,11 @@ class SearchParams:
    v_max: float = 10_000.0
    a_min: float = 0.0
    a_max: float = math.inf
-    ratios: tuple[float, ...] = (1.0, 2.0)  # 1 = 구 · 그 밖 = 축 방향 반축 배수(타원체)
+    ratios: tuple[float, ...] = (1.0, 2.0)  # 구(늘 봄) + 타원체 축 방향 배수(구 임계 노드에서)
    sub: int = 2  # 셀 안 표본 sub×sub(부분 기둥)
    refine_rounds: int = 4  # 가장 고운 간격에서 더 내려가 보는 횟수 상한
-    refine_below: float = 2.0  # 고운 탐색은 이 F 아래로 고친 노드만(높은 곳 지도는 거친 값)
+    # 고운 탐색 · 타원체는 F < max(이 값, 전역 최소 × 1.2) 로 고친 노드만(높은 곳 지도는 거친 값)
+    refine_below: float = 2.0
    groundwater: dict | None = None
    seismic_kh: float = 0.0

@@ -275,7 +276,7 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
        z = strata.z_top[cy, cx] + p.h_min + ih * p.h_step
        return np.column_stack([x, y, z]), cy, cx, p.h_min + ih * p.h_step

-    def run(ids: np.ndarray) -> set[int]:
+    def run(ids: np.ndarray, ratios: tuple[float, ...] = (1.0,)) -> set[int]:
        """노드들의 모든 반지름 · 비 — 셀 최소를 고친 노드 번호를 돌려줌."""
        improved: set[int] = set()
        xyz0, cy0, cx0, _ = node_xyz(ids)
@@ -289,7 +290,7 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
        dirs = np.column_stack([fall_x[cy, cx], fall_y[cy, cx]])
        r_lo = np.maximum(p.dr, np.floor(0.5 * h / p.dr) * p.dr)
        r_hi = h + p.max_depth
-        for ratio in p.ratios:
+        for ratio in ratios:
            done = np.zeros(len(ids), bool)
            for r in np.arange(r_lo.min(), r_hi.max() + 1e-9, p.dr):
                live = np.flatnonzero(~done & (r_lo <= r + 1e-9) & (r <= r_hi + 1e-9))
@@ -330,7 +331,7 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
                        volume,
                        method,
                        node_map,
-                        p.refine_below,
+                        max(p.refine_below, 1.2 * best["governing"]),
                    )
                    k = (
                        int(np.nanargmin(np.where(np.isfinite(f), f, np.inf)))
@@ -388,6 +389,13 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
        seen.update(cand)
        seeds = run(np.fromiter(cand, np.int64))

+    # 타원체 — 구 탐색이 찾은 임계 노드(F < refine_below)에서만(시간 1/3 · 노선 1 km 실측)
+    ell = tuple(r for r in p.ratios if r != 1.0)
+    limit = max(p.refine_below, 1.2 * best["governing"])
+    crit = np.unique(node_map[(node_map >= 0) & (f_min < limit)])
+    if ell and len(crit):
+        run(crit, ell)
+
    f_out = np.where(np.isfinite(f_min), f_min, np.nan)
    return SearchResult(
        f_min=f_out,
@@ -103,3 +103,21 @@ def test_rock_boundary_and_layers():
    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]
@@ -3,6 +3,7 @@
 Reid 외(2015) Scoops3D 표 6.8 · 6.11(CLARA-W 와 같은 구): H 10 m 제방
 - 30° · λ 0.125(γ 19 · φ 20° · c 8.64 kPa) → Bishop 1.32 · 부피 1,124 m³
 - Donald & Giam 1a 2:1 · γ 20 · φ 19.6° · c 3 kPa → 1.04
+- 표 6.7 원뿔 30° · r_u 0.4 → 1.626
 탐색 격자는 시험 시간 때문에 문헌 임계 중심 둘레로 좁힘.
 """

@@ -58,3 +59,29 @@ def test_donald_giam_1a():
        ratios=(1.0,),
    )
    assert search(st, p).best["governing"] == pytest.approx(1.04, rel=0.03)
+
+
+def test_wet_cone_pore_pressure_ratio():
+    """Scoops3D 표 6.7 · 원뿔 30° · λ 0.125(γ 20 · φ 30°) · r_u 0.4 → Bishop 1.626(도표 1.644)."""
+    h = 10.0
+    r0 = h / np.tan(np.radians(30))
+    x = np.arange(-r0 - 15, r0 + 15 + 1e-9, 0.5)
+    xx, yy = np.meshgrid(x, x)
+    z = np.clip(h - np.hypot(xx, yy) * np.tan(np.radians(30)), 0, h)
+    c = 0.125 * 20 * h * np.tan(np.radians(30))
+    mat = {"id": "m", "name": "m", "gamma_t": 20, "gamma_sat": 20, "c": c, "phi": 30}
+    st = Strata(x, x, z, np.zeros((0,) + z.shape), [mat])
+    p = SearchParams(
+        box=(8, r0 + 4, -1, 1),
+        spacing=1.0,
+        h_min=1,
+        h_max=21,
+        h_step=2,
+        coarse=4,
+        max_depth=15,
+        v_min=20,
+        v_max=3e4,
+        ratios=(1.0,),
+        groundwater={"mode": "ru", "ru": 0.4},
+    )
+    assert search(st, p).best["governing"] == pytest.approx(1.626, rel=0.03)