Compound collision and center of mass

[PJani]

Hey.

I have one problem which i don't know how to solve, I am building compound collision from smaller collisions which can have different mass/density than others. I need to calculate center of mass and corresponding inertia for this compound collision.

Does somebody have idea how to do such thing.

[Julio Jerez]

are you doing in the engine or it is something you like to do for something else?
In engine, core 300, handle it automatically, even with scale, here is the code

Code: Select all

void dgCollisionCompound::MassProperties ()
{
   dgVector origin (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));
   dgVector inertia (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));
   dgVector crossInertia (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));

   dgPolyhedraMassProperties localData;
   DebugCollision (dgGetIdentityMatrix(), CalculateInertia, &localData);
   dgFloat32 volume = localData.MassProperties (origin, inertia, crossInertia);
   _ASSERTE (volume > dgFloat32 (0.0f));

   dgFloat32 invVolume = dgFloat32 (1.0f)/volume;
   m_centerOfMass = origin.Scale (invVolume);
   m_centerOfMass.m_w = volume;
   m_inertia = inertia.Scale (invVolume);
   m_crossInertia = crossInertia.Scale(invVolume);
}

if you want to do this for something you are trying you can take the class dgPolyhedraMassProperties and use it.
you can also use the central axis theorem if each shape already have the inertia relative to its center, this is the method that Collision instace use.
This is very fast and can be use in real time each tiem a shape is scale of reposition.

Code: Select all

dgMatrix dgCollisionConvex::CalculateInertiaAndCenterOfMass (const dgMatrix& matrix, const dgVector& scale) const
{
   // using general central theorem, is much faster and more accurate;
   //IImatrix = IIorigin + mass * (displacemnet % displacemnet) * identityMatrix - transpose(displacement) * displacement );
   
   dgVector scaleII (scale.CompProduct(scale));
   dgVector scaleIJ (scale[0] * scale[1], scale[0] * scale[2], scale[1] * scale[2], dgFloat32 (0.0f));

   dgMatrix inertia (dgGetIdentityMatrix());
   inertia[0][0] = m_inertia[0] * scaleII[0];
   inertia[1][1] = m_inertia[1] * scaleII[1];
   inertia[2][2] = m_inertia[2] * scaleII[2];
   inertia[0][1] = m_crossInertia[2] * scaleIJ[0];
   inertia[1][0] = m_crossInertia[2] * scaleIJ[0];
   inertia[0][2] = m_crossInertia[1] * scaleIJ[1];
   inertia[2][0] = m_crossInertia[1] * scaleIJ[1];
   inertia[1][2] = m_crossInertia[0] * scaleIJ[2];
   inertia[2][1] = m_crossInertia[0] * scaleIJ[2];
   inertia = matrix.Inverse() * inertia * matrix;
   
   dgVector origin (m_centerOfMass + matrix.m_posit);
   origin = origin.CompProduct(scale);
   dgFloat32 mag = origin % origin;

   //dgFloat32 unitMass = m_centerOfMass.m_w;
   dgFloat32 unitMass = dgFloat32 (1.0f);
   for (dgInt32 i = 0; i < 3; i ++) {
      inertia[i][i] += unitMass * (mag - origin[i] * origin[i]);
      for (dgInt32 j = i + 1; j < 3; j ++) {
         dgFloat32 crossIJ = - unitMass * origin[i] * origin[j];
         inertia[i][j] += crossIJ;
         inertia[j][i] += crossIJ;
      }
   }

   inertia.m_posit = origin;
   inertia.m_posit.m_w = 1.0f;
   return inertia;
}
 

[Julio Jerez]

The above code assumes uniform density, on a compound, by you can do a very simple trick. to get it do the inertia for different density.
Basically the trick is to prescale the sub shapes to match and shape of a uniform density.
This is how. Say the compound only has one box. Of dimension (X x Y x Z)
The inertia of any object of uniform mass is given by
Sum (Mi * xi * xi)
Since the mass if uniform you can actually factor the mass from the equation and use the sumation of each point particle. thist is a practice used in mechanical and civil
engineering when for specifying the shape of Beam and othe premaded construction materials. For cross sections it is called Area of Inertia, for a volume is called Volume of Inertia.
This simplifies dealing with building material of the same shape but of different marterials.
a typical exercise in mechanic 101 is to calculate the size a Wood Beam would have to be to have the strength of a give steel steal bean of the same shape.
The answerer is you multiple the steel beam cross section size by the square root of the ratio of densities between wood and steel.
You can do the same for a compound collision, here are the steps:

