navaltoolbox 0.9.2

// Copyright (C) 2026 Antoine ANCEAU
//
// This file is part of navaltoolbox.
//
// navaltoolbox 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/>.

//! Vessel module.
//!
//! Provides the Vessel container for hull geometries, tanks, and vessel-level properties.

use crate::appendage::Appendage;
use crate::deckedge::DeckEdge;
use crate::downflooding::DownfloodingOpening;
use crate::hull::Hull;
use crate::mesh::get_bounds;
use crate::silhouette::Silhouette;
use crate::tanks::SharedTank;

/// Pre-computed contact surface between two hulls.
///
/// Stores the face indices from each hull that are in contact,
/// along with the total contact area. This avoids expensive O(N×M)
/// runtime detection in hydrostatic and stability calculations.
#[derive(Clone, Debug)]
pub struct ContactSurface {
    /// Index of the first hull
    pub hull_i: usize,
    /// Index of the second hull
    pub hull_j: usize,
    /// Indices of faces in hull_i that are in contact
    pub face_indices_i: Vec<usize>,
    /// Indices of faces in hull_j that are in contact
    pub face_indices_j: Vec<usize>,
    /// Total pre-computed contact area in m²
    pub total_area: f64,
}

/// Represents a vessel containing hull geometries, tanks, and vessel-level properties.
///
/// The Vessel class serves as a container for hull geometries and manages
/// vessel-level reference positions such as the forward and aft perpendiculars
/// (FP and AP). It supports both single-hull vessels (monohulls) and multi-hull
/// vessels (catamarans, trimarans).
#[derive(Clone)]
pub struct Vessel {
    /// List of hull geometries
    hulls: Vec<Hull>,
    /// List of tanks (shared references)
    tanks: Vec<SharedTank>,
    /// List of appendages
    appendages: Vec<Appendage>,
    /// Deck edges for freeboard calculation
    deck_edges: Vec<DeckEdge>,
    /// Aft Perpendicular position (None = auto from bounds)
    ap: Option<f64>,
    /// Forward Perpendicular position (None = auto from bounds)
    fp: Option<f64>,
    /// Wind silhouette profiles (hull, superstructure, containers, etc.)
    silhouettes: Vec<Silhouette>,
    /// Downflooding openings for θf calculation
    downflooding_openings: Vec<DownfloodingOpening>,
    /// Pre-computed contact surfaces between hull pairs
    contact_surfaces: Vec<ContactSurface>,
}

impl Vessel {
    /// Creates a new Vessel with a single hull.
    pub fn new(hull: Hull) -> Self {
        Self {
            hulls: vec![hull],
            tanks: Vec::new(),
            appendages: Vec::new(),
            deck_edges: Vec::new(),
            ap: None,
            fp: None,
            silhouettes: Vec::new(),
            downflooding_openings: Vec::new(),
            contact_surfaces: Vec::new(),
        }
    }

