eomsangdon/Aislo
1
0
Fork 0
Code Issues Pull Requests 2 Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
eomsangdon
2026-10-26 04:23:48 +09:00
parent ab8013fbf4 e966a37b5f
commit b546569ed1
3 changed files with 176 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files
E04_Stab3D/E04_Stab3D_Engine_Search.py
+139 -58
View File
@@ -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,
resources/tester/test_533_2_stab3d_column.py
+16
View File
@@ -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)
resources/tester/test_533_2_stab3d_search.py
+21
View File
@@ -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(