Inertia confusion

[JoeJ]

inertia.JPG (67.1 KiB) Viewed 6790 times

On the right side of the image there is a long box with inertia visulaized in red about any axis (forming an ellipsoid as expected).

On the left side, the same box is composed from 7 cubes, but the mass properties should be exactly the same.
I calculate inertia about an axis by iterating and summing up over all 7 bodies, code below.

On the very right of the image the numbers show that on principal axis the inarta values match for both,
but the visualization shows a more doughnut like inertia.


Question: Does the doughnutshape indicate a bug of mine?
Or is the ellipsoid shape just an approximization to keep inertia representable by few numbers?

I expect the latter, but if i'm wrong it might be the source of all my problems...

Code: Select all

inline float NonuniformInertiaBodiesGetGlobalInertia (const sVec3 &axis, const sVec3 &com, Body **bodies, const int numBodies) // unproofed
{
   float inertia = 0;
      
   for (int i=0; i<numBodies; i++)
   {
      sVec3 bcom; BodyGetCenterOfMass (bodies[i], bcom);
      float bmass, Ixx, Iyy, Izz;
      BodyGetMassMatrix (bodies[i], bmass, Ixx, Iyy, Izz);
      sVec3 d = bcom - com;   
      inertia += bmass * sVec3(d - axis * axis.Dot (d)).SqL(); // sql = squared length
      
      sMat4 matrix; BodyGetMatrix (bodies[i], matrix);
      sVec3 localAxis = matrix.Unrotate (axis);
      
      localAxis[0] /= Ixx;
      localAxis[1] /= Iyy;
      localAxis[2] /= Izz;

      inertia += 1.0f / localAxis.Length();
   }
   return inertia;
}


[Julio Jerez]

you can not just add inertias, you need to apply the parallel axis theorem.
there is a good explanation in Wikipedia, but if it get complicated you can check out function.

dgMatrix dgCollisionConvex::CalculateInertiaAndCenterOfMass (const dgMatrix& m_alignMatrix, const dgVector& localScale, const dgMatrix& matrix) const

It have a much simpler implementationtioan of the generalized central axis theorem.
The equation in matrix form is expressed in the comment, the Wikipedia version is a lithe confusing

basically the idea is that for you to add the inertias you need to calculate the inertia of each solid part relative to the same point. Since you have many cubes, each one with inertia relative its origin, you need to iterate over each cube calculating the inertia of each cube relative to the center of all cubes, that's where you apply the central or paralleled theorem.
then you can add all the inertias.
in your case since you know the cross inertia will be zero, you can easy do a test.

given the inertia of a box relative to its own origin, calculate the inertia relative to the center of all boxes, that is give by expression. (form the parallel axis theorem)

I = Ixx + m * (P - O) ^2)

p is the position of the box and O is the center of all cubes
as you can see, adding the inertias has tow parts, the local inertial and a compensation factor.
you seem to have the compensation but not the local part Ixx.
for cross inertia is more complex, but is just more algebra, look at the derivation in the newton function.

[JoeJ]

This drives me really crazy. I've read about parallel axis theorem again but i cant't find a bug in my code.

I = Ixx + m * (P - O) ^2)

To simplify, let's talk about point masses who do not have inertia (So there should be no local Ixx to bother, right?).

This still produces the doughnut shape, so the bug should be there, but it's the same formula i see everywhere:

Code: Select all

inline float Inertia_From_PointMassBodies (const sVec3 &axis, const sVec3 &allBodiesCom, Body **bodies, const int numBodies)
{
   float inertia = 0;
      
   for (int i=0; i<numBodies; i++)
   {
      sVec3 bodyCom = BodyGetGlobalCenterOfMass (bodies[i]);
      float bodyMass = BodyGetMass (bodies[i]);

      sVec3 diff = bodyCom - allBodiesCom; // allBodiesCom is precomputed   
      diff -= axis * axis.Dot(diff); // we want sqdist to axis through allBodiesCom
      
      float squareDistToAxis = diff.Dot(diff);
      inertia += bodyMass * squareDistToAxis;
   }
   return inertia;
}


Please look again, i made the code more readable - maybe you can spot something.
But it should be right, same as here:
I'm still not convinced the doughnut is wrong.

---

dgCollisionConvex::CalculateInertiaAndCenterOfMass

That's not the same i think. I need to get Inertia for a set of bodies, not from a convex shape.
To use Newton to verify, it should work to make a compound from the bodies at many rotations.
Tried this at runtime but no luck so far - maybe threading issues - need to do it at setup...

[Julio Jerez]

but why do you do this

Code: Select all

 float squareDistToAxis = diff.Dot(diff);
 inertia += bodyMass * squareDistToAxis;

that not what the formula apply to each component independent
say you have appoint mass,

the inertia of the point is given by

Ixx = mass * sum ((Pix - Po.x) ^ 2)
Iyy = mass * sum ((Piy - Po.y) ^ 2)
Izz = mass * sum ((Piz - Po.z) ^ 2)

it is the component of each point not the dot product, when you project the vector alone a direction that does no work, because the inertia is only valid for axis. you loop will only work when all matrices of all bodies are identity.

if you are dong this calculate the center of a group of bodies will must use the loop of the function newton, I can write a loop for you if you want.

and yes the function if force compound no for bodies, I menat for you to look athe code no to call call the function, the look is the same whet it is for bodies or chanpes.



vector float squareDistToAxis = diff.Dot(diff);
inertia += bodyMass * squareDistToAxis;

[JoeJ]

doughnut proofed.JPG (66.75 KiB) Viewed 6752 times

