""" 구조해석 서버(fea_server.py) — 포트 8097 ========================================= 하는 일: 제관웹(fab.html)에서 보낸 부재 메시 + 하중/구속 조건을 받아 CalculiX(ccx) 로 구조해석을 돌리고, 결과(변위·von Mises 응력·항복 대비 %)를 JSON 으로 돌려준다. [2단계 — 실제 형상 해석] 1단계의 '외곽 상자 근사' 를 버리고, 보내온 삼각형 메시를 그대로 씁니다. ① 정점 병합(weld) → ② 닫힌 솔리드인지 검사 → ③ tetgen 으로 사면체 채우기 → ④ 2차 사면체(C3D10)로 승격 → ⑤ ccx 해석 · 속 빈 각파이프·채널·구멍 뚫린 판재가 실제 단면 그대로 계산됩니다. · 실패(안 닫힌 형상 등) 시에만 1단계 외곽상자 방식으로 물러나며, 결과 JSON 의 method 가 'obb_fallback' 이 되고 화면에 크게 경고가 뜹니다. ★ 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 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. 삼각형 메시 손질 — 정점 병합 · 닫힘 검사 · 두께 추정 # ══════════════════════════════════════════════════════════════════ 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]) def build_tet_mesh(V, F, wd, target_elems_through_wall=1.2, budget_nodes=None, quality_timeout=15.0): # ★ 벽당 요소 수 목표를 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: raise RuntimeError('가장 성긴 설정에서도 절점 %d개로 예산(%d) 초과 — 부재가 너무 복잡합니다' % (len(base['nodes']), int(budget_nodes))) 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 build_weld_equations(nodes, slave_ids, master_faces6, fixed_set, max_dist): """비적합(서로 안 맞는) 메시끼리 용접 접합 = 내가 직접 만드는 MPC(*EQUATION). ★ CalculiX *TIE 를 쓰지 않는 이유(실측): *TIE 는 종속 절점을 상대 면 위로 '실제로 옮겨서' 주변 사면체가 뒤집히고 'nonpositive jacobian' 으로 해석이 죽는다. *EQUATION 은 절점을 옮기지 않고 변위만 묶으므로 그 문제가 없다.""" from scipy.spatial import cKDTree if not len(slave_ids) or not len(master_faces6): return [], 0 ctr = nodes[master_faces6[:, :3]].mean(axis=1) tree = cKDTree(ctr) eqs, used = [], set() for s in slave_ids: s = int(s) if s in fixed_set or s in used: continue cand = tree.query_ball_point(nodes[s], max_dist) if not cand: continue best, bestd = None, 1e30 for fi in cand: P6 = nodes[master_faces6[fi]] w = tri6_weights(nodes[s], P6) if w is None: continue d = np.linalg.norm(nodes[s] - (w[:, None] * P6).sum(axis=0)) if d < bestd: bestd, best = d, (master_faces6[fi], w) if best is None or bestd > max_dist: continue used.add(s) eqs.append((s, best[0], best[1])) return eqs, len(used) def find_weld_pairs(parts, tol): """부재 쌍마다 '서로 맞닿은 곳'(용접선)을 찾는다. parts[i] = {'nodes','elems','n0','e0'} (n0/e0 = 전체 모델에서의 번호 시작점)""" from scipy.spatial import cKDTree welds = [] for i in range(len(parts)): for j in range(i + 1, len(parts)): A, B = parts[i], parts[j] # 겉면 절점만 비교(속 절점은 맞닿을 수 없다) ta = cKDTree(A['surf_pts']) tb = cKDTree(B['surf_pts']) pairs = ta.query_ball_tree(tb, tol) ai = [k for k, v in enumerate(pairs) if v] if len(ai) < 3: continue bj = sorted(set(x for v in pairs for x in v)) an = A['surf_ids'][np.array(ai)] bn = B['surf_ids'][np.array(bj)] zone = np.vstack([A['nodes'][an - A['n0']], B['nodes'][bn - B['n0']]]) ext = zone.max(axis=0) - zone.min(axis=0) welds.append({ 'i': i, 'j': j, 'dep_nodes': an, # A 의 접촉 절점(종속) 'ind_part': j, 'zone_min': zone.min(axis=0), 'zone_max': zone.max(axis=0), 'center': zone.mean(axis=0), 'length_mm': float(np.sort(ext)[-1]), # 용접선 길이 ≈ 접촉영역 최장변 }) return welds 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): 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', '*CLOAD'] else: o += ['*STEP', '*STATIC', '*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])) o += ['*NODE FILE', 'U, S', '*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): disp = np.zeros((nnode, 3)) stress = np.zeros((nnode, 6)) 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 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] else: n, v = _frd_fields(line, 6) if 1 <= n <= nnode: stress[n - 1] = v[:6] return disp, stress 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)) # ══════════════════════════════════════════════════════════════════ # 6. 경계조건 — 클릭 지점 → 절점 집합 # ══════════════════════════════════════════════════════════════════ 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] # ══════════════════════════════════════════════════════════════════ # 7. 해석 본체 # ══════════════════════════════════════════════════════════════════ 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): """positions: 평탄한 [x,y,z,...] 삼각형 수프. fix_points : [{'p':[x,y,z], 'r':반경}, ...] (없으면 자동 외팔보) load_points : [{'p':[x,y,z], 'r':반경, 'dir':[x,y,z], 'N':크기}, ...] """ 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: raise ValueError('고정점만 찍혔습니다 — 하중점도 찍어주세요.') 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 = [], [], [] try: if force_obb: raise RuntimeError('외곽상자 강제 지정') # ★ 부품 경계를 먼저 나눈 뒤, 부품마다 따로 정점을 병합한다. # 전체를 한꺼번에 병합하면 '맞닿은 부재끼리 정점이 붙어' 한 덩어리가 되고, # 그 덩어리는 속에 내부면이 낀 자기교차 형상이라 tetgen 이 거부한다 # (실측: 부재 5개 투입 → 성분 2개로 뭉침 → "A self-intersection was detected"). if part_sizes: comps = [] off = 0 for cnt in part_sizes: rows = int(cnt) // 3 # part_sizes 는 '좌표값 개수'(정점수x3)로 받는다 seg = P[off:off + rows] off += rows if len(seg) < 4: continue Vp, Fp = weld_mesh(seg) comps.extend(split_components(Vp, Fp)) V, F = weld_mesh(P) # 통계 표시용 else: V, F = weld_mesh(P) comps = split_components(V, F) if not comps: raise RuntimeError('삼각형 덩어리를 찾지 못했습니다') # 여러 부재면 절점 예산을 나눠 쓴다 _tot_budget = float(budget_nodes or MAX_C3D10_NODES) _vols = [max(mesh_metrics(vc, fc)['volume_mm3'], 1.0) for vc, fc in comps] _vsum = sum(_vols) # 부품이 많을수록 부품당 품질개선 시간을 짧게 — 전체가 늘어지지 않게 # 부재 1개 목표는 10초 이내 → 품질개선에 8초까지만 준다(검증용 각파이프는 2.8초라 여유). _qto = 8.0 if len(comps) == 1 else max(2.5, 40.0 / len(comps)) nlist, elist, n0 = [], [], 0 walls, esizes, coarse, flips, qual = [], [], False, 0, '' for ci, (Vc, Fc) in enumerate(comps): sub = os.path.join(wd, 'c%d' % ci) os.makedirs(sub, exist_ok=True) per = max(_tot_budget * _vols[ci] / _vsum, 12000) tm = build_tet_mesh(Vc, Fc, sub, budget_nodes=per, quality_timeout=_qto) 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, 'faces_n': _nn, 'faces_ctr': _ctr, 'n1': n0 + len(nd), 'surf_pts': nd[np.unique(tri)], 'wall': tm['metrics']['wall_mm'], }) 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) # 요소 시작번호 다시 계산 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: wtol = max(min(walls) * 1.2, 1.5) 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), 대신 법선이 상대 부재를 향하는지로 옆면을 걸러낸다. sel_f = np.where((dA <= wtol) & (facing < -0.7))[0] if len(sel_mpc): w['face_idx'] = sel_mpc w['force_idx'] = sel_f w['wtol'] = wtol meshinfo = { 'parts': len(parts), 'wall_mm': round(float(np.mean(walls)), 3), 'elem_size_mm': round(float(np.mean(esizes)), 3), 'through_wall': round(float(np.mean(walls) / np.mean(esizes)), 2), 'coarsened': coarse, 'flipped': flips, 'quality': qual, 'welded_verts': len(V), 'tris': len(F), 'welds': 0, } except Exception as ex: # 사면체 실패 → 1단계 외곽상자로 후퇴(크게 경고) method, etype = 'obb_fallback', 'C3D8I' warn = str(ex) 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 = {'fallback_reason': warn} if len(nodes) > (budget_nodes or MAX_C3D10_NODES) * 1.3: shutil.rmtree(wd, ignore_errors=True) raise RuntimeError('메시가 너무 큽니다(절점 %d개)' % len(nodes)) # ── 경계조건 ── 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 = {} bc_info = {'mode': 'click', 'fix_nodes': 0, 'load_nodes': 0, 'radius_mm': round(r_def, 2)} 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)])) for lp in (load_points or []): idx = nodes_near(nodes, lp['p'], lp.get('r') or r_def) if not len(idx): 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.]) per = float(lp.get('N', load_N)) / len(idx) for i in idx: loads[int(i)] = loads.get(int(i), np.zeros(3)) + d * per 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' if len(parts) > 1: # 부재가 여러 개인데 조건을 안 찍으면, 전체 외곽의 한쪽 끝만 고정된다. # 멀리 있는 부재는 용접 몇 점으로만 매달려 말이 안 되는 큰 변위가 나온다. bc_info['warn'] = ('부재 %d개인데 고정점을 안 찍었습니다. 전체 한쪽 끝만 자동 고정되어 ' '멀리 있는 부재가 매달린 상태로 계산됩니다 — 실제 취부점을 찍어주세요.' % len(parts)) if not len(fix): shutil.rmtree(wd, ignore_errors=True) raise RuntimeError('고정점 근처에 절점이 없습니다 — 반경을 키우거나 다른 곳을 찍어주세요') if not loads: shutil.rmtree(wd, ignore_errors=True) raise RuntimeError('하중점 근처에 절점이 없습니다 — 반경을 키우거나 다른 곳을 찍어주세요') bc_info['fix_nodes'] = int(len(fix)) bc_info['load_nodes'] = int(len(loads)) # ── 용접 접합 방정식 생성 ── # ★ 반드시 고정절점(fix)이 정해진 뒤에 만든다. 한 절점이 MPC 종속이면서 동시에 # 구속(SPC)이면 CalculiX 가 "detected on the dependent side of a MPC and a SPC" 로 거부한다. fixed_set = set(int(x) for x in fix) if weld_mode != 'none': for w in weld_list: if 'face_idx' not in w: continue B = parts[w['ind_part']] eqs, cnt = build_weld_equations(nodes, w['dep_nodes'], B['faces6'][w['face_idx']], fixed_set, w['wtol'] * 1.5) for e in eqs: fixed_set.add(e[0]) # 같은 절점을 두 번 종속시키지 않는다 tie_blocks.extend(eqs) w['tied_nodes'] = cnt 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: tied_cnt = {} for (sn, fn6, _w) in tie_blocks: for pi, pp in enumerate(parts): if pp['n0'] <= sn < pp['n1']: tied_cnt[pi] = tied_cnt.get(pi, 0) + 1 if pp['n0'] <= int(fn6[0]) < pp['n1']: 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 fixed_set) if not has_fix and tied_cnt.get(pi, 0) < 3: weak.append(pi + 1) if weak: bc_info['weak_parts'] = weak bc_info['warn'] = (('부재 %s 번이 고정도 안 되고 용접 이음도 거의 없습니다(사실상 공중에 뜬 상태). ' '결과 변위·응력을 믿지 마세요.') % ', '.join(map(str, weak))) 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)) 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 = parse_frd(frd, len(nodes)) vm = von_mises(stress) dmag = np.linalg.norm(disp, axis=1) # 취약점 상위 3곳 — 서로 충분히 떨어진 지점만 고른다 order = np.argsort(vm)[::-1] hots, minsep = [], max(bbox.max() * 0.08, 5.0) for i in order: if len(hots) >= 3: break if all(np.linalg.norm(nodes[i] - nodes[j]) > minsep for j in hots): hots.append(int(i)) imax, idmax = hots[0], 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({ 'label': '부재%d-부재%d 용접 (용접선 %.0fmm, 목두께 a=%.1fmm)' % (w['i'] + 1, w['j'] + 1, 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 '참고값(미검증) — 이 부재가 여러 곳에서 붙어 있어 힘 배분이 정역학만으로 안 정해집니다.', }) res = { 'ok': True, 'method': method, 'method_kr': ('실제 형상 사면체(C3D10)' if method == 'tet_c3d10' else '외곽 상자 근사 — 실제 형상 실패로 후퇴'), 'warn': warn, '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, 'obb': {'L_mm': size[0], 'b_mm': size[1], 't_mm': size[2]}, 'result': { '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, '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_mode': weld_mode, 'solve_sec': round(secs, 2), '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(), 'indices': surf.reshape(-1).astype(int).tolist(), }, } 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')) 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')) 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()