brep_kernel/props/mass_properties/

solid_props.rs

use super::*;

pub fn solid_mass_properties(solid: &BrepSolid) -> Result<MassProperties, String> {
    let mut properties = MassProperties {
        surface_area: 0.0,
        volume: 0.0,
    };
    let mut volume_compensation = 0.0;
    for shell in &solid.shells {
        if shell.faces.is_empty() {
            continue;
        }
        // Use the same local reference as the signed-volume orientation gate.
        // World-origin flux terms can dwarf a small translated solid's volume;
        // even tiny quadrature or boundary errors then dominate their sum.
        let reference = shell_volume_reference(shell)?;
        let kinds = [Integrand::Area, Integrand::VolumeAbout(reference)];
        for face in &shell.faces {
            let (area, volume) = if let Some(values) = biperiodic_band_integral(face, &kinds)? {
                (values[0], values[1] / 3.0)
            } else if !is_affine(&face.surface)? && !is_untrimmed(face)? {
                // Preserve the shared cell decomposition and integration pass.
                let values = integrate_trimmed_multi(face, &kinds)?;
                (values[0], values[1] / 3.0)
            } else {
                (
                    face_area(face)?,
                    face_volume_contribution_about(face, reference)?,
                )
            };
            properties.surface_area += area;
            let next = properties.volume + volume;
            if properties.volume.abs() >= volume.abs() {
                volume_compensation += (properties.volume - next) + volume;
            } else {
                volume_compensation += (volume - next) + properties.volume;
            }
            properties.volume = next;
        }
    }
    properties.volume += volume_compensation;
    Ok(properties)
}

/// Exact signed volume only — same per-face integration paths as
/// `solid_mass_properties` but without the surface-area pass. The boolean
/// assembly orientation gate consumes only the volume sign, and the area
/// integral costs as much again as the volume one.
pub fn solid_signed_volume(solid: &BrepSolid) -> Result<f64, String> {
    let mut volume = 0.0;
    for shell in &solid.shells {
        volume += shell_signed_volume(shell)?;
    }
    Ok(volume)
}

pub(super) fn shell_volume_reference(shell: &ShellRecord) -> Result<Vec3, String> {
    let finite = |point: Vec3| point.x.is_finite() && point.y.is_finite() && point.z.is_finite();
    for face in &shell.faces {
        if let Some(coedge) = face
            .loops
            .first()
            .and_then(|loop_record| loop_record.coedges.first())
        {
            if let Ok(domain) = coedge.pcurve.domain() {
                if let Ok(uv) = coedge.pcurve.evaluate(domain[0]) {
                    if let Ok(point) = face.surface.evaluate(uv.x, uv.y) {
                        if finite(point) {
                            return Ok(point);
                        }
                    }
                }
            }
        }
        if let Ok((u_breaks, v_breaks)) = surface_breaks(&face.surface) {
            let u = 0.5 * (u_breaks[0] + u_breaks[u_breaks.len() - 1]);
            let v = 0.5 * (v_breaks[0] + v_breaks[v_breaks.len() - 1]);
            if let Ok(point) = face.surface.evaluate(u, v) {
                if finite(point) {
                    return Ok(point);
                }
            }
        }
        for row in &face.surface.control_points {
            for control in row {
                if let Ok(point) = control.point() {
                    if finite(point) {
                        return Ok(point);
                    }
                }
            }
        }
    }
    Err("mass_properties: shell has no finite geometric reference".to_string())
}

/// Exact signed volume of one closed shell. Multi-shell tessellation uses this
/// to preserve the authored material orientation: exterior shells contribute
/// positively, while a void boundary contributes negatively.
pub(crate) fn shell_signed_volume(shell: &ShellRecord) -> Result<f64, String> {
    let reference = shell_volume_reference(shell)?;

    // A closed shell's divergence-theorem volume is independent of origin,
    // but evaluating each face about the world origin can cancel enormous
    // translated face terms down to a tiny cavity volume. Anchor the
    // integrand on an authored boundary point and compensate the remaining
    // face sum so the sign is stable for small, far-translated shells.
    let mut volume = 0.0;
    let mut compensation = 0.0;
    for face in &shell.faces {
        let contribution = face_volume_contribution_about(face, reference)?;
        let next = volume + contribution;
        if volume.abs() >= contribution.abs() {
            compensation += (volume - next) + contribution;
        } else {
            compensation += (contribution - next) + volume;
        }
        volume = next;
    }
    Ok(volume + compensation)
}

