Add triangle volume function

src/geom/mod.rs

@@ -11,6 +11,13 @@ pub fn det3x3t(a: &[f64], b: &[f64], c: &[f64]) -> f64 {
     ((a[0] - c[0]) * (b[1] - c[1])) - ((a[1] - c[1]) * (b[0] - c[0]))
 }
 
+pub fn crossproduct(a: &[f64], b: &[f64]) -> Vec<f64> {
+    let i = (a[1] * b[2]) - (a[2] * b[1]);
+    let j = (a[0] * b[2]) - (a[2] * b[0]);
+    let k = (a[0] * b[1]) - (a[1] * b[0]);
+    return vec![i, -j, k];
+}
+
 pub fn area2d_triangle(a: &[f64], b: &[f64], c: &[f64]) -> f64 {
     det3x3t(a, b, c) / 2.0
 }
@@ -24,6 +31,33 @@ pub fn area3d_triangle(a: &[f64], b: &[f64], c: &[f64]) -> f64 {
     return 0.5 * (i * i + j * j + k * k).sqrt();
 }
 
+pub fn volume_triangle_from_base(a: &[f64], b: &[f64], c: &[f64], basez: f64) -> f64 {
+    //-- volume of the prism from the basez to minz
+    let minz = a[2].min(b[2]).min(c[2]);
+    let vol0 = area2d_triangle(&a, &b, &c) * (minz - basez);
+    //-- construct a polyhedron, 1+ of the vertices will be at minz
+    let a_ = vec![a[0], a[1], minz];
+    let b_ = vec![b[0], b[1], minz];
+    let c_ = vec![c[0], c[1], minz];
+    //-- split this polyhedron into 8 faces and sum their signed volume
+    let mut vol = 0.0;
+    vol += signed_volume(&a, &b, &c);
+    vol += signed_volume(&a_, &c_, &b_);
+    vol += signed_volume(&a, &a_, &b);
+    vol += signed_volume(&a_, &b_, &b);
+    vol += signed_volume(&b, &b_, &c);
+    vol += signed_volume(&b_, &c_, &c);
+    vol += signed_volume(&c, &c_, &a);
+    vol += signed_volume(&c_, &a_, &a);
+    return vol0 + vol.abs();
+}
+
+pub fn signed_volume(a: &[f64], b: &[f64], c: &[f64]) -> f64 {
+    let xp = crossproduct(&b, &c);
+    let dp = (a[0] * xp[0]) + (a[1] * xp[1]) + (a[2] * xp[2]);
+    return dp / 6.0;
+}
+
 pub fn normal_triangle(a: &[f64], b: &[f64], c: &[f64], normalise: bool) -> Vec<f64> {
     let v0: Vec<f64> = vec![a[0] - b[0], a[1] - b[1], a[2] - b[2]];
     let v1: Vec<f64> = vec![a[0] - c[0], a[1] - c[1], a[2] - c[2]];

src/lib.rs

@@ -979,6 +979,18 @@ impl Triangulation {
         }
     }
 
+    pub fn volume_triangle(&self, tr: &Triangle, planez: f64) -> Result<f64, StartinError> {
+        match self.is_triangle(tr) {
+            false => Err(StartinError::TriangleNotPresent),
+            true => Ok(geom::volume_triangle_from_base(
+                &self.stars[tr.v[0]].pt,
+                &self.stars[tr.v[1]].pt,
+                &self.stars[tr.v[2]].pt,
+                planez,
+            )),
+        }
+    }
+
     /// Returns the degree of the vertex with ID `vi`.
     pub fn degree(&self, vi: usize) -> Result<usize, StartinError> {
         match self.is_vertex_removed(vi) {