""" 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 invalidate(self) -> None: """Invalida la caché de la superficie NURBS. Llamar siempre que se modifiquen los offsets in-place (p.ej. arrastre interactivo en los visores 2D) para que la próxima llamada a ``surface`` o ``to_mesh`` reconstruya la geometría desde los datos actualizados. """ self._surface = None 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 def waterplane_coefficient(self, draft: Optional[float] = None) -> float: """Coeficiente de plano de flotación Cw = Awp / (Lpp · B). IACS Rec.34 §3.3 — parámetro adimensional de la forma del plano de flotación. """ T = draft if draft is not None else self.draft awp = self.waterplane_area(T) return awp / (self.lpp * self.beam) def it_waterplane(self, draft: Optional[float] = None) -> float: """Segundo momento de área del plano de flotación sobre el eje de crujía IT [m⁴]. IT = (2/3) · ∫₀^L y(x,T)³ dx Rawson & Tupper, "Basic Ship Theory" 5ª ed., Cap. 3. """ 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]) integrand = (2.0 / 3.0) * y_wl ** 3 if len(x) >= 3: return abs(float(simpson(integrand, x=x))) return abs(float(np.trapz(integrand, x))) def il_waterplane(self, draft: Optional[float] = None) -> float: """Segundo momento de área del plano de flotación sobre el centro de flotación IL [m⁴]. IL = ∫₀^L (x − LCF)² · 2y(x,T) dx Rawson & Tupper, "Basic Ship Theory" 5ª ed., Cap. 3. """ 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]) strip = 2.0 * y_wl if len(x) >= 3: awp = float(simpson(strip, x=x)) if awp > 1e-12: lcf = float(simpson(strip * x, x=x)) / awp else: lcf = self.lpp / 2.0 return abs(float(simpson(strip * (x - lcf) ** 2, x=x))) awp = float(np.trapz(strip, x)) lcf = float(np.trapz(strip * x, x)) / awp if awp > 1e-12 else self.lpp / 2.0 return abs(float(np.trapz(strip * (x - lcf) ** 2, x))) def bm_transverse(self, draft: Optional[float] = None) -> float: """Radio metacéntrico transversal BM_T = IT / V [m].""" T = draft if draft is not None else self.draft vol = self.volume_of_displacement(T) return self.it_waterplane(T) / vol if vol > 1e-12 else 0.0 def bm_longitudinal(self, draft: Optional[float] = None) -> float: """Radio metacéntrico longitudinal BM_L = IL / V [m].""" T = draft if draft is not None else self.draft vol = self.volume_of_displacement(T) return self.il_waterplane(T) / vol if vol > 1e-12 else 0.0 def km_transverse(self, draft: Optional[float] = None) -> float: """Altura del metacentro transversal KM_T = KB + BM_T [m]. Rawson & Tupper, "Basic Ship Theory" 5ª ed., §3.2. """ T = draft if draft is not None else self.draft return self.vcb(T) + self.bm_transverse(T) def tpc(self, draft: Optional[float] = None, rho: float = 1025.0) -> float: """Toneladas por centímetro de inmersión TPC [t/cm]. TPC = Awp · ρ / 100 000 Equivale a la masa añadida necesaria para aumentar el calado 1 cm. """ T = draft if draft is not None else self.draft return self.waterplane_area(T) * rho / 100_000.0 def mct1cm( self, draft: Optional[float] = None, rho: float = 1025.0, kg: Optional[float] = None, ) -> float: """Momento para cambiar asiento 1 cm MCT [t·m/cm]. MCT = Δ · GM_L / (100 · Lpp) GM_L = KB + BM_L − KG Si *kg* es None se usa la estimación KG ≈ depth × 0.55 (válida para embarcaciones con DWT vacío sin peso de carga). """ T = draft if draft is not None else self.draft if kg is None: kg = self.depth * 0.55 gml = max(self.vcb(T) + self.bm_longitudinal(T) - kg, 0.0) delta = self.displacement_tonnes(T, rho) return delta * gml / (100.0 * self.lpp) # ------------------------------------------------------------------ # 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() # ------------------------------------------------------------------ # Serialización JSON (.arsd) # ------------------------------------------------------------------ def to_dict(self) -> dict: """Serializa el Hull a un diccionario JSON-serializable. Formato interno: ``hull_v1``. Los arrays numpy se convierten a listas de Python para compatibilidad con json.dumps sin dependencias adicionales. IACS Rec.34 §6 — trazabilidad de datos de entrada (offsets guardados fielmente con la precisión de la tabla original). """ ot = self.offsets return { "format": "hull_v1", "name": self.name, "lpp": self.lpp, "beam": self.beam, "depth": self.depth, "draft": self.draft, "offsets": { "lpp": ot.lpp, "beam": ot.beam, "draft": ot.draft, "x_stations": ot.x_stations.tolist(), "z_waterlines": ot.z_waterlines.tolist(), "station_labels": list(ot.station_labels), "data": ot.data.tolist(), # (n_sta, n_wl) }, } @classmethod def from_dict(cls, data: dict) -> "Hull": """Deserializa un Hull desde un diccionario (leído de un archivo .arsd). Compatible con los formatos ``hull_v1`` y datos heredados sin versión. Parameters ---------- data : dict Diccionario generado por ``Hull.to_dict()``. Raises ------ KeyError Si faltan campos obligatorios. ValueError Si las dimensiones de la tabla son inconsistentes. """ od = data["offsets"] offsets = OffsetsTable( x_stations = np.array(od["x_stations"], dtype=float), z_waterlines = np.array(od["z_waterlines"], dtype=float), data = np.array(od["data"], dtype=float), station_labels = od.get("station_labels", []), lpp = float(od["lpp"]), beam = float(od["beam"]), draft = float(od["draft"]), ) return cls( name = str(data["name"]), lpp = float(data["lpp"]), beam = float(data["beam"]), depth = float(data["depth"]), draft = float(data["draft"]), offsets = offsets, ) # ------------------------------------------------------------------ # Dunder # ------------------------------------------------------------------ def __repr__(self) -> str: return ( f"Hull({self.name!r}, Lpp={self.lpp} m, B={self.beam} m, " f"T={self.draft} m)" )