Files
Aislo/E06_StabRock/E06_StabRock_Model.py
T

465 lines
18 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
"""E06 절리 정본 꼴 · 검사 · 측점 해석 묶음 (요구 345 · 351 · PLAN 536-1).
정본 = `E06_StabRock/joints.json` · 초기값 = `initial_snapshot/E06_StabRock/joints.json`
(절리 없음 + 기본 설정 · 처음 GET 때 한 번). 절리군 한 줄 = 구간(from_m ~ to_m) · 경사방향 · 경사 ·
φj · cj · 지속성(선택). 측점 = 복사본 암 측점(리핑암 · 발파암)의 깎기 쪽만.
사면 = 측점 frame 방향 + 깎기 쪽(경사방향 = 노면 쪽) · 경사 = 깎기(암) 경사비 · 높이 = E03 과 같은
`side_cases` 높이 · 윗면 = 비탈 끝 바깥 원지반 경사. 물성 · 지진 · 기준안전율 = E02 지반 조건.
"""
from __future__ import annotations
import copy
import hashlib
import math
from functools import lru_cache
from pathlib import Path
from typing import Any
import numpy as np
from common_util.common_util_station_key import chainage_key
from E01_StabRegister import E01_StabRegister_Source as src
from E02_StabGround import E02_StabGround_Model as ground_model
from E03_StabSection.E03_StabSection_Engine_Model import ground_from_cross, side_cases
from E06_StabRock import E06_StabRock_Engine as engine
from E06_StabRock.E06_StabRock_Stereonet import PROJECTIONS, stereonet
STAGE_DIR = "E06_StabRock"
JOINTS_NAME = "joints.json"
VERSION = 1
ENGINE_FILES = ("E06_StabRock_Engine.py", "E06_StabRock_Stereonet.py", "E06_StabRock_Model.py")
DEFAULT_GAMMA_R = 25.0 # 암 단위중량 — E02 암 물성이 없을 때(지식DB 07 Hoek & Bray 25~27)
MIN_HEIGHT_M = 1.0 # 이보다 낮은 깎기는 암반 블록 검토 안 함
CRACK_MODES = ("none", "critical", "depth")
CONDITIONS = ("dry", "wet", "seismic")
DEFAULT_SETTINGS: dict[str, Any] = {
"projection": "equal_angle",
"lateral_planar_deg": engine.LATERAL_PLANAR_DEG,
"lateral_toppling_deg": engine.LATERAL_TOPPLING_DEG,
"gamma_r": None, # None = E02 측점 암 물성 γt
"min_block_m3": 0.1, # 이보다 작은 블록은 판정에서 뺌(SWedge 최소 쐐기 크기 · 소규모 = 낙석 관리)
"water_fill_wet": 0.5, # 우기 = 균열 깊이 · 활동면 높이 ½(KDS 11 70 05 표 4.3-3 암반)
"crack": {"mode": "none", "z": None},
"bolt": {"force": 0.0, "plunge": 15.0},
}
# 칸 범위 — 화면 입력 · 서버 검사 같은 값
LIMITS: dict[str, list[float]] = {
"chainage_m": [0.0, 100000.0],
"dip_direction": [0.0, 360.0],
"dip": [0.0, 90.0],
"phi": [0.0, 60.0],
"c": [0.0, 1000.0],
"persistence_m": [0.0, 1000.0],
"lateral_deg": [5.0, 45.0],
"gamma_r": [15.0, 35.0],
"water_fill_wet": [0.0, 1.0],
"min_block_m3": [0.0, 1000.0],
"crack_z": [0.1, 100.0],
"bolt_force": [0.0, 5000.0],
"bolt_plunge": [-30.0, 60.0],
}
# 기본 절리면 값(새 줄) — 지식DB 07 표 7: cj 0~30 · φj 25~35(Barton 기본마찰각) 보수 쪽
NEW_SET = {"phi": 30.0, "c": 0.0}
def joints_path(project_root: Path) -> Path:
return Path(project_root) / STAGE_DIR / JOINTS_NAME
def initial_path(project_root: Path) -> Path:
return Path(project_root) / "initial_snapshot" / STAGE_DIR / JOINTS_NAME
def build_initial() -> dict[str, Any]:
return {"version": VERSION, "sets": [], "settings": copy.deepcopy(DEFAULT_SETTINGS)}
def upgrade(doc: dict[str, Any]) -> dict[str, Any]:
"""빠진 설정 칸만 기본값으로 채운 사본."""
out = copy.deepcopy(doc)
settings = out.setdefault("settings", {})
for key, value in DEFAULT_SETTINGS.items():
settings.setdefault(key, copy.deepcopy(value))
out.setdefault("sets", [])
return out
# ---- 검사 ----------------------------------------------------------------------
def _num(value: Any) -> bool:
return isinstance(value, int | float) and not isinstance(value, bool) and math.isfinite(value)
def _in(value: Any, key: str) -> bool:
lo, hi = LIMITS[key]
return _num(value) and lo <= value <= hi
def check_joints(doc: dict[str, Any]) -> list[dict[str, str]]:
"""틀린 칸 [{path, code}] — code = range · order · name · id_dup · enum · type."""
errors: list[dict[str, str]] = []
def bad(path: str, code: str = "range") -> None:
errors.append({"path": path, "code": code})
sets = doc.get("sets")
if not isinstance(sets, list):
return [{"path": "sets", "code": "type"}]
seen: set[str] = set()
for k, row in enumerate(sets):
base = f"sets.{k}"
if not isinstance(row, dict):
bad(base, "type")
continue
ident = str(row.get("id") or "")
if not ident or ident in seen:
bad(f"{base}.id", "id_dup")
seen.add(ident)
if not str(row.get("name") or "").strip():
bad(f"{base}.name", "name")
for key, limit in (
("from_m", "chainage_m"),
("to_m", "chainage_m"),
("dip_direction", "dip_direction"),
("dip", "dip"),
("phi", "phi"),
("c", "c"),
):
if not _in(row.get(key), limit):
bad(f"{base}.{key}")
if _num(row.get("from_m")) and _num(row.get("to_m")) and row["to_m"] < row["from_m"]:
bad(f"{base}.to_m", "order")
if row.get("persistence_m") is not None and not _in(row["persistence_m"], "persistence_m"):
bad(f"{base}.persistence_m")
settings = doc.get("settings")
if not isinstance(settings, dict):
return errors + [{"path": "settings", "code": "type"}]
if settings.get("projection") not in PROJECTIONS:
bad("settings.projection", "enum")
for key in ("lateral_planar_deg", "lateral_toppling_deg"):
if not _in(settings.get(key), "lateral_deg"):
bad(f"settings.{key}")
if settings.get("gamma_r") is not None and not _in(settings["gamma_r"], "gamma_r"):
bad("settings.gamma_r")
if not _in(settings.get("water_fill_wet"), "water_fill_wet"):
bad("settings.water_fill_wet")
if not _in(settings.get("min_block_m3"), "min_block_m3"):
bad("settings.min_block_m3")
crack = settings.get("crack") if isinstance(settings.get("crack"), dict) else {}
if crack.get("mode") not in CRACK_MODES:
bad("settings.crack.mode", "enum")
elif crack["mode"] == "depth" and not _in(crack.get("z"), "crack_z"):
bad("settings.crack.z")
bolt = settings.get("bolt") if isinstance(settings.get("bolt"), dict) else {}
if not _in(bolt.get("force"), "bolt_force"):
bad("settings.bolt.force")
if not _in(bolt.get("plunge"), "bolt_plunge"):
bad("settings.bolt.plunge")
return errors
def clean(doc: dict[str, Any]) -> dict[str, Any]:
"""정본 칸만 · 숫자는 float."""
sets = []
for row in doc["sets"]:
item = {
"id": str(row["id"]),
"name": str(row["name"]).strip(),
**{
key: float(row[key])
for key in ("from_m", "to_m", "dip_direction", "dip", "phi", "c")
},
"persistence_m": None
if row.get("persistence_m") is None
else float(row["persistence_m"]),
}
item["dip_direction"] %= 360.0
sets.append(item)
s = doc["settings"]
crack = s["crack"]
return {
"version": VERSION,
"sets": sets,
"settings": {
"projection": s["projection"],
"lateral_planar_deg": float(s["lateral_planar_deg"]),
"lateral_toppling_deg": float(s["lateral_toppling_deg"]),
"gamma_r": None if s.get("gamma_r") is None else float(s["gamma_r"]),
"water_fill_wet": float(s["water_fill_wet"]),
"min_block_m3": float(s["min_block_m3"]),
"crack": {
"mode": crack["mode"],
"z": float(crack["z"]) if crack["mode"] == "depth" else None,
},
"bolt": {"force": float(s["bolt"]["force"]), "plunge": float(s["bolt"]["plunge"])},
},
}
# ---- 측점 사면 ------------------------------------------------------------------
def rock_stations(source: Path) -> list[dict[str, Any]]:
"""암 측점(리핑암 · 발파암)."""
return [s for s in src.stations(source) if s.get("ground_type") in ground_model.ROCK_LABELS]
def _upper_dip(ground: list, start: float, sign: float, H: float, face_dip: float) -> float:
"""비탈 끝 바깥 max(H, 5 m) 원지반 평균 경사(오르막만 · 사면 경사 − 1° 아래)."""
gx = np.array([p[0] for p in ground])
gz = np.array([p[1] for p in ground])
end = start + sign * max(H, 5.0)
end = min(max(end, gx[0]), gx[-1])
run = abs(end - start)
if run < 0.5:
return 0.0
rise = float(np.interp(end, gx, gz) - np.interp(start, gx, gz))
return float(min(max(0.0, math.degrees(math.atan2(rise, run))), face_dip - 1.0))
def station_slopes(cross: dict, design: dict) -> list[dict[str, Any]]:
"""측점의 깎기 쪽마다 사면 — {side, dip_direction, dip, height_m, upper_dip}."""
ground = ground_from_cross(cross)
line = sorted(
(float(p["offset_m"]), float(p["elevation_m"])) for p in design.get("design_line") or []
)
ratio = design.get("cut_slope_ratio")
if len(ground) < 2 or len(line) < 2 or not ratio:
return []
face_dip = math.degrees(math.atan2(1.0, float(ratio)))
left = np.array([float(v) for v in cross["frame"]["left_xy"]])
out = []
for case in side_cases(design, ground, line):
if case["role"] != "cut" or case["height_m"] < MIN_HEIGHT_M:
continue
side = case["side"]
sign = 1.0 if side == "left" else -1.0
# 비탈 끝 = 탐색 끝에서 되짚음(side_cases: reach = 비탈 끝 + sign · max(5, 2H))
day = case["limits"][1 if sign > 0 else 0] - sign * max(5.0, 2.0 * case["height_m"])
face = -sign * left # 노면 쪽으로 기움
dip_dir = math.degrees(math.atan2(face[0], face[1])) % 360.0
out.append(
{
"side": side,
"dip_direction": round(dip_dir, 2),
"dip": round(face_dip, 2),
"height_m": round(case["height_m"], 3),
"upper_dip": round(_upper_dip(ground, day, sign, case["height_m"], face_dip), 2),
}
)
return out
# ---- 해석 ----------------------------------------------------------------------
@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 _ground_inputs(ground: dict[str, Any]) -> dict[str, Any]:
"""E02 에서 E06 이 쓰는 칸 — 지문에도 이것만."""
return {
k: ground.get(k)
for k in ("materials", "label_map", "stations", "zones", "safety_factors", "seismic")
}
def _gamma_for(resolved: dict | None, materials: dict, override: float | None) -> float:
if override:
return float(override)
mat = materials.get((resolved or {}).get("lower") or (resolved or {}).get("upper"))
return float(mat["gamma_t"]) if mat else DEFAULT_GAMMA_R
def _fmt(value: float | None) -> float | None:
if value is None or not math.isfinite(value):
return None
return round(value, 4) + 0.0 # -0.0 → 0.0(들림 · 볼트 없음)
def _cases_for(slope: dict, joints: list[dict], kin: dict, settings: dict, gamma: float, conds):
"""위험 모드만 안전율 — 조건 = [(이름, 수위 비, kh)]."""
upper = {"dip_direction": slope["dip_direction"], "dip": slope["upper_dip"]}
by_id = {j["id"]: j for j in joints}
bolt = settings["bolt"]
min_block = settings["min_block_m3"]
crack = settings["crack"] if settings["crack"]["mode"] != "none" else None
out = []
for jid in kin["planar"]:
j = by_id[jid]
for name, fill, kh in conds:
r = engine.plane_fs(
H=slope["height_m"],
slope_dip=slope["dip"],
plane_dip=j["dip"],
c=j["c"],
phi=j["phi"],
gamma=gamma,
upper_dip=slope["upper_dip"],
crack=crack,
water_fill=fill,
kh=kh,
bolt_t=bolt["force"],
bolt_plunge=bolt["plunge"],
)
short = (
r is not None and j.get("persistence_m") is not None and j["persistence_m"] < r["A"]
)
volume = r["W"] / gamma if r else None # 평면 = 폭 1 m 당 면적(m²) → m³
out.append(
{
"mode": "planar",
"ids": [jid],
"condition": name,
"F": _fmt(r["F"]) if r else None,
"formed": r is not None,
"persistence_short": short,
"volume": _fmt(volume),
"small": volume is not None and volume < min_block,
"W": _fmt(r["W"]) if r else None,
"A": _fmt(r["A"]) if r else None,
"z": _fmt(r["z"]) if r else None,
}
)
for item in kin["wedge"]:
a, b = (by_id[k] for k in item["ids"])
for name, fill, kh in conds:
r = engine.wedge_fs(
j1=a,
j2=b,
slope=slope,
upper=upper,
H=slope["height_m"],
gamma=gamma,
water={"fill": fill},
kh=kh,
bolt={
"force": bolt["force"],
"trend": (slope["dip_direction"] + 180.0) % 360.0,
"plunge": bolt["plunge"],
},
)
length = r["line_length"] if r else None
short = any(
j.get("persistence_m") is not None
and length is not None
and j["persistence_m"] < length
for j in (a, b)
)
out.append(
{
"mode": "wedge",
"ids": list(item["ids"]),
"condition": name,
"F": _fmt(r["F"]) if r else None,
"formed": r is not None,
"persistence_short": short,
"sliding": r["mode"] if r else None,
"W": _fmt(r["W"]) if r else None,
"volume": _fmt(r["volume"]) if r else None,
"small": r is not None and r["volume"] < min_block,
}
)
return out
def analyze(source: Path, joints_doc: dict, ground: dict) -> list[dict[str, Any]]:
"""암 측점 × 깎기 쪽 — 사면 · 절리 · Markland · 안전율 · 판정."""
settings = joints_doc["settings"]
sets = joints_doc["sets"]
materials = {m["id"]: m for m in ground.get("materials") or []}
resolved = {chainage_key(r["chainage_m"]): r for r in ground_model.resolve_stations(ground)}
seismic = ground.get("seismic") or {}
factors = ground.get("safety_factors") or {}
conds = [("dry", 0.0, 0.0), ("wet", float(settings["water_fill_wet"]), 0.0)]
if seismic.get("enabled") and seismic.get("kh"):
conds.append(("seismic", 0.0, float(seismic["kh"])))
out = []
for station in rock_stations(source):
ch = float(station["chainage_m"])
key = chainage_key(ch)
joints = [j for j in sets if j["from_m"] <= ch <= j["to_m"]]
cross = src.cross_section(source, ch)
gamma = _gamma_for(resolved.get(key), materials, settings.get("gamma_r"))
for slope in station_slopes(cross, cross["design"]):
kin = engine.kinematic(
slope, joints, settings["lateral_planar_deg"], settings["lateral_toppling_deg"]
)
cases = _cases_for(slope, joints, kin, settings, gamma, conds)
for case in cases:
need = (factors.get(case["condition"]) or {}).get("cut")
case["required"] = need
case["ok"] = (
None
if case["F"] is None or need is None or case["small"]
else case["F"] >= need
)
out.append(
{
"key": key,
"chainage_m": ch,
"label": station.get("label"),
"ground_type": station.get("ground_type"),
"slope": slope,
"gamma_r": gamma,
"joints": [
{k: j[k] for k in ("id", "name", "dip_direction", "dip", "phi", "c")}
for j in joints
],
"kinematic": kin,
"cases": cases,
**_verdict(kin, cases),
}
)
return out
def _verdict(kin: dict, cases: list[dict]) -> dict[str, Any]:
"""판정 — NG(기준 미달 하나라도) · check(전도만 · 블록 안 생김) · OK · none(위험 모드 없음 ·
블록이 모두 소규모). 소규모 블록(최소 블록 크기 미만)은 판정 · 최소 Fs 에서 뺌."""
judged = [c for c in cases if not c["small"]]
values = [c["F"] for c in judged if c["F"] is not None]
min_f = min(values) if values else None
if any(c["ok"] is False for c in judged):
verdict = "NG"
elif kin["toppling"] or (judged and not values):
verdict = "check"
elif judged:
verdict = "OK"
else:
verdict = "none"
small = {(c["mode"], tuple(c["ids"])) for c in cases if c["small"]}
return {
"small_blocks": len(small),
"counts": {
"planar": len(kin["planar"]),
"wedge": len(kin["wedge"]),
"toppling": len(kin["toppling"]),
},
"min_F": min_f,
"verdict": verdict,
}
def station_net(entry: dict, settings: dict) -> dict[str, Any]:
"""한 측점 쪽 스테레오넷(해석 때 절리 그대로) — 두 투영 모두(화면이 고름)."""
joints = entry["joints"]
phi = min((j["phi"] for j in joints), default=NEW_SET["phi"])
nets = {
kind: stereonet(
entry["slope"],
joints,
phi,
kind,
settings["lateral_planar_deg"],
settings["lateral_toppling_deg"],
)
for kind in PROJECTIONS
}
return {**nets, "phi": phi}