""" Isogrid triangulation — generates a proper isogrid rib pattern. Strategy: lay a regular equilateral triangle grid, then adapt it: 1. Generate hex-packed vertex grid at base spacing 2. Scale spacing locally based on density field (growing away from features) 3. Delaunay triangulate the point set 4. Clip to plate boundary, exclude hole keepouts 5. Result: well-proportioned triangles everywhere, denser near holes/edges This is NOT a mesh — it's an isogrid layout. Every triangle should be a viable machining pocket. """ import numpy as np from scipy.spatial import Delaunay from shapely.geometry import Polygon, Point, MultiPoint from shapely.ops import unary_union from .density_field import evaluate_density, density_to_spacing def _generate_hex_grid(bbox, base_spacing): """ Generate a regular hexagonal-packed point grid. This creates the classic isogrid vertex layout: rows of points offset by half-spacing, producing equilateral triangles when Delaunay triangulated. Returns ndarray of shape (N, 2). """ x_min, y_min, x_max, y_max = bbox # Padding pad = base_spacing x_min -= pad y_min -= pad x_max += pad y_max += pad row_height = base_spacing * np.sqrt(3) / 2.0 points = [] row = 0 y = y_min while y <= y_max: # Alternate rows offset by half spacing x_offset = (base_spacing / 2.0) if (row % 2) else 0.0 x = x_min + x_offset while x <= x_max: points.append([x, y]) x += base_spacing y += row_height row += 1 return np.array(points, dtype=np.float64) def _adaptive_hex_grid(geometry, params): """ Generate density-adaptive hex grid. Start with a coarse grid at s_max spacing, then iteratively add points in high-density regions until local spacing matches the density-driven target. """ outer = np.array(geometry['outer_boundary']) x_min, y_min = outer.min(axis=0) x_max, y_max = outer.max(axis=0) bbox = (x_min, y_min, x_max, y_max) s_min = params['s_min'] s_max = params['s_max'] plate_poly = Polygon(geometry['outer_boundary']) # Build hole keepout polygons d_keep = params['d_keep'] keepout_polys = [] 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 boss_radius = hole_radius + d_keep * hole_radius keepout = Point(cx, cy).buffer(boss_radius, resolution=16) keepout_polys.append(keepout) keepout_union = unary_union(keepout_polys) if keepout_polys else Polygon() # Valid region = plate minus keepouts valid_region = plate_poly.difference(keepout_union) # Inset valid region by frame width to keep grid points inside the frame w_frame = params.get('w_frame', 8.0) inner_region = plate_poly.buffer(-w_frame) if inner_region.is_empty: inner_region = plate_poly # Generate base grid at s_max (coarsest) base_pts = _generate_hex_grid(bbox, s_max) # Filter to plate interior (outside keepouts, inside frame) valid_pts = [] for pt in base_pts: p = Point(pt[0], pt[1]) if inner_region.contains(p) and not keepout_union.contains(p): valid_pts.append(pt) valid_pts = np.array(valid_pts, dtype=np.float64) if valid_pts else np.empty((0, 2)) if s_min >= s_max * 0.95: # Uniform mode — no refinement needed return valid_pts, valid_region, keepout_union # Adaptive refinement: add denser points near high-density areas # Generate a fine grid at s_min fine_pts = _generate_hex_grid(bbox, s_min) # For each fine point, check if local density warrants adding it extra_pts = [] for pt in fine_pts: p = Point(pt[0], pt[1]) if not inner_region.contains(p) or keepout_union.contains(p): continue # What spacing does the density field want here? eta = evaluate_density(pt[0], pt[1], geometry, params) target_spacing = density_to_spacing(eta, params) # If target spacing < s_max * 0.8, this is a refinement zone if target_spacing < s_max * 0.8: # Check distance to nearest existing point if len(valid_pts) > 0: dists = np.sqrt(np.sum((valid_pts - pt)**2, axis=1)) min_dist = np.min(dists) # Add point if nearest neighbor is farther than target spacing if min_dist > target_spacing * 0.85: extra_pts.append(pt) # Also add to valid_pts for future distance checks valid_pts = np.vstack([valid_pts, pt]) if extra_pts: all_pts = np.vstack([valid_pts, np.array(extra_pts)]) # Deduplicate (shouldn't be needed but safe) all_pts = np.unique(np.round(all_pts, 6), axis=0) else: all_pts = valid_pts return all_pts, valid_region, keepout_union def _add_boundary_vertices(points, geometry, params, keepout_union): """ Add vertices along the outer boundary and hole keepout boundaries. This ensures triangles conform to boundaries rather than just being clipped. Points are spaced at approximately s_min along boundaries. """ s_min = params['s_min'] w_frame = params.get('w_frame', 8.0) new_pts = list(points) # Add points along inset outer boundary (frame inner edge) plate_poly = Polygon(geometry['outer_boundary']) inner_frame = plate_poly.buffer(-w_frame) if not inner_frame.is_empty and inner_frame.geom_type == 'Polygon': ring = inner_frame.exterior length = ring.length n_pts = max(int(length / s_min), 4) for i in range(n_pts): frac = i / n_pts pt = ring.interpolate(frac, normalized=True) new_pts.append([pt.x, pt.y]) # Add points along hole keepout boundaries if not keepout_union.is_empty: geoms = [keepout_union] if keepout_union.geom_type == 'Polygon' else list(keepout_union.geoms) for geom in geoms: ring = geom.exterior length = ring.length n_pts = max(int(length / (s_min * 0.7)), 6) for i in range(n_pts): frac = i / n_pts pt = ring.interpolate(frac, normalized=True) new_pts.append([pt.x, pt.y]) return np.array(new_pts, dtype=np.float64) def generate_triangulation(geometry, params, max_refinement_passes=3): """ Generate isogrid triangulation: hex-packed grid + Delaunay + clip. Returns dict with 'vertices' (ndarray Nx2) and 'triangles' (ndarray Mx3). """ # Generate adaptive point grid grid_pts, valid_region, keepout_union = _adaptive_hex_grid(geometry, params) if len(grid_pts) < 3: return {'vertices': grid_pts, 'triangles': np.empty((0, 3), dtype=int)} # Add boundary-conforming vertices all_pts = _add_boundary_vertices(grid_pts, geometry, params, keepout_union) # Deduplicate close points all_pts = np.unique(np.round(all_pts, 4), axis=0) if len(all_pts) < 3: return {'vertices': all_pts, 'triangles': np.empty((0, 3), dtype=int)} # Delaunay triangulation of the point set tri = Delaunay(all_pts) triangles = tri.simplices # Filter: keep only triangles whose centroid is inside valid region plate_poly = Polygon(geometry['outer_boundary']) w_frame = params.get('w_frame', 8.0) inner_plate = plate_poly.buffer(-w_frame) if inner_plate.is_empty: inner_plate = plate_poly keep = [] for i, t in enumerate(triangles): cx = np.mean(all_pts[t, 0]) cy = np.mean(all_pts[t, 1]) centroid = Point(cx, cy) # Keep if centroid is inside plate frame and outside keepouts, # AND all 3 vertices are inside the plate boundary (no crossovers) if not inner_plate.contains(centroid): continue if keepout_union.contains(centroid): continue all_inside = True for vi in t: if not plate_poly.contains(Point(all_pts[vi])): # Allow small tolerance (vertex on boundary is OK) if not plate_poly.buffer(0.5).contains(Point(all_pts[vi])): all_inside = False break if all_inside: keep.append(i) triangles = triangles[keep] if keep else np.empty((0, 3), dtype=int) # Filter degenerate / tiny triangles min_area = params.get('min_triangle_area', 20.0) if len(triangles) > 0: v0 = all_pts[triangles[:, 0]] v1 = all_pts[triangles[:, 1]] v2 = all_pts[triangles[:, 2]] areas = 0.5 * np.abs( (v1[:, 0] - v0[:, 0]) * (v2[:, 1] - v0[:, 1]) - (v2[:, 0] - v0[:, 0]) * (v1[:, 1] - v0[:, 1]) ) triangles = triangles[areas >= min_area] return {'vertices': all_pts, 'triangles': triangles}