    /// Creates a new Vessel with multiple hulls (catamaran, trimaran).
    pub fn new_multi(hulls: Vec<Hull>) -> Result<Self, &'static str> {
        if hulls.is_empty() {
            return Err("At least one hull must be provided");
        }
        let mut vessel = Self {
            hulls,
            tanks: Vec::new(),
            appendages: Vec::new(),
            deck_edges: Vec::new(),
            ap: None,
            fp: None,
            silhouettes: Vec::new(),
            downflooding_openings: Vec::new(),
            contact_surfaces: Vec::new(),
        };
        // Auto-compute contact surfaces for multi-hull vessels
        vessel.compute_contact_surfaces();
        Ok(vessel)
    }

    /// Creates a new Vessel with perpendicular positions.
    pub fn with_perpendiculars(hull: Hull, ap: f64, fp: f64) -> Self {
        Self {
            hulls: vec![hull],
            tanks: Vec::new(),
            appendages: Vec::new(),
            deck_edges: Vec::new(),
            ap: Some(ap),
            fp: Some(fp),
            silhouettes: Vec::new(),
            downflooding_openings: Vec::new(),
            contact_surfaces: Vec::new(),
        }
    }

    /// Returns the list of hull geometries.
    pub fn hulls(&self) -> &[Hull] {
        &self.hulls
    }

    /// Returns a mutable reference to the hulls.
    pub fn hulls_mut(&mut self) -> &mut Vec<Hull> {
        &mut self.hulls
    }

    /// Returns true if this is a multi-hull vessel.
    pub fn is_multihull(&self) -> bool {
        self.hulls.len() > 1
    }

    /// Sets the hull plate thickness for a specific hull by index.
    pub fn set_hull_thickness(
        &mut self,
        index: usize,
        thickness: Option<f64>,
    ) -> Result<(), &'static str> {
        if index < self.hulls.len() {
            self.hulls[index].set_thickness(thickness);
            Ok(())
        } else {
            Err("Hull index out of bounds")
        }
    }

    /// Returns the hull plate thickness for a specific hull by index.
    pub fn get_hull_thickness(&self, index: usize) -> Option<f64> {
        if index < self.hulls.len() {
            self.hulls[index].thickness()
        } else {
            None
        }
    }

    /// Returns the list of tanks.
    pub fn tanks(&self) -> &[SharedTank] {
        &self.tanks
    }

    /// Returns a mutable reference to the list of tanks.
    pub fn tanks_mut(&mut self) -> &mut Vec<SharedTank> {
        &mut self.tanks
    }

    /// Returns the Aft Perpendicular position.
    ///
    /// If not explicitly set, returns the minimum X of the combined bounds.
    pub fn ap(&self) -> f64 {
        self.ap.unwrap_or_else(|| self.get_bounds().0)
    }

    /// Returns the Forward Perpendicular position.
    ///
    /// If not explicitly set, returns the maximum X of the combined bounds.
    pub fn fp(&self) -> f64 {
        self.fp.unwrap_or_else(|| self.get_bounds().1)
    }

    /// Sets the Aft Perpendicular position.
    pub fn set_ap(&mut self, ap: f64) {
        self.ap = Some(ap);
    }

    /// Sets the Forward Perpendicular position.
    pub fn set_fp(&mut self, fp: f64) {
        self.fp = Some(fp);
    }

    /// Returns the Length Between Perpendiculars (LBP).
    pub fn lbp(&self) -> f64 {
        self.fp() - self.ap()
    }

    /// Returns the bounding box of all hull geometries combined.
    ///
    /// Returns (xmin, xmax, ymin, ymax, zmin, zmax).
    pub fn get_bounds(&self) -> (f64, f64, f64, f64, f64, f64) {
        if self.hulls.len() == 1 {
            return self.hulls[0].get_bounds();
        }

        let all_bounds: Vec<_> = self.hulls.iter().map(|h| h.get_bounds()).collect();

        let xmin = all_bounds.iter().map(|b| b.0).fold(f64::INFINITY, f64::min);
        let xmax = all_bounds
            .iter()
            .map(|b| b.1)
            .fold(f64::NEG_INFINITY, f64::max);
        let ymin = all_bounds.iter().map(|b| b.2).fold(f64::INFINITY, f64::min);
        let ymax = all_bounds
            .iter()
            .map(|b| b.3)
            .fold(f64::NEG_INFINITY, f64::max);
        let zmin = all_bounds.iter().map(|b| b.4).fold(f64::INFINITY, f64::min);
        let zmax = all_bounds
            .iter()
            .map(|b| b.5)
            .fold(f64::NEG_INFINITY, f64::max);

        (xmin, xmax, ymin, ymax, zmin, zmax)
    }

    // =========================================================================
    // Tank Management
    // =========================================================================

    /// Adds a tank to the vessel (takes shared ownership).
    pub fn add_tank(&mut self, tank: SharedTank) {
        self.tanks.push(tank);
    }

    /// Removes a tank from the vessel by index.
    pub fn remove_tank(&mut self, index: usize) -> Option<SharedTank> {
        if index < self.tanks.len() {
            Some(self.tanks.remove(index))
        } else {
            None
        }
    }

    /// Finds a tank by its name. Returns a shared reference.
    pub fn get_tank_by_name(&self, name: &str) -> Option<SharedTank> {
        self.tanks
            .iter()
            .find(|t| t.read().unwrap().name() == name)
            .cloned()
    }

    /// Finds a tank by its name (mutable access via write lock on returned Arc).
    pub fn get_tank_by_name_mut(&self, name: &str) -> Option<SharedTank> {
        self.get_tank_by_name(name)
    }

    /// Calculates the total mass of all fluid in tanks.
    pub fn get_total_tanks_mass(&self) -> f64 {
        self.tanks
            .iter()
            .map(|t| t.read().unwrap().fluid_mass())
            .sum()
    }

    /// Calculates the combined center of gravity of all tank fluids.
    ///
    /// Returns the mass-weighted average of individual tank CoGs.
    pub fn get_tanks_center_of_gravity(&self) -> [f64; 3] {
        let total_mass = self.get_total_tanks_mass();
        if total_mass <= 0.0 {
            return [0.0, 0.0, 0.0];
        }

        let mut moment = [0.0, 0.0, 0.0];
        for tank_arc in &self.tanks {
            let tank = tank_arc.read().unwrap();
            if tank.fluid_mass() > 0.0 {
                let cog = tank.center_of_gravity();
                moment[0] += tank.fluid_mass() * cog[0];
                moment[1] += tank.fluid_mass() * cog[1];
                moment[2] += tank.fluid_mass() * cog[2];
            }
        }

        [
            moment[0] / total_mass,
            moment[1] / total_mass,
            moment[2] / total_mass,
        ]
    }

    /// Calculates the total free surface moment from all tanks.
    ///
    /// Returns (transverse_moment, longitudinal_moment) in m⁴.
    pub fn get_total_free_surface_moment(&self) -> (f64, f64) {
        let fsm_t: f64 = self
            .tanks
            .iter()
            .map(|t| t.read().unwrap().free_surface_moment_t())
            .sum();
        let fsm_l: f64 = self
            .tanks
            .iter()
            .map(|t| t.read().unwrap().free_surface_moment_l())
            .sum();
        (fsm_t, fsm_l)
    }

    /// Calculates the total free surface correction from all tanks.
    ///
    /// Returns (transverse_correction, longitudinal_correction) in m⁴.
    pub fn get_total_free_surface_correction(&self) -> (f64, f64) {
        let fsc_t: f64 = self
            .tanks
            .iter()
            .map(|t| t.read().unwrap().free_surface_correction_t())
            .sum();
        let fsc_l: f64 = self
            .tanks
            .iter()
            .map(|t| t.read().unwrap().free_surface_correction_l())
            .sum();
        (fsc_t, fsc_l)
    }

    // =========================================================================
    // Silhouette Management
    // =========================================================================

    /// Adds a wind silhouette profile to the vessel.
    pub fn add_silhouette(&mut self, silhouette: Silhouette) {
        self.silhouettes.push(silhouette);
    }

    /// Returns a reference to all wind silhouettes.
    pub fn silhouettes(&self) -> &[Silhouette] {
        &self.silhouettes
    }

    /// Returns a mutable reference to the silhouettes.
    pub fn silhouettes_mut(&mut self) -> &mut Vec<Silhouette> {
        &mut self.silhouettes
    }

    /// Returns the number of silhouettes.
    pub fn num_silhouettes(&self) -> usize {
        self.silhouettes.len()
    }

    /// Returns true if there are any silhouettes.
    pub fn has_silhouettes(&self) -> bool {
        !self.silhouettes.is_empty()
    }

    /// Finds a silhouette by its name.
    pub fn get_silhouette_by_name(&self, name: &str) -> Option<&Silhouette> {
        self.silhouettes.iter().find(|s| s.name() == name)
    }

    /// Removes a silhouette by index.
    pub fn remove_silhouette(&mut self, index: usize) -> Option<Silhouette> {
        if index < self.silhouettes.len() {
            Some(self.silhouettes.remove(index))
        } else {
            None
        }
    }

    /// Removes all silhouettes.
    pub fn clear_silhouettes(&mut self) {
        self.silhouettes.clear();
    }

    /// Calculates the total emerged area from all silhouettes.
    pub fn get_total_emerged_area(&self, waterline_z: f64) -> f64 {
        self.silhouettes
            .iter()
            .map(|s| s.get_emerged_area(waterline_z))
            .sum()
    }

    /// Calculates the total submerged area from all silhouettes.
    pub fn get_total_submerged_area(&self, waterline_z: f64) -> f64 {
        self.silhouettes
            .iter()
            .map(|s| s.get_submerged_area(waterline_z))
            .sum()
    }

    /// Calculates the combined centroid of all emerged areas.
    pub fn get_combined_emerged_centroid(&self, waterline_z: f64) -> [f64; 2] {
        let total_area = self.get_total_emerged_area(waterline_z);
        if total_area < 1e-9 {
            return [0.0, 0.0];
        }

        let mut cx = 0.0;
        let mut cz = 0.0;
        for s in &self.silhouettes {
            let area = s.get_emerged_area(waterline_z);
            if area > 1e-9 {
                let centroid = s.get_emerged_centroid(waterline_z);
                cx += centroid[0] * area;
                cz += centroid[1] * area;
            }
        }

        [cx / total_area, cz / total_area]
    }

    /// Calculates the combined centroid of all submerged lateral areas.
    ///
    /// Returns the area-weighted centroid [x, z] of the underwater portions
    /// of all silhouettes. Used for the exact Z lever computation per
    /// IMO 2008 IS Code §2.3.2: Z = emerged_centroid_z - submerged_centroid_z.
    pub fn get_combined_submerged_centroid(&self, waterline_z: f64) -> [f64; 2] {
        let total_area: f64 = self
            .silhouettes
            .iter()
            .map(|s| s.get_submerged_area(waterline_z))
            .sum();

        if total_area < 1e-9 {
            return [0.0, 0.0];
        }

        let mut cx = 0.0;
        let mut cz = 0.0;
        for s in &self.silhouettes {
            let area = s.get_submerged_area(waterline_z);
            if area > 1e-9 {
                let centroid = s.get_submerged_centroid(waterline_z);
                cx += centroid[0] * area;
                cz += centroid[1] * area;
            }
        }

        [cx / total_area, cz / total_area]
    }

    // =========================================================================
    // Appendage Management
    // =========================================================================

    /// Adds an appendage to the vessel.
    pub fn add_appendage(&mut self, appendage: Appendage) {
        self.appendages.push(appendage);
    }

    /// Returns a reference to all appendages.
    pub fn appendages(&self) -> &[Appendage] {
        &self.appendages
    }

    /// Returns a mutable reference to the appendages.
    pub fn appendages_mut(&mut self) -> &mut Vec<Appendage> {
        &mut self.appendages
    }

    /// Returns the number of appendages.
    pub fn num_appendages(&self) -> usize {
        self.appendages.len()
    }

    /// Removes an appendage by index.
    pub fn delete_appendage(&mut self, index: usize) -> Option<Appendage> {
        if index < self.appendages.len() {
            Some(self.appendages.remove(index))
        } else {
            None
        }
    }

    /// Removes an appendage by name. Returns the removed appendage if found.
    pub fn delete_appendage_by_name(&mut self, name: &str) -> Option<Appendage> {
        if let Some(idx) = self.appendages.iter().position(|a| a.name() == name) {
            Some(self.appendages.remove(idx))
        } else {
            None
        }
    }

    /// Finds an appendage by its name.
    pub fn get_appendage_by_name(&self, name: &str) -> Option<&Appendage> {
        self.appendages.iter().find(|a| a.name() == name)
    }

    /// Finds an appendage by its name (mutable).
    pub fn get_appendage_by_name_mut(&mut self, name: &str) -> Option<&mut Appendage> {
        self.appendages.iter_mut().find(|a| a.name() == name)
    }

    /// Calculates the total volume of all appendages.
    pub fn get_total_appendage_volume(&self) -> f64 {
        self.appendages.iter().map(|a| a.volume()).sum()
    }

    /// Calculates the total wetted surface of all appendages (if specified).
    pub fn get_total_appendage_wetted_surface(&self) -> f64 {
        self.appendages
            .iter()
            .filter_map(|a| a.wetted_surface())
            .sum()
    }

    /// Removes all appendages.
    pub fn clear_appendages(&mut self) {
        self.appendages.clear();
    }

    // =========================================================================
    // Deck Edge Management
    // =========================================================================

    /// Adds a deck edge to the vessel.
    pub fn add_deck_edge(&mut self, deck_edge: DeckEdge) {
        self.deck_edges.push(deck_edge);
    }

    /// Returns a reference to all deck edges.
    pub fn deck_edges(&self) -> &[DeckEdge] {
        &self.deck_edges
    }

    /// Returns a mutable reference to the deck edges.
    pub fn deck_edges_mut(&mut self) -> &mut Vec<DeckEdge> {
        &mut self.deck_edges
    }

    /// Returns the number of deck edges.
    pub fn num_deck_edges(&self) -> usize {
        self.deck_edges.len()
    }

    /// Returns true if any deck edges are defined.
    pub fn has_deck_edges(&self) -> bool {
        !self.deck_edges.is_empty()
    }

    /// Removes a deck edge by index.
    pub fn delete_deck_edge(&mut self, index: usize) -> Option<DeckEdge> {
        if index < self.deck_edges.len() {
            Some(self.deck_edges.remove(index))
        } else {
            None
        }
    }

    /// Finds a deck edge by its name.
    pub fn get_deck_edge_by_name(&self, name: &str) -> Option<&DeckEdge> {
        self.deck_edges.iter().find(|d| d.name() == name)
    }

    /// Removes all deck edges.
    pub fn clear_deck_edges(&mut self) {
        self.deck_edges.clear();
    }

    /// Calculates the minimum freeboard across all deck edges.
    ///
    /// # Arguments
    /// * `heel` - Heel angle in degrees
    /// * `trim` - Trim angle in degrees
    /// * `waterline_z` - Waterline Z coordinate
    ///
    /// # Returns
    /// Minimum freeboard in meters, or None if no deck edges defined
    pub fn get_min_freeboard(&self, heel: f64, trim: f64, waterline_z: f64) -> Option<f64> {
        if self.deck_edges.is_empty() {
            return None;
        }

        let bounds = self.get_bounds();
        let pivot = [
            (self.ap() + self.fp()) / 2.0,
            (bounds.2 + bounds.3) / 2.0,
            waterline_z,
        ];

        Some(
            self.deck_edges
                .iter()
                .map(|de| de.get_freeboard(heel, trim, pivot, waterline_z))
                .fold(f64::INFINITY, f64::min),
        )
    }

    // =========================================================================
    // Downflooding Openings Management
    // =========================================================================

    /// Adds a downflooding opening to the vessel.
    pub fn add_downflooding_opening(&mut self, opening: DownfloodingOpening) {
        self.downflooding_openings.push(opening);
    }

    /// Returns a reference to all downflooding openings.
    pub fn downflooding_openings(&self) -> &[DownfloodingOpening] {
        &self.downflooding_openings
    }

    /// Returns a mutable reference to downflooding openings.
    pub fn downflooding_openings_mut(&mut self) -> &mut Vec<DownfloodingOpening> {
        &mut self.downflooding_openings
    }

    /// Returns the number of downflooding openings.
    pub fn num_downflooding_openings(&self) -> usize {
        self.downflooding_openings.len()
    }

    /// Returns true if any downflooding openings are defined.
    pub fn has_downflooding_openings(&self) -> bool {
        !self.downflooding_openings.is_empty()
    }

    /// Removes all downflooding openings from the vessel.
    pub fn clear_downflooding_openings(&mut self) {
        self.downflooding_openings.clear();
    }

    // =========================================================================
    // Contact Surfaces Management
    // =========================================================================

    /// Pre-computes contact surfaces between all hull pairs.
    ///
    /// Uses an adaptive distance threshold based on the average cell size
    /// in the overlap zone between each hull pair. This makes the detection
    /// scale-independent.
    ///
    /// This is automatically called when creating a multi-hull vessel
    /// via `new_multi()`.
    pub fn compute_contact_surfaces(&mut self) {
        self.contact_surfaces.clear();
        if self.hulls.len() < 2 {
            return;
        }

        for i in 0..self.hulls.len() {
            for j in (i + 1)..self.hulls.len() {
                let mesh_a = self.hulls[i].mesh();
                let mesh_b = self.hulls[j].mesh();

                let bounds_a = get_bounds(mesh_a);
                let bounds_b = get_bounds(mesh_b);

                // Compute adaptive threshold from average edge length
                // in the bounding-box overlap zone
                let threshold = Self::adaptive_threshold(mesh_a, mesh_b, bounds_a, bounds_b);
                if threshold <= 0.0 {
                    continue;
                }

                // Detect contact faces using the adaptive threshold
                let (faces_i, faces_j, total_area) =
                    Self::detect_contact_faces(mesh_a, mesh_b, threshold);

                if !faces_i.is_empty() {
                    self.contact_surfaces.push(ContactSurface {
                        hull_i: i,
                        hull_j: j,
                        face_indices_i: faces_i,
                        face_indices_j: faces_j,
                        total_area,
                    });
                }
            }
        }
    }

    /// Returns the pre-computed contact surfaces.
    pub fn contact_surfaces(&self) -> &[ContactSurface] {
        &self.contact_surfaces
    }

    /// Returns true if contact surfaces have been pre-computed.
    pub fn has_contact_surfaces(&self) -> bool {
        !self.contact_surfaces.is_empty()
    }

    /// Clears all pre-computed contact surfaces.
    pub fn clear_contact_surfaces(&mut self) {
        self.contact_surfaces.clear();
    }

    /// Computes an adaptive distance threshold based on the average cell size
    /// in the overlap zone between two meshes.
    fn adaptive_threshold(
        mesh_a: &parry3d_f64::shape::TriMesh,
        mesh_b: &parry3d_f64::shape::TriMesh,
        bounds_a: (f64, f64, f64, f64, f64, f64),
        bounds_b: (f64, f64, f64, f64, f64, f64),
    ) -> f64 {
        // Find the overlap bounding box (with a small margin)
        let margin = 1.0; // 1m margin to catch faces near the boundary
        let ox_min = bounds_a.0.max(bounds_b.0) - margin;
        let ox_max = bounds_a.1.min(bounds_b.1) + margin;
        let oy_min = bounds_a.2.max(bounds_b.2) - margin;
        let oy_max = bounds_a.3.min(bounds_b.3) + margin;
        let oz_min = bounds_a.4.max(bounds_b.4) - margin;
        let oz_max = bounds_a.5.min(bounds_b.5) + margin;

        // No overlap
        if ox_min > ox_max || oy_min > oy_max || oz_min > oz_max {
            return 0.0;
        }

        // Compute average edge length for faces whose centroids are in the overlap zone
        let mut total_edge_len = 0.0;
        let mut count = 0usize;

        for mesh in [mesh_a, mesh_b] {
            let verts = mesh.vertices();
            for f in mesh.indices() {
                let v0 = verts[f[0] as usize];
                let v1 = verts[f[1] as usize];
                let v2 = verts[f[2] as usize];
                let cx = (v0.x + v1.x + v2.x) / 3.0;
                let cy = (v0.y + v1.y + v2.y) / 3.0;
                let cz = (v0.z + v1.z + v2.z) / 3.0;

                if cx >= ox_min
                    && cx <= ox_max
                    && cy >= oy_min
                    && cy <= oy_max
                    && cz >= oz_min
                    && cz <= oz_max
                {
                    // Average of 3 edge lengths
                    let e0 = nalgebra::distance(&v0, &v1);
                    let e1 = nalgebra::distance(&v1, &v2);
                    let e2 = nalgebra::distance(&v2, &v0);
                    total_edge_len += (e0 + e1 + e2) / 3.0;
                    count += 1;
                }
            }
        }

        if count == 0 {
            return 0.0;
        }

        // Threshold = fraction of the average edge length
        let avg_edge = total_edge_len / count as f64;
        avg_edge * 0.5
    }

    /// Detects contact faces between two meshes.
    ///
    /// Returns (face_indices_a, face_indices_b, total_contact_area).
    fn detect_contact_faces(
        mesh_a: &parry3d_f64::shape::TriMesh,
        mesh_b: &parry3d_f64::shape::TriMesh,
        distance_threshold: f64,
    ) -> (Vec<usize>, Vec<usize>, f64) {
        let verts_a = mesh_a.vertices();
        let faces_a = mesh_a.indices();
        let verts_b = mesh_b.vertices();
        let faces_b = mesh_b.indices();

        struct FaceData {
            centroid: nalgebra::Point3<f64>,
            normal: nalgebra::Vector3<f64>,
            area: f64,
        }

        let b_data: Vec<FaceData> = faces_b
            .iter()
            .map(|f| {
                let v0 = verts_b[f[0] as usize];
                let v1 = verts_b[f[1] as usize];
                let v2 = verts_b[f[2] as usize];
                let cross = (v1 - v0).cross(&(v2 - v0));
                let area = cross.norm() / 2.0;
                let normal = if area > 0.0 {
                    cross / (area * 2.0)
                } else {
                    nalgebra::Vector3::zeros()
                };
                let centroid = nalgebra::Point3::from((v0.coords + v1.coords + v2.coords) / 3.0);
                FaceData {
                    centroid,
                    normal,
                    area,
                }
            })
            .collect();

        let mut contact_faces_a = Vec::new();
        let mut contact_faces_b_set = std::collections::HashSet::new();
        let mut total_area = 0.0;

        for (idx_a, f) in faces_a.iter().enumerate() {
            let v0 = verts_a[f[0] as usize];
            let v1 = verts_a[f[1] as usize];
            let v2 = verts_a[f[2] as usize];
            let cross = (v1 - v0).cross(&(v2 - v0));
            let area = cross.norm() / 2.0;
            if area == 0.0 {
                continue;
            }
            let normal = cross / (area * 2.0);
            let centroid = nalgebra::Point3::from((v0.coords + v1.coords + v2.coords) / 3.0);

            for (idx_b, fb) in b_data.iter().enumerate() {
                if fb.area == 0.0 {
                    continue;
                }
                let plane_dist = ((centroid.coords - fb.centroid.coords).dot(&fb.normal)).abs();
                if plane_dist < distance_threshold && normal.dot(&fb.normal) < -0.9 {
                    let dist_sq = nalgebra::distance_squared(&centroid, &fb.centroid);
                    let allowed_dist = area.sqrt() + fb.area.sqrt() + distance_threshold + 1.0;
                    if dist_sq < allowed_dist * allowed_dist {
                        contact_faces_a.push(idx_a);
                        contact_faces_b_set.insert(idx_b);
                        total_area += area;
                        break;
                    }
                }
            }
        }

        let contact_faces_b: Vec<usize> = contact_faces_b_set.into_iter().collect();
        (contact_faces_a, contact_faces_b, total_area)
    }

    /// Returns the pre-computed contact area for a specific hull pair.
    ///
    /// Returns 0.0 if no contact was found for this pair.
    pub fn get_contact_area_for_pair(&self, hull_i: usize, hull_j: usize) -> f64 {
        let (lo, hi) = if hull_i <= hull_j {
            (hull_i, hull_j)
        } else {
            (hull_j, hull_i)
        };
        self.contact_surfaces
            .iter()
            .find(|cs| cs.hull_i == lo && cs.hull_j == hi)
            .map(|cs| cs.total_area)
            .unwrap_or(0.0)
    }

    /// Returns the contact face indices for hull_i within a specific pair.
    ///
    /// Returns None if no contact was found for this pair.
    pub fn get_contact_face_indices(
        &self,
        hull_i: usize,
        hull_j: usize,
    ) -> Option<(&[usize], &[usize])> {
        let (lo, hi) = if hull_i <= hull_j {
            (hull_i, hull_j)
        } else {
            (hull_j, hull_i)
        };
        self.contact_surfaces
            .iter()
            .find(|cs| cs.hull_i == lo && cs.hull_j == hi)
            .map(|cs| {
                if hull_i <= hull_j {
                    (cs.face_indices_i.as_slice(), cs.face_indices_j.as_slice())
                } else {
                    (cs.face_indices_j.as_slice(), cs.face_indices_i.as_slice())
                }
            })
    }
}