-make your compound collision
-save the scale of each shape.
-take the mass of any of the shape in the compound (the first is just as s good)
-iterate over the rest of all shapes
-For each shape set the scale to be ShapeScale = sqrt (shapeMass / referenceMass)
-calculate the inertia of the scaled compound.
-Multiply the inertia o feth scaled compound by (Numbershapes * referenceMass) his is the total inertia
-Interate again restoring the original size of each sub shape

[PJani]

Thank you julio that the thing i was looking for.

[PJani]

One quick question. Newton uses interia based on center of shape or center of mass?

[Julio Jerez]

from center of mass of the shape. this is why I ussually set the com frpm the calculated inertia of the shape.

Each shape contatin the intertai in teh fore of product of integration that are calculated using teh Green and stok theorem. Then form there is is very fast an eassy to calculet teh intetian using the central thorem.

[PJani]

I have another quick question relating center of mass...

where is being force applied using NewtonBodyAddForce center of mass or center of shape? Because lately i am being so confused because i am working so much with inertia and center of mass

[Julio Jerez]

all forces are applyed at the center of mass of the body, in global space

if you want to apply a force at a point you need to read the global center of mass form the body, and subtract the point form the center of mass.
there is and FAQ for that.

in many case many force call back do not do thes because the origin of teh matrix coicide with eth centr of mass, but any applycation spudl do it right

somethog like this

Code: Select all

   // calculate toruqe generated by a force off centert of mass
   NewtonBodyGetMatrix(body, &matrix[0][0]);
   NewtonBodyGetCentreOfMass (body, &com[0]);
   dVector point (pointOnBodyInGlobalSpace - matrix.RotateVector ( com));
   dVector torque (point * force); // * is a cross product here

   NewtonBodyAddForce (body, &force.m_x);
   NewtonBodyAddTorque (body, &torque.m_x);

[PJani]

I have compound collision with sub collisions each can have different density...

Now...i have problem with inertia and center of mass...

Each time i add sub collision i do(if i am removing i do oposite of this):
where m_mass is current mass of compound, m_center_of_mass is current c.o.m of compound, ci.mass is added sub collision mass, ci.center_of_mass is added sub collision center of mass

Code: Select all

      float mass_sum_inv = 1.0/(m_mass + ci.mass);

      float com_x = (m_center_of_mass.x * m_mass + ci.mass * ci.center_of_mass.x) * mass_sum_inv;
      float com_y = (m_center_of_mass.y * m_mass + ci.mass * ci.center_of_mass.y) * mass_sum_inv;
      float com_z = (m_center_of_mass.z * m_mass + ci.mass * ci.center_of_mass.z) * mass_sum_inv;


I am following this formula for center of mass(sory for tex):
system_{c.o.m.} = (SUM_{i=0}^{N}{objcom_i * objmass_i}) / system_{mass}

Now my result for having all sub collisions same density is(i dont know if this errors(at x and y) are negligible):

Code: Select all

NEWTON_COM: Vector3(1.66308e-008, 7.06808e-009, 0.0285715)
    MY_COM: Vector3(1.06924e-007, 4.23631e-009, 0.0285715)


Now about inertia is even possible to do such thing to calculate inertia at the fly? With minimal iteration over sub collisions at the end? Because i expect lots of adding and removing of sub collisions

[Julio Jerez]

the function to calculate the inertia of a compound colldion is not tha simple
beacisally you nee to use the integrall of inerta no the inerta of each shape.
if you compound do no have too make shape you cna use the newton funtion for that, whi resolve to this.

Code: Select all

