new volume function based on mass props

test/runtests.jl

@@ -7,6 +7,44 @@ using MeshIO
 
 using LinearAlgebra
 
+function _check_volume(triangles; scale=1)
+    """
+    Slow algorithm just used to test the other algorithms
+    Reference:
+        https://people.eecs.berkeley.edu/~wkahan/VtetLang.pdf
+    """
+
+    volume = 0.0
+    for tri in triangles
+        a, b, c = tri .* scale
+
+        W = sqrt(sum((a .- b) .^ 2))
+        V = sqrt(sum((b .- c) .^ 2))
+        U = sqrt(sum((c .- a) .^ 2))
+
+        v = sqrt(sum(a .^ 2))
+        u = sqrt(sum(b .^ 2))
+        w = sqrt(sum(c .^ 2))
+
+        X = (w - U + v) * (U + v + w)
+        Y = (u - V + w) * (V + w + u)
+        Z = (v - W + u) * (W + u + v)
+
+        x = (U - v + w) * (v - w + U)
+        y = (V - w + u) * (w - u + V)
+        z = (W - u + v) * (u - v + W)
+
+        ξ = sqrt(x * Y * Z)
+        η = sqrt(y * Z * X)
+        ζ = sqrt(z * X * Y)
+        λ = sqrt(x * y * z)
+
+        volume += sqrt((ξ + η + ζ - λ) * (λ + ξ + η - ζ) * (η + ζ + λ - ξ) * (ζ + λ + ξ - η)) / (192 * u * v * w)
+    end
+
+    return volume
+end
+
 @testset "simple cube stl" begin
     """
     inertia math:
@@ -25,14 +63,14 @@ using LinearAlgebra
         @test all(eigvals(props.inertia) .≈ (5.0 + 1 / 3))
 
         @testset "get_volume function" begin
-            @test get_volume(stl; scale=1) ≈ props.volume rtol = 0.01
-            @test get_volume(stl; scale=4) ≈ (4 * 2)^3 rtol = 0.01
+            @test _check_volume(stl; scale=1) ≈ props.volume rtol = 0.01
+            @test _check_volume(stl; scale=4) ≈ (4 * 2)^3 rtol = 0.01
         end
     end
     @testset "Compare volumes with scaling" begin
         for scale in 1:5
             props = get_mass_properties(stl; scale=scale)
-            volume = get_volume(stl; scale=scale)
+            volume = _check_volume(stl; scale=scale)
 
             @test props.volume ≈ volume rtol = 0.01
         end