Solid body dims

Project.toml

@@ -1,7 +1,7 @@
 name = "stlProcess"
 uuid = "68914fc9-42cf-4b37-b06f-0b65edf9b8fa"
 authors = ["Anson <anson@ansonbiggs.com>"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"

src/stlProcess.jl

@@ -11,12 +11,14 @@ struct Properties
     inertia::Matrix{Float64}
     surface_area::Float64
     characteristic_length::Float64
+    sb_values::Vector{Float64}
 
-    function Properties(volume, center_of_gravity, inertia, surface_area, characteristic_length)
+    function Properties(volume, center_of_gravity, inertia, surface_area, characteristic_length, sb_values)
         @assert size(center_of_gravity) == (3,)
+        @assert size(sb_values) == (3,)
         @assert size(inertia) == (3, 3)
 
-        return new(volume, center_of_gravity, inertia, surface_area, characteristic_length)
+        return new(volume, center_of_gravity, inertia, surface_area, characteristic_length, sb_values)
     end
 end
 
@@ -29,6 +31,7 @@ function get_mass_properties(triangles; scale=1.0)
     x = reduce(hcat, [[v[1] .* scale for v in tri] for tri in triangles])'
     y = reduce(hcat, [[v[2] .* scale for v in tri] for tri in triangles])'
     z = reduce(hcat, [[v[3] .* scale for v in tri] for tri in triangles])'
+    points = collect(Set(reduce(vcat, [[Array(Float64.(v)) .* scale for v in tri] for tri in triangles])))
 
     function subexpression(x)
         w0, w1, w2 = x[:, 1], x[:, 2], x[:, 3]
@@ -76,16 +79,29 @@ function get_mass_properties(triangles; scale=1.0)
     inertia[1, 2] = inertia[2, 1] = -(intg[8] - volume .* center_of_gravity[1] .* center_of_gravity[2])
     inertia[2, 3] = inertia[3, 2] = -(intg[9] - volume .* center_of_gravity[2] .* center_of_gravity[3])
     inertia[1, 3] = inertia[3, 1] = -(intg[10] - volume .* center_of_gravity[3] .* center_of_gravity[1])
+    inertia = inertia / volume
 
     # https://math.stackexchange.com/questions/128991/how-to-calculate-the-area-of-a-3d-triangle
     surface_area = sum(norm.(eachrow([x0 y0 z0] - [x1 y1 z1]) .× eachrow([x1 y1 z1] - [x2 y2 z2])) / 2)
 
-    characteristic_length = calc_characteristic_length(triangles, inertia, center_of_gravity, scale)
+    characteristic_length = calc_characteristic_length(points, inertia, center_of_gravity)
 
-    return Properties(volume, center_of_gravity, inertia ./ volume, surface_area, characteristic_length)
+    sb_values = solid_body(points)
+
+    return Properties(volume, center_of_gravity, inertia, surface_area, characteristic_length, sb_values)
 end
 
-function calc_characteristic_length(triangles, inertia, center_of_gravity, scale; θtol=0.1)
+function solid_body(points)
+    pts = reduce(hcat, points)
+
+    Xsb = abs(reduce(-, extrema(@view pts[1, :])))
+    Ysb = abs(reduce(-, extrema(@view pts[2, :])))
+    Zsb = abs(reduce(-, extrema(@view pts[3, :])))
+
+    return [Xsb, Ysb, Zsb]
+end
+
+function calc_characteristic_length(points, inertia, center_of_gravity)
     """
     Calculate the Characteristic Length using eigenvectors.
 
@@ -93,27 +109,17 @@ function calc_characteristic_length(triangles, inertia, center_of_gravity, scale
     More information on the geometry used: https://math.stackexchange.com/q/100447
     """
 
-    eigs = eigvecs(inertia)
-    points = collect(Set(reduce(vcat, [[Array(Float64.(v)) .* scale for v in tri] for tri in triangles])))
-
-    function θs_calc(ref, points)
-        """
-        Calculates the angle between a reference vector and an array of points. 
-        """
-        norm_ref = norm(ref)
-        return map(
-            pt ->
-                acos(clamp(dot(ref, pt .- center_of_gravity) / (norm(pt .- center_of_gravity) * norm_ref), -1.0, 1.0)),
-            points,
-        )
-    end
-
     characteristic_points = []
     # find the characteristic points for each eigenvector
-    for eig in eachrow(eigs)
+    for eig in eachrow(eigvecs(inertia))
 
         # Find 3 points for each direction of the eigenvector
-        θs = θs_calc(eig, points)
+        θs = map(
+            point -> acos(
+                clamp(dot(eig, point .- center_of_gravity) / (norm(point .- center_of_gravity) * norm(eig)), -1.0, 1.0),
+            ),
+            points,
+        )
         sort_index = sortperm(θs)
         min_points = points[sort_index[1:3]]
         max_points = points[sort_index[(end - 2):end]]

test/runtests.jl

@@ -53,14 +53,15 @@ end
     models = Dict(
         # Inertia math: https://en.wikipedia.org/wiki/List_of_moments_of_inertia#List_of_3D_inertia_tensors
         # Properties(volume, center_of_gravity, inertia, surface_area, characteristic_length)
-        "cube.stl" => Properties(2.0^3, center, I_mat .* 2^2 / 6, 6 * 2^2, 2), # l = 2
+        "cube.stl" => Properties(2.0^3, center, I_mat .* 2^2 / 6, 6 * 2^2, 2, [2, 2, 2]), # l = 2
-        "sphere.stl" => Properties(4 / 3 * pi, center, I_mat .* 2 / 5, 4π, 2), # r = 1
+        "sphere.stl" => Properties(4 / 3 * pi, center, I_mat .* 2 / 5, 4π, 2, [2, 2, 2]), # r = 1
         "2_4_8_cuboid.stl" => Properties(
             2 * 4 * 8,
             center,
             diagm([2^2 + 4^2, 8^2 + 2^2, 8^2 + 4^2]) ./ 12,
             2(2 * 4 + 2 * 8 + 4 * 8),
             (2 + 4 + 8) / 3,
+            [2, 4, 8],
         ), # h, d, w = 2, 4, 8
         # "slender_y.stl" => Properties(
         #     10 * π * 0.1^2,
@@ -89,6 +90,7 @@ end
             @test eigvals(props.inertia) ≈ eigvals(control.inertia) atol = 0.1
             @test props.surface_area ≈ control.surface_area atol = 0.5
             @test props.characteristic_length ≈ control.characteristic_length atol = 0.5
+            @test props.sb_values ≈ control.sb_values atol = 0.01
         end
         @testset "Compare volumes with scaling for $model" begin
             for scale in 1:4