void dgCollisionCompound::MassProperties ()
{
/*   
   dgVector totalOrigin (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));
   dgVector totalInertia(dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));
   dgVector totalCrossInertia (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));

   dgFloat32 totalVolume = dgFloat32 (0.0f);
   dgTree<dgNodeBase*, dgInt32>::Iterator iter (m_array);

   for (iter.Begin(); iter; iter ++) {
      dgCollisionInstance* const collisionInstance = iter.GetNode()->GetInfo()->GetShape();
      dgCollision* const collision = collisionInstance->GetChildShape();

      dgFloat32 volume = collision->m_centerOfMass.m_w;
      totalVolume += volume;
      dgMatrix inertia (collision->CalculateInertiaAndCenterOfMass (collisionInstance->GetLocalMatrix(), collisionInstance->m_scale));
      totalOrigin += inertia.m_posit.Scale (volume);;
      totalInertia += dgVector (inertia[0][0], inertia[1][1], inertia[2][2], dgFloat32 (0.0f));
      totalCrossInertia += dgVector (inertia[1][2], inertia[0][2], inertia[0][1], dgFloat32 (0.0f));
   }

   if (totalVolume > dgFloat32 (0.0)) {
      m_inertia = totalInertia;
      m_crossInertia = totalCrossInertia;
      m_centerOfMass = totalOrigin.Scale (dgFloat32 (1.0f)/totalVolume);
      m_centerOfMass.m_w = totalVolume;
   }
*/

   dgVector origin (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));
   dgVector inertia (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));
   dgVector crossInertia (dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f), dgFloat32 (0.0f));

   dgPolyhedraMassProperties localData;
   DebugCollision (dgGetIdentityMatrix(), CalculateInertia, &localData);
   dgFloat32 volume = localData.MassProperties (origin, inertia, crossInertia);
   _ASSERTE (volume > dgFloat32 (0.0f));

   dgFloat32 invVolume = dgFloat32 (1.0f)/volume;
   m_centerOfMass = origin.Scale (invVolume);
   m_centerOfMass.m_w = volume;
   m_inertia = inertia.Scale (invVolume);
   m_crossInertia = crossInertia.Scale(invVolume);
}

the commented out code does the calculation in linea time.

Now about inertia is even possible to do such thing to calculate inertia at the fly? With minimal iteration over sub collisions at the end? Because i expect lots of adding and removing of sub collisions

yes It is possible to calculate the inertia opf a compound collision in constant time by keeping the sumations of the chidren withe the compound shape, and updating each time a shape is added, removed, scaled, change the matrix.
whowever I decided to to do that because now the shapes can move and that will requered lots of complexity to teh inetface.

I think that if you just call the functions Newton calculate Inertia, it will be fine if your compound are relativally small.
if it is not then the comented out portion of the code does the same using the precalculated values of the children. and the is one order of magnitud faster.
try those first.

[Julio Jerez]

also notice that function void dgCollisionCompound::EndAddRemove () recompute the inertias of children collision, so if you are adding and removing collisions the inertia are always maintained.
when you call NewtponCalculateIntert that funsion rund in constant time for all shapes, this is what it does for a compound.

Code: Select all

dgMatrix dgCollisionCompound::CalculateInertiaAndCenterOfMass (const dgMatrix& matrix, const dgVector& scale) const
{
   dgVector origin (m_centerOfMass + matrix.m_posit);
   dgFloat32 mag = origin % origin;

   dgMatrix inertia (dgGetIdentityMatrix());
   inertia[0][0] = m_inertia[0];
   inertia[1][1] = m_inertia[1];
   inertia[2][2] = m_inertia[2];
   inertia[0][1] = m_crossInertia[2];
   inertia[1][0] = m_crossInertia[2];
   inertia[0][2] = m_crossInertia[1];
   inertia[2][0] = m_crossInertia[1];
   inertia[1][2] = m_crossInertia[0];
   inertia[2][1] = m_crossInertia[0];
   inertia = matrix.Inverse() * inertia * matrix;

   //dgFloat32 volume = m_centerOfMass.m_w;
   dgFloat32 unitMass = dgFloat32 (1.0f);
   for (dgInt32 i = 0; i < 3; i ++) {
      inertia[i][i] += unitMass * (mag - origin[i] * origin[i]);
      for (dgInt32 j = i + 1; j < 3; j ++) {
         dgFloat32 crossIJ = - unitMass * origin[i] * origin[j];
         inertia[i][j] += crossIJ;
         inertia[j][i] += crossIJ;
      }
   }

   inertia.m_posit = origin;
   inertia.m_posit.m_w = 1.0f;
   return inertia;
}

all this discution only apply if EndAddRemove is used each time a shape is modified. you can batch a bunch of update and them call the update.
I do not think you have any problems. It should all working as it is now.

[PJani]

Hmm i didn't know that 0_0 only thing i am trying to achive is different density for shapes. That why i am trying to compute com and interia at add

EDIT: i do lots of adds removes and then i call end

[Julio Jerez]

if you do may add remove, then you sodul be fine. I will actaully be better that teh amortized solution. In fatc even a reltive small number of add remove will do.

[PJani]

do you know what would be realy great if addsubcollisioncompund would have mass scale or density scale parameter so the newton calculate interia matrix would take different mass/density into account.

[PJani]

ShapeScale = sqrt (shapeMass / referenceMass) why is here sqrt?