#!/usr/bin/env python3 """Task 8: B-rep solids of the main line (본선) 4+745..5+480 -> STEP (Rhino/Autodesk) + GLB (web viewer). Section geometry: standard section C0060303-001 (drawing 03-001), revision rev2 approved 2026-10-21. - lining (라이닝): arch R6982 (inner) / R7282 (outer) about centre C; haunch arcs R3712/R3910 (inner), R4454/R4699 (outer) run below road level as drawn; inner walls and inner chord (road level) as drawn. Lining is symmetric (no superelevation). H/I boxes removed from lining (2D). - excavation band (굴착 여유대): R7402 / R4574 / R4819 minus lining outer, above road level. Symmetric. - road surface / 계획고 (노면, 포장 280 mm): cross slope S per station (superelevation). Only this part tilts. - base / invert (기층·인버트, ASSUMED: not in drawing): region between lining outer and inner chord minus pavement and side drains. Its top follows the tilted pavement bottom. - side drains (측구): box 300 x 836 mm minus pipe r150 (drawing circles). - H/I boxes (공동구, ASSUMED by user decision 2026-10-21). - 시설한계 50 mm: REFERENCE polylines only (not solids), as decided by the user. Coordinates: local section x = lateral along the LEFT normal of the alignment (sign convention flagged), z = height above road centre. Superelevation: z += (S/100)*x (sign flagged). Outputs: model/tunnel_solids.step, model/tunnel_solids.glb, model/tunnel_solids_meta.json, data/verify_solids.json Usage: python3 scripts/build_solids.py "$PWD" """ import csv, json, math, os, sys, tempfile import numpy as np from shapely.geometry import Polygon, box from build123d import Vector, Wire, Solid, Compound, Color, export_step, export_stl, import_step import trimesh ROOT = sys.argv[1] MODE = sys.argv[2] if len(sys.argv) > 2 else 'all' # parts_a | parts_b | assemble (4 GB RAM: run in stages) PARTS = f"{ROOT}/model/parts" RANGE = (int(sys.argv[3]), int(sys.argv[4])) if len(sys.argv) > 4 else (0, 10**6) # run index range for parts_a os.makedirs(PARTS, exist_ok=True) C_CSV = f"{ROOT}/data/alignment_plan_001_centre.csv" P_CSV = f"{ROOT}/data/profile_gl_002_draft.csv" PAT_CSV = f"{ROOT}/data/pattern_stations_dangjin.csv" S0, S1 = 4745, 5480 STEP_M = 10.0 # section spacing (m): ruled loft between sections (STEP size vs curvature) CHUNK_M = 100 PIPE_SEG = 24 CLEAR_OFFSET = 0.05 CLEAR_LINES = 6 CLEAR_RIB_W = 0.02 # ---------- drawing geometry (mm) -> local m ---------- CX, CY, Y0 = 14126.29, 15305.02, 15672.0 def P(x, y): return ((x - CX) / 1000.0, (y - Y0) / 1000.0) def arc(cx, cy, r, a0, a1, n=120): return [(cx + r * math.cos(math.radians(a0 + (a1 - a0) * i / n)), cy + r * math.sin(math.radians(a0 + (a1 - a0) * i / n))) for i in range(n + 1)] def loc(pts): return [P(x, y) for x, y in pts] WALL_L_INNER = [(8771, 15751), (8771, 15826), (8651, 16001), (8621, 16411), (8471, 16411), (8171, 16411), (8045, 16411), (8034, 16446)] WALL_R_INNER = [(19481, 15567), (19481, 15642), (19601, 15817), (19631, 16227), (19781, 16227), (20181, 16227), (20251, 16239)] INNER_MM = (WALL_L_INNER + arc(11621.3, 17406.9, 3712, 195.6, 140) + arc(CX, CY, 6982, 140, 40) + arc(16479.6, 17279.7, 3910, 40, -15.6) + WALL_R_INNER[::-1]) OUTER_MM = (arc(11959.9, 17122.8, 4454, 204.7, 140) + arc(CX, CY, 7282, 140, 40) + arc(16105.0, 16965.3, 4699, 40, -23.7)) EX_MM = (arc(11959.9, 17122.8, 4574, 204.0, 140) + arc(CX, CY, 7402, 140, 40) + arc(16105.0, 16965.3, 4819, 40, -23.0)) inner_poly = Polygon(loc(INNER_MM)).buffer(0) outer_poly = Polygon(loc(OUTER_MM)).buffer(0) ex_poly = Polygon(loc(EX_MM)).buffer(0) def _fill_slivers(g, tol=1e-3): if g.geom_type != 'Polygon': return g return Polygon(g.exterior, [h for h in g.interiors if Polygon(h).area > tol]) ring = _fill_slivers(outer_poly.difference(inner_poly)) # 37 mm wall/haunch join leaves 1e-5 m2 slivers CHORD_A, CHORD_B = P(8771, 15751), P(19481, 15567) _k = (CHORD_B[1] - CHORD_A[1]) / (CHORD_B[0] - CHORD_A[0]) def chord_z(x): return CHORD_A[1] + _k * (x - CHORD_A[0]) below_chord = Polygon([(-9, chord_z(-9)), (9, chord_z(9)), (9, -5), (-9, -5)]) above_chord = Polygon([(-9, chord_z(-9)), (9, chord_z(9)), (9, 9), (-9, 9)]) XL = P(9526, 15764)[0]; XR = P(18726, 15580)[0] T_FIX = -0.02 PAV_T = 0.28 def pav_corners(t): return [(XL, t * XL), (XR, t * XR), (XR, t * XR - PAV_T), (XL, t * XL - PAV_T)] pav_fix = Polygon(pav_corners(T_FIX)) def rect_local(x0_mm, x1_mm, y0_mm, y1_mm): (a, b), (c, d) = P(x0_mm, y0_mm), P(x1_mm, y1_mm) return box(min(a, c), min(b, d), max(a, c), max(b, d)) DRAIN_L = rect_local(9101, 9401, 14864, 15700) DRAIN_R = rect_local(18851, 19151, 14680, 15520) PIPE_L = (P(9251.3, 15064)[0], P(9251.3, 15064)[1], 0.150) PIPE_R = (P(19001.3, 14880)[0], P(19001.3, 14880)[1], 0.150) HBOX = rect_local(8171, 8471, 16051, 16411) IBOX = rect_local(19781, 20181, 15867, 16227) lin_full = _fill_slivers(ring.intersection(above_chord)) assert lin_full.geom_type == 'Polygon' and len(lin_full.interiors) == 0, f"lining topology: {lin_full.geom_type}" LINING_POLY = lin_full.simplify(0.01, preserve_topology=True) # 10 mm: removes 30 um edges (loft invalid) and cuts ~700 -> ~60 section vertices (STEP size) EXCAV_POLY = ex_poly.difference(outer_poly).intersection(box(-9, 0.0, 9, 9)).simplify(0.01, preserve_topology=True) assert EXCAV_POLY.geom_type == 'Polygon' and len(EXCAV_POLY.interiors) == 0, f"excavation topology: {EXCAV_POLY.geom_type}" base_fix = ring.intersection(below_chord).difference(pav_fix).difference(DRAIN_L).difference(DRAIN_R) assert base_fix.geom_type == 'Polygon' and len(base_fix.interiors) == 0, f"base topology: {base_fix.geom_type}" BASE_PTS_FIX = np.array(base_fix.exterior.coords)[:-1] def _find(pt): d = np.hypot(BASE_PTS_FIX[:, 0] - pt[0], BASE_PTS_FIX[:, 1] - pt[1]) i = int(np.argmin(d)); assert d[i] < 1e-6, f"pavement corner not on base boundary: {d[i]}" return i IDX_BL = _find((XL, T_FIX * XL - PAV_T)); IDX_BR = _find((XR, T_FIX * XR - PAV_T)) def base_pts(t): pts = BASE_PTS_FIX.copy() pts[IDX_BL, 1] = t * XL - PAV_T pts[IDX_BR, 1] = t * XR - PAV_T return pts LIN_PTS = np.array(LINING_POLY.exterior.coords)[:-1] EXC_PTS = np.array(EXCAV_POLY.exterior.coords)[:-1] def PIPE_PTS(cx, cz, r): return np.array([(cx + r * math.cos(2 * math.pi * k / PIPE_SEG), cz + r * math.sin(2 * math.pi * k / PIPE_SEG)) for k in range(PIPE_SEG)]) DRAIN_L_PTS = np.array(DRAIN_L.exterior.coords)[:-1] DRAIN_R_PTS = np.array(DRAIN_R.exterior.coords)[:-1] HBOX_PTS = np.array(HBOX.exterior.coords)[:-1] IBOX_PTS = np.array(IBOX.exterior.coords)[:-1] # ---------- alignment, profile ---------- cn = np.array([[int(r["station_m"]), float(r["centre_x"]), float(r["centre_y"])] for r in csv.DictReader(open(C_CSV, encoding="utf-8"))]) dg_d = {} for r in csv.DictReader(open(P_CSV, encoding="utf-8")): if float(r["across_m"]) >= -25: dg_d[int(r["station_anchor"])] = float(r["value"]) 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(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 ORX, ORY = axis(S0)[0], axis(S0)[1] 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 def world(s, x, z): ax, ay, tx, ty = axis(s) nx, ny = -ty, tx zr = float(np.interp(s, dg_s, dg_v)) return (ax - ORX + x * nx, ay - ORY + x * ny, zr + z) # ---------- B-rep helpers ---------- def loft(stations, pts2_fn): wires = [] for s in stations: pts3 = [world(s, float(x), float(z)) for x, z in pts2_fn(s)] wires.append(Wire.make_polygon([Vector(*p) for p in pts3], close=True)) sol = Solid.make_loft(wires, ruled=True) if not sol.is_valid: # end caps of lofts with near-degenerate edges; ShapeFix repairs them sol = sol.fix() return sol def stations_of(a, b, step=STEP_M): n = max(1, int(round((b - a) / step))) return [float(v) for v in np.linspace(a, b, n + 1)] def chunks(): out = []; a = S0 while a < S1: b = min(a + CHUNK_M, S1); out.append((a, b)); a = b return out PAL = {"A": (42, 157, 143), "B-1": (231, 111, 81), "B-2": (244, 162, 97), "C-1": (233, 196, 106), "C-2": (138, 177, 125), "D-1": (38, 70, 83), "D-2": (100, 130, 160), "P-2": (190, 110, 140), "P-3": (150, 90, 60), "RP-1": (80, 170, 200), "WL-1": (200, 200, 120), "WL-2": (120, 90, 170), "개착부": (170, 170, 170)} COL = {"base": (141, 110, 99), "pav": (96, 125, 139), "drain": (30, 136, 229), "hi": (253, 216, 53), "clear": (255, 111, 0)} runs = [] for r in csv.DictReader(open(PAT_CSV, encoding="utf-8")): a, b = int(float(r["sta_from"])), int(float(r["sta_to"])); pat = r["lining_pattern"] a2, b2 = max(a, S0), min(b, S1) if b2 > a2: runs.append((a2, b2, pat)) runs.sort() assert runs[0][0] == S0 and runs[-1][1] == S1, f"pattern coverage {runs[0][0]}..{runs[-1][1]}" for (a1, b1, _), (a2, _, _) in zip(runs, runs[1:]): assert b1 == a2, f"pattern gap/overlap at {b1} vs {a2}" solids, meta_c, clear_wires = [], [], [] def add(sol, cid, label, group, color, pattern=None, a=None, b=None): if not sol.is_valid: sol = sol.fix() # booleans on lofts can leave cap faces invalid assert sol.is_valid, f"invalid solid {cid}" sol.label = cid; sol.color = Color(*[c / 255 for c in color]) solids.append(sol) meta_c.append({"id": cid, "label": label, "group": group, "pattern": pattern, "from_sta": a, "to_sta": b, "volume_m3": round(sol.volume, 4), "color": list(color)}) for i, (a, b, pat) in enumerate(runs): if MODE != 'parts_a' or not (RANGE[0] <= i < RANGE[1]): continue print('lining run', i + 1, pat, a, b, flush=True) lin_run = loft(stations_of(a, b), lambda s: LIN_PTS) lin_run = lin_run - loft(stations_of(a, b), lambda s: HBOX_PTS) - loft(stations_of(a, b), lambda s: IBOX_PTS) add(lin_run, f"LIN_{i+1:02d}_{pat}", f"라이닝 {pat} {a/1000:.3f}~{b/1000:.3f}", "라이닝", PAL[pat], pat, a, b) for i, (a, b, pat) in enumerate(runs): if MODE != 'parts_a' or not (RANGE[0] <= i < RANGE[1]): continue add(loft(stations_of(a, b), lambda s: EXC_PTS), f"EXC_{i+1:02d}_{pat}", f"굴착 여유대 {pat} {a/1000:.3f}~{b/1000:.3f}", "굴착면", (150, 120, 90), pat, a, b) for i, (a, b, pat) in enumerate(runs): # 기타 구성요소도 지보패턴 구간(run) 경계에서 분할 if MODE != 'parts_b' or not (RANGE[0] <= i < RANGE[1]): continue sts = stations_of(a, b) rid = f"{i+1:02d}_{pat}" rng = f"{a/1000:.3f}~{b/1000:.3f}" add(loft(sts, lambda s: np.array(pav_corners(slope_pct(s) / 100.0))), f"PAV_{rid}", f"노면 계획고 포장 {pat} {rng}", "노면(계획고)", COL["pav"], pat, a, b) add(loft(sts, lambda s: base_pts(slope_pct(s) / 100.0)), f"BASE_{rid}", f"기층·인버트(가정) {pat} {rng}", "기층·인버트(가정)", COL["base"], pat, a, b) for side, dpts, pipe in (("L", DRAIN_L_PTS, PIPE_L), ("R", DRAIN_R_PTS, PIPE_R)): bx = loft(sts, lambda s, d=dpts: d) pp = loft(sts, lambda s, p=pipe: PIPE_PTS(p[0], p[1], p[2])) add(bx - pp, f"DRN{side}_{rid}", f"측구 {side} {pat} {rng}", "배수(측구)", COL["drain"], pat, a, b) add(loft(sts, lambda s: HBOX_PTS), f"HBOX_{rid}", f"공동구 H(가정) {pat} {rng}", "공동구(가정)", COL["hi"], pat, a, b) add(loft(sts, lambda s: IBOX_PTS), f"IBOX_{rid}", f"공동구 I(가정) {pat} {rng}", "공동구(가정)", COL["hi"], pat, a, b) # 시설한계 50 mm reference lines (wires, not solids) clr = inner_poly.buffer(CLEAR_OFFSET).exterior clr_pts = [(x, z) for x, z in list(clr.coords)[:-1] if z >= chord_z(x) - 1e-9 and abs(x) <= 5.4] idxs = np.round(np.linspace(0, len(clr_pts) - 1, CLEAR_LINES)).astype(int) clear_lines = [] sts_all = stations_of(S0, S1) for k, ii in enumerate(idxs): x0, z0 = clr_pts[ii] pts3 = [world(s, x0, z0) for s in sts_all] clear_lines.append(pts3) if MODE != 'parts_b' or RANGE[0] != 0: continue w = Wire.make_polygon([Vector(*p) for p in pts3], close=False) w.label = f"REF_CLEAR_{k+1}"; w.color = Color(1.0, 0.44, 0.0) clear_wires.append(w) meta_c.append({"id": f"REF_CLEAR_{k+1}", "label": f"시설한계 50 mm 참조선 {k+1}", "group": "시설한계(참조선)", "pattern": None, "from_sta": S0, "to_sta": S1, "volume_m3": 0, "color": list(COL["clear"]), "reference": True}) if MODE in ('parts_a', 'parts_b'): tag = f'lin_exc_{RANGE[0]:02d}' if MODE == 'parts_a' else f'oth_{RANGE[0]:02d}' export_step(Compound(children=solids + clear_wires, label=tag), f"{PARTS}/{tag}.step") json.dump(meta_c, open(f"{PARTS}/{tag}.json", 'w', encoding='utf-8'), ensure_ascii=False, indent=1) print('parts written', tag, len(solids), 'solids', len(clear_wires), 'ref wires') sys.exit(0) # ---- post-processing modes (each call < 180 s, 4 GB RAM) ---- import glob TAGS = sorted(os.path.basename(x)[:-5] for x in glob.glob(f"{PARTS}/*.json")) from build123d import Solid as _Solid if MODE == 'mesh': # python build_solids.py ROOT mesh tag = sys.argv[3] imp = import_step(f"{PARTS}/{tag}.step") items = [x for x in (imp.children or list(imp.solids())) if isinstance(x, _Solid)] pj = [m for m in json.load(open(f"{PARTS}/{tag}.json", encoding='utf-8')) if not m.get('reference')] assert len(items) == len(pj), f"{tag}: {len(items)} solids vs {len(pj)} meta" sc = trimesh.Scene(); tmp = tempfile.mkdtemp(prefix='stl_') for k, (sol, mm) in enumerate(zip(items, pj)): stl = os.path.join(tmp, f"s{k}.stl") export_stl(sol, stl, tolerance=0.002, angular_tolerance=0.2) m = trimesh.load(stl, force='mesh') m.visual.face_colors = np.tile(np.array(list(mm['color']) + [255], dtype=np.uint8), (len(m.faces), 1)) sc.add_geometry(m, geom_name=mm['id']) sc.export(f"{PARTS}/{tag}.glb") print('mesh', tag, len(sc.geometry), 'geometries'); sys.exit(0) if MODE == 'step': # python build_solids.py ROOT step allitems = [] for tag in TAGS: imp = import_step(f"{PARTS}/{tag}.step") allitems += list(imp.children or list(imp.solids())) os.makedirs(f"{ROOT}/model", exist_ok=True) export_step(Compound(children=allitems, label="본선 4+745~5+480 표준단면 모델"), f"{ROOT}/model/tunnel_solids.step") print('step merged', len(allitems), 'items', os.path.getsize(f"{ROOT}/model/tunnel_solids.step"), 'bytes'); sys.exit(0) if MODE == 'final': # python build_solids.py ROOT final meta_c = [] for tag in TAGS: meta_c += json.load(open(f"{PARTS}/{tag}.json", encoding='utf-8')) scene = trimesh.Scene() for tag in TAGS: g = trimesh.load(f"{PARTS}/{tag}.glb") for name, geom in g.geometry.items(): scene.add_geometry(geom, geom_name=name) glb_faces = {} for k, pts3 in enumerate(clear_lines): # 시설한계 reference ribbons (visual only; exact wires in STEP) P3 = np.array(pts3); V = []; F = [] for i in range(len(P3)): d = P3[min(i + 1, len(P3) - 1)] - P3[max(i - 1, 0)] d = d / (np.linalg.norm(d) + 1e-12) w = np.cross(d, np.array([0, 0, 1.0])); w = w / (np.linalg.norm(w) + 1e-12) * CLEAR_RIB_W V += [P3[i] - w, P3[i] + w] for i in range(len(P3) - 1): a0, b0, c0, d0 = 2 * i, 2 * i + 1, 2 * i + 2, 2 * i + 3 F += [[a0, c0, b0], [b0, c0, d0]] m = trimesh.Trimesh(vertices=np.array(V), faces=np.array(F), process=False) m.visual.face_colors = np.tile(np.array([255, 111, 0, 255], dtype=np.uint8), (len(m.faces), 1)) scene.add_geometry(m, geom_name=f"REF_CLEAR_{k+1}") glb_path = f"{ROOT}/model/tunnel_solids.glb" scene.export(glb_path) step_path = f"{ROOT}/model/tunnel_solids.step" solid_meta = [x for x in meta_c if not x.get('reference')] meta = { "model": ["tunnel_solids.step (STEP B-rep, Rhino/Autodesk)", "tunnel_solids.glb (web viewer, same B-rep tessellated)"], "scope": f"{S0}~{S1} (4+745~5+480)", "units": "m", "origin_plan_xy": [round(ORX, 3), round(ORY, 3)], "sections": {"lining": "rev2: 라이닝 R6982/R7282, 하단 호 도면 끝점까지", "excavation": "굴착 여유대 R7402/R4574/R4819 (노면 위)", "road": "계획고 노면 포장 280 mm, 측점별 편경사 S (노면만)", "base": "기층·인버트: 가정 (도면 미표현)", "drains": "측구 300x836 mm, 관 r150", "hi": "H/I 박스: 공동구 가정(사용자 결정)", "clearance": "시설한계 50 mm 참조선 (솔리드 아님)"}, "assumptions": ["기층·인버트 형상·두께는 가정 (도면 미표현, 사용자 결정 '가정하여 생성')", "H/I 박스 = 공동구 가정 (사용자 결정)", "편경사 부호 z = (S/100)·x, x = 접선 좌측 법선 (부호 확인 필요)", "라이닝·굴착 여유대는 대칭 (편경사 미적용)", "측점 간 축선은 선형 보간, 단면은 10 m 간격 Ruled loft (곡선부 근사)", "라이닝 H/I 박스 제외 부분은 3D 불리언 차감 (공동구 가정)"], "pattern_runs": [{"id": f"LIN_{i+1:02d}_{p}", "pattern": p, "from": a, "to": b} for i, (a, b, p) in enumerate(runs)], "components": meta_c, "open": ["E-1, E-2, RP-2, WL-3 패턴 배치 [미확정] (당진 표 밖)", "NP-5BE-1-WL-4 해석 [확인 필요]", "04 공동구·배수 도면 상세 미반영 (측구·공동구 치수는 03-001 단면 기준)", "편경사 부호(좌·우) 확인 필요"] } json.dump(meta, open(f"{ROOT}/model/tunnel_solids_meta.json", "w", encoding="utf-8"), ensure_ascii=False, indent=2) ver = {"solid_count": len(solid_meta), "runs": len(runs), "pattern_coverage": f"{runs[0][0]}..{runs[-1][1]} contiguous", "lining_runs": sum(1 for x in solid_meta if x['group'] == '라이닝'), "excavation_runs": sum(1 for x in solid_meta if x['group'] == '굴착면'), "volume_sum_m3": round(sum(x['volume_m3'] for x in solid_meta), 3), "step_bytes": os.path.getsize(step_path), "glb_bytes": os.path.getsize(glb_path), "glb_geometries": len(trimesh.load(glb_path).geometry)} def _overlap(sample): base = Polygon(base_pts(slope_pct(sample) / 100.0)); pav = Polygon(pav_corners(slope_pct(sample) / 100.0)) return {"lining_vs_excavation": round(LINING_POLY.intersection(EXCAV_POLY).area, 9), "base_vs_pavement": round(base.intersection(pav).area, 9), "lining_vs_base": round(LINING_POLY.intersection(base).area, 9), "pavement_vs_lining": round(pav.intersection(LINING_POLY).area, 9), "base_vs_drainL": round(base.intersection(DRAIN_L).area, 9), "base_vs_drainR": round(base.intersection(DRAIN_R).area, 9)} ver["section_overlap_m2_by_station"] = {str(st): _overlap(st) for st in [4745, 5105, 5165, 5480]} json.dump(ver, open(f"{ROOT}/data/verify_solids.json", "w", encoding="utf-8"), ensure_ascii=False, indent=2) print(json.dumps(ver, ensure_ascii=False, indent=1)); sys.exit(0)