""" 구조해석 서버(fea_server.py) — 포트 8097 ========================================= 하는 일: 제관웹(fab.html)에서 보낸 부재 메시 + 하중/구속 조건을 받아 CalculiX(ccx) 로 구조해석을 돌리고, 결과(변위·von Mises 응력·항복 대비 %)를 JSON 으로 돌려준다. [2단계 — 실제 형상 해석] 1단계의 '외곽 상자 근사' 를 버리고, 보내온 삼각형 메시를 그대로 씁니다. ① 정점 병합(weld) → ② 닫힌 솔리드인지 검사 → ③ tetgen 으로 사면체 채우기 → ④ 2차 사면체(C3D10)로 승격 → ⑤ ccx 해석 · 속 빈 각파이프·채널·구멍 뚫린 판재가 실제 단면 그대로 계산됩니다. ★ [2026-10-05] 부품 1개 실패가 전체를 죽이지 않는다 — '부품별 격리' 예전에는 사면체 메시를 만드는 전 과정이 try 하나로 묶여 있어서, 293부품 중 단 1개가 깨져도 전체가 '외곽 상자 근사'로 물러났다. 그러면 제품 전체가 속 찬 네모 덩어리가 되고 응력이 0에 가깝게 나온다 — 실패를 결과처럼 보여주는 가장 나쁜 동작이었다. 지금은 ① 부품 1개 = STL 1개 = tetgen 1회 (part_sizes 를 받으면 전부터 그랬다) ② 그 1회를 부품마다 따로 try 로 감싸, 깨진 부품만 빼고 나머지로 해석한다 ③ 빠진 부품은 번호·이름·사유를 결과 JSON(failed_parts)에 반드시 남긴다 ④ 자동 외곽상자 후퇴는 삭제했다. 쓸 수 있는 부품이 하나도 없으면 ok:false 와 "N번 부품(이름)에서 형상 처리 실패" 를 돌려준다. (force_obb=True 는 검증용 비교에만 남아 있다 — selftest 가 쓴다) ★ C3D4(1차 사면체)를 쓰지 않는 이유: 굽힘에서 전단잠김(shear locking)이 심해 얇은 벽 부재가 실제보다 2~5배 뻣뻣하게 나온다. 그래서 중간절점을 직접 만들어 C3D10(2차 사면체)로 올린다(tetgen -o2 의 절점 순서에 의존하지 않는다). 단위: mm / N / MPa 실행: python fea_server.py (서버 시작, 포트 8097) 검증: python fea_server.py --selftest (판재 + 각파이프 이론값 대조) 조사: python fea_server.py --survey (제품 OBJ 들이 닫힌 솔리드인지 전수 조사) """ import glob import json import math import os import pickle import shutil import subprocess import sys import tempfile import time import uuid import numpy as np # ── 외부 실행파일 ──────────────────────────────────────────────────── CCX_EXE = r'E:\도진팩토리\SYSTEMS\calculix\calculix_2.23_4win\ccx_static.exe' TETGEN_EXE = r'E:\도진팩토리\SYSTEMS\calculix\calculix_2.23_4win\tetgen.exe' # ★ ccx·tetgen 은 포트란/C 바이너리라 한글 경로에서 실패할 수 있다. # 그래서 작업폴더는 반드시 ASCII 전용 임시폴더를 쓴다. WORK_ROOT = os.path.join(os.environ.get('LOCALAPPDATA', tempfile.gettempdir()), 'Temp', 'dzw_fea') PORT = 8097 # ── 재질표 (현장에서 실제 쓰는 강종) ───────────────────────────────── # E=탄성계수(MPa), nu=포아송비, fy=항복강도(MPa), rho=밀도(ton/mm^3) MATERIALS = { 'SS275': {'name': 'SS275 (구 SS400 일반구조용강)', 'E': 205000., 'nu': 0.30, 'fy': 275., 'rho': 7.85e-9}, 'SS400': {'name': 'SS400 (구 규격, 두께 16 이하)', 'E': 205000., 'nu': 0.30, 'fy': 245., 'rho': 7.85e-9}, 'SM490': {'name': 'SM490 (용접구조용 고장력강)', 'E': 205000., 'nu': 0.30, 'fy': 325., 'rho': 7.85e-9}, 'S45C': {'name': 'S45C (기계구조용 탄소강)', 'E': 205000., 'nu': 0.30, 'fy': 490., 'rho': 7.85e-9}, 'STKR400': {'name': 'STKR400 (일반구조용 각파이프)', 'E': 205000., 'nu': 0.30, 'fy': 245., 'rho': 7.85e-9}, 'AL6061': {'name': 'AL6061-T6 (알루미늄)', 'E': 68900., 'nu': 0.33, 'fy': 275., 'rho': 2.70e-9}, } DEF_MAT = 'SS275' # 용접 허용전단응력 계수 — 모재 항복강도 대비. 필렛용접 목두께 기준 관용값(0.3fy 수준). WELD_TAU_FACTOR = 0.30 # 해석 규모 상한 (브라우저를 매달아 두지 않기 위한 사전 예산) # 실측 근거: C3D10 절점 394k → ccx 82.5초 / 164k → 약 25초 / 90k → 약 12초. # CEO 요구(30초 이내)를 지키려면 절점 상한을 이 부근에 둬야 한다. MAX_C3D10_NODES = 150000 # tetgen 은 품질조건(q1.4) 때문에 '부피/체적상한' 보다 훨씬 많은 사면체를 만든다. # 실측: 예측 29,000개 → 실제 147,770개(약 5배). 사전 예산 계산에 이 배수를 쓴다. TETGEN_OVERSHOOT = 5.0 # C3D10 절점수 ≈ 사면체수 × 1.9 (실측: 206,144개 → 393,965절점) C3D10_NODES_PER_TET = 1.9 # ── 전체해석(기계 1대 = 100부품 이상) 전용 상한 ────────────────────── # ★ 근거(실측 2026-10-09): 부품당 절점 바닥값이 12,000 이면 146부품 x 12,000 = 1.75M 까지 # 허용돼 '예산' 이 무력화된다(실제 절점 1,077,806개 발생). 바닥을 400 으로 내리고 # 부품당 상한 4,000 을 걸면 146부품이 절점 83,053개로 끝난다. MANY_PARTS = 40 # 이 개수를 넘으면 '전체해석 모드'(수리사다리 + -p 전용 + 부품당 상한) PER_NODE_FLOOR = 400 # 부품당 절점 바닥값 (기존 12,000) PER_CAP = 4000 # 부품당 절점 상한 — 전체해석 모드에서만 걸린다 COARSE_NOTE = '성긴 메시 — 국부 응력값은 참고용' # 부품당 tetgen 품질개선(-q) 제한시간(초). ★ 부품 수와 분리한다 — 근거는 계획 F-8: # 종전 공식(40/부품수)에서 14부품이면 2.86초가 되고, 그 시간을 못 넘긴 부품 8개가 # 조용히 메시 실패로 빠졌다(같은 14부품이 240초에서는 14/14 성공). QUALITY_TIMEOUT_S = 25.0 # 수리 결과 디스크 캐시. 키에 사다리 버전을 넣어 사다리를 고치면 옛 결과가 안 쓰이게 한다. # ★ 실패는 캐시하지 않는다(실측: 껍질 단계를 추가했는데 옛 '실패' 캐시가 돌아와 # 껍질이 한 번도 실행되지 않았다 — 같은 오류가 15.5초에 재현). REPAIR_CACHE = os.path.join(WORK_ROOT, 'repair_cache') LADDER_VER = 'L3-hull400' # ══════════════════════════════════════════════════════════════════ # 1. 삼각형 메시 손질 — 정점 병합 · 닫힘 검사 · 두께 추정 # ══════════════════════════════════════════════════════════════════ def weld_mesh(positions, tol=1e-4): """브라우저의 삼각형 수프(인덱스 없는 position)를 정점 병합해 (V, F) 로 만든다.""" P = np.asarray(positions, dtype=float).reshape(-1, 3) n = len(P) - (len(P) % 3) P = P[:n] key = np.round(P / tol).astype(np.int64) _, uniq, inv = np.unique(key, axis=0, return_index=True, return_inverse=True) V = P[uniq] F = inv.reshape(-1, 3) ok = (F[:, 0] != F[:, 1]) & (F[:, 1] != F[:, 2]) & (F[:, 0] != F[:, 2]) return V, F[ok] def mesh_metrics(V, F): """닫힌 솔리드인지 + 부피·표면적·평균 벽두께 추정. 벽두께 t ≈ 2V/S — 얇은 껍데기 형상에서 성립(각파이프 100x50x4.5 에서 4.5 근처).""" a, b, c = V[F[:, 0]], V[F[:, 1]], V[F[:, 2]] cr = np.cross(b - a, c - a) area = 0.5 * np.linalg.norm(cr, axis=1).sum() vol = abs(np.einsum('ij,ij->i', a, np.cross(b, c)).sum() / 6.0) # 발산정리 # 모든 엣지가 정확히 2번 쓰이면 닫힌 다양체 e = np.vstack([F[:, [0, 1]], F[:, [1, 2]], F[:, [2, 0]]]) e = np.sort(e, axis=1) _, cnt = np.unique(e, axis=0, return_counts=True) closed = bool((cnt == 2).all()) thick = (2.0 * vol / area) if area > 0 else 0.0 return {'closed': closed, 'volume_mm3': vol, 'area_mm2': area, 'wall_mm': thick, 'bad_edges': int((cnt != 2).sum())} def write_stl(path, V, F): with open(path, 'w') as f: f.write('solid s\n') for tri in F: p0, p1, p2 = V[tri[0]], V[tri[1]], V[tri[2]] nrm = np.cross(p1 - p0, p2 - p0) L = np.linalg.norm(nrm) nrm = nrm / L if L > 0 else np.array([0., 0., 1.]) f.write(' facet normal %g %g %g\n outer loop\n' % tuple(nrm)) for p in (p0, p1, p2): f.write(' vertex %.6f %.6f %.6f\n' % tuple(p)) f.write(' endloop\n endfacet\n') f.write('endsolid s\n') # ══════════════════════════════════════════════════════════════════ # 2. tetgen 으로 사면체 채우기 → C3D10 승격 # ══════════════════════════════════════════════════════════════════ def run_tetgen(V, F, wd, maxvol, switches=None, timeout=600): """반환: (절점 (N,3), 사면체 (M,4) 0-base) / 실패 시 예외""" stl = os.path.join(wd, 'shape.stl') write_stl(stl, V, F) for f in glob.glob(os.path.join(wd, 'shape.1.*')): os.remove(f) sw = switches if switches else ('-pq1.4a%.6f' % maxvol) pr = subprocess.run([TETGEN_EXE, sw, 'shape.stl'], cwd=wd, capture_output=True, text=True, errors='replace', timeout=timeout) nf = os.path.join(wd, 'shape.1.node') ef = os.path.join(wd, 'shape.1.ele') if not (os.path.isfile(nf) and os.path.isfile(ef)): raise RuntimeError('tetgen 실패: ' + (pr.stdout or '')[-500:].replace('\n', ' | ')) with open(nf, errors='replace') as f: head = f.readline().split() nn = int(head[0]) pts = np.zeros((nn, 3)) got = 0 for ln in f: if ln.startswith('#') or not ln.strip(): continue p = ln.split() pts[int(p[0]) - 1] = (float(p[1]), float(p[2]), float(p[3])) got += 1 if got >= nn: break with open(ef, errors='replace') as f: head = f.readline().split() ne = int(head[0]) tets = np.zeros((ne, 4), dtype=np.int64) got = 0 for ln in f: if ln.startswith('#') or not ln.strip(): continue p = ln.split() tets[int(p[0]) - 1] = [int(p[1]) - 1, int(p[2]) - 1, int(p[3]) - 1, int(p[4]) - 1] got += 1 if got >= ne: break return pts, tets def fix_orientation(pts, tets): """부피가 음수인 사면체는 절점 2개를 바꿔 양수로. (ccx 는 음의 야코비안을 거부)""" p0, p1, p2, p3 = pts[tets[:, 0]], pts[tets[:, 1]], pts[tets[:, 2]], pts[tets[:, 3]] v = np.einsum('ij,ij->i', p1 - p0, np.cross(p2 - p0, p3 - p0)) bad = v < 0 if bad.any(): tets[bad, 0], tets[bad, 1] = tets[bad, 1].copy(), tets[bad, 0].copy() return tets, int(bad.sum()) # CalculiX C3D10 절점순서: 1~4 꼭짓점, 5=중(1,2) 6=중(2,3) 7=중(3,1) 8=중(1,4) 9=중(2,4) 10=중(3,4) C3D10_EDGES = [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)] def make_c3d10(pts, tets): """1차 사면체 → 2차 사면체. 중간절점을 직접 만들어 ccx 순서를 보장한다 (tetgen -o2 의 내부 순서에 의존하지 않음).""" ea = np.concatenate([tets[:, [a, b]] for a, b in C3D10_EDGES], axis=0) ea = np.sort(ea, axis=1) uniq, inv = np.unique(ea, axis=0, return_inverse=True) mids = 0.5 * (pts[uniq[:, 0]] + pts[uniq[:, 1]]) nodes = np.vstack([pts, mids]) ne = len(tets) mid_idx = (inv.reshape(6, ne).T + len(pts)) return nodes, np.hstack([tets, mid_idx]) # ══════════════════════════════════════════════════════════════════ # 2-b. 수리 사다리 — '열린 형상' 부품을 사면체 메시가 되는 형상으로 고친다 # ══════════════════════════════════════════════════════════════════ # 왜 필요한가: CAD 에서 변환된 부재 146개 중 55개(38%)가 '닫힌 솔리드' 가 아니어서 # tetgen 이 거부한다. 그 55개 중 상당수가 '힘이 지나가는 교량' 이라 빼 버리면 # 사슬이 18곳에서 끊겨 전체해석이 원리적으로 불가능해진다(실측). # # 사다리 순서와 그 근거(전부 실측): # 원본 → 91/146 성공 # A manifold3d → +23. ★ 반드시 '원본' 에 걸어야 한다. 정리(tm_clean)를 먼저 하면 # 면 하나가 지워져 열려서 NotManifold 가 되고 수율이 49%→0% 가 된다. # B 중복·퇴화면 정리 → +2 # C 구멍메움(earcut) → +8 # D 메움+manifold3d # H 볼록껍질 대체 → +21 (교량 복구). 오목 부품이 속찬 덩어리가 되어 실제보다 # 단단해지므로 반드시 '껍질 근사' 로 표시해야 한다. # 정점병합 허용오차 상향 / tetgen -T 상향 은 0% 였다 — 쓰지 않는다. # 복셀 재메시는 쓰지 않는다 — 메모리 보호 루프가 pitch 를 키워 부품을 '큐브 1개' 로 # 바꿔 놓고도 tetgen 을 통과한다(조용히 틀린 형상). 연결 회복에도 무용했다. def _tet_ok_p(V, F, wd, timeout=30): try: pts, tets = run_tetgen(V, F, wd, 0.0, switches='-p', timeout=timeout) return len(tets) > 0 except Exception: return False def _thickness_ok(V, F, emin): """판두께 타당성 게이트. ★ tetgen 통과만으로 수리를 받아들이면 안 된다 — 구멍메움이 부피를 25배로 부풀려도 tetgen 은 통과시킨다(실측 Group73/81: 1,580mm3 → 40,690mm3). 판두께 t=2V/S 가 0.8mm 이상이고 부품 최소변의 1.2배+1 이하여야 물리적으로 가능하다.""" try: t = mesh_metrics(V, F)['wall_mm'] except Exception: return False return (t >= 0.8) and (t <= max(emin, 1.0) * 1.2 + 1.0) def _tm_clean(V, F): import trimesh m = trimesh.Trimesh(vertices=V, faces=F, process=False) m.merge_vertices(merge_tex=True, merge_norm=True) m.update_faces(m.nondegenerate_faces()) m.update_faces(m.unique_faces()) m.remove_unreferenced_vertices() try: m.fix_normals() except Exception: pass return np.asarray(m.vertices), np.asarray(m.faces) def _manifold_pass(V, F): import manifold3d as M mesh = M.Mesh(vert_properties=np.asarray(V, np.float32), tri_verts=np.asarray(F, np.uint32)) o = M.Manifold(mesh).to_mesh() return np.asarray(o.vert_properties[:, :3], float), np.asarray(o.tri_verts, np.int64) def _dir_loops(F): """유향 halfedge 로 경계 고리를 올바른 순서로 추출한다.""" he, cnt = set(), {} for t in F: for a, b in ((t[0], t[1]), (t[1], t[2]), (t[2], t[0])): k = (min(a, b), max(a, b)) cnt[k] = cnt.get(k, 0) + 1 he.add((int(a), int(b))) bnd = [(a, b) for (a, b) in he if cnt[(min(a, b), max(a, b))] == 1 and (b, a) not in he] nxt = {} for a, b in bnd: nxt.setdefault(a, []).append(b) used, loops = set(), [] for a0, _ in bnd: if a0 in used: continue loop, cur = [], a0 while cur in nxt and cur not in used: used.add(cur) loop.append(cur) cand = [x for x in nxt[cur] if x not in used] if not cand: break cur = cand[0] if len(loop) >= 3: loops.append(loop) return loops def _earcut_fill(V, F): """경계 고리를 최적합 평면에 투영해 earcut 으로 삼각분할해 덮는다.""" import mapbox_earcut as ec add = [] for lp in _dir_loops(F): idx = np.array(lp, dtype=np.int64) Q = V[idx] - V[idx].mean(axis=0) try: _, _, Vt = np.linalg.svd(Q, full_matrices=False) except Exception: continue uv = np.stack([Q @ Vt[0], Q @ Vt[1]], axis=1).astype(np.float64) A = 0.5 * np.sum(uv[:, 0] * np.roll(uv[:, 1], -1) - np.roll(uv[:, 0], -1) * uv[:, 1]) if A < 0: uv, idx = uv[::-1], idx[::-1] try: tri = ec.triangulate_float64(uv, np.array([len(uv)])) except Exception: continue if len(tri): add.append(idx[np.asarray(tri, np.int64).reshape(-1, 3)]) if not add: return V, F, 0 return V, np.vstack([F] + add), sum(len(a) for a in add) def _decimate(V, F, target): """삼각형 수를 목표치로 줄인다(open3d 2차오차 감쇄).""" import open3d as o3d m = o3d.geometry.TriangleMesh(o3d.utility.Vector3dVector(np.asarray(V, float)), o3d.utility.Vector3iVector(np.asarray(F, np.int32))) m = m.simplify_quadric_decimation(int(target)) m.remove_degenerate_triangles() m.remove_duplicated_vertices() m.remove_unreferenced_vertices() return np.asarray(m.vertices), np.asarray(m.triangles) HULL_TRI_CAP = 400 # 껍질 삼각형 상한 (실측: p158 껍질 3,260삼각형 → 11,939절점) # ── 두께가 0 인 '면 형상' 을 판두께 추정으로 압출한다 ──────────────────────── # 왜 필요한가 (실측 2026-10-11, 146부재 전체해석): 해석에서 빠진 17개 중 6개가 # '형상 처리 실패' 였고 그 중 5개가 `두께가 0인 면 형상` 이었다 # (Group44 300x0x300 / WINCH_BASE_PLATE_25T x4 28x30x0). # CAD 에서 판재가 '면' 으로만 내려온 것이라 속을 채울 수 없어 사면체가 안 된다. # 수리 사다리(manifold3d / 구멍메움 / 볼록껍질)는 전부 무용하다 — 부피가 0 인 # 형상은 무엇을 해도 부피가 0 이고, 볼록껍질도 납작한 삼각형 판이 된다. # → 판두께를 **추정해 압출**하는 것만이 길이다. 추정값이므로 반드시 라벨을 남긴다. # # ★ 'bbox 최소변' 은 이 경우에 쓸 수 없다 — 바로 그 변이 0 이기 때문이다. # 1순위: 경계상자가 맞닿는 이웃 부재들의 실측 벽두께 중앙값 # 2순위: 이 기계 전체 부재의 실측 벽두께 중앙값 # (부품 이름의 '25T' 같은 문자열은 쓰지 않는다 — 이름은 근거가 아니다) SHEET_FLAT_TOL_MM = 0.05 # bbox 최소변이 이것보다 작으면 '두께 0 면 형상' SHEET_T_MIN_MM = 1.0 # 추정 두께 하한 SHEET_T_MAX_MM = 40.0 # 추정 두께 상한 SHEET_MIN_SPAN_MM = 5.0 # 판으로 볼 최소 '면 안쪽' 크기 def is_flat_sheet(V): """bbox 최소변이 0 에 가까운가 → 두께 0 면 형상. ★ 단, 면 안쪽(나머지 두 변)도 판이라 할 만큼은 커야 한다. 실측: 85번 Group74 는 0 x 0.06 x 2.25mm 삼각형 1개다 — 두 변이 동시에 찌그러진 CAD 찌꺼기이고 부재가 아니다. 이것을 압출해 '살렸다' 고 적으면 계산 안 한 것(있지도 않은 부재)을 표시하는 짓이다. 사유를 그대로 남긴다.""" V = np.asarray(V, float) if len(V) < 3: return False, -1 ext = V.max(axis=0) - V.min(axis=0) ax = int(np.argmin(ext)) if ext[ax] > SHEET_FLAT_TOL_MM: return False, ax inplane = sorted(float(ext[k]) for k in range(3) if k != ax) if inplane[0] < SHEET_MIN_SPAN_MM: return False, ax return True, ax def guess_sheet_t(bb, neigh): """판두께 추정. neigh = [(이웃 bbox(lo,hi), 실측 벽두께), ...] 반환: (두께 mm, 근거 문장)""" lo, hi = bb near = [] for (nlo, nhi), w in neigh: if not (SHEET_T_MIN_MM <= w <= SHEET_T_MAX_MM): continue if np.all(np.asarray(nlo) - 30.0 <= np.asarray(hi)) and \ np.all(np.asarray(lo) - 30.0 <= np.asarray(nhi)): near.append(float(w)) if near: return float(np.median(near)), '이웃 부재 %d개 벽두께 중앙값' % len(near) allw = [float(w) for (_b, w) in neigh if SHEET_T_MIN_MM <= w <= SHEET_T_MAX_MM] if allw: return float(np.median(allw)), '이웃 없음 — 기계 전체 부재 벽두께 중앙값 %d개' % len(allw) return None, '추정 근거가 없습니다(벽두께를 잰 부재가 하나도 없음)' def extrude_flat_sheet(V, F, t, ax): """평면 삼각형 판을 두께 t 의 닫힌 솔리드로(판 중앙면 기준 +-t/2). 위/아래 면 + 경계엣지 측벽 → 모든 엣지가 정확히 2번 쓰이는 닫힌 껍질.""" V = np.asarray(V, float) F = np.asarray(F, np.int64) if len(F) < 1 or not (t > 0): return None d = np.zeros(3) d[int(ax)] = float(t) / 2.0 n = len(V) V2 = np.vstack([V - d, V + d]) cnt = {} for tr in F: for a, b in ((tr[0], tr[1]), (tr[1], tr[2]), (tr[2], tr[0])): k = (min(int(a), int(b)), max(int(a), int(b))) cnt[k] = cnt.get(k, 0) + 1 side = [] for tr in F: for a, b in ((tr[0], tr[1]), (tr[1], tr[2]), (tr[2], tr[0])): a, b = int(a), int(b) if cnt[(min(a, b), max(a, b))] == 1: side.append([a, b, b + n]) side.append([a, b + n, a + n]) if not side: return None # 경계가 없으면(이미 닫힌 면) 압출로 만들 수 없다 F2 = np.vstack([F[:, ::-1], F + n, np.asarray(side, np.int64)]) return V2, F2 def repair_for_tet(V, F, wd): """사다리를 한 단씩 올라가며 'tetgen -p 성공 + 두께 타당' 을 만족하는 형상을 찾는다. 반환: (V, F, 방법이름). 전부 실패하면 (빈배열, 빈배열, '실패'). 디스크 캐시는 (V,F,사다리버전) SHA-1 키. 성공만 캐시한다.""" import hashlib key = hashlib.sha1(np.ascontiguousarray(np.round(V, 4)).tobytes() + np.ascontiguousarray(np.asarray(F, np.int64)).tobytes() + LADDER_VER.encode()).hexdigest()[:16] cp = os.path.join(REPAIR_CACHE, key + '.npz') try: os.makedirs(REPAIR_CACHE, exist_ok=True) if os.path.isfile(cp): z = np.load(cp) return z['V'], z['F'], str(z['m']) except Exception: pass def save(Vx, Fx, meth): try: np.savez_compressed(cp, V=Vx, F=Fx, m=np.array(meth)) except Exception: pass return Vx, Fx, meth emin = float(np.ptp(V, axis=0).min()) if len(V) else 0.0 tries = 0 def accept(Vx, Fx, meth): if not len(Fx) or len(Vx) < 4: return None d = os.path.join(wd, 'rp%d' % tries) os.makedirs(d, exist_ok=True) if not _tet_ok_p(Vx, Fx, d): return None if not _thickness_ok(Vx, Fx, emin): return None return save(Vx, Fx, meth) r = accept(V, F, '원본') if r: return r tries += 1 # A — manifold3d 를 '원본' 에 (순서 중요: 정리를 먼저 하면 수율이 0% 가 된다) try: Vm, Fm = _manifold_pass(V, F) r = accept(Vm, Fm, 'A manifold3d(원본)') if r: return r except Exception: pass tries += 1 try: V1, F1 = _tm_clean(V, F) except Exception: V1, F1 = V, F r = accept(V1, F1, 'B 중복·퇴화면 정리') if r: return r tries += 1 try: V2, F2, na = _earcut_fill(V1, F1) except Exception: V2, F2, na = V1, F1, 0 if na: r = accept(V2, F2, 'C 구멍메움') if r: return r tries += 1 try: V3, F3 = _manifold_pass(V2, F2) r = accept(V3, F3, 'D 메움+manifold3d') if r: return r except Exception: pass tries += 1 # H — 볼록껍질 대체(교량 복구). # ★ 여기에는 두께 게이트를 걸지 않는다. 원본이 닫힌 솔리드가 아니라 기준 부피 자체가 # 쓰레기값이고(껍질이 원본을 포함하는데도 -98% 가 나온다), 게이트를 걸면 교량 21개가 # 전부 탈락해 전체해석이 다시 끊긴다. 대신 '껍질 근사' 라벨을 반드시 남긴다. try: import trimesh as _tm ch = _tm.Trimesh(vertices=V, faces=F, process=False).convex_hull Vh, Fh = np.asarray(ch.vertices), np.asarray(ch.faces) if len(Fh) > HULL_TRI_CAP: try: Vh, Fh = _decimate(Vh, Fh, HULL_TRI_CAP) except Exception: pass d = os.path.join(wd, 'rphull') os.makedirs(d, exist_ok=True) if len(Fh) >= 4 and _tet_ok_p(Vh, Fh, d): return save(Vh, Fh, 'H 볼록껍질대체') except Exception: pass # ★ 실패는 캐시하지 않는다. return np.zeros((0, 3)), np.zeros((0, 3), np.int64), '실패' def build_tet_mesh_repaired(V, F, wd, budget_nodes=None): """전체해석 모드(부품 40개 초과) 전용 메시 생성. ① 수리 사다리로 형상을 고치고 ② tetgen 은 '-p'(표면보존) 전용으로만 쓰고 ③ 부품당 절점 상한을 넘으면 표면 삼각형을 줄여 다시 잰다. ★ 품질조건(-q) 승격을 건너뛰는 근거: -q 가 절점을 51배(최악 1,272배)로 불린다. 얇은 부품일수록 체적상한이 작아져 20삼각형 부품이 27,985절점이 됐다.""" import hashlib cap = min(float(budget_nodes or PER_CAP), float(PER_CAP)) # ★ 완성된 메시까지 캐시한다. '떠 있는 부품 자동 제외' 는 부품을 빼고 다시 푸는데, # 남은 부품의 형상은 하나도 안 바뀌므로 tetgen 을 다시 돌릴 이유가 없다. # 실측: 캐시 없이 146부품 전체해석이 105초였다(2회차 재해석이 절반을 차지). mkey = hashlib.sha1(np.ascontiguousarray(np.round(V, 4)).tobytes() + np.ascontiguousarray(np.asarray(F, np.int64)).tobytes() + ('%s|%d' % (LADDER_VER, int(cap))).encode()).hexdigest()[:16] mcp = os.path.join(REPAIR_CACHE, 'm' + mkey + '.npz') try: os.makedirs(REPAIR_CACHE, exist_ok=True) if os.path.isfile(mcp): z = np.load(mcp) _mm = str(z['m']) return {'nodes': z['nd'], 'elems': z['el'], 'metrics': {'wall_mm': float(z['w']), 'volume_mm3': float(z['v']), 'closed': True, 'bad_edges': 0, 'area_mm2': float(z['a'])}, 'elem_size_mm': float(z['h']), 'through_wall': float(z['tw']), 'coarsened': True, 'quality': str(z['q']), 'flipped': int(z['fl']), 'repair': _mm, 'hull': bool(_mm == 'H 볼록껍질대체')} except Exception: pass Vr, Fr, meth = repair_for_tet(V, F, wd) if not len(Fr): raise RuntimeError('형상 수리 실패(사다리 전부 소진) — 닫힌 솔리드로 만들 수 없습니다') nodes = tets10 = None flipped = 0 tag = '-p전용' for attempt in range(4): pts, tets = run_tetgen(Vr, Fr, wd, 0.0, switches='-p', timeout=30) tets2, flipped = fix_orientation(pts, tets) nodes, tets10 = make_c3d10(pts, tets2) if len(nodes) <= cap or len(Fr) < 80: break tgt = max(60, int(len(Fr) * cap / float(len(nodes)) * 0.85)) try: Vd, Fd = _decimate(Vr, Fr, tgt) except Exception: tag = '-p전용(상한초과 유지)' break d = os.path.join(wd, 'dec%d' % attempt) os.makedirs(d, exist_ok=True) if not len(Fd) or len(Fd) >= len(Fr) or not _tet_ok_p(Vd, Fd, d): tag = '-p전용(상한초과 유지)' break Vr, Fr = Vd, Fd met = mesh_metrics(Vr, Fr) hh = (6.0 * max(met['volume_mm3'], 1.0) / max(len(tets10), 1)) ** (1.0 / 3.0) tw = (met['wall_mm'] / hh) if hh else 0.0 try: np.savez_compressed(mcp, nd=nodes, el=tets10, w=met['wall_mm'], v=met['volume_mm3'], a=met.get('area_mm2', 0.0), h=hh, tw=tw, q=np.array(tag), fl=flipped, m=np.array(meth)) except Exception: pass return {'nodes': nodes, 'elems': tets10, 'metrics': met, 'elem_size_mm': hh, 'through_wall': tw, 'coarsened': True, 'quality': tag, 'flipped': flipped, 'repair': meth, 'hull': (meth == 'H 볼록껍질대체')} # ── 전체해석 '품질 메시' 경로 (0단계 실험으로 확정 — 계획 부록 G) ──────────── # ★ 왜 별도 경로인가: 기본 bulk 경로는 `-p`(표면보존) 전용이라 요소의 43~57% 가 # 수치적으로 쓸 수 없는 납작한 사면체다. 품질조건(`-pq2.0`)을 걸면 그 비율이 # 0.14% 로 떨어지지만 절점이 89만이 되어 메모리 약 20GB·260초가 든다. # 그래서 기본값으로 켜지 않고 '정밀' 선택으로만 켠다(요청의 quality_bulk). # ★ 운전 설정은 실측으로 확정했다(부록 G-7): hs=55mm / 두께게이트 0.80. # hs=75mm 로 더 거칠게 깎으면 절점은 78.6만으로 줄지만 벽이 깎여 # 최대변위가 4.19mm → 55.50mm 로 가짜 유연함이 생긴다. 그 함정에 다시 빠지지 말 것. QBULK_HS_MM = 55.0 # 표면 줄이기 목표 변길이 QBULK_WALL_GATE = 0.80 # 벽두께 보존비 하한(줄인 뒤 벽 ÷ 줄이기 전 벽) QBULK_TETGEN_SW = '-pq2.0' QBULK_TIMEOUT_S = 300.0 # ★ 부품 수와 분리한다(계획 F-8) def build_tet_mesh_quality_bulk(V, F, wd, budget_nodes=None): """전체해석 품질 메시. 수리 → 표면 줄이기(벽두께 게이트) → tetgen 품질조건. 절점 수를 정하는 것은 품질조건이 아니라 벽두께다(계획 F-5 실측). 그래서 손잡이는 '표면 줄이기' 하나이고, 벽이 무너지면 한 단 덜 거칠게 물러난다. 한 번만 시도하면 게이트가 거부한 부품이 원본 밀도를 유지해 상위 5부품이 전체 절점의 39% 를 차지한다(실측) — 그래서 사다리로 물러난다. """ import hashlib mkey = hashlib.sha1(np.ascontiguousarray(np.round(V, 4)).tobytes() + np.ascontiguousarray(np.asarray(F, np.int64)).tobytes() + ('q|%s|%s|%.2f|%.2f' % (LADDER_VER, QBULK_TETGEN_SW, QBULK_HS_MM, QBULK_WALL_GATE)).encode() ).hexdigest()[:16] mcp = os.path.join(REPAIR_CACHE, 'q' + mkey + '.npz') try: os.makedirs(REPAIR_CACHE, exist_ok=True) if os.path.isfile(mcp): z = np.load(mcp, allow_pickle=True) _mm = str(z['m']) return {'nodes': z['nd'], 'elems': z['el'], 'metrics': {'wall_mm': float(z['w']), 'volume_mm3': float(z['v']), 'closed': True, 'bad_edges': 0, 'area_mm2': float(z['a'])}, 'elem_size_mm': float(z['h']), 'through_wall': float(z['tw']), 'coarsened': True, 'quality': str(z['q']), 'flipped': int(z['fl']), 'repair': _mm, 'hull': bool(_mm == 'H 볼록껍질대체')} except Exception: pass Vr, Fr, meth = repair_for_tet(V, F, wd) if not len(Fr): raise RuntimeError('형상 수리 실패(사다리 전부 소진) — 닫힌 솔리드로 만들 수 없습니다') met = mesh_metrics(Vr, Fr) base_tris = len(Fr) tag = QBULK_TETGEN_SW + '(품질)' for k, mult in enumerate((1.0, 0.75, 0.5, 0.35)): hs = QBULK_HS_MM * mult tgt = int(max(met.get('area_mm2', 0.0), 1.0) / (0.433 * hs * hs)) if not (4 <= tgt < int(base_tris * 0.9)): continue # 표면이 이미 그만큼 거칠다 → 줄일 것이 없다 try: Vd, Fd = _decimate(Vr, Fr, tgt) except Exception: continue m1 = mesh_metrics(Vd, Fd) if len(Fd) < 4 or not m1['closed']: continue if m1['wall_mm'] / max(met['wall_mm'], 1e-9) < QBULK_WALL_GATE: continue # 벽이 무너진다 → 거부하고 한 단 덜 거칠게 d = os.path.join(wd, 'qd%d' % k) os.makedirs(d, exist_ok=True) if not _tet_ok_p(Vd, Fd, d): continue Vr, Fr, met = Vd, Fd, m1 tag += '+줄임%.0f%%' % (100.0 * (1.0 - len(Fd) / float(base_tris))) break pts, tets = run_tetgen(Vr, Fr, wd, 0.0, switches=QBULK_TETGEN_SW, timeout=QBULK_TIMEOUT_S) tets2, flipped = fix_orientation(pts, tets) nodes, tets10 = make_c3d10(pts, tets2) hh = (6.0 * max(met['volume_mm3'], 1.0) / max(len(tets10), 1)) ** (1.0 / 3.0) tw = (met['wall_mm'] / hh) if hh else 0.0 try: np.savez_compressed(mcp, nd=nodes, el=tets10, w=met['wall_mm'], v=met['volume_mm3'], a=met.get('area_mm2', 0.0), h=hh, tw=tw, q=np.array(tag), fl=flipped, m=np.array(meth)) except Exception: pass return {'nodes': nodes, 'elems': tets10, 'metrics': met, 'elem_size_mm': hh, 'through_wall': tw, 'coarsened': False, 'quality': tag, 'flipped': flipped, 'repair': meth, 'hull': (meth == 'H 볼록껍질대체')} def build_tet_mesh_local(V, F, wd, budget_nodes=None, quality_timeout=QUALITY_TIMEOUT_S, h_target=None): """국부 재해석(2단계) 전용 — 수리 사다리로 형상을 고친 뒤 **품질 경로**로 메시. 왜 따로 필요한가 (실측 2026-10-10): · build_tet_mesh(품질 경로)만 쓰면 CAD 가 열린 부품이 그대로 실패한다. 국부 영역 25부품 중 20개가 '닫힌 솔리드가 아님' 으로 죽었다. · build_tet_mesh_repaired(전체해석 경로)는 수리는 하지만 부품당 절점이 PER_CAP(4,000)에 하드캡으로 걸려 '요소 크기 절반' 사다리를 만들 수 없다 (두 단계가 같은 메시로 나와 수렴 판정이 무의미해진다). → 수리는 사다리에서, 밀도는 체적상한에서 — 두 경로의 필요한 절반씩만 쓴다. ★ 결함 1 고침 (실측 2026-10-09): 국부 영역에서 부재가 빠지는 원인은 '형상 처리 실패' 가 **아니었다**. 전수 확인한 12부재 × 2단계에서 탈락 사유 10건이 전부 `가장 성긴 설정에서도 절점 N개로 예산(M) 초과` 였다. 국부 경로는 _bulk 가 아니어서 PER_CAP 이 안 걸리는 대신 per = 총예산 × 부피비 로 쪼개지는데, 1단계 총예산 8만을 12부재가 나누면 작은 부재 몫이 PER_NODE_FLOOR(400)까지 내려가 표면보존(-p) 최소 메시조차 못 들어간다. 단계마다 그 경계가 달라 메시 성공 부재 집합이 7→8, 11→12 로 바뀌었다(결함 2 의 원인도 이것). → 국부 경로에서는 부재별 예산 초과를 **탈락 사유로 쓰지 않는다**. 가장 성긴 메시를 그대로 받아들이고(allow_over_budget), 총 절점은 영역 부재 수가 LOCAL_MAX_PARTS 로 묶여 있어 발산하지 않는다. """ Vr, Fr, meth = repair_for_tet(V, F, wd) if not len(Fr): raise RuntimeError('형상 수리 실패(사다리 전부 소진) — 닫힌 솔리드로 만들 수 없습니다') # ── 요소 크기를 **직접** 지정하는 길 (수렴 사다리의 전제) ──────────────── # ★ 왜 절점 예산이 아니라 요소 크기인가 (실측 2026-10-09): 예산을 8배로 올려도 # 부재의 절반 이상은 메시가 한 톨도 안 바뀌었다. build_tet_mesh 에는 # ① `dense_surface(표면삼각형>2000) → -p 결과 그대로 반환` 단축로 # ② `h = max(벽두께/1.2, 예산h)` 하한 # 이 있어서 CAD 부재(표면이 촘촘)는 예산과 무관하게 같은 메시가 나온다. # 그 결과 1단계 요소 5.5mm → 2단계 8.4mm 로 **거꾸로** 굵어졌다. # 수렴 판정은 '요소를 절반으로 줄였다' 를 전제하므로 이 상태로는 뜻이 없다. # → 국부 경로에서는 tetgen 체적상한(-a)을 h 로 직접 걸고, 부재마다 목표 h 와 # 실제 h 를 같이 돌려준다. 못 맞춘 부재는 숨기지 않고 적는다. if h_target and h_target > 0: _a = (float(h_target) ** 3) / 6.0 for sw, tg in (('-pq2.0a%.6f' % _a, '요소지정(품질)'), ('-pa%.6f' % _a, '요소지정')): try: pts, tets = run_tetgen(Vr, Fr, wd, 0.0, switches=sw, timeout=max(quality_timeout, 60.0)) except Exception: continue if not len(tets): continue tets2, flipped = fix_orientation(pts, tets) nd, el = make_c3d10(pts, tets2) met = mesh_metrics(Vr, Fr) hh = (6.0 * max(met['volume_mm3'], 1.0) / max(len(el), 1)) ** (1.0 / 3.0) return {'nodes': nd, 'elems': el, 'metrics': met, 'elem_size_mm': hh, 'through_wall': (met['wall_mm'] / hh if hh else 0.0), 'coarsened': False, 'quality': tg, 'flipped': flipped, 'repair': meth, 'hull': (meth == 'H 볼록껍질대체'), 'h_target_mm': float(h_target), # ★ h_hit 은 'tetgen 이 돌아왔다' 가 아니라 '목표를 실제로 맞췄다' # 여야 한다. 체적상한(-a)은 벽두께가 더 빡빡하면 구속력이 없어 # 목표 19.45mm 를 걸어도 실제가 8.09mm 로 나온다. 그것을 # '전 부재 달성' 이라고 적으면 계산 안 한 것을 표시하는 짓이다. 'h_hit': bool(hh <= float(h_target) * 1.25)} try: tm = build_tet_mesh(Vr, Fr, wd, budget_nodes=budget_nodes, quality_timeout=quality_timeout, allow_over_budget=True) except Exception: # ★ 마지막 후퇴: 전체해석 경로와 **똑같이** tetgen '-p' 만 돌린다. # 수리 사다리가 통과시킨 형상(= tetgen -p 성공이 사다리의 합격 조건)인데 # build_tet_mesh 의 mesh_metrics(closed) 게이트에서 떨어지는 경우가 있다. # 전체해석에서 살아난 부재가 국부에서만 빠지면 후보 자리를 잴 수 없다. pts, tets = run_tetgen(Vr, Fr, wd, 0.0, switches='-p', timeout=60) tets2, flipped = fix_orientation(pts, tets) nd, el = make_c3d10(pts, tets2) met = mesh_metrics(Vr, Fr) hh = (6.0 * max(met['volume_mm3'], 1.0) / max(len(el), 1)) ** (1.0 / 3.0) tm = {'nodes': nd, 'elems': el, 'metrics': met, 'elem_size_mm': hh, 'through_wall': (met['wall_mm'] / hh if hh else 0.0), 'coarsened': True, 'quality': '표면보존만(후퇴)', 'flipped': flipped} tm['repair'] = meth tm['hull'] = (meth == 'H 볼록껍질대체') tm['h_target_mm'] = (float(h_target) if h_target else None) tm['h_hit'] = False if h_target else None return tm def build_tet_mesh(V, F, wd, target_elems_through_wall=1.2, budget_nodes=None, quality_timeout=15.0, allow_over_budget=False): # ★ 벽당 요소 수 목표를 3.0 → 1.2 로 내린 근거(실측): # 각파이프 검증에서 '벽당 0.97겹' 메시가 처짐 -6.99% / 굽힘응력 +2.15% 였다. # 2차 요소(C3D10)는 두께방향 변형률이 1차식이라 굽힘에 한 겹이면 충분하다. # 3.0 은 작은 부재를 쓸데없이 잘게 쪼갠다 — 실측: 186정점 브래킷이 55,561절점(15.8초). """벽두께를 재서 tetgen 체적상한을 정하고 C3D10 메시를 만든다. 절점이 예산을 넘으면 자동으로 성기게 하고 그 사실을 알린다.""" budget_nodes = float(budget_nodes or MAX_C3D10_NODES) met = mesh_metrics(V, F) if not met['closed']: raise RuntimeError('닫힌 솔리드가 아님(열린 엣지 %d개) — 사면체 메시 불가' % met['bad_edges']) t = met['wall_mm'] if not (t > 0): raise RuntimeError('벽두께 추정 실패') # ① 정확도가 원하는 요소크기(벽을 몇 겹으로 나눌지) h_want = t / target_elems_through_wall n_tet_ok = max(budget_nodes / C3D10_NODES_PER_TET, 500.0) h_budget = (6.0 * met['volume_mm3'] * TETGEN_OVERSHOOT / n_tet_ok) ** (1.0 / 3.0) h = max(h_want, h_budget) # ② tetgen 전략 사다리. # ★ 실측으로 확인한 것: 요소 수를 지배하는 것은 체적상한(-a) 이 아니라 품질조건(-q) 이다. # 실제 부재(710x300x9 판재, 표면삼각형 9,084개) 기준 # -pq1.4a102 → 사면체 556,583개 # -pq2.0a102 → 275,993개 # -pa102 → 275,993개 (tetgen 은 -a 를 주면 품질조건을 기본 적용) # -p → 6,870개 (0.2초) # 즉 -a 를 22mm 로 키워도 525,087개라 예산으로 못 줄인다. 그래서 품질조건을 낮추는 # 사다리로 간다: 고품질(작은 형상) → 표준품질 → 표면보존만(-p, CAD 부재처럼 # 표면이 이미 촘촘한 경우 이것만으로도 충분한 밀도가 나온다). # ★ 전략: 먼저 '-p'(표면보존만)로 싸게 재 본다(실측 0.2초). # CAD 에서 변환된 실제 부재는 표면이 이미 촘촘해서(예: 9,084 삼각형) -p 만으로도 # 충분한 밀도가 나오고, 여기에 품질조건을 얹으면 40배로 폭발한다(6,870 → 275,993개). # 반대로 시험용 단순 형상(표면 삼각형 수십 개)은 -p 로는 너무 성기니 품질조건으로 채운다. # → 표면 삼각형 수로 갈라서, 헛되이 비싼 시도를 반복하지 않는다. last = None base = None try: pts, tets = run_tetgen(V, F, wd, 0.0, switches='-p') tets2, flipped = fix_orientation(pts, tets) nodes, tets10 = make_c3d10(pts, tets2) hh = (6.0 * met['volume_mm3'] / max(len(tets2), 1)) ** (1.0 / 3.0) base = {'nodes': nodes, 'elems': tets10, 'metrics': met, 'elem_size_mm': hh, 'through_wall': t / hh, 'coarsened': True, 'quality': '표면보존만(성김)', 'flipped': flipped} except Exception as ex: last = ex dense_surface = len(F) > 2000 if base is not None and len(base['nodes']) > budget_nodes: # ★ allow_over_budget: 국부 재해석 전용. 예산 초과를 탈락 사유로 쓰면 단계마다 # 메시 성공 부재 집합이 달라져 수렴 판정이 불가능해진다(결함 1·2 실측). # 기본값(False)은 종전 그대로 — selftest 의 이론값 대조가 성긴 메시로는 # 통과하지 못하므로 일반 경로의 동작을 건드리지 않는다. if not allow_over_budget: raise RuntimeError('가장 성긴 설정에서도 절점 %d개로 예산(%d) 초과 — 부재가 너무 복잡합니다' % (len(base['nodes']), int(budget_nodes))) base['quality'] = '표면보존만(성김·예산초과 허용)' return base if base is not None and (dense_surface or base['through_wall'] >= 1.2): return base # 표면이 이미 촘촘하다 → 그대로 쓴다(가장 빠름) # ★ 품질 개선(-q)에는 제한시간을 둔다. 이건 '있으면 좋은' 단계지 필수가 아니다. # 실측: RJ102-0020 의 228정점짜리 작은 부재 하나가 -q 단계에서 192.75초를 쓰고도 # 사면체를 342개밖에 못 만들었다(형상이 까다로우면 tetgen 이 헤맨다). # 22부품 조립이 353초 걸린 주범이 이것 — 제한시간을 넘기면 -p 결과로 물러난다. for sw, tag in [('-pq1.4a%.6f' % (h ** 3 / 6.0), '고품질'), ('-pq2.0a%.6f' % (h ** 3 / 6.0), '표준품질')]: try: pts, tets = run_tetgen(V, F, wd, 0.0, switches=sw, timeout=quality_timeout) except subprocess.TimeoutExpired: last = RuntimeError('품질 개선 %.0f초 초과 — 성긴 메시로 진행' % quality_timeout) break except Exception as ex: last = ex continue tets, flipped = fix_orientation(pts, tets) nodes, tets10 = make_c3d10(pts, tets) if len(nodes) > budget_nodes: last = RuntimeError('절점 %d개 > 예산 %d' % (len(nodes), int(budget_nodes))) continue hh = (6.0 * met['volume_mm3'] / max(len(tets), 1)) ** (1.0 / 3.0) return {'nodes': nodes, 'elems': tets10, 'metrics': met, 'elem_size_mm': hh, 'through_wall': t / hh, 'coarsened': tag != '고품질', 'quality': tag, 'flipped': flipped} if base is not None: return base raise RuntimeError('사면체 메시 실패: %s' % (last,)) # ══════════════════════════════════════════════════════════════════ # 3. 외곽 상자(OBB) — 사면체 실패 시 후퇴용 + 자동 하중 방향 결정용 # ══════════════════════════════════════════════════════════════════ def obb_from_points(P): P = np.asarray(P, dtype=float).reshape(-1, 3) P = np.unique(np.round(P, 4), axis=0) c0 = P.mean(axis=0) Q = P - c0 w, Vv = np.linalg.eigh(np.cov(Q.T)) Vv = Vv[:, np.argsort(w)[::-1]] if np.linalg.det(Vv) < 0: Vv[:, 2] *= -1.0 Lp = Q @ Vv lo, hi = Lp.min(axis=0), Lp.max(axis=0) size = hi - lo o2 = np.argsort(size)[::-1] Vv, lo, hi, size = Vv[:, o2], lo[o2], hi[o2], size[o2] if np.linalg.det(Vv) < 0: Vv[:, 2] *= -1.0 lo[2], hi[2] = -hi[2], -lo[2] return c0 + Vv @ ((lo + hi) / 2.0), Vv, size def box_hex_mesh(size, nu, nv, nw): L, b, t = size xs, ys, zs = np.linspace(0, L, nu + 1), np.linspace(0, b, nv + 1), np.linspace(0, t, nw + 1) def nid(i, j, k): return i * (nv + 1) * (nw + 1) + j * (nw + 1) + k + 1 nodes = np.zeros(((nu + 1) * (nv + 1) * (nw + 1), 3)) for i in range(nu + 1): for j in range(nv + 1): for k in range(nw + 1): nodes[nid(i, j, k) - 1] = (xs[i], ys[j], zs[k]) elems = [] for i in range(nu): for j in range(nv): for k in range(nw): elems.append([nid(i, j, k), nid(i + 1, j, k), nid(i + 1, j + 1, k), nid(i, j + 1, k), nid(i, j, k + 1), nid(i + 1, j, k + 1), nid(i + 1, j + 1, k + 1), nid(i, j + 1, k + 1)]) return nodes, np.array(elems, dtype=int) - 1 # ══════════════════════════════════════════════════════════════════ # 3-b. 여러 부재(연결성분) 분리 + 용접 접합 # ══════════════════════════════════════════════════════════════════ def split_components(V, F): """삼각형이 서로 정점을 공유하는 덩어리별로 나눈다 = 부재 1개씩.""" n = len(V) parent = np.arange(n) def find(x): while parent[x] != x: parent[x] = parent[parent[x]] x = parent[x] return x for tri in F: a = find(tri[0]) for k in (1, 2): b = find(tri[k]) if a != b: parent[b] = a roots = np.array([find(i) for i in range(n)]) labels, inv = np.unique(roots, return_inverse=True) comps = [] for ci in range(len(labels)): vmask = inv == ci if not vmask.any(): continue fsel = F[vmask[F[:, 0]]] if len(fsel) < 4: continue idx = np.where(vmask)[0] remap = -np.ones(n, dtype=np.int64) remap[idx] = np.arange(len(idx)) comps.append((V[idx], remap[fsel])) return comps # C3D10 각 면의 (꼭짓점3, 중간절점3) — 중간절점 국부번호 = 4 + C3D10_EDGES 인덱스 # (m01, m12, m20 순서 = 2차 삼각형 형상함수 N4,N5,N6 순서와 일치) C3D10_FACE_NODES = [ (0, 1, 2, 4, 5, 6), (0, 3, 1, 7, 8, 4), (1, 3, 2, 8, 9, 5), (2, 3, 0, 9, 7, 6), ] def boundary_faces6(elems10): """겉면(한 번만 쓰인 면)을 6절점 삼각형으로 반환. 반환: faces(n,6), opp(n,) = 그 면이 속한 사면체의 '반대쪽 꼭짓점'. ★ opp 가 있어야 바깥 방향 법선을 일관되게 정할 수 있다. (감김만으로 정하면 면마다 안팎이 뒤섞여, 용접면 전달력 적분이 상쇄되지 않고 수십 배로 튄다 — 실측)""" blocks, opps = [], [] allc = {0, 1, 2, 3} for f in C3D10_FACE_NODES: blocks.append(elems10[:, list(f)]) o = list(allc - set(f[:3]))[0] opps.append(elems10[:, o]) faces = np.vstack(blocks) opp = np.concatenate(opps) key = np.sort(faces[:, :3], axis=1) _, idx, cnt = np.unique(key, axis=0, return_index=True, return_counts=True) keep = idx[cnt == 1] return faces[keep], opp[keep] def tri6_weights(p, P6): """점 p 를 6절점 삼각형(P6) 에 투영해 2차 형상함수 가중치 6개를 구한다.""" a, b, c = P6[0], P6[1], P6[2] n = np.cross(b - a, c - a) nn = np.dot(n, n) if nn <= 1e-20: return None q = p - n * (np.dot(p - a, n) / nn) # 평면으로 투영 v0, v1, v2 = b - a, c - a, q - a d00, d01, d11 = np.dot(v0, v0), np.dot(v0, v1), np.dot(v1, v1) d20, d21 = np.dot(v2, v0), np.dot(v2, v1) den = d00 * d11 - d01 * d01 if abs(den) < 1e-20: return None l2 = (d11 * d20 - d01 * d21) / den l3 = (d00 * d21 - d01 * d20) / den l1 = 1.0 - l2 - l3 # 삼각형 밖이면 안으로 당긴다(가장자리 절점 대비) L = np.clip(np.array([l1, l2, l3]), 0.0, 1.0) s = L.sum() if s <= 0: return None L = L / s L1, L2, L3 = L return np.array([L1 * (2 * L1 - 1), L2 * (2 * L2 - 1), L3 * (2 * L3 - 1), 4 * L1 * L2, 4 * L2 * L3, 4 * L3 * L1]) def _wloss(n_in=0, spc_drop=0, taken_drop=0, no_face=0, too_far=0, degen=0, eq=0): """용접 씨앗이 '어디서 몇 개 사라졌나' 를 담는 칸. 합이 반드시 닫힌다: n_in = spc_drop + taken_drop + no_face + too_far + degen + eq ★ 왜 taken_drop 과 spc_drop 을 쪼개는가 (2026-10-12): 호출부에서 `fixed_set = set(spc_set)` 로 시작하므로, 버려진 씨앗에는 '원래 구속절점' (고칠 수 없음) 과 '앞 쌍이 먼저 가져감'(배정 결함, 고칠 수 있음) 이 섞인다. 두 수를 한 칸에 담으면 배정을 고쳐도 효과가 측정되지 않는다.""" return {'in': int(n_in), 'spc_drop': int(spc_drop), 'taken_drop': int(taken_drop), 'no_face': int(no_face), 'too_far': int(too_far), 'degen': int(degen), 'eq': int(eq)} def build_weld_equations(nodes, slave_ids, master_faces6, fixed_set, max_dist, fixed_set0=None): """비적합(서로 안 맞는) 메시끼리 용접 접합 = 내가 직접 만드는 MPC(*EQUATION). 반환: (eqs, len(eqs), loss) — loss 는 _wloss 의 분해 표(계측 전용, 해석 불변). ★ CalculiX *TIE 를 쓰지 않는 이유(실측): *TIE 는 종속 절점을 상대 면 위로 '실제로 옮겨서' 주변 사면체가 뒤집히고 'nonpositive jacobian' 으로 해석이 죽는다. *EQUATION 은 절점을 옮기지 않고 변위만 묶으므로 그 문제가 없다. ★★ 후보 상대면을 '면 중심까지의 거리' 로 고르면 안 된다 (실측 2026-10-05) 예전에는 종속절점에서 max_dist 안에 있는 '면 중심' 만 후보로 삼았다. 상대 부재가 성긴 메시면 면 하나가 수백 mm 라서, 종속절점이 그 면 위에 딱 붙어 있어도 면 중심은 수백 mm 떨어져 있다 → 후보가 0개 → 이음이 하나도 안 생긴다. 실측: 500PRO-R 16부품 묶음에서 작은 부재(RE102_0215_1)가 상대면 12개를 가졌는데 종속절점 374+135+46개에 대해 묶인 절점이 0+2+0개였다. 그 부재는 사실상 공중에 떠서 변위가 66,534,314 mm(6,600만 mm)로 나왔다 — 오류 없이 가짜 결과가 나온 것. → 면의 크기(중심~꼭짓점 최대거리)를 반경에 더해 후보를 고르고, 거리는 점↔삼각형 정확거리로 판정한다. 메시 밀도가 달라도 성립한다.""" from scipy.spatial import cKDTree _n_in = len(slave_ids) if not len(slave_ids) or not len(master_faces6): return [], 0, _wloss(n_in=_n_in, no_face=_n_in) _f0 = fixed_set if fixed_set0 is None else fixed_set0 sids, _spcd, _takd = [], 0, 0 for s in slave_ids: s = int(s) if s in fixed_set: if s in _f0: _spcd += 1 # 원래 구속절점 — 배정을 고쳐도 살릴 수 없다 else: _takd += 1 # 앞 쌍이 먼저 가져감 — 배정 결함(채널②) continue sids.append(s) sids = np.array(sids, dtype=np.int64) if not len(sids): return [], 0, _wloss(n_in=_n_in, spc_drop=_spcd, taken_drop=_takd) spts = nodes[sids] tree = cKDTree(spts) T = nodes[master_faces6[:, :3]] # 상대면 꼭짓점 (m,3,3) ctr = T.mean(axis=1) rad = np.linalg.norm(T - ctr[:, None, :], axis=2).max(axis=1) + max_dist cand = tree.query_ball_point(ctr, rad) best_d = np.full(len(sids), np.inf) best_f = -np.ones(len(sids), dtype=np.int64) for fi, idxs in enumerate(cand): if not len(idxs): continue ix = np.asarray(idxs, dtype=np.int64) d = point_tri_dist(spts[ix], T[fi]) upd = d < best_d[ix] if upd.any(): ii = ix[upd] best_d[ii] = d[upd] best_f[ii] = fi eqs = [] _nof = _tfar = _deg = 0 for k in range(len(sids)): fi = int(best_f[k]) if fi < 0: _nof += 1 # 상대면 후보가 아예 없다 continue if best_d[k] > max_dist: _tfar += 1 # 가장 가까운 면이 묶을 거리 한도를 넘는다 continue w = tri6_weights(nodes[sids[k]], nodes[master_faces6[fi]]) if w is None: _deg += 1 # 찌그러진 삼각형 — 가중치를 못 만든다 continue eqs.append((int(sids[k]), master_faces6[fi], w)) return eqs, len(eqs), _wloss(n_in=_n_in, spc_drop=_spcd, taken_drop=_takd, no_face=_nof, too_far=_tfar, degen=_deg, eq=len(eqs)) def point_tri_dist_pairs(P, T): """짝지어진 점↔삼각형 거리를 한꺼번에 계산한다. P:(M,3), T:(M,3,3) → (M,). ★ 왜 이 함수가 따로 필요한가: point_tri_dist 는 '삼각형 1개 대 점 여러 개' 라서 후보 삼각형이 수만 개면 파이썬 반복이 수만 번 돌아 해석보다 오래 걸린다. 용접 전달력 면 선별(B-6)은 '면마다 다른 후보 삼각형' 이라 짝 계산이 맞다. 평면에 내린 발이 삼각형 안이면 수직거리, 밖이면 세 변까지의 거리 중 최소. """ P = np.asarray(P, float) A, B, C = T[:, 0], T[:, 1], T[:, 2] v0, v1 = B - A, C - A n = np.cross(v0, v1) nn = np.einsum('ij,ij->i', n, n) safe = nn > 1e-20 nn_s = np.where(safe, nn, 1.0) w = P - A h = np.einsum('ij,ij->i', w, n) / nn_s q = P - h[:, None] * n wq = q - A d00 = np.einsum('ij,ij->i', v0, v0) d01 = np.einsum('ij,ij->i', v0, v1) d11 = np.einsum('ij,ij->i', v1, v1) den = d00 * d11 - d01 * d01 den_s = np.where(np.abs(den) > 1e-20, den, 1.0) l2 = (d11 * np.einsum('ij,ij->i', wq, v0) - d01 * np.einsum('ij,ij->i', wq, v1)) / den_s l3 = (d00 * np.einsum('ij,ij->i', wq, v1) - d01 * np.einsum('ij,ij->i', wq, v0)) / den_s inside = safe & (np.abs(den) > 1e-20) & (l2 >= -1e-9) & (l3 >= -1e-9) & (l2 + l3 <= 1 + 1e-9) dperp = np.linalg.norm(P - q, axis=1) dedge = np.full(len(P), np.inf) for (U, Vv) in ((A, B), (B, C), (C, A)): e = Vv - U L2 = np.maximum(np.einsum('ij,ij->i', e, e), 1e-20) t = np.clip(np.einsum('ij,ij->i', P - U, e) / L2, 0.0, 1.0) dedge = np.minimum(dedge, np.linalg.norm(P - (U + t[:, None] * e), axis=1)) return np.where(inside, dperp, dedge) def point_tri_dist(p, T): """점 여러 개(p: n,3) → 삼각형 1개(T: 3,3) 까지의 정확한 최단거리. 평면에 내린 발이 삼각형 안이면 그 거리, 밖이면 세 변까지의 거리 중 최소.""" a, b, c = T[0], T[1], T[2] n = np.cross(b - a, c - a) nn = float(np.dot(n, n)) if nn <= 1e-20: # 찌그러진 삼각형 → 변 거리만 d = np.full(len(p), np.inf) for (u, v) in ((a, b), (b, c), (c, a)): e = v - u L2 = float(np.dot(e, e)) t = np.clip(((p - u) @ e) / L2, 0.0, 1.0) if L2 > 1e-20 else np.zeros(len(p)) d = np.minimum(d, np.linalg.norm(p - (u + t[:, None] * e), axis=1)) return d h = ((p - a) @ n) / nn q = p - h[:, None] * n # 평면으로 내린 발 v0, v1 = b - a, c - a d00, d01, d11 = np.dot(v0, v0), np.dot(v0, v1), np.dot(v1, v1) w = q - a d20, d21 = w @ v0, w @ v1 den = d00 * d11 - d01 * d01 if abs(den) < 1e-20: l2 = np.zeros(len(p)); l3 = np.zeros(len(p)) else: l2 = (d11 * d20 - d01 * d21) / den l3 = (d00 * d21 - d01 * d20) / den inside = (l2 >= -1e-9) & (l3 >= -1e-9) & (l2 + l3 <= 1.0 + 1e-9) d = np.linalg.norm(p - q, axis=1) # 평면까지 수직거리 if inside.all(): return d de = np.full(len(p), np.inf) for (u, v) in ((a, b), (b, c), (c, a)): e = v - u L2 = float(np.dot(e, e)) t = np.clip(((p - u) @ e) / L2, 0.0, 1.0) if L2 > 1e-20 else np.zeros(len(p)) de = np.minimum(de, np.linalg.norm(p - (u + t[:, None] * e), axis=1)) return np.where(inside, d, de) def contact_nodes(slave_pts, master_tris, tol, tree_slave=None): """종속 부재의 겉면 절점 중 '상대 부재 표면 위(또는 tol 이내)' 에 있는 것을 고른다. 반환: (절점 인덱스 배열, 쓰인 상대 삼각형 인덱스 집합) ★★ 왜 절점↔절점이 아니라 절점↔표면인가 (실측 근거, 2026-10-05) 예전에는 두 부재의 겉면 '절점끼리' tol 이내인지로 접촉을 판정했다. 그런데 두 부재의 메시는 서로 맞춰 만든 게 아니라서(비적합 메시), 면이 완전히 겹쳐 붙어 있어도 절점은 상대 절점 사이 한가운데에 떨어진다. 절점 간격이 20mm면 최근접 절점까지 10mm 이상 나오므로 tol(=최소벽두께x1.2, 6.6mm) 안에 아무도 안 들어온다. → 실측: 500PRO-R 블럭에서 표면간 거리 1mm 이내로 '진짜 붙어 있는' 쌍 16개 중 절점↔절점 방식은 4개만 찾아냈고, 15부품이 11덩어리로 갈라졌다(하중 전달 불가). 절점이 상대 '표면' 위에 있는지는 메시 밀도와 무관하게 성립한다. 그래서 거리 기준 tol 도 비로소 '용접 가능한 틈(mm)' 이라는 물리적 뜻을 갖는다.""" from scipy.spatial import cKDTree if not len(slave_pts) or not len(master_tris): return np.array([], dtype=np.int64), set() tree = tree_slave if tree_slave is not None else cKDTree(slave_pts) ctr = master_tris.mean(axis=1) # 삼각형 중심에서 꼭짓점까지 최대거리 = 그 삼각형이 영향을 줄 수 있는 반경 rad = np.linalg.norm(master_tris - ctr[:, None, :], axis=2).max(axis=1) + tol cand = tree.query_ball_point(ctr, rad) hit = np.zeros(len(slave_pts), dtype=bool) used = set() for fi, idxs in enumerate(cand): if not len(idxs): continue ix = np.asarray(idxs, dtype=np.int64) ix = ix[~hit[ix]] if hit.any() else ix if not len(ix): continue d = point_tri_dist(slave_pts[ix], master_tris[fi]) ok = d <= tol if ok.any(): hit[ix[ok]] = True used.add(fi) return np.where(hit)[0], used def find_weld_pairs(parts, tol): """부재 쌍마다 '서로 맞닿은 곳'(용접선)을 찾는다. parts[i] = {'nodes','elems','n0','e0'} (n0/e0 = 전체 모델에서의 번호 시작점)""" from scipy.spatial import cKDTree welds = [] # ★ 부품이 많아지면 쌍의 수가 n^2 로 폭발한다(291부품 → 42,195쌍). # 쌍마다 KD트리를 새로 만들면 몇십 분이 걸린다. 그래서 # ① 부품별 KD트리를 한 번만 만들어 재사용하고 # ② 경계상자가 tol 안에서 겹치지도 않는 쌍은 트리를 보기 전에 버린다. # (경계상자가 안 겹치면 표면 절점도 절대 tol 안에 들어올 수 없다 — 안전한 선별) bb = [] trees = [] for pp in parts: sp = pp['surf_pts'] bb.append((sp.min(axis=0), sp.max(axis=0))) trees.append(cKDTree(sp)) for i in range(len(parts)): for j in range(i + 1, len(parts)): if not (np.all(bb[i][0] - tol <= bb[j][1]) and np.all(bb[j][0] - tol <= bb[i][1])): continue # ★★ 종속(slave) 방향은 '양쪽 다 시도해서 접촉절점이 많은 쪽' 으로 고른다. # 왜: 예전에는 '겉면 절점 총수가 많은 쪽' 을 종속으로 고정했다. 그런데 총수는 # 국부 밀도의 대리지표가 못 된다. 1225x635mm 큰 판재는 절점이 844개여도 # 45mm 짜리 작은 브래킷 삼각형 위에 올라가는 절점은 0~2개뿐이라 탈락했고, # 반대 방향(브래킷 8절점이 판재의 큰 삼각형 위)이면 전부 잡힌다. # 실측: 단방향 111쌍 → 양방향 243쌍(2.2배), 공중에 뜬 부품 83개 → 21개. # ★ 요구 접촉절점 수도 3 → 1 로 내렸다. 절점 1~2개 접촉은 '점용접' 이지만 # 그래도 힘을 전달한다. 3개를 요구하면 작은 브래킷이 전부 떨어져 나간다 # (실측: 용접쌍 193 → 243, 공중에 뜬 부품 26 → 21). _got = [] for (_si, _mi) in ((i, j), (j, i)): _A, _B = parts[_si], parts[_mi] _mt = _B['nodes'][_B['faces6'][:, :3] - _B['n0']] _ai, _uf = contact_nodes(_A['surf_pts'], _mt, tol, tree_slave=trees[_si]) if len(_ai) >= 1: _got.append((len(_ai), _si, _mi, _ai, _uf)) if not _got: continue # ★ 재현성: 접촉절점 수가 같으면 종전에는 (i,j) 가 이겨서 **부품 순서**가 # 종속 방향을 정했다. 동수일 때는 부품번호가 아니라 '종속 부재의 접촉 # 절점 좌표합' 으로 가른다 - CAD 값이라 순서와 무관하다. def _tiekey(t): _A = parts[t[1]] _pts = _A['surf_pts'][t[3]] return (-t[0], round(float(_pts.sum()), 3), int(len(_pts))) _got.sort(key=_tiekey) _, si, mi, ai, used_faces = _got[0] # 진 방향도 들고 간다 — 방정식이 0개로 나오면 역할을 뒤집어 재시도하기 위해. _rev = _got[1] if len(_got) > 1 else None A, B = parts[si], parts[mi] mt = B['nodes'][B['faces6'][:, :3] - B['n0']] # 상대 부재 겉면 삼각형 an = A['surf_ids'][ai] zone_pts = [A['nodes'][an - A['n0']]] if used_faces: zone_pts.append(mt[sorted(used_faces)].reshape(-1, 3)) zone = np.vstack(zone_pts) ext = zone.max(axis=0) - zone.min(axis=0) welds.append({ 'i': si, 'j': mi, 'dep_nodes': an, # 종속 부재의 접촉 절점 'ind_part': mi, 'near_faces': np.array(sorted(used_faces), dtype=np.int64), 'zone_min': zone.min(axis=0), 'zone_max': zone.max(axis=0), 'center': zone.mean(axis=0), 'length_mm': float(np.sort(ext)[-1]), # 용접선 길이 ≈ 접촉영역 최장변 'spot': bool(len(ai) < 3), # 접촉절점 1~2개 = 점용접 # 역할 뒤집기 재시도용 (종속=mi 쪽 절점, 상대면=si 쪽 삼각형) 'dep_nodes_rev': (parts[_rev[1]]['surf_ids'][_rev[3]] if _rev is not None else np.array([], dtype=np.int64)), 'near_faces_rev': (np.array(sorted(_rev[4]), dtype=np.int64) if _rev is not None else np.array([], dtype=np.int64)), }) return welds WELD_SEED_PITCH_MM = 25.0 # 용접 씨앗점 격자 간격(mm) — CAD 기준, 메시와 무관 WELD_SEED_MIN = 3 # 접합면 하나당 최소 씨앗점 수(3점이면 면·모멘트 전달) WELD_SEED_PITCH_FLOOR = 4.0 # 격자를 더 못 쪼개는 한계 def _lattice_pick(P, pitch): """점구름을 전역 격자(원점 0, 간격 pitch)로 나누고 칸마다 대표점 1개만 남긴다. ★ 왜 격자인가: CAD 삼각형마다 점을 뽑으면 씨앗 수가 'CAD 분할 방식' 에 따라 달라진다. 전역 격자에 스냅하면 메시도 CAD 분할도 바꿔도 같은 씨앗이 나온다.""" if not len(P): return P key = np.floor(np.asarray(P, float) / pitch).astype(np.int64) _, inv = np.unique(key, axis=0, return_inverse=True) inv = np.asarray(inv).ravel() out = [] for c in range(int(inv.max()) + 1): idx = np.where(inv == c)[0] m = P[idx].mean(axis=0) out.append(int(idx[np.argmin(((P[idx] - m) ** 2).sum(axis=1))])) return P[np.array(sorted(out), dtype=np.int64)] def _pts_near_tris(pts, T, tol, tree_pts=None): """점구름 중 '삼각형 집합의 표면에서 tol 이내' 인 점의 인덱스. 완전 벡터화. contact_nodes 와 판정은 같지만 삼각형마다 파이썬 반복을 돌지 않는다 — CAD 삼각형은 수만 개라 반복이 돌면 씨앗 계산이 해석보다 오래 걸린다.""" from scipy.spatial import cKDTree pts = np.asarray(pts, float) if not len(pts) or not len(T): return np.array([], dtype=np.int64) tree = tree_pts if tree_pts is not None else cKDTree(pts) ctr = T.mean(axis=1) rad = np.linalg.norm(T - ctr[:, None, :], axis=2).max(axis=1) + tol cand = tree.query_ball_point(ctr, rad) pi, ti = [], [] for fi, idxs in enumerate(cand): if not len(idxs): continue pi.append(np.asarray(idxs, dtype=np.int64)) ti.append(np.full(len(idxs), fi, dtype=np.int64)) if not pi: return np.array([], dtype=np.int64) pi = np.concatenate(pi); ti = np.concatenate(ti) out = [] for k in range(0, len(pi), 2000000): # 메모리 상한 — 토막내서 계산 sl = slice(k, k + 2000000) d = point_tri_dist_pairs(pts[pi[sl]], T[ti[sl]]) out.append(pi[sl][d <= tol]) return np.unique(np.concatenate(out)) if out else np.array([], dtype=np.int64) def _assign_unique(pool, pts, kmax=8): """씨앗점마다 서로 겹치지 않는 최근접 절점 1개씩을 배정한다 (결함 2). ★ 왜 필요한가 (실측 2026-10-10): 씨앗 3,425점 중 507점(14.8%)이 다른 씨앗과 같은 절점에 떨어졌다. 종전 코드는 `np.unique` 로 그걸 삼켜 묶인절점이 씨앗 수보다 적어졌고, 메시가 달라지면 삼켜지는 수도 달라져 단계 간 '같은 묶음' 이 깨졌다. → 가까운 순서(k=1..kmax)로 보며 이미 쓴 절점은 건너뛰고 다음 절점에 배정한다. kmax 까지 다 찼으면 어쩔 수 없이 최근접을 공유한다(그 수를 세어 보고한다). 반환: (거리(n,), 절점자리번호(n,)) """ from scipy.spatial import cKDTree pool = np.asarray(pool, float) pts = np.asarray(pts, float).reshape(-1, 3) kk = int(min(max(kmax, 1), len(pool))) d, i = cKDTree(pool).query(pts, k=kk) d = np.asarray(d, float).reshape(len(pts), kk) i = np.asarray(i, np.int64).reshape(len(pts), kk) # 가까운 씨앗부터 배정해, 먼 씨앗이 밀려나게 한다(총 이동거리 최소화에 가깝다) taken = set() outd = d[:, 0].copy() outi = i[:, 0].copy() for s in np.argsort(d[:, 0]): for c in range(kk): if int(i[s, c]) not in taken: taken.add(int(i[s, c])) outd[s], outi[s] = d[s, c], i[s, c] break return outd, outi def weld_seeds_cad(parts, tol, pitch=WELD_SEED_PITCH_MM): """용접 씨앗점을 **메시를 만들기 전의 CAD 표면**에서 1회만 정한다 (1단계). ★ 왜 필요한가 (실측 근거): 종전에는 메시를 만든 뒤 '표면 절점이 상대 표면에 닿았는가' 로 종속절점을 골랐다. 그래서 메시를 촘촘하게 하면 묶이는 절점이 43% 늘어나 '메시를 키웠더니 결과가 바뀌는' 비교가 불가능했다. 또 묶인 절점이 전체의 29.7% 를 덮어, 구속에서 떨어진 순수 구조 후보가 하나도 드러나지 않았다. → 씨앗점은 CAD 삼각형끼리의 거리로만 정하므로 메시 밀도와 무관하다. 어떤 메시든 씨앗점마다 가장 가까운 절점 1개씩만 종속으로 쓰므로 묶음 수가 고정된다. 반환: (seeds, stat) seeds[(i,j)] = (n,3) 씨앗점 좌표 (i= lo_, axis=1) & np.all(tl <= hi_, axis=1))[0] sp = np.where(np.all(samp[si] >= lo_, axis=1) & np.all(samp[si] <= hi_, axis=1))[0] return ft, sp fti, spi = _sel(i, j) # i 의 삼각형·표본점을 j 경계상자로 선별 ftj, spj = _sel(j, i) ai = (spi[_pts_near_tris(samp[i][spi], tris[j][ftj], tol)] if len(spi) and len(ftj) else np.array([], dtype=np.int64)) aj = (spj[_pts_near_tris(samp[j][spj], tris[i][fti], tol)] if len(spj) and len(fti) else np.array([], dtype=np.int64)) if not len(ai) and not len(aj): continue # 양쪽 접촉 표본점을 합쳐 접촉영역을 만든다(한쪽 CAD 가 성길 수 있다) cloud = np.vstack([x for x in (samp[i][ai] if len(ai) else None, samp[j][aj] if len(aj) else None) if x is not None]) # ★ 재현성 [2026-10-11]: cloud 는 i WELD_SEED_PITCH_FLOOR and len(cloud) > len(pts): pit = max(pit / 2.0, WELD_SEED_PITCH_FLOOR) pts = _lattice_pick(cloud, pit) ref += 1 if ref: stat['refined_pairs'] += 1 seeds[(i, j)] = pts stat['pairs'] += 1 stat['seeds'] += len(pts) stat['pair_list'].append((parts[i]['part_no'], parts[j]['part_no'], len(pts), round(pit, 1))) return seeds, stat # ── 떠 있는 덩어리 '다리 잇기' ────────────────────────────────────────────── # 왜 필요한가 (실측 2026-10-11, 머신 BC 6고정/-Z/3톤, 146부품): # 고정 덩어리와 이어지지 않아 자동 제외된 11부품은 서로 떨어진 11개가 아니라 # 3덩어리였고(47·49·50 / 77 / 90·115·117·118·119·120·121), 고정 쪽으로 가는 # 다리는 딱 3군데였다 — 실측 점↔삼각형 최소거리로: # 49/50 ↔ 18번 = 2.920mm : 씨앗 4·5개가 이미 생겼는데 MPC 가 0개였다. # 원인은 간격이 아니다. MPC 생성 루프가 weld_list(메시 접촉쌍)만 돌고, # 이 쌍은 거기에 없어서(CAD 씨앗만 있는 쌍) 통째로 건너뛰어졌다. # 90 ↔ 91번 = 5.635mm : 접촉쌍은 있는데 tol 3.0mm 밖이라 씨앗이 0개. # 77 ↔ 17번 = 8.000mm : 씨앗도 접촉쌍도 없음. # → 이 세 가지를 한 가지 방법으로 고친다: '떠 있는 덩어리' 의 부재에서 실측 # 표면간격이 BRIDGE_MAX_GAP_MM 이내인 쌍만 골라, 그 쌍만 간격+1mm 로 씨앗을 # 다시 뽑아 용접 방정식을 만든다. 전역 tol 은 한 숫자도 건드리지 않는다 # (전역 tol 을 올리면 모델 전체의 묶음이 다시 정해져 기준값과 비교가 깨진다). # ★ 10mm 를 넘는 쌍은 잇지 않는다 — 그것은 CAD 에서 실제로 떨어져 있는 부품이고, # 이으면 없는 용접을 만들어 내는 것이다. 제외 사유에 실측 간격을 적어 남긴다. BRIDGE_MAX_GAP_MM = 10.0 # 이 이상 떨어진 쌍은 '실제 분리' 로 보고 잇지 않는다 BRIDGE_TOL_MARGIN_MM = 1.0 # 쌍별 tol = 실측 간격 + 이 여유(씨앗이 표면에 걸리도록) # ── 부재별 메시 밀도 상향 (연결 실패가 '간격' 이 아니라 '메시' 문제일 때) ──────── # 왜 필요한가 (실측 2026-10-11, 머신 BC): 49·50·77 은 CAD 간격이 각각 2.92·2.92· # 8.00mm 로 이을 범위(10mm) 안인데도 다리가 안 걸렸다. 사유는 간격이 아니었다 — # 그 부재의 겉면 절점이 성겨서 접합면 자리에 절점이 아예 없었다(겉면 절점 ↔ 상대 # 메시 최소거리 11.00 / 11.00 / 22.42mm. 중간절점까지 후보에 넣어도 9.19mm). # → 그 부재 **하나만** 요소 크기를 내려 다시 메시하면 겉면 절점이 접합면에 생긴다. # 실측(_exp_fine_probe.py, 상대 CAD 표면 기준): # 목표 h 실제 h 49번 절점 겉면절점↔상대CAD # 2.0mm 1.58mm 491,957 0.007mm # 4.0 3.07 73,548 0.006 # 6.0 4.45 25,945 0.006 # 12.0 7.79 5,592 2.774 # 요구는 '간격+여유' 안에 절점이 들어오는 것이고, 전체 절점 100만 상한이 있으므로 # **요구를 만족하는 가장 거친 설정** 을 쓴다. 부재당 상한을 두고 사다리로 물러난다. FINE_PART_NODE_CAP = 30000 # 밀도 상향 부재 1개의 절점 상한(100만 예산을 지키려고) FINE_H_LADDER = (1.0, 1.5, 2.0, 3.0, 4.0) # 목표 h 배수 — 상한을 넘으면 거칠게 물러난다 def build_tet_mesh_fine(V, F, wd, h_req, node_cap=FINE_PART_NODE_CAP, logf=None): """부재 1개를 요소크기 h_req 로 다시 메시한다. 절점 상한을 넘으면 거칠게 물러난다. ★ 숫자를 맞추려는 조정이 아니다 — h_req 는 '이 부재가 상대 표면에서 간격+여유 안에 절점을 가져야 한다' 는 요구에서 나온 값(= 실측 간격 + 여유)이다. 상한은 전체 절점 100만을 넘기지 않기 위한 예산이고, 물러난 사실은 적어 남긴다. """ last = None for mult in FINE_H_LADDER: h = float(h_req) * mult d = os.path.join(wd, 'f%.0f' % (mult * 10)) os.makedirs(d, exist_ok=True) try: tm = build_tet_mesh_local(V, F, d, budget_nodes=1e9, quality_timeout=180.0, h_target=h) except Exception: continue tm['fine_h_req_mm'] = float(h_req) tm['fine_h_mult'] = float(mult) tm['fine_capped'] = bool(len(tm['nodes']) > node_cap) last = tm if len(tm['nodes']) <= node_cap: if logf: logf(' [밀도상향] 요구 h %.2fmm × %.1f = %.2fmm → 실제 %.2fmm / 절점 %s' % (h_req, mult, h, tm['elem_size_mm'], '{:,}'.format(len(tm['nodes'])))) return tm return last def cad_min_gap(Vi, Fi, Vj, Fj, k=24): """부품 i 표면 ↔ 부품 j 표면의 최소 간격(mm). 점↔삼각형 정확식. 표본점은 i 의 꼭짓점 + 삼각형 무게중심. 표본점마다 j 의 삼각형 중 무게중심이 가까운 k개만 정확거리를 재고 그 최소를 취한다 — 전수(수만×수만)는 돌지 않는다. k 를 키우면 값은 단조 감소하므로 이 값은 참 최소거리의 상한이다. """ from scipy.spatial import cKDTree Ti, Tj = Vi[Fi], Vj[Fj] S = np.vstack([Vi, Ti.mean(axis=1)]) C = Tj.mean(axis=1) kk = int(min(k, len(C))) if not len(S) or not kk: return float('inf') _d, idx = cKDTree(C).query(S, k=kk) idx = np.asarray(idx).reshape(len(S), kk) best = float('inf') for col in range(kk): best = min(best, float(point_tri_dist_pairs(S, Tj[idx[:, col]]).min())) return best def bridge_floating(parts, nodes, grp, fix_groups, fixed_set, tie_blocks, base_tol, logf=None, max_gap=BRIDGE_MAX_GAP_MM, fixed_set0=None): """떠 있는 덩어리를 실측 간격 근거로만 이어 붙인다. tie_blocks 에 바로 더한다. 반환: 이은 다리 목록 [{'a','b','gap_mm','tol_mm','seeds','eqs','kind'}] + 못 이은 쌍의 실측 간격(사유 표기용) ★ 순서 의존성: 이 함수가 만드는 방정식도 fixed_set('한 절점은 한 번만 종속')을 누적하므로 도는 순서가 결과를 바꿀 수 있다. 성분은 '부품번호 최솟값', 쌍은 '(부품번호, 부품번호)' 로 정렬해 입력 순서와 무관하게 만든다. """ lo = np.array([p['lo'] for p in parts], float) hi = np.array([p['hi'] for p in parts], float) comps = {} for i, g in enumerate(grp): comps.setdefault(int(g), []).append(i) flo = sorted((g for g in comps if g not in fix_groups), key=lambda g: min(parts[i]['part_no'] for i in comps[g])) bridges, rejected = [], [] for g in flo: mem = sorted(comps[g], key=lambda i: parts[i]['part_no']) memset = set(mem) cand = [] for a in mem: d = np.maximum(np.maximum(lo - hi[a], lo[a] - hi), 0.0).max(axis=1) for b in np.where(d <= max_gap)[0]: b = int(b) if b in memset: continue cand.append((int(parts[a]['part_no']), int(parts[b]['part_no']), a, b)) cand.sort() for (na, nb, a, b) in cand: A, B = parts[a], parts[b] gap = cad_min_gap(np.asarray(A['cadV'], float), np.asarray(A['cadF'], np.int64), np.asarray(B['cadV'], float), np.asarray(B['cadF'], np.int64)) if not np.isfinite(gap) or gap > max_gap: rejected.append({'a': na, 'b': nb, 'gap_mm': (round(gap, 3) if np.isfinite(gap) else None)}) continue tol = max(float(base_tol), float(gap) + BRIDGE_TOL_MARGIN_MM) sd, _st = weld_seeds_cad([A, B], tol) pts = sd.get((0, 1)) if pts is None or not len(pts): rejected.append({'a': na, 'b': nb, 'gap_mm': round(gap, 3), 'why': '씨앗 0개'}) continue zlo, zhi = pts.min(axis=0) - tol, pts.max(axis=0) + tol eqs, cnt, roles = [], 0, [] _bl, _bside = None, None for (S, M) in ((A, B), (B, A)): # 한쪽 절점이 이미 소비됐으면 역할을 뒤집는다 # ★ 상대면 선별에 '면 중심이 상자 안' 을 쓰면 안 된다(실측으로 당했다): # 성긴 메시는 면 하나가 수백 mm 라, 씨앗 바로 위에 덮인 면이라도 # 중심은 상자 밖이다 → 후보 0개 → 방정식 0개. # (2026-10-11 실측: 49↔18·50↔18·77↔17 세 쌍이 전부 이것 때문에 0개였다) # → 면마다 '중심~꼭짓점 최대거리' 를 상자에 더해 고른다. 그 거리보다 # 가까운 면은 전부 담기므로 빠뜨리지 않는다. _mt = nodes[M['faces6'][:, :3]] _mc = _mt.mean(axis=1) _mr = np.linalg.norm(_mt - _mc[:, None, :], axis=2).max(axis=1)[:, None] sel = np.where(np.all(_mc >= zlo - _mr, axis=1) & np.all(_mc <= zhi + _mr, axis=1))[0] # ★ 실측값을 남긴다: 겉면 절점이 상대 '메시 표면' 에서 얼마나 # 떨어져 있는지. CAD 간격(gap)과 이 값이 다르면 원인은 간격이 # 아니라 '메시가 성겨 접합면에 절점이 없다' 는 것이다. _sp = nodes[np.asarray(S['surf_ids'], np.int64)] _ng = (cad_min_gap(_sp, np.arange(len(_sp)).reshape(-1, 1).repeat(3, 1), nodes, M['faces6'][:, :3]) if len(_sp) else float('inf')) roles.append({'slave_n': int(len(_sp)), 'sel_faces': int(len(sel)), 'node_gap_mm': (round(_ng, 3) if np.isfinite(_ng) else None)}) if not len(sel): continue # ★ 종속절점을 '씨앗점의 최근접 절점' 으로 고르면 안 된다(실측 2026-10-11): # 전체해석 메시는 요소가 수십 mm 라 씨앗에서 가장 가까운 절점도 # 수십 mm 떨어져 있다. 그 절점은 상대 표면에서 간격(2.9~8mm)보다 # 훨씬 멀어, 묶을 거리 한도를 넘겨 방정식이 0개가 됐다 # (49↔18 / 50↔18 / 77↔17 세 쌍 전부 이것 때문이었다). # → 거꾸로 한다: 겉면 절점 전부를 후보로 주고 '상대 표면에서 # 간격+여유 안' 인 절점만 묶게 한다(거리 판정은 점↔삼각형 정확식). # 그러면 묶는 거리가 실측 간격 수준으로 작아진다. _all, _n, _bl0 = build_weld_equations( nodes, S['surf_ids'], M['faces6'][sel], fixed_set, float(gap) + BRIDGE_TOL_MARGIN_MM, fixed_set0) if not _n: continue # 메시 무관성 유지: 씨앗점마다 그 후보 중 1개씩만 남긴다 # (전 접합면을 다 묶으면 메시가 촘촘해질 때 묶음 수가 늘어난다) _cid = np.array([int(e[0]) for e in _all], dtype=np.int64) _dk, _ik = _assign_unique(nodes[_cid], pts) _keep = set(int(_cid[k]) for k in _ik) eqs = [e for e in _all if int(e[0]) in _keep] cnt = len(eqs) if cnt: # ★ 종속절점이 어느 부품 쪽인가(=힘의 부호 기준). 1단계에서 # 전체/국부 접합면 힘을 비교할 때 이 값이 없으면 부호를 # 맞출 수 없다. _bl, _bside = _bl0, int(S['part_no']) break if not cnt: rejected.append({'a': na, 'b': nb, 'gap_mm': round(gap, 3), 'roles': roles, 'why': ('씨앗 %d개인데 간격+%.0fmm 안에 묶을 절점이 없음 ' '(겉면 절점↔상대 메시 최소거리 %s mm)' % (len(pts), BRIDGE_TOL_MARGIN_MM, '/'.join(str(r['node_gap_mm']) for r in roles) or '?'))}) continue for e in eqs: fixed_set.add(e[0]) tie_blocks.extend(eqs) bridges.append({ 'a': na, 'b': nb, 'gap_mm': round(float(gap), 3), 'tol_mm': round(tol, 3), 'seeds': int(len(pts)), 'eqs': int(cnt), # ★ 두 원인을 한 라벨로 뭉개지 않는다. 간격이 기준 tol 안이면 # '간격 추정' 이 아니라 '접촉쌍 누락' 이 사실이다. 'kind': ('접촉쌍 누락 복구' if gap <= base_tol else '간격 추정 용접'), # 계측(해석 불변): 손실 분해 + 종속절점이 어느 부품 쪽인가. # assign_cut = 방정식은 만들었지만 '씨앗 1개당 절점 1개' 규칙으로 # 되돌려 버린 수(메시 무관성 유지 비용) — loss['eq'] 와 cnt 의 차. 'loss': _bl, 'assign_cut': int((_bl['eq'] if _bl else 0) - cnt), 'side_part_no': _bside, # ★ 1단계: 다리도 '접합면 힘' 에 들어가야 한다. 종전 _iface_rf 는 # weld_list 만 돌아 다리 방정식이 통째로 빠져 있었고, 그래서 # 국부 재해석이 '기준값 없는 쌍' 을 강제하고 있었다. 'tied_ids': np.array([int(e[0]) for e in eqs], dtype=np.int64), }) if logf: logf(' [다리] %d번↔%d번 간격 %.3fmm → tol %.2fmm, 씨앗 %d, 방정식 %d (%s)' % (na, nb, gap, tol, len(pts), cnt, bridges[-1]['kind'])) return bridges, rejected def interface_perimeter(nodes, faces3): """접합면 삼각형 집합의 바깥 둘레 길이 = 용접선 길이.""" e = np.vstack([faces3[:, [0, 1]], faces3[:, [1, 2]], faces3[:, [2, 0]]]) k = np.sort(e, axis=1) u, cnt = np.unique(k, axis=0, return_counts=True) b = u[cnt == 1] if not len(b): return 0.0 return float(np.linalg.norm(nodes[b[:, 0]] - nodes[b[:, 1]], axis=1).sum()) def weld_interface_force(nodes, stress, part, welds_w, tol): """용접 접합면을 지나 전달되는 힘 = 접촉면 삼각형에 작용하는 트랙션 t=σ·n 의 적분. (반력 추출이 아니라 응력장에서 직접 구하므로 CalculiX 옵션에 의존하지 않는다)""" pts = nodes sel_idx = welds_w.get('force_idx') if sel_idx is None or not len(sel_idx): return None g = part['faces6'][sel_idx][:, :3] opp = part['faces_opp'][sel_idx] p0, p1, p2 = pts[g[:, 0]], pts[g[:, 1]], pts[g[:, 2]] cr = np.cross(p1 - p0, p2 - p0) ln = np.linalg.norm(cr, axis=1) area = 0.5 * ln nrm = cr / np.maximum(ln[:, None], 1e-12) ctr = (p0 + p1 + p2) / 3.0 # 반대쪽 꼭짓점에서 면 쪽으로 = 바깥 방향. 어긋난 면은 뒤집는다. flip = np.einsum('ij,ij->i', nrm, ctr - pts[opp]) < 0 nrm[flip] *= -1.0 sel = np.ones(len(g), dtype=bool) # ★ 6절점(2차) 삼각형의 면적분은 '중간절점 3개의 평균' 이 정확하다. # 2차 형상함수를 삼각형 위에서 적분하면 꼭짓점 가중치는 0, 중간절점은 각각 1/3 이다. # 꼭짓점 평균을 쓰면 굽힘처럼 기울기가 센 곳에서 적분이 틀어진다. gm = part['faces6'][sel_idx][:, 3:] S = stress[gm].mean(axis=1) sx, sy, sz, txy, tyz, tzx = S.T n1, n2, n3 = nrm[:, 0], nrm[:, 1], nrm[:, 2] t = np.stack([sx * n1 + txy * n2 + tzx * n3, txy * n1 + sy * n2 + tyz * n3, tzx * n1 + tyz * n2 + sz * n3], axis=1) F = (t[sel] * area[sel, None]).sum(axis=0) return {'force_N': F, 'area_mm2': float(area[sel].sum())} # ══════════════════════════════════════════════════════════════════ # 4. 표면 삼각형 추출 (화면 표시용) # ══════════════════════════════════════════════════════════════════ TET_FACES = [(0, 2, 1), (0, 1, 3), (1, 2, 3), (2, 0, 3)] HEX_FACES = [(0, 3, 2, 1), (4, 5, 6, 7), (0, 1, 5, 4), (1, 2, 6, 5), (2, 3, 7, 6), (3, 0, 4, 7)] def surface_tris_tet(tets): """한 번만 쓰인 면 = 겉면.""" faces = np.vstack([tets[:, list(f)] for f in TET_FACES]) key = np.sort(faces, axis=1) _, idx, cnt = np.unique(key, axis=0, return_index=True, return_counts=True) return faces[idx[cnt == 1]] def surface_tris_hex(hexes): quads = np.vstack([hexes[:, list(f)] for f in HEX_FACES]) key = np.sort(quads, axis=1) _, idx, cnt = np.unique(key, axis=0, return_index=True, return_counts=True) q = quads[idx[cnt == 1]] return np.vstack([q[:, [0, 1, 2]], q[:, [0, 2, 3]]]) # ══════════════════════════════════════════════════════════════════ # 5. .inp 작성 / .frd 파싱 # ══════════════════════════════════════════════════════════════════ def write_inp(path, nodes, elems, etype, fix_nodes, loads, E, nu, ties=None, plastic_fy=None, prescribed=None): o = ['*NODE, NSET=Nall'] for i, p in enumerate(nodes, 1): o.append('%d, %.6f, %.6f, %.6f' % (i, p[0], p[1], p[2])) o.append('*ELEMENT, TYPE=%s, ELSET=Eall' % etype) for i, e in enumerate(elems, 1): o.append('%d, %s' % (i, ', '.join(str(int(x) + 1) for x in e))) # ── 용접 접합 = *EQUATION (내가 만든 MPC). 절점을 옮기지 않아 요소가 안 뒤집힌다 ── for (sn, fnodes, w) in (ties or []): for dof in (1, 2, 3): terms = [(sn + 1, dof, 1.0)] for k in range(6): if abs(w[k]) > 1e-9: terms.append((int(fnodes[k]) + 1, dof, -float(w[k]))) o.append('*EQUATION') o.append('%d' % len(terms)) line = [] for ti, (nn, dd, cc) in enumerate(terms): line.append('%d, %d, %.10g' % (nn, dd, cc)) if len(line) == 4 or ti == len(terms) - 1: o.append(', '.join(line) + (',' if ti != len(terms) - 1 else '')) line = [] o += ['*MATERIAL, NAME=MAT', '*ELASTIC', '%.1f, %.3f' % (E, nu)] if plastic_fy: # 완전탄소성에 가까운 약간의 가공경화(수렴 안정용) o += ['*PLASTIC', '%.1f, 0.0' % plastic_fy, '%.1f, 0.20' % (plastic_fy * 1.15)] o += ['*SOLID SECTION, ELSET=Eall, MATERIAL=MAT', '*BOUNDARY'] for n in sorted(fix_nodes): o.append('%d, 1, 3' % (n + 1)) if plastic_fy: o += ['*STEP, INC=200', '*STATIC', '0.1, 1.0, 1e-5, 0.5'] else: o += ['*STEP', '*STATIC'] # ── 강제변위(prescribed displacement) — 스텝 안에서만 적용되는 *BOUNDARY ── # 기존 전체고정(SPC, 3DOF)과 달리 1개 DOF만 지정값으로 미는 것이라 # 별도 블록으로 둔다(용접 MPC 충돌 방지는 호출부에서 fixed_set 처리로 해결). if prescribed: o.append('*BOUNDARY') for n in sorted(prescribed): for dof, val in prescribed[n]: o.append('%d, %d, %d, %.8f' % (n + 1, dof, dof, val)) if loads: o.append('*CLOAD') for n, vec in loads.items(): for d in range(3): if abs(vec[d]) > 1e-12: o.append('%d, %d, %.8f' % (n + 1, d + 1, vec[d])) # ★ RF(반력)를 함께 출력한다. .frd 에 FORC 블록이 생기지만 parse_frd 는 블록 이름으로 # 분기하고 그 외는 mode=None 으로 건너뛰므로 기존 파싱이 깨지지 않는다(출력 전용). o += ['*NODE FILE', 'U, S, RF', '*END STEP'] with open(path, 'w', encoding='ascii') as f: f.write('\n'.join(o) + '\n') def _frd_fields(line, count): vals = [] for i in range(count): s = line[13 + 12 * i: 13 + 12 * (i + 1)] if not s.strip(): break vals.append(float(s)) return int(line[3:13]), vals def parse_frd(path, nnode): """.frd 읽기 — 변위(DISP) / 응력(STRESS) / 반력(FORC). ★ FORC(반력)를 읽는 이유 = 구조해석 실무의 가장 값싸고 강한 1차 검증이 'Σ반력 = Σ하중' 이기 때문이다(F-1). write_inp 가 '*NODE FILE U,S,RF' 로 이미 출력하고 있었는데 읽지 않아 검증 장치가 비어 있었다. ★ 블록 이름·형식은 실측 확인(job.frd): ' -4 DISP 4 1' / ' -4 STRESS 6 1' / ' -4 FORC 4 1' / ' -4 ERROR 1 1'. FORC 는 DISP 와 같은 12자 필드 3개다. 모르는 이름은 mode=None 으로 건너뛰므로 기존 파싱이 깨지지 않는다. """ disp = np.zeros((nnode, 3)) stress = np.zeros((nnode, 6)) forc = np.zeros((nnode, 3)) err = np.full(nnode, np.nan) mode = None with open(path, 'r', errors='replace') as f: for line in f: if line.startswith(' -4'): tok = line[5:].split() nm = tok[0] if tok else '' mode = ('U' if nm == 'DISP' else 'S' if nm == 'STRESS' else 'R' if nm == 'FORC' else 'E' if nm == 'ERROR' else None) elif line.startswith(' -3'): mode = None elif line.startswith(' -1') and mode: if mode == 'U': n, v = _frd_fields(line, 3) if 1 <= n <= nnode: disp[n - 1] = v[:3] elif mode == 'R': n, v = _frd_fields(line, 3) if 1 <= n <= nnode: forc[n - 1] = v[:3] elif mode == 'E': # ★ ERROR / STR(%) = ccx 가 스스로 추정한 절점별 응력오차(%). # 요청하지 않아도 ccx 가 써 준다 — 추가 계산비용 0. # 종전에는 블록 이름을 모른다고 건너뛰어 이 값을 버렸다(계획 A-6). n, v = _frd_fields(line, 1) if 1 <= n <= nnode and v: err[n - 1] = v[0] else: n, v = _frd_fields(line, 6) if 1 <= n <= nnode: stress[n - 1] = v[:6] return disp, stress, forc, err def _grow_one_layer(seed, elems, nnode): """씨앗 절점이 들어 있는 요소의 모든 절점으로 1겹 넓힌다(F-2 제외 마스크용). 왜 '요소 1겹' 인가: 용접 MPC 종속절점에 생기는 가짜 첨두는 그 절점 하나가 아니라 그 절점이 속한 요소 전체에 퍼진다(2차 요소라 중간절점까지). 반경을 지어내는 대신 메시 자체의 연결로 넓히는 것이 메시 밀도가 바뀌어도 일관된다. """ out = np.zeros(nnode, dtype=bool) if not len(elems) or not seed.any(): return out E = np.asarray(elems, dtype=np.int64) hit = seed[E].any(axis=1) if hit.any(): out[np.unique(E[hit].ravel())] = True return out def tet_quality(nodes, elems): """요소 품질을 전수 계산한다(계획 B-7 / D-1-1). 사면체 꼭짓점 4개만 쓴다. ★ 왜 이것이 필요한가: 전체해석 메시의 43~57% 가 상용 해석기라면 거부당할 납작한 사면체(슬리버)인데, 코드에 재는 곳이 없어 지금까지 보이지 않았다. 그 상태에서 나온 응력·안전율을 '판정' 으로 쓰면 거짓말이 된다. · sj(스케일드 야코비안) = sqrt(2)·det(e1,e2,e3)/(|e1||e2||e3|) 의 네 꼭짓점 최소값. 정사면체 = 1, 완전 납작 = 0. 상용 해석기 경고선 ≈ 0.1, 거부선 ≈ 0.02. · ar(종횡비) = 최장변/최단변. 좋음 < 3. """ E = np.asarray(elems) if not len(E): return {} E = E[:, :4].astype(np.int64) P = np.asarray(nodes, float)[E] # (M,4,3) sj = np.full(len(E), np.inf) for a, b, c, d in ((0, 1, 2, 3), (1, 2, 0, 3), (2, 0, 1, 3), (3, 2, 1, 0)): e1, e2, e3 = P[:, b] - P[:, a], P[:, c] - P[:, a], P[:, d] - P[:, a] det = np.einsum('ij,ij->i', np.cross(e1, e2), e3) nrm = (np.linalg.norm(e1, axis=1) * np.linalg.norm(e2, axis=1) * np.linalg.norm(e3, axis=1)) sj = np.minimum(sj, np.sqrt(2.0) * np.abs(det) / np.maximum(nrm, 1e-30)) pairs = ((0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)) L = np.stack([np.linalg.norm(P[:, i] - P[:, j], axis=1) for i, j in pairs], axis=1) ar = L.max(axis=1) / np.maximum(L.min(axis=1), 1e-30) vol6 = np.einsum('ij,ij->i', np.cross(P[:, 1] - P[:, 0], P[:, 2] - P[:, 0]), P[:, 3] - P[:, 0]) # ★ 요소별 실측 크기. 종전 meshinfo 의 elem_size_mm 은 '부품별 평균의 평균' # 이라 실제 요소크기의 2.6배로 나왔다(계획 F-1 — 세분화 사다리가 무력해진 원인). h_e = np.abs(vol6) ** (1.0 / 3.0) return {'elems': int(len(E)), 'sj_median': round(float(np.median(sj)), 4), 'sj_min': float('%.3g' % float(sj.min())), 'sj_bad_pct': round(float(100.0 * np.mean(sj < 0.02)), 4), 'ar_median': round(float(np.median(ar)), 2), 'ar_p99': round(float(np.percentile(ar, 99)), 1), 'ar_max': round(float(ar.max()), 1), 'neg_vol': int(np.sum(vol6 <= 0.0)), 'elem_size_median_mm': round(float(np.median(h_e)), 3), '_sj': sj, '_ar': ar} def von_mises(S): sx, sy, sz, txy, tyz, tzx = S.T return np.sqrt(0.5 * ((sx - sy) ** 2 + (sy - sz) ** 2 + (sz - sx) ** 2) + 3.0 * (txy ** 2 + tyz ** 2 + tzx ** 2)) # ══════════════════════════════════════════════════════════════════ # 5-b. 취약 후보 자동 탐지 (계획 D-2 / CEO 요구 3) # ══════════════════════════════════════════════════════════════════ # ★ 왜 이 함수가 필요한가: 종전 `hotspots` 는 '마스크 후 von Mises 상위 3점을 # 서로 떨어뜨려 고른 것' 이다 = CEO 가 금지한 "최대 응력값 하나로 위험 확정". # 여기서는 여러 지표의 '합의' 로 순위를 정하고, 각 후보에 왜 뽑혔는지와 # 믿을 수 있는지를 숫자로 붙여 내려보낸다. 값을 '위험' 으로 확정하지 않는다. # # 쓰는 지표 6개 (CEO 지시): # ① von Mises — 크기. 단, 특이점 자리는 따로 떼어 낸다 # ② 응력 구배 |∇σ| — 첨두가 '한 점' 인지 '영역' 인지 가른다. 노치·용접토우에서 크다 # ③ 변형에너지밀도 SED — 응력·변형의 곱. 응력 성분 전부를 쓰므로 단축 편향이 없다 # ④ ERR% — 낮을수록 그 자리 결과를 믿을 수 있다(순위에 섞지 않는다) # ⑤ 용접 이음부 전달력 — 적분량이라 메시에 둔감. 하중경로를 보여준다 # ⑥ 피로(FAT80) 순위 — 용접 이음부 기준. 기존 화면 로직과 같은 기준을 서버에서도 쓴다 STRESS_SINGULAR_NEAR = 2.0 # 특이점(고정·하중·MPC) 근처 판정 배수 — 요소크기의 몇 배까지 MESH_DEFECT_PCT = 5.0 # 후보 영역의 불량요소 비율이 이 값을 넘으면 '메시 결함' def strain_energy_density(S, E, nu): """응력만으로 변형에너지밀도(N·mm/mm³)를 계산한다. ★ 선형 등방 탄성에서는 변형률을 따로 출력받지 않아도 응력에서 정확히 나온다: u = 1/(2E)[σx²+σy²+σz² − 2ν(σxσy+σyσz+σzσx) + 2(1+ν)(τxy²+τyz²+τzx²)] (`*EL FILE, E` 를 추가한 것이 아니다 — 응력에서 유도한 값이다.) """ sx, sy, sz, txy, tyz, tzx = np.asarray(S, float).T return (1.0 / (2.0 * E)) * (sx * sx + sy * sy + sz * sz - 2.0 * nu * (sx * sy + sy * sz + sz * sx) + 2.0 * (1.0 + nu) * (txy * txy + tyz * tyz + tzx * tzx)) def _nodal_gradient(nodes, elems, val): """요소 안 절점값의 최대 차이 ÷ 요소 크기 = 요소별 구배. 절점에는 최대값을 뿌린다. ★ 좌표 반경으로 이웃을 찾지 않고 '메시 연결' 로만 계산한다 — 반경을 지어내면 메시 밀도가 바뀔 때 값이 따라 바뀌어 비교가 안 된다. """ E4 = np.asarray(elems)[:, :4].astype(np.int64) out = np.zeros(len(nodes)) if not len(E4): return out v = val[E4] rng = v.max(axis=1) - v.min(axis=1) P = nodes[E4] L = np.stack([np.linalg.norm(P[:, i] - P[:, j], axis=1) for i, j in ((0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3))], axis=1) g = rng / np.maximum(L.mean(axis=1), 1e-9) np.maximum.at(out, E4.ravel(), np.repeat(g, 4)) return out def _pct_rank(v, mask): """mask 인 절점들 안에서의 백분위 순위(0~1). mask 밖은 0.""" out = np.zeros(len(v)) idx = np.where(mask)[0] if not len(idx): return out o = np.argsort(np.argsort(v[idx])) out[idx] = o / max(len(idx) - 1, 1) return out def _cluster_seeds(elems, seed, nnode): """씨앗 절점을 '요소를 공유하면 한 덩어리' 로 묶는다(메시 연결 기반 클러스터링). 반환: 절점수 많은 순으로 정렬된 [절점번호 배열, ...]""" from scipy.sparse import coo_matrix from scipy.sparse.csgraph import connected_components E4 = np.asarray(elems)[:, :4].astype(np.int64) if not len(E4) or not seed.any(): return [] keep = seed[E4] rows, cols = [], [] for a in range(4): for b in range(a + 1, 4): m = keep[:, a] & keep[:, b] if m.any(): rows.append(E4[m, a]); cols.append(E4[m, b]) sid = np.where(seed)[0] if not rows: return [np.array([i]) for i in sid] r = np.concatenate(rows); c = np.concatenate(cols) g = coo_matrix((np.ones(len(r)), (r, c)), shape=(nnode, nnode)) n, lab = connected_components(g, directed=False) groups = {} for i in sid: groups.setdefault(int(lab[i]), []).append(int(i)) out = [np.array(v) for v in groups.values()] out.sort(key=len, reverse=True) return out def detect_candidates(nodes, elems, vm, stress, errpct, excl, sj_elem, E, nu, fy, parts, bc_spheres, mpc_seed, welds, h_elem, top_n=8, seed_frac=0.004): """취약 후보 영역을 사용자 입력 없이 뽑는다. 반환: (구조후보목록, 특이점후보목록, 요약) 순위는 '단일 최대응력' 이 아니라 ①②③ 백분위 순위의 평균(지표 합의)이다. ④ERR% 와 요소품질은 순위에 섞지 않고 '신뢰도' 로만 붙인다 — 섞으면 "오차가 커서 순위가 올라갔다" 같은 뜻 없는 결과가 나온다. """ nn = len(nodes) sed = strain_energy_density(stress, E, nu) grad = _nodal_gradient(nodes, elems, vm) E4 = np.asarray(elems)[:, :4].astype(np.int64) # 절점별 '이웃 요소 중 불량(sj<0.02) 비율' — 후보가 메시 결함인지 가린다 bad_cnt = np.zeros(nn); all_cnt = np.zeros(nn) if sj_elem is not None and len(sj_elem) == len(E4): np.add.at(all_cnt, E4.ravel(), 1.0) np.add.at(bad_cnt, E4.ravel(), np.repeat((sj_elem < 0.02).astype(float), 4)) bad_frac = bad_cnt / np.maximum(all_cnt, 1.0) # 절점 → 부품 대응표 pid = np.full(nn, -1, dtype=int) pnm = {} for pp in parts: pid[pp['n0']:pp['n1']] = int(pp.get('part_no', -1)) pnm[int(pp.get('part_no', -1))] = pp.get('part_name', '') # 특이점(고정·하중·MPC)까지의 거리 — '위험' 으로 올리지 않고 숫자로 적어 보여준다 d_bc = np.full(nn, np.inf) for (_p, _r) in (bc_spheres or []): d_bc = np.minimum(d_bc, np.linalg.norm(nodes - np.asarray(_p, float), axis=1) - float(_r)) d_mpc = np.full(nn, np.inf) if mpc_seed is not None and mpc_seed.any(): from scipy.spatial import cKDTree d_mpc, _ = cKDTree(nodes[mpc_seed]).query(nodes) def run(pop, kind): if not pop.any(): return [] r1 = _pct_rank(vm, pop); r2 = _pct_rank(grad, pop); r3 = _pct_rank(sed, pop) cons = (r1 + r2 + r3) / 3.0 cons[~pop] = -1.0 k = max(60, int(pop.sum() * seed_frac)) thr = float(np.partition(cons[pop], -min(k, int(pop.sum())))[-min(k, int(pop.sum()))]) seed = pop & (cons >= thr) out = [] for g in _cluster_seeds(elems, seed, nn)[:top_n * 3]: if len(g) < 3: continue P = nodes[g] ctr = P.mean(axis=0) rad = float(np.linalg.norm(P - ctr, axis=1).max()) ps = sorted(set(int(x) for x in pid[g] if x >= 0)) ev = errpct[g][np.isfinite(errpct[g])] bf = float(100.0 * bad_frac[g].mean()) dbc = float(np.min(d_bc[g])) if np.isfinite(d_bc[g]).any() else None dmp = float(np.min(d_mpc[g])) if np.isfinite(d_mpc[g]).any() else None near = STRESS_SINGULAR_NEAR * max(h_elem, 1.0) if bf > MESH_DEFECT_PCT: cls, why_cls = '메시 결함', ('이 영역 요소의 %.1f%% 가 쓸 수 없는 품질(sj<0.02)이라 ' '응력값이 형상 때문인지 메시 때문인지 가릴 수 없습니다' % bf) elif kind == 'singular': cls, why_cls = '특이점 의심', ('고정·하중·용접 구속 자리(또는 그 1요소 이내)입니다 — ' '여기 응력은 계산상 인공 첨두일 수 있습니다') elif (dbc is not None and dbc <= near) or (dmp is not None and dmp <= near): cls, why_cls = '특이점 의심', ('가장 가까운 구속까지 %s mm (요소 %.1fmm 의 %.1f배 이내)' % (('%.1f' % min([x for x in (dbc, dmp) if x is not None])) if (dbc is not None or dmp is not None) else '-', h_elem, STRESS_SINGULAR_NEAR)) else: cls, why_cls = '구조 후보', '구속·껍질 자리에서 떨어진 본체이고 요소 품질도 정상입니다' # 가장 가까운 용접 이음부와 그 전달력 (⑤ 하중경로) wnear = None if welds: wd_ = [(float(np.linalg.norm(np.asarray(w['at'], float) - ctr)), w) for w in welds if w.get('at')] if wd_: dd, w = min(wd_, key=lambda t: t[0]) wnear = {'dist_mm': round(dd, 1), 'label': w.get('label'), 'force_N': round(float(w.get('force_N') or 0.0), 1), 'fatigue_pct': round(float(w.get('pct') or 0.0), 1)} out.append({ 'kind': cls, 'at': [round(float(x), 2) for x in ctr], 'radius_mm': round(rad, 1), 'nodes': int(len(g)), 'parts': ps[:12], 'part_names': [pnm.get(x, '') for x in ps[:12]], 'score': round(float(cons[g].max()), 4), 'vonmises_max_MPa': round(float(vm[g].max()), 2), 'vonmises_med_MPa': round(float(np.median(vm[g])), 2), 'grad_max_MPa_per_mm': round(float(grad[g].max()), 3), 'sed_max_Nmm_per_mm3': float('%.4g' % float(sed[g].max())), 'errpct_median': round(float(np.median(ev)), 1) if len(ev) else None, 'bad_elem_pct': round(bf, 2), 'dist_to_bc_mm': (round(dbc, 1) if dbc is not None and np.isfinite(dbc) else None), 'dist_to_weld_mpc_mm': (round(dmp, 1) if dmp is not None and np.isfinite(dmp) else None), 'nearest_weld': wnear, 'yield_ratio_pct': round(float(vm[g].max() / fy * 100.0), 1), 'class_reason': why_cls, # ★ '근거' 는 문장으로 남긴다 — 숫자만 주면 왜 뽑혔는지 아무도 모른다. 'why': ('von Mises 상위 %.0f%% · 응력구배 상위 %.0f%% · 변형에너지밀도 상위 %.0f%% ' '(세 지표 합의 %.2f). 솔버 응력오차 %s%%, 이 영역 불량요소 %.2f%%.' % (100 * (1 - r1[g].max()), 100 * (1 - r2[g].max()), 100 * (1 - r3[g].max()), float(cons[g].max()), ('%.1f' % np.median(ev)) if len(ev) else '없음', bf)), 'verdict': '판정 보류', # 수렴 확인 전에는 안전/위험을 확정하지 않는다 }) out.sort(key=lambda d: d['score'], reverse=True) return out[:top_n] main = run(~excl, 'struct') sing = run(excl, 'singular') summary = { 'method': ('지표 합의(von Mises 순위 + 응력구배 순위 + 변형에너지밀도 순위의 평균) ' '+ 메시 연결 기반 클러스터링. 단일 최대응력으로 순위를 정하지 않는다.'), 'sed_note': '변형에너지밀도는 응력에서 유도했다(변형률 출력은 요청하지 않는다)', 'seed_frac': seed_frac, 'elem_size_mm': round(float(h_elem), 3), 'pop_struct': int((~excl).sum()), 'pop_singular': int(excl.sum()), 'singular_near_factor': STRESS_SINGULAR_NEAR, 'mesh_defect_gate_pct': MESH_DEFECT_PCT, 'verdict_policy': ('후보에는 안전율·위험을 붙이지 않는다. 국부 재해석으로 수렴을 ' '확인한 뒤에만 값을 말한다(CEO 지시 7).'), } return main, sing, summary # ══════════════════════════════════════════════════════════════════ # 6. 경계조건 — 클릭 지점 → 절점 집합 # ══════════════════════════════════════════════════════════════════ def interp_disp(gnodes, gelems, gdisp, pts, k=24): """전체 해석의 변위장을 임의 좌표에서 보간한다 (2단계 국부 재해석의 경계 강제용). ★ 왜 보간인가: 국부 메시는 절점이 새로 생기므로 전체 해석의 절점번호를 쓸 수 없다. 점을 품은 사면체를 찾아 꼭짓점 4개의 무게중심좌표로 섞는다(C3D10 의 모서리절점은 쓰지 않는 1차 보간 — 경계면에서의 오차는 요소 크기의 제곱에 비례해 작다). 사면체를 못 찾으면 최근접 절점 값으로 대체하고 **그 개수를 따로 보고한다** (어떤 방식이 몇 개에 쓰였는지 숨기면 결과를 검증할 수 없다). 반환: (disp(m,3), mode(m,)) mode 1=사면체 보간, 0=최근접 절점 대체""" from scipy.spatial import cKDTree pts = np.asarray(pts, float) corn = np.asarray(gelems, np.int64)[:, :4] gnodes = np.asarray(gnodes, float) gdisp = np.asarray(gdisp, float) out = np.zeros((len(pts), 3)) mode = np.zeros(len(pts), dtype=np.int8) if not len(pts): return out, mode ctr = gnodes[corn].mean(axis=1) kk = int(min(k, len(ctr))) _, ii = cKDTree(ctr).query(pts, k=kk) ii = np.asarray(ii).reshape(len(pts), kk) todo = np.ones(len(pts), dtype=bool) for c in range(kk): idx = np.where(todo)[0] if not len(idx): break e = ii[idx, c] A = gnodes[corn[e]] v0, v1, v2 = A[:, 1] - A[:, 0], A[:, 2] - A[:, 0], A[:, 3] - A[:, 0] q = pts[idx] - A[:, 0] det = np.einsum('ij,ij->i', v0, np.cross(v1, v2)) safe = np.abs(det) > 1e-12 dd = np.where(safe, det, 1.0) b1 = np.einsum('ij,ij->i', q, np.cross(v1, v2)) / dd b2 = np.einsum('ij,ij->i', v0, np.cross(q, v2)) / dd b3 = np.einsum('ij,ij->i', v0, np.cross(v1, q)) / dd b0 = 1.0 - b1 - b2 - b3 tb = -1e-6 ins = safe & (b0 >= tb) & (b1 >= tb) & (b2 >= tb) & (b3 >= tb) if ins.any(): jj = idx[ins] W = np.stack([b0[ins], b1[ins], b2[ins], b3[ins]], axis=1) out[jj] = np.einsum('ik,ikj->ij', W, gdisp[corn[e[ins]]]) mode[jj] = 1 todo[jj] = False if todo.any(): # ★ 결함 3 고침: 종전에는 '최근접 절점 값 복사' 로 때웠다(실측 2~27%). # 절점 하나의 값을 그대로 쓰면 경계 변위가 요소 한 칸만큼 계단처럼 틀어져 # 국부 단면력이 전체와 어긋난다. 씨앗점은 CAD 표면 위에 있어 사면체 # '안' 에 안 들어가는 경우가 많은데, 그래도 가장 가까운 사면체의 **면 위** # 에는 있다. 그래서 무게중심좌표를 단체(simplex)로 투영해 그 사면체의 # 1차 보간식으로 값을 만든다 — 절점 복사가 아니라 면 위 보간이다. idx = np.where(todo)[0] e = ii[idx, 0] A = gnodes[corn[e]] v0, v1, v2 = A[:, 1] - A[:, 0], A[:, 2] - A[:, 0], A[:, 3] - A[:, 0] q = pts[idx] - A[:, 0] det = np.einsum('ij,ij->i', v0, np.cross(v1, v2)) dd = np.where(np.abs(det) > 1e-12, det, 1.0) b1 = np.einsum('ij,ij->i', q, np.cross(v1, v2)) / dd b2 = np.einsum('ij,ij->i', v0, np.cross(q, v2)) / dd b3 = np.einsum('ij,ij->i', v0, np.cross(v1, q)) / dd W = np.clip(np.stack([1.0 - b1 - b2 - b3, b1, b2, b3], axis=1), 0.0, 1.0) W = W / np.maximum(W.sum(axis=1, keepdims=True), 1e-12) out[idx] = np.einsum('ik,ikj->ij', W, gdisp[corn[e]]) mode[idx] = 2 return out, mode def nodes_near(nodes, pt, radius): d = np.linalg.norm(nodes - np.asarray(pt, float), axis=1) return np.where(d <= radius)[0] def auto_cantilever(nodes, center, R, size): """클릭 조건이 없을 때의 기본값: 긴 축 한쪽 끝 고정 + 반대 끝에 두께방향 하중.""" local = (nodes - center) @ R u = local[:, 0] lo, hi = u.min(), u.max() band = max((hi - lo) * 0.02, size[2] * 0.5) fix = np.where(u <= lo + band)[0] tip = np.where(u >= hi - band)[0] return fix, tip, R[:, 2] # ══════════════════════════════════════════════════════════════════ # 6-b. 실패한 부품을 사람 말로 설명하기 # ══════════════════════════════════════════════════════════════════ def _pname(part_names, pi): """부품 이름. 이름을 안 받았으면 번호만.""" try: nm = (part_names or [])[pi] except Exception: nm = None return str(nm) if nm else ('%d번 부품' % (pi + 1)) def _why_mesh_failed(ex, met, Vc): """tetgen/메시 실패 사유를 현장에서 알아들을 한국어로 바꾼다. ★ 사유를 뭉개면 '어느 도면을 고쳐야 하는지' 를 알 수 없다. 형상 자체의 결함과 규모 초과를 구분해서 적는다.""" s = str(ex) span = (Vc.max(axis=0) - Vc.min(axis=0)) if len(Vc) else np.zeros(3) if float(span.min()) <= 1e-6 or met.get('volume_mm3', 0.0) <= 1e-6: return ('두께가 0인 면 형상입니다(크기 %.0f x %.0f x %.0f mm) — 속을 채울 수 없어 해석 대상이 아닙니다' % tuple(span)) if 'self-intersection' in s or 'vertex lies' in s or 'PLC Error' in s or 'Two facets' in s: return '형상이 스스로 교차합니다(CAD 변환 결함) — 그 부품만 모델을 고쳐야 합니다' if '닫힌 솔리드가 아님' in s: return s + ' — 면이 빠진 형상입니다' if '예산' in s or '절점' in s: return s if '벽두께 추정 실패' in s: return '벽두께를 재지 못했습니다(부피 또는 표면적이 0)' return s def _failed_text(failed_parts): """화면에 그대로 띄울 한 줄 요약.""" if not failed_parts: return '' head = ', '.join('%d번(%s): %s' % (f['no'], f['name'], f['why']) for f in failed_parts[:6]) more = (' …외 %d개' % (len(failed_parts) - 6)) if len(failed_parts) > 6 else '' return '형상 처리 실패 %d개 — %s%s' % (len(failed_parts), head, more) def _part_bbox(P, part_sizes, pi): """원래 보내온 삼각형 수프에서 pi번(0-base) 부품의 경계상자.""" off = 0 for k, cnt in enumerate(part_sizes or []): rows = int(cnt) // 3 if k == pi: seg = P[off:off + rows] if not len(seg): return None return seg.min(axis=0), seg.max(axis=0) off += rows return None def _no_node_msg(kind, pts, failed_parts, part_sizes, P): """'그 자리에 절점이 없다' 를 원인별로 갈라 말한다. ★ 같은 증상에 원인이 둘이고 운전자가 할 일이 다르다. ① 찍은 자리의 부품이 형상 처리에 실패해 아예 메시에서 빠졌다 → 그 부품 이름을 말한다 ② 그 자리에 형상 자체가 없다(프리셋이 외곽상자 모서리를 짚었다) → 직접 찍으라고 한다""" hit = [] if failed_parts and part_sizes and pts: for f in failed_parts: bb = _part_bbox(P, part_sizes, f['no'] - 1) if bb is None: continue lo, hi = bb mg = max(float(np.max(hi - lo)) * 0.05, 5.0) for q in pts: p = np.asarray(q.get('p'), float) if p.size == 3 and np.all(p >= lo - mg) and np.all(p <= hi + mg): hit.append(f) break if hit: return ('%s이 있는 %s에서 형상 처리가 실패해 그 부품이 해석에서 빠졌습니다 (%s). ' '다른 부재에 %s을 찍어주세요.' % (kind, ', '.join('%d번 부품(%s)' % (f['no'], f['name']) for f in hit), hit[0]['why'], kind)) base = ('이 자리에는 형상이 없습니다 — %s을 부재 표면에 직접 찍어주세요. ' '(각파이프·채널·L형강은 바깥 치수 모서리가 허공입니다)' % kind) if failed_parts: base += ' ※ ' + _failed_text(failed_parts) return base # ══════════════════════════════════════════════════════════════════ # 7. 해석 본체 # ══════════════════════════════════════════════════════════════════ LOCAL_MAX_PARTS = 12 # 국부 영역 부품 수 상한 (반경 R 을 이 수로 자동 결정) # ★ 25 → 12 로 줄인 근거(실측 2026-10-10): 25부품이면 절점 예산을 다 써도 요소가 # 33~45mm 로 씨앗 격자(25mm)보다 굵어진다. 그러면 씨앗 여러 개가 한 절점에 # 겹쳐 떨어지고, 메시를 조이면 그 겹침이 풀려 묶인절점이 70 → 89개로 늘었다 # (= 경계조건이 달라져 수렴 판정 불가). 영역을 줄여 두 단계 모두 요소를 # 25mm 보다 작게 만든다 — 숫자를 맞추려는 조정이 아니라 비교가 성립하게 하는 조건. LOCAL_DRIVE_R_MM = 0.0 # 0 = 씨앗점마다 최근접 절점 1개만 강제(메시 무관) LOCAL_BUDGETS = (80000, 640000) # (구) 절점 예산 사다리 — 아래 이유로 쓰지 않는다 # ★ 수렴 사다리는 '요소 크기' 로 올린다(실측 2026-10-09). # 절점 예산을 8배로 올려도 CAD 부재 절반은 메시가 안 바뀌어 요소크기가 # 5.5mm → 8.4mm 로 거꾸로 굵어졌다. 수렴 판정의 전제('요소를 절반으로')가 # 깨지므로 tetgen 체적상한으로 요소 크기를 직접 지정한다. # 배수는 전체 해석 요소크기(h_global) 기준: 2h → h. # 왜 h/2 까지 안 내려가는가: 12부재를 h/2 로 깎으면 절점이 160만을 넘어 # '100만 절점 초과 금지' 를 깬다. 2h → h 도 체적당 요소수는 8배(=변 1/2)다. # ★★ 실측(2026-10-09)으로 한 번 더 고쳤다: 영역 12부재 **전체**를 2h→h 로 # 깎으면 절점이 83만 → 234만이 되어 '100만 절점 초과 금지' 를 정면으로 깬다. # → 세분화는 **후보가 있는 부재에만** 건다. 둘레 부재는 두 단계 모두 # 똑같은 메시를 쓴다. 이것이 오히려 더 깨끗한 수렴 시험이다: # 경계조건·둘레 강성이 글자 하나 안 바뀌고 후보 자리 요소만 절반이 되므로 # 응력 변화가 '메시 때문' 임이 분명해진다(종전에는 둘레까지 같이 바뀌어 # 무엇 때문에 값이 변했는지 가릴 수 없었다). LOCAL_H_MULT = (2.0, 1.0) # ★ 고칠 것 2 [2026-10-11] — 세분화 비율을 '보장' 하는 방식으로 바꾼다. # 왜 바꾸는가(실측): 1단계 목표 19.45mm 를 걸었는데 실제는 8.4mm 가 나왔다. # tetgen 체적상한(-a)은 **상한**이라, 표면 삼각형과 품질조건이 더 빡빡하면 # 구속력이 전혀 없다. 그래서 2단계 목표를 9.73mm(목표의 절반)로 걸어도 이미 # 8.4mm 인 메시에는 아무 영향이 없어 실제 축소비가 1.09배밖에 안 됐다. # → 목표를 '1단계 목표의 절반' 이 아니라 **'1단계 실측 h 의 절반'** 으로 잡는다. # -a 는 상한이므로 실측보다 작은 목표는 반드시 구속력을 가진다. # ※ 숫자를 맞추려는 조정이 아니다 — 비율이 보장되지 않는 원인(상한 vs 실측)을 # 없애는 것이다. LOCAL_MAX_STAGES = 3 # 1 -> 1/2 -> 1/4 (절점 상한에 걸리면 거기서 멈춘다) LOCAL_NODE_CEIL = 950000 # 국부 모델 절점 상한 — 100만 절점 금지선을 넘지 않는다 # ★ 변을 1/2 로 줄이면 이상적으로는 절점이 8배다. 실측은 14.2~17.5배였다 # (후보1 77,768 -> 1,359,721 = 17.5배 / 후보2·3 46,944 -> 665,758 = 14.2배). # tetgen 의 품질조건이 체적상한 위로 요소를 더 넣기 때문이다. 예측에는 실측 # **하한**(8배의 1.8배 = 14.4배)을 쓴다 — 과하게 보수적으로 잡으면 들어갈 수 # 있는 단계를 포기해 버리고, 예측이 어긋나 상한을 넘은 단계는 판정에서 뺀다. LOCAL_REFINE_OVERSHOOT = 1.8 # 시도할 축소비 — 큰 쪽부터. 1.8 은 수렴 판정에 필요한 하한이다. LOCAL_RATIOS = (2.0, 1.9, 1.8) LOCAL_SURROUND_BUDGET = 600000 # 둘레 부재 절점 예산 — 두 단계 모두 같은 값(고정) CONV_GATE_PCT = 5.0 # 변화율이 이 이하면 '수렴' DIVERGE_PCT = 15.0 # 세분화할수록 응력이 이 이상 계속 커지면 '특이점 확정' def _bbox_dist(ctr, lo, hi): """점에서 경계상자까지의 거리(안에 있으면 0).""" return float(np.linalg.norm(np.maximum(np.maximum(lo - ctr, ctr - hi), 0.0))) def local_refine(cands, g, base_kw, top_n=3, budgets=None, logf=None): """2·3단계 — 상위 후보마다 국부 영역을 잘라 2단계 메시로 다시 풀고 수렴을 판정한다. 왜 '부품 단위' 로 자르는가: 부품 가운데를 잘라 그 단면에 변위를 강제하면 잘린 면 전체가 구속처럼 굳어 그 근처 응력이 가짜가 된다. 부품 경계(=용접 접합부)에서 자르면 영역 밖 구조가 **원래 연결돼 있던 그 자리(용접 씨앗점)** 로만 힘을 주고받으므로, 그 씨앗점에만 전체 해석의 변위를 강제하면 된다. 영역 안의 자유 표면은 그대로 자유다. """ def _log(m): if logf: logf(m) # 요소 크기 사다리 — 1단계만 미리 정한다(전체 해석 요소크기의 2배). # 2단계부터는 '이번 단계에서 후보 부재가 실제로 나온 h' 의 1/2 을 목표로 # 삼는다. 미리 정한 목표는 체적상한이 구속력을 못 가져 비율이 보장되지 않는다. hsteps = [float(g.get('h_global') or 10.0) * LOCAL_H_MULT[0]] # ★ 둘레 부재 예산은 두 단계 **같은 값**이다. 단계마다 바꾸면 둘레 메시가 # 달라져 경계조건이 변하고 수렴 판정이 또 무의미해진다(결함 2 와 같은 함정). budgets = budgets or [LOCAL_SURROUND_BUDGET] * len(hsteps) def _next_h(st, c): """다음 단계 목표 h = 이번 단계 '후보 부재 실측 h' 의 1/2. 반환 (목표 h 또는 None, 멈춘 사유). 절점 상한을 넘길 것 같으면 올라가지 않고 사유를 남긴다 — 100만 절점 금지.""" hs = [float(x) for x in (st.get('cand_h_mm') or []) if x] if not hs: return None, '후보 부재의 실측 요소크기를 읽지 못했습니다' cur = min(hs) # 후보 부재가 여러 개면 가장 가는 쪽 기준 cp = set(int(x) for x in (c.get('parts') or [])) cn = sum(int(x.get('nodes') or 0) for x in (st.get('per_part_h') or []) if int(x['no']) in cp) base = max(float(st.get('nodes') or 0) - cn, 0.0) # 둘레(안 깎는 쪽) # 들어가는 가장 큰 축소비를 고른다. 1.8 조차 안 들어가면 올라가지 않는다. for rr in LOCAL_RATIOS: proj = base + cn * (rr ** 3) * LOCAL_REFINE_OVERSHOOT if proj <= LOCAL_NODE_CEIL: return cur / rr, '' _p18 = base + cn * (LOCAL_RATIOS[-1] ** 3) * LOCAL_REFINE_OVERSHOOT return None, ('요소를 %.1f배 가늘게 하는 것조차 절점 약 %s개로 상한 %s개를 ' '넘습니다(후보 부재 %s절점 / 둘레 %s절점) - 여기서 멈춥니다' % (LOCAL_RATIOS[-1], '{:,}'.format(int(_p18)), '{:,}'.format(LOCAL_NODE_CEIL), '{:,}'.format(int(cn)), '{:,}'.format(int(base)))) out = [] parts = g['parts'] pno = [int(pp['no']) for pp in parts] seed_pairs = g.get('seed_pairs') or {} # ── 결함 4: 같은 자리를 두 번 계산하지 않는다 ───────────────────────── # ★ 반드시 top_n 으로 자르기 **전에** 합친다. 자른 뒤에 합치면 후보 1·2 가 # 하나로 줄어 3개가 아니라 2개 영역만 보게 된다. # 문턱은 전체 해석 요소크기의 3배 — 그 안이면 같은 요소 이웃이라 국부 영역 # 선택(부재 12개)도 같은 집합이 나온다(실측: 후보 1·2 가 9mm 간격, h=9.7mm). _hm = float(g.get('h_global') or 10.0) _mrg = max(3.0 * _hm, 1.0) uniq, merged = [], [] for c in cands: _p = np.asarray(c['at'], float) hit = None for ui, u in enumerate(uniq): if float(np.linalg.norm(_p - np.asarray(u['at'], float))) <= _mrg: hit = ui break if hit is None: uniq.append(c) merged.append([]) else: merged[hit].append({'at': c['at'], 'kind': c.get('kind'), 'dist_mm': round(float(np.linalg.norm( _p - np.asarray(uniq[hit]['at'], float))), 2)}) _log('후보 중복 제거: %d개 → %d개 (문턱 %.1fmm = 요소 %.1fmm x3)' % (len(cands), len(uniq), _mrg, _hm)) for ci, c in enumerate(uniq[:top_n]): ctr = np.asarray(c['at'], float) d = np.array([_bbox_dist(ctr, pp['lo'], pp['hi']) for pp in parts]) order = np.argsort(d)[:LOCAL_MAX_PARTS] keep = sorted(int(k) for k in order) keep_no = set(pno[k] for k in keep) R = float(d[order].max()) # ★ 영역 부재 수를 '메시 전 추정' 으로 줄이려 했으나 **되돌렸다** # (실측 2026-10-11). 추정식을 세울 수가 없다 — 같은 영역 12부재에서 # 표면삼각형 기준 추정은 실측의 0.17~1.44배, 부피·벽두께 기준 추정은 # 0.17~156배로 틀렸다. 표면이 성긴 부재(53번 24삼각형 -> 28,819절점)는 # 품질조건(-q)이 절점을 결정하고, 촘촘한 부재는 '-p' 가 결정하기 때문이다. # 맞지도 않는 식으로 영역을 깎으면 '원인 없이 숫자 바꾸기' 가 된다. # -> 영역은 종전대로 가까운 12부재로 두고, 2단계가 95만 절점에 들어가는지는 # **1단계 실측 절점 수**로만 판정한다(_next_h). # ★ 결함 3: 씨앗점마다 '어느 접합면에서 온 점인지' 를 들고 간다. 경계 힘을 # 알짜 합으로 비교하면 평형 때문에 언제나 맞아 검증이 공허하다 — # 접합면 하나하나의 힘을 비교해야 경계조건의 타당성을 잴 수 있다. dp, dsrc = [], [] # ★ 재현성: 접합면을 보는 순서가 drive/dsrc 순서를 정하고, 그것이 # '씨앗점마다 최근접 절점 1개(충돌이면 두 번째)' 배정을 바꾼다. # 부품번호 쌍으로 정렬해 순서를 고정한다. for (pa, pb), v in sorted(seed_pairs.items(), key=lambda kv: (int(kv[0][0]), int(kv[0][1]))): if (pa in keep_no) != (pb in keep_no) and len(v): dp.append(v) dsrc.extend([(int(pa), int(pb))] * len(v)) drive = np.vstack(dp) if dp else np.zeros((0, 3)) sep = (float(np.linalg.norm(drive - ctr, axis=1).min()) if len(drive) else None) # 전체 해석에서 '잘린 접합면' 들을 지나던 힘의 합(벡터) — 국부 해석의 # 경계 강제변위 절점 절점력 합과 맞아야 경계 조건이 타당하다. # ★ 1단계 (2026-10-12): iface_rf 의 값이 '힘 벡터' 에서 # {'F','side_part_no','n','kind'} 로 바뀌었다. 옛 캐시(벡터)도 읽을 수 # 있게 둘 다 받는다 — 캐시를 다시 풀라고 사람을 시키지 않는다. _irf = g.get('iface_rf') or {} _cutd, _cutside, _cutkind = {}, {}, {} for (pa, pb), v in _irf.items(): if (pa in keep_no) == (pb in keep_no): continue _kk = (int(pa), int(pb)) if isinstance(v, dict): _cutd[_kk] = np.asarray(v.get('F'), float) _cutside[_kk] = v.get('side_part_no') _cutkind[_kk] = v.get('kind') else: _cutd[_kk] = np.asarray(v, float) _cutside[_kk] = None _cutkind[_kk] = '접합면' _cut = list(_cutd.values()) gcut = (np.sum(_cut, axis=0) if _cut else None) soup, sizes, names = [], [], [] for k in keep: T = parts[k]['cadV'][parts[k]['cadF']].reshape(-1, 3) soup.append(T) sizes.append(int(T.size)) names.append(parts[k]['name'] or str(parts[k]['no'])) Pn = np.vstack(soup) lo = Pn.min(axis=0) - 1.0 hi = Pn.max(axis=0) + 1.0 lps = [lp for lp in (base_kw.get('load_points') or []) if np.all(np.asarray(lp['p'], float) >= lo) and np.all(np.asarray(lp['p'], float) <= hi)] rec = {'cand': ci + 1, 'at': c['at'], 'kind': c['kind'], 'parts_n': len(keep), 'R_mm': round(R, 1), 'drive_pts': int(len(drive)), 'cand_to_drive_mm': (round(sep, 1) if sep is not None else None), 'load_points_inside': len(lps), 'cut_ifaces': len(_cut), 'merged_n': len(merged[ci]), 'merged': merged[ci], 'merge_tol_mm': round(_mrg, 1), 'global_cut_N': ([round(float(x), 2) for x in gcut] if gcut is not None else None), 'stage': []} if not len(drive): rec['verdict'] = '판정 불가' rec['why'] = '영역 밖과 잇던 용접 씨앗점이 없습니다 - 경계 변위를 줄 수 없습니다' out.append(rec) continue # ★ 결함 3: 영역 안에 들어온 **실제 지지점은 살려 둔다**. 종전에는 전부 # 버렸다. 후보 3 의 영역에는 이 기계의 지지점 6개가 전부 들어 있어서, # 지지를 버리면 국부 모델이 강제변위만으로 떠 있는 전혀 다른 구조가 된다. # ※ 종전에 여기 적었던 '경계 힘 어긋남 1,764% 의 정체가 이것' 은 **취소한다**. # 통제실험(2026-10-11)에서 그 어긋남은 강제 절점 수에 딸린 값임이 드러났다 # (207->58 절점에서 20,525N->24N, 형상은 그대로). 지지점을 살리는 것은 # 물리적으로 옳아서 유지하는 것이지, 그 수치의 원인이어서가 아니다. fps = [fp for fp in (base_kw.get('fix_points') or []) if np.all(np.asarray(fp['p'], float) >= lo) and np.all(np.asarray(fp['p'], float) <= hi)] rec['fix_points_inside'] = len(fps) kw = dict(base_kw) kw.pop('load_points', None) kw.pop('fix_points', None) rad = float(c.get('radius_mm') or 0.0) # ★ 후보 자리 응력을 재는 반경은 '전체 해석의 요소크기' 로 후보마다 한 번만 # 정한다. 단계마다 그 단계의 요소크기로 정하면 두 단계가 서로 다른 영역을 # 재게 되어 변화율이 뜻을 잃는다. selr = max(rad + 3.0 * float(g.get('h_global') or 10.0), 20.0) rec['cand_sel_r_mm'] = round(selr, 1) hsteps = list(hsteps[:1]) # 후보마다 1단계부터 다시 쌓는다 si = -1 while True: si += 1 if si >= len(hsteps) or si >= LOCAL_MAX_STAGES: break _ht = hsteps[si] # ★ 둘레 부재 예산은 모든 단계에서 **같은 값**이다(결함 2 와 같은 함정). # 세분화하는 것은 후보 부재뿐이고 둘레 메시는 단계마다 동일해야 # 경계조건과 둘레 강성이 변하지 않는다. # ★ run_fea 의 하드 벽은 budget_nodes x 1.3 이다. 국부 모델의 상한은 # LOCAL_NODE_CEIL(95만) 이므로 그 벽이 95만에 오도록 역산해 넘긴다. # 종전 600,000 은 벽이 78만이어서, 후보 부재를 절반으로 줄인 단계가 # '부품이 너무 많습니다' 로 통째로 실패했다(실측 87만·85만 절점). bud = LOCAL_NODE_CEIL / 1.3 t0 = time.time() try: r = run_fea(Pn.reshape(-1).tolist(), part_sizes=sizes, part_names=names, load_points=(lps or None), fix_points=(fps or None), budget_nodes=bud, quality_bulk=False, weld_seed=True, keep_arrays=True, # ★ 결함 1 고침: 이걸 안 넘기면 국부 해석의 부품번호가 # 1..12 로 새로 붙어, 후보 부재의 '전체 해석 부품번호' # (예: 51번)와 **절대 일치할 수 없다**. 그래서 # cand_parts_in_mesh 가 항상 비고 '형상 처리 실패'라는 # 틀린 사유로 전 후보가 판정 불가가 됐다(실측 확인: # 후보 부재 51·28번은 두 단계 모두 메시에 성공한다). _orig_nos=[int(parts[k]['no']) for k in keep], submodel={'nodes': g['nodes'], 'elems': g['elems'], 'disp': g['disp'], 'drive_pts': drive, 'drive_src': dsrc, 'drive_r': LOCAL_DRIVE_R_MM, 'h_target': _ht, # ★ 둘레 부재에는 h_target 을 걸지 않는다(실측 # 2026-10-11에 걸어 보고 되돌렸다): tetgen 은 # 체적상한 -a 를 주면 품질조건을 기본 적용해서 # '성기게' 가 아니라 거꾸로 폭발한다 — # 둘레가 305,806 -> 761,728절점, 요소 8.4 -> 1.2mm. # 예산은 둘레를 깎는 대신 **영역 부재 수**로 조절한다 # (아래 영역 선택 참조). 'h_parts': [int(x) for x in (c.get('parts') or [])]}, **kw) except Exception as ex: rec['stage'].append({'n': si + 1, 'h_target_mm': round(_ht, 2), 'error': '%s: %s' % (type(ex).__name__, ex), 'sec': round(time.time() - t0, 1)}) _log(' cand%d stage%d FAIL: %s' % (ci + 1, si + 1, ex)) continue A = r['_arrays'] _h = r['mesh'].get('elem_size_real_mm') or 10.0 _dn = np.linalg.norm(A['nodes'] - ctr, axis=1) # ★ 반경을 넓혀 가며 찾지 않는다. 넓히면 '후보에서 1m 떨어진 어딘가의 # 최대응력' 을 후보 응력이라고 적게 된다 — 계산 안 한 것을 표시하는 짓이다. okn = (_dn <= selr) & (A['excl'] == 0) # ★★ 후보 자리 응력은 **후보가 있는 부재 위에서만** 재야 한다. # 실측(2026-10-09): 반경만으로 고르면 둘레의 세분화 안 된 부재가 # 최대값을 차지해, 후보 부재를 절반으로 깎아도 cand_MPa 가 소수 # 둘째자리까지 똑같이 나왔다(157.82 → 157.82) → '수렴' 이라는 # 거짓 판정이 나온다. 안 깎은 부재의 값으로 수렴을 말하는 것은 # 계산 안 한 것을 표시하는 짓이다. _cp = set(int(x) for x in (c.get('parts') or [])) if _cp: _pm = np.zeros(len(A['nodes']), dtype=bool) for _pp in (A['parts'] or []): if int(_pp['no']) in _cp: _pm[int(_pp['n0']):int(_pp['n1'])] = True if _pm.any(): okn = okn & _pm st = { 'n': si + 1, 'h_target_mm': round(_ht, 2), 'h_miss_parts': r['mesh'].get('h_miss_parts'), 'per_part_h': r['mesh'].get('per_part_h'), 'nodes': r['mesh']['nodes'], 'elems': r['mesh']['elements'], 'parts': r['mesh']['parts'], 'elem_mm': _h, 'welds_tied': r['mesh'].get('welds'), 'seeds': (r['mesh'].get('weld_seed') or {}).get('seeds'), 'disp_mm': r['result']['max_disp_mm'], 'struct_MPa': r['result']['max_vonmises_MPa'], 'cand_MPa': (round(float(A['vm'][okn].max()), 2) if okn.any() else None), 'cand_nodes': int(okn.sum()), 'cand_parts_in_mesh': sorted(set( int(x) for x in (c.get('parts') or []) if x in set(int(pp['no']) for pp in (A['parts'] or [])))), 'parts_meshed': r['mesh']['parts'], 'cand_h_mm': ([x['h_mm'] for x in (r['mesh'].get('per_part_h') or []) if int(x['no']) in set(int(y) for y in (c.get('parts') or []))] or None), 'compliance_Nmm': r.get('compliance_Nmm'), 'driven': (r.get('submodel') or {}).get('driven_nodes'), 'interp_in_tet': (r.get('submodel') or {}).get('interp_in_tet'), 'interp_clamped': (r.get('submodel') or {}).get('interp_clamped'), 'interp_nearest': (r.get('submodel') or {}).get('interp_nearest'), 'rf_driven_N': r.get('rf_driven_N'), 'rf_by_iface_N': r.get('rf_driven_by_iface_N'), 'rf_by_iface_side': r.get('rf_driven_by_iface_side'), 'sec': round(time.time() - t0, 1), 'solve_sec': r['solve_sec'], } # 경계 단면력 일치율 — 전체 해석의 '잘린 접합면 힘 합' 대비 몇 % 어긋났나. # 두 값이 모두 작으면 비율이 뜻을 잃으므로 크기도 같이 적는다. if st['rf_driven_N'] and gcut is not None: _lv = np.asarray(st['rf_driven_N'], float) _gm = float(np.linalg.norm(gcut)) st['cut_force_global_N'] = round(_gm, 2) st['cut_force_local_N'] = round(float(np.linalg.norm(_lv)), 2) st['cut_force_err_pct'] = (round(float(np.linalg.norm(_lv - gcut)) / _gm * 100.0, 2) if _gm > 1e-6 else None) # ★ 결함 3 의 본체 — 접합면별 비교. # 알짜 합(위 cut_force_err_pct)은 솔버가 평형을 강제하므로 경계조건의 # 좋고 나쁨을 재지 못한다(계산이 틀려도 총합은 맞는다). 접합면 # 하나하나의 힘은 평형이 정해 주지 않으므로 이것이 진짜 검증이다. _lbi = st.get('rf_by_iface_N') or {} if _lbi and _cutd: _rows, _errs = [], [] for _k, _gv in sorted(_cutd.items()): _lvv = _lbi.get('%d-%d' % _k) or _lbi.get('%d-%d' % (_k[1], _k[0])) if _lvv is None: _rows.append({'iface': '%d-%d' % _k, 'global_N': round(float(np.linalg.norm(_gv)), 1), 'local_N': None, 'err_pct': None, 'note': '이 접합면 씨앗이 국부 메시에 없습니다'}) continue _lvv = np.asarray(_lvv, float) # ★★ 1단계 고침 (2026-10-12) — 부호를 맞춘다. # 전체 해석의 힘은 '종속절점이 있는 부품' 쪽 절점력 합이고, # 국부의 힘은 '영역 안 부품' 쪽 경계 반력 합이다. 두 쪽이 # 서로 반대 부품이면 같은 힘이 부호만 뒤집혀 나오고, 그러면 # 어긋남이 정의상 200% 로 찍힌다(물리적 오차가 아니다). # → 어느 쪽인지 기록해 둔 side_part_no 로 판정해 뒤집는다. _gs = _cutside.get(_k) _ls = (st.get('rf_by_iface_side') or {}).get('%d-%d' % _k) if _ls is None: _ls = (st.get('rf_by_iface_side') or {}).get('%d-%d' % (_k[1], _k[0])) _flip = (_gs is not None and _ls is not None and int(_gs) != int(_ls)) if _flip: _lvv = -_lvv _gm2 = float(np.linalg.norm(_gv)) # ★ 힘이 거의 0 인 접합면의 비율은 뜻이 없다(실측 8.9e11%). # 가장 큰 접합면의 1% 미만이면 비율을 적지 않는다 — # 0 으로 나눈 숫자를 '어긋남' 으로 보고하면 거짓 경보다. _gmax = max(float(np.linalg.norm(v)) for v in _cutd.values()) _ep = (float(np.linalg.norm(_lvv - _gv)) / _gm2 * 100.0 if _gm2 > max(1e-6, 0.01 * _gmax) else None) _rows.append({'iface': '%d-%d' % _k, 'global_N': round(_gm2, 1), 'local_N': round(float(np.linalg.norm(_lvv)), 1), 'err_pct': (round(_ep, 1) if _ep is not None else None), 'kind': _cutkind.get(_k), 'side_global': _gs, 'side_local': _ls, 'flipped': bool(_flip), 'note': (None if (_gs is not None and _ls is not None) else '종속 쪽을 몰라 부호를 못 맞췄습니다')}) if _ep is not None: _errs.append(_ep) # ★ 1단계: '기준값이 없는 쌍' 을 표에서 숨기지 않는다. # 국부가 강제한 접합면인데 전체 해석에 그 접합면 힘이 없으면, # 그 자리는 **맞는지 틀리는지 확인할 방법이 없는 경계**다. # 비율에는 넣지 않고(기준이 없으니) 표에 그대로 남긴다. _gkeys = set('%d-%d' % _k for _k in _cutd) | set('%d-%d' % (_k[1], _k[0]) for _k in _cutd) for _lk, _lvv0 in sorted(_lbi.items()): if _lk in _gkeys: continue _rows.append({'iface': _lk, 'global_N': None, 'local_N': round(float(np.linalg.norm( np.asarray(_lvv0, float))), 1), 'err_pct': None, 'kind': None, 'side_global': None, 'side_local': (st.get('rf_by_iface_side') or {}).get(_lk), 'flipped': False, 'note': '전체 해석에 이 접합면 힘이 없습니다 — ' '기준값 없이 국부가 강제한 경계입니다'}) st['iface_force'] = _rows st['iface_force_no_ref'] = int(sum(1 for x in _rows if x['global_N'] is None)) # 힘이 큰 접합면에 가중을 준다 — 거의 0N 인 접합면의 비율은 뜻이 없다. _wg = np.array([x['global_N'] for x in _rows if x['err_pct'] is not None]) _we = np.array([x['err_pct'] for x in _rows if x['err_pct'] is not None]) if len(_we): st['iface_force_err_wpct'] = round(float((_wg * _we).sum() / max(_wg.sum(), 1e-9)), 2) st['iface_force_err_max_pct'] = round(float(_we.max()), 2) st['iface_force_n'] = int(len(_we)) st['iface_force_missing'] = int(sum(1 for x in _rows if x['err_pct'] is None)) rec['stage'].append(st) _log(' cand%d stage%d: nodes %s / h %.1fmm / disp %.4f / candMPa %s / %.0fs' % (ci + 1, si + 1, '{:,}'.format(st['nodes']), _h, st['disp_mm'], st['cand_MPa'], st['sec'])) _log(' h목표 %.1fmm / 못맞춘부재 %s / 접합면힘 가중오차 %s%% (최대 %s%%)' % (_ht, st.get('h_miss_parts'), st.get('iface_force_err_wpct'), st.get('iface_force_err_max_pct'))) # ── 다음 단계 목표 = 이번 단계 후보 부재 실측 h 의 1/2 ────────── if si + 1 < LOCAL_MAX_STAGES: _nh, _stopwhy = _next_h(st, c) if _nh: hsteps.append(_nh) _log(' 다음 단계 목표 %.2fmm (이번 실측 %.2fmm 의 1/2)' % (_nh, _nh * 2.0)) else: rec['ladder_stop'] = _stopwhy _log(' 사다리 멈춤: %s' % _stopwhy) del r, A # ★ 100만 절점 금지선을 넘긴 단계는 판정 근거로 쓰지 않는다. # h_target 경로는 budget_nodes 를 무시하므로(요소 크기를 직접 걸기 때문) # _next_h 의 사전 예측이 유일한 방어선이다. 예측이 어긋나 상한을 넘은 # 단계가 생기면 그 단계 위에 판정을 올리지 않는다. for _x in rec['stage']: if 'error' not in _x and int(_x.get('nodes') or 0) > LOCAL_NODE_CEIL: _x['over_ceiling'] = True ok = [x for x in rec['stage'] if 'error' not in x and not x.get('over_ceiling')] _nocand = [x for x in ok if not x.get('cand_parts_in_mesh')] if len(ok) < 2: rec['verdict'] = '판정 불가' rec['why'] = ('판정에 쓸 수 있는 단계가 %d개뿐입니다' % len(ok) + (' — ' + rec['ladder_stop'] if rec.get('ladder_stop') else '') + (' (절점 상한을 넘겨 제외한 단계 %d개)' % sum(1 for x in rec['stage'] if x.get('over_ceiling')) if any(x.get('over_ceiling') for x in rec['stage']) else '')) elif _nocand: rec['verdict'] = '판정 불가' rec['why'] = ('후보 부재 %s번이 국부 메시에 없습니다(형상 처리 실패) — ' '이 자리 응력을 잴 수 없어 수렴을 판정하지 않습니다' % ', '.join(str(x) for x in (c.get('parts') or []))) else: a, b = ok[-2], ok[-1] def ch(k): va, vb = a.get(k), b.get(k) if not va or vb is None: return None return abs(vb - va) / abs(va) * 100.0 cd, cs, cc = ch('disp_mm'), ch('cand_MPa'), ch('compliance_Nmm') # ★ 최대변위는 판정에 쓰지 않는다. 서브모델에서는 최대변위 절점이 # 거의 항상 '경계 강제변위 절점' 이라, 전체 해석에서 그대로 받아 쓴 # 값이 다시 나온다(실측: 두 단계가 소수 4자리까지 같다). 그것을 # '변위가 0% 변했으니 수렴' 의 근거로 적으면 계산 안 한 것을 # 표시하는 짓이다. 참고로만 남기고 게이트는 응력·컴플라이언스로 본다. vals = [x for x in (cs, cc) if x is not None] conv = max(vals) if vals else None # ★ 실측 축소비 — '절반으로 줄였다' 고 쓰지 않고 실제로 몇 배인지 적는다. _ha0 = min([float(x) for x in (a.get('cand_h_mm') or []) if x] or [0]) or None _hb0 = min([float(x) for x in (b.get('cand_h_mm') or []) if x] or [0]) or None rec['h_ratio'] = (round(_ha0 / _hb0, 2) if (_ha0 and _hb0) else None) rec['h_pair_mm'] = [_ha0, _hb0] rec['stages_done'] = len(ok) rec['conv'] = {'disp_pct_ref_only': (round(cd, 2) if cd is not None else None), 'disp_note': '경계 강제값 — 판정에 쓰지 않음', 'cand_MPa_pct': (round(cs, 2) if cs is not None else None), 'compliance_pct': (round(cc, 2) if cc is not None else None), 'worst_pct': (round(conv, 2) if conv is not None else None)} grow = (b.get('cand_MPa') or 0) - (a.get('cand_MPa') or 0) if cs is None: rec['verdict'] = '판정 불가' rec['why'] = '후보 자리 응력을 두 단계 모두에서 재지 못했습니다' elif conv is not None and conv <= CONV_GATE_PCT: rec['verdict'] = '구조 후보(수렴)' if c['kind'] == '구조 후보' else '수렴' rec['why'] = ('후보자리 응력·컴플라이언스 변화율이 모두 %.1f%% 이하입니다' ' (최대변위는 경계 강제값이라 판정에 쓰지 않음)' % CONV_GATE_PCT) elif cs is not None and grow > 0 and cs > DIVERGE_PCT: rec['verdict'] = '특이점 확정' # ★ '절반' 이라고 쓰지 않는다. 목표는 2배였지만 벽두께가 더 빡빡해 # 실제로는 1.3~1.8배만 가늘어진다 — 실측한 그 비율을 그대로 적는다. _ha = _ha0 _hb = _hb0 _hw = ('요소를 %.1fmm → %.1fmm (%.2f배) 로 줄였는데' % (_ha, _hb, _ha / _hb) if (_ha and _hb and _hb > 0) else '요소를 더 가늘게 했는데') rec['why'] = ('%s 이 자리 응력이 %.1f%% 더 커졌습니다 ' '- 값이 발산하는 특이점입니다(실제 위험이 아님)' % (_hw, cs)) else: rec['verdict'] = '미수렴(판정 유보)' rec['why'] = ('변화율 %.1f%% - %.1f%% 게이트를 못 넘겼습니다' % (conv or 0, CONV_GATE_PCT)) out.append(rec) return out def run_fea(positions, load_N=1000.0, material=DEF_MAT, E=None, nu=None, fix_points=None, load_points=None, pick_radius=None, keep=False, force_obb=False, budget_nodes=None, weld_mode='bond', weld_throat=0.0, part_sizes=None, plastic=False, part_names=None, _drop_rounds=2, _auto_dropped=None, _orig_nos=None, quality_bulk=False, weld_seed=True, submodel=None, keep_arrays=False, local_refine_n=0, _logf=None, _fine_parts=None, _fine_rounds=1): """positions: 평탄한 [x,y,z,...] 삼각형 수프. fix_points : [{'p':[x,y,z], 'r':반경}, ...] (없으면 자동 외팔보) load_points : [{'p':[x,y,z], 'r':반경, 'dir':[x,y,z], 'N':크기}, ...] part_names : part_sizes 와 같은 순서의 부품 이름. 실패한 부품을 이름으로 알려주기 위함. """ mat = MATERIALS.get(material, MATERIALS[DEF_MAT]) E = float(E) if E else mat['E'] nu = float(nu) if nu else mat['nu'] fy = mat['fy'] P = np.asarray(positions, dtype=float).reshape(-1, 3) if len(P) < 4: raise ValueError('정점이 너무 적습니다(%d개)' % len(P)) # 조건을 찍기 시작했는데 하중점이 없으면, 메시를 만들기 전에 바로 알려준다 # (메시부터 만들면 수십 초 버리고 실패한다) if fix_points and not load_points and not submodel: raise ValueError('고정점만 찍혔습니다 — 하중점도 찍어주세요.') # ★ submodel(국부 재해석)에서는 하중이 없어도 된다 — 경계 강제변위가 # 하중 역할을 한다. 영역 안에 들어온 실제 지지점은 살려야 하므로(결함 3) # 이 안내문이 국부 해석을 막으면 안 된다. # ★ 자동 제외 재귀 때문에 부품 자리번호가 밀린다. 화면에 보여줄 번호는 반드시 # '원래 블럭에서의 부품번호' 여야 하므로 그 대응표를 들고 다닌다. def _onum(k): try: return int(_orig_nos[k]) if _orig_nos else k + 1 except Exception: return k + 1 job = uuid.uuid4().hex[:10] wd = os.path.join(WORK_ROOT, job) os.makedirs(wd, exist_ok=True) center, R, size = obb_from_points(P) method, etype, warn = 'tet_c3d10', 'C3D10', None meshinfo = {} parts, weld_list, tie_blocks = [], [], [] _seed_stat = {} # 1단계 용접 씨앗점 통계 _seed_pairs_out = {} # (부품번호쌍) -> 씨앗점 좌표 — 국부 재해석이 그대로 쓴다 _sjarr = None # 요소별 스케일드 야코비안 — 후보가 '메시 결함' 인지 가리는 데 쓴다 failed_parts = [] # [{'no':1-base 부품번호, 'name':이름, 'why':한국어 사유}] if not force_obb: # ★ 부품 경계를 먼저 나눈 뒤, 부품마다 따로 정점을 병합한다. # 전체를 한꺼번에 병합하면 '맞닿은 부재끼리 정점이 붙어' 한 덩어리가 되고, # 그 덩어리는 속에 내부면이 낀 자기교차 형상이라 tetgen 이 거부한다 # (실측: 부재 5개 투입 → 성분 2개로 뭉침 → "A self-intersection was detected"). # ★ comps 는 (정점, 삼각형, 원래 부품번호) 로 들고 다닌다. # 원래 부품번호를 잃어버리면 '부재 7번이 공중에 떴다' 같은 경고가 # 실패한 부품만큼 밀려서 엉뚱한 부품을 가리킨다. if part_sizes: comps = [] off = 0 for pi, cnt in enumerate(part_sizes): rows = int(cnt) // 3 # part_sizes 는 '좌표값 개수'(정점수x3)로 받는다 seg = P[off:off + rows] off += rows if len(seg) < 4: failed_parts.append({'no': _onum(pi), 'name': _pname(part_names, pi), 'why': '정점이 %d개뿐 — 형상이 없습니다' % len(seg)}) continue Vp, Fp = weld_mesh(seg) sub = split_components(Vp, Fp) if not sub: failed_parts.append({'no': _onum(pi), 'name': _pname(part_names, pi), 'why': '삼각형이 너무 적어 덩어리를 못 만듭니다'}) continue comps.extend([(a, b, pi) for (a, b) in sub]) V, F = weld_mesh(P) # 통계 표시용 else: V, F = weld_mesh(P) comps = [(a, b, ci) for ci, (a, b) in enumerate(split_components(V, F))] if not comps: raise RuntimeError('삼각형 덩어리를 찾지 못했습니다') # 여러 부재면 절점 예산을 나눠 쓴다 _tot_budget = float(budget_nodes or MAX_C3D10_NODES) _mets = [mesh_metrics(vc, fc) for vc, fc, _pi in comps] _vols = [max(m['volume_mm3'], 1.0) for m in _mets] _vsum = sum(_vols) # 부품이 많을수록 부품당 품질개선 시간을 짧게 — 전체가 늘어지지 않게 # 부재 1개 목표는 10초 이내 → 품질개선에 8초까지만 준다(검증용 각파이프는 2.8초라 여유). # ★ 계획 F-8 수정: 종전 공식 `max(2.5, 40.0/len(comps))` 은 부품당 시간을 # 부품 수에 반비례시켰다. 실측으로 이것이 '연결성 버그' 로 보였던 현상의 # 진짜 원인이었다 — 같은 14부품이 2.86초에서는 8개가 메시 실패로 조용히 # 빠지고, 240초에서는 14/14 전부 성공했다. 국부 재해석의 품질이 영역 크기에 # 따라 달라지면 수렴 판정이 무의미해지므로 부품 수와 분리한다. _qto = QUALITY_TIMEOUT_S # ★ 전체해석 모드 — 부품이 MANY_PARTS(40) 를 넘으면 # ① 수리 사다리로 열린 형상을 고치고 ② 품질조건(-q) 승격을 건너뛰고 # ③ 부품당 절점 상한(PER_CAP)을 건다. # 부품 1~40개(검증·일상 작업)는 기존 경로를 그대로 쓴다 — selftest 이론값 대조가 # 성긴 메시로는 통과하지 못하기 때문이다. _bulk = len(comps) > MANY_PARTS _repair_stat = {} _hull_parts = [] # ★ 고칠 것 3 — '두께 0 면 형상' 을 판두께 추정 압출로 살린다. # 반드시 repair_for_tet **전에** 한다: repair_for_tet 은 입력 V 의 # bbox 최소변(emin)으로 두께 게이트 상한(emin*1.2+1)을 잡는데, # 납작한 판은 emin=0 이라 상한이 2.2mm 가 되어 어떤 두께의 압출도 # 거부된다. 압출한 형상을 '원본' 으로 넣어야 게이트가 맞는다. _sheet_est = [] # [{'no','name','t_mm','basis'}] — 화면 라벨용 _nb = [] for _k, (_vk, _fk, _pk) in enumerate(comps): _vk = np.asarray(_vk, float) _nb.append(((_vk.min(axis=0), _vk.max(axis=0)), float(_mets[_k].get('wall_mm') or 0.0))) comps2 = [] for _k, (_vk, _fk, _pk) in enumerate(comps): _isflat, _ax = is_flat_sheet(_vk) if not _isflat: comps2.append((_vk, _fk, _pk)) continue _t, _basis = guess_sheet_t(_nb[_k][0], [_nb[j] for j in range(len(_nb)) if j != _k]) _ok = None if _t: _ok = extrude_flat_sheet(_vk, _fk, _t, _ax) if _ok is None: comps2.append((_vk, _fk, _pk)) # 못 살렸다 — 종전대로 실패 사유를 남긴다 continue comps2.append((_ok[0], _ok[1], _pk)) _sheet_est.append({'no': _onum(_pk), 'name': _pname(part_names, _pk), 't_mm': round(float(_t), 2), 'basis': _basis}) if _logf: _logf(' 두께 0 면 형상 %d번(%s) → 두께 추정 %.2fmm 압출 (%s)' % (_onum(_pk), _pname(part_names, _pk), _t, _basis)) if _sheet_est: comps = comps2 _mets = [mesh_metrics(vc, fc) for vc, fc, _pi in comps] _vols = [max(m['volume_mm3'], 1.0) for m in _mets] _vsum = sum(_vols) _hstat = [] # 국부 재해석: 부재별 목표/실제 요소크기(단조 세분화 확인용) _fine_stat = [] # 밀도 상향 부재: 요구/실제 요소크기와 절점 수(실측 표기용) nlist, elist, n0 = [], [], 0 walls, esizes, coarse, flips, qual = [], [], False, 0, '' for ci, (Vc, Fc, _pi) in enumerate(comps): sub = os.path.join(wd, 'c%d' % ci) os.makedirs(sub, exist_ok=True) per = max(_tot_budget * _vols[ci] / _vsum, PER_NODE_FLOOR) if _bulk: per = min(per, PER_CAP) # ★★ 여기가 '부품별 격리' 의 핵심. # 이 try 가 없으면 부품 1개가 깨질 때 전체가 외곽상자로 물러났다. # 깨진 부품은 사유를 적어 두고 건너뛰고, 나머지로 계속 간다. try: if _bulk: # ★ 품질 경로는 선택식이다(계획 부록 G). 기본값은 종전 `-p` 전용 — # 품질 경로는 절점 89만·메모리 약 20GB·260초가 들어서 # 운전자가 모르게 켜면 기계가 죽을 수 있다. # ★ 밀도 상향 부재(_fine_parts)는 이 부재만 요소크기를 지정해 # 다시 메시한다. 다른 부재는 한 톨도 안 바뀐다. _fh = (_fine_parts or {}).get(_onum(_pi)) if _fh and quality_bulk: tm = build_tet_mesh_fine(Vc, Fc, sub, float(_fh), logf=_logf) if tm is None: raise RuntimeError('밀도 상향 메시 실패(요구 %.2fmm)' % float(_fh)) _fine_stat.append({'no': _onum(_pi), 'h_req_mm': round(float(_fh), 3), 'h_mm': round(float(tm['elem_size_mm']), 2), 'nodes': int(len(tm['nodes'])), 'capped': bool(tm.get('fine_capped'))}) else: tm = (build_tet_mesh_quality_bulk(Vc, Fc, sub) if quality_bulk else build_tet_mesh_repaired(Vc, Fc, sub, budget_nodes=per)) _repair_stat[tm['repair']] = _repair_stat.get(tm['repair'], 0) + 1 if tm.get('hull'): _hull_parts.append(_onum(_pi)) elif submodel: # 국부 재해석 — 수리 + 품질 경로(위 build_tet_mesh_local 주석 참조) tm = build_tet_mesh_local(Vc, Fc, sub, budget_nodes=per, quality_timeout=_qto, h_target=( submodel.get('h_target') if (not submodel.get('h_parts') or _onum(_pi) in set( submodel.get('h_parts'))) else None)) _hstat.append({'no': _onum(_pi), 'h_mm': round(float(tm['elem_size_mm']), 2), 'h_target_mm': tm.get('h_target_mm'), 'h_hit': tm.get('h_hit'), 'nodes': int(len(tm['nodes']))}) _repair_stat[tm['repair']] = _repair_stat.get(tm['repair'], 0) + 1 if tm.get('hull'): _hull_parts.append(_onum(_pi)) else: tm = build_tet_mesh(Vc, Fc, sub, budget_nodes=per, quality_timeout=_qto) except Exception as _ex: failed_parts.append({'no': _onum(_pi), 'name': _pname(part_names, _pi), 'why': _why_mesh_failed(_ex, _mets[ci], Vc)}) continue nd, el = tm['nodes'], tm['elems'] f6, fopp = boundary_faces6(el) f6 = f6 + n0 fopp = fopp + n0 tri = f6[:, :3] - n0 _g = f6[:, :3] - n0 _cr = np.cross(nd[_g[:, 1]] - nd[_g[:, 0]], nd[_g[:, 2]] - nd[_g[:, 0]]) _ln = np.maximum(np.linalg.norm(_cr, axis=1)[:, None], 1e-12) _nn = _cr / _ln _ctr = nd[_g].mean(axis=1) _fl = np.einsum('ij,ij->i', _nn, _ctr - nd[fopp - n0]) < 0 _nn[_fl] *= -1.0 sids = np.unique(tri) + n0 parts.append({ 'nodes': nd, 'elems': el, 'n0': n0, 'e0': len(elist) and sum(len(x) for x in elist) or 0, 'faces6': f6, 'faces_opp': fopp, 'surf_ids': sids, # ★★ 2a (2026-10-12): 씨앗을 받을 후보 절점 풀을 **중간절점까지** # 넓힌 것. surf_ids 는 꼭짓점(코너)만이라, 성긴 부재에서는 씨앗에서 # 가장 가까운 코너도 수십 mm 떨어져 있고 그 절점은 상대 표면에서 # 묶을 거리 한도(wtol x1.5 = 4.5mm)를 넘는다 → 손실 채널 # too_far 668개(전체 손실 883의 75.6%, 0단계 실측)의 정체다. # C3D10 은 변마다 중간절점이 있으므로 풀을 넓히면 같은 메시에서 # 씨앗까지의 거리가 절반으로 줄어든다. 메시를 바꾸지 않는다. # ※ 기존 surf_ids / surf_pts 키는 그대로 둔다 — 면 선별·떠있음 # 점검·진단 덤프가 그 키를 쓰고 있어 건드리면 딴 게 바뀐다. 'surf_ids_all': np.unique(f6), 'surf_pts_all': nd[np.unique(f6) - n0], 'faces_n': _nn, 'faces_ctr': _ctr, 'n1': n0 + len(nd), 'surf_pts': nd[np.unique(tri)], 'wall': tm['metrics']['wall_mm'], # ★ 원래 부품번호/이름을 반드시 들고 간다. 실패한 부품을 건너뛰면 # parts 안의 자리번호가 밀리므로, 화면에 보여줄 번호는 이것만 쓴다. 'part_no': _onum(_pi), 'part_name': _pname(part_names, _pi), 'lo': nd.min(axis=0), 'hi': nd.max(axis=0), # ★ 1단계(용접 씨앗점 고정)용 — 메시 전 CAD 표면을 그대로 들고 간다. # 씨앗점은 이 CAD 로만 정해지므로 메시를 바꿔도 같은 씨앗이 나온다. 'cadV': Vc, 'cadF': Fc, }) nlist.append(nd) elist.append(el + n0) walls.append(tm['metrics']['wall_mm']) esizes.append(tm['elem_size_mm']) coarse = coarse or tm['coarsened'] qual = tm.get('quality', '') flips += tm['flipped'] n0 += len(nd) # ★ 쓸 수 있는 부품이 하나도 없으면 외곽상자로 둔갑시키지 않는다. # 어느 부품이 왜 깨졌는지 그대로 돌려준다(요구 3). if not parts: shutil.rmtree(wd, ignore_errors=True) raise RuntimeError('형상 처리에 성공한 부품이 없습니다 — ' + _failed_text(failed_parts)) # 요소 시작번호 다시 계산 acc = 0 for pi, e in enumerate(elist): parts[pi]['e0'] = acc acc += len(e) nodes = np.vstack(nlist) elems = np.vstack(elist) surf = surface_tris_tet(elems[:, :4]) # ── 용접 접합면 찾기 (방정식 생성은 경계조건 확정 뒤에) ── if len(parts) > 1: # ★ 접촉 허용오차 3.0mm 의 근거(실측, 점↔삼각형 정확거리 481쌍): # 부품 간 실제 간격 중앙값 1.72mm, <=1.5mm 49% / <=3.0mm 53%. # 1.5mm 면 그래프가 5덩어리로 갈라지고(최대 131/146), 3.0mm 면 3덩어리(138/146). wtol = max(min(walls) * 1.2, 3.0) weld_list = find_weld_pairs(parts, wtol) from scipy.spatial import cKDTree as _KD for w in weld_list: A = parts[w['i']] B = parts[w['ind_part']] f6 = B['faces6'] ctr = nodes[f6[:, :3]].mean(axis=1) # ★ 접합면은 '상대 부재 표면에 실제로 붙어 있는 면' 만. # 넉넉한 상자로 고르면 옆면까지 딸려 들어와 전달력이 수십 배로 튄다(실측 36,251N vs 실제 1,000N). # ★ 맞닿은 면 판정은 '두 부재의 바깥법선이 서로 마주보는가' 로 한다. # 앞서 쓴 '상대 절점 쪽 방향' 은 두 면이 겹쳐 있을 때 그 방향이 면 안쪽(접선)으로 # 나와 버려서 무의미하다 — 실측으로 접합면의 5%(24.8/480mm2)만 잡혔다. dA, iA = _KD(A['faces_ctr']).query(ctr) facing = np.einsum('ij,ij->i', B['faces_n'], A['faces_n'][iA]) # ① 접합(MPC)용 — 넉넉하게. 좁히면 결합이 성겨져 가짜 응력첨두가 생긴다(실측 8,028MPa) lo, hi = w['zone_min'] - wtol, w['zone_max'] + wtol sel_mpc = np.where(np.all((ctr >= lo) & (ctr <= hi), axis=1))[0] # ② 전달력 적분용 — 진짜 맞닿은 면만(가깝고 + 법선이 상대 부재를 향함) # ★ 거리 기준을 절점간격보다 작게(2.9mm) 잡으면 안 된다. # A 표면 '절점' 까지의 거리라서, 접합면 위에 딱 붙어 있어도 절점 사이 # 한가운데 있으면 탈락한다. 실측: 480mm2 접합면 중 24.8mm2(5%)만 잡혀 # 부분면 적분값을 전체 전달력으로 오인했다(2,256N). # → 거리는 넉넉히(wtol), 대신 법선이 상대 부재를 향하는지로 옆면을 걸러낸다. # ★ B-6 버그 수정: dA 는 '두 부재 면중심끼리의 거리' 다. 성긴 메시에서 # 면 하나가 수백 mm 라 맞닿은 면도 중심거리가 wtol 을 넘어 전부 탈락했고, # 그래서 전달력이 0건으로 비어 있었다. 점↔삼각형 정확거리로 바꾼다 # (point_tri_dist 는 이미 코드에 있다). 후보 삼각형은 중심 KD 로 좁힌다. # 비용을 억제하려고 '법선이 마주보고 + 중심거리가 wtol 의 6배 안' 인 # 면만 정확거리를 잰다. # ★ 상한은 '어림수 몇 배' 가 아니라 **수학적으로 안전한 경계**를 쓴다: # 삼각형 T 위의 어떤 점까지의 거리가 wtol 이하라면, 중심까지의 거리는 # 반드시 wtol + (중심에서 꼭짓점까지의 최대거리) 이하다. 그러므로 # 이 경계로 거르면 **참인 접촉면을 하나도 놓치지 않는다.** # 종전에 쓰려던 'wtol 의 6배'(약 18mm)는 B-6 이 지적한 그 함정을 # 그대로 되살린다 — 성긴 메시의 면 하나가 수백 mm 라 참인 접촉면이 # 정확거리를 재 보기도 전에 탈락한다. _Avtx = nodes[A['faces6'][:, :3]] # (nA,3,3) _Arad = float(np.linalg.norm( _Avtx - A['faces_ctr'][:, None, :], axis=2).max()) if len(_Avtx) else 0.0 _pre = np.where((facing < -0.7) & (dA <= wtol + _Arad))[0] dA_tri = np.full(len(ctr), np.inf) if len(_pre): _nc = min(8, len(A['faces_ctr'])) _, _ii = _KD(A['faces_ctr']).query(ctr[_pre], k=_nc) _ii = _ii.reshape(len(_pre), _nc) _Atri = _Avtx _pp = ctr[_pre] _best = np.full(len(_pre), np.inf) for _c in range(_nc): # 반복 8회, 각 회는 완전 벡터화 _best = np.minimum(_best, point_tri_dist_pairs(_pp, _Atri[_ii[:, _c]])) dA_tri[_pre] = _best sel_f = np.where(dA_tri <= wtol)[0] # ★ 접촉 판정에서 '실제로 절점이 닿은 면'(near_faces)은 반드시 넣는다. # 상자 필터만 쓰면 접합면이 상자 경계에 걸쳐 있을 때 정작 종속절점이 # 올라탈 면이 빠져 MPC 가 0개가 된다. nf = w.get('near_faces') if nf is not None and len(nf): sel_mpc = np.union1d(sel_mpc, nf) if len(sel_mpc): w['face_idx'] = sel_mpc w['force_idx'] = sel_f w['wtol'] = wtol # ★ 접합면 면적을 실제로 잰다. 종전 결과의 'iface_area_mm2' 는 # **어디에서도 채워지지 않는 키**라 언제나 0.0 이었다(화면이 빈 # 숫자를 보여 주고 있었다). 선별이 잘 되고 있는지 확인할 수 있는 # 유일한 숫자이므로 반드시 채운다. if len(sel_f): _ft = nodes[f6[sel_f][:, :3]] w['iface_area'] = float(0.5 * np.linalg.norm( np.cross(_ft[:, 1] - _ft[:, 0], _ft[:, 2] - _ft[:, 0]), axis=1).sum()) else: w['iface_area'] = 0.0 # ── 1단계: 용접 씨앗점 고정 (메시와 무관한 묶음) ────────────── # 여기서 w['dep_nodes'] 를 '씨앗점마다 최근접 절점 1개' 로 **갈아끼운다**. # 면(face) 선별은 그대로 둔다 — 종속절점은 여전히 상대 '면' 에 tri6 # 가중치로 묶이므로(절점↔절점 강체막대 금지) 기하적으로 올바르다. if weld_seed: _seeds, _sstat = weld_seeds_cad(parts, wtol) _mpair = set((min(w['i'], w['ind_part']), max(w['i'], w['ind_part'])) for w in weld_list) _spair = set(k for k, v in _seeds.items() if len(v)) _sstat['pairs_mesh'] = len(_mpair) _sstat['pairs_cad'] = len(_spair) _sstat['pairs_both'] = len(_mpair & _spair) _sstat['pairs_cad_only'] = len(_spair - _mpair) _sstat['pairs_mesh_only'] = len(_mpair - _spair) _noseed = 0 for w in weld_list: _k = (min(w['i'], w['ind_part']), max(w['i'], w['ind_part'])) _sp = _seeds.get(_k) if _sp is None or not len(_sp): w['dep_nodes'] = np.array([], dtype=np.int64) w['dep_nodes_rev'] = np.array([], dtype=np.int64) _noseed += 1 continue _Ap, _Bp = parts[w['i']], parts[w['ind_part']] # ★ 2a: 후보 풀을 중간절점까지 넓힌 'surf_*_all' 로 바꾼다. # (옛 키가 없는 옛 캐시·경로도 깨지지 않게 폴백을 둔다) _Apt = _Ap.get('surf_pts_all', _Ap['surf_pts']) _Bpt = _Bp.get('surf_pts_all', _Bp['surf_pts']) _Aid = _Ap.get('surf_ids_all', _Ap['surf_ids']) _Bid = _Bp.get('surf_ids_all', _Bp['surf_ids']) _dA2, _iA2 = _assign_unique(_Apt, _sp) _dB2, _iB2 = _assign_unique(_Bpt, _sp) w['dep_nodes'] = np.unique(_Aid[_iA2]) w['dep_nodes_rev'] = np.unique(_Bid[_iB2]) w['seed_n'] = int(len(_sp)) # ★ 결함 2 고침: 두 씨앗이 같은 절점에 떨어지면 종전에는 unique 가 # 하나를 조용히 삼켜 묶음 수가 줄었다(507건, 14.8%). 이제는 # 두 번째·세 번째 최근접 절점으로 밀어 **씨앗 하나당 절점 하나**를 # 보장한다 — 묶음 수가 씨앗 수로 고정되어야 단계 간 비교가 된다. w['seed_collide'] = int(len(_sp) - len(w['dep_nodes'])) w['seed_snap_max_mm'] = round(float(max(_dA2.max(), _dB2.max())), 2) _sstat['welds_without_seed'] = _noseed _sstat['collisions'] = int(sum(w.get('seed_collide', 0) for w in weld_list)) _seed_stat = _sstat _seed_pairs_out = {(parts[k[0]]['part_no'], parts[k[1]]['part_no']): v for k, v in _seeds.items() if len(v)} # ★ 요소 품질 전수 측정(계획 D-1-1). 이 값이 없으면 '쓸 수 없는 메시' 에서 # 나온 응력을 판정으로 내보내게 된다. _tq = tet_quality(nodes, elems) if etype == 'C3D10' else {} _sjarr = _tq.pop('_sj', None) _tq.pop('_ar', None) # ★ 실측 요소크기로 '벽당 겹수' 를 다시 잰다. 종전 값은 평균의 평균이라 # 실제의 2.6배로 나왔다(계획 F-1). _h_real = _tq.get('elem_size_median_mm') or float(np.mean(esizes)) meshinfo = { 'parts': len(parts), 'wall_mm': round(float(np.mean(walls)), 3), 'elem_size_mm': round(float(np.mean(esizes)), 3), # 종전 정의(부품별 평균의 평균) — 호환용으로 남긴다. 판단에는 쓰지 말 것. 'through_wall_meanmean': round(float(np.mean(walls) / np.mean(esizes)), 2), 'elem_size_real_mm': round(float(_h_real), 3), 'through_wall': round(float(np.median(walls) / max(_h_real, 1e-9)), 2), 'quality_metrics': _tq, 'coarsened': coarse, 'flipped': flips, 'quality': qual, 'welded_verts': len(V), 'tris': len(F), 'welds': 0, 'parts_in': len(part_sizes) if part_sizes else len(comps), 'parts_failed': len(failed_parts), # ★ 두께를 '추정' 해 압출한 부재 — 가정이므로 숨기지 않고 화면에 적는다 # (볼록껍질 대체를 '껍질 근사' 로 적는 것과 같은 이유). 'sheet_est': _sheet_est, # 화면(fab_fea.js)이 쓰는 기계판독용 자료 — 문장만 주면 회색처리·라벨을 못 붙인다 'repair': _repair_stat, # 수리 방법별 부품 수 'hull_parts': _hull_parts, # '껍질 근사' 로 대체된 원래 부품번호 'bulk': bool(_bulk), 'coarse_note': (COARSE_NOTE if (_bulk and not quality_bulk) else ''), 'quality_bulk': bool(_bulk and quality_bulk), 'weld_seed': _seed_stat, # 국부 재해석 전용 — 부재마다 목표 요소크기를 맞췄는지. 집합체 중앙값만 # 보면 '일부 부재는 안 가늘어졌다' 를 놓친다(실측으로 당한 함정). 'per_part_h': _hstat, 'h_target_mm': (submodel.get('h_target') if submodel else None), 'h_miss_parts': [x['no'] for x in _hstat if x.get('h_hit') is False], # 밀도 상향 부재 — 화면에 '왜 이 부재만 촘촘한가' 를 실측으로 적는다 'fine_parts': _fine_stat, } if failed_parts: warn = _failed_text(failed_parts) # ── 외곽상자 모드는 '검증용 비교' 로만 남아 있다(force_obb=True). ─────────── # 자동 후퇴는 삭제했다 — 실패를 네모 덩어리 결과로 둔갑시켜 보여주는 것이 # 아무 결과도 안 주는 것보다 나쁘기 때문이다. if force_obb: method, etype = 'obb_forced', 'C3D8I' parts, weld_list, tie_blocks, failed_parts = [], [], [], [] nodes_l, hexes = box_hex_mesh(size, 20, 4, 4) nodes = center + (nodes_l - size / 2.0) @ R.T elems = hexes surf = surface_tris_hex(elems) meshinfo = {'parts': 1, 'welds': 0, 'forced_obb': True} _nbudget = (1.1e6 if quality_bulk else (budget_nodes or MAX_C3D10_NODES) * 1.3) if len(nodes) > _nbudget: shutil.rmtree(wd, ignore_errors=True) # ★ 왜 못 푸는지, 지금 상태가 어떤지 숫자로 말한다. '너무 큽니다' 한 마디로는 # 운전자가 다음에 뭘 해야 할지 모른다(실측: 291부품 → 절점 1,077,806개). raise RuntimeError( '선택한 부품이 너무 많습니다 — 절점 %s개로 한 번에 풀 수 있는 양(%s개)의 %.1f배입니다. ' '부품 %d개 중 %d개는 형상 처리에 성공했고 %d개는 실패했습니다. ' '부재를 10~30개씩 나눠서 해석하세요(왼쪽 툴바에서 일부만 선택).%s' % ('{:,}'.format(len(nodes)), '{:,}'.format(int(budget_nodes or MAX_C3D10_NODES)), len(nodes) / float(budget_nodes or MAX_C3D10_NODES), (meshinfo.get('parts_in') or len(parts)), len(parts), len(failed_parts), (' ※ ' + _failed_text(failed_parts)) if failed_parts else '')) # ── 경계조건 ── bbox = nodes.max(axis=0) - nodes.min(axis=0) r_def = pick_radius or max(bbox.max() * 0.05, meshinfo.get('elem_size_mm', 5) * 3, 3.0) loads = {} prescribed = {} # node0idx -> [(dof(1~3), mm값), ...] 강제변위 모드 disp_node_set = set() # 강제변위가 걸린 절점 — 용접 MPC 제외 + 하중점 취급용 bc_info = {'mode': 'click', 'fix_nodes': 0, 'load_nodes': 0, 'radius_mm': round(r_def, 2)} req_sum = 0.0 # 운전자가 요청한 하중 총합(N) drop_sum = 0.0 # 근처 절점이 없어 실리지 못하고 사라진 하중(N) overlap_nodes = 0 # 고정점이면서 하중점이기도 했던 절점 수(하중에서 배제한 수) # ★ F-2 — '실제로 쓴' 경계조건 반경 목록. 종전엔 bc_info 에 기본 반경(r_def) 하나만 # 담겨, 화면의 특이점 제외가 실제 반경(하중점 60mm)과 다른 값으로 돌아갔다. # (보고서 C-6: 구조 본체까지 과하게 지워져 '구조 최대' 가 실제보다 작게 보임) bc_spheres = [] # [([x,y,z], r_mm), ...] # ── 2단계: 서브모델 경계 강제변위 ──────────────────────────────────── # 전체 해석에서 '잘린 면' 에 해당하는 좌표(CAD 용접 씨앗점)를 받아, 그 근처 # 국부 절점에 전체 해석의 변위를 보간해 3방향 모두 *BOUNDARY 로 지정한다. # 이렇게 하면 국부 영역 밖 구조의 강성을 '변위 이력' 으로 대신 받는다. sub_info = {} _drv_src_nodes = {} # 접합면키 -> 그 접합면 씨앗이 강제한 절점들(결함 3 검증용) # ★ 1단계: 강제된 절점이 **어느 부품 쪽**인가. 전체 해석의 접합면 힘은 # '종속절점이 있는 부품' 의 절점력 합이므로, 국부 쪽 절점이 반대 부품에 # 있으면 두 값은 부호가 뒤집힌 채 비교된다(정의상 200% 어긋남). _drv_src_part = {} # 접합면키 -> {부품번호: 강제절점 수} if submodel: _dp = np.asarray(submodel['drive_pts'], float).reshape(-1, 3) _dr = float(submodel.get('drive_r') or 0.0) # ★★ 왜 '반경 안의 모든 절점' 이 아니라 '씨앗점마다 1개' 인가 (실측 2026-10-10) # 반경 15mm 로 잡으면 강제된 절점이 252개 → 771개(메시를 조이면 3배)로 늘었고, # 외부가 해 준 일(Σ반력·변위)이 89,165 → 171,486 Nmm 로 92% 튀었다. 강제된 # 절점 수가 메시에 딸려 변하면 '경계조건이 같은 모델' 이 아니라서 수렴을 # 판정할 수 없다 — 1단계 용접 씨앗점과 똑같은 함정이다. # → 씨앗점마다 최근접 절점 1개만 강제한다(용접 MPC 와 같은 이산화). if _dr > 0: _hl = [nodes_near(nodes, _p, _dr) for _p in _dp] _hl = [x for x in _hl if len(x)] _hit = np.unique(np.concatenate(_hl)) if _hl else np.array([], dtype=int) else: # ★ 결함 2·3: 두 씨앗이 같은 절점에 떨어지면 종전 unique 가 하나를 # 삼켜 강제 절점 수가 메시에 따라 달라졌다. 겹치지 않게 배정한다. _dd2, _ii2 = _assign_unique(nodes, _dp) _hit = np.asarray(_ii2, np.int64) if len(_hit): _du, _dmode = interp_disp(submodel['nodes'], submodel['elems'], submodel['disp'], nodes[_hit]) for _q, _ii0 in enumerate(_hit): _ii = int(_ii0) for _dof in (1, 2, 3): prescribed.setdefault(_ii, []).append((_dof, float(_du[_q, _dof - 1]))) disp_node_set.add(_ii) # 씨앗 1개 ↔ 절점 1개 대응이 살아 있을 때만 접합면별로 묶는다. _dsrc = submodel.get('drive_src') if _dsrc is not None and _dr <= 0 and len(_dsrc) == len(_hit): _pbnd = np.array([pp['n0'] for pp in parts] + [parts[-1]['n1']], dtype=np.int64) for _q, _ii0 in enumerate(_hit): _kk0 = tuple(int(x) for x in _dsrc[_q]) _drv_src_nodes.setdefault(_kk0, []).append(int(_ii0)) _pk = int(np.searchsorted(_pbnd, int(_ii0), side='right')) - 1 if 0 <= _pk < len(parts): _d0 = _drv_src_part.setdefault(_kk0, {}) _pn = int(parts[_pk]['part_no']) _d0[_pn] = _d0.get(_pn, 0) + 1 sub_info = {'drive_pts': int(len(_dp)), 'drive_r_mm': round(_dr, 2), 'driven_nodes': int(len(_hit)), 'driven_dup': int(len(_dp) - len(_hit)) if _dr <= 0 else None, 'interp_in_tet': int((_dmode == 1).sum()), # 사면체 '면 위' 로 투영해 보간한 수(종전 '최근접 절점 복사' 대체) 'interp_clamped': int((_dmode == 2).sum()), 'interp_nearest': int((_dmode == 0).sum()), 'snap_max_mm': round(float(np.max(_dd2)), 2) if _dr <= 0 else None, 'driven_disp_max_mm': round(float(np.linalg.norm(_du, axis=1).max()), 4)} else: sub_info = {'drive_pts': int(len(_dp)), 'drive_r_mm': round(_dr, 2), 'driven_nodes': 0} if fix_points or load_points: fix = np.unique(np.concatenate( [nodes_near(nodes, f['p'], f.get('r') or r_def) for f in (fix_points or [])] or [np.array([], dtype=int)])) # ★ 고정 절점은 하중 집합에서 뺀다(F-1). # 왜: 고정구(球)와 하중구가 겹치면 그 절점에 *BOUNDARY 와 *CLOAD 가 동시에 # 붙는다. CalculiX 는 구속을 적용하므로 그 하중은 반력으로 바로 흡수되어 # 구조에 전혀 실리지 않는데, 종전 코드는 applied_N 에 그대로 세고 # drop_sum 은 0 을 유지했다 → "하중 소멸 0%" 가 거짓이 되는 결함. # ※ 이 결함은 반력 검증(Σ반력=Σ하중)으로는 안 잡힌다 — 흡수된 하중도 # 고정 절점의 반력에 그대로 나타나 검증은 통과해 버린다. 그래서 둘 다 필요. _fixset = set(int(x) for x in fix) for f in (fix_points or []): bc_spheres.append([[float(x) for x in f['p']], float(f.get('r') or r_def)]) for lp in (load_points or []): bc_spheres.append([[float(x) for x in lp['p']], float(lp.get('r') or r_def)]) for lp in (load_points or []): idx = nodes_near(nodes, lp['p'], lp.get('r') or r_def) if lp.get('mm') is None: req_sum += abs(float(lp.get('N', load_N))) if not len(idx): # 하중점 근처에 절점이 없어 조용히 버려지던 지점 — 누락분으로 집계한다 if lp.get('mm') is None: drop_sum += abs(float(lp.get('N', load_N))) continue d = np.asarray(lp.get('dir') or [0, 0, -1], float) nrm = np.linalg.norm(d) d = d / nrm if nrm > 0 else np.array([0., 0., -1.]) if lp.get('mm') is not None: # ── 강제변위 모드: N 대신 mm(몇 mm 밀기)가 들어옴 ── # dir의 주된 축 1개 DOF에만 부호 포함 값(mm*sign)을 지정한다. # 3DOF 모두 고정하면 옆방향까지 구속돼 가짜 응력이 생기므로 1DOF만 건다. axis = int(np.argmax(np.abs(d))) sign = 1.0 if d[axis] >= 0 else -1.0 val = float(lp['mm']) * sign dof = axis + 1 for i in idx: ii = int(i) prescribed.setdefault(ii, []).append((dof, val)) disp_node_set.add(ii) else: # per 는 '원래 절점 수' 로 나눈다 — 살아남은 절점에 재분배하면 # 결함(하중이 고정점에 흡수됨)이 숫자에서 사라져 버린다. per = float(lp.get('N', load_N)) / len(idx) for i in idx: ii = int(i) if ii in _fixset: overlap_nodes += 1 drop_sum += abs(per) continue loads[ii] = loads.get(ii, np.zeros(3)) + d * per elif submodel: # 경계 강제변위가 구속을 대신한다 — 자동 외팔보를 걸면 안 된다 fix = np.array([], dtype=int) bc_info['mode'] = 'submodel' else: fix, tip, d = auto_cantilever(nodes, center, R, size) per = load_N / max(len(tip), 1) for i in tip: loads[int(i)] = -d * per bc_info['mode'] = 'auto_cantilever' req_sum = abs(float(load_N)) if len(parts) > 1: # 부재가 여러 개인데 조건을 안 찍으면, 전체 외곽의 한쪽 끝만 고정된다. # 멀리 있는 부재는 용접 몇 점으로만 매달려 말이 안 되는 큰 변위가 나온다. bc_info['warn'] = ('부재 %d개인데 고정점을 안 찍었습니다. 전체 한쪽 끝만 자동 고정되어 ' '멀리 있는 부재가 매달린 상태로 계산됩니다 — 실제 취부점을 찍어주세요.' % len(parts)) # ★ "절점이 없다" 를 뭉개지 않는다(요구 5). # 같은 증상에 원인이 둘이고 운전자가 해야 할 일이 서로 다르다. # ① 그 자리에 있던 부품의 형상 처리가 실패했다 → 그 부품 이름을 말해줘야 한다 # ② 애초에 그 자리에 형상이 없다(외곽상자 모서리 = 허공) → 직접 찍으라고 해야 한다 if not len(fix) and not prescribed: shutil.rmtree(wd, ignore_errors=True) if submodel: # 경계 씨앗점 주변에 절점이 하나도 없다 = 경계 부재들이 전부 메시 실패했다. # '고정점을 찍어주세요' 로 뭉개면 원인을 영원히 못 찾는다. raise RuntimeError( '국부 영역 경계에 절점이 없습니다 — 경계 씨앗점 %d개 주변 %.0fmm 안에 ' '절점이 하나도 없습니다(형상 처리 실패 %d개). %s' % (len(np.asarray(submodel['drive_pts'], float).reshape(-1, 3)), float(submodel.get('drive_r') or 0.0), len(failed_parts), _failed_text(failed_parts))) raise RuntimeError(_no_node_msg('고정점', fix_points, failed_parts, part_sizes, P)) if not loads and not prescribed: shutil.rmtree(wd, ignore_errors=True) raise RuntimeError(_no_node_msg('하중점', load_points, failed_parts, part_sizes, P)) bc_info['fix_nodes'] = int(len(fix)) bc_info['load_nodes'] = int(len(loads)) if prescribed: bc_info['mode'] = 'disp' bc_info['prescribed_nodes'] = int(len(disp_node_set)) # ── 용접 접합 방정식 생성 ── # ★ 반드시 고정절점(fix)이 정해진 뒤에 만든다. 한 절점이 MPC 종속이면서 동시에 # 구속(SPC)이면 CalculiX 가 "detected on the dependent side of a MPC and a SPC" 로 거부한다. # 강제변위 노드도 SPC와 같은 성격(지정 변위)이라 MPC 종속절점에서 제외해야 한다. # ★★ spc_set 과 fixed_set 을 반드시 구분한다. # fixed_set 은 아래 용접 루프에서 '이미 종속으로 쓴 절점' 을 계속 더해 나간다 # (한 절점을 두 번 종속시키면 CalculiX 가 거부하므로 필요한 장치다). # 그래서 루프가 끝난 fixed_set 에는 용접 종속절점 수천 개가 섞여 있다. 이것을 # '고정된 절점' 으로 오인하면 # · 고정도 안 된 부재가 '고정됨' 으로 보여 공중에 뜬 것을 못 잡고 # · 하중-고정 연결성 검사도 통째로 무력해진다(아무 부재나 고정 덩어리로 보임). # → 진짜 구속(SPC/강제변위)만 담은 spc_set 을 따로 떠 두고 점검에는 그것만 쓴다. spc_set = set(int(x) for x in fix) | disp_node_set fixed_set = set(spc_set) # ★ 계측 전용 스냅샷(해석 불변). 용접 루프가 fixed_set 에 종속절점을 쌓으므로, # '버려진 씨앗' 안에서 **원래 구속절점**(고칠 수 없음)과 **앞 쌍이 먼저 가져감** # (배정 결함)을 가르려면 루프 진입 전 집합이 필요하다. fixed_set0 = frozenset(spc_set) if weld_mode != 'none': # ★ 재현성 [2026-10-11 통제실험]: 아래 루프는 fixed_set 을 누적하며 # '한 절점은 최대 한 번만 종속' 규칙을 적용한다. 그래서 **접합면을 보는 # 순서가 어느 절점이 용접 MPC 종속이 되는지를 결정한다**. weld_list 는 # find_weld_pairs 의 i 464.2MPa, +5.0%. 절점·묶인절점·씨앗은 동일). # -> 접합면을 **기하 키**(접촉영역 중심 좌표)로 정렬해 부품 순서와 # 무관하게 같은 종속절점 배정이 나오게 한다. 좌표는 CAD 값이므로 # 부품을 어떤 순서로 보내도 같다. # ※ 숫자를 맞추려는 조정이 아니다 - 순서 의존성 자체를 없앤다. def _weld_order_key(w): c = np.asarray(w.get('center'), float).reshape(-1) dn = w.get('dep_nodes') return (round(float(c[0]), 3), round(float(c[1]), 3), round(float(c[2]), 3), round(float(w.get('length_mm') or 0.0), 3), int(0 if dn is None else len(dn))) weld_list = sorted(weld_list, key=_weld_order_key) for w in weld_list: if 'face_idx' not in w: continue B = parts[w['ind_part']] eqs, cnt, _wl = build_weld_equations(nodes, w['dep_nodes'], B['faces6'][w['face_idx']], fixed_set, w['wtol'] * 1.5, fixed_set0) _wside = int(parts[w['i']]['part_no']) # ★ 방정식이 0개면 역할을 뒤집어 한 번만 재시도한다. # 원인: '한 절점은 최대 한 번만 종속'(fixed_set 누적) 규칙 때문에, # 표면절점 8개짜리 작은 부품이 이웃 3개와 닿으면 1번만 종속될 수 있고 # 나머지 쌍이 방정식 0개가 된다. 상대 부재의 절점은 아직 소비되지 않았으므로 # 역할을 뒤집으면 살아난다. (방정식이 이미 생긴 쌍은 건드리지 않는다) if cnt == 0: Aw = parts[w['i']] nf2 = w.get('near_faces_rev') if nf2 is not None and len(nf2): eqs2, cnt2, _wl2 = build_weld_equations( nodes, w['dep_nodes_rev'], Aw['faces6'][nf2], fixed_set, w['wtol'] * 1.5, fixed_set0) if cnt2 > 0: eqs, cnt, _wl = eqs2, cnt2, _wl2 w['flipped_role'] = True # 역할을 뒤집었으면 종속절점은 상대 부품 쪽이다(부호 기준) _wside = int(parts[w['ind_part']]['part_no']) for e in eqs: fixed_set.add(e[0]) # 같은 절점을 두 번 종속시키지 않는다 tie_blocks.extend(eqs) w['tied_nodes'] = cnt w['loss'] = _wl # 계측(해석 불변) — 손실 분해 표 w['side_part_no'] = _wside # 종속절점이 어느 부품 쪽인가(부호 기준) # ★ 이 접합면을 '실제로' 묶은 종속절점 번호. *EQUATION 종속절점의 절점력은 # 그 구속이 전달하는 힘 그 자체다 → 접합면별 단면력을 잴 수 있다 # (국부 재해석의 경계 강제가 타당한지 검증하는 유일한 숫자). w['tied_ids'] = np.array([int(e[0]) for e in eqs], dtype=np.int64) # ★ 4단계: 화면이 '씨앗 몇 개 중 몇 개를 못 묶었나' 를 사유별로 말할 수 # 있게 쌍별 loss 를 합쳐 둔다. 계산 안 한 것을 숨기지 않기 위한 값이다. _ls = {'seed': 0, 'tied': 0, 'compete': 0, 'too_far': 0, 'no_face': 0, 'spc': 0, 'degen': 0, 'no_seed': 0} for w in weld_list: _ls['seed'] += int(w.get('seed_n') or 0) _ls['tied'] += int(w.get('tied_nodes') or 0) if 'face_idx' in w and not int(w.get('seed_n') or 0): _ls['no_seed'] += 1 _L = w.get('loss') if _L: _ls['compete'] += _L['taken_drop'] _ls['too_far'] += _L['too_far'] _ls['no_face'] += _L['no_face'] _ls['spc'] += _L['spc_drop'] _ls['degen'] += _L['degen'] _ls['lost'] = int(_ls['seed'] - _ls['tied']) meshinfo['weld_loss'] = _ls if meshinfo.get('parts', 1) > 1: meshinfo['welds'] = sum(w.get('tied_nodes', 0) for w in weld_list) # ── 부재가 제대로 붙어 있는지 점검 ── # 살짝 스치기만 한 부재는 용접 이음이 몇 점밖에 안 생겨, 사실상 매달린 상태가 된다. # 이러면 CalculiX 가 오류를 내지 않고 '엄청나게 큰 변위' 로 그럴듯하게 풀어 버린다. # (실측: 부재 8개 조립품에서 변위 1,891mm 가 나옴) → 미리 잡아서 경고한다. if len(parts) > 1: # 절점번호 → 부품 자리번호. 부품이 수백 개면 선형탐색이 비싸니 경계배열로 찾는다. _bnd = np.array([pp['n0'] for pp in parts] + [parts[-1]['n1']], dtype=np.int64) def _part_of(nid): k = int(np.searchsorted(_bnd, nid, side='right')) - 1 return k if 0 <= k < len(parts) else None # ── 고정까지 이어지지 않은 '떠 있는 덩어리' 점검 ─────────────────────── # ★ 용접으로 이어진 덩어리에 구속점이 하나도 없으면 그 덩어리는 공중에서 # 자유롭게 움직인다(강체운동 6자유도) → 전체 강성행렬이 특이해지고, # CalculiX 는 오류 없이 터무니없는 변위를 그럴듯하게 돌려준다. # 실측: 작은 부재 하나가 이음 2개로만 붙어 있어 변위 66,534,314 mm. # (부품 1개가 형상 실패로 빠졌는데 그게 '다리' 였던 경우에도 똑같이 갈라진다) # → 풀기 전에 잡아서, 실패를 결과로 둔갑시키지 않는다. # ★ 함수로 둔 이유: 아래에서 '다리 잇기' 로 방정식을 더한 뒤 **같은 방법으로 # 다시** 세어야 한다. 두 번 베껴 쓰면 한쪽만 고치는 사고가 난다. def _conn(): pof = list(range(len(parts))) def find(x): while pof[x] != x: pof[x] = pof[pof[x]] x = pof[x] return x for (sn, fn6, _w) in tie_blocks: a, b = _part_of(int(sn)), _part_of(int(fn6[0])) if a is None or b is None: continue ra, rb = find(a), find(b) if ra != rb: pof[rb] = ra fixg = set() for f in spc_set: # ★ fixed_set 이 아니라 진짜 구속절점만 pi = _part_of(int(f)) if pi is not None: fixg.add(find(pi)) loadg = set() for nid in list(loads.keys()) + list(disp_node_set): pi = _part_of(int(nid)) if pi is not None: loadg.add(find(pi)) flo, flo_loaded = [], False for pi, pp in enumerate(parts): g = find(pi) if g not in fixg: flo.append((pp['part_no'], pp['part_name'])) if g in loadg: flo_loaded = True return flo, flo_loaded, [find(i) for i in range(len(parts))], fixg, loadg floating, float_loaded, _grp, fix_groups, load_groups = _conn() # ── 다리 잇기: 떠 있는 덩어리를 '실측 간격' 근거로만 이어 본다 ────────── _bridges, _bridge_rej = [], [] # ★ 통제실험 전용 스위치(2026-10-11). 환경변수 DZW_FEA_NOBRIDGE=1 이면 다리 # 잇기를 건너뛴다. 기본값(변수 없음)에서는 한 줄도 달라지지 않는다 — # '변위가 +43% 커진 원인' 을 가르려면 다리만 끄고 같은 요청을 돌려야 하는데 # 그 스위치가 없으면 코드를 임시로 고쳐야 하고, 그러면 실험값과 수정 사고를 # 구별할 수 없게 된다. 숫자(tol·간격)는 하나도 건드리지 않는다. _nobridge = os.environ.get('DZW_FEA_NOBRIDGE') == '1' if _nobridge and floating and _logf: _logf(' [다리] DZW_FEA_NOBRIDGE=1 — 다리 잇기를 건너뛴다(통제실험)') if floating and weld_mode != 'none' and not _nobridge: _bridges, _bridge_rej = bridge_floating( parts, nodes, _grp, fix_groups, fixed_set, tie_blocks, float(locals().get('wtol') or 3.0), logf=_logf, fixed_set0=fixed_set0) if _bridges: floating, float_loaded, _grp, fix_groups, load_groups = _conn() if _logf: _logf(' [다리] %d군데를 이어 떠 있는 부재 %d개가 남았다' % (len(_bridges), len(floating))) # ★ tied_ids 는 numpy 배열이라 JSON 으로 못 나간다 — 화면용에서는 뺀다. # (내부에서는 _bridges 가 그대로 들고 있어 접합면 힘 계산에 쓰인다) meshinfo['weld_bridges'] = [{k: v for k, v in b.items() if k != 'tied_ids'} for b in _bridges] # ★ 다리로 만든 방정식도 '묶인절점' 에 더한다. 안 더하면 화면의 묶인절점 수가 # 실제로 쓴 구속보다 적게 나와, 전·후 비교에서 사라진 숫자가 생긴다. if _bridges: meshinfo['weld_bridge_eqs'] = int(sum(b['eqs'] for b in _bridges)) meshinfo['welds'] = int(meshinfo.get('welds') or 0) + meshinfo['weld_bridge_eqs'] # 못 이은 쌍 중 가장 가까운 간격 — 제외 사유에 실측값으로 적는다 _gap_of = {} for r in _bridge_rej: for (x, y) in ((r['a'], r['b']), (r['b'], r['a'])): gv = r.get('gap_mm') if gv is not None and (x not in _gap_of or gv < _gap_of[x][0]): _gap_of[x] = (gv, y, r.get('why') or '') # ── 이음이 거의 없는 부재 경고 (다리 잇기 뒤의 상태로 센다) ────────── tied_cnt = {} for (sn, fn6, _w) in tie_blocks: for nid in (int(sn), int(fn6[0])): pi = _part_of(nid) if pi is not None: tied_cnt[pi] = tied_cnt.get(pi, 0) + 1 weak = [] for pi, pp in enumerate(parts): has_fix = any(pp['n0'] <= f < pp['n1'] for f in spc_set) if not has_fix and tied_cnt.get(pi, 0) < 3: # ★ pi+1 이 아니라 원래 부품번호. 실패한 부품을 건너뛰면 자리번호가 밀린다. weak.append((pp['part_no'], pp['part_name'])) if weak: bc_info['weak_parts'] = [w[0] for w in weak] bc_info['warn'] = (('부재 %s 가 고정도 안 되고 용접 이음도 거의 없습니다(사실상 공중에 뜬 상태). ' '결과 변위·응력을 믿지 마세요.') % ', '.join('%d번(%s)' % w for w in weak)) # ── 접합면별 묶인절점 덤프(환경변수 DZW_FEA_PAIRDUMP 가 있을 때만) ── # 왜 필요한가 (2026-10-11): '묶인절점 2,647 → 2,212' 처럼 합계만 보면 # 줄어든 435개가 '어느 쌍에서' 줄었는지 알 수 없다. 용접 루프는 # fixed_set('한 절점은 한 번만 종속')을 누적하므로 부품 구성이 바뀌면 # 종속절점 쟁탈 결과가 재배분된다. 두 실행에 **공통인 쌍만** 골라 빼면 # '재배분 때문인가, 새 부품의 쌍 때문인가' 가 추측 없이 갈린다. # ※ 자동 제외 뒤 재귀 실행에서도 덮어쓴다 — 마지막으로 **실제 푼** 모델의 # 배정이 알고 싶은 값이기 때문이다(DIAG 와 의도가 반대다). if os.environ.get('DZW_FEA_PAIRDUMP'): try: _pd = [] for _w in weld_list: if 'face_idx' not in _w: continue _a = int(parts[_w['i']]['part_no']) _b = int(parts[_w['ind_part']]['part_no']) _c = np.asarray(_w.get('center'), float).reshape(-1) _dn = _w.get('dep_nodes') _pd.append({'a': min(_a, _b), 'b': max(_a, _b), 'center': [round(float(x), 3) for x in _c[:3]], 'tied': int(_w.get('tied_nodes') or 0), 'flipped': bool(_w.get('flipped_role')), # ★ 0단계 계측: 씨앗이 어디서 사라졌나를 쌍마다 남긴다. # 합계만 보면 '28.4% 손실' 의 채널을 가를 수 없다. 'seed_n': int(_w.get('seed_n') or 0), 'dep_n': int(0 if _dn is None else len(_dn)), 'side_part_no': _w.get('side_part_no'), 'loss': _w.get('loss'), 'kind': '접합면'}) for _b2 in _bridges: _pd.append({'a': min(_b2['a'], _b2['b']), 'b': max(_b2['a'], _b2['b']), 'center': None, 'tied': int(_b2['eqs']), 'flipped': False, 'seed_n': int(_b2.get('seeds') or 0), 'dep_n': None, 'side_part_no': _b2.get('side_part_no'), 'loss': _b2.get('loss'), 'assign_cut': _b2.get('assign_cut'), 'kind': '다리'}) json.dump({'parts_meshed': int(len(parts)), 'parts_nos': [int(pp['part_no']) for pp in parts], 'welds_tied_total': int(meshinfo.get('welds') or 0), 'bridge_eqs': int(meshinfo.get('weld_bridge_eqs') or 0), 'floating': [int(t[0]) for t in floating], 'pairs': _pd}, open(os.environ['DZW_FEA_PAIRDUMP'], 'w', encoding='utf-8'), ensure_ascii=False, indent=1) if _logf: _logf(' [쌍덤프] 부품 %d개 / 쌍 %d개 → %s' % (len(parts), len(_pd), os.environ['DZW_FEA_PAIRDUMP'])) except Exception as _e: if _logf: _logf(' [쌍덤프] 실패: %s' % _e) # ── 진단 덤프(환경변수 DZW_FEA_DIAG 가 있을 때만) ──────────────────── # 왜 여기인가: '연결 불가로 자동 제외' 판정이 바로 이 지점에서 끝난다. # 여기서 한 번 받아 적어 두면 400초짜리 전체 해석을 다시 돌리지 않고 # 부품별 원인(표면 간격·씨앗·MPC)을 몇 번이든 다시 따질 수 있다. # 숫자는 전부 이 실행에서 실제로 쓴 값이다(지어낸 값 없음). # ※ floating 이 있을 때만 쓴다. 자동 제외 뒤 재귀 실행에서는 floating 이 # 비므로, 덮어쓰면 '왜 빠졌는지' 가 담긴 1회차 덤프가 사라진다. if os.environ.get('DZW_FEA_DIAG') and floating: try: _L = locals() _dg = { 'wtol_mm': float(_L.get('wtol') or 0.0), 'parts': [{'no': pp['part_no'], 'name': pp['part_name'], 'lo': pp['lo'].tolist(), 'hi': pp['hi'].tolist(), 'wall_mm': float(pp.get('wall') or 0.0), 'cadV': np.asarray(pp['cadV'], np.float32), 'cadF': np.asarray(pp['cadF'], np.int32)} for pp in parts], 'floating': [int(t[0]) for t in floating], 'group': [int(x) for x in _grp], 'bridges': list(_bridges), 'bridge_rejected': list(_bridge_rej), # ★ 메시까지 담는다(떠 있는 부재 + 그 10mm 이웃만). 이게 없으면 # '절점이 접합면에 있나' 같은 질문 하나에 400초짜리 전체 해석을 # 다시 돌려야 한다. 담아 두면 _fea_연결진단.py replay 로 몇 초다. 'nodes': nodes.astype(np.float32), 'mesh': {}, 'fix_groups': sorted(int(x) for x in fix_groups), 'spc_parts': sorted(set(int(_part_of(int(f))) for f in spc_set if _part_of(int(f)) is not None)), 'sheet_est': list(_L.get('_sheet_est') or []), 'failed': list(failed_parts), 'welds': [{'i': int(w['i']), 'j': int(w['ind_part']), 'wtol': float(w.get('wtol') or 0.0), 'seed_n': int(w.get('seed_n') or 0), 'seed_collide': int(w.get('seed_collide') or 0), 'tied': int(w.get('tied_nodes') or 0), 'flipped': bool(w.get('flipped_role')), 'dep_n': int(len(w.get('dep_nodes') if w.get('dep_nodes') is not None else [])), 'iface_area': float(w.get('iface_area') or 0.0)} for w in (_L.get('weld_list') or [])], 'seeds': {(int(k[0]), int(k[1])): int(len(v)) for k, v in (_L.get('_seeds') or {}).items()}, 'tie_pair': {}, } for (sn, fn6, _w) in tie_blocks: _a, _b = _part_of(int(sn)), _part_of(int(fn6[0])) if _a is None or _b is None: continue _kk = (min(_a, _b), max(_a, _b)) _dg['tie_pair'][_kk] = _dg['tie_pair'].get(_kk, 0) + 1 _flo_i = [i for i, pp in enumerate(parts) if pp['part_no'] in set(int(t[0]) for t in floating)] _lo = np.array([pp['lo'] for pp in parts], float) _hi = np.array([pp['hi'] for pp in parts], float) _need = set(_flo_i) for _a in _flo_i: _d = np.maximum(np.maximum(_lo - _hi[_a], _lo[_a] - _hi), 0.0).max(axis=1) _need |= set(int(x) for x in np.where(_d <= BRIDGE_MAX_GAP_MM)[0]) for _i in sorted(_need): _pp = parts[_i] _dg['mesh'][_i] = {'faces6': np.asarray(_pp['faces6'], np.int32), 'surf_ids': np.asarray(_pp['surf_ids'], np.int32), 'n0': int(_pp['n0']), 'n1': int(_pp['n1'])} with open(os.environ['DZW_FEA_DIAG'], 'wb') as _fh: pickle.dump(_dg, _fh, protocol=4) if _logf: _logf(' [진단] 연결 진단 덤프 저장 — 부품 %d / 떠있음 %d / %s' % (len(parts), len(floating), os.environ['DZW_FEA_DIAG'])) except Exception as _dex: if _logf: _logf(' [진단] 덤프 실패: %r' % (_dex,)) # ── 밀도 상향 한 바퀴: '간격' 이 아니라 '메시' 때문에 못 이은 부재만 ────── # ★ 왜 여기인가: 다리 잇기가 끝난 바로 이 지점에서 '왜 못 이었는지' 가 # 실측값으로 손에 있다. 사유가 '간격이 멀다' 면 손댈 게 없지만, # 'CAD 간격은 이을 범위 안인데 그 거리 안에 묶을 절점이 없다' 면 그것은 # 메시 밀도 문제다 — 그 부재만 다시 촘촘히 메시해 한 번 더 시도한다. # 부품번호를 코드에 박지 않는다. 판단 근거는 전부 이 실행의 실측값이다. if floating and _fine_rounds > 0 and part_sizes and quality_bulk: _floset = set(int(t[0]) for t in floating) _fineq = {} for r in (_bridge_rej or []): gv = r.get('gap_mm') if gv is None or gv > BRIDGE_MAX_GAP_MM: continue # 실제로 떨어진 쌍 — 메시 문제 아니다 _ngs = [x.get('node_gap_mm') for x in (r.get('roles') or []) if x.get('node_gap_mm') is not None] if not _ngs or min(_ngs) <= gv + BRIDGE_TOL_MARGIN_MM: continue # 절점은 이미 가깝다 — 다른 원인이다 for _no in (r['a'], r['b']): if int(_no) not in _floset: continue # 떠 있는 쪽만 손댄다 _need = float(gv) + BRIDGE_TOL_MARGIN_MM if _no not in _fineq or _need < _fineq[int(_no)]: _fineq[int(_no)] = _need if _fineq: _fp = dict(_fine_parts or {}) _fp.update(_fineq) if _fp != dict(_fine_parts or {}): shutil.rmtree(wd, ignore_errors=True) if _logf: _logf(' [밀도상향] 메시가 성겨 못 이은 부재 %s → 요소크기 %s 로 ' '다시 메시하고 한 번 더 잇는다' % (sorted(_fineq), ', '.join('%d번 %.2fmm' % (k, v) for k, v in sorted(_fineq.items())))) return run_fea(positions, load_N=load_N, material=material, E=E, nu=nu, fix_points=fix_points, load_points=load_points, pick_radius=pick_radius, keep=keep, force_obb=force_obb, budget_nodes=budget_nodes, weld_mode=weld_mode, weld_throat=weld_throat, part_sizes=part_sizes, plastic=plastic, part_names=part_names, _drop_rounds=_drop_rounds, _auto_dropped=_auto_dropped, _orig_nos=_orig_nos, quality_bulk=quality_bulk, weld_seed=weld_seed, submodel=submodel, keep_arrays=keep_arrays, local_refine_n=local_refine_n, _logf=_logf, _fine_parts=_fp, _fine_rounds=_fine_rounds - 1) if floating and _drop_rounds > 0 and part_sizes and len(floating) < len(parts) - 4: # ★ 전체해석(기계 1대)에서는 '떠 있는 부품' 을 자동으로 빼고 다시 푼다. # 근거(실측): 146부품 투입 → 21개가 공중에 뜸 → 자동 제외 → 125부품으로 완주. # 한 부품이라도 떠 있으면 강성행렬이 특이해져 가짜 변위가 나오므로, 운전자에게 # '빼고 다시 찍으세요' 라고 되돌려 보내는 것보다 빼고 풀어 주는 편이 낫다. # 제외된 부품은 결과에 목록으로 돌려줘 화면에서 회색으로 표시한다. _dropno = set(int(t[0]) for t in floating) _off, _keepP, _keepS, _keepN, _keepO = 0, [], [], [], [] for _pi, _cnt in enumerate(part_sizes): _rows = int(_cnt) // 3 if _onum(_pi) not in _dropno: _keepO.append(_onum(_pi)) _keepP.append(P[_off:_off + _rows]) _keepS.append(int(_cnt)) _keepN.append(_pname(part_names, _pi)) _off += _rows shutil.rmtree(wd, ignore_errors=True) if _keepP and len(_keepS) < len(part_sizes): # ★ 사유에 '실측 간격' 을 적는다. '이어지지 않았다' 만 적으면 운전자가 # CAD 에서 실제로 떨어진 부품인지, 프로그램이 못 찾은 것인지 구분할 수 없다. def _why_float(no): gv = _gap_of.get(int(no)) if gv is None: return ('고정된 부재와 용접으로 이어지지 않아 자동 제외 ' '— 실측 %.0fmm 안에 이을 수 있는 이웃 부재가 없습니다' % BRIDGE_MAX_GAP_MM) if gv[0] > BRIDGE_MAX_GAP_MM: return ('실제 분리 — 가장 가까운 %d번 부재와 실측 간격 %.2fmm ' '(%.0fmm 넘게 떨어져 있어 용접으로 잇지 않았습니다)' % (int(gv[1]), gv[0], BRIDGE_MAX_GAP_MM)) # 간격은 이을 범위 안인데 이어지지 않았다 = 프로그램이 못 이은 것이다. # '실제 분리' 로 적으면 거짓이 된다 — 사실대로 적는다. return ('고정된 부재와 이어지지 않아 자동 제외 — 가장 가까운 %d번 부재와 ' '실측 간격 %.2fmm (이을 범위 안인데 접합을 못 만들었습니다: %s)' % (int(gv[1]), gv[0], gv[2] or '사유 미기록')) _acc = list(_auto_dropped or []) + [{'no': int(t[0]), 'name': t[1], 'why': _why_float(t[0])} for t in floating] return run_fea(np.vstack(_keepP).ravel().tolist(), load_N=load_N, material=material, E=E, nu=nu, fix_points=fix_points, load_points=load_points, pick_radius=pick_radius, keep=keep, force_obb=force_obb, budget_nodes=budget_nodes, weld_mode=weld_mode, weld_throat=weld_throat, part_sizes=_keepS, plastic=plastic, part_names=_keepN, _drop_rounds=_drop_rounds - 1, _auto_dropped=_acc, _orig_nos=_keepO, quality_bulk=quality_bulk, weld_seed=weld_seed, submodel=submodel, keep_arrays=keep_arrays, local_refine_n=local_refine_n, _logf=_logf, _fine_parts=_fine_parts, _fine_rounds=_fine_rounds) if floating: shutil.rmtree(wd, ignore_errors=True) lst = ', '.join('%d번(%s)' % t for t in floating[:8]) if len(floating) > 8: lst += ' …외 %d개' % (len(floating) - 8) raise RuntimeError( ('%s 가 고정된 부재와 용접으로 이어져 있지 않습니다(공중에 뜬 상태). ' '그대로 풀면 변위가 수천만 mm 로 나오는 가짜 결과가 됩니다.%s ' '고정점을 그 덩어리에도 찍거나, 사이를 잇는 부품을 함께 선택하거나, ' '그 부품을 선택에서 빼주세요.%s') % (lst, ' 특히 하중이 그 덩어리에 걸려 있습니다.' if float_loaded else '', (' ※ ' + _failed_text(failed_parts)) if failed_parts else '')) if method == 'obb_forced': bc_info['warn'] = ('검증용 외곽상자 모드입니다(바깥 치수 상자). ' '부재별 고정·용접 점검은 하지 않았습니다.') inp = os.path.join(wd, 'job.inp') write_inp(inp, nodes, elems, etype, set(int(x) for x in fix), loads, E, nu, ties=tie_blocks, plastic_fy=(fy if plastic else None), prescribed=prescribed) t0 = time.time() pr = subprocess.run([CCX_EXE, '-i', 'job'], cwd=wd, capture_output=True, text=True, errors='replace', timeout=900) secs = time.time() - t0 frd = os.path.join(wd, 'job.frd') if not os.path.isfile(frd) or 'Job finished' not in (pr.stdout or ''): tail = (pr.stdout or '')[-1200:] + (pr.stderr or '')[-600:] shutil.rmtree(wd, ignore_errors=True) raise RuntimeError('CalculiX 해석 실패:\n' + tail) disp, stress, rfall, errpct = parse_frd(frd, len(nodes)) vm = von_mises(stress) dmag = np.linalg.norm(disp, axis=1) # ══ ERR%(솔버 응력오차 추정) 요약 — 공짜 지표(계획 A-6 / D-1-2) ═════════ # 읽는 법: 낮을수록 그 자리의 응력을 믿을 수 있다. 실측으로 메시 품질에 # 반응함을 확인했다(전체 -p 30.7% → 국부 품질 8.68% → 더 조인 국부 9.17%). _ev = errpct[np.isfinite(errpct)] err_info = {'available': bool(len(_ev))} if len(_ev): err_info.update({ 'median_pct': round(float(np.median(_ev)), 2), 'p90_pct': round(float(np.percentile(_ev, 90)), 2), 'max_pct': round(float(_ev.max()), 2), 'nodes': int(len(_ev)), }) # ══ 컴플라이언스 ΣF·u — 전역 적분량 수렴지표(계획 A-3 / D-1-3) ══════════ # 왜 최대변위가 아니라 이것인가: 최대변위는 절점 1개 값이라 요소 1개가 # 망가지면 흔들린다. ΣF·u 는 전체 적분량이라 둔감하다. # 단서: 단조증가는 '구속조건 고정' 일 때만 성립한다 — 용접 묶음이 단계마다 # 달라지면 이 값의 단조성으로 이산화를 판정할 수 없다(계획 F-9). compliance = 0.0 for _ni, _fv in (loads or {}).items(): compliance += float(np.dot(np.asarray(_fv, float), disp[int(_ni)])) # ★ 서브모델(경계 강제변위) 모드에서는 외력 절점이 없어 위 합이 0 이 된다. # 그때는 강제변위 절점에서 외부가 해 준 일(Σ반력·변위)을 쓴다 — 같은 뜻의 # 전역 적분량이라 수렴지표로 쓸 수 있다. for _ni in sorted(disp_node_set): compliance += float(np.dot(rfall[int(_ni)], disp[int(_ni)])) # ══ F-1 반력 검증 — Σ반력 + Σ하중 ≒ 0 ══════════════════════════════ # 왜 '고정 절점만' 더하는가: 전체 절점으로 더하면 솔버가 전역 평형을 강제하므로 # 언제나 0 이 나와 검증이 공허해진다. 지지점이 실제로 받는 힘을 봐야 한다. # ★ CalculiX 의 FORC 는 '절점 총 힘'이다 — 하중 절점에서는 가한 외력이 그대로 # 나온다(실측: 판재 1,000N 하중 절점 합 = +1,000N). 그래서 '그 밖의 절점' # 합계를 볼 때는 고정 절점과 하중·강제변위 절점을 모두 빼야 한다. # 그러고도 남는 힘이 있으면 MPC(용접) 같은 다른 구속이 하중을 흡수한 것이다. _fixarr = np.asarray(sorted(set(int(x) for x in fix)), dtype=int) rf_fix = np.sum(rfall[_fixarr], axis=0) if len(_fixarr) else np.zeros(3) _other = np.ones(len(nodes), dtype=bool) if len(_fixarr): _other[_fixarr] = False for _ii in list(loads.keys()) + list(disp_node_set): _other[int(_ii)] = False rf_other = np.sum(rfall[_other], axis=0) f_sum = np.sum(list(loads.values()), axis=0) if loads else np.zeros(3) f_mag = float(np.linalg.norm(f_sum)) react_mag = float(np.linalg.norm(rf_fix)) # ★ 진짜 평형 점검 = '모든 절점의 절점력 합'. 실측(기계 1대)으로 확정했다: # 모든 절점 합 = 7.45 N (하중 29,400N 의 0.025%) → 솔버는 평형이다. # 그런데 고정 절점만 더하면 29,544N 으로 1.05% 어긋난다. # 차이는 용접 MPC(*EQUATION) 종속절점에 1,729N 이 절점력으로 잡히기 때문이다. # → 그러므로 "Σ고정반력 ≠ Σ하중" 을 '평형이 깨졌다' 로 말하면 용접 조립품마다 # 거짓 경보가 나고 운전자가 경고를 무시하게 된다. 둘을 나눠 보고한다: # ① equilibrium_err_pct : 전체 절점력 합 / 하중 ← 이것이 1% 게이트 # ② mpc_share_pct : 지지점이 아닌 곳(용접 MPC)으로 빠진 힘의 비율 rf_all = np.sum(rfall, axis=0) if f_mag > 1e-9: equil_err = float(np.linalg.norm(rf_all) / f_mag * 100.0) react_err = float(np.linalg.norm(rf_fix + f_sum) / f_mag * 100.0) mpc_share = float(np.linalg.norm(rf_other) / f_mag * 100.0) else: equil_err = react_err = mpc_share = None # 강제변위 모드 — 판정 안 함 react_info = { 'reaction_N': react_mag, 'reaction_xyz': [float(x) for x in rf_fix], 'load_xyz': [float(x) for x in f_sum], # 지지점 반력이 하중과 얼마나 다른가(용접으로 샌 몫이 여기 들어온다) 'reaction_err_pct': react_err, # 솔버 평형 점검 — 이것이 통과 기준이다 'equilibrium_N': float(np.linalg.norm(rf_all)), 'equilibrium_err_pct': equil_err, 'reaction_ok': (equil_err is not None and equil_err <= 1.0), # ★ 진단용 원시값. '하중의 몇 %가 용접으로 갔다' 로 읽지 말 것 — # 실측에서 이 값이 하중의 114% 로 나왔다(29,400N 하중에 33,628N). # 용접 MPC 종속절점의 절점력은 개별로 크고 서로 상쇄되는 내부력이라, # 벡터합의 크기가 '빠져나간 힘' 을 뜻하지 않는다. 뜻을 확정하지 못했으므로 # 화면에 비율로 띄우지 않는다(확인 못 한 것을 숫자로 말하지 않는다). 'reaction_other_N': float(np.linalg.norm(rf_other)), 'other_raw_pct': mpc_share, 'fix_nodes': int(len(_fixarr)), } # ══ F-2 구조 판정 마스크 — 대표 숫자의 기준을 하나로 통일 ═══════════════ # 왜 서버에서 만드는가: 종전에는 화면(bcMask)이 '찍은 좌표 + 기본 반경' 으로 # 따로 계산해, 안전율은 전역 최대 기준인데 색칠·배너는 제외 기준이라 # 같은 결과에 두 기준이 섞였다(보고서 C-2). 마스크를 서버가 한 번만 만들어 # 내려보내면 안전율·Max 라벨·컬러바 상한이 구조적으로 같은 값에서 나온다. # # 제외 대상 셋 — 전부 '계산상 인공 첨두' 가 생기는 자리다: # ① 고정·하중 경계조건 반경 안(실제 사용 반경 × 1.5) + 실제 구속·하중 절점 # ② 용접 MPC(*EQUATION) 종속·주절점과 그에 닿은 요소 1겹(±1요소) # ③ 볼록껍질로 대체된 부품 전체 — 오목 단면이 속 찬 덩어리가 되어 # 강성이 과대하므로 그 부품의 응력은 구조 판정에 쓸 수 없다 excl = np.zeros(len(nodes), dtype=bool) EXCL_R_FACTOR = 1.5 # ① 경계조건 for (_p, _r) in bc_spheres: if _r > 0: excl[np.linalg.norm(nodes - np.asarray(_p, float), axis=1) <= _r * EXCL_R_FACTOR] = True bc_seed = np.zeros(len(nodes), dtype=bool) if len(_fixarr): bc_seed[_fixarr] = True for _ii in list(loads.keys()) + list(disp_node_set): bc_seed[int(_ii)] = True if not bc_spheres: # 자동 외팔보 등 '찍은 좌표' 가 없는 경로 — 실제 구속·하중 절점에서 # 요소 1겹을 넓혀 제외한다(반경을 지어내지 않는다). excl |= _grow_one_layer(bc_seed, elems, len(nodes)) excl |= bc_seed n_bc = int(excl.sum()) # ② 용접 MPC 경계 ±1 요소 mpc_seed = np.zeros(len(nodes), dtype=bool) for (_sn, _fn6, _w) in (tie_blocks or []): mpc_seed[int(_sn)] = True for _k in range(6): mpc_seed[int(_fn6[_k])] = True n_mpc_seed = int(mpc_seed.sum()) if n_mpc_seed: excl |= _grow_one_layer(mpc_seed, elems, len(nodes)) | mpc_seed n_mpc = int(excl.sum()) - n_bc # ③ 껍질 대체 부품 전체 _hullno = set(int(x) for x in (meshinfo.get('hull_parts') or [])) n_hull_parts = 0 if _hullno: for pp in parts: if int(pp.get('part_no', -1)) in _hullno: excl[pp['n0']:pp['n1']] = True n_hull_parts += 1 n_hull = int(excl.sum()) - n_bc - n_mpc kept = int((~excl).sum()) kept_pct = kept / float(len(nodes)) * 100.0 # ★ 조용히 '제외 안 함' 으로 되돌리지 않는다. 남은 절점이 너무 적으면 그 사실을 # 숫자로 말하고 struct_usable=False 로 내린다 — 되돌리면 화면이 '구조 기준' # 이라고 말하면서 실제로는 전역 최대를 보여주는 거짓말이 된다. struct_usable = kept >= max(50, len(nodes) * 0.02) if struct_usable: _vmk = np.where(excl, -1.0, vm) simax = int(np.argmax(_vmk)) else: simax = int(np.argmax(vm)) mask_info = { 'excl_total': int(excl.sum()), 'kept': kept, 'kept_pct': round(kept_pct, 2), 'by_bc': n_bc, 'by_mpc': n_mpc, 'by_hull': n_hull, 'mpc_seed_nodes': n_mpc_seed, 'hull_parts_n': n_hull_parts, 'radius_factor': EXCL_R_FACTOR, 'radii_mm': sorted(set(round(float(r), 2) for (_p, r) in bc_spheres)), 'usable': bool(struct_usable), } # 취약점 상위 3곳 — 서로 충분히 떨어진 지점만 고른다 # ★ F-2: 순위도 '구조 판정 마스크' 로 거른다. 머리글 숫자만 구조 기준으로 # 바꾸고 순위표를 전역 최대로 남겨 두면, 안전율과 '먼저 항복하는 순서' 가 # 또 서로 다른 기준이 되어 C-2 를 되살린다. _vmrank = np.where(excl, -1.0, vm) if struct_usable else vm order = np.argsort(_vmrank)[::-1] hots, minsep = [], max(bbox.max() * 0.08, 5.0) for i in order: if len(hots) >= 3: break if _vmrank[i] < 0: break if all(np.linalg.norm(nodes[i] - nodes[j]) > minsep for j in hots): hots.append(int(i)) if not hots: hots = [int(np.argmax(vm))] # imax = 구조 최대(대표 숫자) / ipeak = 전역 최대(계산상 첨두, 참고용) imax = int(simax) ipeak = int(np.argmax(vm)) idmax = int(np.argmax(dmag)) ratio = vm / fy # ── 용접부 검토 ── # ★ 전달력을 '접합면 응력 적분' 으로 구하지 않는다(원인 규명 결과). # 접합면은 용접토(re-entrant corner) 응력 특이점이라, 절점평균 응력이 메시에 따라 # 요동치고 평형을 만족하지 않는다. 실측: 접합면 면적을 480mm2 로 정확히 잡아도 # 전달력이 5,313N(참값 1,000N, +431%), 축력은 0이어야 하는데 -1,317N 이 나왔다. # → 용접이 전달하는 힘은 '용접 건너편 부재에 걸린 외력의 합' 이며 이는 정역학으로 정확하다. # (부재가 그 용접 하나로만 붙어 있을 때. 여러 곳에 붙었으면 부정정이라 미검증으로 표시) weld_out = [] if weld_list and weld_mode != 'none' and weld_throat: allow = mat['fy'] * WELD_TAU_FACTOR for w in weld_list: if 'face_idx' not in w: continue B = parts[w['ind_part']] f3 = B['faces6'][w['force_idx']][:, :3] if len(w.get('force_idx', [])) else None if f3 is None or not len(f3): continue Lw = interface_perimeter(nodes, f3) if Lw <= 0: continue ctr_w = nodes[np.unique(f3)].mean(axis=0) # 이 용접 건너편(하중이 걸린 쪽) 부재에 실린 외력 합 + 용접중심에 대한 모멘트 lo_n, hi_n = B['n0'], B['n1'] Fv = np.zeros(3) Mv = np.zeros(3) for nidx, vec in loads.items(): if lo_n <= nidx < hi_n: Fv += vec Mv += np.cross(nodes[nidx] - ctr_w, vec) n_conn = sum(1 for x in weld_list if x['ind_part'] == w['ind_part'] or x['i'] == w['ind_part']) if np.linalg.norm(Fv) < 1e-9: continue Aw = weld_throat * Lw # 용접 목단면적 Zw = weld_throat * Lw ** 2 / 6.0 # 선용접군 단면계수 tau = float(np.linalg.norm(Fv) / Aw) sig = float(np.linalg.norm(Mv) / Zw) if Zw > 0 else 0.0 comb = float(math.sqrt(sig ** 2 + 3.0 * tau ** 2)) ok_stat = (n_conn == 1) weld_out.append({ # ★ 자리번호(w['i']+1)가 아니라 원래 부품번호 — 실패 부품을 건너뛰면 밀린다 'label': '부재%d(%s)-부재%d(%s) 용접 (용접선 %.0fmm, 목두께 a=%.1fmm)' % (parts[w['i']]['part_no'], parts[w['i']]['part_name'], parts[w['j']]['part_no'], parts[w['j']]['part_name'], Lw, weld_throat), 'at': ctr_w.tolist(), 'force_N': float(np.linalg.norm(Fv)), 'force_xyz': [float(x) for x in Fv], 'moment_Nmm': float(np.linalg.norm(Mv)), 'iface_area_mm2': float(w.get('iface_area', 0) or 0), 'length_mm': Lw, 'throat_mm': weld_throat, 'tau_MPa': tau, 'sigma_MPa': sig, 'combined_MPa': comb, 'allow_MPa': float(allow), 'pct': float(comb / allow * 100.0), 'verified': bool(ok_stat), 'note': None if ok_stat else '참고값(미검증) — 이 부재가 여러 곳에서 붙어 있어 힘 배분이 정역학만으로 안 정해집니다.', }) # ★ 용접부 검토를 요청했는데 결과가 0건이면 조용히 비워 두지 않는다. # 원인(실측 2026-10-05): 접합면 선별(force_idx)이 '두 부재 면 중심끼리의 거리' 로 # 후보를 고른다. 실제 제품처럼 메시가 성기면 면 하나가 수백 mm 라서 면 중심끼리는 # 멀어 전부 탈락한다 → 0건. # ★ 2026-10-09 고쳤다(B-6): force_idx 선별을 점↔삼각형 정확거리 # (point_tri_dist_pairs)로 바꿨다. 그래도 0건이면 사유는 다른 데 있다. weld_note = None if weld_mode != 'none' and weld_throat and not weld_out: weld_note = ('용접부 검토가 0건입니다 — 맞닿은 접합면을 찾지 못했습니다. ' '응력·변위 결과는 유효하지만 용접목 응력은 이 화면에서 판정하지 마세요.') # ══ 접합면별 단면력 (용접 MPC 종속절점의 절점력 합) ════════════════════ # 국부 재해석에서 '잘린 접합면' 을 지나는 힘이 전체 해석과 맞는지 보려면 # 알짜 합(≈0)이 아니라 **접합면 하나하나의 힘**을 봐야 한다. # ★★ 1단계 고침 (2026-10-12): 종전에는 값이 '힘 벡터' 뿐이었고 다리(bridge) # 방정식이 통째로 빠져 있었다. 그래서 # · 국부 재해석과 비교할 때 **어느 부품 쪽 절점력인지**(부호 기준)를 몰라, # 종속 쪽이 서로 반대면 정의상 어긋남이 200% 로 찍혔다. # · 다리로 이은 쌍은 기준값이 아예 없는데 국부가 그 자리를 강제했다. # → {키: {'F','side_part_no','n','kind'}} 로 바꾸고 다리도 넣는다. _iface_rf = {} for _w in (weld_list or []): _ti = _w.get('tied_ids') if _ti is None or not len(_ti): continue _k = (parts[_w['i']]['part_no'], parts[_w['ind_part']]['part_no']) _iface_rf[_k] = {'F': np.sum(rfall[_ti], axis=0), 'side_part_no': _w.get('side_part_no'), 'tied_ids': np.asarray(_ti, np.int64), 'n': int(len(_ti)), 'kind': '접합면'} for _b3 in (locals().get('_bridges') or []): _ti = _b3.get('tied_ids') if _ti is None or not len(_ti): continue _k = (int(_b3['a']), int(_b3['b'])) if _k in _iface_rf or (_k[1], _k[0]) in _iface_rf: continue _iface_rf[_k] = {'F': np.sum(rfall[_ti], axis=0), 'side_part_no': _b3.get('side_part_no'), 'tied_ids': np.asarray(_ti, np.int64), 'n': int(len(_ti)), 'kind': '다리'} # ══ 취약 후보 자동 탐지 (CEO 요구 3 / 계획 D-2) ═══════════════════════ # ★ 사용자는 위치를 입력하지 않는다. 전체 구조에서 지표 합의로 후보를 뽑는다. try: _hreal = (meshinfo.get('elem_size_real_mm') or meshinfo.get('elem_size_mm') or 5.0) cand_main, cand_sing, cand_sum = detect_candidates( nodes, elems, vm, stress, errpct, excl, _sjarr, E, nu, fy, parts, bc_spheres, mpc_seed, weld_out, float(_hreal)) except Exception as _cex: cand_main, cand_sing = [], [] cand_sum = {'error': '후보 탐지 실패: %s: %s' % (type(_cex).__name__, _cex)} res = { 'ok': True, 'method': method, 'method_kr': ('실제 형상 사면체(C3D10)' if method == 'tet_c3d10' else '외곽 상자 근사(검증용 강제 모드)'), 'warn': warn, # ★ 빠진 부품을 반드시 결과에 남긴다 — 조용히 빼면 '전체를 풀었다'는 착각을 준다. 'failed_parts': failed_parts, 'failed_text': _failed_text(failed_parts), # ★ 해석에서 빠진 부품 목록(형상 실패 + 자동 제외). 화면에서 반투명 회색으로 표시한다. 'excluded_parts': ([{'no': x['no'], 'name': x['name'], 'why': x['why']} for x in failed_parts] + list(_auto_dropped or [])), 'auto_dropped': list(_auto_dropped or []), 'material': {'code': material, 'name': mat['name'], 'E_MPa': E, 'nu': nu, 'fy_MPa': fy}, 'mesh': {'nodes': int(len(nodes)), 'elements': int(len(elems)), 'type': etype, **meshinfo}, 'bc': bc_info, 'applied_N': float(np.linalg.norm(np.sum(list(loads.values()), axis=0))) if loads else 0.0, # 요청 하중 총합 / 절점을 못 찾아 사라진 하중 — 화면에 "요청 ○○N 중 ○○N만 실림"을 띄우기 위함 'requested_N': float(req_sum), 'dropped_N': float(drop_sum), # ★ F-1 반력 검증 — 구조해석 1차 게이트. 1% 넘으면 화면에 빨간 경고. 'reaction': react_info, 'reaction_N': react_info['reaction_N'], 'reaction_err_pct': react_info['reaction_err_pct'], # 고정점이면서 하중점이기도 해서 하중에서 배제한 절점 수(배제분은 dropped_N 에 포함) 'overlap_nodes': int(overlap_nodes), 'obb': {'L_mm': size[0], 'b_mm': size[1], 't_mm': size[2]}, # ★ F-2 구조 판정 마스크 정보 — 무엇을 왜 몇 개 뺐는지 숫자로 남긴다. 'struct_mask': mask_info, # ★ 공짜 지표 3종(계획 D-1). 이 숫자가 없으면 '쓸 수 없는 메시' 에서 나온 # 응력을 판정으로 내보내게 된다. 'stress_error': err_info, # 솔버 자신의 응력오차 추정 ERR% 'compliance_Nmm': round(float(compliance), 3), # ΣF·u (변형에너지의 2배) # ★ CEO 요구 3 — 사용자 입력 없이 뽑은 취약 후보. 순위는 지표 합의이고 # 각 후보에 '근거' 문장과 '특이점 의심 / 구조 후보 / 메시 결함' 분류가 붙는다. 'candidates': cand_main, # 구속·껍질 자리(마스크로 제외된 영역)의 후보는 따로 담는다 — 용접 이음부는 # 피로 균열이 실제로 시작되는 자리이므로 숨기면 안 되고, 그렇다고 '위험' 으로 # 올려도 안 된다. 거리를 숫자로 적어 운전자가 직접 보게 한다(계획 D-2 금지항). 'candidates_singular': cand_sing, 'candidates_info': cand_sum, 'result': { # ★ F-2: 대표 숫자(최대 응력·안전율)는 '구조 최대' 기준이다. # = 고정·하중 반경 + 용접 MPC 경계 ±1요소 + 껍질 대체 부품 전부 제외. 'max_vonmises_MPa': float(vm[imax]), 'max_vonmises_at': nodes[imax].tolist(), 'yield_ratio_pct': float(ratio[imax] * 100.0), 'safety_factor': float(fy / vm[imax]) if vm[imax] > 0 else None, 'basis': 'struct' if struct_usable else 'global', # 전역 최대 = '계산상 첨두(참고)'. 메시를 잘게 쪼갤수록 커지며 수렴하지 않는다. 'peak_vonmises_MPa': float(vm[ipeak]), 'peak_vonmises_at': nodes[ipeak].tolist(), 'peak_yield_ratio_pct': float(ratio[ipeak] * 100.0), 'peak_safety_factor': float(fy / vm[ipeak]) if vm[ipeak] > 0 else None, 'max_disp_mm': float(dmag[idmax]), 'max_disp_at': nodes[idmax].tolist(), 'hotspots': [{'at': nodes[i].tolist(), 'MPa': float(vm[i]), 'pct': float(ratio[i] * 100.0)} for i in hots], }, 'welds': weld_out, 'weld_note': weld_note, 'weld_mode': weld_mode, 'solve_sec': round(secs, 2), 'submodel': sub_info, # 경계 강제변위 절점의 절점력 합 = 국부 영역이 밖에서 받는 알짜 힘(N). # 전체 해석에서 같은 부품들의 절점력 합과 맞아야 경계 강제가 타당하다. 'rf_driven_N': ([round(float(x), 2) for x in np.sum(rfall[sorted(disp_node_set)], axis=0)] if disp_node_set else None), # ★ 결함 3: 접합면별 경계 반력 합. 알짜 합은 평형 때문에 언제나 맞아서 # 검증이 공허하다 — 접합면 하나하나를 전체 해석과 맞춰 봐야 한다. 'rf_driven_by_iface_N': {('%d-%d' % k): [round(float(x), 2) for x in np.sum(rfall[sorted(set(v))], axis=0)] for k, v in _drv_src_nodes.items() if v}, # ★ 1단계: 그 힘이 '어느 부품 쪽 절점' 의 합인가. 전체 해석 쪽과 다르면 # 비교 전에 부호를 뒤집어야 한다(안 뒤집으면 정의상 200% 어긋남). 'rf_driven_by_iface_side': {('%d-%d' % k): int(max(v.items(), key=lambda x: x[1])[0]) for k, v in _drv_src_part.items() if v}, 'viz': { 'positions': nodes.round(4).reshape(-1).tolist(), 'ratio': (ratio * 100.0).round(2).tolist(), # 변형 애니메이션용 절점 변위(mm). 화면에서 과장배율을 곱해 쓴다. 'disp': disp.round(5).reshape(-1).tolist(), 'vonmises': vm.round(3).tolist(), # ★ F-2: 절점별 '구조 판정 제외' 표(1=제외). 화면이 따로 계산하지 않고 # 이것만 쓰게 해서 안전율·Max 라벨·컬러바 상한이 같은 기준에서 나오게 한다. 'excl': excl.astype(np.uint8).tolist(), # ★ 절점별 솔버 응력오차(%). 후보에 '신뢰도' 를 붙이는 근거(계획 D-2). # 값이 없는 절점(ccx 가 안 써 준 경우)은 -1 로 보낸다. 'errpct': np.where(np.isfinite(errpct), errpct, -1.0).round(1).tolist(), 'indices': surf.reshape(-1).astype(int).tolist(), }, } # ── 2·3단계: 상위 후보 자동 국부 재해석 + 수렴 판정 ────────────────── # 전체 해석이 끝난 뒤에만 돌린다(전체 변위장을 경계조건으로 쓰므로). if local_refine_n and cand_main and not submodel: _g = {'nodes': nodes, 'elems': elems, 'disp': disp, 'parts': [{'no': pp['part_no'], 'name': pp['part_name'], 'lo': pp['lo'], 'hi': pp['hi'], 'cadV': pp['cadV'], 'cadF': pp['cadF']} for pp in parts], 'seed_pairs': _seed_pairs_out, 'iface_rf': _iface_rf, 'h_global': float(_hreal)} _bk = {'material': material, 'E': E, 'nu': nu, 'load_N': load_N, 'pick_radius': pick_radius, 'weld_mode': weld_mode, 'weld_throat': weld_throat, 'plastic': plastic, 'load_points': load_points} try: _loc = local_refine(cand_main, _g, _bk, top_n=int(local_refine_n), logf=_logf) except Exception as _lex: _loc = [] res['local_error'] = '%s: %s' % (type(_lex).__name__, _lex) res['local'] = _loc # ★ 중복 제거(결함 4) 뒤에는 국부 결과의 순번이 후보 목록의 순번과 다르다. # 자리번호로 붙이면 엉뚱한 후보에 남의 판정이 붙는다 → 좌표로 붙인다. # 합쳐진 후보에게도 같은 판정을 붙이고 '합쳐짐' 을 표시한다. for _lr in _loc: _ats = [_lr['at']] + [m['at'] for m in (_lr.get('merged') or [])] for _cd in res['candidates']: if any(np.allclose(np.asarray(_cd['at'], float), np.asarray(_a, float), atol=1e-6) for _a in _ats): _cd['local'] = _lr _cd['merged_into'] = (None if np.allclose( np.asarray(_cd['at'], float), np.asarray(_lr['at'], float), atol=1e-6) else [round(float(x), 1) for x in _lr['at']]) if keep_arrays: # 국부 재해석(2단계)이 쓰는 날배열 — 같은 프로세스 안에서만 쓴다(JSON 금지) res['_arrays'] = {'nodes': nodes, 'elems': elems, 'disp': disp, 'rfall': rfall, 'vm': vm, 'excl': excl, 'parts': [{'no': pp['part_no'], 'name': pp['part_name'], 'n0': pp['n0'], 'n1': pp['n1'], 'lo': pp['lo'], 'hi': pp['hi'], 'cadV': pp['cadV'], 'cadF': pp['cadF']} for pp in parts], 'seed_pairs': _seed_pairs_out, 'iface_rf': _iface_rf, 'h_global': float(_hreal)} if keep: res['workdir'] = wd else: shutil.rmtree(wd, ignore_errors=True) return res # ══════════════════════════════════════════════════════════════════ # 8. 시험 형상 만들기 (검증용) — 삼각형 수프로 반환 # ══════════════════════════════════════════════════════════════════ def _quad(pts, a, b, c, d, out): for t in ((a, b, c), (a, c, d)): for k in t: out.extend(pts[k]) def make_plate(L, b, t): c = [[0, 0, 0], [L, 0, 0], [L, b, 0], [0, b, 0], [0, 0, t], [L, 0, t], [L, b, t], [0, b, t]] o = [] _quad(c, 0, 3, 2, 1, o); _quad(c, 4, 5, 6, 7, o) _quad(c, 0, 1, 5, 4, o); _quad(c, 1, 2, 6, 5, o) _quad(c, 2, 3, 7, 6, o); _quad(c, 3, 0, 4, 7, o) return o def _emit_quad(o, p0, p1, p2, p3): for tri in ((p0, p1, p2), (p0, p2, p3)): for p in tri: o.extend(p) def make_rhs_tube(L, H, B, t): """각파이프(직각 모서리, 속 빔). H=춤(하중 방향 z), B=폭(y), t=두께. 구성: 바깥 옆벽 4장 + 안쪽 옆벽 4장 + 양 끝 사각 링(도넛) — 끝면은 막지 않는다. ※ 바깥은 법선이 바깥, 안쪽은 법선이 안쪽(=구멍 쪽)을 향해야 닫힌 솔리드가 된다.""" hi, bi = H - 2.0 * t, B - 2.0 * t x0, x1 = 0.0, float(L) def rect(b, h): # (y,z) 네 모서리, 반시계 return [(-b / 2, -h / 2), (b / 2, -h / 2), (b / 2, h / 2), (-b / 2, h / 2)] O, I = rect(B, H), rect(bi, hi) o = [] for i in range(4): j = (i + 1) % 4 # 바깥 옆벽 (법선 바깥) _emit_quad(o, [x0, O[i][0], O[i][1]], [x0, O[j][0], O[j][1]], [x1, O[j][0], O[j][1]], [x1, O[i][0], O[i][1]]) # 안쪽 옆벽 (감김 반대 = 법선이 구멍 쪽) _emit_quad(o, [x1, I[i][0], I[i][1]], [x1, I[j][0], I[j][1]], [x0, I[j][0], I[j][1]], [x0, I[i][0], I[i][1]]) # 끝 링: x=L 은 법선 +x, x=0 은 법선 -x _emit_quad(o, [x1, O[i][0], O[i][1]], [x1, O[j][0], O[j][1]], [x1, I[j][0], I[j][1]], [x1, I[i][0], I[i][1]]) _emit_quad(o, [x0, O[i][0], O[i][1]], [x0, I[i][0], I[i][1]], [x0, I[j][0], I[j][1]], [x0, O[j][0], O[j][1]]) return o # ══════════════════════════════════════════════════════════════════ # 9. 자체 검증 # ══════════════════════════════════════════════════════════════════ def _section_stress(r, x_at, band): """길이축(x) 특정 단면 근처 절점들의 최대 von Mises. ★ 고정단 바로 옆은 구속에 의한 응력집중(수치 특이점)이라 단면 굽힘응력과 비교하면 안 된다. 그래서 고정단에서 충분히 떨어진 단면에서 비교한다.""" p = np.array(r['viz']['positions']).reshape(-1, 3) vm = np.array(r['viz']['vonmises']) sel = np.abs(p[:, 0] - x_at) <= band return float(vm[sel].max()) if sel.any() else float('nan') def _check(title, r, I, c, L, P, E, x_chk): """처짐(전체) + 단면 굽힘응력(고정단에서 떨어진 곳) 두 가지로 이론 대조.""" m = r['mesh'] th_d = P * L ** 3 / (3.0 * E * I) fd = r['result']['max_disp_mm'] ed = (fd - th_d) / th_d * 100.0 th_s = P * (L - x_chk) * c / I # M(x)=P(L-x), σ=Mc/I fsig = _section_stress(r, x_chk, max(L * 0.02, 8.0)) es = (fsig - th_s) / th_s * 100.0 print(' %s' % title) print(' 방식 : %s' % r['method_kr']) print(' 메시 : %s %d개 / 절점 %d개 (벽두께추정 %.2fmm, 벽당 %.2f요소, 요소 %.2fmm)' % (m['type'], m['elements'], m['nodes'], m.get('wall_mm', 0), m.get('through_wall', 0), m.get('elem_size_mm', 0))) print(' 단면2차 I : %.0f mm^4 (손계산)' % I) print(' 처짐 : FEM %.4f / 이론 %.4f mm → %+.2f%%' % (fd, th_d, ed)) print(' 굽힘응력 : FEM %.2f / 이론 %.2f MPa (x=%gmm 단면) → %+.2f%%' % (fsig, th_s, x_chk, es)) print(' 풀이시간 : %.1f초 (전체 %.1f초)' % (r['solve_sec'], r.get('_total', 0))) return abs(ed), abs(es) def _timed(**kw): t0 = time.time() r = run_fea(**kw) r['_total'] = round(time.time() - t0, 1) return r def selftest(): fails = [] P = 1000. # ── 검증 1: 속 찬 판재 ── print('\n[검증 1] 판재 외팔보 300 x 40 x 6 mm, 끝단 %gN (자동 외팔보 조건)' % P) L, b, t = 300., 40., 6. r = _timed(positions=make_plate(L, b, t), load_N=P, material='SS275') E = r['material']['E_MPa'] # 자동 조건은 '가장 얇은 축' 으로 민다 → 두께 t 방향 굽힘 ed, es = _check('판재(속 찬 단면)', r, b * t ** 3 / 12., t / 2, L, P, E, L * 0.5) if ed > 10 or r['method'] != 'tet_c3d10': fails.append('판재 처짐 %.1f%%' % ed) # ── 검증 2: 속 빈 각파이프 (2단계의 핵심) ── L2, H, B, t2 = 1000., 100., 50., 4.5 hi, bi = H - 2 * t2, B - 2 * t2 I_weak = (H * B ** 3 - hi * bi ** 3) / 12. # 자동조건이 미는 축(가장 얇은 축=B) I_strong = (B * H ** 3 - bi * hi ** 3) / 12. print('\n[검증 2] 각파이프(속 빔) %gx%gx%gT, 길이 %gmm, 끝단 %gN' % (H, B, t2, L2, P)) print(' ※ 1단계 외곽상자 방식이 가장 크게 틀리던 형상') print(' 손계산 I: 약축 %.0f / 강축 %.0f mm^4 (속 찬 상자였다면 %.0f)' % (I_weak, I_strong, B * H ** 3 / 12.)) r2 = _timed(positions=make_rhs_tube(L2, H, B, t2), load_N=P, material='STKR400') ed2, es2 = _check('각파이프(약축 굽힘)', r2, I_weak, B / 2, L2, P, r2['material']['E_MPa'], 400.) if ed2 > 10 or r2['method'] != 'tet_c3d10': fails.append('각파이프 처짐 %.1f%%' % ed2) # ── 1단계(외곽상자)와 직접 비교: 얼마나 위험하게 틀렸었나 ── print('\n[비교] 같은 각파이프를 1단계(외곽상자=속 찬 것으로 봄) 방식으로 풀면') r3 = _timed(positions=make_rhs_tube(L2, H, B, t2), load_N=P, material='STKR400', force_obb=True) th_hollow = P * L2 ** 3 / (3 * 205000. * I_weak) print(' 처짐 %.4f mm ←→ 실제 각파이프 이론 %.4f mm' % (r3['result']['max_disp_mm'], th_hollow)) print(' → 1단계는 이 부재를 %.1f배 튼튼한 것으로 잘못 봤다(응력 %.1f → 실제 %.1f MPa).' % (th_hollow / r3['result']['max_disp_mm'], r3['result']['max_vonmises_MPa'], r2['result']['max_vonmises_MPa'])) print('\n판정:', 'PASS' if not fails else 'FAIL ' + str(fails)) return 0 if not fails else 1 def survey(): """제품 OBJ 전수 조사 — 닫힌 솔리드 비율과 사면체 메시 성공률.""" import collections base = os.path.dirname(os.path.abspath(__file__)) files = sorted(glob.glob(os.path.join(base, 'fab_models', '**', '*.obj'), recursive=True)) print('조사 대상 %d개\n' % len(files)) stat = collections.Counter() bad = [] for fp in files: Vl, Fl = [], [] try: for ln in open(fp, errors='replace'): p = ln.split() if not p: continue if p[0] == 'v': Vl.append([float(p[1]), float(p[2]), float(p[3])]) elif p[0] == 'f': ix = [int(x.split('/')[0]) - 1 for x in p[1:]] for k in range(1, len(ix) - 1): Fl.append([ix[0], ix[k], ix[k + 1]]) except Exception as e: stat['읽기실패'] += 1 bad.append((os.path.basename(fp), 'read: %s' % e)) continue if len(Fl) < 4: stat['면부족'] += 1 bad.append((os.path.basename(fp), '면 %d개' % len(Fl))) continue V = np.array(Vl) soup = V[np.array(Fl, dtype=int)].reshape(-1) Vw, Fw = weld_mesh(soup) m = mesh_metrics(Vw, Fw) if m['closed']: stat['닫힘'] += 1 else: stat['열림'] += 1 bad.append((os.path.basename(fp), '열린엣지 %d개' % m['bad_edges'])) print('결과:', dict(stat)) if bad: print('\n문제 형상 (최대 15개):') for n, w in bad[:15]: print(' %-46s %s' % (n[:46], w)) return 0 # ══════════════════════════════════════════════════════════════════ # 10. HTTP 서버 # ══════════════════════════════════════════════════════════════════ def serve(): from http.server import BaseHTTPRequestHandler, ThreadingHTTPServer class H(BaseHTTPRequestHandler): def _cors(self): self.send_header('Access-Control-Allow-Origin', '*') self.send_header('Access-Control-Allow-Headers', 'Content-Type') self.send_header('Access-Control-Allow-Methods', 'POST, GET, OPTIONS') def do_OPTIONS(self): self.send_response(204); self._cors(); self.end_headers() def _json(self, code, obj): body = json.dumps(obj, ensure_ascii=False).encode('utf-8') self.send_response(code) self.send_header('Content-Type', 'application/json; charset=utf-8') self.send_header('Content-Length', str(len(body))) self._cors(); self.end_headers() self.wfile.write(body) def do_GET(self): if self.path.startswith('/api/fea/ping'): self._json(200, {'ok': True, 'ccx': os.path.isfile(CCX_EXE), 'tetgen': os.path.isfile(TETGEN_EXE), 'port': PORT}) elif self.path.startswith('/api/fea/materials'): self._json(200, {'ok': True, 'default': DEF_MAT, 'materials': MATERIALS}) else: self._json(404, {'ok': False, 'error': 'not found'}) def do_POST(self): if not self.path.startswith('/api/fea'): self._json(404, {'ok': False, 'error': 'not found'}); return try: n = int(self.headers.get('Content-Length') or 0) q = json.loads(self.rfile.read(n).decode('utf-8')) # ★ 마지막 요청을 그대로 저장한다(F-3 메시 수렴 시험용). # 화면과 '완전히 같은 모델·같은 경계조건' 으로 헤드리스 재실행을 # 하려면 형상·고정·하중을 손으로 다시 만들면 안 된다 — 조금만 달라도 # '같은 모델에서 메시만 바꾼 비교' 가 아니게 되어 수렴 판단이 무의미해진다. try: os.makedirs(WORK_ROOT, exist_ok=True) with open(os.path.join(WORK_ROOT, 'last_request.json'), 'w', encoding='utf-8') as _lf: json.dump(q, _lf) except Exception: pass res = run_fea(q.get('positions') or [], load_N=float(q.get('load_N', 1000.0)), material=q.get('material', DEF_MAT), E=q.get('E'), nu=q.get('nu'), fix_points=q.get('fix_points'), load_points=q.get('load_points'), pick_radius=q.get('pick_radius'), force_obb=bool(q.get('force_obb')), weld_mode=q.get('weld_mode') or 'bond', weld_throat=float(q.get('weld_throat') or 0.0), part_sizes=q.get('part_sizes'), part_names=q.get('part_names'), quality_bulk=bool(q.get('quality_bulk')), weld_seed=bool(q.get('weld_seed', True)), local_refine_n=int(q.get('local_refine') or 0)) self._json(200, res) except Exception as e: self._json(200, {'ok': False, 'error': '%s: %s' % (type(e).__name__, e)}) def log_message(self, *a): pass os.makedirs(WORK_ROOT, exist_ok=True) print('구조해석 서버 시작 → http://localhost:%d/api/fea' % PORT) print(' ccx :', CCX_EXE) print(' tetgen :', TETGEN_EXE) print(' 작업폴더(ASCII 전용):', WORK_ROOT) ThreadingHTTPServer(('127.0.0.1', PORT), H).serve_forever() if __name__ == '__main__': os.makedirs(WORK_ROOT, exist_ok=True) if '--selftest' in sys.argv: sys.exit(selftest()) if '--survey' in sys.argv: sys.exit(survey()) serve()