From e966a37b5f65c1c3e7d3416f8947dbd8501d2b2c Mon Sep 17 00:00:00 2001 From: umsangdon Date: Mon, 26 Oct 2026 04:22:26 +0900 Subject: [PATCH] =?UTF-8?q?fix(533-2=202=EC=B0=A8):=20E04=20=EA=B8=89?= =?UTF-8?q?=EA=B2=BD=EC=82=AC=20=EC=A0=84=EC=97=AD=20=EC=B5=9C=EC=86=8C=20?= =?UTF-8?q?=E2=80=94=20=EC=A0=90=20=EB=8B=A8=EC=9C=84=20=EB=8D=A9=EC=9D=B4?= =?UTF-8?q?=20=EB=82=98=EB=88=94=20=C2=B7=20=EC=A0=84=EC=97=AD=20=EC=B5=9C?= =?UTF-8?q?=EC=86=8C=20=EB=8B=A4=EB=93=AC=EA=B8=B0=20=C2=B7=20Zhang=20?= =?UTF-8?q?=EA=B2=BD=EC=9A=B0=203=20=EC=8B=9C=ED=97=98?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Search: 활동면이 지표를 자른 점을 셀 안 점 단위로 덩이 나눔 — 발끝을 살짝 비켜 지나는 구의 사면 몸통 · 발끝 앞 접시가 따로(셀 중심 줄 꺾임만 검사 · 지표 쌍선형 + 활동면 볼록) - Search: 전역 최소 면을 찾은 노드 · 반지름 한 눈금 안에서 Nelder-Mead 로 다듬음(box · node_mask 안) — 셀 지도도 그 면으로 고침 - Scoops3D 60° 0.771~0.776 → 0.754(0.75 · +0.5 %) · 30° 1.325 → 1.3197(1.32) - 검증_sub5 1 km 띠: 전역 최소 1.185 → 1.174(−0.9 %) · p1/p5/p50 같음 · 셀 9.7만 중 168 이 1 % 넘게 바뀜 · 12 프로세스 28 s → 45 s - 시험: 60° 제방 ±1 % · Zhang 경우 3(약층 + 수위선 · FK 1977 그림 8) Bishop 1.524 = Lam & Fredlund 1.481~1.54 · Gharti 1.49 ±3 % · 경우3/2 비 0.943 Co-Authored-By: Claude Opus 5.5 --- E04_Stab3D/E04_Stab3D_Engine_Search.py | 197 +++++++++++++------ resources/tester/test_533_2_stab3d_column.py | 16 ++ resources/tester/test_533_2_stab3d_search.py | 21 ++ 3 files changed, 176 insertions(+), 58 deletions(-) diff --git a/E04_Stab3D/E04_Stab3D_Engine_Search.py b/E04_Stab3D/E04_Stab3D_Engine_Search.py index 3ef40f964..5a85460b6 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Search.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Search.py @@ -7,7 +7,9 @@ - 타원체 = 활동 단면(d−z)은 원 · 회전축 방향 반축만 ratio 배(Hungr 1989 꼴) — 수직력이 축을 지나 Bishop 그대로. 방향 = 노드 지표 낙하 방향(5×5 평활). 구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1). -- 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈(ndimage.label · 8 이웃). +- 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈 — 셀 안 점 단위로 가름(_split_pieces · + 발끝을 살짝 비켜 지나는 구의 사면 몸통 · 발끝 앞 접시가 따로). +- 전역 최소 면은 마지막에 중심 · 반지름을 눈금 사이로 다듬음(_polish). - 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음). - 후보 덩이는 프로세스로 나눠 셈 · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값. @@ -172,9 +174,9 @@ def trial_columns( # 2) 부분 기둥 — sub×sub 점 fr = ((np.arange(sub) + 0.5) / sub - 0.5) if sub > 1 else np.zeros(1) - nf = nw * sub # 점 격자 한 변 - win_y, win_x = (ci // nw) % nw, ci % nw - pts: dict[str, list] = {k: [] for k in ("c", "f", "depth", "zs", "gx", "gy", "px", "py")} + pts: dict[str, list] = {k: [] for k in ("c", "depth", "zs", "gx", "gy", "px", "py")} + slot = np.full((sub, sub, len(ci)), -1, np.int32) # (ka, kb, 셀) → 점 번호 + n_pt = 0 for ka, fy in enumerate(fr): for kb, fx in enumerate(fr): xs = strata.x[0] + (cx + fx) * dx @@ -184,7 +186,8 @@ def trial_columns( z_b = np.maximum(z_s, hard) k = np.flatnonzero((rem > 0) & (z_b < zg)) pts["c"].append(k) - pts["f"].append((tt[k] * nf + win_y[k] * sub + ka) * nf + win_x[k] * sub + kb) + slot[ka, kb, k] = n_pt + np.arange(len(k)) + n_pt += len(k) for key, val in ( ("depth", zg - z_b), ("zs", z_s), @@ -195,10 +198,11 @@ def trial_columns( ): pts[key].append(val[k]) pt = {key: np.concatenate(v) for key, v in pts.items()} - n_pt = len(pt["c"]) if not n_pt: return empty - lab_pt, n_lab = _split_pieces(strata, pt, tt, cy, cx, hard, n_t * nf * nf, nf, sub, surface) + lab_pt, n_lab = _split_pieces( + strata, slot, ci, n_t * nw * nw, nw, tt, cy, cx, hard, sub, surface, n_pt + ) # 기둥 = (셀, 덩이) 짝 — 두 덩이가 한 셀을 나누면 기둥 둘 key, inv = np.unique(pt["c"] * n_lab + lab_pt, return_inverse=True) cell_k, lab = key // n_lab, key % n_lab + 1 # 덩이 번호 1.. (0 = 없음) @@ -279,44 +283,54 @@ def _part(job: tuple, k: dict | None = None) -> tuple: return light, trial, big, res -def _split_pieces(strata, pt, tt, cy, cx, hard, n_fine, nf, sub, surface): - """활동면이 지표를 자른 점들을 덩이로 — 이웃 점(4 방향) 사이 가운데 점도 땅속일 때만 이음. +def _split_pieces(strata, slot, ci, n_win, nw, tt, cy, cx, hard, sub, surface, n_pt): + """활동면이 지표를 자른 점들을 덩이로 — 이웃 점(4 방향)이 둘 다 땅속이고 사이도 땅속이면 이음. 한 셀 안에서도 덩이가 갈림(발끝을 살짝 비켜 지나는 구 = 사면 몸통 · 발끝 앞 접시 따로). - 반환 (점별 덩이 번호 0.., 덩이 수). + 지표는 셀 중심 사이 쌍선형(축 방향 직선) · 활동면은 아래로 볼록 → 두 점 사이 깊이는 + 지표 꺾임(셀 중심 줄)에서만 끝점보다 낮아짐 → 셀 중심 줄을 사이에 둔 짝만 그 줄에서 봄. + slot = (ka, kb, 셀) → 점 번호(−1 없음) · ci = 창 칸 번호 · n_win = 창 칸 수. + 반환 (점별 덩이 번호, 수). """ from scipy.sparse import coo_matrix from scipy.sparse.csgraph import connected_components - n_pt = len(pt["c"]) - where = np.full(n_fine, -1, np.int32) # 점 격자 칸 → 점 번호 - where[pt["f"]] = np.arange(n_pt, dtype=np.int32) fr = ((np.arange(sub) + 0.5) / sub - 0.5) if sub > 1 else np.zeros(1) - half = 0.5 / sub rows, cols_ = [], [] - for step, axis in ((1, "x"), (nf, "y")): - pos = pt["f"] % nf if axis == "x" else (pt["f"] // nf) % nf - i = np.flatnonzero(pos < nf - 1) # 줄 · 창 끝 넘김 막음 - q = where[pt["f"][i] + step] - i, q = i[q >= 0], q[q >= 0] - c = pt["c"][i] - # 가운데 점 — 점 i 의 셀에서 (fy, fx) + half - ka = ((pt["f"][i] // nf) % nf) % sub - kb = (pt["f"][i] % nf) % sub - mx = pt["px"][i] + (half * strata.dx if axis == "x" else 0.0) - my = pt["py"][i] + (half * strata.dx if axis == "y" else 0.0) - rem, z_s, _, _ = surface(tt[c], mx, my) - zg = np.empty(len(i)) - for va in range(sub): - for vb in range(sub): - m = (ka == va) & (kb == vb) - if m.any(): - fy = fr[va] + (half if axis == "y" else 0.0) - fx = fr[vb] + (half if axis == "x" else 0.0) - zg[m] = strata.ground_at(cy[c[m]], cx[c[m]], fy, fx) - ok = (rem > 0) & (np.maximum(z_s, hard[c]) < zg) - rows.append(i[ok]) - cols_.append(q[ok]) + + def link(a, b): + both = (a >= 0) & (b >= 0) + rows.append(a[both]) + cols_.append(b[both]) + + # 셀 안 — 꺾임(f = 0)을 사이에 둔 짝은 꺾임 점에서 땅속인지 봄 + for k in range(sub - 1): + for m in range(sub): + for axis in ("x", "y"): + a = slot[m, k] if axis == "x" else slot[k, m] + b = slot[m, k + 1] if axis == "x" else slot[k + 1, m] + c = np.flatnonzero((a >= 0) & (b >= 0)) + if fr[k] < 0 < fr[k + 1]: + fy, fx = (fr[m], 0.0) if axis == "x" else (0.0, fr[m]) + xs = strata.x[0] + (cx[c] + fx) * strata.dx + ys = strata.y[0] + (cy[c] + fy) * strata.dx + rem, z_s, _, _ = surface(tt[c], xs, ys) + zg = strata.ground_at(cy[c], cx[c], fy, fx) + c = c[(rem > 0) & (np.maximum(z_s, hard[c]) < zg)] + rows.append(a[c]) + cols_.append(b[c]) + # 셀 사이 — 이웃 셀(창 안 · 줄 넘김 막음)의 맞닿은 점끼리 + cell_at = np.full(n_win + nw, -1, np.int32) # 창 칸 → 셀 번호(끝 넘침 칸 −1) + cell_at[ci] = np.arange(len(ci), dtype=np.int32) + for step, pos in ((1, ci % nw), (nw, (ci // nw) % nw)): + j = cell_at[ci + step] + nb = np.flatnonzero((j >= 0) & (pos < nw - 1)) + j = j[nb] + for m in range(sub): + if step == 1: + link(slot[m, sub - 1, nb], slot[m, 0, j]) + else: + link(slot[sub - 1, m, nb], slot[0, m, j]) r_ = np.concatenate(rows) c_ = np.concatenate(cols_) graph = coo_matrix((np.ones(len(r_), np.int8), (r_, c_)), shape=(n_pt, n_pt)) @@ -379,14 +393,13 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: method = np.full((ny, nx), -1, np.int8) node_map = np.full((ny, nx), -1, np.int64) best = {"governing": np.inf} + at_best: dict = {} # 전역 최소의 노드 번호 · 방향(다듬기용 · 결과에 안 냄) seen: set[int] = set() counters = {"trials": 0, "evaluated": 0} stage = [0, 40, "거친 탐색"] # 진행률 — 시작 % · 폭 % · 이름 def tick(done: float) -> None: if p.should_stop and p.should_stop(): - if pool: - pool.shutdown(cancel_futures=True) raise InterruptedError if p.on_progress: p.on_progress(int(stage[0] + stage[1] * done), str(stage[2])) @@ -463,25 +476,11 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: thr, ref, ) - k = ( - int(np.nanargmin(np.where(np.isfinite(f), f, np.inf))) - if np.isfinite(f).any() - else -1 - ) - if k >= 0 and f[k] < best["governing"]: - best.clear() - best.update( - { - key: (v[k].tolist() if isinstance(v, np.ndarray) else v) - for key, v in res.items() - } - ) - best.update( - governing=float(f[k]), - method=int(meth[k]), - center=xyz[part][trial[k]].tolist(), - radius=float(r), - ratio=float(ratio), + k = _set_best(best, res, f, meth, xyz[part][trial], r, ratio) + if k >= 0: + at_best.update( + node=int(ids[part][trial[k]]), + direction=dirs[part][trial[k]] if ratio != 1.0 else None, ) return improved @@ -533,6 +532,11 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: stage[:] = [85, 14, "타원체"] run(crit, ell) + if at_best: + stage[:] = [99, 1, "다듬기"] + tick(0.0) + _polish(strata, p, best, at_best, f_min, volume, method, node_map) + finally: if pool: # 끊기(InterruptedError) · 오류에도 일꾼을 닫음 pool.shutdown(cancel_futures=True) @@ -549,6 +553,83 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: ) +def _set_best(best: dict, res: dict, f, meth, centers, r: float, ratio: float) -> int: + """덩이들 가운데 best 보다 낮은 것이 있으면 best 를 바꾸고 그 번호 · 없으면 −1.""" + if not np.isfinite(f).any(): + return -1 + k = int(np.nanargmin(np.where(np.isfinite(f), f, np.inf))) + if not f[k] < best["governing"]: + return -1 + best.clear() + best.update( + {key: (v[k].tolist() if isinstance(v, np.ndarray) else v) for key, v in res.items()} + ) + best.update( + governing=float(f[k]), + method=int(meth[k]), + center=np.asarray(centers[k], float).tolist(), + radius=float(r), + ratio=float(ratio), + ) + return k + + +def _polish(strata, p, best, at_best, f_min, volume, method, node_map) -> None: + """전역 최소 다듬기 — 그 면의 중심 · 반지름을 노드 · 반지름 눈금 사이로 옮겨 봄(Nelder-Mead). + + 급경사 전역 최소는 발끝에 거의 닿는 면 — 눈금(노드 간격 · dr)에 걸린 값이 연속 최소보다 + 1~2 % 높음(60° 제방 1 m 격자 0.7595 → 0.7519 · Scoops3D 0.75). 비 · 방향은 그대로. + """ + from scipy.optimize import minimize + + ratio = best["ratio"] + d = at_best["direction"] + d = None if d is None else np.asarray(d, float)[None] + + def trial(v): + cols, tr, _ = trial_columns( + strata, v[None, :3], v[3], ratio, d, p.v_min, p.v_max, p.a_min, p.a_max, p.sub + ) + if not len(tr): + return None + res = solve(strata, cols, p.groundwater, p.seismic_kh) + return (cols, tr, res, *governing(res)) + + def objective(v): + if p.node_mask is not None: # 둘레 네 셀이 모두 노드 칸이어야(가장자리 밖으로 안 샘) + fy = (v[1] - strata.y[0]) / strata.dx + fx = (v[0] - strata.x[0]) / strata.dx + iy = np.clip([math.floor(fy), math.ceil(fy)], 0, len(strata.y) - 1) + ix = np.clip([math.floor(fx), math.ceil(fx)], 0, len(strata.x) - 1) + if not p.node_mask[np.ix_(iy, ix)].all(): + return np.inf + out = trial(v) + if out is None or not np.isfinite(out[3]).any(): + return np.inf + return float(np.nanmin(out[3])) + + v0 = np.r_[best["center"], best["radius"]] + step = np.array([p.spacing, p.spacing, p.h_step, p.dr]) / 2 + simplex = np.vstack([v0, v0 + np.diag(step)]) + # 찾은 노드 · 반지름에서 한 눈금 안(그 사이를 놓친 것만 고침) · box 안 + lo, hi = v0 - 2 * step, v0 + 2 * step + if p.box is not None: + lo[:2] = np.maximum(lo[:2], (p.box[0], p.box[2])) + hi[:2] = np.minimum(hi[:2], (p.box[1], p.box[3])) + o = minimize( + objective, + v0, + method="Nelder-Mead", + bounds=list(zip(lo, hi)), + options={"initial_simplex": simplex, "xatol": 1e-3, "fatol": 1e-5, "maxiter": 400}, + ) + if not o.fun < best["governing"] * (1 - 1e-6): + return + cols, tr, res, f, meth = trial(o.x) + _update(f, meth, res, cols, tr, np.array([at_best["node"]]), f_min, volume, method, node_map) + _set_best(best, res, f, meth, np.repeat(o.x[None, :3], len(tr), 0), float(o.x[3]), ratio) + + def _update( f, meth, diff --git a/resources/tester/test_533_2_stab3d_column.py b/resources/tester/test_533_2_stab3d_column.py index b85de1337..ceec2043c 100644 --- a/resources/tester/test_533_2_stab3d_column.py +++ b/resources/tester/test_533_2_stab3d_column.py @@ -6,6 +6,9 @@ · 타원체 중심 (6.1, 0, 21.3) Rx=Rz=24.4 · Ry=78.2 · c 29 kPa φ 20° γ 18.8 → Hungr 3D Bishop 2.167 · Zhang 2.122 / 약층(c 0 · φ 10°, 발끝 아래 0.94 m — 그림 화소로 잼) → Hungr 1.620 · Zhang 1.533. +- 경우 3 = 경우 2 + 수위선(Fredlund & Krahn 1977 그림 8 · Slide3 2D Extruded #16 그림 16: + 뒤 경계 42.7 m 에서 발끝 위 6.1 m → 발끝 → 발끝 앞 지표) → Lam & Fredlund 1993 1.481~1.54 + · Gharti 2011 1.49 · Zhang 1.441(경우 2 처럼 Zhang 은 다른 풀이보다 4~6 % 낮음). """ import numpy as np @@ -76,6 +79,19 @@ def test_zhang_case2_weak_layer_truncated(): assert 1.45 < r["janbu_corrected"] < 1.80 # 문헌 폭 1.53~1.67 둘레 +def test_zhang_case3_weak_layer_piezometric_line(): + st = _zhang(True) + xx, _ = np.meshgrid(st.x, st.y) + st.z_water = np.clip(xx * 6.1 / 42.67, 0.0, None) + args = ([6.1, 0.0, 21.3], 24.4, 78.2 / 24.4, [-1.0, 0.0]) + r = evaluate_surface(st, *args, groundwater={"mode": "line"}) + assert 1.481 * 0.99 < r["bishop"] < 1.54 * 1.01 # Lam & Fredlund 1993 네 값 폭 + assert r["bishop"] == pytest.approx(1.49, rel=0.03) # Gharti 2011 + # 수위선이 낮추는 비 = Lam & Fredlund 경우 3/경우 2(0.940~0.951) + dry = evaluate_surface(_zhang(True), *args)["bishop"] + assert 0.93 < r["bishop"] / dry < 0.96 + + def test_pore_pressure_and_seismic(): st, c = _gens(0.04) dry = evaluate_surface(st, c, 1.0) diff --git a/resources/tester/test_533_2_stab3d_search.py b/resources/tester/test_533_2_stab3d_search.py index aed8858d9..7f0ad434e 100644 --- a/resources/tester/test_533_2_stab3d_search.py +++ b/resources/tester/test_533_2_stab3d_search.py @@ -2,6 +2,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³ +- 60° · λ 0.125 → 0.75 — 발끝에 거의 닿는 면(발끝 앞 접시는 따로 덩이 · 마지막 다듬기) ±1 % - Donald & Giam 1a 2:1 · γ 20 · φ 19.6° · c 3 kPa → 1.04 - 표 6.7 원뿔 30° · r_u 0.4 → 1.626 탐색 격자는 시험 시간 때문에 문헌 임계 중심 둘레로 좁힘. @@ -45,6 +46,26 @@ def test_embankment_30deg_global_minimum(): assert np.isfinite(r.volume[np.isfinite(f)]).all() +def test_embankment_60deg_toe_minimum(): + st = _embankment(10 / np.tan(np.radians(60)), 19, 20, 8.64) + p = SearchParams( + box=(-9, -5, -1, 1), + h_min=13, + h_max=17, + h_step=1, + spacing=1.0, + coarse=2, + dr=0.5, + max_depth=20, + v_min=50, + v_max=2e4, + ratios=(1.0,), + ) + r = search(st, p) + assert r.best["governing"] == pytest.approx(0.75, rel=0.01) + assert np.nanmin(r.f_min) == pytest.approx(r.best["governing"]) + + def test_donald_giam_1a(): st = _embankment(20, 20, 19.6, 3.0) p = SearchParams(