add tests

src/composite.jl

@@ -23,7 +23,7 @@ sum_inertia(I...) = sum(I)
 # centre of mass
 centre_of_mass(ms, rs) = reduce(+, ms .* rs) / sum(ms)
 
-centre_of_mass(solids::Vector{AbstractSolid}, rs) = centre_of_mass(getfield.(solids, :mass), rs)
+centre_of_mass(solids::Vector{T}, rs) where {T<:AbstractSolid} = centre_of_mass(getfield.(solids, :mass), rs)
 
 # calculate total inertia of collection of solids about a point p
 function calculate_inertia(Is, Rs, ms, rs, p)
@@ -32,7 +32,7 @@ function calculate_inertia(Is, Rs, ms, rs, p)
 
     for (I, R, m, r) in zip(Is, Rs, ms, rs)
         # inertia about axes parallel to global frame
-        I′ = R * I * R'
+        I′ = rotated_inertia(I, R')
 
         # inertia about com
         I′cm = parallel_axis(I′, m, p - r)
@@ -46,6 +46,18 @@ end
 # total inertia about CoM
 calculate_inertia(Is, Rs, ms, rs) = calculate_inertia(Is, Rs, ms, rs, centre_of_mass(ms, rs))
 
-calculate_inertia(solids::Vector{AbstractSolid}, Rs, rs, p) = begin
+"""
+Calculate the combined inertia of a collection of rigid bodies about a point.
+
+- `solids`: vector of solids.
+- `Rs`: vector of rotation matricies for the rotation from the global frame to the body-fixed frame.
+- `rs`: vector of positions of each solids centre of mass.
+- `p`: point about which inertia is calculated. If omitted defaults to the centre of mass of the system.
+"""
+calculate_inertia(solids::Vector{T}, Rs, rs, p) where {T<:AbstractSolid} = begin
     calculate_inertia(getfield.(solids, :moi), Rs, getfield.(solids, :mass), rs, p)
 end
+
+calculate_inertia(solids::Vector{T}, Rs, rs) where {T<:AbstractSolid} = begin
+    calculate_inertia(getfield.(solids, :moi), Rs, getfield.(solids, :mass), rs)
+end

src/solids.jl

@@ -7,47 +7,49 @@ struct Cuboid{T} <: AbstractSolid
     width::T
     height::T
     mass::T
-    moi::SMatrix{3,T} # principle inertia tensor
+    moi::SMatrix{3,3,T,9} # principle inertia tensor
+end
 
-    # determine moment of inertia from dimensions and mass
-    function Cuboid{T}(l, w, h, m) where {T}
-        c = 1 / 12 * m
-        I = @SMatrix [c*(w^2+h^2) 0 0; 0 c*(l^2+h^2) 0; 0 0 c*(l^2+w^2)]
+# determine moment of inertia from dimensions and mass
+function Cuboid(l, w, h, m)
+    c = 1 / 12 * m
+    I = @SMatrix [c*(w^2+h^2) 0 0; 0 c*(l^2+h^2) 0; 0 0 c*(l^2+w^2)]
 
-        return new{T}(l, w, h, m, I)
-    end
+    return Cuboid(l, w, h, m, I)
 end
 
+
 ## sphere
 struct Sphere{T} <: AbstractSolid
     radius::T
     mass::T
-    moi::SMatrix{3,T} # principle inertia tensor
+    moi::SMatrix{3,3,T,9} # principle inertia tensor
+end
 
-    # determine moment of inertia from dimensions and mass
-    function Sphere{T}(r, m) where {T}
-        c = 2 / 5 * m
-        I = @SMatrix [c*r^2 0 0; 0 c*r^2 0; 0 0 c*r^2]
+# determine moment of inertia from dimensions and mass
+function Sphere(r, m)
+    c = 2 / 5 * m
+    I = @SMatrix [c*r^2 0 0; 0 c*r^2 0; 0 0 c*r^2]
 
-        return new{T}(r, m, I)
-    end
+    return Sphere(r, m, I)
 end
 
+
 ## cylinder
 struct Cylinder{T} <: AbstractSolid
     radius::T
     length::T
     mass::T
-    moi::SMatrix{3,T} # principle inertia tensor
+    moi::SMatrix{3,3,T,9} # principle inertia tensor
+end
 
-    # determine moment of inertia from dimensions and mass
-    function Cylinder{T}(r, l, m) where {T}
-        c1 = 1 / 12 * m
-        c2 = 1 / 2 * m
-        I = @SMatrix [c2*r^2 0 0; 0 c1*(3r^2+l^2) 0; 0 0 c1*(3r^2+l^2)]
+# determine moment of inertia from dimensions and mass
+function Cylinder(r, l, m)
+    c1 = 1 / 12 * m
+    c2 = 1 / 2 * m
+    I = @SMatrix [c2*r^2 0 0; 0 c1*(3r^2+l^2) 0; 0 0 c1*(3r^2+l^2)]
 
-        return new{T}(r, l, m, I)
-    end
+    return Cylinder(r, l, m, I)
 end
 
 

test/Manifest.toml

@@ -0,0 +1,78 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.8.2"
+manifest_format = "2.0"
+project_hash = "9e1473ed2095b1bff92020688eb7650655ed2d6c"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "0.5.2+0"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.20+0"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
+git-tree-sha1 = "f86b3a049e5d05227b10e15dbb315c5b90f14988"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.5.9"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.0"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.1.1+0"

test/Project.toml

@@ -0,0 +1,3 @@
+[deps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -0,0 +1,75 @@
+using Test, Inertia
+using StaticArrays
+
+@testset "solids" verbose = true begin
+
+    @testset "cuboid" verbose = true begin
+        cuboid = Cuboid(1.0, 2.0, 3.0, 0.1)
+        @test isapprox(cuboid.moi, @SMatrix [0.10833 0 0; 0 0.08333 0; 0 0 0.04167]; atol=1e-5)
+
+    end
+
+    @testset "sphere" verbose = true begin
+        sphere = Sphere(1.0, 1.0)
+        @test isapprox(sphere.moi, @SMatrix [0.4 0 0; 0 0.4 0; 0 0 0.4]; atol=1e-5)
+
+    end
+
+
+    @testset "cylinder" verbose = true begin
+        cylinder = Cylinder(0.1, 2.3, 3.5)
+        @test isapprox(cylinder.moi, @SMatrix [0.0175 0 0; 0 1.5517 0; 0 0 1.5517]; atol=1e-4)
+
+    end
+
+end
+
+
+# test that two cuboids next to each other such that their inertia about combined centre of mass is the same as what is would be for a larger cuboid
+@testset "composite shapes" verbose = true begin
+    # cuboids of same dimensions and mass
+    l = 1.0
+    w = 2.0
+    h = 3.0
+    m = 0.1
+    cuboid1 = Cuboid(l, w, h, m)
+    cuboid2 = Cuboid(l, w, h, m)
+
+    # dispalcement from origin of cuboids com
+    r1 = [-l / 2, 0, 0]
+    r2 = [l / 2, 0, 0]
+
+    solids = [cuboid1, cuboid2]
+    rs = [r1, r2]
+
+    @testset "centre of mass" verbose = true begin
+        com = centre_of_mass(solids, rs)
+        @test com ≈ [0.0, 0.0, 0.0]
+
+    end
+
+    @testset "parallel axis" verbose = true begin
+        @test Inertia.parallel_axis(cuboid1, [0, 0, 0]) == cuboid1.moi
+        @test Inertia.parallel_axis(cuboid1, [1, 1, 1]) == cuboid1.moi + m * [2 -1 -1; -1 2 -1; -1 -1 2]
+
+    end
+
+    @testset "rotate inertia" verbose = true begin
+        R = [1 0 0; 0 0 -1; 0 1 0] # rotates 90 degrees about x axis
+        @test Inertia.rotated_inertia(cuboid1.moi, R') == Cuboid(l, h, w, m).moi
+
+    end
+
+    @testset "moment of inertia" verbose = true begin
+        Rs = fill([1 0 0; 0 1 0; 0 0 1], 2)
+        moi = calculate_inertia(solids, Rs, rs)
+        @test moi == Cuboid(2l, w, h, 2m).moi
+
+        R = [1 0 0; 0 0 -1; 0 1 0] # rotates 90 degrees about x axis
+        Rs = fill(R, 2)
+        moi = calculate_inertia(solids, Rs, rs)
+        @test moi == Cuboid(2l, h, w, 2m).moi
+
+
+    end
+end