brepkit-check 4.6.8

//! Geometric properties: volume, area, center of mass, inertia tensor.

pub mod accumulator;
pub mod analytic;
pub mod bbox;
mod boundary;
pub mod face_integrator;

pub use accumulator::GProps;

use brepkit_math::aabb::Aabb3;
use brepkit_math::vec::{Point3, Vec3};
use brepkit_topology::Topology;
use brepkit_topology::face::FaceId;
use brepkit_topology::solid::SolidId;

use crate::CheckError;

/// Options for property computation.
#[derive(Debug, Clone)]
pub struct PropertiesOptions {
    /// Gauss quadrature order (default 5).
    pub gauss_order: usize,
    /// Adaptive integration tolerance (default 1e-6).
    pub adaptive_eps: f64,
    /// Maximum adaptive subdivision depth (default 8).
    pub max_depth: usize,
}

impl Default for PropertiesOptions {
    fn default() -> Self {
        Self {
            gauss_order: 5,
            adaptive_eps: 1e-6,
            max_depth: 8,
        }
    }
}

/// Compute the bounding box of a solid.
///
/// # Errors
///
/// Returns an error if any topology entity is missing or the solid has no vertices.
pub fn bounding_box(topo: &Topology, solid: SolidId) -> Result<Aabb3, CheckError> {
    bbox::bounding_box(topo, solid)
}

/// Compute the volume of a solid via face integration.
///
/// Uses the divergence theorem: V = (1/3) sum of integral (P - C) dot N dA
/// over every face of the solid, cavity shells included, with C the origin
/// for a solid near it and the centre of its box for one far away.
///
/// # Errors
///
/// Returns an error if any topology entity is missing or integration fails.
pub fn solid_volume(
    topo: &Topology,
    solid: SolidId,
    options: &PropertiesOptions,
) -> Result<f64, CheckError> {
    let about = volume_anchor(topo, solid);
    let mut total_volume = 0.0;
    for fid in brepkit_topology::explorer::solid_faces(topo, solid)? {
        let contrib = face_integrator::integrate_face_about(topo, fid, options.gauss_order, about)?;
        total_volume += contrib.volume;
    }
    Ok(total_volume)
}

/// [`solid_volume`] when every face integrates exactly, else `None`.
///
/// A plane face integrates along its edges, a cylinder, cone or sphere face
/// along its wires. `None` for a solid with a torus or NURBS face, or a face
/// neither path reads (a wire that does not chain, a cone face across its
/// apex), for a caller that would otherwise measure it from a mesh.
///
/// # Errors
///
/// Returns an error if any topology entity is missing.
pub fn exact_solid_volume(topo: &Topology, solid: SolidId) -> Result<Option<f64>, CheckError> {
    use brepkit_topology::face::FaceSurface;
    let about = volume_anchor(topo, solid);
    let mut total = 0.0;
    for fid in brepkit_topology::explorer::solid_faces(topo, solid)? {
        let face = topo.face(fid)?;
        let sign = if face.is_reversed() { -1.0 } else { 1.0 };
        let contribution = match face.surface() {
            FaceSurface::Plane { normal, .. } => {
                face_integrator::planar_face_by_edges(topo, fid, *normal * sign, about)?
            }
            FaceSurface::Cylinder(_) | FaceSurface::Cone(_) | FaceSurface::Sphere(_) => {
                boundary::curved_face_by_boundary(topo, face, sign, about)?
            }
            FaceSurface::Nurbs(_) => boundary::nurbs_face_by_boundary(topo, face, sign, about)?,
            FaceSurface::Torus(_) => None,
        };
        let Some(contribution) = contribution else {
            return Ok(None);
        };
        total += contribution.volume;
    }
    Ok(Some(total))
}

/// The point [`solid_volume`] sums about: the origin, as it always has, unless
/// the solid sits more than ten of its own half-diagonals from it, where the
/// far-off coordinates would swamp the flux. A solid without vertices (an
/// empty boolean result) has no box and sums about the origin.
fn volume_anchor(topo: &Topology, solid: SolidId) -> Vec3 {
    let Ok(b) = bbox::bounding_box(topo, solid) else {
        return Vec3::new(0.0, 0.0, 0.0);
    };
    let centre = Vec3::new(
        f64::midpoint(b.min.x(), b.max.x()),
        f64::midpoint(b.min.y(), b.max.y()),
        f64::midpoint(b.min.z(), b.max.z()),
    );
    let reach = (b.max - b.min).length() * 0.5;
    if centre.length() <= 10.0 * reach.max(1.0) {
        Vec3::new(0.0, 0.0, 0.0)
    } else {
        centre
    }
}

