AR-Shipdesign/arshipdesign/core/hull.py

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