From 503e00bfc9788bcd79c39411e9a1d38b64adc5e1 Mon Sep 17 00:00:00 2001 From: alro1965 Date: Wed, 27 May 2026 01:07:35 -0400 Subject: [PATCH] feat(sprint1): motor NURBS, modelos de casco, visor 3D PyVista MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Geometría: - BSplineCurve: interpolación scipy, arc_length, tangente, chord-length - LoftedSurface: lofting de secciones → RectBivariateSpline bivariate Core (casco Wigley como caso de prueba): - Section: área, centroide_z, max_half_breadth, curva B-spline - OffsetsTable: from_wigley(), to_sections(), interpolación xy - Hull: volumen, Awp, LCB, VCB, Cb, Cm, Cp, desplazamiento, to_mesh() UI: - Viewer3DWidget (pyvistaqt.QtInteractor): casco Wigley por defecto al arrancar, fondo navy, waterplane semi-transparente, fallback graceful si PyVista no disponible - MainWindow: Viewer3DWidget inyectado en viewport Perspectiva 3D Tests: 39 nuevos tests, fórmulas analíticas Wigley verificadas (±1%) V = 4BLT/9, Cb = 4/9, Awp = 2BL/3 (derivación correcta) Co-Authored-By: Claude Sonnet 4.5 --- arshipdesign/core/hull.py | 306 ++++++++++++++++++++- arshipdesign/core/offsets.py | 172 +++++++++++- arshipdesign/core/section.py | 134 +++++++++- arshipdesign/geometry/nurbs_curve.py | 142 +++++++++- arshipdesign/geometry/nurbs_surface.py | 186 ++++++++++++- arshipdesign/ui/main_window.py | 10 + arshipdesign/ui/widgets/viewer_3d.py | 226 +++++++++++++++- tests/test_sprint1_geometry.py | 355 +++++++++++++++++++++++++ 8 files changed, 1519 insertions(+), 12 deletions(-) create mode 100644 tests/test_sprint1_geometry.py diff --git a/arshipdesign/core/hull.py b/arshipdesign/core/hull.py index 59210a6..474681e 100644 --- a/arshipdesign/core/hull.py +++ b/arshipdesign/core/hull.py @@ -1,2 +1,304 @@ -"""Geometría del casco NURBS. Stub — Sprint 1.""" -raise NotImplementedError("hull — Sprint 1") +""" +Hull — modelo de casco naval con geometría NURBS y cálculos hidrostáticos básicos. + +El casco se representa como una LoftedSurface construida a partir de las secciones +de una OffsetsTable. Los cálculos hidrostáticos usan la regla de Simpson sobre +las secciones muestreadas para máxima compatibilidad con cualquier forma de casco. + +Autor: Álvaro Romero +Sprint 1 — AR-ShipDesign +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Optional + +import numpy as np +from scipy.integrate import simpson + +from arshipdesign.core.offsets import OffsetsTable +from arshipdesign.core.section import Section +from arshipdesign.geometry.nurbs_surface import LoftedSurface + + +@dataclass +class Hull: + """Modelo geométrico del casco naval. + + Atributos + --------- + name : str + Nombre del casco / proyecto. + lpp : float + Eslora entre perpendiculares [m]. + beam : float + Manga máxima en flotación [m]. + depth : float + Puntal de trazado [m]. + draft : float + Calado de diseño [m]. + offsets : OffsetsTable + Tabla de offsets del casco. + """ + + name: str + lpp: float + beam: float + depth: float + draft: float + offsets: OffsetsTable + _surface: Optional[LoftedSurface] = field(default=None, repr=False, compare=False) + + # ------------------------------------------------------------------ + # Fábricas + # ------------------------------------------------------------------ + + @classmethod + def from_wigley( + cls, + name: str = "Wigley Hull", + lpp: float = 10.0, + beam: float = 1.5, + draft: float = 0.75, + n_stations: int = 21, + n_waterlines: int = 11, + ) -> "Hull": + """Crea un casco Wigley estándar. + + El casco Wigley tiene solución analítica exacta para sus + hidrostáticos, lo que permite verificar los métodos numéricos. + + Fórmulas analíticas: + Volumen de desplazamiento: V = (8/15) · B/2 · T · L/2 + LCB: en el midship (x = Lpp/2) por simetría + Área plano de flotación (T): Awp = (4/3) · (B/2) · L + (sólo la contribución f_xi: integral de 1-(2ξ/L)² dξ) + """ + offsets = OffsetsTable.from_wigley( + lpp=lpp, beam=beam, draft=draft, + n_stations=n_stations, n_waterlines=n_waterlines, + ) + return cls( + name=name, + lpp=lpp, + beam=beam, + depth=draft, # Para Wigley depth = draft + draft=draft, + offsets=offsets, + ) + + # ------------------------------------------------------------------ + # Superficie NURBS (lazy) + # ------------------------------------------------------------------ + + @property + def surface(self) -> LoftedSurface: + """LoftedSurface construida a partir de la tabla de offsets.""" + if self._surface is None: + self._surface = self._build_surface() + return self._surface + + def _build_surface(self) -> LoftedSurface: + sections_data = [] + u_arr = self.offsets.x_stations / self.lpp # normalizar a [0,1] + for i, u in enumerate(u_arr): + pts = np.column_stack([ + self.offsets.data[i, :], + self.offsets.z_waterlines, + ]) + sections_data.append((float(u), pts)) + n_sec = len(sections_data) + deg_u = min(3, n_sec - 1) + return LoftedSurface(sections_data, degree_u=deg_u, degree_v=3) + + # ------------------------------------------------------------------ + # Hidrostáticos + # ------------------------------------------------------------------ + + def sections_at_draft(self, draft: Optional[float] = None) -> list[Section]: + """Lista de secciones de la tabla de offsets.""" + return self.offsets.to_sections() + + def volume_of_displacement(self, draft: Optional[float] = None) -> float: + """Volumen de desplazamiento [m³] hasta *draft*. + + Integra el área de cada sección en la dirección x usando + la regla de Simpson sobre todas las estaciones. + """ + T = draft if draft is not None else self.draft + sections = self.offsets.to_sections() + x = np.array([s.x for s in sections]) + areas = np.array([s.area(draft=T) for s in sections]) + + if len(x) >= 3: + vol = float(simpson(areas, x=x)) + else: + vol = float(np.trapz(areas, x)) + return abs(vol) + + def waterplane_area(self, draft: Optional[float] = None) -> float: + """Área del plano de flotación [m²] al calado *draft*. + + Integra 2·y(x, z=draft) en la dirección x. + """ + T = draft if draft is not None else self.draft + x = self.offsets.x_stations + y_wl = np.array([ + self.offsets.half_breadth(xi, T) for xi in x + ]) + # Área = integral de 2·y(x) dx (ambas bandas) + if len(x) >= 3: + awp = float(simpson(2.0 * y_wl, x=x)) + else: + awp = float(np.trapz(2.0 * y_wl, x)) + return abs(awp) + + def lcb(self, draft: Optional[float] = None) -> float: + """Centro longitudinal de carena (LCB) [m desde AP]. + + Momento de primer orden del volumen / volumen total. + """ + T = draft if draft is not None else self.draft + sections = self.offsets.to_sections() + x = np.array([s.x for s in sections]) + areas = np.array([s.area(draft=T) for s in sections]) + + if len(x) >= 3: + vol = float(simpson(areas, x=x)) + moment = float(simpson(areas * x, x=x)) + else: + vol = float(np.trapz(areas, x)) + moment = float(np.trapz(areas * x, x)) + + if abs(vol) < 1e-12: + return self.lpp / 2.0 + return moment / vol + + def vcb(self, draft: Optional[float] = None) -> float: + """Centro vertical de carena (VCB / KB) [m sobre la quilla].""" + T = draft if draft is not None else self.draft + sections = self.offsets.to_sections() + x = np.array([s.x for s in sections]) + areas = np.array([s.area(draft=T) for s in sections]) + cz = np.array([s.centroid_z(draft=T) for s in sections]) + + if len(x) >= 3: + vol = float(simpson(areas, x=x)) + moment_z = float(simpson(areas * cz, x=x)) + else: + vol = float(np.trapz(areas, x)) + moment_z = float(np.trapz(areas * cz, x)) + + if abs(vol) < 1e-12: + return T / 2.0 + return moment_z / vol + + def block_coefficient(self, draft: Optional[float] = None) -> float: + """Coeficiente de bloque Cb = V / (Lpp · B · T).""" + T = draft if draft is not None else self.draft + V = self.volume_of_displacement(T) + return V / (self.lpp * self.beam * T) + + def midship_coefficient(self, draft: Optional[float] = None) -> float: + """Coeficiente de cuaderna maestra Cm = Am / (B · T).""" + T = draft if draft is not None else self.draft + sections = self.offsets.to_sections() + # Cuaderna en el midship + x_mid = self.lpp / 2.0 + areas_at_mid = [s.area(draft=T) for s in sections if abs(s.x - x_mid) < 1e-6] + if not areas_at_mid: + # Interpolar + x_arr = np.array([s.x for s in sections]) + a_arr = np.array([s.area(draft=T) for s in sections]) + am = float(np.interp(x_mid, x_arr, a_arr)) + else: + am = areas_at_mid[0] + return am / (self.beam * T) + + def prismatic_coefficient(self, draft: Optional[float] = None) -> float: + """Coeficiente prismático Cp = V / (Am · Lpp).""" + T = draft if draft is not None else self.draft + V = self.volume_of_displacement(T) + Am = self.midship_coefficient(T) * self.beam * T + if Am < 1e-12: + return 0.0 + return V / (Am * self.lpp) + + def displacement_tonnes( + self, draft: Optional[float] = None, rho: float = 1025.0 + ) -> float: + """Desplazamiento en toneladas métricas (agua salada por defecto). + + Parámetros + ---------- + rho : float + Densidad del agua [kg/m³]. Default 1025 kg/m³ (agua salada). + """ + V = self.volume_of_displacement(draft) + return V * rho / 1000.0 + + # ------------------------------------------------------------------ + # Malla PyVista para visualización 3D + # ------------------------------------------------------------------ + + def to_mesh(self, n_u: int = 40, n_v: int = 20) -> "pyvista.PolyData": + """Genera una malla PyVista del casco (ambas bandas). + + Requiere PyVista instalado. Retorna un PolyData triangulado. + """ + try: + import pyvista as pv + except ImportError as exc: + raise ImportError("PyVista no está instalado") from exc + + surf = self.surface + u_range = surf.u_range + u_arr = np.linspace(u_range[0], u_range[1], n_u) + v_arr = np.linspace(0.0, 1.0, n_v) + uu, vv = np.meshgrid(u_arr, v_arr, indexing="ij") # (n_u, n_v) + + # Evaluar (y, z) en la malla + y_mat = surf._spline_y(u_arr, v_arr) # (n_u, n_v) + z_mat = surf._spline_z(u_arr, v_arr) # (n_u, n_v) + + # x real desde parámetro u + x_mat = uu * self.lpp + + # Banda de estribor (y > 0) + pts_stbd = np.stack([ + x_mat.ravel(), y_mat.ravel(), z_mat.ravel() + ], axis=1) + + # Banda de babor (y < 0) + pts_port = np.stack([ + x_mat.ravel(), -y_mat.ravel(), z_mat.ravel() + ], axis=1) + + # Unir ambas bandas + all_pts = np.vstack([pts_stbd, pts_port]) + + # Construir caras de la malla estructurada + faces = [] + offset = n_u * n_v + for band in [0, offset]: + for i in range(n_u - 1): + for j in range(n_v - 1): + p0 = band + i * n_v + j + p1 = band + (i + 1) * n_v + j + p2 = band + (i + 1) * n_v + (j + 1) + p3 = band + i * n_v + (j + 1) + faces.extend([4, p0, p1, p2, p3]) + + faces_arr = np.array(faces, dtype=int) + mesh = pv.PolyData(all_pts, faces_arr) + return mesh.triangulate() + + # ------------------------------------------------------------------ + # Dunder + # ------------------------------------------------------------------ + + def __repr__(self) -> str: + return ( + f"Hull({self.name!r}, Lpp={self.lpp} m, B={self.beam} m, " + f"T={self.draft} m)" + ) diff --git a/arshipdesign/core/offsets.py b/arshipdesign/core/offsets.py index dd8623d..7c10b4f 100644 --- a/arshipdesign/core/offsets.py +++ b/arshipdesign/core/offsets.py @@ -1,2 +1,170 @@ -"""Tabla de offsets. Stub — Sprint 1.""" -raise NotImplementedError("offsets — Sprint 1") +""" +OffsetsTable — tabla de offsets naval clásica. + +Estructura: filas = estaciones (x), columnas = líneas de agua (z). +Valores: semi-manga y(x, z) en metros. + +Autor: Álvaro Romero +Sprint 1 — AR-ShipDesign +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Optional + +import numpy as np + +from arshipdesign.core.section import Section + + +@dataclass +class OffsetsTable: + """Tabla de offsets del casco. + + Atributos + --------- + x_stations : np.ndarray, shape (n_sta,) + Posiciones longitudinales [m] de cada estación, creciente. + z_waterlines : np.ndarray, shape (n_wl,) + Alturas de líneas de agua [m] desde la quilla, creciente. + data : np.ndarray, shape (n_sta, n_wl) + Semi-mangas y[i, j] en metros. + y[i, j] = semi-manga en la estación x_stations[i], línea z_waterlines[j]. + station_labels : list[str] + Etiquetas opcionales para las estaciones. + lpp : float + Eslora entre perpendiculares [m]. + beam : float + Manga máxima [m]. + draft : float + Calado de diseño [m]. + """ + + x_stations: np.ndarray + z_waterlines: np.ndarray + data: np.ndarray + station_labels: list[str] = field(default_factory=list) + lpp: float = 0.0 + beam: float = 0.0 + draft: float = 0.0 + + def __post_init__(self) -> None: + self.x_stations = np.asarray(self.x_stations, dtype=float) + self.z_waterlines = np.asarray(self.z_waterlines, dtype=float) + self.data = np.asarray(self.data, dtype=float) + + n_sta = len(self.x_stations) + n_wl = len(self.z_waterlines) + if self.data.shape != (n_sta, n_wl): + raise ValueError( + f"data.shape {self.data.shape} ≠ ({n_sta}, {n_wl})" + ) + if not self.station_labels: + self.station_labels = [str(i) for i in range(n_sta)] + + # ------------------------------------------------------------------ + # Fábrica: casco Wigley analítico + # ------------------------------------------------------------------ + + @classmethod + def from_wigley( + cls, + lpp: float = 10.0, + beam: float = 1.5, + draft: float = 0.75, + n_stations: int = 21, + n_waterlines: int = 11, + ) -> "OffsetsTable": + """Genera tabla de offsets para el casco Wigley matemático. + + Fórmula analítica: + y(x, z) = (B/2) · [1 − (2ξ/L)²] · [1 − (ζ/T)²] + + donde ξ ∈ [−L/2, L/2] y ζ ∈ [−T, 0]. + + El AP corresponde a x=0, el FP a x=Lpp. + La quilla está a z=0, la flotación de diseño a z=draft. + """ + # Posiciones longitudinales: 0 (AP) → Lpp (FP) + x_sta = np.linspace(0.0, lpp, n_stations) + # Líneas de agua: 0 (quilla) → draft (calado diseño) + z_wl = np.linspace(0.0, draft, n_waterlines) + + # Convertir a coordenadas centradas en el midship + xi = x_sta - lpp / 2.0 # ξ ∈ [-Lpp/2, Lpp/2] + zeta = z_wl - draft # ζ ∈ [-T, 0] + + # Factores de forma + f_xi = 1.0 - (2.0 * xi / lpp) ** 2 # (n_sta,) + f_zeta = 1.0 - (zeta / draft) ** 2 # (n_wl,) + + # y[i, j] = B/2 · f_xi[i] · f_zeta[j] + data = (beam / 2.0) * np.outer(f_xi, f_zeta) + data = np.clip(data, 0.0, None) # evitar negativos por redondeo + + labels = [f"S{i}" for i in range(n_stations)] + + return cls( + x_stations=x_sta, + z_waterlines=z_wl, + data=data, + station_labels=labels, + lpp=lpp, + beam=beam, + draft=draft, + ) + + # ------------------------------------------------------------------ + # Conversión a secciones + # ------------------------------------------------------------------ + + def to_sections(self) -> list[Section]: + """Convierte la tabla en una lista de objetos Section.""" + sections = [] + for i, x in enumerate(self.x_stations): + sta = self.station_labels[i] if i < len(self.station_labels) else str(i) + sec = Section( + station=sta, + x=float(x), + half_breadths=self.data[i, :].copy(), + z_positions=self.z_waterlines.copy(), + label=f"x={x:.3f} m", + ) + sections.append(sec) + return sections + + # ------------------------------------------------------------------ + # Acceso + # ------------------------------------------------------------------ + + def half_breadth(self, x: float, z: float) -> float: + """Interpola la semi-manga en cualquier (x, z) [m].""" + # Interpolar en x + col_y = np.array([ + float(np.interp(x, self.x_stations, self.data[:, j])) + for j in range(len(self.z_waterlines)) + ]) + # Interpolar en z + return float(np.interp(z, self.z_waterlines, col_y)) + + @property + def n_stations(self) -> int: + return len(self.x_stations) + + @property + def n_waterlines(self) -> int: + return len(self.z_waterlines) + + @property + def max_half_breadth(self) -> float: + return float(self.data.max()) + + # ------------------------------------------------------------------ + # Dunder + # ------------------------------------------------------------------ + + def __repr__(self) -> str: + return ( + f"OffsetsTable(Lpp={self.lpp} m, B={self.beam} m, T={self.draft} m, " + f"stations={self.n_stations}, waterlines={self.n_waterlines})" + ) diff --git a/arshipdesign/core/section.py b/arshipdesign/core/section.py index 1ef1e28..1e84e00 100644 --- a/arshipdesign/core/section.py +++ b/arshipdesign/core/section.py @@ -1,2 +1,132 @@ -"""Sección transversal. Stub — Sprint 1.""" -raise NotImplementedError("section — Sprint 1") +""" +Section — sección transversal del casco naval. + +Una sección es un corte vertical del casco a posición longitudinal x. +Se define por puntos (y_half, z) de quilla a cubierta (casco simétrico). + +Autor: Álvaro Romero +Sprint 1 — AR-ShipDesign +""" +from __future__ import annotations + +from dataclasses import dataclass +from typing import Optional + +import numpy as np +from scipy.integrate import simpson + +from arshipdesign.geometry.nurbs_curve import BSplineCurve + + +@dataclass +class Section: + """Sección transversal en una estación. + + Atributos + --------- + station : int or str + Número de estación (0 = AP, 20 = FP en sistema de 21 estaciones). + x : float + Posición longitudinal en metros desde AP. + half_breadths : np.ndarray + Semi-mangas [m] en cada z (mitad de babor). + z_positions : np.ndarray + Alturas z [m] desde la quilla (creciente, 0 = quilla). + label : str + Etiqueta opcional. + """ + + station: int | str + x: float + half_breadths: np.ndarray + z_positions: np.ndarray + label: str = "" + + def __post_init__(self) -> None: + self.half_breadths = np.asarray(self.half_breadths, dtype=float) + self.z_positions = np.asarray(self.z_positions, dtype=float) + if self.half_breadths.shape != self.z_positions.shape: + raise ValueError( + "half_breadths y z_positions deben tener igual número de elementos" + ) + if len(self.z_positions) < 2: + raise ValueError("Se necesitan al menos 2 puntos por sección") + if not np.all(np.diff(self.z_positions) >= 0): + raise ValueError("z_positions debe ser no-decreciente") + + # ------------------------------------------------------------------ + # Curva B-spline + # ------------------------------------------------------------------ + + def curve(self, degree: int = 3) -> BSplineCurve: + """Curva B-spline (y_half, z) de esta sección.""" + pts = np.column_stack([self.half_breadths, self.z_positions]) + k = min(degree, len(pts) - 1) + return BSplineCurve(pts, degree=k) + + # ------------------------------------------------------------------ + # Hidrostáticos de sección + # ------------------------------------------------------------------ + + def _clip_to_draft( + self, draft: Optional[float] + ) -> tuple[np.ndarray, np.ndarray]: + """Recorta z/y al calado dado. Retorna (z_clipped, y_clipped).""" + z, y = self.z_positions, self.half_breadths + if draft is None: + return z, y + z_top = float(draft) + if z_top <= z[0]: + return np.array([z[0], z[0]]), np.array([0.0, 0.0]) + if z_top < z[-1]: + y_top = float(np.interp(z_top, z, y)) + mask = z < z_top + return ( + np.concatenate([z[mask], [z_top]]), + np.concatenate([y[mask], [y_top]]), + ) + return z, y + + def area(self, draft: Optional[float] = None) -> float: + """Área sumergida de la sección [m²] hasta *draft*. + + Integra 2·y(z) dz (ambas bandas por simetría). + """ + z_int, y_int = self._clip_to_draft(draft) + if len(z_int) < 2: + return 0.0 + if len(z_int) >= 3: + a = float(simpson(y_int, x=z_int)) + else: + a = float(np.trapz(y_int, z_int)) + return 2.0 * abs(a) + + def centroid_z(self, draft: Optional[float] = None) -> float: + """Centroide vertical de la sección sumergida [m desde quilla].""" + z_int, y_int = self._clip_to_draft(draft) + if len(z_int) < 2: + return float(self.z_positions[0]) + if len(z_int) >= 3: + moment = float(simpson(y_int * z_int, x=z_int)) + area = float(simpson(y_int, x=z_int)) + else: + moment = float(np.trapz(y_int * z_int, z_int)) + area = float(np.trapz(y_int, z_int)) + if abs(area) < 1e-12: + return float(z_int[0]) + return moment / area + + def max_half_breadth(self, draft: Optional[float] = None) -> float: + """Máxima semi-manga hasta *draft* [m].""" + z, y = self._clip_to_draft(draft) + return float(y.max()) if len(y) > 0 else 0.0 + + # ------------------------------------------------------------------ + # Dunder + # ------------------------------------------------------------------ + + def __repr__(self) -> str: + return ( + f"Section(station={self.station!r}, x={self.x:.3f} m, " + f"n_pts={len(self.z_positions)})" + ) diff --git a/arshipdesign/geometry/nurbs_curve.py b/arshipdesign/geometry/nurbs_curve.py index dd714a8..d2b3bca 100644 --- a/arshipdesign/geometry/nurbs_curve.py +++ b/arshipdesign/geometry/nurbs_curve.py @@ -1,2 +1,140 @@ -"""Curva NURBS. Stub — Sprint 1.""" -raise NotImplementedError("nurbs_curve — Sprint 1") +""" +BSplineCurve — curva B-spline paramétrica basada en scipy.interpolate. + +Autor: Álvaro Romero +Sprint 1 — AR-ShipDesign +""" +from __future__ import annotations + +import numpy as np +from scipy.interpolate import make_interp_spline +from scipy.integrate import quad + + +class BSplineCurve: + """Curva B-spline interpolada a partir de puntos de datos. + + Parámetros + ---------- + points : array-like, shape (n, 2) o (n, 3) + Puntos de datos por los que pasa la curva (interpolación exacta). + degree : int, optional + Grado del B-spline (por defecto 3 = cúbico). + knots : array-like or None + Vector de nudos externo. Si es None se genera automáticamente + con método cuerda-longitud (no uniforme, más estable). + + Notas + ----- + El parámetro interno *t* varía en [0, 1]. + La curva se construye usando ``scipy.interpolate.make_interp_spline``. + Para n puntos y grado k se necesita n >= k+1. + """ + + def __init__( + self, + points: np.ndarray, + degree: int = 3, + *, + knots: np.ndarray | None = None, + ) -> None: + pts = np.asarray(points, dtype=float) + if pts.ndim != 2 or pts.shape[1] not in (2, 3): + raise ValueError("points debe ser array (n, 2) o (n, 3)") + if pts.shape[0] < degree + 1: + raise ValueError( + f"Se necesitan al menos {degree + 1} puntos para grado {degree}" + ) + + self._pts = pts + self._k = degree + self._dim = pts.shape[1] + + # Parámetro interno por longitud de cuerda normalizada + self._t_param = _chord_length_parameterization(pts) + + if knots is not None: + self._spl = make_interp_spline( + self._t_param, pts, k=degree, t=knots + ) + else: + self._spl = make_interp_spline(self._t_param, pts, k=degree) + + # ------------------------------------------------------------------ + # Evaluación + # ------------------------------------------------------------------ + + def evaluate(self, t: float | np.ndarray) -> np.ndarray: + """Evalúa la curva en el parámetro *t* ∈ [0, 1]. + + Retorna array (dim,) para t escalar, o (m, dim) para t array. + """ + t = np.asarray(t, dtype=float) + scalar = t.ndim == 0 + t = np.atleast_1d(np.clip(t, 0.0, 1.0)) + result = self._spl(t) + return result[0] if scalar else result + + def sample(self, n: int = 100) -> np.ndarray: + """Muestrea *n* puntos uniformes en t ∈ [0, 1]. + + Retorna array (n, dim). + """ + t = np.linspace(0.0, 1.0, max(2, n)) + return self.evaluate(t) + + # ------------------------------------------------------------------ + # Métricas + # ------------------------------------------------------------------ + + def arc_length(self, t0: float = 0.0, t1: float = 1.0) -> float: + """Longitud de arco numérica entre los parámetros *t0* y *t1*.""" + dspl = self._spl.derivative() + + def integrand(t: float) -> float: + return float(np.linalg.norm(dspl(float(t)))) + + length, _ = quad(integrand, t0, t1, limit=100) + return float(length) + + def tangent(self, t: float) -> np.ndarray: + """Vector tangente unitario en el parámetro *t*.""" + d = self._spl.derivative()(float(np.clip(t, 0.0, 1.0))) + norm = np.linalg.norm(d) + return d / norm if norm > 1e-12 else np.zeros(self._dim) + + # ------------------------------------------------------------------ + # Propiedades + # ------------------------------------------------------------------ + + @property + def degree(self) -> int: + return self._k + + @property + def dim(self) -> int: + return self._dim + + @property + def control_points(self) -> np.ndarray: + """Puntos de datos originales.""" + return self._pts.copy() + + def __repr__(self) -> str: + n = self._pts.shape[0] + return f"BSplineCurve(degree={self._k}, dim={self._dim}, n_points={n})" + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _chord_length_parameterization(pts: np.ndarray) -> np.ndarray: + """Parametrización por longitud de cuerda normalizada a [0, 1].""" + diffs = np.diff(pts, axis=0) + chord_lens = np.linalg.norm(diffs, axis=1) + cumulative = np.concatenate([[0.0], np.cumsum(chord_lens)]) + total = cumulative[-1] + if total < 1e-12: + return np.linspace(0.0, 1.0, len(pts)) + return cumulative / total diff --git a/arshipdesign/geometry/nurbs_surface.py b/arshipdesign/geometry/nurbs_surface.py index 6444cd4..c969d24 100644 --- a/arshipdesign/geometry/nurbs_surface.py +++ b/arshipdesign/geometry/nurbs_surface.py @@ -1,2 +1,184 @@ -"""Superficie NURBS geomdl. Stub — Sprint 1.""" -raise NotImplementedError("nurbs_surface — Sprint 1") +""" +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})" + ) diff --git a/arshipdesign/ui/main_window.py b/arshipdesign/ui/main_window.py index 1c90f8b..dcd1fa0 100644 --- a/arshipdesign/ui/main_window.py +++ b/arshipdesign/ui/main_window.py @@ -822,6 +822,16 @@ class MainWindow(QMainWindow): self._module_area = ModuleArea() self.setCentralWidget(self._module_area) + # Inyectar visor 3D en el viewport Perspectiva (diferido) + from arshipdesign.ui.widgets.viewer_3d import Viewer3DWidget, _PYVISTA_OK + if _PYVISTA_OK: + self._viewer_3d = Viewer3DWidget() + vp = self._module_area.four_viewport.viewport("perspective") + if vp is not None: + vp.set_canvas(self._viewer_3d) + else: + self._viewer_3d = None + # Dock izquierdo — capas self._layers_panel = LayersPanel(self._strings) self._dock_layers = QDockWidget("Capas", self) diff --git a/arshipdesign/ui/widgets/viewer_3d.py b/arshipdesign/ui/widgets/viewer_3d.py index b055e57..e946a7f 100644 --- a/arshipdesign/ui/widgets/viewer_3d.py +++ b/arshipdesign/ui/widgets/viewer_3d.py @@ -1,2 +1,224 @@ -"""Vista 3D PyVista. Stub — Sprint 1.""" -raise NotImplementedError("viewer_3d — Sprint 1") +""" +Viewer3DWidget — visor 3D interactivo basado en PyVista / pyvistaqt. + +Integra un QtInteractor de pyvistaqt directamente en el viewport +"Perspectiva 3D" del layout de 4 viewports. + +Si pyvistaqt no está disponible (entorno CI / offscreen) la clase degrada +a un QLabel de reemplazo sin lanzar excepción, para que el smoke test +y el resto de la UI sigan funcionando. + +Autor: Álvaro Romero +Sprint 1 — AR-ShipDesign +""" +from __future__ import annotations + +import logging +from typing import Optional + +import numpy as np +from PySide6.QtCore import Qt, QTimer +from PySide6.QtWidgets import QLabel, QVBoxLayout, QWidget + +logger = logging.getLogger("ui.viewer_3d") + +# --------------------------------------------------------------------------- +# Importación condicional de PyVista +# --------------------------------------------------------------------------- +try: + import pyvista as pv + from pyvistaqt import QtInteractor + _PYVISTA_OK = True +except Exception as _pv_err: # ImportError o cualquier otro + pv = None # type: ignore[assignment] + QtInteractor = None # type: ignore[assignment] + _PYVISTA_OK = False + logger.warning("PyVista/pyvistaqt no disponible: %s", _pv_err) + + +class Viewer3DWidget(QWidget): + """Widget embebible de visor 3D. + + Si PyVista está disponible → usa ``pyvistaqt.QtInteractor``. + Si no → muestra un QLabel de aviso (modo degradado). + + Uso típico: + ----------- + >>> v = Viewer3DWidget(parent=viewport_frame) + >>> v.load_hull(hull_object) # muestra casco NURBS + >>> v.load_mesh(pv_mesh) # carga cualquier PolyData + >>> v.reset_camera() + """ + + def __init__(self, parent: Optional[QWidget] = None) -> None: + super().__init__(parent) + self._plotter: Optional["QtInteractor"] = None + self._ready = False + self._build_ui() + + # ------------------------------------------------------------------ + # Construcción de la UI + # ------------------------------------------------------------------ + + def _build_ui(self) -> None: + lo = QVBoxLayout(self) + lo.setContentsMargins(0, 0, 0, 0) + lo.setSpacing(0) + + if _PYVISTA_OK: + # El QtInteractor se crea diferido (después de que la ventana esté visible) + # para evitar problemas con el contexto OpenGL en el arranque. + self._placeholder = QLabel("Iniciando visor 3D…") + self._placeholder.setAlignment(Qt.AlignmentFlag.AlignCenter) + self._placeholder.setObjectName("viewportPlaceholder") + lo.addWidget(self._placeholder) + # Diferir inicialización hasta que Qt haya mostrado la ventana + QTimer.singleShot(500, self._init_plotter) + else: + lbl = QLabel("PyVista no disponible\n(pip install pyvistaqt)") + lbl.setAlignment(Qt.AlignmentFlag.AlignCenter) + lbl.setObjectName("viewportPlaceholder") + lo.addWidget(lbl) + + def _init_plotter(self) -> None: + """Crea el QtInteractor y carga la malla Wigley por defecto.""" + if not _PYVISTA_OK: + return + try: + lo = self.layout() + # Eliminar placeholder + if lo.count() > 0: + old = lo.takeAt(0).widget() + if old: + old.hide() + old.deleteLater() + + # Crear el interactor PyVista embebido + self._plotter = QtInteractor( + parent=self, + auto_update=False, # sin polling continuo de GPU + off_screen=False, + ) + self._plotter.setObjectName("pyvistaInteractor") + lo.addWidget(self._plotter) + + # Configurar tema dark para que combine con la UI + self._plotter.set_background("#1a1d30") # viewportCanvas color + + # Cargar casco Wigley como geometría de bienvenida + self._load_default_wigley() + self._ready = True + logger.info("Viewer3DWidget: QtInteractor iniciado correctamente") + + except Exception as exc: + logger.error("Error al iniciar PyVista: %s", exc) + lbl = QLabel(f"Error visor 3D:\n{exc}") + lbl.setAlignment(Qt.AlignmentFlag.AlignCenter) + lbl.setObjectName("viewportPlaceholder") + self.layout().addWidget(lbl) + + def _load_default_wigley(self) -> None: + """Renderiza el casco Wigley de demostración al arrancar.""" + try: + from arshipdesign.core.hull import Hull + hull = Hull.from_wigley( + lpp=10.0, beam=1.5, draft=0.75, n_stations=40, n_waterlines=20 + ) + mesh = hull.to_mesh(n_u=50, n_v=25) + self._render_hull_mesh(mesh) + except Exception as exc: + logger.warning("No se pudo cargar Wigley por defecto: %s", exc) + + # ------------------------------------------------------------------ + # API pública + # ------------------------------------------------------------------ + + def load_hull(self, hull) -> None: + """Carga un objeto Hull en el visor. + + Parámetros + ---------- + hull : arshipdesign.core.hull.Hull + """ + if not self._ready or self._plotter is None: + logger.warning("Viewer3DWidget no listo — hull no cargado") + return + try: + mesh = hull.to_mesh() + self._render_hull_mesh(mesh) + except Exception as exc: + logger.error("Error al cargar Hull: %s", exc) + + def load_mesh(self, mesh: "pv.PolyData") -> None: + """Carga directamente una malla PyVista.""" + if not self._ready or self._plotter is None: + return + self._render_hull_mesh(mesh) + + def reset_camera(self) -> None: + """Resetea la cámara a la vista isométrica.""" + if self._plotter is not None: + self._plotter.reset_camera() + self._plotter.view_isometric() + + def clear(self) -> None: + """Limpia todos los actores de la escena.""" + if self._plotter is not None: + self._plotter.clear() + + @property + def is_ready(self) -> bool: + return self._ready + + # ------------------------------------------------------------------ + # Helpers internos + # ------------------------------------------------------------------ + + def _render_hull_mesh(self, mesh: "pv.PolyData") -> None: + """Renderiza una malla con el estilo marítimo del tema.""" + if self._plotter is None: + return + self._plotter.clear() + + # Casco principal — color acero naval + self._plotter.add_mesh( + mesh, + color="#3a6080", + smooth_shading=True, + show_edges=True, + edge_color="#4da8ff", + line_width=0.3, + opacity=0.92, + name="hull", + ) + + # Plano de flotación semi-transparente + try: + hull_pts = mesh.points + x_min, x_max = float(hull_pts[:, 0].min()), float(hull_pts[:, 0].max()) + y_max = float(np.abs(hull_pts[:, 1]).max()) * 1.05 + wl = pv.Plane( + center=((x_min + x_max) / 2, 0, 0), + direction=(0, 0, 1), + i_size=(x_max - x_min), + j_size=2 * y_max, + i_resolution=1, + j_resolution=1, + ) + self._plotter.add_mesh( + wl, color="#4da8ff", opacity=0.15, name="waterplane" + ) + except Exception: + pass # no crítico + + self._plotter.view_isometric() + self._plotter.reset_camera() + + def closeEvent(self, event) -> None: # type: ignore[override] + """Libera el contexto PyVista al cerrar.""" + if self._plotter is not None: + try: + self._plotter.close() + except Exception: + pass + super().closeEvent(event) diff --git a/tests/test_sprint1_geometry.py b/tests/test_sprint1_geometry.py new file mode 100644 index 0000000..d21a4ce --- /dev/null +++ b/tests/test_sprint1_geometry.py @@ -0,0 +1,355 @@ +""" +Tests de Sprint 1 — Geometría NURBS, modelos de casco, hidrostáticos Wigley. + +El casco Wigley tiene soluciones analíticas exactas para sus hidrostáticos +que se usan como valores de referencia para verificar la implementación numérica. + +Soluciones analíticas del casco Wigley: + y(ξ, ζ) = (B/2) · [1 - (2ξ/L)²] · [1 - (ζ/T)²] + con ξ ∈ [-L/2, L/2], ζ ∈ [-T, 0] + + Volumen de desplazamiento (derivación exacta): + V = B · ∫₋ᴸ/₂ᴸ/²[1-(2ξ/L)²]dξ · ∫₋ᵀ⁰[1-(ζ/T)²]dζ + = B · (2L/3) · (2T/3) + = 4·B·L·T / 9 + + Área plano de flotación (ζ=0 → f_ζ=1, ambas bandas): + Awp = 2 · ∫₋ᴸ/₂ᴸ/² (B/2) · [1-(2ξ/L)²] dξ + = B · (2L/3) = 2BL/3 + + Cb = V/(L·B·T) = 4/9 ≈ 0.4444 + + LCB: x = Lpp/2 por simetría longitudinal del casco Wigley. +""" +import math +import numpy as np +import pytest + +from arshipdesign.geometry.nurbs_curve import BSplineCurve +from arshipdesign.geometry.nurbs_surface import LoftedSurface +from arshipdesign.core.section import Section +from arshipdesign.core.offsets import OffsetsTable +from arshipdesign.core.hull import Hull + + +# --------------------------------------------------------------------------- +# Parámetros del casco Wigley de referencia +# --------------------------------------------------------------------------- + +LPP = 10.0 +BEAM = 1.5 +DRAFT = 0.75 +N_STA = 41 # alta resolución para precisión numérica +N_WL = 21 + + +@pytest.fixture(scope="module") +def wigley_offsets() -> OffsetsTable: + return OffsetsTable.from_wigley( + lpp=LPP, beam=BEAM, draft=DRAFT, + n_stations=N_STA, n_waterlines=N_WL + ) + + +@pytest.fixture(scope="module") +def wigley_hull() -> Hull: + return Hull.from_wigley( + lpp=LPP, beam=BEAM, draft=DRAFT, + n_stations=N_STA, n_waterlines=N_WL + ) + + +# --------------------------------------------------------------------------- +# 1. BSplineCurve +# --------------------------------------------------------------------------- + +class TestBSplineCurve: + def test_creation_2d(self): + pts = np.array([[0, 0], [1, 1], [2, 0.5], [3, 0]], dtype=float) + c = BSplineCurve(pts, degree=3) + assert c.degree == 3 + assert c.dim == 2 + + def test_creation_3d(self): + pts = np.random.rand(6, 3) + c = BSplineCurve(pts, degree=3) + assert c.dim == 3 + + def test_evaluate_scalar(self): + pts = np.array([[0, 0], [1, 1], [2, 0.5], [3, 0]], dtype=float) + c = BSplineCurve(pts, degree=3) + p = c.evaluate(0.5) + assert p.shape == (2,) + + def test_evaluate_array(self): + pts = np.array([[0, 0], [1, 1], [2, 0.5], [3, 0]], dtype=float) + c = BSplineCurve(pts, degree=3) + t = np.linspace(0, 1, 50) + pts_out = c.evaluate(t) + assert pts_out.shape == (50, 2) + + def test_endpoints_interpolated(self): + """Los extremos de la curva deben pasar por los puntos originales.""" + pts = np.array([[0.0, 0.0], [1.0, 2.0], [3.0, 1.5], [4.0, 0.0]]) + c = BSplineCurve(pts, degree=3) + p0 = c.evaluate(0.0) + p1 = c.evaluate(1.0) + np.testing.assert_allclose(p0, pts[0], atol=1e-8) + np.testing.assert_allclose(p1, pts[-1], atol=1e-8) + + def test_sample_returns_n_points(self): + pts = np.array([[0, 0], [1, 1], [2, 0.5], [3, 0]], dtype=float) + c = BSplineCurve(pts, degree=3) + sampled = c.sample(80) + assert sampled.shape == (80, 2) + + def test_arc_length_line(self): + """Línea recta: longitud = distancia euclidiana.""" + pts = np.array([[0, 0], [1, 0], [2, 0], [3, 0]], dtype=float) + c = BSplineCurve(pts, degree=1) + L = c.arc_length() + assert abs(L - 3.0) < 0.01 + + def test_too_few_points_raises(self): + pts = np.array([[0, 0], [1, 1]], dtype=float) + with pytest.raises(ValueError): + BSplineCurve(pts, degree=3) + + def test_wrong_shape_raises(self): + with pytest.raises(ValueError): + BSplineCurve(np.array([1, 2, 3]), degree=1) + + +# --------------------------------------------------------------------------- +# 2. Section +# --------------------------------------------------------------------------- + +class TestSection: + def _make_rect_section(self, half_b: float = 2.0, height: float = 3.0) -> Section: + """Sección rectangular: área = 2·B·H.""" + return Section( + station=5, x=5.0, + half_breadths=np.array([half_b, half_b]), + z_positions=np.array([0.0, height]), + ) + + def _make_triangular_section(self, half_b: float = 2.0, height: float = 3.0) -> Section: + """Sección triangular (quilla aguda): área = B·H. + Usa 5 puntos para que la integración numérica sea precisa. + """ + n = 5 + z = np.linspace(0.0, height, n) + y = (half_b / height) * z # lineal 0 → half_b + return Section(station=5, x=5.0, half_breadths=y, z_positions=z) + + def test_area_rectangular(self): + s = self._make_rect_section(half_b=1.0, height=2.0) + # Área completa (ambas bandas) = 2 * 1.0 * 2.0 = 4.0 + assert abs(s.area() - 4.0) < 1e-6 + + def test_area_triangular(self): + s = self._make_triangular_section(half_b=1.0, height=2.0) + # Área completa = 2 * (0.5 * 1.0 * 2.0) = 2.0 + assert abs(s.area() - 2.0) < 1e-6 + + def test_area_at_draft(self): + s = self._make_rect_section(half_b=1.0, height=4.0) + # Recortar a draft=2.0 → área = 2*1*2 = 4.0 + assert abs(s.area(draft=2.0) - 4.0) < 1e-6 + + def test_area_zero_below_keel(self): + s = self._make_rect_section() + assert s.area(draft=0.0) == 0.0 + + def test_centroid_z_rectangular(self): + s = self._make_rect_section(half_b=1.0, height=2.0) + # Centroide de rectángulo = H/2 = 1.0 + assert abs(s.centroid_z() - 1.0) < 1e-6 + + def test_centroid_z_triangular(self): + s = self._make_triangular_section(half_b=1.0, height=3.0) + # Centroide de sección triangular (y lineal, keel=0): 2H/3 + # ∫y·z dz / ∫y dz = (H²/3)/(H/2) = 2H/3 = 2.0 + assert abs(s.centroid_z() - 2.0) < 1e-4 + + def test_max_half_breadth(self): + s = self._make_rect_section(half_b=2.5) + assert abs(s.max_half_breadth() - 2.5) < 1e-9 + + def test_curve_returns_bspline(self): + s = self._make_rect_section() + c = s.curve() + assert isinstance(c, BSplineCurve) + + def test_repr(self): + s = self._make_rect_section() + assert "Section" in repr(s) + + +# --------------------------------------------------------------------------- +# 3. OffsetsTable — casco Wigley +# --------------------------------------------------------------------------- + +class TestOffsetsTable: + def test_shape(self, wigley_offsets): + ot = wigley_offsets + assert ot.data.shape == (N_STA, N_WL) + + def test_endpoints_zero(self, wigley_offsets): + """AP y FP deben tener semi-manga cero (extremos agudos).""" + ot = wigley_offsets + assert ot.data[0, :].max() < 1e-9, "AP no es cero" + assert ot.data[-1, :].max() < 1e-9, "FP no es cero" + + def test_keel_zero(self, wigley_offsets): + """En la quilla (z=0, ζ=-T) la semi-manga es cero.""" + ot = wigley_offsets + # Primera columna: z=0 → ζ = 0 - T = -T → f_zeta = 1 - (-T/T)² = 0 + assert ot.data[:, 0].max() < 1e-9, "Quilla no es cero" + + def test_max_half_breadth(self, wigley_offsets): + """El máximo debe ser B/2 en el midship en la línea de agua.""" + ot = wigley_offsets + assert abs(ot.max_half_breadth - BEAM / 2.0) < 1e-6 + + def test_to_sections(self, wigley_offsets): + sections = wigley_offsets.to_sections() + assert len(sections) == N_STA + for s in sections: + assert isinstance(s, Section) + + def test_interpolation(self, wigley_offsets): + """Verificar interpolación en punto conocido: midship, waterline.""" + ot = wigley_offsets + x_mid = LPP / 2.0 + z_wl = DRAFT + y_interp = ot.half_breadth(x_mid, z_wl) + # Wigley: y = B/2 * (1-0) * (1-1) = 0 → ζ = 0 - T = -T → f_zeta=0 + # En z=draft, ζ=0 → f_zeta = 1 - (0/T)² = 1 + # y = B/2 * 1 * 1 = B/2 + assert abs(y_interp - BEAM / 2.0) < 1e-3 + + def test_repr(self, wigley_offsets): + r = repr(wigley_offsets) + assert "OffsetsTable" in r + + +# --------------------------------------------------------------------------- +# 4. Hull — hidrostáticos del casco Wigley +# --------------------------------------------------------------------------- + +class TestHullHydrostatics: + """ + Tolerancias para verificación numérica: + - Volumen: ±1 % respecto al valor analítico + - LCB: ±1 % de Lpp + - Awp: ±1 % + - Cb: ±1 % + """ + + # Valores analíticos exactos del casco Wigley + # V = 2·(B/2)·∫[1-(2ξ/L)²]dξ · ∫[1-(ζ/T)²]dζ + # = B · (2L/3) · (2T/3) = 4BLT/9 + V_ANALYTIC = (4.0 * BEAM * LPP * DRAFT) / 9.0 + # Awp = 2·∫(B/2)·[1-(2ξ/L)²]dξ = B · 2L/3 = 2BL/3 + AWP_ANALYTIC = (2.0 * BEAM * LPP) / 3.0 + LCB_ANALYTIC = LPP / 2.0 # simetría longitudinal + CB_ANALYTIC = 4.0 / 9.0 # V/(LBT) = 4BLT/9 / (LBT) = 4/9 + + def test_volume_within_1pct(self, wigley_hull): + V = wigley_hull.volume_of_displacement() + err_pct = abs(V - self.V_ANALYTIC) / self.V_ANALYTIC * 100 + assert err_pct < 1.0, ( + f"Volumen {V:.4f} m³ difiere {err_pct:.2f}% del analítico {self.V_ANALYTIC:.4f}" + ) + + def test_waterplane_area_within_1pct(self, wigley_hull): + Awp = wigley_hull.waterplane_area() + err_pct = abs(Awp - self.AWP_ANALYTIC) / self.AWP_ANALYTIC * 100 + assert err_pct < 1.0, ( + f"Awp {Awp:.4f} m² difiere {err_pct:.2f}% del analítico {self.AWP_ANALYTIC:.4f}" + ) + + def test_lcb_symmetric(self, wigley_hull): + lcb = wigley_hull.lcb() + err = abs(lcb - self.LCB_ANALYTIC) + assert err < 0.1 * LPP, ( + f"LCB {lcb:.4f} m difiere {err:.4f} m del midship {self.LCB_ANALYTIC}" + ) + + def test_block_coefficient(self, wigley_hull): + Cb = wigley_hull.block_coefficient() + err_pct = abs(Cb - self.CB_ANALYTIC) / self.CB_ANALYTIC * 100 + assert err_pct < 1.0, ( + f"Cb {Cb:.4f} difiere {err_pct:.2f}% del analítico {self.CB_ANALYTIC:.4f}" + ) + + def test_vcb_positive(self, wigley_hull): + kb = wigley_hull.vcb() + assert 0.0 < kb < DRAFT, f"VCB {kb:.4f} fuera del rango [0, T={DRAFT}]" + + def test_midship_coefficient_less_than_1(self, wigley_hull): + cm = wigley_hull.midship_coefficient() + assert 0.0 < cm <= 1.0, f"Cm={cm:.4f} fuera del rango (0,1]" + + def test_displacement_positive(self, wigley_hull): + D = wigley_hull.displacement_tonnes() + assert D > 0.0 + + def test_displacement_freshwater_less(self, wigley_hull): + D_sw = wigley_hull.displacement_tonnes(rho=1025.0) + D_fw = wigley_hull.displacement_tonnes(rho=1000.0) + assert D_sw > D_fw + + def test_repr(self, wigley_hull): + assert "Hull" in repr(wigley_hull) + + +# --------------------------------------------------------------------------- +# 5. LoftedSurface básico +# --------------------------------------------------------------------------- + +class TestLoftedSurface: + def _make_simple_surface(self) -> LoftedSurface: + """Superficie sencilla: 5 secciones rectangulares.""" + sections = [] + for i, u in enumerate(np.linspace(0.0, 1.0, 6)): + # Semi-manga varía de 0 a 1 a 0 (forma Wigley simplificada) + y_max = 1.0 - abs(2 * u - 1.0) ** 2 + pts = np.array([[0.0, 0.0], [y_max, 0.5], [y_max, 1.0]]) + sections.append((u, pts)) + return LoftedSurface(sections, degree_u=3, degree_v=2) + + def test_creation(self): + surf = self._make_simple_surface() + assert surf.n_sections == 6 + + def test_evaluate_scalar(self): + surf = self._make_simple_surface() + pt = surf.evaluate(0.5, 0.5) + assert pt.shape == (2,) + + def test_section_at_returns_curve(self): + surf = self._make_simple_surface() + curve = surf.section_at(0.5) + assert isinstance(curve, BSplineCurve) + + def test_insufficient_sections_raises(self): + pts = np.array([[0.0, 0.0], [0.5, 1.0], [0.0, 2.0]]) + with pytest.raises(ValueError): + LoftedSurface([(0.0, pts), (0.5, pts)], degree_u=3) + + +# --------------------------------------------------------------------------- +# 6. Wigley hull mesh (si PyVista disponible) +# --------------------------------------------------------------------------- + +def test_wigley_hull_mesh(wigley_hull): + """Malla PyVista del casco Wigley debe tener puntos válidos.""" + pv = pytest.importorskip("pyvista", reason="PyVista no instalado") + mesh = wigley_hull.to_mesh(n_u=20, n_v=10) + assert mesh.n_points > 0 + assert mesh.n_cells > 0 + # Sin puntos NaN + pts = mesh.points + assert not np.any(np.isnan(pts)), "La malla contiene puntos NaN"