765 lines
32 KiB
Python
765 lines
32 KiB
Python
"""
|
||
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)
|
||
# Altura del sheer (cubierta) por estación [m].
|
||
# Permite líneas de cubierta con arrufo/caída sin alterar el puntal escalar.
|
||
# Default vacío → se usa hull.depth uniforme en todas las estaciones.
|
||
sheer_z: np.ndarray = field(default_factory=lambda: np.array([]), repr=False, compare=False)
|
||
# Puntos de control de la roda (stem) en el plano X-Z — shape (n, 2).
|
||
# Default vacío → 3 puntos colineales de quilla-FP a cubierta-FP.
|
||
stem_ctrl: np.ndarray = field(default_factory=lambda: np.zeros((0, 2)), repr=False, compare=False)
|
||
# Puntos de control del contorno del espejo (AP) en el plano X-Z — shape (n, 2).
|
||
# Default vacío → 3 puntos colineales de quilla-AP a cubierta-AP.
|
||
transom_ctrl: np.ndarray = field(default_factory=lambda: np.zeros((0, 2)), repr=False, compare=False)
|
||
# Desviaciones X per-estación para los nodos de quilla y sheer [m].
|
||
# X efectiva quilla(i) = offsets.x_stations[i] + keel_x_offsets[i].
|
||
keel_x_offsets: np.ndarray = field(default_factory=lambda: np.array([]), repr=False, compare=False)
|
||
sheer_x_offsets: np.ndarray = field(default_factory=lambda: np.array([]), repr=False, compare=False)
|
||
# Posiciones X de los planos de estación para la visualización [m].
|
||
# Independiente de la malla paramétrica (x_stations).
|
||
# Default vacío → se generan n=11 planos uniformes entre 0 y Lpp.
|
||
station_planes: np.ndarray = field(default_factory=lambda: np.array([]), repr=False, compare=False)
|
||
# Nodos marcados como esquina: rompen la suavidad B-spline en ese punto.
|
||
# Cada elemento es [i, j] donde j puede ser _KEEL_IDX(-1), _SHEER_IDX(-2)
|
||
# o un índice de línea de agua (0 .. n_wl-1).
|
||
corner_nodes: list = field(default_factory=list, 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 get_stem_ctrl(self) -> np.ndarray:
|
||
"""Puntos de control de la roda en X-Z.
|
||
|
||
El primer punto siempre coincide con la quilla en FP y el último
|
||
con el sheer en FP, garantizando continuidad del contorno del perfil.
|
||
Solo los puntos intermedios son libremente editables.
|
||
"""
|
||
x_fp = float(self.offsets.x_stations[-1])
|
||
z0 = float(self.offsets.keel_z[-1]) # quilla en FP
|
||
z1 = float(self.get_sheer_z()[-1]) # sheer en FP
|
||
if self.stem_ctrl.shape[0] >= 3:
|
||
ctrl = self.stem_ctrl.copy()
|
||
ctrl[0] = [x_fp, z0] # snap inferior → quilla-FP
|
||
ctrl[-1] = [x_fp, z1] # snap superior → sheer-FP
|
||
return ctrl
|
||
# Default: roda vertical con 3 puntos
|
||
return np.array([[x_fp, z0], [x_fp, (z0 + z1) * 0.5], [x_fp, z1]])
|
||
|
||
def get_transom_ctrl(self) -> np.ndarray:
|
||
"""Puntos de control del contorno del espejo (AP) en X-Z.
|
||
|
||
Orden: [0] = quilla-AP, [-1] = sheer-AP.
|
||
El primer y último punto siempre están fijados a quilla/sheer en AP.
|
||
Solo los puntos intermedios son libremente editables.
|
||
"""
|
||
x_ap = float(self.offsets.x_stations[0])
|
||
z0 = float(self.offsets.keel_z[0]) # quilla en AP
|
||
z1 = float(self.get_sheer_z()[0]) # sheer en AP
|
||
if self.transom_ctrl.shape[0] >= 3:
|
||
ctrl = self.transom_ctrl.copy()
|
||
ctrl[0] = [x_ap, z0] # snap inferior → quilla-AP
|
||
ctrl[-1] = [x_ap, z1] # snap superior → sheer-AP
|
||
return ctrl
|
||
return np.array([[x_ap, z0], [x_ap, (z0 + z1) * 0.5], [x_ap, z1]])
|
||
|
||
def get_sheer_z(self) -> np.ndarray:
|
||
"""Alturas de cubierta (sheer) por estación [m].
|
||
|
||
Si ``sheer_z`` no está inicializado o tiene dimensión incorrecta,
|
||
devuelve un array uniforme con el valor del puntal (``depth``).
|
||
"""
|
||
n = self.offsets.n_stations
|
||
if len(self.sheer_z) == n:
|
||
return self.sheer_z
|
||
return np.full(n, self.depth)
|
||
|
||
def get_keel_x_offsets(self) -> np.ndarray:
|
||
"""Desviaciones X per-estación para los nodos de quilla [m]."""
|
||
n = self.offsets.n_stations
|
||
if len(self.keel_x_offsets) == n:
|
||
return self.keel_x_offsets
|
||
return np.zeros(n)
|
||
|
||
def get_sheer_x_offsets(self) -> np.ndarray:
|
||
"""Desviaciones X per-estación para los nodos de sheer [m]."""
|
||
n = self.offsets.n_stations
|
||
if len(self.sheer_x_offsets) == n:
|
||
return self.sheer_x_offsets
|
||
return np.zeros(n)
|
||
|
||
def is_corner(self, i: int, j: int) -> bool:
|
||
"""Indica si el nodo (i, j) está marcado como esquina."""
|
||
return any(int(e[0]) == i and int(e[1]) == j for e in self.corner_nodes)
|
||
|
||
def toggle_corner(self, i: int, j: int) -> None:
|
||
"""Alterna el estado de esquina del nodo (i, j)."""
|
||
if self.is_corner(i, j):
|
||
self.corner_nodes = [e for e in self.corner_nodes
|
||
if not (int(e[0]) == i and int(e[1]) == j)]
|
||
else:
|
||
self.corner_nodes.append([i, j])
|
||
|
||
def get_station_planes(self, n: int = 11) -> np.ndarray:
|
||
"""Posiciones X [m] de los planos de estación para visualización.
|
||
|
||
Si no están configurados, devuelve n planos uniformes entre 0 y Lpp.
|
||
"""
|
||
if len(self.station_planes) >= 2:
|
||
return self.station_planes
|
||
return np.linspace(0.0, self.lpp, n)
|
||
|
||
def insert_station(self, x: float) -> None:
|
||
"""Inserta una nueva estación en la posición x [m], interpolando todos los arrays.
|
||
|
||
Actualiza OffsetsTable (x_stations, data, keel_z) y Hull (sheer_z).
|
||
AP y FP no se pueden sobrepasar.
|
||
"""
|
||
ot = self.offsets
|
||
x = float(np.clip(x, float(ot.x_stations[0]) + 1e-3, float(ot.x_stations[-1]) - 1e-3))
|
||
idx = int(np.searchsorted(ot.x_stations, x))
|
||
|
||
old_x = ot.x_stations.copy()
|
||
old_sheer = self.get_sheer_z().copy()
|
||
|
||
# Interpolar semi-mangas y desviaciones para la nueva estación
|
||
new_y = np.array([
|
||
float(np.interp(x, old_x, ot.data[:, j]))
|
||
for j in range(ot.n_waterlines)
|
||
])
|
||
new_keel_z = float(np.interp(x, old_x, ot.keel_z))
|
||
new_sheer_z = float(np.interp(x, old_x, old_sheer))
|
||
new_z_offsets = np.array([
|
||
float(np.interp(x, old_x, ot.z_offsets[:, j]))
|
||
for j in range(ot.n_waterlines)
|
||
])
|
||
new_x_offsets = np.array([
|
||
float(np.interp(x, old_x, ot.x_offsets[:, j]))
|
||
for j in range(ot.n_waterlines)
|
||
])
|
||
new_keel_x_off = float(np.interp(x, old_x, self.get_keel_x_offsets()))
|
||
new_sheer_x_off = float(np.interp(x, old_x, self.get_sheer_x_offsets()))
|
||
|
||
ot.x_stations = np.insert(old_x, idx, x)
|
||
ot.data = np.insert(ot.data, idx, new_y, axis=0)
|
||
ot.keel_z = np.insert(ot.keel_z, idx, new_keel_z)
|
||
ot.z_offsets = np.insert(ot.z_offsets, idx, new_z_offsets, axis=0)
|
||
ot.x_offsets = np.insert(ot.x_offsets, idx, new_x_offsets, axis=0)
|
||
lbl = f"S{idx}"
|
||
ot.station_labels.insert(idx, lbl)
|
||
|
||
self.sheer_z = np.insert(old_sheer, idx, new_sheer_z)
|
||
self.keel_x_offsets = np.insert(self.get_keel_x_offsets(), idx, new_keel_x_off)
|
||
self.sheer_x_offsets = np.insert(self.get_sheer_x_offsets(), idx, new_sheer_x_off)
|
||
self.invalidate()
|
||
|
||
def insert_waterline(self, z: float) -> None:
|
||
"""Inserta una nueva línea de agua en altura z [m], interpolando semi-mangas.
|
||
|
||
No afecta keel_z ni sheer_z (son arrays por estación, no por LdA).
|
||
"""
|
||
ot = self.offsets
|
||
z = float(np.clip(z, float(ot.z_waterlines[0]) + 1e-3, float(ot.z_waterlines[-1]) - 1e-3))
|
||
idx = int(np.searchsorted(ot.z_waterlines, z))
|
||
new_y = np.array([
|
||
float(np.interp(z, ot.z_waterlines, ot.data[i, :]))
|
||
for i in range(ot.n_stations)
|
||
])
|
||
ot.z_waterlines = np.insert(ot.z_waterlines, idx, z)
|
||
ot.data = np.insert(ot.data, idx, new_y, axis=1)
|
||
ot.z_offsets = np.insert(ot.z_offsets, idx, 0.0, axis=1)
|
||
ot.x_offsets = np.insert(ot.x_offsets, idx, 0.0, axis=1)
|
||
self.invalidate()
|
||
|
||
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 snap_boundary_nodes_to_contours(self) -> None:
|
||
"""Enclava los nodos extremos de cada línea de agua sobre las
|
||
aristas de terminación del casco: roda (FP, índice -1) y espejo
|
||
(AP, índice 0).
|
||
|
||
El stem y el transom son **aristas de terminación** (boundary edges)
|
||
— no son curvas interiores de suavizado. Cada línea de agua debe
|
||
terminar exactamente donde intersecta esa arista. Este método
|
||
calcula la coordenada X de esa intersección para cada altura Z y
|
||
actualiza x_offsets[0,:] y x_offsets[-1,:] en consecuencia.
|
||
|
||
Llamar tras:
|
||
- Cargar un casco desde disco (``from_dict``).
|
||
- Arrastrar un punto de control del stem o del transom.
|
||
"""
|
||
from arshipdesign.geometry.nurbs_curve import BSplineCurve
|
||
|
||
ot = self.offsets
|
||
n_wl = ot.n_waterlines
|
||
NSAMP = max(300, n_wl * 30)
|
||
|
||
def _sample_edge(ctrl: np.ndarray) -> np.ndarray:
|
||
"""Muestrea la curva spline de una arista de terminación.
|
||
|
||
Devuelve (NSAMP, 2) [X, Z] ordenado por Z creciente para
|
||
poder usar np.interp(z, ...) de forma segura.
|
||
"""
|
||
m = len(ctrl)
|
||
if m < 2:
|
||
return ctrl[np.argsort(ctrl[:, 1])]
|
||
try:
|
||
deg = min(3, m - 1)
|
||
curve = BSplineCurve(ctrl, degree=deg)
|
||
pts = curve.sample(NSAMP) # (NSAMP, 2): X, Z
|
||
except Exception:
|
||
t_raw = np.linspace(0.0, 1.0, m)
|
||
t_new = np.linspace(0.0, 1.0, NSAMP)
|
||
pts = np.column_stack([
|
||
np.interp(t_new, t_raw, ctrl[:, 0]),
|
||
np.interp(t_new, t_raw, ctrl[:, 1]),
|
||
])
|
||
return pts[np.argsort(pts[:, 1])] # ordenar por Z
|
||
|
||
stem_samp = _sample_edge(self.get_stem_ctrl()) # FP (i = -1)
|
||
trans_samp = _sample_edge(self.get_transom_ctrl()) # AP (i = 0)
|
||
|
||
# Asegurar forma correcta
|
||
if ot.x_offsets.shape != (ot.n_stations, n_wl):
|
||
ot.x_offsets = np.zeros((ot.n_stations, n_wl))
|
||
|
||
z_stem_min, z_stem_max = float(stem_samp[:, 1].min()), float(stem_samp[:, 1].max())
|
||
z_trans_min, z_trans_max = float(trans_samp[:, 1].min()), float(trans_samp[:, 1].max())
|
||
|
||
for j in range(n_wl):
|
||
z_fp = float(ot.z_waterlines[j]) + float(ot.z_offsets[-1, j])
|
||
z_ap = float(ot.z_waterlines[j]) + float(ot.z_offsets[0, j])
|
||
|
||
# Clamp dentro del rango de la arista (evita extrapolación)
|
||
z_fp = float(np.clip(z_fp, z_stem_min, z_stem_max))
|
||
z_ap = float(np.clip(z_ap, z_trans_min, z_trans_max))
|
||
|
||
x_on_stem = float(np.interp(z_fp, stem_samp[:, 1], stem_samp[:, 0]))
|
||
x_on_trans = float(np.interp(z_ap, trans_samp[:, 1], trans_samp[:, 0]))
|
||
|
||
# x_offsets = X_efectiva − X_referencia_paramétrica
|
||
ot.x_offsets[-1, j] = x_on_stem - float(ot.x_stations[-1])
|
||
ot.x_offsets[0, j] = x_on_trans - float(ot.x_stations[0])
|
||
|
||
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 + self.offsets.z_offsets[i, :],
|
||
])
|
||
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).
|
||
|
||
Cada sección se construye desde la quilla (keel_z per-estación)
|
||
hasta la cubierta (sheer_z per-estación), igual que BodyPlanViewer.
|
||
Requiere PyVista instalado. Retorna un PolyData triangulado.
|
||
"""
|
||
try:
|
||
import pyvista as pv
|
||
except ImportError as exc:
|
||
raise ImportError("PyVista no está instalado") from exc
|
||
|
||
from arshipdesign.geometry.nurbs_curve import BSplineCurve
|
||
|
||
ot = self.offsets
|
||
sheer = self.get_sheer_z()
|
||
|
||
# ── Paso 1: perfiles (y, z) a las estaciones originales ──────────
|
||
def _section_yz(i: int, n_pts: int) -> np.ndarray:
|
||
"""Muestrea n_pts puntos del perfil (y, z) de la sección i."""
|
||
kz = float(ot.keel_z[i])
|
||
sz = float(sheer[i])
|
||
y_arr = ot.data[i, :]
|
||
z_arr = ot.z_waterlines + ot.z_offsets[i, :]
|
||
keel_pt = np.array([[0.0, kz]])
|
||
sheer_pt = np.array([[float(y_arr[-1]), sz]])
|
||
raw = np.vstack([keel_pt, np.column_stack([y_arr, z_arr]), sheer_pt])
|
||
m = len(raw)
|
||
try:
|
||
k = min(3, max(m - 1, 1))
|
||
curve = BSplineCurve(raw, degree=k)
|
||
return curve.sample(n_pts) # (n_pts, 2)
|
||
except Exception:
|
||
# Fallback: linear re-sample of the raw polyline
|
||
t_raw = np.linspace(0, 1, m)
|
||
t_new = np.linspace(0, 1, n_pts)
|
||
return np.column_stack([
|
||
np.interp(t_new, t_raw, raw[:, 0]),
|
||
np.interp(t_new, t_raw, raw[:, 1]),
|
||
])
|
||
|
||
n_sta = ot.n_stations
|
||
profiles = np.array([_section_yz(i, n_v) for i in range(n_sta)])
|
||
# profiles shape: (n_sta, n_v, 2) — (y, z) at each station
|
||
|
||
# ── Paso 2: interpolar a n_u estaciones uniformes ────────────────
|
||
x_orig = ot.x_stations / self.lpp # normalized [0,1]
|
||
x_new = np.linspace(0.0, 1.0, n_u)
|
||
|
||
yz_grid = np.zeros((n_u, n_v, 2))
|
||
for v_idx in range(n_v):
|
||
for c in range(2): # y, z
|
||
yz_grid[:, v_idx, c] = np.interp(x_new, x_orig, profiles[:, v_idx, c])
|
||
|
||
# ── Paso 3: construir vértices 3D ─────────────────────────────────
|
||
x_grid = (x_new[:, None] * self.lpp) * np.ones((n_u, n_v))
|
||
y_grid = yz_grid[:, :, 0] # semi-manga (estribor +)
|
||
z_grid = yz_grid[:, :, 1]
|
||
|
||
pts_stbd = np.stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()], axis=1)
|
||
pts_port = np.stack([x_grid.ravel(), -y_grid.ravel(), z_grid.ravel()], axis=1)
|
||
all_pts = np.vstack([pts_stbd, pts_port])
|
||
|
||
# ── Paso 4: 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])
|
||
|
||
mesh = pv.PolyData(all_pts, np.array(faces, dtype=int))
|
||
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,
|
||
"sheer_z": self.get_sheer_z().tolist(),
|
||
"stem_ctrl": self.get_stem_ctrl().tolist(),
|
||
"transom_ctrl": self.get_transom_ctrl().tolist(),
|
||
"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)
|
||
"keel_z": ot.keel_z.tolist(),
|
||
"z_offsets": ot.z_offsets.tolist(), # (n_sta, n_wl)
|
||
"x_offsets": ot.x_offsets.tolist(), # (n_sta, n_wl)
|
||
},
|
||
"keel_x_offsets": self.get_keel_x_offsets().tolist(),
|
||
"sheer_x_offsets": self.get_sheer_x_offsets().tolist(),
|
||
"station_planes": self.get_station_planes().tolist(),
|
||
"corner_nodes": self.corner_nodes,
|
||
}
|
||
|
||
@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"]),
|
||
keel_z = np.array(od.get("keel_z", []), dtype=float),
|
||
z_offsets = np.array(od.get("z_offsets", np.zeros((0, 0))), dtype=float),
|
||
x_offsets = np.array(od.get("x_offsets", np.zeros((0, 0))), dtype=float),
|
||
)
|
||
hull = cls(
|
||
name = str(data["name"]),
|
||
lpp = float(data["lpp"]),
|
||
beam = float(data["beam"]),
|
||
depth = float(data["depth"]),
|
||
draft = float(data["draft"]),
|
||
offsets = offsets,
|
||
sheer_z = np.array(data.get("sheer_z", []), dtype=float),
|
||
stem_ctrl = np.array(data.get("stem_ctrl", []), dtype=float).reshape(-1, 2),
|
||
transom_ctrl = np.array(data.get("transom_ctrl", []), dtype=float).reshape(-1, 2),
|
||
keel_x_offsets = np.array(data.get("keel_x_offsets", []), dtype=float),
|
||
sheer_x_offsets = np.array(data.get("sheer_x_offsets", []), dtype=float),
|
||
station_planes = np.array(data.get("station_planes", []), dtype=float),
|
||
corner_nodes = list(data.get("corner_nodes", [])),
|
||
)
|
||
# NOTA: NO snap_boundary_nodes_to_contours() al cargar.
|
||
# Los x_offsets guardados son la posición que el usuario definió
|
||
# en el Visor Perfil (eje X libre) y se restauran tal cual.
|
||
# El snap solo aplica en la creación inicial del proyecto (wizard).
|
||
return hull
|
||
|
||
# ------------------------------------------------------------------
|
||
# Dunder
|
||
# ------------------------------------------------------------------
|
||
|
||
def __repr__(self) -> str:
|
||
return (
|
||
f"Hull({self.name!r}, Lpp={self.lpp} m, B={self.beam} m, "
|
||
f"T={self.draft} m)"
|
||
)
|