diff --git a/Project.toml b/Project.toml index eee994f..b675e89 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "stlProcess" uuid = "68914fc9-42cf-4b37-b06f-0b65edf9b8fa" authors = ["Anson "] -version = "0.2.1" +version = "0.2.2" [deps] FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" diff --git a/src/stlProcess.jl b/src/stlProcess.jl index bfd3798..e1f7137 100644 --- a/src/stlProcess.jl +++ b/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]] diff --git a/test/runtests.jl b/test/runtests.jl index afb5ace..91117f7 100644 --- a/test/runtests.jl +++ b/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 - "sphere.stl" => Properties(4 / 3 * pi, center, I_mat .* 2 / 5, 4π, 2), # r = 1 + "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, [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