/// Compute the total surface area of a solid.
///
/// Sums the area of every face of the solid, cavity shells included.
///
/// # Errors
///
/// Returns an error if any topology entity is missing or integration fails.
pub fn solid_area(
    topo: &Topology,
    solid: SolidId,
    options: &PropertiesOptions,
) -> Result<f64, CheckError> {
    let mut total_area = 0.0;
    for fid in brepkit_topology::explorer::solid_faces(topo, solid)? {
        let contrib = face_integrator::integrate_face(topo, fid, options.gauss_order)?;
        total_area += contrib.area;
    }
    Ok(total_area)
}

/// Compute the center of mass of a solid.
///
/// Uses the divergence theorem: for each coordinate axis, integrates
/// `(1/2) x_i^2 * n_i` over the solid's boundary, then divides by total
/// volume to obtain the volumetric centroid (solid CoM).
///
/// # Errors
///
/// Returns an error if any topology entity is missing, integration fails,
/// or the solid has zero volume.
pub fn center_of_mass(
    topo: &Topology,
    solid: SolidId,
    options: &PropertiesOptions,
) -> Result<Point3, CheckError> {
    let mut total_volume = 0.0;
    let mut mx = 0.0;
    let mut my = 0.0;
    let mut mz = 0.0;

    for fid in brepkit_topology::explorer::solid_faces(topo, solid)? {
        let contrib = face_integrator::integrate_face(topo, fid, options.gauss_order)?;
        total_volume += contrib.volume;
        mx += contrib.volume_moment_x;
        my += contrib.volume_moment_y;
        mz += contrib.volume_moment_z;
    }

    if total_volume.abs() < 1e-30 {
        return Err(CheckError::IntegrationFailed(
            "solid has zero volume".into(),
        ));
    }

    Ok(Point3::new(
        mx / total_volume,
        my / total_volume,
        mz / total_volume,
    ))
}

