diff --git a/tools/adaptive-isogrid/src/brain/__main__.py b/tools/adaptive-isogrid/src/brain/__main__.py index 8e852e8d..dd82fafd 100644 --- a/tools/adaptive-isogrid/src/brain/__main__.py +++ b/tools/adaptive-isogrid/src/brain/__main__.py @@ -48,7 +48,39 @@ def _plot_boundary_polyline(geometry: Dict[str, Any], arc_pts: int = 64) -> np.n pts = typed_segments_to_polyline(typed, arc_pts=arc_pts) if len(pts) >= 3: return np.asarray(pts, dtype=float) - return np.asarray(geometry["outer_boundary"], dtype=float) + + outer = np.asarray(geometry["outer_boundary"], dtype=float) + if len(outer) < 4: + return outer + # v1 fallback: dense polyline boundaries may encode fillets. + # Resample uniformly for smoother plotting while preserving the polygon path. + if np.allclose(outer[0], outer[-1]): + ring = outer + else: + ring = np.vstack([outer, outer[0]]) + seg = np.diff(ring, axis=0) + seg_len = np.hypot(seg[:, 0], seg[:, 1]) + nonzero = seg_len[seg_len > 1e-9] + if len(nonzero) == 0: + return outer + step = max(float(np.median(nonzero) * 0.5), 0.5) + cum = np.r_[0.0, np.cumsum(seg_len)] + total = float(cum[-1]) + if total <= step: + return outer + samples = np.arange(0.0, total, step, dtype=float) + if samples[-1] < total: + samples = np.r_[samples, total] + out = [] + j = 0 + for s in samples: + while j < len(cum) - 2 and cum[j + 1] < s: + j += 1 + den = max(cum[j + 1] - cum[j], 1e-12) + t = (s - cum[j]) / den + p = ring[j] + t * (ring[j + 1] - ring[j]) + out.append([float(p[0]), float(p[1])]) + return np.asarray(out, dtype=float) def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Path, resolution: float = 3.0) -> None: @@ -62,7 +94,13 @@ def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Pa ax.plot(np.r_[outer[:, 0], outer[0, 0]], np.r_[outer[:, 1], outer[0, 1]], "w-", lw=1.5) for hole in geometry.get("holes", []): - hb = np.asarray(hole["boundary"]) + if hole.get("is_circular", False) and "center" in hole and "diameter" in hole: + cx, cy = hole["center"] + r = float(hole["diameter"]) * 0.5 + a = np.linspace(0.0, 2.0 * np.pi, 90, endpoint=True) + hb = np.column_stack([cx + r * np.cos(a), cy + r * np.sin(a)]) + else: + hb = np.asarray(hole["boundary"]) ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]], "k-", lw=0.9) ax.set_aspect("equal", adjustable="box") @@ -123,7 +161,13 @@ def _plot_triangulation(geometry: Dict[str, Any], triangulation: Dict[str, Any], # Blue hole circles for hole in geometry.get("holes", []): - hb = np.asarray(hole["boundary"]) + if hole.get("is_circular", False) and "center" in hole and "diameter" in hole: + cx, cy = hole["center"] + r = float(hole["diameter"]) * 0.5 + a = np.linspace(0.0, 2.0 * np.pi, 90, endpoint=True) + hb = np.column_stack([cx + r * np.cos(a), cy + r * np.sin(a)]) + else: + hb = np.asarray(hole["boundary"]) ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]], "b-", lw=0.8, zorder=3) @@ -182,7 +226,13 @@ def _plot_final_profile(geometry, pockets, ribbed_plate, out_path: Path, params: # Blue hole circles for hole in geometry.get("holes", []): - hb = np.asarray(hole["boundary"]) + if hole.get("is_circular", False) and "center" in hole and "diameter" in hole: + cx, cy = hole["center"] + r = float(hole["diameter"]) * 0.5 + a = np.linspace(0.0, 2.0 * np.pi, 90, endpoint=True) + hb = np.column_stack([cx + r * np.cos(a), cy + r * np.sin(a)]) + else: + hb = np.asarray(hole["boundary"]) ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]], "b-", lw=0.8, zorder=3) diff --git a/tools/adaptive-isogrid/src/brain/triangulation.py b/tools/adaptive-isogrid/src/brain/triangulation.py index 7e77c837..c2667aa7 100644 --- a/tools/adaptive-isogrid/src/brain/triangulation.py +++ b/tools/adaptive-isogrid/src/brain/triangulation.py @@ -14,68 +14,45 @@ a viable machining pocket. import numpy as np from scipy.spatial import Delaunay -from shapely.geometry import Polygon, Point, MultiPoint, LinearRing +from shapely.geometry import Polygon, Point, LinearRing +from shapely.geometry.base import BaseGeometry from shapely.ops import unary_union from src.shared.arc_utils import inset_arc, typed_segments_to_polyline, typed_segments_to_sparse from .density_field import evaluate_density, density_to_spacing -def _boundary_layer_offset_for_segment(mid_pt, geometry, params): - """Choose inward offset for boundary seed row.""" - explicit = params.get('boundary_layer_offset', None) - if explicit is not None: - return max(float(explicit), 0.0) - eta = evaluate_density(mid_pt[0], mid_pt[1], geometry, params) - return max(float(density_to_spacing(eta, params)), 1e-3) +def _geometry_to_list(geom: BaseGeometry) -> list[BaseGeometry]: + if geom.is_empty: + return [] + return [geom] if geom.geom_type == 'Polygon' else list(getattr(geom, 'geoms', [])) -def _add_boundary_layer_seed_points(points, geometry, params, plate_poly, keepout_union): - """Add a structured point row offset inward from each straight outer edge.""" +def _add_boundary_layer_seed_points(points, geometry, params, inner_plate, keepout_union): + """Add a row of points just inside the inset boundary to align boundary triangles.""" + offset = float(params.get('boundary_layer_offset', 0.0) or 0.0) + if offset <= 0.0: + offset = max(params.get('s_min', 10.0) * 0.6, 0.5) + + layer_geom = inner_plate.buffer(-offset) + if layer_geom.is_empty: + return points + boundary_pts = [] - ring = LinearRing(geometry['outer_boundary']) - is_ccw = bool(ring.is_ccw) - - # Prefer typed segments to avoid treating discretized arcs as straight edges - typed = geometry.get('outer_boundary_typed') - if typed: - segments = [seg for seg in typed if seg.get('type', 'line') == 'line'] - edge_pairs = [(_np(seg['start']), _np(seg['end'])) for seg in segments] - else: - coords = np.asarray(geometry['outer_boundary'], dtype=float) - if len(coords) >= 2 and np.allclose(coords[0], coords[-1]): - coords = coords[:-1] - edge_pairs = [] - for i in range(len(coords)): - edge_pairs.append((coords[i], coords[(i + 1) % len(coords)])) - - for a, b in edge_pairs: - dx, dy = b[0] - a[0], b[1] - a[1] - edge_len = float(np.hypot(dx, dy)) - if edge_len < 1e-9: + for poly in _geometry_to_list(layer_geom): + ring = poly.exterior + length = float(ring.length) + if length <= 1e-9: continue - - mid = np.array([(a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5], dtype=float) - spacing = float(density_to_spacing(evaluate_density(mid[0], mid[1], geometry, params), params)) - spacing = max(spacing, 1e-3) - offset = _boundary_layer_offset_for_segment(mid, geometry, params) - - nx_l, ny_l = (-dy / edge_len), (dx / edge_len) - nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l) - - n_pts = max(int(np.floor(edge_len / spacing)), 1) - for k in range(1, n_pts + 1): - t = k / (n_pts + 1) - bx = a[0] + t * dx - by = a[1] + t * dy - px = bx + offset * nx - py = by + offset * ny - p = Point(px, py) - if not plate_poly.buffer(1e-6).contains(p): + n_pts = max(int(np.ceil(length / max(params.get('s_min', 10.0), 1e-3))), 8) + for i in range(n_pts): + frac = i / n_pts + p = ring.interpolate(frac, normalized=True) + if not inner_plate.buffer(1e-6).covers(p): continue - if not keepout_union.is_empty and keepout_union.contains(p): + if not keepout_union.is_empty and keepout_union.buffer(1e-6).covers(p): continue - boundary_pts.append([px, py]) + boundary_pts.append([p.x, p.y]) if boundary_pts: return np.vstack([points, np.asarray(boundary_pts, dtype=np.float64)]) @@ -86,6 +63,40 @@ def _np(pt): return np.asarray([float(pt[0]), float(pt[1])], dtype=float) +def _hole_keepout_seed_points(geometry, params): + """Return sparse keepout seed points: 3 for circular holes, coarse ring for others.""" + d_keep = float(params['d_keep']) + pts = [] + for hole in geometry.get('holes', []): + if hole.get('is_circular', False) and 'center' in hole: + cx, cy = hole['center'] + diameter = float(hole.get('diameter', 10.0) or 10.0) + hole_radius = diameter / 2.0 + keepout_radius = hole_radius * (1.0 + d_keep) + for k in range(3): + a = (2.0 * np.pi * k) / 3.0 + pts.append([cx + keepout_radius * np.cos(a), cy + keepout_radius * np.sin(a)]) + continue + + hb = np.asarray(hole.get('boundary', []), dtype=float) + if len(hb) < 3: + continue + ring = LinearRing(hb) + keepout = Polygon(ring).buffer(max(d_keep, 0.0)) + if keepout.is_empty: + continue + for geom in _geometry_to_list(keepout): + ext = geom.exterior + n = max(int(np.ceil(ext.length / max(params.get('s_min', 10.0), 1e-3))), 6) + for i in range(n): + p = ext.interpolate(i / n, normalized=True) + pts.append([p.x, p.y]) + + if not pts: + return np.empty((0, 2), dtype=np.float64) + return np.asarray(pts, dtype=np.float64) + + def _generate_hex_grid(bbox, base_spacing): """ Generate a regular hexagonal-packed point grid. @@ -266,11 +277,7 @@ def _add_boundary_vertices(points, geometry, params, keepout_union): # Even spacing along inset boundary (as before) if not inner_plate.is_empty and inner_plate.is_valid: - if inner_plate.geom_type == 'MultiPolygon': - inner_polys = list(inner_plate.geoms) - else: - inner_polys = [inner_plate] - for inner_poly in inner_polys: + for inner_poly in _geometry_to_list(inner_plate): ext = inner_poly.exterior length = ext.length n_pts = max(int(length / s_min), 4) @@ -279,17 +286,16 @@ def _add_boundary_vertices(points, geometry, params, keepout_union): pt = ext.interpolate(frac, normalized=True) new_pts.append([pt.x, pt.y]) else: + # v1 geometry fallback: inset from polygon and force all inset ring vertices. inner_plate = plate_poly.buffer(-w_frame) if not inner_plate.is_empty: - if inner_plate.geom_type == 'MultiPolygon': - inner_polys = list(inner_plate.geoms) - else: - inner_polys = [inner_plate] - for inner_poly in inner_polys: + for inner_poly in _geometry_to_list(inner_plate): ext = inner_poly.exterior + # Always keep true inset corners/vertices from the offset polygon. coords = list(ext.coords)[:-1] for cx, cy in coords: new_pts.append([cx, cy]) + # Plus evenly spaced boundary points for perimeter alignment. length = ext.length n_pts = max(int(length / s_min), 4) for i in range(n_pts): @@ -297,19 +303,10 @@ def _add_boundary_vertices(points, geometry, params, keepout_union): pt = ext.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 - for cx, cy in list(ring.coords)[:-1]: - new_pts.append([cx, cy]) - 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]) + # Add sparse points for hole keepouts (3 points per circular keepout). + keepout_seeds = _hole_keepout_seed_points(geometry, params) + if len(keepout_seeds) > 0: + new_pts.extend(keepout_seeds.tolist()) if inner_plate.is_empty or not inner_plate.is_valid: inner_plate = plate_poly @@ -334,7 +331,7 @@ def generate_triangulation(geometry, params, max_refinement_passes=3): # Add structured boundary-layer seed row along straight edges plate_poly = Polygon(geometry['outer_boundary']) - all_pts = _add_boundary_layer_seed_points(all_pts, geometry, params, plate_poly, keepout_union) + all_pts = _add_boundary_layer_seed_points(all_pts, geometry, params, inner_plate, keepout_union) # Deduplicate close points all_pts = np.unique(np.round(all_pts, 4), axis=0) @@ -351,6 +348,10 @@ def generate_triangulation(geometry, params, max_refinement_passes=3): if inner_plate.is_empty: inner_plate = plate_poly + inner_valid_region = inner_plate + if not keepout_union.is_empty: + inner_valid_region = inner_plate.difference(keepout_union) + keep = [] for i, t in enumerate(triangles): cx = np.mean(all_pts[t, 0]) @@ -369,8 +370,15 @@ def generate_triangulation(geometry, params, max_refinement_passes=3): if not plate_poly.buffer(0.5).contains(Point(all_pts[vi])): all_inside = False break - if all_inside: - keep.append(i) + if not all_inside: + continue + + tri_poly = Polygon(all_pts[t]) + if tri_poly.is_empty or not tri_poly.is_valid: + continue + if not inner_valid_region.buffer(1e-6).covers(tri_poly): + continue + keep.append(i) triangles = triangles[keep] if keep else np.empty((0, 3), dtype=int)