""" 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")