diff --git a/tools/adaptive-isogrid/src/atomizer_study.py b/tools/adaptive-isogrid/src/atomizer_study.py
index 3a348e6d..36922b17 100644
--- a/tools/adaptive-isogrid/src/atomizer_study.py
+++ b/tools/adaptive-isogrid/src/atomizer_study.py
@@ -37,8 +37,8 @@ DEFAULT_PARAMS = {
     'p': 2.0,
     'beta': 0.3,
     'R_edge': 15.0,
-    's_min': 45.0,
-    's_max': 55.0,
+    's_min': 20.0,
+    's_max': 80.0,
     't_min': 2.5,
     't_0': 3.0,
     'gamma': 1.0,
diff --git a/tools/adaptive-isogrid/src/brain/__main__.py b/tools/adaptive-isogrid/src/brain/__main__.py
index dd82fafd..1cf7ce6f 100644
--- a/tools/adaptive-isogrid/src/brain/__main__.py
+++ b/tools/adaptive-isogrid/src/brain/__main__.py
@@ -43,44 +43,21 @@ def _merge_params(geometry: Dict[str, Any], params_file: Path | None) -> Dict[st
 
 
 def _plot_boundary_polyline(geometry: Dict[str, Any], arc_pts: int = 64) -> np.ndarray:
+    """Get the sandbox boundary as a clean polyline for plotting.
+
+    For v2 typed segments: densify arcs, keep lines exact.
+    For v1: just use the raw outer_boundary vertices — no resampling needed
+    since these are the true polygon corners from NX.
+    """
     typed = geometry.get("outer_boundary_typed")
     if typed:
         pts = typed_segments_to_polyline(typed, arc_pts=arc_pts)
         if len(pts) >= 3:
             return np.asarray(pts, dtype=float)
 
+    # v1 fallback: use raw boundary vertices directly — they ARE the geometry.
     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)
+    return outer
 
 
 def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Path, resolution: float = 3.0) -> None:
@@ -114,65 +91,83 @@ def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Pa
 
 def _plot_triangulation(geometry: Dict[str, Any], triangulation: Dict[str, Any], out_path: Path, params: Dict[str, Any] = None) -> None:
     """
-    Plot the Delaunay triangulation clipped to the inner frame (w_frame inset).
-    Green boundary, blue rib edges, blue hole circles.
+    Plot the triangulation with:
+    - Green: original sandbox boundary
+    - Red dashed: inset contour (w_frame offset) — this is what the mesh fills
+    - Blue: triangle edges
+    - Blue circles: hole boundaries
+    - Orange dashed: hole keepout rings
     """
-    from shapely.geometry import Polygon as ShapelyPolygon, LineString, Point
+    from shapely.geometry import Polygon as ShapelyPolygon, Point as ShapelyPoint
 
     verts = triangulation["vertices"]
     tris = triangulation["triangles"]
 
     outer = _plot_boundary_polyline(geometry)
     plate_poly = ShapelyPolygon(outer)
-    w_frame = (params or {}).get("w_frame", 2.0)
+    w_frame = (params or {}).get("w_frame", 8.0)
+    d_keep = (params or {}).get("d_keep", 1.5)
     inner_plate = plate_poly.buffer(-w_frame)
-    if inner_plate.is_empty or not inner_plate.is_valid:
-        inner_plate = plate_poly
 
-    fig, ax = plt.subplots(figsize=(8, 6), dpi=160)
+    fig, ax = plt.subplots(figsize=(10, 8), dpi=160)
 
-    # Draw triangle edges clipped to inner plate
+    # Draw all triangle edges
     drawn_edges = set()
     for tri in tris:
-        centroid = Point(verts[tri].mean(axis=0))
-        if not inner_plate.contains(centroid):
-            continue
         for i in range(3):
             edge = tuple(sorted([tri[i], tri[(i + 1) % 3]]))
             if edge in drawn_edges:
                 continue
             drawn_edges.add(edge)
             p1, p2 = verts[edge[0]], verts[edge[1]]
