[ODE] friction anisotropy solution
Graham Fyffe
gfyffe at gmail.com
Mon May 2 21:38:41 MST 2005
Hi again, everyone. I'm not sure if this got lost in the shuffle,
because my original post was too long. So here it is again:
I think I have a solution to the frictiin anisotropy problem. It
works well in all my tests. Here are my main comments, and then the
whole code follows:
If we just let ODE compute the friction directions, we get anisotropy
(due to the friction pyramid axes), which results in objects
"catching" along the coordinate axes (X, Y, Z). This is especially
noticible for sliding objects, as an object sliding just off the X
direction will change direction until it is sliding on X, etc. This
is a Bad Thing.
My solution is to align the friction direction with the body's tangent
velocity, in the hopes of minimizing the problem.
Further notes: I'm actually aligning the friction direction to be
perpendicular to both the contact normal and the relative tangent
velocity between the two colliding bodies, measured at the collision
position. Then, ODE will calculate the second friction direction
perpendicular to mine and to the collision normal, guaranteeing that
the two friction directions are orthogonal to the normal, with one of
them aligned as much as possible to the relative tangent velocity.
The main benefit of this method is that the friction errors are
restricted to the cases where there are external forces acting on the
bodies. In this case, the solution can only be as bad as the original
friction pyramid. Lacking significant external forces, a sliding /
spinning body will behave as expected, with nice, isotropic friction
behavior. Boxes sliding in a straight line will continue to do so
until they come to rest, regardless of the direction or the
orientation :)
One final note: I have not covered here how to deal with the case
where you actually WANT anisotropic friction. Basically, what you
should do is calculate the desired friction mu that exists in the
direction that is computed for fdir1 (e.g. using a friction ellipse)
and also a nice mu2 in the perpendicular direction. Same goes if you
want different values for slip1 and slip2. Basically just sample them
on the directions that you get.
The code follows. Take special note of 'nearCallback'. Cheers!
- Graham
// header comments stripped for message brevity //
#include <ode/ode.h>
#include <drawstuff/drawstuff.h>
#ifdef _MSC_VER
#pragma warning(disable:4244 4305) // for VC++, no precision loss complaints
#endif
// select correct drawing functions
#ifdef dDOUBLE
#define dsDrawBox dsDrawBoxD
#define dsDrawSphere dsDrawSphereD
#define dsDrawCylinder dsDrawCylinderD
#define dsDrawCappedCylinder dsDrawCappedCylinderD
#endif
// some constants
#define LENGTH 0.2 // box length & width
#define HEIGHT 0.05 // box height
#define MASS 0.2 // mass of box[i][j] = (i+1) * MASS
#define FORCE 0.05 // force applied to box[i][j] = (j+1) * FORCE
#define MU 0.5 // the global mu to use
#define GRAVITY 0.5 // the global gravity to use
#define N1 10 // number of different forces to try
#define N2 10 // number of different masses to try
// dynamics and collision objects
static dWorldID world;
static dSpaceID space;
static dBodyID body[N1][N2];
static dJointGroupID contactgroup;
static dGeomID ground;
static dGeomID box[N1][N2];
static void nearCallback (void *data, dGeomID o1, dGeomID o2)
{
int i;
// only collide things with the ground
// [NOTE] even though we do this here, the friction stuff does
handle the case of 2 moving objects.
bool g1 = (o1 == ground);
bool g2 = (o2 == ground);
if (g1 == g2) return;
dBodyID b1 = dGeomGetBody(o1);
dBodyID b2 = dGeomGetBody(o2);
if (g1)
{
dBodyID t = b1;
b1 = b2;
b2 = t;
}
dContact contact[3]; // up to 3 contacts per box
for (i=0; i<3; i++) {
contact[i].surface.mode = dContactSoftCFM | dContactApprox1;
contact[i].surface.mu = MU;
contact[i].surface.soft_cfm = 0.01;
// If we just let ODE compute the friction directions, we get anisotropy,
// which results in objects "catching" along the coordinate
axes (X, Y, Z).
// This is especially noticible for sliding objects, as an object sliding
// just off the X direction will change direction until it is
sliding on X, etc.
// This is a Bad Thing.
// This first test just uses the body's local X:
//dBodyVectorToWorld(b1, 1, 0, 0, contact[i].fdir1);
// Unfortunately this just changes the axes of the problem, and
the problem
// itself is still present in full force.
// This next test tries to align the friction direction with
the body's tangent velocity,
// in the hopes of minimizing the problem:
dBodyGetPointVel(b1,
contact[i].geom.pos[0], contact[i].geom.pos[1],
contact[i].geom.pos[2],
contact[i].fdir1);
dVector3 fdir1_b2;
if (b2)
{
dBodyGetPointVel(b2,
contact[i].geom.pos[0], contact[i].geom.pos[1],
contact[i].geom.pos[2],
fdir1_b2);
contact[i].fdir1[0] -= fdir1_b2[0];
contact[i].fdir1[1] -= fdir1_b2[1];
contact[i].fdir1[2] -= fdir1_b2[2];
}
// at this point, contact[i].fdir1 is the relative tangent
velocity of the two bodies.
dCROSS(contact[i].fdir1, =, contact[i].fdir1, contact[i].geom.normal);
// now, contact[i].fdir1 is perpendicular to both the normal and
// the relative tangent velocity.
double length = sqrt(contact[i].fdir1[0] * contact[i].fdir1[0]
+ contact[i].fdir1[1] * contact[i].fdir1[1]
+ contact[i].fdir1[2] * contact[i].fdir1[2]);
if (length > 1e-12)
{
// we only use our calculated direction if it has
enough precision
contact[i].fdir1[0] /= length;
contact[i].fdir1[1] /= length;
contact[i].fdir1[2] /= length;
dNormalize3(contact[i].fdir1);
contact[i].surface.mode |= dContactFDir1;
}
// This system performs very well, on sliding box tests and on
// sliding, spinning box tests!
}
if (int numc = dCollide (o1,o2,3,&contact[0].geom,sizeof(dContact))) {
for (i=0; i<numc; i++) {
dJointID c = dJointCreateContact (world,contactgroup,contact+i);
dJointAttach (c,b1,b2);
}
}
}
static void start()
{
static float xyz[3] = {1.7772,-0.7924,2.7600};
static float hpr[3] = {90.0000,-54.0000,0.0000};
dsSetViewpoint (xyz,hpr);
}
static void simLoop (int pause)
{
int i;
if (!pause) {
int const maxIt = 10;
for (int t = 0; t < maxIt; t ++)
{
dSpaceCollide (space,0,&nearCallback);
dWorldStep (world,0.05 / maxIt);
dJointGroupEmpty (contactgroup);
}
}
dsSetColor (1,0,1);
dReal sides[3] = {LENGTH,LENGTH,HEIGHT};
for (i=0; i<N1; i++) {
for (int j=0; j<N2; j++) {
dsDrawBox (dGeomGetPosition(box[i][j]),dGeomGetRotation(box[i][j]),
sides);
}
}
}
static void reset()
{
int i,j;
// bodies
for (i=0; i<N1; i++)
{
for (j=0; j<N2; j++)
{
dBodySetPosition(body[i][j], i*2*LENGTH,
j*2*LENGTH, HEIGHT*0.5);
// this first scenerio tests sliding while spinning
dMatrix3 R;
dRFromAxisAndAngle(R, 0, 0, 1, 0);
dBodySetRotation(body[i][j], R);
dBodySetAngularVel(body[i][j], 0, 0, i * 2);
dBodySetLinearVel(body[i][j], 0, j * 0.1, 0);
// this next one tests friction at various angles
// dMatrix3 R;
// dRFromAxisAndAngle(R, 0, 0, 1, i / (N1 - 1.0)
* 3.14159 / 2.0);
// dBodySetRotation(body[i][j], R);
// dBodySetLinearVel(body[i][j], 0, j * 0.5, 0);
// dBodySetAngularVel(body[i][j], 0, 0, 0);
}
}
}
static void command (int cmd)
{
switch (cmd)
{
case 'r':
reset();
break;
}
}
int main (int argc, char **argv)
{
int i,j;
dMass m;
// setup pointers to drawstuff callback functions
dsFunctions fn;
fn.version = DS_VERSION;
fn.start = &start;
fn.step = &simLoop;
fn.command = command;
fn.stop = 0;
fn.path_to_textures = "../../drawstuff/textures";
// create world
world = dWorldCreate();
space = dHashSpaceCreate (0);
contactgroup = dJointGroupCreate (0);
dWorldSetGravity (world,0,0,-GRAVITY);
ground = dCreatePlane (space,0,0,1,0);
// bodies
for (i=0; i<N1; i++) {
for (j=0; j<N2; j++) {
body[i][j] = dBodyCreate (world);
dMassSetBox (&m,1,LENGTH,LENGTH,HEIGHT);
//dMassAdjust (&m,MASS*(j+1));
dBodySetMass (body[i][j],&m);
box[i][j] = dCreateBox (space,LENGTH,LENGTH,HEIGHT);
dGeomSetBody (box[i][j],body[i][j]);
}
}
reset();
// run simulation
dsSimulationLoop (argc,argv,352,288,&fn);
dJointGroupDestroy (contactgroup);
dSpaceDestroy (space);
dWorldDestroy (world);
return 0;
}
More information about the ODE
mailing list