diff --git a/E04_Stab3D/E04_Stab3D_Engine_Search.py b/E04_Stab3D/E04_Stab3D_Engine_Search.py index 99ae971fd..3e95f6a82 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Search.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Search.py @@ -167,30 +167,43 @@ def trial_columns( 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 @@ -262,6 +275,51 @@ 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 방향) 사이 가운데 점도 땅속일 때만 이음. + + 한 셀 안에서도 덩이가 갈림(발끝을 살짝 비켜 지나는 구 = 사면 몸통 · 발끝 앞 접시 따로). + 반환 (점별 덩이 번호 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)