-            line = LineString([p1, p2])
-            clipped = inner_plate.intersection(line)
-            if clipped.is_empty:
-                continue
-            if clipped.geom_type == "LineString":
-                cx, cy = clipped.xy
-                ax.plot(cx, cy, color="#1f77b4", lw=0.5, alpha=0.85)
-            elif clipped.geom_type == "MultiLineString":
-                for seg in clipped.geoms:
-                    cx, cy = seg.xy
-                    ax.plot(cx, cy, color="#1f77b4", lw=0.5, alpha=0.85)
+            ax.plot([p1[0], p2[0]], [p1[1], p2[1]],
+                    color="#1f77b4", lw=0.5, alpha=0.85)
 
-    # Green sandbox boundary
+    # Green: original sandbox boundary
     ax.plot(np.r_[outer[:, 0], outer[0, 0]], np.r_[outer[:, 1], outer[0, 1]],
-            color="#228833", lw=1.8, zorder=5)
+            color="#228833", lw=1.8, zorder=5, label="Sandbox boundary")
 
-    # Blue hole circles
+    # Red dashed: inset contour (w_frame)
+    if not inner_plate.is_empty:
+        polys = [inner_plate] if inner_plate.geom_type == 'Polygon' else list(inner_plate.geoms)
+        for poly in polys:
+            ix, iy = poly.exterior.xy
+            ax.plot(ix, iy, color="#cc3333", lw=1.2, ls="--", zorder=4,
+                    label=f"Inset contour (w_frame={w_frame}mm)")
+
+    # Blue hole circles + Orange keepout rings
     for hole in geometry.get("holes", []):
         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)
+            # Hole circle
             hb = np.column_stack([cx + r * np.cos(a), cy + r * np.sin(a)])
+            ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]],
+                    "b-", lw=0.8, zorder=3)
+            # Keepout ring
+            kr = r * (1.0 + d_keep)
+            kb = np.column_stack([cx + kr * np.cos(a), cy + kr * np.sin(a)])
+            ax.plot(np.r_[kb[:, 0], kb[0, 0]], np.r_[kb[:, 1], kb[0, 1]],
+                    color="#ff8800", lw=0.8, ls="--", zorder=3)
         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)
+            ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]],
+                    "b-", lw=0.8, zorder=3)
+
+    # De-duplicate legend entries
+    handles, labels = ax.get_legend_handles_labels()
+    seen = set()
+    unique_handles, unique_labels = [], []
+    for h, l in zip(handles, labels):
+        if l not in seen:
+            seen.add(l)
+            unique_handles.append(h)
+            unique_labels.append(l)
+    ax.legend(unique_handles, unique_labels, loc="upper right", fontsize=8)
 
     ax.set_aspect("equal", adjustable="box")
-    ax.set_title("Constrained Delaunay Triangulation / Rib Pattern")
+    ax.set_title(f"Isogrid Triangulation ({len(tris)} triangles, {len(drawn_edges)} edges)")
     ax.set_xlabel("x [mm]")
     ax.set_ylabel("y [mm]")
     fig.tight_layout()
diff --git a/tools/adaptive-isogrid/src/brain/triangulation.py b/tools/adaptive-isogrid/src/brain/triangulation.py
index c2667aa7..8601b706 100644
--- a/tools/adaptive-isogrid/src/brain/triangulation.py
+++ b/tools/adaptive-isogrid/src/brain/triangulation.py
@@ -1,397 +1,365 @@
 """
 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
+Uses Shewchuk's Triangle library (constrained Delaunay + quality refinement)
+for boundary-conforming, clean triangulations.
 
-This is NOT a mesh — it's an isogrid layout. Every triangle should be
-a viable machining pocket.
+Strategy:
+1. Build PSLG (Planar Straight Line Graph):
+   - Outer contour = inset boundary (w_frame offset inward)
+   - Holes = keepout circles (d_keep offset outward from holes)
+2. Triangulate with quality constraints:
+   - Minimum angle ~25° (no slivers)
+   - Maximum area varies by density heatmap (smaller near holes, larger away)
+3. Result: clean mesh that fills 100% of valid area, with density-graded
+   triangle sizes and edges parallel to boundary/hole contours.
 """
 
 import numpy as np
-from scipy.spatial import Delaunay
-from shapely.geometry import Polygon, Point, LinearRing
+import triangle as tr
+from shapely.geometry import Polygon, Point, LinearRing, MultiPolygon
 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 src.shared.arc_utils import inset_arc, typed_segments_to_polyline
 from .density_field import evaluate_density, density_to_spacing
 
 
-def _geometry_to_list(geom: BaseGeometry) -> list[BaseGeometry]:
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _geometry_to_list(geom: BaseGeometry) -> list:
     if geom.is_empty:
         return []
