Files
Aislo/E04_Stab3D/E04_Stab3D_Run.py
T

434 lines
18 KiB
Python

"""E04 3D 해석 한 번 — 입력 모으기 · 입력 지문 · 엔진 부르기 · 파생 저장본 (PLAN 533-3).
계산 자리 = ③ 서버 단독. 엔진(탐색 · 기둥 솔버 · 합성 지형) 주인 = 533-1 · 2 — 여기는 부르기만.
저장본 = E 프로젝트 `E04_Stab3D/results/<구간 열쇠>.json` · 열쇠 = `시작_끝_조건`
(예 `0.00_100.00_wet` · 옛 `시작_끝` 파일은 그 안 조건으로 읽음) · 입력 지문(복사본 지문 ·
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 functools import lru_cache
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_FILES = (
"E04_Stab3D_Engine_Column.py",
"E04_Stab3D_Engine_Grid.py",
"E04_Stab3D_Engine_Search.py",
"E04_Stab3D_Run.py",
)
CONDITIONS = ("dry", "wet")
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"
@lru_cache(maxsize=1)
def engine_hash() -> str:
"""엔진 판 = 엔진 코드 지문 — 코드가 바뀌면 지문이 달라져 다시 셈."""
digest = hashlib.sha1()
for name in ENGINE_FILES:
digest.update((Path(__file__).parent / name).read_bytes())
return digest.hexdigest()[:12]
def section_key(start_m: float, end_m: float, condition: str = "dry") -> str:
return f"{start_m:.2f}_{end_m:.2f}_{condition}"
def _doc_key(doc: dict[str, Any]) -> str:
summary = doc["summary"]
condition = summary.get("condition") or doc["settings"].get("condition") or "dry"
return section_key(float(summary["start_m"]), float(summary["end_m"]), condition)
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_hash(),
},
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"]],
"columns": int(best.get("columns") or len(best["cells"])),
"few_columns": bool(best.get("few_columns")),
},
"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_hash(),
"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()),
"columns": best.get("columns"),
"few_columns": bool(best.get("few_columns")),
"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:
"""열쇠 파일 → 없으면 옛 이름(`시작_끝`) 파일이 같은 조건일 때."""
folder = results_dir(project_root)
path = folder / f"{key}.json"
if path.is_file():
return json.loads(path.read_text(encoding="utf-8"))
base, _, condition = key.rpartition("_")
legacy = folder / f"{base}.json"
if condition in CONDITIONS and legacy.is_file():
doc = json.loads(legacy.read_text(encoding="utf-8"))
if _doc_key(doc) == key:
return doc
return None
def saved_summaries(project_root: Path) -> list[dict[str, Any]]:
"""저장본 요약(래스터 빼고) — 구간 목록 색 · 판정용.
열쇠 = 구간 + 조건 · 같은 열쇠면 새 이름 파일이 옛 이름 파일을 이김."""
found: dict[str, dict[str, Any]] = {}
for path in sorted(results_dir(project_root).glob("*.json")):
doc = json.loads(path.read_text(encoding="utf-8"))
key = _doc_key(doc)
if key in found and path.stem != key:
continue
found[key] = {
"key": key,
"fingerprint": doc["fingerprint"],
"settings": doc["settings"],
**doc["summary"],
}
return [found[k] for k in sorted(found)]
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"], settings["condition"]),
result,
)
return result["summary"]