Calculate inertial parameters from a mesh

There is some package that can be used to calculate inertial parameters (centers of masses, moments of inertial, etc.) from a triangular mesh? E.g. implementation of [1] or similar.

I’ve been looking at the packages in, and there are packages to deal with meshes, but I have not found that kind of functions in particular.

[1] Kallay, M. (2006). Computing the Moment of Inertia of a Solid Defined by a Triangle Mesh. Journal of Graphics Tools, 11(2), 51-57.

That paper is behind a paywall.
The code however is on a Github repository
The code is 74 lines

I don’t understand this in AddTriangleContribution

// Contribution to the mass
_m += v;

Surely equating mass to volume assumes a density of 1
And hell’s bells guys, you are not writing in Fortran. Single letter variable names are a bit too old school.

Also everything is double precision. Really? The coordinates of some points on a mesh? How darned big is your mesh? Yes I know double precision is simply a default, I’m just questioning whether it is needed for some dimensions.

The final mass is calculated using m = _m / 6 on line 55 which implicitly assumes a density of 1/6th but density could easily be changed at this point (as long as you assume uniform density).

Also the GetResults method isn’t idempotent which seems silly for what it does. The author even includes a warning against running it twice.

Yes, there is the C++ code, and it is easy to rewrite it in Julia.
But for that very reason, I wondered if that had already been done by the people who created packages like Meshing, etc.

Thankyou @Jordan_Cluts Love that two double floats are being divided by an integer.
I know it works fine, but… densities are floating point values (hopefully positive …)
What if at a later stage someone adds the density as a constant somewhere in the code? Yeah, it will work fine. But assumptions like that lead to moon shots missing by a country mile…
I guess I am not getting it that this is an example code, intended to show you the method, so is a bit like pseudo code.

I am not sure why you think that dividing floats by integers is problematic. They are promoted automatically in most languages, and using a constant can be folded, eg in Julia:

julia> @code_warntype 12.0 / 6
 316 1 ─ %1 = Base.sitofp(Float64, %%y)::Float64                              │╻╷╷╷ promote
     │   %2 = Base.div_float(%%x, %1)::Float64                                │╻    /
     └──      return %2                                                       │    

julia> @code_warntype 12.0 / 6.0
 401 1 ─ %1 = Base.div_float(%%x, %%y)::Float64                                           │
     └──      return %1                                                                   │

julia> div6(x) = x/6
div6 (generic function with 1 method)

julia> @code_warntype div6(12.0)
 1 1 ─ %1 = Base.div_float(%%x, 6.0)::Float64                                         │╻╷ /
   └──      return %1                                                                 │  

I don’t think it is problematic. I’m just thinking out loud that using an integer for a quantity which may be replaced at some stage with a floating point value is the sort of thing that leads to bugs later in the life of software. I stand to be corrected.

Sorry - I should explain that I spent a few months as a software test engineer. I still have the facial tics.

My (admittedly superficial) understanding is that the 6 there comes from the tetrahedron volume formula. The only scenario I can imagine it changing is a fundamental redesign of the algorithm, but then the whole context would change so I would not worry about this too much.

1 Like

Thankyou. My bad. If it is a factor of 6 then it of course will never change. Not unless you have some freaky non-linear universe.

Ah you are totally right that the 6 is from the tetrahedron volume formula. I had assumed without checking that that was all accounted for in the // Signed volume of this tetrahedron (an understandable oversight I personally think) and that the modification later would be for density. So this formula does assume a density of 1 (which would obviously be trivial to change by simply multiplying by the density in the correct units of ones coordinate system).

Check out FinEtools. Examples of these calculations are e.g. in test_acoustics.jl (module mmiintegrationmm).


Thanks for answering!
The previous conversation about the bad c++ code published by JGT might be fun, but this is more useful. :slightly_smiling_face:

By the way: are you interested in the inertial properties of solids (represented by polyhedra), or “inertial” properties of surfaces (represented by polygons, presumably with some thickness)?

I’m interested in the inertial parameters of the solid enclosed by the mesh.

Sorry all. I fell short of Bill and Teds advice to “Be excellent to each other” in this discussion. Note to self - place ears in listening mode.

1 Like

How about JuliaFEM to calculate this kind of properties?

using JuliaFEM
using JuliaFEM.Preprocess
add_elements! = JuliaFEM.add_elements!


mesh = abaqus_read_mesh(joinpath("cross_section_properties", "beam.inp"))

volume_elements = create_elements(mesh, "BEAM")
left_elements = create_surface_elements(mesh, "LEFT")

info("Number of volume elements: ", length(volume_elements))
info("Number of surface elements in LEFT", length(left_elements))

update!(volume_elements, "density", 5.0)
volume = 0.0
mass = 0.0

# Calculate volume and mass of body

time = 0.0
for element in volume_elements
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        volume += ip.weight * detJ
        density = element("density", ip, time)
        mass += ip.weight * density * detJ

info("Volume of mesh: ", round(volume, 2))
info("Mass of mesh: ", round(mass, 2))

b = 10
h = 20
L = 100

using Base.Test
@test isapprox(volume, b*h*L)
@test isapprox(mass, 5*volume)

# Calculate first moment of inertia with respect to origin

Qy = 0.0
Qz = 0.0
area = 0.0
for element in left_elements
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        x, y, z = element("geometry", ip, time)
        Qy += ip.weight * z * detJ
        Qz += ip.weight * y * detJ
        area += ip.weight * detJ

info("Qy: ", round(Qy, 2))
info("Qz: ", round(Qz, 2))

# Calculate centroid

yc = Qz/area
zc = Qy/area

info("yc: ", round(yc, 2))
info("zc: ", round(zc, 2))

@test isapprox(yc, b/2)
@test isapprox(zc, h/2)

# Calculate second moment of inertia with respect to centroid

Iy = 0.0
Iz = 0.0
for element in left_elements
    for ip in get_integration_points(element, 1)
        detJ = element(ip, time, Val{:detJ})
        x, y, z = element("geometry", ip, time)
        Iy += ip.weight * (y-yc)^2 * detJ
        Iz += ip.weight * (z-zc)^2 * detJ
        area += ip.weight * detJ

info("Iy: ", round(Iy, 2))
info("Iz: ", round(Iz, 2))

@test isapprox(Iy, h*b^3/12)
@test isapprox(Iz, b*h^3/12)

1 Like