gistfile1.txt

import numpy as np
from shapely.geometry import Polygon, Point
from geometry_utils import is_inside_polygon

def parameterize_boundary(poly, min_points=50):
    """Parameterize a polygon boundary by arc length, ensuring sufficient points."""
    if not poly.is_valid:
        raise ValueError("Invalid polygon provided for parameterization.")
    coords = list(poly.exterior.coords)[:-1]
    if len(coords) < 3:
        raise ValueError(f"Polygon has too few points: {len(coords)}")
    
    # Remove duplicate or near-duplicate points
    cleaned_coords = [coords[0]]
    for i in range(1, len(coords)):
        dist = np.sqrt((coords[i][0] - cleaned_coords[-1][0])**2 + 
                       (coords[i][1] - cleaned_coords[-1][1])**2)
        if dist > 1e-6:
            cleaned_coords.append(coords[i])
    coords = cleaned_coords
    if len(coords) < 3:
        raise ValueError(f"After removing duplicates, too few points remain: {len(coords)}")

    if len(coords) < min_points:
        from scipy.interpolate import splev, splprep
        x_coords, y_coords = np.array(coords)[:, 0], np.array(coords)[:, 1]
        if np.any(np.isnan(x_coords)) or np.any(np.isnan(y_coords)) or \
           np.any(np.isinf(x_coords)) or np.any(np.isinf(y_coords)):
            raise ValueError("Invalid coordinates (NaN or inf) detected in polygon points")
        
        print(f"Before splprep: {len(coords)} points, coords={coords}")
        try:
            tck, u = splprep([x_coords, y_coords], s=0, per=True)
        except Exception as e:
            print(f"splprep failed with error: {e}")
            u = np.linspace(0, 1, len(coords))
            u_fine = np.linspace(0, 1, min_points)
            x_fine = np.interp(u_fine, u, x_coords)
            y_fine = np.interp(u_fine, u, y_coords)
            coords = list(zip(x_fine, y_fine))
        else:
            u_fine = np.linspace(0, 1, min_points)
            x_fine, y_fine = splev(u_fine, tck)
            coords = list(zip(x_fine, y_fine))
    
    # Recheck for duplicates after interpolation
    cleaned_coords = [coords[0]]
    for i in range(1, len(coords)):
        dist = np.sqrt((coords[i][0] - cleaned_coords[-1][0])**2 + 
                       (coords[i][1] - cleaned_coords[-1][1])**2)
        if dist > 1e-6:
            cleaned_coords.append(coords[i])
    coords = cleaned_coords
    if len(coords) < 3:
        raise ValueError(f"After interpolation and cleaning, too few points remain: {len(coords)}")

    lengths = [np.sqrt((coords[i][0] - coords[i-1][0])**2 + 
                       (coords[i][1] - coords[i-1][1])**2) for i in range(len(coords))]
    total_length = sum(lengths)
    if total_length < 1e-6:
        raise ValueError("Polygon has near-zero perimeter.")
    eta = np.cumsum([0] + lengths) / total_length
    print(f"Parameterized boundary with {len(coords)} points")
    return coords, eta

def generate_harmonic_grid(hull, core, desired_cell_area=1.0):
    """Generate a grid using transfinite interpolation for an elliptical core."""
    # Validate and compute floor area
    if not hull.contains(core):
        raise ValueError("Core must be fully contained within hull.")
    floor = hull.difference(core.buffer(0.05))
    if not floor.is_valid:
        raise ValueError("Invalid floor geometry (check hull/core overlap or geometry).")
    area = floor.area
    print(f"Floor area: {area:.2f} square meters")

    # Parameterize boundaries
    core_coords, core_eta = parameterize_boundary(core, min_points=87)
    N = len(core_coords)
    hull_coords, hull_eta = parameterize_boundary(hull, min_points=87)
    print(f"Hull points: {len(hull_coords)}, Core points: {len(core_coords)}")

    # Determine grid resolution with adjusted scaling
    M = max(int(np.round(area / (N * 0.5))), 10)
    if M < 4 or N < 4:
        raise ValueError(f"Floor area too small for grid (M={M}, N={N}).")
    print(f"Grid resolution: {M}x{N}")

    # Computational grid
    xi = np.linspace(0, 1, M)
    eta = np.linspace(0, 1, N)

    # Initialize physical coordinates using transfinite interpolation
    x = np.zeros((M, N))
    y = np.zeros((M, N))
    hull_indices = []
    core_indices = []
    for j in range(N):
        t = eta[j]
        idx_hull = min(int(t * len(hull_coords)), len(hull_coords) - 1)
        idx_core = min(int(t * len(core_coords)), len(core_coords) - 1)
        hull_indices.append(idx_hull)
        core_indices.append(idx_core)
        hull_point = hull_coords[idx_hull]
        core_point = core_coords[idx_core]
        for i in range(M):
            if i == 0:
                x[i, j] = core_point[0]
                y[i, j] = core_point[1]
            else:
                x[i, j] = (1 - xi[i]) * core_point[0] + xi[i] * hull_point[0]
                y[i, j] = (1 - xi[i]) * core_point[1] + xi[i] * hull_point[1]

    # Verify boundary alignment
    for j in range(N):
        idx_core = core_indices[j]
        idx_hull = hull_indices[j]
        if not (abs(x[0, j] - core_coords[idx_core][0]) < 1e-6 and abs(y[0, j] - core_coords[idx_core][1]) < 1e-6):
            print(f"Warning: Core boundary point at j={j} misaligned: ({x[0, j]}, {y[0, j]}) vs ({core_coords[idx_core][0]}, {core_coords[idx_core][1]})")
        if not (abs(x[M-1, j] - hull_coords[idx_hull][0]) < 1e-6 and abs(y[M-1, j] - hull_coords[idx_hull][1]) < 1e-6):
            print(f"Warning: Hull boundary point at j={j} misaligned: ({x[M-1, j]}, {y[M-1, j]}) vs ({hull_coords[idx_hull][0]}, {hull_coords[idx_hull][1]})")

    # Generate grid cells, skipping i=0 to avoid core intrusion
    cells = []
    for i in range(1, M-1):  # Start from i=1 to skip the core boundary
        for j in range(N):
            j_next = (j + 1) % N
            cell_points = [
                (x[i, j], y[i, j]),
                (x[i+1, j], y[i+1, j]),
                (x[i+1, j_next], y[i+1, j_next]),
                (x[i, j_next], y[i, j_next])
            ]
            cell_poly = Polygon(cell_points)
            if not cell_poly.is_valid:
                print(f"Warning: Invalid cell at (i={i}, j={j})")
                continue
            # Check if the cell intersects the core
            intersects_core = core.intersects(cell_poly)
            if intersects_core:
                print(f"Debug: Cell at (i={i}, j={j}) intrudes into core. Cell points: {cell_points}")
                cells.append({'points': cell_points, 'status': 'void'})
                continue
            # Check if cell is in the floor area
            centroid = cell_poly.centroid
            in_floor = floor.contains(centroid)
            if in_floor:
                cells.append({'points': cell_points, 'status': 'usable'})
            else:
                cells.append({'points': cell_points, 'status': 'void'})
    return cells, x, y, M, N