From 906037f974629cb4f3d9e408a0accbb8becbd2c1 Mon Sep 17 00:00:00 2001 From: Antoine Date: Tue, 17 Feb 2026 21:48:55 +0000 Subject: [PATCH] feat(adaptive-isogrid): add Gmsh Frontal-Delaunay triangulation - Replaces scipy/Triangle iterative refinement with single-pass Gmsh - Separate distance fields for holes (I(x)) and edges (E(x)) - Frontal-Delaunay produces boundary-conforming, quasi-structured mesh - Better triangle quality for manufacturing (more equilateral) - Drop-in replacement: same signature as generate_triangulation() --- .../src/brain/triangulation_gmsh.py | 337 ++++++++++++++++++ 1 file changed, 337 insertions(+) create mode 100644 tools/adaptive-isogrid/src/brain/triangulation_gmsh.py diff --git a/tools/adaptive-isogrid/src/brain/triangulation_gmsh.py b/tools/adaptive-isogrid/src/brain/triangulation_gmsh.py new file mode 100644 index 00000000..0e763222 --- /dev/null +++ b/tools/adaptive-isogrid/src/brain/triangulation_gmsh.py @@ -0,0 +1,337 @@ +""" +Isogrid triangulation using Gmsh Frontal-Delaunay algorithm. + +Replaces scipy/Triangle-based approach with Gmsh's more structured meshing: +- Frontal-Delaunay advances from boundaries inward → natural boundary conformance +- Background size field handles density variation in ONE pass (no iterative refinement) +- Boolean hole handling → cleaner than PSLG workarounds +- More equilateral triangles → better rib widths + pocket shapes for machining + +Reference: Gmsh documentation, https://gmsh.info/doc/texinfo/gmsh.html +""" + +import numpy as np +import gmsh +from shapely.geometry import Polygon, Point, LinearRing +from shapely.ops import unary_union + +from .density_field import evaluate_density, density_to_spacing + + +def _ensure_gmsh_initialized(): + """Initialize Gmsh if not already done.""" + try: + gmsh.isInitialized() + except: + gmsh.initialize() + else: + if not gmsh.isInitialized(): + gmsh.initialize() + + +def _add_polygon_to_gmsh(coords: np.ndarray, lc: float = 1.0) -> tuple: + """ + Add a polygon to Gmsh geometry. + + Args: + coords: Nx2 array of vertices (NOT closed — first != last) + lc: characteristic length (mesh size at vertices) + + Returns: + (point_tags, curve_loop_tag) + """ + geo = gmsh.model.geo + + # Add points + point_tags = [] + for x, y in coords: + tag = geo.addPoint(float(x), float(y), 0.0, lc) + point_tags.append(tag) + + # Add lines connecting points + line_tags = [] + n = len(point_tags) + for i in range(n): + line = geo.addLine(point_tags[i], point_tags[(i + 1) % n]) + line_tags.append(line) + + # Create curve loop + curve_loop = geo.addCurveLoop(line_tags) + + return point_tags, line_tags, curve_loop + + +def _add_circle_to_gmsh(cx: float, cy: float, radius: float, lc: float = 1.0) -> tuple: + """ + Add a circle to Gmsh geometry using 4 arcs. + + Returns: + (center_tag, curve_loop_tag, arc_tags) + """ + geo = gmsh.model.geo + + # Center point (not meshed, just for arc definition) + center = geo.addPoint(cx, cy, 0.0, lc) + + # 4 points on circle (cardinal directions) + p_e = geo.addPoint(cx + radius, cy, 0.0, lc) + p_n = geo.addPoint(cx, cy + radius, 0.0, lc) + p_w = geo.addPoint(cx - radius, cy, 0.0, lc) + p_s = geo.addPoint(cx, cy - radius, 0.0, lc) + + # 4 arcs (counter-clockwise) + arc1 = geo.addCircleArc(p_e, center, p_n) + arc2 = geo.addCircleArc(p_n, center, p_w) + arc3 = geo.addCircleArc(p_w, center, p_s) + arc4 = geo.addCircleArc(p_s, center, p_e) + + arc_tags = [arc1, arc2, arc3, arc4] + curve_loop = geo.addCurveLoop(arc_tags) + + return center, curve_loop, arc_tags + + +def generate_triangulation_gmsh(geometry, params): + """ + Generate isogrid triangulation using Gmsh Frontal-Delaunay. + + Args: + geometry: dict with 'outer_boundary', 'holes', 'typed_boundary' (optional) + params: dict with s_min, s_max, w_frame, d_keep, etc. + + Returns: + dict with 'vertices' (Nx2), 'triangles' (Mx3), 'inner_plate' (Shapely Polygon) + """ + s_min = float(params['s_min']) + s_max = float(params['s_max']) + w_frame = float(params.get('w_frame', 8.0)) + d_keep = float(params.get('d_keep', 0.5)) + + # Build inner plate (inset by w_frame) + outer_poly = Polygon(geometry['outer_boundary']) + inner_plate = outer_poly.buffer(-w_frame, resolution=16) + + if inner_plate.is_empty or not inner_plate.is_valid: + return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)} + + # Handle MultiPolygon (take largest piece) + if inner_plate.geom_type == 'MultiPolygon': + inner_plate = max(inner_plate.geoms, key=lambda p: p.area) + + # Initialize Gmsh + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", 0) # Suppress output + gmsh.model.add("isogrid") + + geo = gmsh.model.geo + field = gmsh.model.mesh.field + + try: + # --- Step 1: Add outer boundary --- + outer_coords = np.array(inner_plate.exterior.coords)[:-1] # Remove closing point + _, outer_lines, outer_loop = _add_polygon_to_gmsh(outer_coords, s_max) + + # Create plane surface + surface = geo.addPlaneSurface([outer_loop]) + + # Track curves separately for distance fields + hole_curves = [] # For hole distance field + hole_loops = [] + hole_centers = [] # For diagnostics + + # --- Step 2: Add hole keepouts --- + for hole in geometry.get('holes', []): + diameter = float(hole.get('diameter', 10.0) or 10.0) + + if hole.get('is_circular', False) and 'center' in hole: + cx, cy = hole['center'] + hole_radius = diameter / 2.0 + keepout_radius = hole_radius * (1.0 + d_keep) + + # Check if keepout intersects inner plate + keepout_poly = Point(cx, cy).buffer(keepout_radius) + if not inner_plate.intersects(keepout_poly): + continue + + # Add circle + center_tag, hole_loop, arc_tags = _add_circle_to_gmsh( + cx, cy, keepout_radius, s_min + ) + hole_loops.append(hole_loop) + hole_curves.extend(arc_tags) + hole_centers.append((cx, cy, keepout_radius)) + + else: + # Non-circular hole — use polygon boundary + buffer + hb = np.asarray(hole.get('boundary', []), dtype=float) + if len(hb) < 3: + continue + + hole_poly = Polygon(hb).buffer(d_keep * diameter / 2.0) + if hole_poly.is_empty or not inner_plate.intersects(hole_poly): + continue + + # Clip to inner plate + clipped = inner_plate.intersection(hole_poly) + if clipped.is_empty or clipped.area < 1.0: + continue + + if clipped.geom_type == 'MultiPolygon': + clipped = max(clipped.geoms, key=lambda p: p.area) + + hole_coords = np.array(clipped.exterior.coords)[:-1] + _, hole_lines, hole_loop = _add_polygon_to_gmsh(hole_coords, s_min) + hole_loops.append(hole_loop) + hole_curves.extend(hole_lines) + + # Store centroid for diagnostics + c = clipped.centroid + hole_centers.append((c.x, c.y, np.sqrt(clipped.area / np.pi))) + + # Subtract holes from surface + if hole_loops: + # Recreate surface with holes + geo.remove([(2, surface)]) + surface = geo.addPlaneSurface([outer_loop] + [-h for h in hole_loops]) + + geo.synchronize() + + # --- Step 3: Set up adaptive size field --- + # Separate distance fields for holes vs edges (maps to η(x) = η₀ + α·I(x) + β·E(x)) + + # Get density params with defaults + R_0 = float(params.get('R_0', 30.0)) # Hole influence radius + R_edge = float(params.get('R_edge', 20.0)) # Edge influence radius + + size_fields = [] + + # --- Hole influence field: I(x) --- + if hole_curves: + hole_dist = field.add("Distance") + field.setNumbers(hole_dist, "CurvesList", hole_curves) + field.setNumber(hole_dist, "Sampling", 100) + + hole_thresh = field.add("Threshold") + field.setNumber(hole_thresh, "InField", hole_dist) + field.setNumber(hole_thresh, "SizeMin", s_min) # Dense near holes + field.setNumber(hole_thresh, "SizeMax", s_max) # Sparse far away + field.setNumber(hole_thresh, "DistMin", R_0 * 0.3) # Influence starts + field.setNumber(hole_thresh, "DistMax", R_0 * 2.0) # Influence ends + size_fields.append(hole_thresh) + + # --- Edge/boundary influence field: E(x) --- + edge_dist = field.add("Distance") + field.setNumbers(edge_dist, "CurvesList", outer_lines) + field.setNumber(edge_dist, "Sampling", 100) + + edge_thresh = field.add("Threshold") + field.setNumber(edge_thresh, "InField", edge_dist) + field.setNumber(edge_thresh, "SizeMin", s_min) # Dense near edges + field.setNumber(edge_thresh, "SizeMax", s_max) # Sparse in interior + field.setNumber(edge_thresh, "DistMin", 0.0) # Starts at boundary + field.setNumber(edge_thresh, "DistMax", R_edge) # Influence ends + size_fields.append(edge_thresh) + + # --- Combine fields: take minimum size at each point --- + if len(size_fields) > 1: + min_field = field.add("Min") + field.setNumbers(min_field, "FieldsList", size_fields) + field.setAsBackgroundMesh(min_field) + else: + field.setAsBackgroundMesh(size_fields[0]) + + # --- Step 4: Meshing options --- + # Algorithm 6 = Frontal-Delaunay (2D) + gmsh.option.setNumber("Mesh.Algorithm", 6) + + # Quality settings + gmsh.option.setNumber("Mesh.MeshSizeMin", s_min * 0.8) + gmsh.option.setNumber("Mesh.MeshSizeMax", s_max * 1.2) + gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + + # Angle constraints for quality triangles + gmsh.option.setNumber("Mesh.MinimumCircleNodes", 12) + gmsh.option.setNumber("Mesh.MinimumCurveNodes", 3) + + # --- Step 5: Generate mesh --- + gmsh.model.mesh.generate(2) + + # Optimize for better triangle quality + gmsh.model.mesh.optimize("Laplace2D") + + # --- Step 6: Extract mesh --- + node_tags, node_coords, _ = gmsh.model.mesh.getNodes() + elem_types, elem_tags, elem_node_tags = gmsh.model.mesh.getElements(2) + + # Vertices: reshape to Nx3, take only x,y + vertices = np.array(node_coords).reshape(-1, 3)[:, :2] + + # Build node tag → index mapping (tags aren't necessarily 0-indexed) + tag_to_idx = {tag: i for i, tag in enumerate(node_tags)} + + # Triangles: type 2 = 3-node triangles + triangles = [] + for etype, etags, enodes in zip(elem_types, elem_tags, elem_node_tags): + if etype == 2: # Triangle + tri_nodes = np.array(enodes).reshape(-1, 3) + for tri in tri_nodes: + idx_tri = [tag_to_idx[t] for t in tri] + triangles.append(idx_tri) + + triangles = np.array(triangles, dtype=np.int32) if triangles else np.empty((0, 3), dtype=np.int32) + + finally: + gmsh.finalize() + + return { + 'vertices': vertices, + 'triangles': triangles, + 'inner_plate': inner_plate, + } + + +# Alias for drop-in replacement +generate_triangulation = generate_triangulation_gmsh + + +if __name__ == "__main__": + # Quick test with sample geometry + import matplotlib.pyplot as plt + + # Simple rectangular plate with 3 holes + geometry = { + 'outer_boundary': [ + [0, 0], [200, 0], [200, 150], [0, 150] + ], + 'holes': [ + {'center': [50, 75], 'diameter': 20, 'is_circular': True}, + {'center': [100, 75], 'diameter': 15, 'is_circular': True}, + {'center': [150, 75], 'diameter': 25, 'is_circular': True}, + ] + } + + params = { + 's_min': 8.0, + 's_max': 25.0, + 'w_frame': 6.0, + 'd_keep': 0.5, + } + + print("Generating mesh with Gmsh Frontal-Delaunay...") + result = generate_triangulation_gmsh(geometry, params) + + verts = result['vertices'] + tris = result['triangles'] + + print(f"Vertices: {len(verts)}") + print(f"Triangles: {len(tris)}") + + # Plot + fig, ax = plt.subplots(figsize=(12, 9)) + ax.triplot(verts[:, 0], verts[:, 1], tris, 'b-', linewidth=0.5) + ax.set_aspect('equal') + ax.set_title(f'Gmsh Frontal-Delaunay: {len(tris)} triangles') + plt.savefig('/tmp/gmsh_isogrid_test.png', dpi=150) + print("Saved to /tmp/gmsh_isogrid_test.png")