1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

use crate::mass_properties::MassProperties;
use crate::math::{Point, Real};
use crate::shape::Triangle;

impl MassProperties {
    /// Computes the mass properties of a triangle-mesh.
    pub fn from_trimesh(
        density: Real,
        vertices: &[Point<Real>],
        indices: &[[u32; 3]],
    ) -> MassProperties {
        let (area, com) = trimesh_area_and_center_of_mass(vertices, indices);

        if area == 0.0 {
            return MassProperties::new(com, 0.0, 0.0);
        }

        let mut itot = 0.0;

        for idx in indices {
            let triangle = Triangle::new(
                vertices[idx[0] as usize],
                vertices[idx[1] as usize],
                vertices[idx[2] as usize],
            );

            // TODO: is the parallel axis theorem correctly applied here?
            let area = triangle.area();
            let ipart = triangle.unit_angular_inertia();
            itot += ipart * area;
        }

        Self::new(com, area * density, itot * density)
    }
}

/// Computes the area and center-of-mass of a triangle-mesh.
pub fn trimesh_area_and_center_of_mass(
    vertices: &[Point<Real>],
    indices: &[[u32; 3]],
) -> (Real, Point<Real>) {
    let mut res = Point::origin();
    let mut areasum = 0.0;

    for idx in indices {
        let triangle = Triangle::new(
            vertices[idx[0] as usize],
            vertices[idx[1] as usize],
            vertices[idx[2] as usize],
        );
        let area = triangle.area();
        let center = triangle.center();

        res += center.coords * area;
        areasum += area;
    }

    if areasum == 0.0 {
        (areasum, res)
    } else {
        (areasum, res / areasum)
    }
}