I just made the test using Newton Compound - the doughnut is in deed correct.
I rotated the bodies, but made compounds always using identity orientation,
then caching the x axis times Ixx for display - same doughnut shape.

So inertia matrix seems a bad approx for the real thing.
I always expected this to be the case but i'm surprised it has such a large effect.

I don't know how much using cross inertia could improve this.
also i don't request to use a more accurate inertia representation, most probably stability would decrease for a small gain in realism.

but why do you do this

Code: Select all
float squareDistToAxis = diff.Dot(diff);
inertia += bodyMass * squareDistToAxis;



that not what the formula apply to each component independent
say you have appoint mass,

the inertia of the point is given by

Ixx = mass * sum ((Pix - Po.x) ^ 2)
Iyy = mass * sum ((Piy - Po.y) ^ 2)
Izz = mass * sum ((Piz - Po.z) ^ 2)

I do this because my axis of interest is no principial axis like you are used to.

If my axis would be (1,0,0), then it's Ixx = mass * sum ((Pix - Po.x) ^ 2) like you say, but that's an opimized special case.

E. g. If axis is (0.71, 0.71, 0), then we need:
InertiaAboutAxisThroughPo = sum (mass[i] * sqdist(P[i], axisThroughPo))


The reason i need to know inertia over an arbitary axis is this:
The ankle joint inverted pendulum model may map one body to the base (foot),
but 19 bodies to the upper pendulum.
The inertia of the pendulum will be very non-uniform, so i do not approximate it with a point mass (like most papers suggest).
Also, there is no known best fit orientation for a precomputed Inertia matrix.

That's why i always loop over all bodies to get accurate inertia for a given axis.
Later i'll try to se what happens when i reduce accurary for speed. But first it needs to work...

if you are dong this calculate the center of a group of bodies will must use the loop of the function newton, I can write a loop for you if you want.

You may need to do so anyways for the same reasons.
For now i assume i'm correct but i'd trust you more than me

[JoeJ]

Edit:

Notice the white rays show only the cross section of the doughnut, like i've pointed in blue on the first posts picture.

Here the 'proof test' code:

Code: Select all

#define INERTIA_TEST_NUM_SAMPLES 64
   sVec3 inertiaVis[INERTIA_TEST_NUM_SAMPLES];

   void VisInertiaTest ()
   {
      for (int i=0; i<INERTIA_TEST_NUM_SAMPLES; i++)
         RenderVector (sVec3 (50, 10, 6), inertiaVis[i], 1,1,1);
   }

   void CreateInertiaTest ()
   {
      for (int a=0; a<INERTIA_TEST_NUM_SAMPLES; a++)
      {
         float angle = float (a) / float(INERTIA_TEST_NUM_SAMPLES) * PI*2.0;

         NewtonCollision *compoundShape = NewtonCreateCompoundCollision (world->world, 10000+a);
         NewtonCompoundCollisionBeginAddRemove (compoundShape);

         sVec3 origin (20 + 2.2 * (a&7), 10 + 2.2 * (a>>3), 6);
         sVec3 com (0,0,0);
         for (int i=0; i<MP_NUM_BODIES; i++)
         {
            sVec3 offset = sVec3 (0, mpd.baseDim[1] * 0.5 + mpd.segDist * i, 0);
            sMat4 matrix; matrix.Rotation()->FromEulerZYX (sVec3 (0,0,angle));
            inertiaVis[a] = matrix[0];
            matrix[3] = origin + matrix.Rotate(offset);
            matrix.Reset4thColumn();
            
            sVec3 dim = mpd.boxDim;
            Shape *boxShape = world->CreateBoxShape (dim);
            NewtonCollisionSetMatrix (boxShape, (float*)&matrix);
            NewtonCompoundCollisionAddSubCollision (compoundShape, boxShape);
            world->ReleaseShape (boxShape);

            com += matrix[3];
         }
         NewtonCompoundCollisionEndAddRemove (compoundShape);

         com /= float(MP_NUM_BODIES);
         sMat4 matrix; //
         matrix.Identity();
         matrix.Reset4thColumn();
         matrix[3] = com;

         Body *tempBody = world->CreateRigidBody (mpd.boxMass * MP_NUM_BODIES, matrix, compoundShape);
         
         float mass, Ixx, Iyy, Izz;
         BodyGetMassMatrix (tempBody, mass,  Ixx, Iyy, Izz);
         base_debug::logA->Print ("Ixx = ", Ixx);
         NewtonDestroyCollision (compoundShape);
         
         inertiaVis[a] *= Ixx;
      }
   }


[Julio Jerez]

I do this because my axis of interest is no principal axis like you are used to.

yes and that is part the problem.

I don't know how much using cross inertia could improve this.

Imagine a body rotate 90 degree, in the case the x become y, and so on.
for a body with two axis of symmetry the cross inertia add to zero, however when you brake a body into parts, the individual part has cross contributions, It is the sum of the contributions that add to zero not the parts. also the angular orientation add thing that are counter intuitive.
when doing this things the best is to abandoned the intuition and trust the Algebra.
usually when the algebra is right, all make sense but it is not until is all right and not before.

I am going to write a utility function that will do just what you want.
I need it for two reason, the multi body vehicle and pending feature that will be a balancing played controller similar to what you are doing.

I need it immediately because for multitier body vehicle, there is a place that need to calculate the tire lumped mass (the part for the vehicle the residing on top of the tire when then vehicle is at rest). I was using the weight of the chassis ignoring the rest, but the engine is a significant part of the vehicle weight so it need to be considered that too, so I may as well add the functionality.
Later I will need for the same reasons you need it now, the player pendulum model.

[JoeJ]

Ok, i've high hope on your function...