Merge remote-tracking branch 'origin/dev' into sub_laptop_4

This commit is contained in:
2026-10-26 03:42:02 +09:00
12 changed files with 410 additions and 143 deletions
+60 -23
View File
@@ -244,6 +244,27 @@ def _finish(res: np.ndarray, drive: np.ndarray) -> np.ndarray:
return out
def _fixed_point(st: dict, f: np.ndarray, coef: dict, step) -> tuple[np.ndarray, np.ndarray]:
"""F = step(기둥 묶음, F) 고정점 — 수렴한 활동체의 기둥은 빼며 셈(값은 빼지 않을 때와 같음).
coef = 기둥별 상수 배열(group 포함) · 반환 (F, 수렴 여부).
"""
n = st["n"]
sub = coef
conv = np.zeros(n, bool)
for _ in range(MAX_ITER):
f_new = step(sub, f)
done = np.abs(f_new - f) < TOL
f = np.where(conv, f, f_new)
conv |= done
if conv.all():
break
live = ~conv[sub["group"]]
if live.mean() < 0.6: # 남은 기둥만 — 매번 줄이면 줄이는 값이 더 듦
sub = {k: v[live] for k, v in sub.items()}
return f, conv
def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
grp, n = st["group"], st["n"]
wq = st["w"] + st["q"]
@@ -252,19 +273,29 @@ def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
lift = np.zeros(n)
np.maximum.at(lift, grp, neg)
f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0) + lift
conv = np.zeros(n, bool)
for _ in range(MAX_ITER):
s_res, _ = _s_res(st, f)
s = s_res / np.maximum(f[grp], 1e-3)
nn = (wq - s * st["sin_a"]) / st["cos_psi"] # 연직 평형에서 N
res = _gsum(grp, st["arm_s"] * s_res, n)
drive = _gsum(grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + nn * st["m_n"], n)
f_new = _finish(res, drive)
done = np.abs(f_new - f) < TOL
f = np.where(conv, f, f_new)
conv |= done
if conv.all():
break
# S_res = num / mα · N = (W+Q − S sinα)/cosψ → 구동 = d1 − Σ S_res·k_n / F (기둥별 상수 미리)
num = st["c"] * st["a_h"] + np.maximum(wq - st["u_ah"], 0.0) * st["tan_phi"]
d1 = _gsum(
grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + wq * st["m_n"] / st["cos_psi"], n
)
coef = {
"group": grp,
"num": num,
"cos_psi": st["cos_psi"],
"k": st["sin_a"] * st["tan_phi"],
"arm_s": st["arm_s"],
"k_n": st["sin_a"] * st["m_n"] / st["cos_psi"],
}
def step(c: dict, f: np.ndarray) -> np.ndarray:
g = c["group"]
fg = np.maximum(f, 1e-3)
s_res = c["num"] / (c["cos_psi"] + c["k"] / fg[g])
res = _gsum(g, c["arm_s"] * s_res, n)
drive = d1 - _gsum(g, s_res * c["k_n"], n) / fg
return _finish(res, drive)
f, conv = _fixed_point(st, f, coef, step)
_, m_a = _s_res(st, f)
ma_min = _gmin(grp, m_a, n)
# 못 쓰는 면 — 수렴 안 함 · mα ≤ 0(바닥 수직력이 음) · F ≤ 0 → nan(후보에서 빠짐)
@@ -280,17 +311,23 @@ def _gmin(grp: np.ndarray, values: np.ndarray, n: int) -> np.ndarray:
def _janbu(st: dict, f0: np.ndarray) -> np.ndarray:
grp, n = st["group"], st["n"]
drive = _gsum(grp, (st["w"] + st["q"]) * st["tan_a"] + st["kh"] * st["w"], n)
wq = st["w"] + st["q"]
drive = _gsum(grp, wq * st["tan_a"] + st["kh"] * st["w"], n)
f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0)
conv = np.zeros(n, bool)
for _ in range(MAX_ITER):
s_res, _ = _s_res(st, f)
f_new = _finish(_gsum(grp, s_res / st["cos_a"], n), drive)
done = np.abs(f_new - f) < TOL
f = np.where(conv, f, f_new)
conv |= done
if conv.all():
break
num = st["c"] * st["a_h"] + np.maximum(wq - st["u_ah"], 0.0) * st["tan_phi"]
coef = {
"group": grp,
"num": num / st["cos_a"],
"cos_psi": st["cos_psi"],
"k": st["sin_a"] * st["tan_phi"],
}
def step(c: dict, f: np.ndarray) -> np.ndarray:
g = c["group"]
s_res = c["num"] / (c["cos_psi"] + c["k"] / np.maximum(f, 1e-3)[g])
return _finish(_gsum(g, s_res, n), drive)
f, conv = _fixed_point(st, f, coef, step)
f[~conv | (f <= 0)] = np.nan
return f
+5 -1
View File
@@ -34,6 +34,7 @@ ROCK_KINDS = ("ripping_rock", "blasting_rock")
#: 층 → 물성 id(E03 기본표) · 못 뚫는 층
LAYER_MATERIALS = ("fill", "weathered_soil", "weathered_rock", "soft_rock")
LAYER_HARD = (False, False, False, True)
LAYER_LABELS = ("fill", "soil", "ripping_rock", "blasting_rock") # E02 label_map 열쇠
def station_from_b06(cross: dict, design: dict) -> dict:
@@ -255,12 +256,15 @@ def build_grid(
def strata_from_grid(
grid: dict,
materials: list[dict] | None = None,
label_map: dict[str, str] | None = None,
*,
water_depth_m: float | None = None,
with_surcharge: bool = True,
) -> Strata:
"""격자 → 기둥 솔버 지층(층 차례 = 머리 주석) · 물성은 id 로 고름(E03 꼴 · 빠지면 기본표).
label_map — E02 층 이름 → 물성 id(soil · ripping_rock · blasting_rock · fill) · 빠지면 기본 짝.
water_depth_m — 지하수위면 = 합성 지표 − 깊이(0 = 지표까지 포화 · KDS 우기 규칙)
· None = 지하수 없음.
"""
@@ -276,7 +280,7 @@ def strata_from_grid(
grid["y"],
z,
np.stack([fill_bot, soil_bot, rip_bot]),
[by_id[k] for k in LAYER_MATERIALS],
[by_id[(label_map or {}).get(lab, k)] for lab, k in zip(LAYER_LABELS, LAYER_MATERIALS)],
hard=list(LAYER_HARD),
z_water=None if water_depth_m is None else z - water_depth_m,
q=grid["q"] if with_surcharge else None,
+233 -93
View File
@@ -9,6 +9,7 @@
구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1).
- 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈(ndimage.label · 8 이웃).
- 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음).
- 후보 덩이는 프로세스로 나눠 셈 · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값.
다스리는 안전율(governing):
- 못 뚫는 층(기본 발파암 = 연암)에 잘린 바닥 면적 비 ≥ 10 % → Janbu 보정으로 전환
@@ -21,9 +22,15 @@
from __future__ import annotations
import math
import os
import pickle
import tempfile
import time
import weakref
from collections.abc import Callable
from concurrent.futures import ProcessPoolExecutor
from dataclasses import dataclass, field
from types import SimpleNamespace
import numpy as np
from scipy import ndimage
@@ -57,6 +64,7 @@ class SearchParams:
a_max: float = math.inf
ratios: tuple[float, ...] = (1.0, 2.0) # 구(늘 봄) + 타원체 축 방향 배수(구 임계 노드에서)
sub: int = 2 # 셀 안 표본 sub×sub(부분 기둥)
workers: int = 0 # 프로세스 수 · 0 = 자동(격자 5만 칸 이상이면 CPU − 2 · 12 까지)
refine_rounds: int = 4 # 가장 고운 간격에서 더 내려가 보는 횟수 상한
# 고운 탐색 · 타원체는 F < max(이 값, 전역 최소 × 1.2) 로 고친 노드만(높은 곳 지도는 거친 값)
refine_below: float = 2.0
@@ -107,13 +115,23 @@ def trial_columns(
ixg = ix0[:, None, None] + off[None, None, :]
iyg = iy0[:, None, None] + off[None, :, None]
inb = (ixg >= 0) & (ixg < nx) & (iyg >= 0) & (iyg < ny)
ixc, iyc, inb = np.broadcast_arrays(np.clip(ixg, 0, nx - 1), np.clip(iyg, 0, ny - 1), inb)
shape = ixc.shape
nw = 2 * w + 1
shape = (n_t, nw, nw)
d = None if ratio == 1.0 else np.asarray(directions, float)
k2 = 1.0 / (ratio * ratio)
def remainder(t, xs, ys):
"""r² − (중심에서 수평 거리 꼴)² — 0 보다 크면 활동면이 그 점 위를 지남."""
dxp = xs - centers[t, 0]
dyp = ys - centers[t, 1]
if d is None:
return r * r - dxp * dxp - dyp * dyp
pp = dxp * d[t, 0] + dyp * d[t, 1]
qq = -dxp * d[t, 1] + dyp * d[t, 0]
return r * r - pp * pp - k2 * qq * qq
def surface(t, xs, ys):
"""후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작 · 창 모두)."""
"""후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작)."""
dxp = xs - centers[t, 0]
dyp = ys - centers[t, 1]
if d is None:
@@ -130,13 +148,16 @@ def trial_columns(
return rem, centers[t, 2] - root, gx, gy
# 1) 셀 중심으로 거르기 — 걸친 셀과 그 이웃만 셀 안 표본을 봄(창 대부분은 허공 · 땅속 밖)
t_grid = np.broadcast_to(np.arange(n_t)[:, None, None], shape)
xs0 = strata.x[0] + ixg * dx
ys0 = strata.y[0] + iyg * dx
rem0, zs0, _, _ = surface(t_grid, xs0, ys0)
# 창 칸 번호 → (후보, 셀) — 큰 창 배열을 펴지 않고 셈
def cell_of(flat):
t = flat // (nw * nw)
return t, iy0[t] + off[(flat // nw) % nw], ix0[t] + off[flat % nw]
rem0 = remainder(np.arange(n_t)[:, None, None], strata.x[0] + ixg * dx, strata.y[0] + iyg * dx)
in0 = np.flatnonzero(inb & (rem0 > 0))
yy0, xx0 = iyc.ravel()[in0], ixc.ravel()[in0]
zb0 = np.maximum(zs0.ravel()[in0], strata.hard_top[yy0, xx0])
t0, yy0, xx0 = cell_of(in0)
zs0 = centers[t0, 2] - np.sqrt(rem0.ravel()[in0])
zb0 = np.maximum(zs0, strata.hard_top[yy0, xx0])
act0 = np.zeros(shape, bool)
act0.ravel()[in0[zb0 < strata.z_top[yy0, xx0]]] = True
cand = act0 if sub == 1 else ndimage.binary_dilation(act0, structure=_STRUCT) & inb
@@ -146,36 +167,47 @@ def trial_columns(
empty = (Columns(e, e, e, np.zeros(0), np.zeros(0), np.zeros(0), np.zeros((0, 3))), e, big)
if not len(ci):
return empty
tt = t_grid.ravel()[ci]
cy = iyc.ravel()[ci]
cx = ixc.ravel()[ci]
tt, cy, cx = cell_of(ci)
hard = strata.hard_top[cy, cx]
# 2) 부분 기둥 — sub×sub 점
n_act = np.zeros(len(ci))
acc = {k: np.zeros(len(ci)) for k in ("depth", "zs", "gx", "gy", "px", "py")}
fr = ((np.arange(sub) + 0.5) / sub - 0.5) if sub > 1 else np.zeros(1)
for fy in fr:
for fx in fr:
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")}
for ka, fy in enumerate(fr):
for kb, fx in enumerate(fr):
xs = strata.x[0] + (cx + fx) * dx
ys = strata.y[0] + (cy + fy) * dx
rem, z_s, gx, gy = surface(tt, xs, ys)
zg = strata.ground_at(cy, cx, fy, fx)
z_b = np.maximum(z_s, hard)
act = (rem > 0) & (z_b < zg)
n_act += act
acc["depth"] += np.where(act, zg - z_b, 0.0)
for key, val in (("zs", z_s), ("gx", gx), ("gy", gy), ("px", xs), ("py", ys)):
acc[key] += np.where(act, val, 0.0)
keep = n_act > 0
ci, tt, cy, cx, n_act = ci[keep], tt[keep], cy[keep], cx[keep], n_act[keep]
acc = {k: v[keep] for k, v in acc.items()}
active = np.zeros(shape, bool)
active.ravel()[ci] = True
labels, n_lab = ndimage.label(active, structure=_STRUCT)
if n_lab == 0:
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)
for key, val in (
("depth", zg - z_b),
("zs", z_s),
("gx", gx),
("gy", gy),
("px", xs),
("py", ys),
):
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 = labels.ravel()[ci]
lab_pt, n_lab = _split_pieces(strata, pt, tt, cy, cx, hard, n_t * nf * nf, nf, sub, surface)
# 기둥 = (셀, 덩이) 짝 — 두 덩이가 한 셀을 나누면 기둥 둘
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 = 없음)
n_act = np.bincount(inv, minlength=len(key)).astype(float)
acc = {
k: np.bincount(inv, weights=pt[k], minlength=len(key))
for k in ("depth", "zs", "gx", "gy", "px", "py")
}
tt, cy, cx = tt[cell_k], cy[cell_k], cx[cell_k]
a_h = dx * dx
depth_cell = acc["depth"] / (sub * sub)
vol = np.bincount(lab, weights=depth_cell, minlength=n_lab + 1) * a_h
@@ -211,6 +243,87 @@ def trial_columns(
return cols, trial, big
_PART_OPTS = ("v_min", "v_max", "a_min", "a_max", "sub", "groundwater", "seismic_kh")
_PART: dict = {} # 일꾼 프로세스의 지층 · 조건
def _part_init(strata: Strata | str, opts: dict) -> None:
if isinstance(strata, str):
with open(strata, "rb") as fh:
strata = pickle.load(fh)
_PART.update(strata=strata, **opts)
def _open_pool(strata: Strata, opts: dict, n_workers: int) -> ProcessPoolExecutor:
"""일꾼 풀 — 지층은 임시 파일로 넘김(spawn 관으로 넘기면 일꾼이 차례로 뜸 · 12 개 6 s)."""
fd, path = tempfile.mkstemp(suffix=".e04.pkl")
with os.fdopen(fd, "wb") as fh:
pickle.dump(strata, fh, protocol=pickle.HIGHEST_PROTOCOL)
pool = ProcessPoolExecutor(n_workers, initializer=_part_init, initargs=(path, opts))
weakref.finalize(pool, os.remove, path) # 멈춤 · 오류로 빠져나가도 지움
return pool
def _part(job: tuple, k: dict | None = None) -> tuple:
"""후보 덩이 하나 — 기둥 + 솔버. 셀 합치기에 쓰는 칸(iy · ix · group)만 돌려줌(옮길 양 줄임)."""
xyz, dirs, r, ratio = job
k = k or _PART
st = k["strata"]
cols, trial, big = trial_columns(
st, xyz, r, ratio, dirs, k["v_min"], k["v_max"], k["a_min"], k["a_max"], k["sub"]
)
res = solve(st, cols, k["groundwater"], k["seismic_kh"]) if len(trial) else None
light = SimpleNamespace(
iy=cols.iy.astype(np.int32), ix=cols.ix.astype(np.int32), group=cols.group.astype(np.int32)
)
return light, trial, big, res
def _split_pieces(strata, pt, tt, cy, cx, hard, n_fine, nf, sub, surface):
"""활동면이 지표를 자른 점들을 덩이로 — 이웃 점(4 방향) 사이 가운데 점도 땅속일 때만 이음.
한 셀 안에서도 덩이가 갈림(발끝을 살짝 비켜 지나는 구 = 사면 몸통 · 발끝 앞 접시 따로).
반환 (점별 덩이 번호 0.., 덩이 수).
"""
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])
r_ = np.concatenate(rows)
c_ = np.concatenate(cols_)
graph = coo_matrix((np.ones(len(r_), np.int8), (r_, c_)), shape=(n_pt, n_pt))
n_lab, lab = connected_components(graph, directed=False)
return lab, n_lab
def governing(res: dict) -> tuple[np.ndarray, np.ndarray]:
"""다스리는 안전율 · 방법(0 Bishop · 1 Janbu 보정) — 머리 주석 규칙."""
jc = np.where(np.isfinite(res["janbu_corrected"]), res["janbu_corrected"], np.inf)
@@ -272,6 +385,8 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
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]))
@@ -287,6 +402,14 @@ 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
# 후보 덩이 = 프로세스로 나눠 셈(_part) · 결과는 낸 차례대로 합침 → 일꾼 수와 무관한 같은 값
n_workers = p.workers or (
min(12, max(1, (os.cpu_count() or 2) - 2)) if ny * nx >= 50_000 else 1
)
opts = {k: getattr(p, k) for k in _PART_OPTS}
pool = _open_pool(strata, opts, n_workers) if n_workers > 1 else None
local = {"strata": strata, **opts}
def run(ids: np.ndarray, ratios: tuple[float, ...] = (1.0,)) -> set[int]:
"""노드들의 모든 반지름 · 비 — 셀 최소를 고친 노드 번호를 돌려줌."""
improved: set[int] = set()
@@ -310,27 +433,20 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
if not len(live):
continue
w_now = (2 * math.ceil(r * max(ratio, 1) / strata.dx) + 3) ** 2
chunk = max(1, BATCH_CELLS // max(w_now, 1))
for s in range(0, len(live), chunk):
part = live[s : s + chunk]
d_part = dirs[part] if ratio != 1.0 else None
cols, trial, big = trial_columns(
strata,
xyz[part],
r,
ratio,
d_part,
p.v_min,
p.v_max,
p.a_min,
p.a_max,
p.sub,
)
# 메모리 = 일꾼 수 합쳐 BATCH_CELLS · 일꾼마다 한 덩이 이상(작은 고운 탐색도 나눔)
chunk = max(1, BATCH_CELLS // (w_now * n_workers))
chunk = min(chunk, -(-len(live) // n_workers))
parts = [live[s : s + chunk] for s in range(0, len(live), chunk)]
jobs = [(xyz[q], dirs[q] if ratio != 1.0 else None, r, ratio) for q in parts]
outs = pool.map(_part, jobs) if pool else (_part(j, local) for j in jobs)
# 고친 노드 판정 기준 = 이 반지름 시작 값(덩이 나눔과 무관)
ref = f_min.copy()
thr = max(p.refine_below, 1.2 * best["governing"])
for part, (cols, trial, big, res) in zip(parts, outs):
done[part[big > p.v_max]] = True
counters["trials"] += len(part)
if not len(trial):
continue
res = solve(strata, cols, p.groundwater, p.seismic_kh)
counters["evaluated"] += len(trial)
f, meth = governing(res)
improved |= _update(
@@ -344,7 +460,8 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
volume,
method,
node_map,
max(p.refine_below, 1.2 * best["governing"]),
thr,
ref,
)
k = (
int(np.nanargmin(np.where(np.isfinite(f), f, np.inf)))
@@ -368,50 +485,57 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
)
return improved
step = max(1, int(p.coarse))
grid = np.stack(
np.meshgrid(
np.arange(0, nhs, step), np.arange(0, nys, step), np.arange(0, nxs, step), indexing="ij"
),
-1,
).reshape(-1, 3)
ids = (grid[:, 0] * nys + grid[:, 1]) * nxs + grid[:, 2]
seen.update(ids.tolist())
seeds = run(ids)
rounds = passes = 0
while seeds and (step > 1 or rounds < p.refine_rounds):
if step == 1:
rounds += 1
step = max(1, step // 2)
passes += 1
stage[:] = [min(40 + 8 * (passes - 1), 77), 8, f"고운 탐색 {passes}"]
cand = set()
for nid in seeds:
ih, rem = divmod(nid, nys * nxs)
iy, ix = divmod(rem, nxs)
for dh in (-step, 0, step):
for dy in (-step, 0, step):
for dxs in (-step, 0, step):
h2, y2, x2 = ih + dh, iy + dy, ix + dxs
if 0 <= h2 < nhs and 0 <= y2 < nys and 0 <= x2 < nxs:
cid = (h2 * nys + y2) * nxs + x2
if cid not in seen:
cand.add(cid)
if not cand:
try:
step = max(1, int(p.coarse))
grid = np.stack(
np.meshgrid(
np.arange(0, nhs, step),
np.arange(0, nys, step),
np.arange(0, nxs, step),
indexing="ij",
),
-1,
).reshape(-1, 3)
ids = (grid[:, 0] * nys + grid[:, 1]) * nxs + grid[:, 2]
seen.update(ids.tolist())
seeds = run(ids)
rounds = passes = 0
while seeds and (step > 1 or rounds < p.refine_rounds):
if step == 1:
break
continue
seen.update(cand)
seeds = run(np.fromiter(cand, np.int64))
rounds += 1
step = max(1, step // 2)
passes += 1
stage[:] = [min(40 + 8 * (passes - 1), 77), 8, f"고운 탐색 {passes}"]
cand = set()
for nid in seeds:
ih, rem = divmod(nid, nys * nxs)
iy, ix = divmod(rem, nxs)
for dh in (-step, 0, step):
for dy in (-step, 0, step):
for dxs in (-step, 0, step):
h2, y2, x2 = ih + dh, iy + dy, ix + dxs
if 0 <= h2 < nhs and 0 <= y2 < nys and 0 <= x2 < nxs:
cid = (h2 * nys + y2) * nxs + x2
if cid not in seen:
cand.add(cid)
if not cand:
if step == 1:
break
continue
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):
stage[:] = [85, 14, "타원체"]
run(crit, ell)
# 타원체 — 구 탐색이 찾은 임계 노드(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):
stage[:] = [85, 14, "타원체"]
run(crit, ell)
finally:
if pool: # 끊기(InterruptedError) · 오류에도 일꾼을 닫음
pool.shutdown(cancel_futures=True)
f_out = np.where(np.isfinite(f_min), f_min, np.nan)
return SearchResult(
f_min=f_out,
@@ -426,19 +550,35 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
def _update(
f, meth, res, cols, trial, node_ids, f_min, volume, method, node_map, refine_below=np.inf
f,
meth,
res,
cols,
trial,
node_ids,
f_min,
volume,
method,
node_map,
refine_below=np.inf,
ref=None,
) -> set[int]:
"""셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김)."""
"""셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김).
돌려주는 노드 = ref(없으면 f_min) 보다 낮춘 후보의 노드 중 F < refine_below.
"""
fc = f[cols.group]
ok = np.isfinite(fc)
iy, ix, grp, fc = cols.iy[ok], cols.ix[ok], cols.group[ok], fc[ok]
order = np.argsort(-fc, kind="stable")
iy, ix, grp, fc = iy[order], ix[order], grp[order], fc[order]
ref = f_min if ref is None else ref
hit = (fc < ref[iy, ix] * (1 - 1e-6)) & (fc < refine_below)
seeds = set(np.unique(node_ids[trial[grp[hit]]]).tolist())
better = fc < f_min[iy, ix] * (1 - 1e-6)
iy, ix, grp, fc = iy[better], ix[better], grp[better], fc[better]
f_min[iy, ix] = fc
volume[iy, ix] = res["volume"][grp]
method[iy, ix] = meth[grp]
nid = node_ids[trial[grp]]
node_map[iy, ix] = nid
return set(np.unique(nid[fc < refine_below]).tolist())
node_map[iy, ix] = node_ids[trial[grp]]
return seeds
+12 -8
View File
@@ -28,14 +28,13 @@ import numpy as np
from scipy.spatial import cKDTree
from E01_StabRegister import E01_StabRegister_Source as source
from E02_StabGround.E02_StabGround_Model import load_defaults as load_ground_defaults
from E02_StabGround import E02_StabGround_Model as ground_model
from E03_StabSection.E03_StabSection_Types import DEFAULT_ROAD_SURCHARGE_KPA
from E04_Stab3D import E04_Stab3D_Engine_Grid as grid_engine
from E04_Stab3D import E04_Stab3D_Engine_Search as search_engine
ENGINE_VERSION = "e04-3" # 엔진 판 — 바뀌면 지문이 달라져 다시 셈
STAGE_DIR = "E04_Stab3D"
GROUND_FILE = Path("E02_StabGround") / "ground.json"
SECTION_LENGTH_M = 100.0 # 자동 구간 길이
MAX_SECTION_M = 500.0 # 한 번에 셀 구간 상한(시간 · 메모리)
METHOD_NAMES = ("bishop", "janbu_corrected")
@@ -69,10 +68,12 @@ def auto_sections(stations: list[dict[str, Any]]) -> list[dict[str, float]]:
return out
def read_ground(project_root: Path) -> dict[str, Any] | None:
"""E02 작업본(없으면 None = 기본값)."""
path = Path(project_root) / GROUND_FILE
return json.loads(path.read_text(encoding="utf-8")) if path.is_file() else None
def read_ground(project_root: Path) -> dict[str, Any]:
"""E02 작업본 → 없으면 초기값 → 없으면 마스터 기본값(E05 와 같은 차례)."""
for path in (ground_model.ground_path(project_root), ground_model.initial_path(project_root)):
if path.is_file():
return json.loads(path.read_text(encoding="utf-8"))
return ground_model.pick(ground_model.load_defaults())
def fingerprint(
@@ -171,7 +172,7 @@ def _floor(strata: Any, best: dict[str, Any]) -> list[list[float]]:
def reference_fs(ground: dict[str, Any] | None, condition: str) -> float:
"""기준안전율(E02 · 없으면 마스터 기본) — 깎기 · 쌓기 중 큰 값(보수)."""
factors = (ground or load_ground_defaults()).get("safety_factors") or {}
factors = (ground or ground_model.load_defaults()).get("safety_factors") or {}
row = factors.get(condition) or {}
return float(max(row.get("cut", 1.5), row.get("fill", 1.5)))
@@ -202,7 +203,10 @@ def run_engine(
)
wet = settings["condition"] == "wet"
strata = grid_engine.strata_from_grid(
grid, (ground or {}).get("materials"), water_depth_m=0.0 if wet else None
grid,
(ground or {}).get("materials"),
(ground or {}).get("label_map"),
water_depth_m=0.0 if wet else None,
)
mask = _node_mask(stations, grid, settings["start_m"], settings["end_m"], half)
ratios = (1.0,) if settings["shape"] == "sphere" else (1.0, float(settings["ellipsoid_ratio"]))
+14 -14
View File
@@ -3,7 +3,7 @@
* E04 안정성평가 3D 해석(요구 347-1 · 요구서 4-2 · 4-3 · PLAN 533-3) — B05~B09 틀(group_workflow).
* 좌 = 구간 목록 · 탐색 설정(기본 · 고급 접힘) · [해석] · 진행률 · [끊기]
* 가운데 = 평면 지도(셀별 최소 안전율 · 임계 체적) · [3D 보기](지표 + 노선면 + 임계 활동면)
* 우 = 구간 결과(3D 최소 · 2D 최소(같은 구간 E03 측점 fs_min 최소) · 3D/2D 비 · 체적 · 방법)
* 우 = 구간 결과(3D 최소 · 2D 최소(같은 구간 E03 측점 · 3D 와 같은 조건 안전율 최소) · 3D/2D 비 · 체적 · 방법)
* 계산 = 서버 단독(POST …/stab/analyze3d · 뒤에서) · 화면은 값 · 그림만. 탐색 설정 = sessionStorage.
* ========================================================================== */
@@ -21,6 +21,7 @@ import {
FS_STEPS,
VOLUME_STEPS,
legendRows,
fs2dOf,
min2d,
ratio3d2d,
stepColor,
@@ -73,18 +74,13 @@ async function api<T>(path: string, init?: RequestInit): Promise<T> {
return body;
}
/** E03 2D 측점 요약(GET …/stab/result2d · stations[].fs_min · 측점 하나에 깎기 · 쌓기 두 줄) — 없으면 빈 목록 */
async function fetch2d(
projectId: string,
): Promise<{ chainage_m: number; fs_min: number | null }[]> {
/** E03 2D 측점 줄(GET …/stab/result2d · fs_by_condition · min = 조건별) — 없으면 빈 목록 */
type Station2d = Parameters<typeof fs2dOf>[0] & { chainage_m: number };
async function fetch2d(projectId: string): Promise<Station2d[]> {
try {
const body = await api<{ stations?: Record<string, unknown>[] }>(
`/projects/${projectId}/stab/result2d`,
);
return (body.stations ?? []).map((s) => ({
chainage_m: Number(s.chainage_m),
fs_min: typeof s.fs_min === "number" ? s.fs_min : null,
}));
const body = await api<{ stations?: Station2d[] }>(`/projects/${projectId}/stab/result2d`);
return (body.stations ?? []).map((s) => ({ ...s, chainage_m: Number(s.chainage_m) }));
} catch {
return [];
}
@@ -227,7 +223,7 @@ export function renderE04Stab3D(root: HTMLElement): void {
// ---- 상태 ----
let overview: Overview | null = null;
let result: Result3d | null = null;
let stations2d: { chainage_m: number; fs_min: number | null }[] = [];
let stations2d: Station2d[] = [];
let running = false;
const remember = (): void =>
@@ -290,7 +286,11 @@ export function renderE04Stab3D(root: HTMLElement): void {
return;
}
const s = result.summary;
const f2d = min2d(stations2d, s.start_m, s.end_m);
const f2d = min2d(
stations2d.map((x) => ({ chainage_m: x.chainage_m, fs_min: fs2dOf(x, s.condition) })),
s.start_m,
s.end_m,
);
const ratio = ratio3d2d(s.f_min, f2d);
const station = overview?.stations.reduce<StationPoint | null>(
(best, x) =>
+12
View File
@@ -180,6 +180,18 @@ export function rasterRgba(
return out;
}
/** 2D 측점 한 줄의 3D 조건 안전율 — 지배 비탈 fs_by_condition · 없으면 조건별 최소 min · 없으면 null */
export function fs2dOf(
row: {
fs_by_condition?: Record<string, number | null>;
min?: Record<string, number | null>;
},
condition: string,
): number | null {
const v = row.fs_by_condition?.[condition] ?? row.min?.[condition];
return typeof v === "number" && Number.isFinite(v) ? v : null;
}
/** 구간 안 2D 측점 최소 안전율(없으면 null) */
export function min2d(
stations: readonly { chainage_m: number; fs_min: number | null }[],
+1 -1
View File
@@ -40,7 +40,7 @@ DEFAULT_SETTINGS: dict[str, Any] = {
"lateral_planar_deg": engine.LATERAL_PLANAR_DEG,
"lateral_toppling_deg": engine.LATERAL_TOPPLING_DEG,
"gamma_r": None, # None = E02 측점 암 물성 γt
"water_fill_wet": 0.5, # 우기 수위 = 비탈 높이 ½(KDS 11 70 05 표 4.3-3 암반)
"water_fill_wet": 0.5, # 우기 = 균열 깊이 · 활동면 높이 ½(KDS 11 70 05 표 4.3-3 암반)
"crack": {"mode": "none", "z": None},
"bolt": {"force": 0.0, "plunge": 15.0},
}
@@ -85,3 +85,15 @@ def test_wet_cone_pore_pressure_ratio():
groundwater={"mode": "ru", "ru": 0.4},
)
assert search(st, p).best["governing"] == pytest.approx(1.626, rel=0.03)
def test_workers_same_result():
"""프로세스 나눔 — 일꾼 수와 무관하게 셀별 최소 F · 전역 최소가 같음."""
st = _embankment(10 / np.tan(np.radians(30)), 19, 20, 8.64)
kw = dict(box=(-2, 10, -4, 4), h_min=4, h_max=20, spacing=2.0, coarse=2, max_depth=16)
kw.update(v_min=50, v_max=2e4, ratios=(1.0, 2.0))
one = search(st, SearchParams(workers=1, **kw))
two = search(st, SearchParams(workers=2, **kw))
np.testing.assert_array_equal(one.f_min, two.f_min)
assert one.best["governing"] == two.best["governing"]
assert one.evaluated == two.evaluated
+32
View File
@@ -182,3 +182,35 @@ def test_node_mask_section_and_route_ends(root):
assert mask[(np.abs(yy) <= 12.0) & (xx >= 1) & (xx <= 99)].all()
part = run._node_mask(STATIONS, grid, 20.0, 60.0, 12.0)
assert not part[(xx < 19.5) | (xx > 60.5)].any() and part[(yy == 0) & (xx == 40)].all()
def test_ground_order_and_label_map(tmp_path):
"""물성 = E02 작업본 → 초기값 → 마스터 기본값(성토 φ 28) · 층 물성 = label_map 짝."""
import json
from E04_Stab3D import E04_Stab3D_Engine_Grid as grid_engine
assert {m["id"]: m["phi"] for m in run.read_ground(tmp_path)["materials"]}["fill"] == 28.0
init = tmp_path / "initial_snapshot" / "E02_StabGround" / "ground.json"
init.parent.mkdir(parents=True)
init.write_text(json.dumps({"tag": "initial"}), encoding="utf-8")
assert run.read_ground(tmp_path) == {"tag": "initial"}
work = tmp_path / "E02_StabGround" / "ground.json"
work.parent.mkdir(parents=True)
work.write_text(json.dumps({"tag": "work"}), encoding="utf-8")
assert run.read_ground(tmp_path) == {"tag": "work"}
z = np.full((2, 2), 10.0)
grid = {
"x": np.arange(2.0),
"y": np.arange(2.0),
"z": z,
"z_orig": z,
"rock_top": z - 3,
"rock_kind": np.ones((2, 2), int),
"q": np.zeros((2, 2)),
}
ids = [m["id"] for m in grid_engine.strata_from_grid(grid).materials]
assert ids == ["fill", "weathered_soil", "weathered_rock", "soft_rock"]
mapped = grid_engine.strata_from_grid(grid, label_map={"soil": "colluvium", "fill": "fill"})
assert [m["id"] for m in mapped.materials][:2] == ["fill", "colluvium"]
+4
View File
@@ -56,6 +56,9 @@ writeFileSync(process.argv[3], JSON.stringify({
min2d: R.min2d([{ chainage_m: 0, fs_min: 1.1 }, { chainage_m: 20, fs_min: 1.4 },
{ chainage_m: 40, fs_min: null }, { chainage_m: 60, fs_min: 0.9 }], 10, 50),
ratio: [R.ratio3d2d(1.68, 1.4), R.ratio3d2d(1.5, null), R.ratio3d2d(1.5, 0)],
fs2d: [R.fs2dOf({ fs_by_condition: { wet: 1.2 }, min: { wet: 1.1 } }, "wet"),
R.fs2dOf({ fs_by_condition: { dry: 1.6 }, min: { wet: 1.1 } }, "wet"),
R.fs2dOf({ fs_by_condition: { wet: null }, min: {} }, "wet"), R.fs2dOf({}, "dry")],
terrain: mesh(M.terrainMesh(c.grid, c.z)),
sphere: mesh(M.slipMesh(c.grid, c.z, c.sphere)),
ellipsoid: mesh(M.slipMesh(c.grid, c.z, c.ellipsoid)),
@@ -120,6 +123,7 @@ def test_raster_min_pixels_ratio(ts):
assert ts["pixelEmpty"] == 0 # 빈 셀 = 투명
assert ts["min2d"] == 1.4 # 10~50 m 안 = 20 m(1.4) · null 뺌
assert ts["ratio"][0] == pytest.approx(1.2) and ts["ratio"][1:] == [None, None]
assert ts["fs2d"] == [1.2, 1.1, None, None] # 지배 비탈 → 조건별 최소 → 없음
def _py_slip(best, x, y):
+19
View File
@@ -48,6 +48,25 @@ def test_plane_hoek_tension_crack():
close(r["crack_b"], 15.3357, 1e-4) # 마루 뒤 거리(RocPlane 검증 #1)
def test_wet_default_kds_half_depth():
"""KDS 11 70 05 표 4.3-3 암반 우기 = 균열 · 활동면 높이 ½ 심도 포화 · 기준 1.2."""
import json
from pathlib import Path
fill = model.DEFAULT_SETTINGS["water_fill_wet"]
assert fill == 0.5
z = 10.0
r = eng.plane_fs(**SAU_MAU, water_fill=fill, crack={"mode": "depth", "z": z})
close(r["V"], 0.5 * 10 * (z / 2) ** 2, 1e-9) # 균열 안 수위 = z/2
dry = eng.plane_fs(**SAU_MAU) # 균열 없음 = 활동면 오름 60 · 가운데 수두 = (½ · 60)/2
close(eng.plane_fs(**SAU_MAU, water_fill=fill)["U"], 0.5 * 10 * 15 * dry["A"], 1e-9)
root = Path(__file__).resolve().parents[2]
ground = json.loads(
(root / "resources/master_template/stability/ground_defaults.json").read_text("utf-8")
)
assert ground["safety_factors"]["wet"]["cut"] == 1.2
def test_plane_critical_crack_search_matches_hoek():
"""건조 · 윗면 수평 — 수치 탐색 임계 자리 = Hoek-Bray 식 z · b."""
r = eng.plane_fs(
+6 -3
View File
@@ -36,7 +36,10 @@ export const ui_locales_e06 = {
E06_Set_LateralPlanar: ["평면 측방 한계 (°)", "Planar lateral limit (°)"],
E06_Set_LateralTopple: ["전도 측방 한계 (°)", "Toppling lateral limit (°)"],
E06_Set_Gamma: ["암 단위중량 (kN/m³ · 비우면 E02)", "Rock unit weight (kN/m³, blank = E02)"],
E06_Set_WaterWet: ["우기 수위 비 (0~1)", "Wet water ratio (0–1)"],
E06_Set_WaterWet: [
"우기 수위 비 (0~1 · 균열 · 활동면 높이 대비)",
"Wet water ratio (0–1, of crack / sliding-plane height)",
],
E06_Set_Crack: ["인장균열", "Tension crack"],
E06_Crack_None: ["없음", "None"],
E06_Crack_Critical: ["임계(최소 안전율 자리)", "Critical (minimum FS)"],
@@ -45,8 +48,8 @@ export const ui_locales_e06 = {
E06_Set_BoltForce: ["볼트 힘 (평면 kN/m · 쐐기 kN)", "Bolt force (plane kN/m · wedge kN)"],
E06_Set_BoltPlunge: ["볼트 경사 (수평 아래 °)", "Bolt plunge (below horizontal °)"],
E06_SettingsHint: [
"측방 한계 기본 20° (Dips 기본 · 지식DB 03 전도는 30°) · 우기 = 비탈 높이 ½ 수위(KDS 11 70 05) · 볼트 = 수동 지지력.",
"Lateral limit default 20° (Dips default; knowledge DB uses 30° for toppling) · wet = water at half height (KDS 11 70 05) · bolts = passive support.",
"측방 한계 기본 20° (Dips 기본 · 지식DB 03 전도는 30°) · 우기 = 인장균열 깊이 또는 활동면 높이의 ½ 까지 지하수 · 기준 안전율 1.2(KDS 11 70 05 표 4.3-3 암반) · 볼트 = 수동 지지력.",
"Lateral limit default 20° (Dips default; knowledge DB uses 30° for toppling) · wet = groundwater to ½ of tension-crack depth or sliding-plane height, required FS 1.2 (KDS 11 70 05 Table 4.3-3, rock) · bolts = passive support.",
],
E06_Analyze: ["해석", "Analyze"],
E06_Running: ["해석 중…", "Analyzing…"],