[ODE] Bug? Different behavior with different mass representations

[Matthias Baas]

Hi,

I have a problem with bodies that have mass properties where the center 
of gravity is not the origin and where the inertia tensor is not diagonal.

Below I have a small sample program that shows the problem. It's a 
simulation of a cube of 1 kg that lies between [0,1] in the x-, y- and 
z-direction. A force is applied at the global position (0.5, 0.5, 0) 
(=the center of one face) in z direction, so the cube should just start 
moving along the Z axis without spinning (there's no gravity or other 
external forces).
In the sample program you can switch between two different 
representations of the above scenario (by changing the "if (1)" to "if 
(0)").
In the first representation the origin of the local coordinate system of 
the cube is at the center of gravity and the cube is moved at position 
(0.5, 0.5, 0.5). In this case the center of gravity is (0,0,0) (in the 
body frame) and the inertia tensor is diagonal.
In the second representation the origin of the local coordinate system 
lies at a corner of the cube, so that the center of gravity is at (0.5, 
0.5, 0.5) and the inertia tensor is not diagonal anymore.

As I see it, both representations should result in the same simulation 
(barring rounding errors). However, that's not the case. The first 
representation shows the expected behavior where the rotation of the 
cube remains the identity. But in the second case, the cube starts spinning.

So can anybody explain where my mistake is? Or is it really a bug in 
ODE? By the way, I'm using the 0.5 binary release for Windows (it 
doesn't matter if it's the "single-trimesh" or the "double-no-trimesh" 
version). And it also doesn't matter if I'm using dWorldStep() or 
dWorldQuickStep()).

- Matthias -


// ODE bug?

#include <stdio.h>
#include "ode/ode.h"

int main()
{
   dWorldID world;
   dBodyID body;
   dMass mass;
   const dReal* pos;
   const dReal* rot;
   int i;

   // Create a simulation with a cube of 1 kg that lies in the
   // range [0,1] on all axes. At the beginning an initial force
   // is applied at point (0.5, 0.5, 0) which is the center of one
   // face of the cube. The direction of the force is in negative
   // normal direction. As a result the cube should just move in
   // Z direction without spinning.
   //
   // You can choose between two different representations of this
   // scenario which leads to two different results.

   world = dWorldCreate();
   body = dBodyCreate(world);

   // Set the mass parameters and the initial position of the cube.
   // Here you can select which representation should be used.
   if (1)
   {
     // Representation 1:
     // Set the mass parameters of a cube that's defined between
     // [-0.5, 0.5], so the center of gravity is identical with
     // the origin and the inertia tensor is diagonal. The cube is
     // then moved to position (0.5, 0.5, 0.5) so that it lies in
     // the range [0,1].
     dMassSetParameters(&mass, 1.0, 0,0,0, 0.5/3, 0.5/3, 0.5/3, 0,0,0);
     dBodySetPosition(body, 0.5, 0.5, 0.5);
   }
   else
   {
     // Representation 2:
     // Set the mass parameters of a cube that's defined between [0,1].
     // The center of gravity at (0.5, 0.5, 0.5) is now different
     // from the origin and the inertia tensor is not diagonal anymore.
     // The cube's position now remains at the origin.
     dMassSetParameters(&mass, 1.0, 0.5,0.5,0.5, 2.0/3, 2.0/3, 2.0/3, 
-0.25,-0.25,-0.25);
     dBodySetPosition(body, 0.0, 0.0, 0.0);
   }

   dBodySetMass(body, &mass);

   // Apply a force at a position given in the global coordinate system
   dBodyAddForceAtPos(body, 0,0,10, 0.5,0.5,0);

   // Do a couple of simulation steps and print the position and
   // orientation of the cube...
   // Using representation 1 the rotation matrix remains the identity
   // which is the expected result. However, using representation 2
   // the cube starts spinning.
   for(i=0; i<20; i+=1)
   {
     pos = dBodyGetPosition(body);
     rot = dBodyGetRotation(body);
     printf("Step %2d: pos: %f, %f, %f\n", i, pos[0], pos[1], pos[2]);
     printf("  %f  %f  %f\n", rot[0], rot[1], rot[2]);
     printf("  %f  %f  %f\n", rot[4], rot[5], rot[6]);
     printf("  %f  %f  %f\n", rot[8], rot[9], rot[10]);

     dWorldStep(world, 0.04);
     //dWorldQuickStep(world, 0.04);
   }

   dBodyDestroy(body);
   dWorldDestroy(world);
   return 0;
}