use nalgebra::Point3;
use parry3d_f64::shape::{Shape, TriMesh};
use std::path::Path;
use crate::hull::Hull;
use crate::mesh::{
calculate_waterplane_properties, clip_at_waterline, clip_by_axis_aligned_plane, load_stl,
load_vtk, Axis,
};
#[derive(Clone)]
pub struct Tank {
name: String,
mesh: TriMesh,
total_volume: f64,
fluid_density: f64,
fill_level: f64,
water_density: f64,
bounds: (f64, f64, f64, f64, f64, f64),
fsm_mode: FSMMode,
max_fsm_cache: Option<(f64, f64)>,
permeability: f64,
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum FSMMode {
Actual,
Maximum,
Fixed { t: f64, l: f64 },
}
impl Tank {
pub fn new(name: &str, mesh: TriMesh, fluid_density: f64) -> Self {
let mass_props = mesh.mass_properties(1.0);
let total_volume = mass_props.mass().abs();
let aabb = mesh.local_aabb();
let bounds = (
aabb.mins.x,
aabb.maxs.x,
aabb.mins.y,
aabb.maxs.y,
aabb.mins.z,
aabb.maxs.z,
);
Self {
name: name.to_string(),
mesh,
total_volume,
fluid_density,
fill_level: 0.0,
water_density: 1025.0,
bounds,
fsm_mode: FSMMode::Actual,
max_fsm_cache: None,
permeability: 1.0,
}
}
pub fn from_file<P: AsRef<Path>>(path: P, fluid_density: f64) -> Result<Self, String> {
let path = path.as_ref();
if !path.exists() {
return Err(format!("File not found: {}", path.display()));
}
let ext = path
.extension()
.and_then(|e| e.to_str())
.map(|e| e.to_lowercase())
.unwrap_or_default();
let mesh = match ext.as_str() {
"stl" => load_stl(path).map_err(|e| format!("Failed to load STL: {}", e))?,
"vtk" | "vtp" | "vtu" => {
load_vtk(path).map_err(|e| format!("Failed to load VTK: {}", e))?
}
_ => return Err(format!("Unsupported file format: {}", ext)),
};
if mesh.vertices().is_empty() {
return Err("Loaded mesh has no vertices".to_string());
}
let name = path.file_stem().and_then(|s| s.to_str()).unwrap_or("Tank");
Ok(Self::new(name, mesh, fluid_density))
}
#[allow(clippy::too_many_arguments)]
pub fn from_box_hull_intersection(
name: &str,
hull: &Hull,
x_min: f64,
x_max: f64,
y_min: f64,
y_max: f64,
z_min: f64,
z_max: f64,
fluid_density: f64,
) -> Result<Self, String> {
let mut mesh = hull.mesh().clone();
if x_min >= x_max || y_min >= y_max || z_min >= z_max {
return Err("Invalid box dimensions: min must be less than max".to_string());
}
mesh = clip_by_axis_aligned_plane(&mesh, Axis::X, x_min, false).0
.ok_or_else(|| format!("Tank '{}': No geometry remaining after X_min ({:.2}) clip. Box is fully aft of hull?", name, x_min))?;
mesh = clip_by_axis_aligned_plane(&mesh, Axis::X, x_max, true).0
.ok_or_else(|| format!("Tank '{}': No geometry remaining after X_max ({:.2}) clip. Box is fully fwd of hull?", name, x_max))?;
mesh = clip_by_axis_aligned_plane(&mesh, Axis::Y, y_min, false).0
.ok_or_else(|| format!("Tank '{}': No geometry remaining after Y_min ({:.2}) clip. Box is fully stbd of hull?", name, y_min))?;
mesh = clip_by_axis_aligned_plane(&mesh, Axis::Y, y_max, true).0
.ok_or_else(|| format!("Tank '{}': No geometry remaining after Y_max ({:.2}) clip. Box is fully port of hull?", name, y_max))?;
mesh = clip_by_axis_aligned_plane(&mesh, Axis::Z, z_min, false).0
.ok_or_else(|| format!("Tank '{}': No geometry remaining after Z_min ({:.2}) clip. Box is fully below hull?", name, z_min))?;
mesh = clip_by_axis_aligned_plane(&mesh, Axis::Z, z_max, true).0
.ok_or_else(|| format!("Tank '{}': No geometry remaining after Z_max ({:.2}) clip. Box is fully above hull?", name, z_max))?;
let mass_props = mesh.mass_properties(1.0);
let volume = mass_props.mass().abs();
if volume < 1e-6 {
return Err(format!("Tank '{}': Intersection resulted in near-zero volume ({:.2e} m³). Check that box overlaps hull interior.", name, volume));
}
Ok(Self::new(name, mesh, fluid_density))
}
#[allow(clippy::too_many_arguments)]
pub fn from_box(
name: &str,
x_min: f64,
x_max: f64,
y_min: f64,
y_max: f64,
z_min: f64,
z_max: f64,
fluid_density: f64,
) -> Self {
let vertices = vec![
Point3::new(x_min, y_min, z_min),
Point3::new(x_max, y_min, z_min),
Point3::new(x_max, y_max, z_min),
Point3::new(x_min, y_max, z_min),
Point3::new(x_min, y_min, z_max),
Point3::new(x_max, y_min, z_max),
Point3::new(x_max, y_max, z_max),
Point3::new(x_min, y_max, z_max),
];
let indices = vec![
[0, 2, 1],
[0, 3, 2],
[4, 5, 6],
[4, 6, 7],
[0, 1, 5],
[0, 5, 4],
[2, 3, 7],
[2, 7, 6],
[0, 4, 7],
[0, 7, 3],
[1, 2, 6],
[1, 6, 5],
];
let mesh = TriMesh::new(vertices, indices).expect("Failed to create tank mesh");
Self::new(name, mesh, fluid_density)
}
pub fn name(&self) -> &str {
&self.name
}
pub fn get_bounds(&self) -> (f64, f64, f64, f64, f64, f64) {
self.bounds
}
pub fn set_name(&mut self, name: &str) {
self.name = name.to_string();
}
pub fn total_volume(&self) -> f64 {
self.total_volume
}
pub fn fill_level(&self) -> f64 {
self.fill_level
}
pub fn set_fill_level(&mut self, level: f64) {
self.fill_level = level.clamp(0.0, 1.0);
}
pub fn permeability(&self) -> f64 {
self.permeability
}
pub fn set_permeability(&mut self, permeability: f64) {
self.permeability = permeability.clamp(0.0, 1.0);
if self.fsm_mode == FSMMode::Maximum {
self.max_fsm_cache = None;
self.max_fsm_cache = Some(self.calculate_max_fsm());
}
}
pub fn fill_percent(&self) -> f64 {
self.fill_level * 100.0
}
pub fn set_fill_percent(&mut self, percent: f64) {
self.fill_level = (percent / 100.0).clamp(0.0, 1.0);
}
pub fn fluid_density(&self) -> f64 {
self.fluid_density
}
pub fn set_fluid_density(&mut self, density: f64) {
self.fluid_density = density;
}
pub fn fill_volume(&self) -> f64 {
self.total_volume * self.fill_level * self.permeability
}
pub fn fluid_mass(&self) -> f64 {
self.fill_volume() * self.fluid_density
}
pub fn get_fluid_mesh(&self) -> Option<TriMesh> {
self.get_fluid_mesh_at(0.0, 0.0)
}
pub fn get_fluid_mesh_at(&self, heel: f64, trim: f64) -> Option<TriMesh> {
if self.fill_level <= 1e-6 {
return None;
}
if self.fill_level >= 1.0 - 1e-6 {
return Some(self.mesh.clone());
}
let pivot = nalgebra::Point3::new(0.0, 0.0, 0.0);
let transformed_mesh = crate::mesh::transform_mesh(&self.mesh, heel, trim, pivot);
let z = self.find_z_for_mesh(&transformed_mesh, self.total_volume * self.fill_level);
if let Some(clipped) = clip_at_waterline(&transformed_mesh, z).0 {
use nalgebra::Rotation3;
let roll = heel.to_radians();
let pitch = trim.to_radians();
let rotation = Rotation3::from_euler_angles(roll, pitch, 0.0);
let pivot_pt = nalgebra::Point3::new(0.0, 0.0, 0.0);
let inverse_rotation = rotation.inverse();
let new_vertices: Vec<Point3<f64>> = clipped
.vertices()
.iter()
.map(|v| {
let p = Point3::new(v.x, v.y, v.z);
pivot_pt + inverse_rotation * (p - pivot_pt)
})
.collect();
Some(TriMesh::new(new_vertices, clipped.indices().to_vec()).unwrap())
} else {
None
}
}
pub fn mesh(&self) -> &TriMesh {
&self.mesh
}
pub fn center_of_gravity(&self) -> [f64; 3] {
self.center_of_gravity_at(0.0, 0.0)
}
pub fn center_of_gravity_at(&self, heel: f64, trim: f64) -> [f64; 3] {
if self.fill_level <= 0.0 {
return [0.0, 0.0, 0.0];
}
let pivot = nalgebra::Point3::new(0.0, 0.0, 0.0);
let transformed_mesh = crate::mesh::transform_mesh(&self.mesh, heel, trim, pivot);
let z = self.find_z_for_mesh(&transformed_mesh, self.total_volume * self.fill_level);
let transformed_cog = if let Some(clipped) = clip_at_waterline(&transformed_mesh, z).0 {
let mass_props = clipped.mass_properties(1.0);
let com = mass_props.local_com;
nalgebra::Point3::new(com.x, com.y, com.z)
} else {
let mass_props = transformed_mesh.mass_properties(1.0);
let com = mass_props.local_com;
nalgebra::Point3::new(com.x, com.y, com.z)
};
use nalgebra::Rotation3;
let roll = heel.to_radians();
let pitch = trim.to_radians();
let rotation = Rotation3::from_euler_angles(roll, pitch, 0.0);
let inverse_rotation = rotation.inverse();
let original_cog = inverse_rotation * (transformed_cog - pivot) + pivot.coords;
[original_cog.x, original_cog.y, original_cog.z]
}
fn find_z_for_mesh(&self, mesh: &TriMesh, target_volume: f64) -> f64 {
let aabb = mesh.local_aabb();
let z_min = aabb.mins.z;
let z_max = aabb.maxs.z;
let tolerance = target_volume.max(0.001) * 1e-4;
let max_iter = 20;
let mut low = z_min;
let mut high = z_max;
for _ in 0..max_iter {
let mid = (low + high) / 2.0;
let volume = if let Some(clipped) = clip_at_waterline(mesh, mid).0 {
clipped.mass_properties(1.0).mass().abs()
} else {
0.0
};
let diff = volume - target_volume;
if diff.abs() < tolerance {
return mid;
}
if diff > 0.0 {
high = mid;
} else {
low = mid;
}
}
(low + high) / 2.0
}
pub fn free_surface_moment_t(&self) -> f64 {
match self.fsm_mode {
FSMMode::Actual => {
if self.fill_level <= 0.0 || self.fill_level >= 0.98 {
if self.fill_level >= 1.0 - 1e-6 {
return 0.0;
}
if self.fill_level <= 1e-6 {
return 0.0;
}
}
let target_volume = self.total_volume * self.fill_level;
let z = self.find_z_for_mesh(&self.mesh, target_volume);
if let Some(wp) = calculate_waterplane_properties(&self.mesh, z) {
wp.i_transverse * self.permeability
} else {
0.0
}
}
FSMMode::Maximum => self.max_fsm_cache.map(|(t, _)| t).unwrap_or(0.0),
FSMMode::Fixed { t, .. } => t,
}
}
pub fn free_surface_moment_l(&self) -> f64 {
match self.fsm_mode {
FSMMode::Actual => {
if self.fill_level <= 1e-6 || self.fill_level >= 1.0 - 1e-6 {
return 0.0;
}
let target_volume = self.total_volume * self.fill_level;
let z = self.find_z_for_mesh(&self.mesh, target_volume);
if let Some(wp) = calculate_waterplane_properties(&self.mesh, z) {
wp.i_longitudinal * self.permeability
} else {
0.0
}
}
FSMMode::Maximum => self.max_fsm_cache.map(|(_, l)| l).unwrap_or(0.0),
FSMMode::Fixed { l, .. } => l,
}
}
pub fn set_fsm_mode(&mut self, mode: FSMMode) {
self.fsm_mode = mode;
if self.fsm_mode == FSMMode::Maximum && self.max_fsm_cache.is_none() {
self.max_fsm_cache = Some(self.calculate_max_fsm());
}
}
fn calculate_max_fsm(&self) -> (f64, f64) {
let levels = [
0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99,
];
let mut max_t = 0.0;
let mut max_l = 0.0;
for level in levels {
let target_volume = self.total_volume * level;
let z = self.find_z_for_mesh(&self.mesh, target_volume);
if let Some(wp) = calculate_waterplane_properties(&self.mesh, z) {
let wp_it = wp.i_transverse * self.permeability;
if wp_it > max_t {
max_t = wp_it;
}
let wp_il = wp.i_longitudinal * self.permeability;
if wp_il > max_l {
max_l = wp_il;
}
}
}
(max_t, max_l)
}
pub fn fsm_mode(&self) -> FSMMode {
self.fsm_mode
}
pub fn free_surface_correction_t(&self) -> f64 {
self.free_surface_moment_t() * (self.fluid_density / self.water_density)
}
pub fn free_surface_correction_l(&self) -> f64 {
self.free_surface_moment_l() * (self.fluid_density / self.water_density)
}
}
impl std::fmt::Debug for Tank {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Tank")
.field("name", &self.name)
.field("total_volume", &self.total_volume)
.field("fill_level", &self.fill_level)
.field("fluid_density", &self.fluid_density)
.finish()
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_box_tank_volume() {
let tank = Tank::from_box("Test", 0.0, 10.0, 0.0, 5.0, 0.0, 2.0, 1000.0);
assert!(
(tank.total_volume() - 100.0).abs() < 0.1,
"Volume was {}",
tank.total_volume()
);
}
#[test]
fn test_tank_fill() {
let mut tank = Tank::from_box("Test", 0.0, 10.0, 0.0, 5.0, 0.0, 2.0, 1000.0);
tank.set_fill_percent(50.0);
assert!((tank.fill_level() - 0.5).abs() < 1e-6);
assert!((tank.fill_volume() - 50.0).abs() < 0.1);
assert!((tank.fluid_mass() - 50000.0).abs() < 100.0);
}
#[test]
fn test_from_box_hull_intersection() {
use crate::hull::Hull;
let hull = Hull::from_box(20.0, 10.0, 10.0);
let tank = Tank::from_box_hull_intersection(
"InnerTank",
&hull,
5.0,
15.0,
-2.0,
2.0,
0.0,
5.0,
1025.0,
)
.expect("Intersection failed");
assert!(
(tank.total_volume() - 200.0).abs() < 1e-3,
"Volume was {}",
tank.total_volume()
);
let res = Tank::from_box_hull_intersection(
"Outside", &hull, 30.0, 40.0, 0.0, 1.0, 0.0, 1.0, 1025.0,
);
assert!(res.is_err());
}
#[test]
fn test_tank_fsm() {
let tank = Tank::from_box("FSM_Test", 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 1000.0);
let mut tank = tank;
tank.set_fill_percent(50.0);
let i_t = tank.free_surface_moment_t();
let i_l = tank.free_surface_moment_l();
assert!((i_t - 833.333).abs() < 1.0, "I_t was {}", i_t);
assert!((i_l - 833.333).abs() < 1.0, "I_l was {}", i_l);
tank.set_fill_level(1.0);
assert_eq!(tank.free_surface_moment_t(), 0.0);
tank.set_fill_level(0.0);
assert_eq!(tank.free_surface_moment_t(), 0.0);
}
#[test]
fn test_tank_permeability() {
let mut tank = Tank::from_box("Perm_Test", 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 1000.0);
tank.set_fill_percent(50.0);
let orig_volume = tank.fill_volume();
let orig_mass = tank.fluid_mass();
let orig_fsm_t = tank.free_surface_moment_t();
tank.set_permeability(0.9);
assert!((tank.fill_volume() - orig_volume * 0.9).abs() < 1e-6);
assert!((tank.fluid_mass() - orig_mass * 0.9).abs() < 1e-6);
assert!((tank.free_surface_moment_t() - orig_fsm_t * 0.9).abs() < 1e-6);
}
#[test]
fn test_cog_shift_direction() {
let tank = Tank::from_box("ShiftTest", 0.0, 2.0, -5.0, 5.0, 0.0, 2.0, 1000.0);
let mut tank = tank;
tank.set_fill_percent(50.0);
let cog_upright = tank.center_of_gravity();
assert!((cog_upright[1] - 0.0).abs() < 1e-6);
let cog_heeled = tank.center_of_gravity_at(10.0, 0.0);
assert!(
cog_heeled[1] < -0.01,
"COG Y should be negative (Starboard) for positive heel. Got: {}",
cog_heeled[1]
);
}
}