/// Exact moments for affine faces via Green's theorem over the trim
/// pcurves: ∬ g du dv = ∮ G dv with G(u,v) = ∫ g dt.  On an affine carrier
/// every moment integrand is a low-degree polynomial, so the inner Gauss
/// antiderivative is exact and the boundary quadrature has the same quality
/// as `parameter_space_area` — no trim-polygon sampling error.  The loop
/// winding supplies the orientation sign, matching the affine volume path.
pub(super) fn affine_moment(face: &FaceRecord, kind: Integrand) -> Result<f64, String> {
    let ku = crate::KnotVector::new(face.surface.knots_u.clone(), face.surface.degree_u)?;
    let kv = crate::KnotVector::new(face.surface.knots_v.clone(), face.surface.degree_v)?;
    let [u0, u1] = ku.domain();
    let [v0, v1] = kv.domain();
    let points = &face.surface.control_points;
    let p00 = points[0][0].point()?;
    let p10 = points[1][0].point()?;
    let p01 = points[0][1].point()?;
    let du = p10.sub(p00).scale(1.0 / (u1 - u0));
    let dv = p01.sub(p00).scale(1.0 / (v1 - v0));
    let weighted_normal = du.cross(dv);
    let g = |u: f64, v: f64| {
        let point = p00.add(du.scale(u - u0)).add(dv.scale(v - v0));
        integrand_value(kind, point, weighted_normal)
    };
    let inner = |u: f64, v: f64| {
        let half = (u - u0) * 0.5;
        let middle = (u + u0) * 0.5;
        let mut sum = 0.0;
        for index in 0..GAUSS_X.len() {
            sum += GAUSS_W[index] * g(middle + half * GAUSS_X[index], v);
        }
        sum * half
    };
    let mut total = 0.0;
    for loop_record in &face.loops {
        for coedge in &loop_record.coedges {
            for pair in curve_breaks(&coedge.pcurve)?.windows(2) {
                let half = (pair[1] - pair[0]) * 0.5;
                let middle = (pair[1] + pair[0]) * 0.5;
                for index in 0..GAUSS_X.len() {
                    let parameter = middle + half * GAUSS_X[index];
                    let (point, tangent) = coedge.pcurve.deriv1(parameter)?;
                    total += GAUSS_W[index] * half * inner(point.x, point.y) * tangent.y;
                }
            }
        }
    }
    Ok(total)
}

pub(super) fn face_moment(face: &FaceRecord, kind: Integrand) -> Result<f64, String> {
    if is_affine(&face.surface)? {
        return affine_moment(face, kind);
    }
    if let Some(values) = biperiodic_band_integral(face, &[kind])? {
        return Ok(values[0]);
    }
    if is_untrimmed(face)? {
        integrate_untrimmed(face, kind)
    } else {
        integrate_trimmed(face, kind)
    }
}

/// Symmetric-matrix product `A[i][j] = M[k][i] * M[k][j]` reused by the Jacobi
/// sweep — kept tiny and explicit rather than pulling in a matrix crate.
pub(super) fn mat3_mul(a: [[f64; 3]; 3], b: [[f64; 3]; 3]) -> [[f64; 3]; 3] {
    let mut r = [[0.0f64; 3]; 3];
    for i in 0..3 {
        for j in 0..3 {
            for k in 0..3 {
                r[i][j] += a[i][k] * b[k][j];
            }
        }
    }
    r
}

pub(super) fn mat3_transpose(a: [[f64; 3]; 3]) -> [[f64; 3]; 3] {
    let mut r = [[0.0f64; 3]; 3];
    for i in 0..3 {
        for j in 0..3 {
            r[i][j] = a[j][i];
        }
    }
    r
}

