Files
Aislo/E04_Stab3D/E04_Stab3D_Run.py
T

390 lines
16 KiB
Python

"""E04 3D 해석 한 번 — 입력 모으기 · 입력 지문 · 엔진 부르기 · 파생 저장본 (PLAN 533-3).
계산 자리 = ③ 서버 단독. 엔진(탐색 · 기둥 솔버 · 합성 지형) 주인 = 533-1 · 2 — 여기는 부르기만.
저장본 = E 프로젝트 `E04_Stab3D/results/<구간 열쇠>.json` · 입력 지문(복사본 지문 · E02 지반 ·
요청 · 엔진 판)이 같으면 다시 셈하지 않음.
래스터 꼴(화면 `logic/E04_Stab3D_Raster.ts` 와 같은 열쇠):
grid = {x0, y0, dx, nx, ny} — 셀 (iy, ix) 중심 = (x0 + ix·dx, y0 + iy·dx) · 줄 우선 평탄 배열
z = 합성 지표(원지형 + 노선면) · f = 셀별 최소 안전율 · v = 그 임계 체적 m³ · m = 방법(0 Bishop
· 1 Janbu 보정) · 닿지 않은 셀 = null
best = 전체 최소 활동면 {x, y, z(중심), r, ratio, azimuth_deg(활동 방향), f, volume, method, ix, iy,
hovland · bishop · janbu_corrected(같은 면 방법별 F) · truncated_ratio · floor([셀, 못 뚫는 층 윗면])
· cells(활동체 덩이 셀 번호 — 한 구가 자른 여러 덩이 중 이 덩이)}
타원체 = 회전축(활동 방향에 수직) 반축만 ratio·r — 533-2 `trial_columns` 와 같은 꼴.
"""
from __future__ import annotations
import hashlib
import json
import math
import time
from collections.abc import Callable
from datetime import datetime
from pathlib import Path
from typing import Any
import numpy as np
from scipy.spatial import cKDTree
from E01_StabRegister import E01_StabRegister_Source as source
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-4" # 엔진 판 — 바뀌면 지문이 달라져 다시 셈
STAGE_DIR = "E04_Stab3D"
SECTION_LENGTH_M = 100.0 # 자동 구간 길이
MAX_SECTION_M = 500.0 # 한 번에 셀 구간 상한(시간 · 메모리)
METHOD_NAMES = ("bishop", "janbu_corrected")
Progress = Callable[[int, str], None]
StopFlag = Callable[[], bool]
def results_dir(project_root: Path) -> Path:
return Path(project_root) / STAGE_DIR / "results"
def section_key(start_m: float, end_m: float) -> str:
return f"{start_m:.2f}_{end_m:.2f}"
def auto_sections(stations: list[dict[str, Any]]) -> list[dict[str, float]]:
"""측점을 앞에서부터 SECTION_LENGTH_M 이상 묶음 — 경계 = 측점 · 꼬리가 짧으면 앞 구간에 붙임."""
chainages = sorted({float(s["chainage_m"]) for s in stations})
out: list[dict[str, float]] = []
start = chainages[0] if chainages else 0.0
for ch in chainages[1:]:
if ch - start >= SECTION_LENGTH_M - 1e-6:
out.append({"start_m": start, "end_m": ch})
start = ch
if chainages and start < chainages[-1]:
if out and chainages[-1] - start < SECTION_LENGTH_M / 2:
out[-1]["end_m"] = chainages[-1]
else:
out.append({"start_m": start, "end_m": chainages[-1]})
return out
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(
source_fingerprint: str, ground: dict[str, Any] | None, settings: dict[str, Any]
) -> str:
body = json.dumps(
{
"source": source_fingerprint,
"ground": ground,
"settings": settings,
"engine": ENGINE_VERSION,
},
ensure_ascii=False,
sort_keys=True,
)
return hashlib.sha1(body.encode("utf-8")).hexdigest()
def section_stations(stations: list[dict[str, Any]], start_m: float, end_m: float) -> list[dict]:
return [s for s in stations if start_m - 1e-6 <= float(s["chainage_m"]) <= end_m + 1e-6]
def bbox_of(
stations: list[dict[str, Any]], half_width_m: float
) -> tuple[float, float, float, float]:
xs = [float(s["center_x"]) for s in stations]
ys = [float(s["center_y"]) for s in stations]
return (
min(xs) - half_width_m,
max(xs) + half_width_m,
min(ys) - half_width_m,
max(ys) + half_width_m,
)
# ---------------------------------------------------------------------------
# 엔진 — 합성 DEM · 지층(533-1 Grid) → 구 · 타원체 탐색 + 기둥 솔버(533-2 Search · Column)
# ---------------------------------------------------------------------------
def _engine_station(root: Path, station: dict[str, Any]) -> dict[str, Any]:
"""복사본 측점 → 533-1 측점 한 줄(중심 · 좌표틀 · 지반선 · B06 설계선 · 암 깊이 · 노면)."""
cross = source.cross_section(root, station["chainage_m"])
row = {
"chainage_m": station["chainage_m"],
"center_x": station["center_x"],
"center_y": station["center_y"],
"frame": cross["frame"],
"samples": cross["samples"],
}
return grid_engine.station_from_b06(row, cross["design"])
def _node_mask(
stations: list[dict[str, Any]],
grid: dict[str, Any],
start_m: float,
end_m: float,
half_width_m: float,
) -> np.ndarray:
"""탐색 노드 자리 = 중심선에서 half_width 안 · 가장 가까운 중심선 점의 누가가 구간 안
· 가장 가까운 점이 중심선 양 끝점인 셀(시점 앞 · 종점 뒤)은 뺌.
중심선 = 측점 중심을 이은 선(B06 누가 그대로 · route_points 는 접선 꺾은선이라 누가가 다름)."""
line = np.array(
[[float(s["center_x"]), float(s["center_y"]), float(s["chainage_m"])] for s in stations]
)
# 0.5 m 로 촘촘히 — 꺾인 점 사이 셀도 바른 누가
dense = [line[:1]]
for a, b in zip(line[:-1], line[1:]):
n = max(1, int(math.ceil(math.hypot(*(b[:2] - a[:2])) / 0.5)))
t = np.linspace(0, 1, n + 1)[1:, None]
dense.append(a + (b - a) * t)
line = np.vstack(dense)
xx, yy = np.meshgrid(grid["x"], grid["y"])
dist, at = cKDTree(line[:, :2]).query(np.column_stack([xx.ravel(), yy.ravel()]))
chain = line[at, 2]
keep = (dist <= half_width_m) & (chain >= start_m - 1e-6) & (chain <= end_m + 1e-6)
keep &= (at > 0) & (at < len(line) - 1)
return keep.reshape(xx.shape)
def _floor(strata: Any, best: dict[str, Any]) -> list[list[float]]:
"""임계면이 못 뚫는 층에 잘린 셀 — [셀 번호, 층 윗면] (화면 활동면 메쉬가 그 높이로 자름)."""
xx, yy = np.meshgrid(strata.x, strata.y)
cx, cy, cz = best["center"]
d = np.asarray(best.get("direction") or [0.0, 1.0], float)
pp = (xx - cx) * d[0] + (yy - cy) * d[1]
qq = (-(xx - cx) * d[1] + (yy - cy) * d[0]) / max(float(best["ratio"]), 1e-9)
rest = float(best["radius"]) ** 2 - pp**2 - qq**2
slip = np.where(rest >= 0, cz - np.sqrt(np.maximum(rest, 0)), np.nan)
hard = strata.hard_top
hit = np.isfinite(slip) & (slip < strata.z_top) & (hard > slip) & (hard < strata.z_top)
if "cells" in best: # 활동체 덩이 셀만
hit &= np.isin(np.arange(hit.size), best["cells"]).reshape(hit.shape)
return [[int(i), round(float(hard.ravel()[i]), 2)] for i in np.flatnonzero(hit.ravel())]
def reference_fs(ground: dict[str, Any] | None, condition: str) -> float:
"""기준안전율(E02 · 없으면 마스터 기본) — 깎기 · 쌓기 중 큰 값(보수)."""
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)))
def run_engine(
root: Path,
stations: list[dict[str, Any]],
ground: dict[str, Any] | None,
settings: dict[str, Any],
on_progress: Progress,
should_stop: StopFlag,
) -> dict[str, Any]:
"""엔진 한 번 → {grid, z, f, v, m(넘파이 (ny, nx)), best, trials}.
stations = 구간 측점 + 앞뒤 이웃 하나씩(구간 끝 합성 지형이 이어지게).
"""
on_progress(1, "합성 지형")
half = float(settings["box_half_width_m"])
margin = half + float(settings["max_depth_m"]) + 10.0 # 활동체가 격자 끝에 닿지 않게
with np.load(source.terrain_grid_path(root)) as npz:
dem = {k: npz[k] for k in ("x", "y", "z", "valid_mask")}
grid = grid_engine.build_grid(
dem,
[_engine_station(root, s) for s in stations],
bbox=bbox_of(stations, margin),
road_q_kpa=float((ground or {}).get("road_surcharge_kpa", DEFAULT_ROAD_SURCHARGE_KPA)),
cell_m=float(settings["cell_m"]),
)
wet = settings["condition"] == "wet"
strata = grid_engine.strata_from_grid(
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"]))
params = search_engine.SearchParams(
node_mask=mask,
spacing=float(settings["spacing_m"]),
h_min=float(settings["h_min_m"]),
h_max=float(settings["h_max_m"]),
h_step=float(settings["h_step_m"]),
dr=float(settings["dr_m"]),
max_depth=float(settings["max_depth_m"]),
v_min=float(settings["v_min_m3"]),
v_max=float(settings["v_max_m3"]),
ratios=ratios,
groundwater={"mode": "line"} if wet else None, # 지하수위면 = strata.z_water
on_progress=lambda pct, msg: on_progress(2 + int(pct * 0.97), msg),
should_stop=should_stop,
)
res = search_engine.search(strata, params)
best = res.best
if not math.isfinite(float(best.get("governing", math.inf))):
raise ValueError("조건에 맞는 활동면이 없습니다 — 부피 · 높이 범위를 넓히세요.")
d = best.get("direction") or [0.0, 1.0]
cx, cy, cz = best["center"]
return {
"grid": {
"x0": float(strata.x[0]),
"y0": float(strata.y[0]),
"dx": float(strata.dx),
"nx": len(strata.x),
"ny": len(strata.y),
},
"z": strata.z_top,
"f": res.f_min,
"v": res.volume,
"m": np.where(res.method >= 0, res.method, np.nan),
"best": {
"x": float(cx),
"y": float(cy),
"z": float(cz),
"r": float(best["radius"]),
"ratio": float(best["ratio"]),
"azimuth_deg": math.degrees(math.atan2(float(d[0]), float(d[1]))) % 360.0,
"f": float(best["governing"]),
"volume": float(best["volume"]),
"method": int(best["method"]),
"hovland": float(best["hovland"]),
"bishop": float(best["bishop"]),
"janbu_corrected": float(best["janbu_corrected"]),
"truncated_ratio": float(best["truncated_ratio"]),
"floor": _floor(strata, best),
"cells": [int(i) for i in best["cells"]],
},
"trials": int(res.trials),
"reference_fs": reference_fs(ground, settings["condition"]),
}
# ---------------------------------------------------------------------------
# 저장본
# ---------------------------------------------------------------------------
def _flat(values: np.ndarray, digits: int) -> list[float | None]:
return [None if not math.isfinite(v) else round(float(v), digits) for v in values.ravel()]
def _nearest_chainage(stations: list[dict[str, Any]], x: float, y: float) -> float:
return min(
stations, key=lambda s: (float(s["center_x"]) - x) ** 2 + (float(s["center_y"]) - y) ** 2
)["chainage_m"]
def pack(
raw: dict[str, Any],
stations: list[dict[str, Any]],
settings: dict[str, Any],
finger: str,
elapsed_s: float,
) -> dict[str, Any]:
"""엔진 결과 → 저장본(JSON) — 요약은 best 에서."""
grid = raw["grid"]
best = dict(raw["best"])
best["ix"] = int(round((best["x"] - grid["x0"]) / grid["dx"]))
best["iy"] = int(round((best["y"] - grid["y0"]) / grid["dx"]))
f = np.asarray(raw["f"], dtype=float)
cells = np.asarray(best.get("cells") or [], int)
# 가까운 측점 = 활동체 평면 중심(없으면 구 중심) — 화면 최소 표시와 같은 자리
where = (
(
grid["x0"] + (cells % grid["nx"]).mean() * grid["dx"],
grid["y0"] + (cells // grid["nx"]).mean() * grid["dx"],
)
if len(cells)
else (best["x"], best["y"])
)
return {
"fingerprint": finger,
"engine": ENGINE_VERSION,
"settings": settings,
"computed_at": datetime.now().isoformat(timespec="seconds"),
"elapsed_s": round(elapsed_s, 2),
"grid": grid,
"z": _flat(np.asarray(raw["z"], dtype=float), 2),
"f": _flat(f, 3),
"v": _flat(np.asarray(raw["v"], dtype=float), 1),
"m": _flat(np.asarray(raw["m"], dtype=float), 0),
"best": best,
"summary": {
"start_m": settings["start_m"],
"end_m": settings["end_m"],
"f_min": round(float(best["f"]), 3),
"volume": round(float(best["volume"]), 1),
"method": METHOD_NAMES[int(best.get("method") or 0)],
"condition": settings["condition"],
"reference_fs": raw["reference_fs"],
"verdict": "ok" if best["f"] >= raw["reference_fs"] else "ng",
"by_method": {
k: round(float(best[k]), 3) for k in ("hovland", "bishop", "janbu_corrected")
},
"chainage_m": float(_nearest_chainage(stations, *where)),
"cells": int(np.isfinite(f).sum()),
"trials": int(raw.get("trials") or 0),
},
}
def save(project_root: Path, key: str, result: dict[str, Any]) -> None:
folder = results_dir(project_root)
folder.mkdir(parents=True, exist_ok=True)
tmp = folder / f"{key}.json.tmp"
tmp.write_text(json.dumps(result, ensure_ascii=False, separators=(",", ":")), encoding="utf-8")
tmp.replace(folder / f"{key}.json")
def load(project_root: Path, key: str) -> dict[str, Any] | None:
path = results_dir(project_root) / f"{key}.json"
return json.loads(path.read_text(encoding="utf-8")) if path.is_file() else None
def saved_summaries(project_root: Path) -> list[dict[str, Any]]:
"""저장본 요약(래스터 빼고) — 구간 목록 색 · 판정용."""
out = []
for path in sorted(results_dir(project_root).glob("*.json")):
doc = json.loads(path.read_text(encoding="utf-8"))
out.append(
{
"key": path.stem,
"fingerprint": doc["fingerprint"],
"settings": doc["settings"],
**doc["summary"],
}
)
return out
def analyze(
project_root: Path,
root: Path,
settings: dict[str, Any],
finger: str,
on_progress: Progress,
should_stop: StopFlag,
) -> dict[str, Any]:
"""한 구간 3D 해석 → 저장본 쓰고 요약 돌려줌."""
started = time.perf_counter()
stations = sorted(source.stations(root), key=lambda s: float(s["chainage_m"]))
picked = section_stations(stations, settings["start_m"], settings["end_m"])
if len(picked) < 2:
raise ValueError("구간 안 측점이 둘 미만입니다.")
first, last = stations.index(picked[0]), stations.index(picked[-1])
around = stations[max(first - 1, 0) : last + 2] # 앞뒤 이웃 — 구간 끝 합성 지형 이음
ground = read_ground(project_root)
raw = run_engine(root, around, ground, settings, on_progress, should_stop)
result = pack(raw, picked, settings, finger, time.perf_counter() - started)
save(project_root, section_key(settings["start_m"], settings["end_m"]), result)
return result["summary"]