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

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 sample_circle(center, radius, num_points=32):
    """Sample a circle as a polygon with `num_points` vertices."""
    cx, cy = center
    angles = np.linspace(0.0, 2.0 * np.pi, num_points, endpoint=False)
    return [[cx + radius * np.cos(a), cy + radius * np.sin(a)] for a in angles]


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 = list(geometry['outer_boundary'])
    # Strip closing duplicate if last == first
    if len(outer) > 2:
        d = np.linalg.norm(np.array(outer[0]) - np.array(outer[-1]))
        if d < 0.01:
            outer = outer[:-1]
    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 boss keepout reservation
    for hole in geometry['holes']:
        diameter = float(hole.get('diameter', 10.0) or 10.0)
        keepout_dist = d_keep * diameter / 2.0

        if hole.get('is_circular', False) and 'center' in hole:
            # Circular boss reservation around hole:
            # r_boss = r_hole + d_keep * hole_diameter / 2
            hole_radius = diameter / 2.0
            boss_radius = hole_radius + keepout_dist
            keepout_boundary = sample_circle(hole['center'], boss_radius, num_points=12)
        else:
            # Fallback for non-circular holes
            keepout_boundary = offset_polygon(hole['boundary'], keepout_dist, inward=False)

        # Strip closing duplicate if present
        if len(keepout_boundary) > 2:
            d = np.linalg.norm(np.array(keepout_boundary[0]) - np.array(keepout_boundary[-1]))
            if d < 0.01:
                keepout_boundary = keepout_boundary[:-1]
        v_start = len(vertices)
        vertices.extend(keepout_boundary)
        n_h = len(keepout_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 this keepout region 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 filter_small_triangles(result, min_triangle_area):
    """Remove triangles smaller than the manufacturing threshold."""
    triangles = result.get('triangles')
    vertices = result.get('vertices')
    if triangles is None or vertices is None or len(triangles) == 0:
        return result

    areas = compute_triangle_areas(vertices, triangles)
    keep_mask = areas >= float(min_triangle_area)

    result['triangle_areas'] = areas
    result['small_triangle_mask'] = ~keep_mask
    result['triangles'] = triangles[keep_mask]
    return result


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)
    
    # Use s_min as the uniform target spacing for initial/non-adaptive pass.
    # Density-adaptive refinement only kicks in when stress results are
    # available (future iterations). For now, uniform triangles everywhere.
    s_target = params['s_min']
    use_adaptive = params.get('adaptive_density', False)

    # Target equilateral triangle area for the chosen spacing
    target_area = (np.sqrt(3) / 4.0) * s_target**2

    # Triangle options: p=PSLG, q30=min angle 30°, a=area constraint, D=conforming
    result = tr.triangulate(pslg, f'pq30Da{target_area:.1f}')

    if use_adaptive:
        # Iterative density-adaptive refinement (for stress-informed passes)
        for iteration in range(max_refinement_passes):
            verts = result['vertices']
            tris = result['triangles']

            areas = compute_triangle_areas(verts, tris)
            centroids = compute_centroids(verts, tris)

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

            if np.all(areas <= target_areas * 1.2):
                break

            result['triangle_max_area'] = target_areas
            result = tr.triangulate(result, 'rpq30D')

    min_triangle_area = params.get('min_triangle_area', 20.0)
    result = filter_small_triangles(result, min_triangle_area)
    return result