/// Compute the v-range for an analytic surface by projecting face wire
/// vertices onto the given axis.
///
/// Iterates over all wires (outer + inner) of `face_id`, projects each vertex
/// position onto `axis` relative to `origin`, and returns `(v_min, v_max)`.
/// If the face has no distinguishable range (e.g. a single vertex),
/// returns `(-1.0, 1.0)` as a fallback.
///
/// # Errors
///
/// Returns an error if any topology entity is missing.
pub fn axial_v_range(
    topo: &Topology,
    face_id: FaceId,
    origin: Point3,
    axis: Vec3,
) -> Result<(f64, f64), CheckError> {
    let face_data = topo.face(face_id)?;
    let outer = topo.wire(face_data.outer_wire())?;

    let mut v_min = f64::MAX;
    let mut v_max = f64::MIN;

    // Chain outer wire and inner wires.
    let inner_wires: Vec<_> = face_data
        .inner_wires()
        .iter()
        .filter_map(|&wid| topo.wire(wid).ok())
        .collect();

    for wire in std::iter::once(outer).chain(inner_wires.iter().copied()) {
        for oe in wire.edges() {
            let edge = topo.edge(oe.edge())?;
            for vid in [oe.oriented_start(edge), oe.oriented_end(edge)] {
                let pt = topo.vertex(vid)?.point();
                let to_pt = Vec3::new(
                    pt.x() - origin.x(),
                    pt.y() - origin.y(),
                    pt.z() - origin.z(),
                );
                let v = axis.dot(to_pt);
                v_min = v_min.min(v);
                v_max = v_max.max(v);
            }
        }
    }

    if v_min < v_max {
        Ok((v_min, v_max))
    } else {
        Ok((-1.0, 1.0))
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod tests {
    use super::*;
    use brepkit_math::vec::Point3;
    use brepkit_topology::Topology;
    use brepkit_topology::test_utils::make_unit_cube_manifold;

    #[test]
    fn gprops_accumulator_two_cubes() {
        // Two unit cubes side by side along x-axis
        let a = analytic::box_props(1.0, 1.0, 1.0);
        let mut b = analytic::box_props(1.0, 1.0, 1.0);
        // Shift b's center to (1.5, 0.5, 0.5) — as if placed at x=1
        b.center = Point3::new(1.5, 0.5, 0.5);

        let mut combined = a;
        combined.add(&b);

        // Total volume = 2
        assert!((combined.mass - 2.0).abs() < 1e-12);
        // Combined center = (1.0, 0.5, 0.5)
        assert!((combined.center.x() - 1.0).abs() < 1e-12);
        assert!((combined.center.y() - 0.5).abs() < 1e-12);
        assert!((combined.center.z() - 0.5).abs() < 1e-12);
    }

    #[test]
    fn empty_solid_has_no_volume() {
        let mut topo = Topology::new();
        let empty = topo.add_empty_solid();
        let volume = solid_volume(&topo, empty, &PropertiesOptions::default()).unwrap();
        assert!(volume.abs() < f64::EPSILON);
    }

    #[test]
    fn box_props_volume_and_com() {
        let props = analytic::box_props(2.0, 3.0, 4.0);
        assert!((props.mass - 24.0).abs() < 1e-12);
        assert!((props.center.x() - 1.0).abs() < 1e-12);
        assert!((props.center.y() - 1.5).abs() < 1e-12);
        assert!((props.center.z() - 2.0).abs() < 1e-12);
        // Ixx = 24/12 * (9 + 16) = 50
        assert!((props.inertia[0] - 50.0).abs() < 1e-12);
    }

    #[test]
    fn sphere_props_volume() {
        let props = analytic::sphere_props(1.0);
        let expected = 4.0 / 3.0 * std::f64::consts::PI;
        assert!((props.mass - expected).abs() < 1e-12);
        assert!((props.center.x()).abs() < 1e-12);
        assert!((props.center.y()).abs() < 1e-12);
        assert!((props.center.z()).abs() < 1e-12);
    }

    #[test]
    fn cylinder_props_volume_and_com() {
        let props = analytic::cylinder_props(1.0, 2.0);
        let expected_v = std::f64::consts::PI * 2.0;
        assert!((props.mass - expected_v).abs() < 1e-12);
        assert!((props.center.z() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn cone_full_volume() {
        let props = analytic::cone_props(1.0, 0.0, 3.0);
        let expected_v = std::f64::consts::PI * 3.0 / 3.0; // pi * h/3 * r^2
        assert!((props.mass - expected_v).abs() < 1e-12);
        // CoM of full cone at h/4 from base
        assert!((props.center.z() - 0.75).abs() < 1e-12);
    }

    #[test]
    fn torus_props_volume() {
        let props = analytic::torus_props(3.0, 1.0);
        let expected_v = 2.0 * std::f64::consts::PI * std::f64::consts::PI * 3.0;
        assert!((props.mass - expected_v).abs() < 1e-12);
    }

    #[test]
    fn box_surface_area() {
        let area = analytic::box_area(2.0, 3.0, 4.0);
        // 2*(6 + 12 + 8) = 52
        assert!((area - 52.0).abs() < 1e-12);
    }

    #[test]
    fn sphere_surface_area() {
        let area = analytic::sphere_area(2.0);
        let expected = 4.0 * std::f64::consts::PI * 4.0;
        assert!((area - expected).abs() < 1e-12);
    }

    #[test]
    fn inertia_matrix_symmetric() {
        let mut props = GProps::new();
        props.inertia = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0];
        let mat = props.matrix_of_inertia();
        // Off-diagonal symmetry
        assert!((mat[0][1] - mat[1][0]).abs() < 1e-15);
        assert!((mat[0][2] - mat[2][0]).abs() < 1e-15);
        assert!((mat[1][2] - mat[2][1]).abs() < 1e-15);
        // Diagonal values
        assert!((mat[0][0] - 10.0).abs() < 1e-15);
        assert!((mat[1][1] - 20.0).abs() < 1e-15);
        assert!((mat[2][2] - 30.0).abs() < 1e-15);
    }

    #[test]
    fn bounding_box_unit_cube() {
        let mut topo = Topology::new();
        let solid = make_unit_cube_manifold(&mut topo);
        let aabb = bounding_box(&topo, solid).unwrap();
        // Unit cube at origin: min=(0,0,0), max=(1,1,1)
        assert!((aabb.min.x()).abs() < 1e-12);
        assert!((aabb.min.y()).abs() < 1e-12);
        assert!((aabb.min.z()).abs() < 1e-12);
        assert!((aabb.max.x() - 1.0).abs() < 1e-12);
        assert!((aabb.max.y() - 1.0).abs() < 1e-12);
        assert!((aabb.max.z() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn gauss_volume_matches_analytic() {
        let mut topo = Topology::new();
        let solid = make_unit_cube_manifold(&mut topo);
        let options = PropertiesOptions::default();
        let vol = solid_volume(&topo, solid, &options).unwrap();
        // Unit cube volume = 1.0
        assert!((vol - 1.0).abs() < 1e-10, "expected volume 1.0, got {vol}");
    }

    #[test]
    fn gauss_area_matches_analytic() {
        let mut topo = Topology::new();
        let solid = make_unit_cube_manifold(&mut topo);
        let options = PropertiesOptions::default();
        let area = solid_area(&topo, solid, &options).unwrap();
        // Unit cube surface area = 6.0
        assert!((area - 6.0).abs() < 1e-10, "expected area 6.0, got {area}");
    }

    #[test]
    fn gauss_com_matches_analytic() {
        let mut topo = Topology::new();
        let solid = make_unit_cube_manifold(&mut topo);
        let options = PropertiesOptions::default();
        let com = center_of_mass(&topo, solid, &options).unwrap();
        // Unit cube CoM at (0.5, 0.5, 0.5)
        assert!((com.x() - 0.5).abs() < 1e-10, "com.x = {}", com.x());
        assert!((com.y() - 0.5).abs() < 1e-10, "com.y = {}", com.y());
        assert!((com.z() - 0.5).abs() < 1e-10, "com.z = {}", com.z());
    }

    #[test]
    fn accumulator_default_is_zero() {
        let props = GProps::default();
        assert!((props.mass).abs() < 1e-15);
        assert!((props.center.x()).abs() < 1e-15);
        assert!((props.center.y()).abs() < 1e-15);
        assert!((props.center.z()).abs() < 1e-15);
        for &c in &props.inertia {
            assert!(c.abs() < 1e-15);
        }
    }

    /// The unit cube with its top face's surface swapped for `top`, a
    /// NURBS surface over the same unit square at z = 1.
    fn cube_topped_by(top: brepkit_math::nurbs::surface::NurbsSurface) -> (Topology, SolidId) {
        use brepkit_topology::face::FaceSurface;
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);
        let shell = topo.solid(cube).unwrap().outer_shell();
        let lid = topo
            .shell(shell)
            .unwrap()
            .faces()
            .iter()
            .copied()
            .find(|&f| {
                matches!(topo.face(f).unwrap().surface(), FaceSurface::Plane { normal, d }
                    if normal.z() > 0.5 && (d - 1.0).abs() < 1e-9)
            })
            .unwrap();
        topo.face_mut(lid)
            .unwrap()
            .set_surface(FaceSurface::Nurbs(top));
        (topo, cube)
    }

    /// A top `(100u / (1 + 99u), v, 1)`, a unit square run at a hundredth of
    /// its pace at one side, still holds the cube's 1: the integrals are
    /// halved until they settle, where a fixed rule gave 0.9445.
    #[test]
    fn a_cube_under_a_crowded_rational_top_holds_its_volume() {
        let p = Point3::new;
        let top = brepkit_math::nurbs::surface::NurbsSurface::new(
            1,
            1,
            vec![0.0, 0.0, 1.0, 1.0],
            vec![0.0, 0.0, 1.0, 1.0],
            vec![
                vec![p(0.0, 0.0, 1.0), p(0.0, 1.0, 1.0)],
                vec![p(1.0, 0.0, 1.0), p(1.0, 1.0, 1.0)],
            ],
            vec![vec![1.0, 1.0], vec![100.0, 100.0]],
        )
        .unwrap();
        let (topo, cube) = cube_topped_by(top);
        let volume = exact_solid_volume(&topo, cube).unwrap().unwrap();
        assert!((volume.abs() - 1.0).abs() < 1e-9, "{volume}");
    }

    /// A flat top whose `v` runs twice as fast past its knot at 0.25: the
    /// cube's side edges, read at eighths, land on that knot line at their
    /// middles, and the cube still holds 1.
    #[test]
    fn a_cube_under_a_top_with_a_sampled_knot_holds_its_volume() {
        let p = Point3::new;
        let top = brepkit_math::nurbs::surface::NurbsSurface::new(
            1,
            1,
            vec![0.0, 0.0, 1.0, 1.0],
            vec![0.0, 0.0, 0.25, 1.0, 1.0],
            vec![
                vec![p(0.0, 0.0, 1.0), p(0.0, 0.5, 1.0), p(0.0, 1.0, 1.0)],
                vec![p(1.0, 0.0, 1.0), p(1.0, 0.5, 1.0), p(1.0, 1.0, 1.0)],
            ],
            vec![vec![1.0; 3]; 2],
        )
        .unwrap();
        let (topo, cube) = cube_topped_by(top);
        let volume = exact_solid_volume(&topo, cube).unwrap().unwrap();
        assert!((volume.abs() - 1.0).abs() < 1e-12, "{volume}");
    }

    /// The flat NURBS patch `(u, v, 1)` over the unit square.
    fn flat_top() -> brepkit_math::nurbs::surface::NurbsSurface {
        let p = Point3::new;
        brepkit_math::nurbs::surface::NurbsSurface::new(
            1,
            1,
            vec![0.0, 0.0, 1.0, 1.0],
            vec![0.0, 0.0, 1.0, 1.0],
            vec![
                vec![p(0.0, 0.0, 1.0), p(0.0, 1.0, 1.0)],
                vec![p(1.0, 0.0, 1.0), p(1.0, 1.0, 1.0)],
            ],
            vec![vec![1.0; 2]; 2],
        )
        .unwrap()
    }

    /// A closed wire at z = 1 through `corners` in turn, each edge a line
    /// unless `curve` gives the curve of the edge from that corner.
    fn wire_through(
        topo: &mut Topology,
        corners: &[(f64, f64)],
        curve: impl Fn(usize) -> Option<brepkit_topology::edge::EdgeCurve>,
    ) -> brepkit_topology::wire::WireId {
        use brepkit_topology::edge::{Edge, EdgeCurve};
        use brepkit_topology::vertex::Vertex;
        use brepkit_topology::wire::{OrientedEdge, Wire};
        let vertices: Vec<_> = corners
            .iter()
            .map(|&(x, y)| topo.add_vertex(Vertex::new(Point3::new(x, y, 1.0), 1e-7)))
            .collect();
        let edges = (0..vertices.len())
            .map(|k| {
                let (a, b) = (vertices[k], vertices[(k + 1) % vertices.len()]);
                let e = topo.add_edge(Edge::new(a, b, curve(k).unwrap_or(EdgeCurve::Line)));
                OrientedEdge::new(e, true)
            })
            .collect();
        topo.add_wire(Wire::new(edges, true).unwrap())
    }

    /// The flat patch with a centred half-unit square hole holds three
    /// quarters of the square, and a third of that for its volume term,
    /// whichever way each of its wires runs.
    #[test]
    fn a_nurbs_face_takes_its_hole_away_whichever_way_it_runs() {
        use brepkit_topology::face::{Face, FaceSurface};
        let square = |lo: f64, hi: f64, ccw: bool| {
            let mut corners = vec![(lo, lo), (hi, lo), (hi, hi), (lo, hi)];
            if !ccw {
                corners.reverse();
            }
            corners
        };
        for (outer_ccw, hole_ccw) in [(true, false), (true, true), (false, true), (false, false)] {
            let mut topo = Topology::new();
            let outer = wire_through(&mut topo, &square(0.0, 1.0, outer_ccw), |_| None);
            let hole = wire_through(&mut topo, &square(0.25, 0.75, hole_ccw), |_| None);
            let face = Face::new(outer, vec![hole], FaceSurface::Nurbs(flat_top()));
            let c = boundary::nurbs_face_by_boundary(&topo, &face, 1.0, Vec3::new(0.0, 0.0, 0.0))
                .unwrap()
                .unwrap();
            assert!(
                (c.area - 0.75).abs() < 1e-12 && (c.volume - 0.25).abs() < 1e-12,
                "outer ccw {outer_ccw}, hole ccw {hole_ccw}: area {}, volume {}",
                c.area,
                c.volume
            );
        }
    }

    /// The flat patch with its side at u = 1 the rational line
    /// `(1, 100t / (1 + 99t), 1)`, weights 1 and 100, which runs at a
    /// hundredth of its pace at one end: a third for its volume term, read
    /// from the curve's own velocity.
    #[test]
    fn a_nurbs_face_reads_a_crowded_edge_by_its_velocity() {
        use brepkit_topology::edge::EdgeCurve;
        use brepkit_topology::face::{Face, FaceSurface};
        let rational = brepkit_math::nurbs::curve::NurbsCurve::new(
            1,
            vec![0.0, 0.0, 1.0, 1.0],
            vec![Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, 1.0)],
            vec![1.0, 100.0],
        )
        .unwrap();
        let mut topo = Topology::new();
        let corners = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)];
        let outer = wire_through(&mut topo, &corners, |k| {
            (k == 1).then(|| EdgeCurve::NurbsCurve(rational.clone()))
        });
        let face = Face::new(outer, vec![], FaceSurface::Nurbs(flat_top()));
        let c = boundary::nurbs_face_by_boundary(&topo, &face, 1.0, Vec3::new(0.0, 0.0, 0.0))
            .unwrap()
            .unwrap();
        assert!(
            (c.area - 1.0).abs() < 1e-12 && (c.volume - 1.0 / 3.0).abs() < 1e-12,
            "area {}, volume {}",
            c.area,
            c.volume
        );
    }
}