-    return [geom] if geom.geom_type == 'Polygon' else list(getattr(geom, 'geoms', []))
+    if geom.geom_type == 'Polygon':
+        return [geom]
+    if geom.geom_type == 'MultiPolygon':
+        return list(geom.geoms)
+    return list(getattr(geom, 'geoms', []))
 
 
-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 = []
-    for poly in _geometry_to_list(layer_geom):
-        ring = poly.exterior
-        length = float(ring.length)
-        if length <= 1e-9:
-            continue
-        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.buffer(1e-6).covers(p):
-                continue
-            boundary_pts.append([p.x, p.y])
-
-    if boundary_pts:
-        return np.vstack([points, np.asarray(boundary_pts, dtype=np.float64)])
-    return points
+def _ring_to_segments(coords: np.ndarray, start_idx: int):
+    """Convert a ring (Nx2 array, NOT closed) to vertex array + segment index pairs."""
+    n = len(coords)
+    segments = []
+    for i in range(n):
+        segments.append([start_idx + i, start_idx + (i + 1) % n])
+    return segments
 
 
-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:
+def _sample_ring(ring, spacing: float) -> np.ndarray:
+    """Sample points along a Shapely ring at given spacing, returning Nx2 array (not closed)."""
+    length = float(ring.length)
+    if length < 1e-9:
         return np.empty((0, 2), dtype=np.float64)
-    return np.asarray(pts, dtype=np.float64)
+    n = max(int(np.ceil(length / max(spacing, 1e-3))), 8)
+    pts = []
+    for i in range(n):
+        p = ring.interpolate(i / n, normalized=True)
+        pts.append([p.x, p.y])
+    return np.array(pts, dtype=np.float64)
 
 
-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
+# ---------------------------------------------------------------------------
+# Step 1 — Build the inset plate polygon (frame inner edge)
+# ---------------------------------------------------------------------------
 