/// Classic Jacobi eigenvalue iteration for a SYMMETRIC 3×3 matrix. Repeatedly
/// applies a Givens rotation in the plane of the largest off-diagonal element,
/// each rotation zeroing that element, until the matrix is diagonal to machine
/// precision. Returns `(eigenvalues, vectors)` where `vectors` holds the
/// eigenvectors as COLUMNS (`vectors[i][k]` is component `i` of eigenvector
/// `k`), unsorted. For a symmetric matrix Jacobi is unconditionally
/// convergent and the accumulated rotations stay orthonormal, so the columns
/// are mutually orthonormal by construction.
pub(super) fn jacobi_eigen_symmetric_3x3(matrix: [[f64; 3]; 3]) -> ([f64; 3], [[f64; 3]; 3]) {
    // Symmetrize defensively against tiny asymmetry from round-off upstream.
    let mut a = [
        [matrix[0][0], 0.0, 0.0],
        [0.0, matrix[1][1], 0.0],
        [0.0, 0.0, matrix[2][2]],
    ];
    a[0][1] = 0.5 * (matrix[0][1] + matrix[1][0]);
    a[1][0] = a[0][1];
    a[0][2] = 0.5 * (matrix[0][2] + matrix[2][0]);
    a[2][0] = a[0][2];
    a[1][2] = 0.5 * (matrix[1][2] + matrix[2][1]);
    a[2][1] = a[1][2];

    let mut v = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
    let scale = a[0][0].abs() + a[1][1].abs() + a[2][2].abs() + 1.0;
    for _sweep in 0..64 {
        // Pick the largest off-diagonal magnitude.
        let pairs = [(0usize, 1usize), (0, 2), (1, 2)];
        let (mut p, mut q, mut best) = (0usize, 1usize, 0.0f64);
        for &(i, j) in &pairs {
            if a[i][j].abs() > best {
                best = a[i][j].abs();
                p = i;
                q = j;
            }
        }
        if best <= 1e-18 * scale {
            break;
        }
        // Angle that zeroes a[p][q]: t = tan(theta) is the smaller root of
        // t² + 2·theta·t − 1 = 0 with theta = (a_qq − a_pp)/(2 a_pq).
        let theta = (a[q][q] - a[p][p]) / (2.0 * a[p][q]);
        let t = if theta == 0.0 {
            1.0
        } else {
            theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt())
        };
        let c = 1.0 / (t * t + 1.0).sqrt();
        let s = t * c;
        // Givens rotation J: J[p][p]=J[q][q]=c, J[p][q]=s, J[q][p]=−s.
        let mut j = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
        j[p][p] = c;
        j[q][q] = c;
        j[p][q] = s;
        j[q][p] = -s;
        // A ← Jᵀ A J   (drives a[p][q] to zero), V ← V J.
        a = mat3_mul(mat3_transpose(j), mat3_mul(a, j));
        v = mat3_mul(v, j);
    }
    ([a[0][0], a[1][1], a[2][2]], v)
}

