# Copyright (C) 2025 Antoine ANCEAU
#
# This file is part of pynavaltoolbox.
#
# pynavaltoolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
Hydrostatic calculation tools.
This module provides utility functions for hydrostatic calculations,
including volume and surface area computations.
"""
import vtk
import numpy as np
from typing import Tuple, Callable
[docs]
def calculate_submerged_properties(
polydata: vtk.vtkPolyData,
surface_function: Callable[[float, float], float],
surface_type: str = 'planar',
draft: float = 0.0,
**kwargs
) -> Tuple[float, np.ndarray, float, Tuple[float, float, float, float, float, float]]:
"""
Calculate submerged volume, center of buoyancy, and wetted surface area.
This is a generic function that can handle different types of free surfaces:
- Planar surface (calm water)
- Sinusoidal surface (waves)
- Custom surface defined by a function
Args:
polydata: The hull geometry as vtkPolyData.
surface_function: Function f(x, y) -> z defining the free surface.
For planar: lambda x, y: draft.
For sinusoidal: lambda x, y: draft + A*sin(k*x + phi).
surface_type: Type of surface ('planar', 'sinusoidal', 'custom').
draft: Reference draft level (used for planar surfaces).
**kwargs: Additional parameters (e.g., wave amplitude, wavelength).
Returns:
Tuple of (volume, center_of_buoyancy, wetted_surface_area, bounds) where:
volume: Submerged volume in m³.
center_of_buoyancy: numpy array [x, y, z] in meters.
wetted_surface_area: Wetted surface area in m².
bounds: (xmin, xmax, ymin, ymax, zmin, zmax) of submerged part.
"""
if surface_type == 'planar':
return _calculate_planar_submersion(polydata, draft)
elif surface_type == 'sinusoidal':
# Future implementation for wave surfaces
raise NotImplementedError("Sinusoidal surfaces not yet implemented")
else:
# Future implementation for custom surfaces
raise NotImplementedError("Custom surfaces not yet implemented")
def _calculate_planar_submersion(
polydata: vtk.vtkPolyData,
draft: float
) -> Tuple[float, np.ndarray, float, Tuple[float, float, float, float, float, float]]:
"""
Calculate properties for a planar free surface (calm water).
Args:
polydata: The hull geometry.
draft: Draft level (Z coordinate of waterline).
Returns:
Tuple of (volume, center_of_buoyancy, wetted_surface_area)
"""
# Get hull bounds to check if it's fully submerged or fully emerged
hull_bounds = polydata.GetBounds()
z_min = hull_bounds[4]
z_max = hull_bounds[5]
# Check if hull is completely above water (fully emerged)
if z_min >= draft:
# Hull is completely out of water - no submerged volume
return 0.0, np.array([0.0, 0.0, 0.0]), 0.0, (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
# Check if hull is completely below water (fully submerged)
if z_max <= draft:
# Hull is completely under water - use full hull volume
# Triangulate the hull
triangulate = vtk.vtkTriangleFilter()
triangulate.SetInputData(polydata)
triangulate.Update()
triangulated = triangulate.GetOutput()
# Clean polydata - only polys, no lines or verts
clean_polys = vtk.vtkPolyData()
clean_polys.SetPoints(triangulated.GetPoints())
clean_polys.SetPolys(triangulated.GetPolys())
# Calculate total volume
mass_props = vtk.vtkMassProperties()
mass_props.SetInputData(clean_polys)
mass_props.Update()
volume = abs(mass_props.GetVolume())
# Calculate total wetted surface area (entire hull surface)
wetted_area = _calculate_surface_area(polydata)
# Calculate center of buoyancy (centroid of entire hull)
center = _calculate_volumetric_centroid(clean_polys)
# Bounds are hull bounds
bounds = hull_bounds
return volume, center, wetted_area, bounds
# Hull is partially submerged - proceed with clipping
# Create a plane at the waterline
plane = vtk.vtkPlane()
plane.SetOrigin(0, 0, draft)
plane.SetNormal(0, 0, 1)
# Clip the hull at the waterline
clipper = vtk.vtkClipPolyData()
clipper.SetInputData(polydata)
clipper.SetClipFunction(plane)
clipper.InsideOutOn() # Keep the part below the plane
clipper.Update()
clipped = clipper.GetOutput()
# Get the waterline contour
cutter = vtk.vtkCutter()
cutter.SetInputData(polydata)
cutter.SetCutFunction(plane)
cutter.Update()
waterline = cutter.GetOutput()
# Create a cap at the waterline
if waterline.GetNumberOfPoints() > 2:
# Use vtkStripper to create polylines from the cutter output
stripper = vtk.vtkStripper()
stripper.SetInputData(waterline)
stripper.Update()
# internal triangulation of the contour
# This requires the contour to be a closed loop (or loops)
# vtkContourTriangulator is robust for this
triangulator = vtk.vtkContourTriangulator()
triangulator.SetInputData(stripper.GetOutput())
triangulator.Update()
cap_polydata = triangulator.GetOutput()
# Orient the cap normals to point downwards (into the water/hull)
# vtkContourTriangulator usually produces upward normals for CCW contours
# We need to check and reverse if needed, but for volume calc,
# usually outward normals are needed.
# If the clipped hull has outward normals (pointing into water),
# the cap should have outward normals (pointing up, out of water).
# Wait, Clipping keeps the 'inside' (Implicit function < 0).
# For a plane with normal (0,0,1), 'inside' is below the plane. This is correct.
# The clipped hull surface normals point OUT of the hull (into water).
# To close the volume, the cap normals should point UP (out of the hull).
# Verify orientation or ensure consistency.
# vtkMassProperties expects valid closed surface with consistent outward normals.
# We'll rely on vtkCleanPolyData and maybe vtkReverseSense if volume comes out negative.
# Combine the clipped hull with the cap
append = vtk.vtkAppendPolyData()
append.AddInputData(clipped)
append.AddInputData(cap_polydata)
append.Update()
combined = append.GetOutput()
# Clean the mesh to remove duplicate points and merge edges
cleaner = vtk.vtkCleanPolyData()
cleaner.SetInputData(combined)
cleaner.SetTolerance(1e-4) # Small tolerance to merge matching points
cleaner.Update()
closed_mesh = cleaner.GetOutput()
# Ensure consistent normals
normals = vtk.vtkPolyDataNormals()
normals.SetInputData(closed_mesh)
normals.ConsistencyOn()
normals.SplittingOff()
normals.AutoOrientNormalsOn() # Ensure outward normals
normals.Update()
closed_mesh = normals.GetOutput()
else:
# If no waterline contour is generated, it means the hull is either fully submerged
# or fully emerged, which should have been caught by the initial bounds check.
# In this specific path (partially submerged), this 'else' block should ideally not be reached.
# However, if it is, it implies the clipped mesh itself forms a closed volume or is the only relevant part.
closed_mesh = clipped
# Triangulate
triangulate = vtk.vtkTriangleFilter()
triangulate.SetInputData(closed_mesh)
triangulate.Update()
triangulated = triangulate.GetOutput()
# Filter out lines and vertices, keep only strips and polys
# vtkTriangleFilter outputs lines if input has lines
# We need a clean polydata with only triangles
clean_polys = vtk.vtkPolyData()
clean_polys.SetPoints(triangulated.GetPoints())
clean_polys.SetPolys(triangulated.GetPolys())
# Do NOT set lines or verts
# Calculate volume using mass properties
mass_props = vtk.vtkMassProperties()
mass_props.SetInputData(clean_polys)
mass_props.Update()
volume = abs(mass_props.GetVolume())
# Calculate wetted surface area (excluding the cap)
# Use the clipped mesh before adding the cap
wetted_area = _calculate_surface_area(clipped)
# Calculate center of buoyancy using volumetric moment integration
# vtkCenterOfMass only calculates the average of points, which is wrong for non-uniform meshes
center = _calculate_volumetric_centroid(clean_polys)
# Bounds of submerged part (clipped)
bounds = clipped.GetBounds()
return volume, center, wetted_area, bounds
def _calculate_volumetric_centroid(polydata: vtk.vtkPolyData) -> np.ndarray:
"""
Calculate the centroid of a closed polyhedron using tetrahedral decomposition.
For a closed mesh, the volume centroid is:
C = (1 / 4V) * sum(V_i * (P1_i + P2_i + P3_i + Origin))
where V_i is the signed volume of the tetrahedron formed by triangle i and the Origin.
Args:
polydata: Closed triangulated mesh.
Returns:
Numpy array [x, y, z] of the centroid.
"""
points = polydata.GetPoints()
polys = polydata.GetPolys()
# Iterate through all triangles
polys.InitTraversal()
pt_ids = vtk.vtkIdList()
total_volume = 0.0
moment = np.zeros(3)
# Pre-allocate for speed
p0 = np.zeros(3)
p1 = np.zeros(3)
p2 = np.zeros(3)
while polys.GetNextCell(pt_ids):
if pt_ids.GetNumberOfIds() != 3:
continue
points.GetPoint(pt_ids.GetId(0), p0)
points.GetPoint(pt_ids.GetId(1), p1)
points.GetPoint(pt_ids.GetId(2), p2)
# Tetrahedron volume (signed)
# V = (1/6) * det(p0, p1, p2)
# = (1/6) * p0 . (p1 x p2)
cross_prod = np.cross(p1, p2)
tet_vol = np.dot(p0, cross_prod) / 6.0
# Tetrahedron centroid (assuming Origin is the 4th point)
# C = (p0 + p1 + p2 + 0) / 4
tet_center = (p0 + p1 + p2) / 4.0
total_volume += tet_vol
moment += tet_vol * tet_center
if abs(total_volume) < 1e-9:
return np.zeros(3)
return moment / total_volume
def _calculate_surface_area(polydata: vtk.vtkPolyData) -> float:
"""
Calculate the surface area of a polydata mesh.
Args:
polydata: The mesh to measure.
Returns:
Surface area in m².
"""
# Triangulate first to ensure all cells are triangles
triangulate = vtk.vtkTriangleFilter()
triangulate.SetInputData(polydata)
triangulate.Update()
triangulated = triangulate.GetOutput()
# Calculate surface area
mass_props = vtk.vtkMassProperties()
mass_props.SetInputData(triangulated)
mass_props.Update()
return mass_props.GetSurfaceArea()
[docs]
def calculate_shell_volume(wetted_area: float, thickness: float) -> float:
"""
Calculate the volume of the shell based on wetted surface area and thickness.
This is an approximation assuming the shell follows the outer surface.
For thin shells, Volume ≈ Area × Thickness
Args:
wetted_area: Wetted surface area in m².
thickness: Shell thickness in meters.
Returns:
Shell volume in m³.
"""
return wetted_area * thickness
[docs]
def calculate_waterplane_geometric_properties(waterline_polydata: vtk.vtkPolyData) -> Tuple[float, np.ndarray, float, float]:
"""
Calculate geometric properties of the waterplane (Area, Centroid, I_t, I_l).
Uses Green's theorem for polygon properties.
Assumes the waterplane is in the XY plane (or parallel to it).
Args:
waterline_polydata: vtkPolyData containing the waterline contour(s).
Returns:
Tuple (Area, Centroid, I_t, I_l, lwl, bwl)
- Area: Waterplane area in m²
- Centroid: Cylinder center [x, y, z] (z is from input data)
- I_t: Transverse moment of inertia (about longitudinal axis through centroid) in m^4
- I_l: Longitudinal moment of inertia (about transverse axis through centroid) in m^4
- lwl: Length on waterline in m
- bwl: Beam on waterline in m
"""
# Check if waterline has any points
if waterline_polydata is None or waterline_polydata.GetNumberOfPoints() == 0:
return 0.0, np.array([0.0, 0.0, 0.0]), 0.0, 0.0, 0.0, 0.0
# Calculate geometric dimensions (Lwl, Bwl) from bounds
bounds = waterline_polydata.GetBounds()
lwl = bounds[1] - bounds[0]
bwl = bounds[3] - bounds[2]
# Use vtkStripper to ensure we have ordered lines
stripper = vtk.vtkStripper()
stripper.SetInputData(waterline_polydata)
stripper.Update()
stripped = stripper.GetOutput()
lines = stripped.GetLines()
points = stripped.GetPoints()
# Check if stripper produced valid output
if points is None or points.GetNumberOfPoints() == 0:
return 0.0, np.array([0.0, 0.0, 0.0]), 0.0, 0.0, lwl, bwl
if lines is None or lines.GetNumberOfCells() == 0:
return 0.0, np.array([0.0, 0.0, 0.0]), 0.0, 0.0, lwl, bwl
lines.InitTraversal()
pt_ids = vtk.vtkIdList()
# Accumulators
area_accum = 0.0
Cx_accum = 0.0
Cy_accum = 0.0
Ixx_accum = 0.0 # About global X axis
Iyy_accum = 0.0 # About global Y axis
z_coord = points.GetPoint(0)[2]
# Iterate through all contours
while lines.GetNextCell(pt_ids):
n = pt_ids.GetNumberOfIds()
if n < 3:
continue
# Iterate through points in the contour
for i in range(n):
# Current and next point indices (wrapping around)
id1 = pt_ids.GetId(i)
id2 = pt_ids.GetId((i + 1) % n)
p1 = points.GetPoint(id1)
p2 = points.GetPoint(id2)
x1, y1 = p1[0], p1[1]
x2, y2 = p2[0], p2[1]
# Common term: cross product z-component (2 * triangle area)
cross_z = x1 * y2 - x2 * y1
# Area
area_accum += cross_z
# Centroid moments (A * x, A * y)
Cx_accum += (x1 + x2) * cross_z
Cy_accum += (y1 + y2) * cross_z
# Second moments of area about global axes
Ixx_accum += (y1**2 + y1 * y2 + y2**2) * cross_z
Iyy_accum += (x1**2 + x1 * x2 + x2**2) * cross_z
area = area_accum / 2.0
if abs(area) < 1e-9:
return 0.0, np.array([0.0, 0.0, z_coord]), 0.0, 0.0, lwl, bwl
# Final Centroid coordinates
# Formula: Cx = (1/6A) * sum((x1+x2)*(x1y2-x2y1))
Cx = Cx_accum / 6.0
Cy = Cy_accum / 6.0
# Actual centroid
centroid_x = Cx / area
centroid_y = Cy / area
centroid = np.array([centroid_x, centroid_y, z_coord])
# Final Moments of Inertia about Global Axes
I_xx_global = Ixx_accum / 12.0
I_yy_global = Iyy_accum / 12.0
# Steiner's Theorem to get I about axes through centroid
I_t = abs(I_xx_global - area * (centroid_y**2))
I_l = abs(I_yy_global - area * (centroid_x**2))
return abs(area), centroid, I_t, I_l, lwl, bwl
[docs]
def calculate_section_area(hull_polydata: vtk.vtkPolyData, x_location: float, draft: float) -> float:
"""
Calculate the submerged area of a transverse section at x_location.
Args:
hull_polydata: Hull geometry.
x_location: Longitudinal position of the section.
draft: Waterline Z-coordinate (area is calculated below this).
Returns:
Area in m².
"""
# 1. Cut the hull at x_location
plane = vtk.vtkPlane()
plane.SetOrigin(x_location, 0, 0)
plane.SetNormal(1, 0, 0)
cutter = vtk.vtkCutter()
cutter.SetInputData(hull_polydata)
cutter.SetCutFunction(plane)
cutter.Update()
section_poly = cutter.GetOutput()
if section_poly.GetNumberOfPoints() == 0:
return 0.0
# 2. Clip the section at the waterline (keep part below draft)
wl_plane = vtk.vtkPlane()
wl_plane.SetOrigin(0, 0, draft)
wl_plane.SetNormal(0, 0, 1) # Normal Up
clipper = vtk.vtkClipPolyData()
clipper.SetInputData(section_poly)
clipper.SetClipFunction(wl_plane)
clipper.InsideOutOn() # Implicit function < 0 (Below plane)
clipper.Update()
submerged_section = clipper.GetOutput()
# 3. Calculate Area Robustly (Centerline Integration)
# For a monohull centered at Y=0, Section Area is the sum of (width_from_CL * height) slices.
# We iterate all line segments.
# Contribution = abs(Y_avg) * abs(dZ)
# This works because:
# - Port side (Y>0) traces the area between shell and CL.
# - Stbd side (Y<0) traces the area between shell and CL.
# - Sum gives total area.
# - Horizontal segments (Bottom/Deck) have dZ=0, contrib 0 (Correct).
# - Assumes no internal voids or offset hulls (catamarans).
total_area = 0.0
lines = submerged_section.GetLines()
points = submerged_section.GetPoints()
lines.InitTraversal()
pt_ids = vtk.vtkIdList()
while lines.GetNextCell(pt_ids):
n = pt_ids.GetNumberOfIds()
if n < 2:
continue
for i in range(n-1):
id1 = pt_ids.GetId(i)
id2 = pt_ids.GetId(i+1)
p1 = points.GetPoint(id1)
p2 = points.GetPoint(id2)
y1, z1 = p1[1], p1[2]
y2, z2 = p2[1], p2[2]
y_avg = (y1 + y2) / 2.0
dz = z2 - z1
# Robust integration (Monohull assumption)
total_area += abs(y_avg) * abs(dz)
return total_area