Files

207 lines
10 KiB
Python

#!/usr/bin/env python3
"""Task 7: axis + standard-section sweep skeleton of the main line (본선), 4+745..5+480.
Inputs (all confirmed by the user, 2026-10-21):
- centreline nodes data/alignment_plan_001_centre.csv (midpoint of side ticks)
- profile data/profile_gl_002_draft.csv (A = 지반고 GL, B = 계획고 = 노면)
- section table DOC/03_설계/설계_표준단면치수표_2026-10-21_v0.1.00.md
R1 6,982 mm (lining inner), lining 300 mm -> outer R 7,282, excavation R 7,402
crown A = 6.615 m above road (F = 0.000 = road reference)
- cross slope: 4+745~5+105 -2 %, 5+105~5+165 linear -2 -> +2 %, 5+165~5+480 +2 %
Model (skeleton, not as-built):
- Section = arch (concentric centre c = 6.615 - 6.982) + vertical walls down to the road.
Walls and invert are APPROXIMATIONS and are flagged in the metadata. No invert.
- Superelevation: z += (S/100) * x across the section (sign convention ASSUMED, flagged).
- Coordinates: plan XY in m (station spacing 100 m matches plan units), origin at 4+745.
Outputs: model/tunnel_skeleton.glb, model/tunnel_skeleton_meta.json, data/verify_skeleton.json
"""
import csv
import json
import math
import sys
import numpy as np
import trimesh
ROOT = sys.argv[1] # Code/3d-drawing
C_CSV = f"{ROOT}/data/alignment_plan_001_centre.csv"
P_CSV = f"{ROOT}/data/profile_gl_002_draft.csv"
R1, R_OUT, R_EX = 6.982, 7.282, 7.402 # m (mm / 1000)
CROWN_A = 6.615 # m above road
CENTRE_C = CROWN_A - R1 # arch centre height above road (m)
S0, S1 = 4745, 5480 # modelled range (fully determined)
DS = 1.0 # station step for rings (m); pattern segments are 2 m
ANG = np.radians(np.linspace(40.0, 140.0, 41)) # arch arc, right springline -> left
def read_centre():
rows = [r for r in csv.DictReader(open(C_CSV, encoding="utf-8"))]
return np.array([[int(r["station_m"]), float(r["centre_x"]), float(r["centre_y"])] for r in rows])
def read_profile():
gl, dg = {}, {}
for r in csv.DictReader(open(P_CSV, encoding="utf-8")):
sta = int(r["station_anchor"])
(gl if float(r["across_m"]) < -25 else dg)[sta] = float(r["value"])
return gl, dg
def slope_pct(s):
if s < 5105: return -2.0
if s < 5165: return -2.0 + 4.0 * (s - 5105) / 60.0
return 2.0
cn = read_centre()
gl_d, dg_d = read_profile()
gl_s = np.array(sorted(gl_d)); gl_v = np.array([gl_d[k] for k in gl_s])
dg_s = np.array(sorted(dg_d)); dg_v = np.array([dg_d[k] for k in dg_s])
cn_s = cn[:, 0]
def axis_xy(s):
x = np.interp(s, cn_s, cn[:, 1]); y = np.interp(s, cn_s, cn[:, 2])
i = min(max(np.searchsorted(cn_s, s) - 1, 0), len(cn_s) - 2)
dx, dy = cn[i + 1, 1] - cn[i, 1], cn[i + 1, 2] - cn[i, 2]
L = math.hypot(dx, dy)
return x, y, dx / L, dy / L
ORIGIN = np.array(axis_xy(S0)[:2])
stations = np.arange(S0, S1 + 0.001, DS)
def section(R, s, slope):
"""ring (x lateral m, z above road m) for radius R, including superelevation."""
xs = np.concatenate([[R * math.cos(ANG[0])], R * np.cos(ANG), [R * math.cos(ANG[-1])]])
zs = np.concatenate([[0.0], CENTRE_C + R * np.sin(ANG), [0.0]])
zs = zs + slope * xs
return xs, zs
def ring_world(R, s):
x, y, tx, ty = axis_xy(s)
nx, ny = -ty, tx # left normal (horizontal)
zr = np.interp(s, dg_s, dg_v) # road (계획고 = 노면)
xs, zs = section(R, s, slope_pct(s) / 100.0)
P = np.column_stack([x - ORIGIN[0] + xs * nx, y - ORIGIN[1] + xs * ny, zr + zs])
return P
PAL = {"A": (42, 157, 143, 255), "B-1": (231, 111, 81, 255), "B-2": (244, 162, 97, 255),
"C-1": (233, 196, 106, 255), "C-2": (138, 177, 125, 255), "D-1": (38, 70, 83, 255),
"D-2": (100, 130, 160, 255), "P-2": (190, 110, 140, 255), "P-3": (150, 90, 60, 255),
"RP-1": (80, 170, 200, 255), "WL-1": (200, 200, 120, 255), "WL-2": (120, 90, 170, 255),
"개착부": (170, 170, 170, 255)}
ZONES = [(int(float(r["sta_from"])), int(float(r["sta_to"])), r["lining_pattern"], "엑셀 지보패턴표(당진방향)")
for r in csv.DictReader(open(f"{ROOT}/data/pattern_stations_dangjin.csv", encoding="utf-8"))]
def pattern_of(sta):
for a, b, pat, _ in ZONES:
if a <= sta < b:
return pat
return ZONES[-1][2]
def make_band_pattern(R, name):
rings = [ring_world(R, s) for s in stations]
n = len(rings[0]); m = len(rings)
V = np.vstack(rings)
F, C = [], []
for i in range(m - 1):
col = PAL[pattern_of(stations[i])]
for j in range(n - 1):
a = i * n + j; b = a + 1; c = a + n; d = c + 1
F += [[a, b, c], [b, d, c]]
C += [col, col]
mesh = trimesh.Trimesh(vertices=V, faces=np.array(F), process=False)
mesh.visual.face_colors = np.array(C, dtype=np.uint8)
mesh.metadata["name"] = name
mesh.metadata["ring_points"] = n
return mesh
def make_ribbon(z_fn, half_w, color, name):
V, F = [], []
for k, s in enumerate(stations):
x, y, tx, ty = axis_xy(s)
nx, ny = -ty, tx
zr = z_fn(s)
V.append([x - ORIGIN[0] - half_w * nx, y - ORIGIN[1] - half_w * ny, zr])
V.append([x - ORIGIN[0] + half_w * nx, y - ORIGIN[1] + half_w * ny, zr])
V = np.array(V)
for k in range(len(stations) - 1):
a = 2 * k; b = a + 1; c = a + 2; d = a + 3
F += [[a, c, b], [b, c, d]]
mesh = trimesh.Trimesh(vertices=V, faces=np.array(F), process=False)
mesh.visual.face_colors = np.tile(np.array(color, dtype=np.uint8), (len(F), 1))
mesh.metadata["name"] = name
return mesh
def make_band(R, color, name):
"""uniform-colour band (excavation outline)."""
rings = [ring_world(R, s) for s in stations]
n = len(rings[0]); m = len(rings)
V = np.vstack(rings)
F = []
for i in range(m - 1):
for j in range(n - 1):
a = i * n + j; b = a + 1; c = a + n; d = c + 1
F += [[a, b, c], [b, d, c]]
mesh = trimesh.Trimesh(vertices=V, faces=np.array(F), process=False)
mesh.visual.face_colors = np.tile(np.array(color, dtype=np.uint8), (len(F), 1))
mesh.metadata["name"] = name
return mesh
parts = [
make_band_pattern(R_OUT, "lining_outer_R7282"),
make_band_pattern(R1, "lining_inner_R6982"),
make_band(R_EX, (150, 120, 90, 110), "excavation_R7402"),
make_ribbon(lambda s: float(np.interp(s, dg_s, dg_v)), 5.0, (70, 70, 75, 255), "road_surface_계획고"),
make_ribbon(lambda s: float(np.interp(s, gl_s, gl_v)), 20.0, (110, 150, 90, 200), "ground_GL_지반고"),
]
scene = trimesh.Scene()
for p in parts:
scene.add_geometry(p, geom_name=p.metadata["name"])
glb = f"{ROOT}/model/tunnel_skeleton.glb"
scene.export(glb)
# ---- verification (assertions recorded, not forced) ----
checks = {}
# radius check in section plane: distance from arch centre (lateral=0, z=centre_c)
xs, zs = section(R1, S0, 0.0)
rad = np.hypot(xs[1:-1], zs[1:-1] - CENTRE_C)
checks["inner_radius_max_dev_m"] = float(np.max(np.abs(rad - R1)))
checks["lining_thickness_crown_m"] = round(R_OUT - R1, 3)
crown_z = ring_world(R1, 5000)[:, 2].max() - np.interp(5000, dg_s, dg_v)
checks["crown_above_road_m (expect 6.615 at S=0 station 5000)"] = round(float(crown_z), 4)
checks["road_z_4745_m"] = round(float(np.interp(S0, dg_s, dg_v)), 3)
checks["road_z_5480_m"] = round(float(np.interp(S1, dg_s, dg_v)), 3)
L_chord = float(np.sum([math.dist(cn[i, 1:], cn[i + 1, 1:]) for i in range(len(cn) - 1)]))
checks["plan_chord_length_total_m (4+700..5+500)"] = round(L_chord, 3)
checks["station_range_modelled"] = f"{S0} .. {S1} (ring step {DS} m, {len(stations)} rings)"
checks["pattern_zones (지보패턴표 당진방향)"] = len(ZONES)
checks["pattern_coverage_all_stations"] = bool(all(pattern_of(x) in PAL for x in stations))
checks["pattern_counts"] = {k: sum(1 for z in ZONES if z[2] == k) for k in PAL}
checks["glb_bytes"] = int(__import__("os").path.getsize(glb))
reload = trimesh.load(glb)
checks["glb_reload_geometries"] = len(reload.geometry)
json.dump(checks, open(f"{ROOT}/data/verify_skeleton.json", "w", encoding="utf-8"), ensure_ascii=False, indent=2)
meta = {
"model": "tunnel_skeleton.glb", "scope": f"{S0}~{S1} (4+745~5+480)",
"units": "m", "origin_plan_xy": [round(ORIGIN[0], 3), round(ORIGIN[1], 3)],
"sources": {"centreline": "C0060302-001 CA-ALGN-STSM (side-tick midpoint)",
"road_grade_B": "C0060302-002 CR-TABL-TEX2 (계획고=노면)",
"GL_A": "C0060302-002 CR-TABL-TEX2 (지반고)",
"section": "C0060303-001 CS-TABL-TEX1/2, CS-CONC-MAJR, CS-EXCV"},
"approximations": ["walls vertical down to road; invert not modelled",
"lining band thickness at walls is radial-consistent but not as-built",
"centreline and grade linearly interpolated between nodes",
"superelevation sign convention assumed (z += S/100 * x, x = left of tangent)"],
"layers": [p.metadata["name"] for p in parts],
"ring_step_m": DS, "ring_points": int(parts[0].metadata["ring_points"]), "first_station": S0,
"pattern_status": "지보패턴: 엑셀 참고_지보패턴현황 당진방향 표(4+740~5+500, 42행, 빈틈 없음) 적용. 출처 data/pattern_stations_dangjin.csv (사용자 결정 2026-10-21)",
"patterns": [{"from": a, "to": b, "pattern": pat} for a, b, pat, _ in ZONES],
"open": ["4+700~4+745 and 5+480~5+500 excluded ([미확정])", "지보패턴 측점표 적용(당진). 도면 E-1~2, RP-2, WL-3~4 는 당진 표에 없어 [미확정]. 대산 표 WL-4(4+758.85~4+768.85)는 참고만 [확인 필요]",
"lining reference recorded (meta.lining_reference, 04-004~009); sheet-to-station mapping [미확정]",
"[확정 2026-10-21] R1 6,982 mm 유지 (사용자 결정). 라이닝 도면 6장의 R1 6.892 m 표기는 확인 대상 오기로 기록", "pattern-group to excavation-pattern mapping: user explanation conflicts with filename (WL-4) - re-confirm [확인 필요]"]
}
_ls = json.load(open(f"{ROOT}/data/lining_spec.json", encoding="utf-8"))
meta["lining_reference"] = {
"source": _ls["source"], "concrete": "fck 27 MPa", "rebar": "SD400 (fy 400 MPa)", "aggregate_mm": 25,
"summary": "C0060304-004~009 | fck 27 MPa | SD400 | 입출구부 50 m 취약구간 철근 상세 (측점 대응 미확정)",
"sheets": [{"sheet": r["file"][:14], "groups": [[g["group"], g["kind"]] for g in r["pattern_details"]]} for r in _ls["sheets"]]
}
json.dump(meta, open(f"{ROOT}/model/tunnel_skeleton_meta.json", "w", encoding="utf-8"), ensure_ascii=False, indent=2)
print(json.dumps(checks, ensure_ascii=False, indent=1))