/// Principal axes/moments from a centroidal inertia tensor (Golovanov §8.11).
/// Diagonalizes the symmetric tensor, sorts the eigenpairs by moment
/// ASCENDING, returns each eigenvector as a ROW (`axes[i]` pairs with
/// `moments[i]`) normalized to unit length, and fixes the sign of the third
/// axis so the frame is right-handed (determinant +1).
pub(super) fn principal_frame(inertia: [[f64; 3]; 3]) -> ([f64; 3], [[f64; 3]; 3]) {
    let (values, vectors) = jacobi_eigen_symmetric_3x3(inertia);
    // Column k of `vectors` is the eigenvector for `values[k]`.
    let mut order = [0usize, 1, 2];
    order.sort_by(|&a, &b| values[a].total_cmp(&values[b]));
    let mut moments = [0.0f64; 3];
    let mut axes = [[0.0f64; 3]; 3];
    for (slot, &k) in order.iter().enumerate() {
        moments[slot] = values[k];
        let mut axis = [vectors[0][k], vectors[1][k], vectors[2][k]];
        let length = (axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]).sqrt();
        if length > 0.0 {
            axis = [axis[0] / length, axis[1] / length, axis[2] / length];
        }
        axes[slot] = axis;
    }
    // Right-hand the frame: if axis0 × axis1 points opposite axis2, flip axis2.
    let cross = [
        axes[0][1] * axes[1][2] - axes[0][2] * axes[1][1],
        axes[0][2] * axes[1][0] - axes[0][0] * axes[1][2],
        axes[0][0] * axes[1][1] - axes[0][1] * axes[1][0],
    ];
    let det = cross[0] * axes[2][0] + cross[1] * axes[2][1] + cross[2] * axes[2][2];
    if det < 0.0 {
        axes[2] = [-axes[2][0], -axes[2][1], -axes[2][2]];
    }
    (moments, axes)
}

/// Area, volume, centroid, and centroidal inertia (unit density).  Area and
/// volume use the same exact paths as `solid_mass_properties`; the moment
/// integrals use divergence-theorem surface quadrature (exact for untrimmed
/// spans, trim-polygon scanline accuracy for trimmed faces).
pub fn solid_mass_properties_full(solid: &BrepSolid) -> Result<FullMassProperties, String> {
    let base = solid_mass_properties(solid)?;
    const MOMENT_KINDS: [Integrand; 9] = [
        Integrand::MomentX,
        Integrand::MomentY,
        Integrand::MomentZ,
        Integrand::SecondXX,
        Integrand::SecondYY,
        Integrand::SecondZZ,
        Integrand::ProductXY,
        Integrand::ProductXZ,
        Integrand::ProductYZ,
    ];
    let mut moments = [0.0f64; 3];
    let mut seconds = [0.0f64; 3];
    let mut products = [0.0f64; 3];
    for shell in &solid.shells {
        for face in &shell.faces {
            let results = if !is_affine(&face.surface)? && !is_untrimmed(face)? {
                // One cell decomposition, all nine moment integrands per
                // station, instead of nine full passes per face.
                integrate_trimmed_multi(face, &MOMENT_KINDS)?
            } else {
                MOMENT_KINDS
                    .iter()
                    .map(|kind| face_moment(face, *kind))
                    .collect::<Result<Vec<_>, _>>()?
            };
            moments[0] += results[0];
            moments[1] += results[1];
            moments[2] += results[2];
            seconds[0] += results[3];
            seconds[1] += results[4];
            seconds[2] += results[5];
            products[0] += results[6];
            products[1] += results[7];
            products[2] += results[8];
        }
    }
    let volume = base.volume;
    if volume.abs() <= 1e-30 {
        return Err("solid_mass_properties_full: non-positive volume".into());
    }
    let centroid = Vec3::new(
        moments[0] / volume,
        moments[1] / volume,
        moments[2] / volume,
    );
    // Inertia about the origin, then parallel-axis down to the centroid.
    let ixx =
        seconds[1] + seconds[2] - volume * (centroid.y * centroid.y + centroid.z * centroid.z);
    let iyy =
        seconds[0] + seconds[2] - volume * (centroid.x * centroid.x + centroid.z * centroid.z);
    let izz =
        seconds[0] + seconds[1] - volume * (centroid.x * centroid.x + centroid.y * centroid.y);
    let ixy = -(products[0] - volume * centroid.x * centroid.y);
    let ixz = -(products[1] - volume * centroid.x * centroid.z);
    let iyz = -(products[2] - volume * centroid.y * centroid.z);
    let inertia = [[ixx, ixy, ixz], [ixy, iyy, iyz], [ixz, iyz, izz]];
    let (principal_moments, principal_axes) = principal_frame(inertia);
    Ok(FullMassProperties {
        surface_area: base.surface_area,
        volume,
        centroid,
        inertia,
        principal_moments,
        principal_axes,
    })
}