[ODE] rolling friction

[Nagymathe Denes]

Hi,

Using angular damping may look correct in most situations, but physically it 
isn't (at least you have to check whether the body is touching some surface 
or not to avoid the damping in mid-air). Some while ago i have implemented 
an estimation of rolling resistance by changing the contact point a bit 
further to the direction of the movement. It works pretty well. This code 
fraction comes into nearcallback, just before creating the contact joint.
I have tried two versions; their results are quite similar, but the second 
one is simpler.

#ifdef ROLLING_RESISTANCE
// The basic method for rolling resistance: shift contact point a bit 
forward to get some torque against rolling
// The normal should also be changed; we don't do that yet, but that seems 
unnecessary
// it could be limited to shapes with curved surfaces (sphere, capsule, 
cylinder)
// and/or the factor should be set in the geom or body
// Unsolved issue: spheres still keep their spin around the 'vertical' axis, 
althought they do stop rolling
{
dReal rrfactor = 0.01f;
int j;
//method a:
/* dReal bq[4], dq[4];
dVector3 p1,p2,p3;
dMatrix3 bR;
dWtoDQ (b1->avel,b1->q,dq);
for (j=0; j<4; j++) bq[j] = b1->q[j] + 0.5*dq[j]*rrfactor;
dNormalize4 (bq);
dQtoR (bq,bR);
for (j=0; j<3; j++) p1[j] = nc->geom.pos[j] - b1->pos[j];
dMULTIPLY1_331 (p2,b1->R,p1);
dMULTIPLY0_331 (p3,bR,p2);
for (j=0; j<3; j++) nc->geom.pos[j] -= p3[j] - p1[j];
if (b2)
{
dWtoDQ (b2->avel,b2->q,dq);
for (j=0; j<4; j++) bq[j] = b2->q[j] + 0.5*dq[j]*rrfactor;
dNormalize4 (bq);
dQtoR (bq,bR);
for (j=0; j<3; j++) p1[j] = nc->geom.pos[j] - b2->pos[j];
dMULTIPLY1_331 (p2,b2->R,p1);
dMULTIPLY0_331 (p3,bR,p2);
for (j=0; j<3; j++) nc->geom.pos[j] -= p3[j] - p1[j];
}
*/
//method b:
dVector3 p1,p2;
const dReal *bp = dBodyGetPosition(b1);
const dReal *av = dBodyGetAngularVel(b1);
for (j=0; j<3; j++) p1[j] = nc->geom.pos[j] - bp[j];
dCROSS(p2, =, av, p1);
for (j=0; j<3; j++) nc->geom.pos[j] -= p2[j]*rrfactor;
if (b2)
{
bp = dBodyGetPosition(b2);
for (j=0; j<3; j++) p1[j] = nc->geom.pos[j] - bp[j];
av = dBodyGetAngularVel(b2);
dCROSS(p2, =, av, p1);
for (j=0; j<3; j++) nc->geom.pos[j] -= p2[j]*rrfactor;
}
}
#endif //ROLLING_RESISTANCE