Atomizer/tools/adaptive-isogrid/src/brain/triangulation.py

"""
Constrained Delaunay triangulation with density-adaptive refinement.

Uses Shewchuk's Triangle library to generate adaptive mesh that
respects plate boundary, hole keepouts, and density-driven spacing.
"""

import numpy as np
import triangle as tr
from shapely.geometry import Polygon, LinearRing

from .density_field import evaluate_density, density_to_spacing


def offset_polygon(coords, distance, inward=True):
    """
    Offset a polygon boundary by `distance`.
    inward=True shrinks, inward=False expands.
    
    Uses Shapely buffer (negative = inward for exterior ring).
    """
    poly = Polygon(coords)
    if inward:
        buffered = poly.buffer(-distance)
    else:
        buffered = poly.buffer(distance)
    
    if buffered.is_empty:
        return coords  # can't offset that much, return original
    
    if hasattr(buffered, 'exterior'):
        return list(buffered.exterior.coords)[:-1]  # remove closing duplicate
    return coords


def build_pslg(geometry, params):
    """
    Build Planar Straight Line Graph for Triangle library.
    
    Parameters
    ----------
    geometry : dict
        Plate geometry (outer_boundary, holes).
    params : dict
        Must contain d_keep (hole keepout multiplier).
    
    Returns
    -------
    dict : Triangle-compatible PSLG with vertices, segments, holes.
    """
    d_keep = params['d_keep']
    
    vertices = []
    segments = []
    hole_markers = []
    
    # Outer boundary (no offset needed — frame is handled in profile assembly)
    outer = geometry['outer_boundary']
    v_start = len(vertices)
    vertices.extend(outer)
    n = len(outer)
    for i in range(n):
        segments.append([v_start + i, v_start + (i + 1) % n])
    
    # Each hole with keepout offset
    for hole in geometry['holes']:
        keepout_dist = d_keep * (hole.get('diameter', 10.0) or 10.0) / 2.0
        hole_boundary = offset_polygon(hole['boundary'], keepout_dist, inward=False)
        
        v_start = len(vertices)
        vertices.extend(hole_boundary)
        n_h = len(hole_boundary)
        for i in range(n_h):
            segments.append([v_start + i, v_start + (i + 1) % n_h])
        
        # Marker inside hole tells Triangle to leave it empty
        hole_markers.append(hole['center'])
    
    result = {
        'vertices': np.array(vertices, dtype=np.float64),
        'segments': np.array(segments, dtype=np.int32),
    }
    if hole_markers:
        result['holes'] = np.array(hole_markers, dtype=np.float64)
    
    return result


def compute_triangle_areas(vertices, triangles):
    """Compute area of each triangle."""
    v0 = vertices[triangles[:, 0]]
    v1 = vertices[triangles[:, 1]]
    v2 = vertices[triangles[:, 2]]
    # Cross product / 2
    areas = 0.5 * np.abs(
        (v1[:, 0] - v0[:, 0]) * (v2[:, 1] - v0[:, 1]) -
        (v2[:, 0] - v0[:, 0]) * (v1[:, 1] - v0[:, 1])
    )
    return areas


def compute_centroids(vertices, triangles):
    """Compute centroid of each triangle."""
    v0 = vertices[triangles[:, 0]]
    v1 = vertices[triangles[:, 1]]
    v2 = vertices[triangles[:, 2]]
    return (v0 + v1 + v2) / 3.0


def generate_triangulation(geometry, params, max_refinement_passes=3):
    """
    Generate density-adaptive constrained Delaunay triangulation.
    
    Parameters
    ----------
    geometry : dict
        Plate geometry.
    params : dict
        Full parameter set (density field + spacing + manufacturing).
    max_refinement_passes : int
        Number of iterative refinement passes.
    
    Returns
    -------
    dict : Triangle result with 'vertices' and 'triangles'.
    """
    pslg = build_pslg(geometry, params)
    
    # Initial triangulation with global max area
    s_max = params['s_max']
    global_max_area = (np.sqrt(3) / 4.0) * s_max**2
    
    # Triangle options: p=PSLG, q30=min angle 30°, a=area constraint, D=conforming Delaunay
    result = tr.triangulate(pslg, f'pq30Da{global_max_area:.1f}')
    
    # Iterative refinement based on density field
    for iteration in range(max_refinement_passes):
        verts = result['vertices']
        tris = result['triangles']
        
        areas = compute_triangle_areas(verts, tris)
        centroids = compute_centroids(verts, tris)
        
        # Compute target area for each triangle based on density at centroid
        target_areas = np.array([
            (np.sqrt(3) / 4.0) * density_to_spacing(
                evaluate_density(cx, cy, geometry, params), params
            )**2
            for cx, cy in centroids
        ])
        
        # Check if all triangles satisfy constraints (20% tolerance)
        if np.all(areas <= target_areas * 1.2):
            break
        
        # Set per-triangle max area and refine
        result['triangle_max_area'] = target_areas
        result = tr.triangulate(result, 'rpq30D')
    
    return result