-    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 inset boundary vertices (arc-aware) and keepout boundary points."""
-    s_min = params['s_min']
-    w_frame = params.get('w_frame', 8.0)
-
-    new_pts = list(points)
+def _build_inner_plate(geometry, params) -> Polygon:
+    """Offset sandbox boundary inward by w_frame."""
+    w_frame = float(params.get('w_frame', 8.0))
     plate_poly = Polygon(geometry['outer_boundary'])
 
     typed_segments = geometry.get('outer_boundary_typed')
     if typed_segments:
         ring = LinearRing(geometry['outer_boundary'])
         is_ccw = bool(ring.is_ccw)
-        inset_segments = []
 
+        inset_segments = []
         for seg in typed_segments:
             stype = seg.get('type', 'line')
             if stype == 'arc':
                 center_inside = plate_poly.contains(Point(seg['center']))
                 inset_segments.append(inset_arc({**seg, 'center_inside': center_inside}, w_frame))
-                continue
+            else:
+                x1, y1 = seg['start']
+                x2, y2 = seg['end']
+                dx, dy = (x2 - x1), (y2 - y1)
+                ln = np.hypot(dx, dy)
+                if ln < 1e-12:
+                    continue
+                nx_l, ny_l = (-dy / ln), (dx / ln)
+                nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l)
+                inset_segments.append({
+                    'type': 'line',
+                    'start': [x1 + w_frame * nx, y1 + w_frame * ny],
+                    'end': [x2 + w_frame * nx, y2 + w_frame * ny],
+                })
 
-            x1, y1 = seg['start']
-            x2, y2 = seg['end']
-            dx, dy = (x2 - x1), (y2 - y1)
-            ln = np.hypot(dx, dy)
-            if ln < 1e-12:
-                continue
-            nx_l, ny_l = (-dy / ln), (dx / ln)
-            nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l)
-            inset_segments.append({
-                'type': 'line',
-                'start': [x1 + w_frame * nx, y1 + w_frame * ny],
-                'end': [x2 + w_frame * nx, y2 + w_frame * ny],
-            })
-
-        # Sparse points forced into Delaunay (3/arc, 2/line)
-        sparse_pts = typed_segments_to_sparse(inset_segments)
-        new_pts.extend(sparse_pts)
-
-        # Dense inset polygon for containment checks
-        dense = typed_segments_to_polyline(inset_segments, arc_pts=16)
+        dense = typed_segments_to_polyline(inset_segments, arc_pts=32)
         if len(dense) >= 3:
             inner_plate = Polygon(dense)
             if not inner_plate.is_valid:
                 inner_plate = inner_plate.buffer(0)
-            if inner_plate.is_empty:
-                inner_plate = plate_poly.buffer(-w_frame)
-        else:
-            inner_plate = plate_poly.buffer(-w_frame)
-
-        # Even spacing along inset boundary (as before)
-        if not inner_plate.is_empty and inner_plate.is_valid:
-            for inner_poly in _geometry_to_list(inner_plate):
-                ext = inner_poly.exterior
-                length = ext.length
-                n_pts = max(int(length / s_min), 4)
-                for i in range(n_pts):
-                    frac = i / n_pts
-                    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:
-            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):
-                    frac = i / n_pts
-                    pt = ext.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 not inner_plate.is_empty:
+                return inner_plate
 
+    inner_plate = plate_poly.buffer(-w_frame)
     if inner_plate.is_empty or not inner_plate.is_valid:
         inner_plate = plate_poly
+    return inner_plate
 
-    return np.array(new_pts, dtype=np.float64), inner_plate
 
+# ---------------------------------------------------------------------------
+# Step 2 — Build hole keepout polygons
+# ---------------------------------------------------------------------------
+
+def _build_keepouts(geometry, params) -> list:
+    """Build list of keepout Polygons (one per hole)."""
+    d_keep = float(params['d_keep'])
+    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 + d_keep * hole_radius
+            keepout = Point(cx, cy).buffer(keepout_radius, resolution=32)
+            keepouts.append(keepout)
+        else:
+            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 not keepout.is_empty:
+                keepouts.append(keepout)
+
+    return keepouts
+
+
+# ---------------------------------------------------------------------------
+# Step 3 — Build PSLG for Triangle library
+# ---------------------------------------------------------------------------
+
+def _build_pslg(inner_plate: Polygon, keepouts: list, boundary_spacing: float):
+    """
+    Build a Planar Straight Line Graph (PSLG) for the Triangle library.
+
+    Returns dict with:
+      'vertices': Nx2 array
+      'segments': Mx2 array of vertex index pairs
+      'holes': Hx2 array of hole marker points
+    """
+    vertices = []
+    segments = []
+    hole_markers = []
+
+    # Outer boundary of the inner plate
+    # Handle MultiPolygon (take largest piece)
+    polys = _geometry_to_list(inner_plate)
+    if not polys:
+        return None
+    main_poly = max(polys, key=lambda p: p.area)
+
+    # Sample the outer ring
+    outer_pts = _sample_ring(main_poly.exterior, boundary_spacing)
+    if len(outer_pts) < 3:
+        return None
+
+    start_idx = 0
+    vertices.extend(outer_pts.tolist())
+    segs = _ring_to_segments(outer_pts, start_idx)
+    segments.extend(segs)
+    start_idx += len(outer_pts)
+
+    # Handle any interior rings of the inner plate (shouldn't happen normally)
+    for interior_ring in main_poly.interiors:
+        ring_pts = _sample_ring(interior_ring, boundary_spacing)
+        if len(ring_pts) < 3:
+            continue
+        vertices.extend(ring_pts.tolist())
+        segs = _ring_to_segments(ring_pts, start_idx)
+        segments.extend(segs)
+        # Hole marker inside this interior ring
+        centroid = Polygon(interior_ring).centroid
+        hole_markers.append([centroid.x, centroid.y])
+        start_idx += len(ring_pts)
+
+    # Keepout holes — clip to inner plate
+    for keepout in keepouts:
+        clipped = inner_plate.intersection(keepout)
+        if clipped.is_empty:
+            continue
+        for kpoly in _geometry_to_list(clipped):
+            if kpoly.area < 1.0:
+                continue
+            ring_pts = _sample_ring(kpoly.exterior, boundary_spacing)
+            if len(ring_pts) < 3:
+                continue
+            vertices.extend(ring_pts.tolist())
+            segs = _ring_to_segments(ring_pts, start_idx)
+            segments.extend(segs)
+            # Hole marker inside the keepout
+            centroid = kpoly.centroid
+            hole_markers.append([centroid.x, centroid.y])
+            start_idx += len(ring_pts)
+
+    return {
+        'vertices': np.array(vertices, dtype=np.float64),
+        'segments': np.array(segments, dtype=np.int32),
+        'holes': np.array(hole_markers, dtype=np.float64) if hole_markers else np.empty((0, 2)),
+    }
+
+
+# ---------------------------------------------------------------------------
+# Step 4 — Compute area constraint from density field
+# ---------------------------------------------------------------------------
+
+def _max_area_for_spacing(spacing: float) -> float:
+    """Area of an equilateral triangle with given edge length."""
+    return (np.sqrt(3) / 4.0) * spacing ** 2
+
+
+def _refine_with_density(tri_result, geometry, params):
+    """
+    Iterative refinement: check each triangle's area against the density-based
+    max area at its centroid. If too large, re-triangulate with tighter constraint.
+
+    Uses Triangle's region-based area constraints via the 'a' switch.
+    """
+    s_min = float(params['s_min'])
+    s_max = float(params['s_max'])
+
+    verts = tri_result['vertices']
+    tris = tri_result['triangles']
+
+    # Compute max allowed area for each triangle based on density at centroid
+    max_areas = np.empty(len(tris))
+    for i, t in enumerate(tris):
+        cx = np.mean(verts[t, 0])
+        cy = np.mean(verts[t, 1])
+        eta = evaluate_density(cx, cy, geometry, params)
+        target_spacing = density_to_spacing(eta, params)
+        max_areas[i] = _max_area_for_spacing(target_spacing)
+
+    # Use the global maximum area as a single constraint
+    # (Triangle library doesn't support per-triangle area easily in basic mode)
+    # Instead, use a moderate constraint and let the quality refinement handle it
+    global_max_area = _max_area_for_spacing(s_max)
+    return global_max_area
+
+
+# ---------------------------------------------------------------------------
+# Main entry point
+# ---------------------------------------------------------------------------
 
 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 isogrid triangulation using Triangle library.
+
+    Returns dict with 'vertices' (Nx2) and 'triangles' (Mx3).
     """
