# Calculate inertial parameters from a mesh

**URL:** <https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633>\
**Category:** General Usage\
**Tags:** package\
**Created:** [July 24, 2018, 10:45am UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633 "2018-07-24T10:45:47Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [July 24, 2018, 10:45am UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/1 "2018-07-24T10:45:47Z")

</div>

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 [https://github.com/JuliaGeometry](https://github.com/JuliaGeometry), 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. [https://doi.org/10.1080/2151237X.2006.10129220](https://doi.org/10.1080/2151237X.2006.10129220)

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [July 24, 2018, 2:21pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/2 "2018-07-24T14:21:19Z")

</div>

That paper is behind a paywall.  
The code however is on a Github repository [https://github.com/erich666/jgt-code](https://github.com/erich666/jgt-code)  
The code is 74 lines  
[https://github.com/erich666/jgt-code/blob/master/Volume\_11/Number\_2/Kallay2006/Moment\_of\_Inertia.cpp](https://github.com/erich666/jgt-code/blob/master/Volume_11/Number_2/Kallay2006/Moment_of_Inertia.cpp)

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [July 24, 2018, 2:45pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/3 "2018-07-24T14:45:44Z")

</div>

I don’t understand this in AddTriangleContribution

```julia
// 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.

---

<div class="post-metadata">

**Author:** ![Jordan\_Cluts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordan_cluts/32/13753_2.png) [@Jordan\_Cluts](https://discourse.julialang.org/u/Jordan_Cluts)\
**Post date:** [July 24, 2018, 2:48pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/4 "2018-07-24T14:48:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [July 24, 2018, 2:49pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/5 "2018-07-24T14:49:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [July 24, 2018, 2:56pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/6 "2018-07-24T14:56:30Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [July 24, 2018, 3:24pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/7 "2018-07-24T15:24:19Z")

</div>

> [@johnh](#):
>
> Love that two double floats are being divided by an integer.  
> I know it works fine, but… densities are floating point values (hopefully positive …)

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
julia> @code_warntype 12.0 / 6
Body::Float64
 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
Body::Float64
 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)
Body::Float64
 1 1 ─ %1 = Base.div_float(%%x, 6.0)::Float64 │╻╷ /
   └── return %1 │  

```

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [July 24, 2018, 3:29pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/8 "2018-07-24T15:29:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [July 24, 2018, 3:40pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/9 "2018-07-24T15:40:55Z")

</div>

> [@johnh](#):
>
> using an integer for a quantity which may be replaced at some stage with a floating point value

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.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [July 24, 2018, 3:43pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/10 "2018-07-24T15:43:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Jordan\_Cluts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordan_cluts/32/13753_2.png) [@Jordan\_Cluts](https://discourse.julialang.org/u/Jordan_Cluts)\
**Post date:** [July 24, 2018, 3:45pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/11 "2018-07-24T15:45:36Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [July 24, 2018, 5:05pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/12 "2018-07-24T17:05:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [July 24, 2018, 5:57pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/13 "2018-07-24T17:57:00Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [July 24, 2018, 6:13pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/14 "2018-07-24T18:13:50Z")

</div>

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)?

---

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [July 24, 2018, 6:18pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/15 "2018-07-24T18:18:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [July 25, 2018, 8:32am UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/16 "2018-07-25T08:32:12Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ahojukka5](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahojukka5/32/4376_2.png) [@ahojukka5](https://discourse.julialang.org/u/ahojukka5)\
**Post date:** [July 25, 2018, 1:02pm UTC](https://discourse.julialang.org/t/calculate-inertial-parameters-from-a-mesh/12633/17 "2018-07-25T13:02:43Z")

</div>

How about JuliaFEM to calculate this kind of properties?

 ![beam](https://global.discourse-cdn.com/julialang/original/3X/e/d/ed516edc885eef276ba2c050398802904bdc8943.png)

```julia
using JuliaFEM
using JuliaFEM.Preprocess
add_elements! = JuliaFEM.add_elements!

# https://mechanicalc.com/reference/cross-sections

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
    end
end

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
    end
end

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
    end
end

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

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

```