From 69c683b5098a85fcd1f6de23b3cffb277c410066 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Sun, 25 Oct 2026 20:45:19 +0900 Subject: [PATCH] =?UTF-8?q?fix(533-1=20=C2=B7=20533-2):=20=EC=95=88=20?= =?UTF-8?q?=EB=8B=AB=ED=9E=8C=20=EC=82=AC=EB=A9=B4=20=EC=97=B0=EC=9E=A5=20?= =?UTF-8?q?=C2=B7=20=CE=B4=20=EB=8A=94=20=EC=9B=90=20DEM=20=EC=97=90?= =?UTF-8?q?=EC=84=9C=20=C2=B7=20=ED=83=80=EC=9B=90=EC=B2=B4=EB=8A=94=20?= =?UTF-8?q?=EA=B5=AC=20=EC=9E=84=EA=B3=84=20=EB=85=B8=EB=93=9C=EC=97=90?= =?UTF-8?q?=EC=84=9C=20=C2=B7=20=EC=9B=90=EB=BF=94=20r=5Fu=20=EC=8B=9C?= =?UTF-8?q?=ED=97=98?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Grid: B06 사면이 반폭 안에서 안 닫히면 끝 기울기로 원지형까지 늘림(자료 끝이면 끝 값 이어 감) — 띠 끝 7 m 가짜 절벽 없앰 - Grid: δ = 설계 − 원 DEM(다시 뜬 1 m 격자의 cm 잡음이 끝 판정을 흔들던 것) - Search: 타원체는 구 탐색 임계 노드에서만(1 km 822 s → 287 s) · 다듬기 문턱 = max(2.0, 전역 최소 × 1.2) - 시험: 안 닫힌 사면 연장 · Scoops3D 원뿔 r_u 0.4(1.626) Co-Authored-By: Claude Opus 5.5 --- E04_Stab3D/E04_Stab3D_Engine_Grid.py | 53 ++++++++++++++++---- E04_Stab3D/E04_Stab3D_Engine_Search.py | 18 +++++-- resources/tester/test_533_1_stab3d_grid.py | 18 +++++++ resources/tester/test_533_2_stab3d_search.py | 27 ++++++++++ 4 files changed, 101 insertions(+), 15 deletions(-) diff --git a/E04_Stab3D/E04_Stab3D_Engine_Grid.py b/E04_Stab3D/E04_Stab3D_Engine_Grid.py index 6c5217624..18277cb89 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Grid.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Grid.py @@ -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 diff --git a/E04_Stab3D/E04_Stab3D_Engine_Search.py b/E04_Stab3D/E04_Stab3D_Engine_Search.py index 2ce1b88be..8a2c1426d 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Search.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Search.py @@ -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, diff --git a/resources/tester/test_533_1_stab3d_grid.py b/resources/tester/test_533_1_stab3d_grid.py index 8e4ced245..639c2aefd 100644 --- a/resources/tester/test_533_1_stab3d_grid.py +++ b/resources/tester/test_533_1_stab3d_grid.py @@ -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] diff --git a/resources/tester/test_533_2_stab3d_search.py b/resources/tester/test_533_2_stab3d_search.py index 1ef54b604..1e3c408bd 100644 --- a/resources/tester/test_533_2_stab3d_search.py +++ b/resources/tester/test_533_2_stab3d_search.py @@ -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)