-    # Generate adaptive point grid
-    grid_pts, valid_region, keepout_union = _adaptive_hex_grid(geometry, params)
+    s_min = float(params['s_min'])
+    s_max = float(params['s_max'])
 
-    if len(grid_pts) < 3:
-        return {'vertices': grid_pts, 'triangles': np.empty((0, 3), dtype=int)}
-
-    # Add boundary-conforming vertices and get inset plate polygon for clipping
-    all_pts, inner_plate = _add_boundary_vertices(grid_pts, geometry, params, keepout_union)
-
-    # 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, inner_plate, 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'])
+    # Step 1: Build inner plate
+    inner_plate = _build_inner_plate(geometry, params)
     if inner_plate.is_empty:
-        inner_plate = plate_poly
+        return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
 
-    inner_valid_region = inner_plate
-    if not keepout_union.is_empty:
-        inner_valid_region = inner_plate.difference(keepout_union)
+    # Step 2: Build keepouts
+    keepouts = _build_keepouts(geometry, params)
+    keepout_union = unary_union(keepouts) if keepouts else Polygon()
 
-    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 not all_inside:
-            continue
+    # Step 3: Build PSLG
+    # Boundary sampling at intermediate spacing for clean boundary conformance
+    boundary_spacing = max(s_min, min(s_max * 0.5, 30.0))
+    pslg = _build_pslg(inner_plate, keepouts, boundary_spacing)
 
-        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)
+    if pslg is None or len(pslg['vertices']) < 3:
+        return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
 
-    triangles = triangles[keep] if keep else np.empty((0, 3), dtype=int)
+    # Step 4: Initial triangulation with coarse area constraint
+    # 'p' = use PSLG (segments + holes)
+    # 'q25' = minimum angle 25° (no slivers)
+    # 'a{area}' = max triangle area (coarse — s_max)
+    max_area = _max_area_for_spacing(s_max)
+    tri_switches = f'pq25a{max_area:.1f}'
 
