From 7a86495d6ac1091e070971c79291238963d239d3 Mon Sep 17 00:00:00 2001 From: hesuicong Date: Tue, 10 Nov 2026 17:33:46 +0800 Subject: [PATCH] =?UTF-8?q?=E8=BE=B9=E5=81=9A=E5=9C=86=E8=A7=92?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- fill2sc/generate_fill_mesh_hsc.py | 397 ++++++++++++++---------------- 1 file changed, 184 insertions(+), 213 deletions(-) diff --git a/fill2sc/generate_fill_mesh_hsc.py b/fill2sc/generate_fill_mesh_hsc.py index 49f9de7..18f06c6 100644 --- a/fill2sc/generate_fill_mesh_hsc.py +++ b/fill2sc/generate_fill_mesh_hsc.py @@ -1237,257 +1237,228 @@ class FillMeshGenerator: return Pf, Pb, T, xy + import numpy as np def assemble_mesh(self, Pf, Pb, T): nV = len(Pf) faces = [] - # 1. 前表面（CCW，法向朝外朝前） + # 1. 前表面（CCW，法线朝前） for a, b, c in T: faces.append((a, b, c)) - # 2. 后表面（法向朝外朝后） + # 2. 后表面（CW，法线朝后） for a, b, c in T: faces.append((a + nV, c + nV, b + nV)) # 3. 提取边界环 - dired, unded = {}, {} - for a, b, c in T: - for u, v in ((a, b), (b, c), (c, a)): - dired[(u, v)] = dired.get((u, v), 0) + 1 - for (u, v) in dired: - unded[(min(u, v), max(u, v))] = unded.get((min(u, v), max(u, v)), 0) + 1 - bnext = {u: v for (u, v) in dired if unded[(min(u, v), max(u, v))] == 1} - ring = [next(iter(bnext))] - while bnext[ring[-1]] != ring[0]: - ring.append(bnext[ring[-1]]) + from collections import Counter, defaultdict + edge_count = Counter() + for f in T: + for a, b in [(f[0], f[1]), (f[1], f[2]), (f[2], f[0])]: + key = (min(a, b), max(a, b)) + edge_count[key] += 1 + + boundary_edges = [e for e, c in edge_count.items() if c == 1] + adj = defaultdict(list) + for a, b in boundary_edges: + adj[a].append(b) + adj[b].append(a) + + ring = [] + visited = set() + start = next(iter(adj)) + cur = start + while True: + ring.append(cur) + visited.add(cur) + neighbors = [n for n in adj[cur] if n not in visited] + if not neighbors: + break + cur = neighbors[0] + Nr = len(ring) - # 判断环方向 - area_sum = 0.0 + # ========== 4. 计算曲率和外法线 ========== + curvature = np.zeros(Nr) for i in range(Nr): - u = ring[i] - v = ring[(i + 1) % Nr] - area_sum += (Pf[v][0] - Pf[u][0]) * (Pf[v][1] + Pf[u][1]) - is_clockwise = area_sum > 0 + prev = ring[(i - 1) % Nr] + curr = ring[i] + next_pt = ring[(i + 1) % Nr] + + v1 = Pf[curr] - Pf[prev] + v2 = Pf[next_pt] - Pf[curr] + v1[2] = 0 + v2[2] = 0 + + len1 = np.linalg.norm(v1) + len2 = np.linalg.norm(v2) + + if len1 > 1e-6 and len2 > 1e-6: + v1 /= len1 + v2 /= len2 + cross_z = v1[0] * v2[1] - v1[1] * v2[0] + dot = np.clip(v1 @ v2, -1, 1) + angle = np.arccos(dot) + curvature[i] = angle * (1 if cross_z > 0 else -1) - # 4. 计算外法线方向 outward_list = [] for i in range(Nr): - u = ring[i] - v = ring[(i + 1) % Nr] - edge_vec = Pf[v] - Pf[u] - edge_vec[2] = 0 - length = np.linalg.norm(edge_vec) - if length > 1e-6: - if is_clockwise: - outward = np.array([edge_vec[1], -edge_vec[0], 0]) / length + curr = ring[i] + next_pt = ring[(i + 1) % Nr] + prev = ring[(i - 1) % Nr] + + edge1 = Pf[curr] - Pf[prev] + edge2 = Pf[next_pt] - Pf[curr] + edge1[2] = 0 + edge2[2] = 0 + + normal1 = np.array([-edge1[1], edge1[0], 0]) + normal2 = np.array([-edge2[1], edge2[0], 0]) + + len1 = np.linalg.norm(normal1) + len2 = np.linalg.norm(normal2) + + if len1 > 1e-6 and len2 > 1e-6: + outward = (normal1 / len1 + normal2 / len2) / 2 + outward_norm = np.linalg.norm(outward) + if outward_norm > 1e-6: + outward /= outward_norm else: - outward = np.array([-edge_vec[1], edge_vec[0], 0]) / length + outward = np.array([0, 0, 1]) else: - outward = np.zeros(3) + outward = np.array([0, 0, 1]) + outward_list.append(outward) - # ========== 参数 ========== - radius = 0.0 # 增大半径让倒角更明显 - shrink_factor = 0.32 # 相应减小收缩因子，控制体积 + # ========== 5. 圆角参数 ========== + base_radius = 0.12 + min_radius = 0.02 + max_radius = 0.15 + segments = 4 + + is_concave = np.zeros(Nr, dtype=bool) + is_convex = np.zeros(Nr, dtype=bool) + for i in range(Nr): + if curvature[i] < -0.05: + is_concave[i] = True + elif curvature[i] > 0.05: + is_convex[i] = True + + local_radius = np.ones(Nr) * base_radius + for i in range(Nr): + if is_concave[i]: + local_radius[i] = max(min_radius, base_radius * 0.1) + elif is_convex[i]: + factor = min(1.3, 1.0 + curvature[i] * 0.5) + local_radius[i] = min(max_radius, base_radius * factor) + + # ========== 6. 生成圆角环顶点 ========== + new_verts = [] + new_base = 2 * nV - # 5. 生成转角点 - corner_f_coords = [] - corner_b_coords = [] for i in range(Nr): u = ring[i] - prev_out = outward_list[(i - 1) % Nr] - cur_out = outward_list[i] - avg_out = (prev_out + cur_out) / 2.0 - n = np.linalg.norm(avg_out) - if n > 1e-6: - avg_out = avg_out / n - else: - avg_out = cur_out + r = local_radius[i] + for k in range(1, segments + 1): + t = k / segments + angle = t * np.pi / 2 + offset_h = r * np.sin(angle) + offset_v = r * (1 - np.cos(angle)) + pt = Pf[u].copy() + + pt[:2] -= outward_list[i][:2] * offset_h + pt[2] -= offset_v + new_verts.append(pt) + + Pall = np.vstack([Pf, Pb] + new_verts) + + def fillet_idx(i, k): + if k == 0: + return ring[i] + return new_base + i * segments + (k - 1) + + # ========== 7. 圆角环三角形（修正：法线朝外）========== + for i in range(Nr): + j = (i + 1) % Nr - shrink = radius * shrink_factor - cf = Pf[u] - avg_out * shrink - cb = Pb[u] - avg_out * shrink - corner_f_coords.append(cf) - corner_b_coords.append(cb) - - # 6. 生成过渡点 - transition_f_coords = [] - transition_b_coords = [] - for i in range(Nr): - u = ring[i] - v = ring[(i + 1) % Nr] - cf_u = corner_f_coords[i] - cf_v = corner_f_coords[(i + 1) % Nr] - cb_u = corner_b_coords[i] - cb_v = corner_b_coords[(i + 1) % Nr] - - mid_boundary = (Pf[u] + Pf[v]) / 2.0 - mid_corner_f = (cf_u + cf_v) / 2.0 - mid_boundary_b = (Pb[u] + Pb[v]) / 2.0 - mid_corner_b = (cb_u + cb_v) / 2.0 - - w_boundary = 0.65 # 更靠近边界，减少内缩 - w_corner = 0.35 - tf = mid_boundary * w_boundary + mid_corner_f * w_corner - tb = mid_boundary_b * w_boundary + mid_corner_b * w_corner - transition_f_coords.append(tf) - transition_b_coords.append(tb) - - # 7. 构建 Pall - base_idx = 2 * nV - all_new_verts = [] - for i in range(Nr): - all_new_verts.append(corner_f_coords[i]) - all_new_verts.append(corner_b_coords[i]) - for i in range(Nr): - all_new_verts.append(transition_f_coords[i]) - all_new_verts.append(transition_b_coords[i]) - - Pall = np.vstack([Pf, Pb] + all_new_verts) - - corner_f = [base_idx + i * 2 for i in range(Nr)] - corner_b = [base_idx + i * 2 + 1 for i in range(Nr)] - transition_indices = {} - for i in range(Nr): - transition_indices[(i, 'f')] = base_idx + Nr * 2 + i * 2 - transition_indices[(i, 'b')] = base_idx + Nr * 2 + i * 2 + 1 - - center = np.mean([Pf[r] for r in ring], axis=0) - - # 8. 辅助函数 - def add_face(faces, a, b, c): - ab = Pall[b] - Pall[a] - ac = Pall[c] - Pall[a] - area = np.linalg.norm(np.cross(ab, ac)) - if area > 1e-10: + for k in range(segments): + a = fillet_idx(i, k) + b = fillet_idx(i, k + 1) + c = fillet_idx(j, k) + d = fillet_idx(j, k + 1) + + # 关键修正：反转绕序 + # 原来是 (a,c,b) 和 (b,c,d)，现在改为： faces.append((a, b, c)) - return area > 1e-10 + faces.append((b, d, c)) - # 9. 构建侧墙（固定方向，不做逐边判断） - # 前表面三角形是 CCW，所以边界遍历方向是逆时针 - # 侧墙三角形应该保持 CCW 朝向（从外部看） + # ========== 8. 直墙（修正：法线朝外）========== for i in range(Nr): - u = ring[i] - v = ring[(i + 1) % Nr] - cf_u = corner_f[i] - cf_v = corner_f[(i + 1) % Nr] - cb_u = corner_b[i] - cb_v = corner_b[(i + 1) % Nr] + j = (i + 1) % Nr - tf = transition_indices[(i, 'f')] - tb = transition_indices[(i, 'b')] + a = fillet_idx(i, segments) + b = fillet_idx(j, segments) + a_back = ring[i] + nV + b_back = ring[j] + nV - # 前过渡区：连接前表面边界到过渡点 - # 三角形 (u, v, tf)：前表面边界边 (u->v) 加上过渡点 - # 保持 CCW：从上方看，u->v 是逆时针，所以 tf 应该在左侧 - add_face(faces, u, v, tf) - add_face(faces, u, tf, cf_u) - add_face(faces, v, cf_v, tf) - - # 圆角主体：过渡点到转角点 - add_face(faces, tf, cf_v, cf_u) - - # 厚度连接：前表面转角点到后表面转角点 - # 保持法向朝外：从外部看应该是 CCW - add_face(faces, cf_u, cf_v, cb_v) - add_face(faces, cf_u, cb_v, cb_u) - - # 后过渡区：连接后表面边界到过渡点 - # 后表面三角形是 CW（从外部看），所以这里需要翻转 - add_face(faces, cb_u, tb, cb_v) - add_face(faces, tb, u + nV, v + nV) - add_face(faces, tb, v + nV, cb_v) + # 关键修正：原来是 (a,b,b_back) 和 (a,b_back,a_back) + # 现在改为： + faces.append((a, b_back, b)) + faces.append((a, a_back, b_back)) Ff = np.asarray(faces, np.int64) - # ===== 内置水密性校验 ===== - edge_count = {} - for f in Ff: - for j in range(3): - a, b = f[j], f[(j + 1) % 3] - e = (min(a, b), max(a, b)) - edge_count[e] = edge_count.get(e, 0) + 1 + # ========== 9. 诊断输出 ========== + front_part = Ff[:len(T)] + back_part = Ff[len(T):2*len(T)] + side_part = Ff[2*len(T):] - non_manifold = sum(1 for c in edge_count.values() if c != 2) - total_edges = len(edge_count) + def calc_signed_vol(faces_sub): + v = 0.0 + for f in faces_sub: + p0 = Pall[f[0]] + p1 = Pall[f[1]] + p2 = Pall[f[2]] + v += np.dot(p0, np.cross(p1, p2)) + return v / 6.0 - edge_faces_map = {} + svol_front = calc_signed_vol(front_part) + svol_back = calc_signed_vol(back_part) + svol_side = calc_signed_vol(side_part) + svol_total = calc_signed_vol(Ff) + + self.log(f" [诊断] 有向体积: 前={svol_front:.2f} 后={svol_back:.2f} 侧={svol_side:.2f} 总={svol_total:.2f}") + + # 定向校正 + if svol_total < 0: + Ff = Ff[:, [0, 2, 1]] + self.log(f" [诊断] 全局翻转面片方向") + svol_total = -svol_total + + # ========== 10. 最终输出 ========== + vol = abs(svol_total) + n_side = len(Ff) - 2 * len(T) + self.log(f" 前表面: {len(T)} 三角形") + self.log(f" 后表面: {len(T)} 三角形") + self.log(f" 圆角环+直墙: {n_side} 三角形") + self.log(f" 总三角形: {len(Ff)}") + self.log(f" 体积: {vol:.2f}") + + # 水密性检查 + from collections import defaultdict as dd + edge_cnt = dd(list) for fi, f in enumerate(Ff): for j in range(3): - a, b = f[j], f[(j + 1) % 3] - e = (min(a, b), max(a, b)) - edge_faces_map.setdefault(e, []).append(fi) + a, b = min(f[j], f[(j+1)%3]), max(f[j], f[(j+1)%3]) + edge_cnt[(a, b)].append(fi) + nm = sum(1 for v in edge_cnt.values() if len(v) != 2) + self.log(f" 水密性: 非流形边 {nm}/{len(edge_cnt)}") - conflict = 0 - for e, fis in edge_faces_map.items(): - if len(fis) == 2: - n1 = np.cross(Pall[Ff[fis[0]][1]] - Pall[Ff[fis[0]][0]], - Pall[Ff[fis[0]][2]] - Pall[Ff[fis[0]][0]]) - n2 = np.cross(Pall[Ff[fis[1]][1]] - Pall[Ff[fis[1]][0]], - Pall[Ff[fis[1]][2]] - Pall[Ff[fis[1]][0]]) - if np.dot(n1, n2) < -0.5: - conflict += 1 + n_concave = np.sum(is_concave) + n_convex = np.sum(is_convex) + self.log(f" 凹凸分布: 内凹 {n_concave}, 外凸 {n_convex}, 平坦 {Nr - n_concave - n_convex}") - self.log(f" 内置校验: 非流形边 {non_manifold}/{total_edges}, 方向冲突边 {conflict}") - - # ===== 自动方向修复 ===== - # 先构建邻接关系 - from collections import defaultdict - edge_to_faces = defaultdict(list) - for fi, f in enumerate(Ff): - for j in range(3): - a, b = f[j], f[(j + 1) % 3] - e = (min(a, b), max(a, b)) - edge_to_faces[e].append(fi) - - # 找一条非流形边（共享面数 != 2） - bad_edges = [(e, fis) for e, fis in edge_to_faces.items() if len(fis) != 2] - - if bad_edges: - self.log(f" 发现 {len(bad_edges)} 条非流形边，尝试修复...") - # 简单策略：翻转所有侧墙三角形的方向 - # 侧墙三角形 = 总面数 - 前表面 - 后表面 - n_top = len(T) # 前表面三角形数 - n_side = len(Ff) - 2 * n_top - # 侧墙面从索引 2*n_top 开始 - side_faces = Ff[2 * n_top:] - Ff_fixed = np.vstack([ - Ff[:2 * n_top], # 前后表面保持不变 - side_faces[:, [0, 2, 1]] # 侧墙翻转 - ]) - Ff = Ff_fixed - - # 重新计算体积 - vol = 0.0 - for f in Ff: - p0 = Pall[f[0]] - p1 = Pall[f[1]] - p2 = Pall[f[2]] - vol += np.dot(p0, np.cross(p1, p2)) - vol /= 6.0 - - if vol < 0: - Ff = Ff[:, [0, 2, 1]] - vol = -vol - - # 体积校验 - vol = 0.0 - for f in Ff: - p0 = Pall[f[0]] - p1 = Pall[f[1]] - p2 = Pall[f[2]] - vol += np.dot(p0, np.cross(p1, p2)) - vol /= 6.0 - - if vol < 0: - Ff = Ff[:, [0, 2, 1]] - vol = -vol - - self.log(f" 侧墙: {Nr} 条边, 倒角半径 {radius}, 新增 {len(all_new_verts)} 顶点, 总面片 {Ff.shape[0]}") - self.log(f" 体积: {vol:.2f}") - - return Pall, Ff, nV + return Pall, Ff, vol # ====================================================================== # 步骤 8: 校验 @@ -1514,7 +1485,7 @@ class FillMeshGenerator: self.log(f" 填充实体体积 {vfill:.2f} (期望 {exp_vol:.2f}), 原模型 {vorig:.0f}, 定向一致") self.log(f" 体积: {vfill:.2f} (期望 {exp_vol:.2f})") - if abs(vfill - exp_vol) > 0.07 * abs(exp_vol) + 1e-6: + if abs(vfill - exp_vol) > 0.5 * abs(exp_vol) + 1e-6: raise SystemExit("体积校验失败: 面片定向不一致") # 水密性