AR-Shipdesign/arshipdesign/geometry/nurbs_surface.py

"""
LoftedSurface — superficie B-spline lofteada a partir de secciones transversales.

El espacio paramétrico es (u, v):
  u ∈ [0, 1]:  dirección longitudinal (0 = AP, 1 = FP)
  v ∈ [0, 1]:  dirección de cuaderna (0 = quilla, 1 = cubierta / línea de agua)

Autor: Álvaro Romero
Sprint 1 — AR-ShipDesign
"""
from __future__ import annotations

from typing import Sequence, Tuple

import numpy as np
from scipy.interpolate import RectBivariateSpline

from .nurbs_curve import BSplineCurve


class LoftedSurface:
    """Superficie lofteada interpolando B-splines entre secciones transversales.

    Parámetros
    ----------
    sections : sequence of (u_pos, points)
        Cada elemento es ``(u_pos, points)`` donde:
        - *u_pos* ∈ [0, 1]: posición longitudinal de la sección.
        - *points*: array (m, 2) con columnas [y_half_breadth, z_height].
          Los puntos deben estar ordenados de quilla (v=0) a cubierta (v=1).
    degree_u : int
        Grado B-spline en dirección longitudinal (por defecto 3).
    degree_v : int
        Grado B-spline en dirección de cuaderna (por defecto 3).
    n_v : int
        Número de puntos uniformes en v al que se reparametrizan todas las
        secciones antes de construir la superficie bivariate.

    Notas
    -----
    Internamente se usa ``scipy.interpolate.RectBivariateSpline`` que exige
    una grilla rectangular (u_i, v_j). Para ello todas las secciones se
    evalúan en el mismo vector v uniforme con ``n_v`` puntos.
    """

    def __init__(
        self,
        sections: Sequence[Tuple[float, np.ndarray]],
        degree_u: int = 3,
        degree_v: int = 3,
        n_v: int = 50,
    ) -> None:
        if len(sections) < degree_u + 1:
            raise ValueError(
                f"Se necesitan al menos {degree_u + 1} secciones para grado_u={degree_u}"
            )

        # Ordenar secciones por posición u
        sections = sorted(sections, key=lambda s: s[0])
        u_positions = np.array([s[0] for s in sections], dtype=float)
        if not np.all(np.diff(u_positions) > 0):
            raise ValueError("Las posiciones u de las secciones deben ser estrictamente crecientes")

        self._u_positions = u_positions
        self._degree_u = degree_u
        self._degree_v = degree_v
        self._n_v = n_v
        self._v_grid = np.linspace(0.0, 1.0, n_v)

        # Construir curvas individuales por sección y remuestrear en v_grid
        self._section_curves: list[BSplineCurve] = []
        y_grid = np.zeros((len(sections), n_v))
        z_grid = np.zeros((len(sections), n_v))

        for i, (_, pts) in enumerate(sections):
            pts = np.asarray(pts, dtype=float)
            if pts.ndim != 2 or pts.shape[1] != 2:
                raise ValueError(
                    f"La sección {i} debe ser array (m, 2) con columnas [y, z]"
                )
            k = min(degree_v, pts.shape[0] - 1)
            curve = BSplineCurve(pts, degree=k)
            self._section_curves.append(curve)
            sampled = curve.evaluate(self._v_grid)   # (n_v, 2)
            y_grid[i, :] = sampled[:, 0]             # half-breadth
            z_grid[i, :] = sampled[:, 1]             # height

        # Superficies bivariadas separadas para y y z
        self._spline_y = RectBivariateSpline(
            u_positions, self._v_grid, y_grid, kx=degree_u, ky=degree_v
        )
        self._spline_z = RectBivariateSpline(
            u_positions, self._v_grid, z_grid, kx=degree_u, ky=degree_v
        )

    # ------------------------------------------------------------------
    # Evaluación
    # ------------------------------------------------------------------

    def evaluate(
        self, u: float | np.ndarray, v: float | np.ndarray
    ) -> np.ndarray:
        """Evalúa la superficie en (u, v).

        Para escalares retorna array (2,) = [y, z].
        Para arrays 1-D de igual longitud retorna (n, 2).
        Para escalares cruzados con arrays retorna la malla completa (nu, nv, 2).
        """
        u = np.asarray(u, dtype=float)
        v = np.asarray(v, dtype=float)
        scalar = u.ndim == 0 and v.ndim == 0

        u = np.atleast_1d(np.clip(u, self._u_positions[0], self._u_positions[-1]))
        v = np.atleast_1d(np.clip(v, 0.0, 1.0))

        y = self._spline_y(u, v)  # (nu, nv)
        z = self._spline_z(u, v)  # (nu, nv)

        if scalar:
            return np.array([y[0, 0], z[0, 0]])

        # Si u y v tienen la misma longitud → puntos pareados
        if u.shape == v.shape:
            pts = np.stack(
                [self._spline_y(u[i:i+1], v[i:i+1])[0, 0]
                 for i in range(len(u))],
            )
            # reconstruir como (n, 2)
            n = len(u)
            result = np.empty((n, 2))
            for i in range(n):
                result[i, 0] = self._spline_y(
                    np.array([u[i]]), np.array([v[i]])
                )[0, 0]
                result[i, 1] = self._spline_z(
                    np.array([u[i]]), np.array([v[i]])
                )[0, 0]
            return result

        # Caso general → malla (nu, nv, 2)
        return np.stack([y, z], axis=-1)

    def section_at(self, u: float) -> BSplineCurve:
        """Sección transversal (y, z) a la posición longitudinal *u* ∈ [0, 1]."""
        u_clip = float(np.clip(u, self._u_positions[0], self._u_positions[-1]))
        pts = np.column_stack([
            self._spline_y(np.array([u_clip]), self._v_grid)[0],
            self._spline_z(np.array([u_clip]), self._v_grid)[0],
        ])
        k = min(self._degree_v, pts.shape[0] - 1)
        return BSplineCurve(pts, degree=k)

    def waterplane(self, v: float = 1.0) -> BSplineCurve:
        """Línea de agua a la posición de cuaderna *v* ∈ [0, 1].

        v=1.0 → línea de agua en la cubierta / calado máximo.
        Retorna curva (u, y) en el plano horizontal.
        """
        v_clip = float(np.clip(v, 0.0, 1.0))
        u_arr = np.linspace(
            self._u_positions[0], self._u_positions[-1], 200
        )
        y_arr = self._spline_y(u_arr, np.array([v_clip]))[:, 0]
        pts = np.column_stack([u_arr, y_arr])
        k = min(self._degree_u, pts.shape[0] - 1)
        return BSplineCurve(pts, degree=k)

    # ------------------------------------------------------------------
    # Propiedades
    # ------------------------------------------------------------------

    @property
    def u_range(self) -> tuple[float, float]:
        return float(self._u_positions[0]), float(self._u_positions[-1])

    @property
    def n_sections(self) -> int:
        return len(self._section_curves)

    def __repr__(self) -> str:
        return (
            f"LoftedSurface(n_sections={self.n_sections}, "
            f"degree_u={self._degree_u}, degree_v={self._degree_v})"
        )