-    # 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]]
+    tri_input = {
+        'vertices': pslg['vertices'],
+        'segments': pslg['segments'],
+    }
+    if len(pslg['holes']) > 0:
+        tri_input['holes'] = pslg['holes']
+
+    try:
+        result = tr.triangulate(tri_input, tri_switches)
+    except Exception:
+        return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
+
+    verts = result.get('vertices', np.empty((0, 2)))
+    tris = result.get('triangles', np.empty((0, 3), dtype=int))
+
+    if len(tris) == 0:
+        return {'vertices': verts, 'triangles': tris}
+
+    # Step 5: Density-based refinement using per-triangle area constraints
+    # Triangle's 'r' switch refines an existing triangulation.
+    # We set 'triangle_max_area' per triangle based on density at centroid.
+    for pass_num in range(max_refinement_passes):
+        # Compute per-triangle max area from density field
+        per_tri_max = np.empty(len(tris))
+        needs_refinement = False
+
+        for i, t in enumerate(tris):
+            cx = np.mean(verts[t, 0])
+            cy = np.mean(verts[t, 1])
+            eta = evaluate_density(cx, cy, geometry, params)
+            target_spacing = density_to_spacing(eta, params)
+            per_tri_max[i] = _max_area_for_spacing(target_spacing)
+
+            # Check if this triangle actually needs refinement
+            v0, v1, v2 = verts[t[0]], verts[t[1]], verts[t[2]]
+            actual_area = 0.5 * abs(
+                (v1[0] - v0[0]) * (v2[1] - v0[1]) -
+                (v2[0] - v0[0]) * (v1[1] - v0[1])
+            )
+            if actual_area > per_tri_max[i] * 1.3:
+                needs_refinement = True
+
+        if not needs_refinement:
+            break
+
+        # Build refinement input with per-triangle area constraints
+        refine_input = {
+            'vertices': verts,
+            'triangles': tris,
+            'segments': pslg['segments'] if len(pslg['segments']) <= len(verts) else np.empty((0, 2), dtype=np.int32),
+            'triangle_max_area': per_tri_max,
+        }
+        if len(pslg['holes']) > 0:
+            refine_input['holes'] = pslg['holes']
+
+        try:
+            result = tr.triangulate(refine_input, 'rpq25')
+            verts = result.get('vertices', verts)
+            tris = result.get('triangles', tris)
+        except Exception:
+            break
+
+    # Step 6: Filter degenerate / tiny triangles
+    min_area_filter = float(params.get('min_triangle_area', 20.0))
+    if len(tris) > 0:
+        v0 = verts[tris[:, 0]]
+        v1 = verts[tris[:, 1]]
+        v2 = verts[tris[:, 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]
+        tris = tris[areas >= min_area_filter]
 
-    return {'vertices': all_pts, 'triangles': triangles}
+    return {'vertices': verts, 'triangles': tris}
diff --git a/tools/adaptive-isogrid/tests/sandbox2_brain_input.json b/tools/adaptive-isogrid/tests/sandbox2_brain_input.json
new file mode 100644
index 00000000..9566c5e8
--- /dev/null
+++ b/tools/adaptive-isogrid/tests/sandbox2_brain_input.json
@@ -0,0 +1,100 @@
+{
+  "plate_id": "sandbox_2",
+  "units": "mm",
+  "thickness": 12.7,
+  "outer_boundary": [
+    [
+      0.0,
+      0.0
+    ],
+    [
+      7.5,
+      -7.5
+    ],
+    [
+      7.5,
+      -22.6
+    ],
+    [
+      22.5,
+      -22.6
+    ],
+    [
+      22.5,
+      -13.496098
+    ],
+    [
+      74.5,
+      -13.496098
+    ],
+    [
+      74.5,
+      -22.6
+    ],
+    [
+      102.5,
+      -22.6
+    ],
+    [
+      102.5,
+      -7.5
+    ],
+    [
+      117.5,
+      -7.5
+    ],
+    [
+      117.5,
+      -22.6
+    ],
+    [
+      140.748693,
+      -22.6
+    ],
+    [
+      140.748693,
+      124.4
+    ],
+    [
+      117.5,
+      124.4
+    ],
+    [
+      117.5,
+      102.5
+    ],
+    [
+      102.5,
+      102.5
+    ],
+    [
+      102.5,
+      124.4
+    ],
+    [
+      7.5,
+      124.4
+    ],
+    [
+      7.5,
+      102.5
+    ],
+    [
+      0.0,
+      95.0
+    ],
+    [
+      -13.5,
+      95.0
+    ],
+    [
+      -13.5,
+      0.0
+    ],
+    [
+      0.0,
+      0.0
+    ]
+  ],
+  "holes": []
+}
\ No newline at end of file