""" 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})" )