impl std::fmt::Debug for Vessel {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        let bounds = self.get_bounds();
        f.debug_struct("Vessel")
            .field("hulls", &self.hulls.len())
            .field("tanks", &self.tanks.len())
            .field("ap", &self.ap())
            .field("fp", &self.fp())
            .field("lbp", &self.lbp())
            .field("bounds", &bounds)
            .finish()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use nalgebra::Point3;
    use parry3d_f64::shape::TriMesh;

    fn create_test_hull() -> Hull {
        let vertices = vec![
            Point3::new(0.0, -5.0, 0.0),
            Point3::new(100.0, -5.0, 0.0),
            Point3::new(100.0, 5.0, 0.0),
            Point3::new(0.0, 5.0, 0.0),
            Point3::new(0.0, -5.0, 10.0),
            Point3::new(100.0, -5.0, 10.0),
            Point3::new(100.0, 5.0, 10.0),
            Point3::new(0.0, 5.0, 10.0),
        ];
        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).unwrap();
        Hull::from_mesh(mesh)
    }

    #[test]
    fn test_vessel_bounds() {
        let hull = create_test_hull();
        let vessel = Vessel::new(hull);

        let bounds = vessel.get_bounds();
        assert!((bounds.0 - 0.0).abs() < 1e-6);
        assert!((bounds.1 - 100.0).abs() < 1e-6);
    }

    #[test]
    fn test_vessel_perpendiculars() {
        let hull = create_test_hull();
        let vessel = Vessel::new(hull);

        assert!((vessel.ap() - 0.0).abs() < 1e-6);
        assert!((vessel.fp() - 100.0).abs() < 1e-6);
        assert!((vessel.lbp() - 100.0).abs() < 1e-6);
    }
}