Jump to content

  • Log In with Google      Sign In   
  • Create Account


Calculating impulse due to rigid body collision with friction


Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

  • You cannot reply to this topic
16 replies to this topic

#1 BirdiePeep   Members   -  Reputation: 134

Like
0Likes
Like

Posted 22 September 2007 - 06:36 AM

I am programming a simple physics simulation and have been following along with some books about the subject. Everything has gone well so far, although at the moment I'm stuck on how to calculate the impulse when two rigid bodies collide. I'll explain the algorithms that I'm using below. In a frictionless environment I would use this equation to find the impulse that would be applied to both rigid bodies in the collision. J = -Vr(1+e) / {1/m1 + 1/m2 + n . [(r1 x n)/I1] x r1 + n . [(r2 x n)/I2] x r2} J - The scalar impulse Vr - The relative closing velocity of the two bodies e - The coefficient of restitution m1 - The mass of rigid body 1 m2 - The mass of rigid body 2 n - The normal of the collision r1 - The vector from the center of mass of rigid body 1 to the point of collision r2 - The vector from the center of mass of rigid body 2 to the point of collision I1 - Inertia tensor for rigid body 1 I2 - Inertia tensor for rigid body 2 I apply the impulse to the objects as such. v1 = v1 + (Jn)/m1 v2 = v2 + (-Jn)/m2 w1 = w1 + (r1 x Jn)/I1 w2 = w2 + (r2 x -Jn)/I2 v1 - Velocity of rigid body 1 v2 - Velocity of rigid body 2 w1 - Angular velocity of rigid body 1 w2 - Angular velocity of rigid body 2 I have tested this and everything works correctly. The thing that I'm clueless about is when I then decide to throw dynamic/kinetic friction into the equation. So far from the books I'm reading this is what I have been able to figure out. This is the updated equations to update the rigid bodies. v1 = v1 + [Jn + (uJ)t]/m1 v2 = v2 + [-Jn + (uJ)t]/m2 w1 = w1 + {r1 x [Jn + (uJ)t]}/I1 w2 = w2 + {r2 x [-Jn + (uJ)t]}/I2 u - The coefficient of friction t - The tangent normal, this is a unit vector that is perpendicular to the contact normal. The tangent normal is calculated with this equation. t = [(n x Vr) x n] t = t/|t| I can understand the concept of what his happening in the above equations. I also know that I'm going to need to change the impulse calculation equation because of friction, although I'm not sure how this is done. The books I have only skim over this aspect and I'm a little unsure how to continue. I would imagine it is not difficult, although my searching has not turned up any useful information. So my question is, when two rigid bodies collide in a frictional envionrment how does one go about calculating the impulse? Thank for reading. :p

Sponsor:

#2 xissburg   Members   -  Reputation: 200

Like
0Likes
Like

Posted 22 September 2007 - 07:06 AM

I've learned this looking at the source code of the MrRowl's physics simulators(thanks for that :) ) and that works very nicely. Here's what I do(this is from the rigid sphere dynamics simulator I did(it has no resting contact :( ):


//Solves a collision between 2 spheres
void PHYSICALWORLD::SolveCollision(const COLLISIONINFO& c)
{
SPHERE *sA, *sB;

//get pointers to the spheres so that sA points to the sphere with mass
//and sB points to the sphere with 0 mass or if both have 0 mass return
if(sphere[c.indexA]->GetMass() > 0.f)
{
sA = sphere[c.indexA];
sB = sphere[c.indexB];
}
else if(sphere[c.indexB]->GetMass() > 0.f)
{
sA = sphere[c.indexB];
sB = sphere[c.indexA];
}
else
return;

bool massB = sB->GetMass() > 0.f;

//put them at the 'timeOfImpact position' and also linearly interpolate the velocities
sA->SetPosition(sA->GetPrevPosition() + (sA->GetPosition() - sA->GetPrevPosition())*c.timeOfImpact, false);
sA->SetOrientation(sA->GetPrevOrientation() + (sA->GetOrientation() - sA->GetPrevOrientation())*c.timeOfImpact,
false, false);
sA->SetLinVelocity((sA->GetPosition() - sA->GetPrevPosition())/(sA->GetDeltaTime()*c.timeOfImpact));
sA->SetAngVelocity((sA->GetOrientation() - sA->GetPrevOrientation())/(sA->GetDeltaTime()*c.timeOfImpact));

if(massB)
{
sB->SetPosition(sB->GetPrevPosition() + (sB->GetPosition() - sB->GetPrevPosition())*c.timeOfImpact, false);
sB->SetOrientation(sB->GetPrevOrientation() + (sB->GetOrientation() - sB->GetPrevOrientation())*c.timeOfImpact,
false, false);

sB->SetLinVelocity((sB->GetPosition() - sB->GetPrevPosition())/(sB->GetDeltaTime()*c.timeOfImpact));
sB->SetAngVelocity((sB->GetOrientation() - sB->GetPrevOrientation())/(sB->GetDeltaTime()*c.timeOfImpact));
}

//compute collision normal
VECTOR3 n = (sA->GetPosition() - sB->GetPosition())/(sA->GetRadius() + sB->GetRadius());
VECTOR3 rA = -n*sA->GetRadius();
VECTOR3 rB = n*sB->GetRadius();

//relative velocity at contact point
VECTOR3 vRel = sA->GetLinVelocity() + (sA->GetAngVelocity() ^ rA);

if(massB)
vRel -= sB->GetLinVelocity() + (sB->GetAngVelocity() ^ rB);

//normal relative velocity
float nvRel = n*vRel;

//the bodies are going to penetrate
if(nvRel < 0.f)
{
//compute and apply normal impulse. there's no need to include
//the inertia tensor here because we're working with spheres only
float numerator = -(1.f+(sA->GetElasticity() + sB->GetElasticity())*0.5f)*nvRel;
float denominator = sA->GetInvMass() + sB->GetInvMass();
float j = numerator/denominator;
VECTOR3 jn = n*j;
sA->ApplyImpulse(jn, rA);
if(massB) sB->ApplyImpulse(-jn, rB);


//FRICTION IMPULSE
//compute the relative tangential velocity at the contact point
VECTOR3 tVel = vRel - n*nvRel;
//tangential speed
float tSpd = tVel.LengthSq();

//the tangential speed is greater than 0 then we must apply friction impulse
if(tSpd > 0.0000001f)
{
tSpd = sqrtf(tSpd);
//friction normal to use in the impulse equation. Use the
//minus sign so that it points in the opposite direction
//of the tangential velocity vector
VECTOR3 fn = -tVel/tSpd;
//compute the denominator of the impulse equation the
//same way you do when computing for normal impulse
//but you use the friction normal instead
//the inertia is included here
denominator = sA->GetInvMass() + fn*((sA->GetInvInertia()*(rA^fn))^rA);
if(massB) denominator += sB->GetInvMass() + fn*((sB->GetInvInertia()*(rB^fn))^rB);

//the friction impulse. the tangential speed is the numerator
float fj = tSpd/denominator;

//as when working with forces, if the friction impulse
//is greater than the normal impulse times the static
//friction coefficient, which means a slide, the
//friction impulse is the normal impulse times the
//dynamic friction coefficient else there's no change
//to do there
if(fj > j*(sA->GetStaticFriction() + sB->GetStaticFriction())*0.5f)
fj = j*(sA->GetDynamicFriction() + sB->GetDynamicFriction())*0.5f;

//compute the friction impulse vector
VECTOR3 fjv = fn * fj;

//apply the impulses that same way
sA->ApplyImpulse(fjv, rA);
if(massB) sB->ApplyImpulse(-fjv, rB);
}

sA->ApplyForce(gravity*sA->GetMass());
sA->ApplyDragForce(fluidDensity);
sA->Integrate(sA->GetDeltaTime()*(1.f - c.timeOfImpact));
sA->ZeroOutForces();

if(massB)
{
sB->ApplyForce(gravity*sB->GetMass());
sB->ApplyDragForce(fluidDensity);
sB->Integrate(sB->GetDeltaTime()*(1.f - c.timeOfImpact));
sB->ZeroOutForces();
}
}
}



Thats all, Hope this helps.

#3 Dirk Gregorius   Members   -  Reputation: 730

Like
0Likes
Like

Posted 22 September 2007 - 11:19 AM

IIRC a good derivation of a contact with friction can be found e.g. in the PhD thesis of B. Mirtich:

http://www.kuffner.org/james/software/dynamics/mirtich/index.html


The method in MrRowl's Jiglib comes IIRC from the Guendelman paper. His PhD and the original paper can be found here:

http://graphics.stanford.edu/~erang/papers/thesis/thesis.pdf
http://graphics.stanford.edu/papers/rigid_bodies-sig03/


Both good lectures...


HTH,
-Dirk

#4 chand81   Members   -  Reputation: 181

Like
0Likes
Like

Posted 22 September 2007 - 07:03 PM

Quote:
Original post by KiroNeem

J = -Vr(1+e) / {1/m1 + 1/m2 + n . [(r1 x n)/I1] x r1 + n . [(r2 x n)/I2] x r2}

J - The scalar impulse
Vr - The relative closing velocity of the two bodies
e - The coefficient of restitution
m1 - The mass of rigid body 1
m2 - The mass of rigid body 2
n - The normal of the collision
r1 - The vector from the center of mass of rigid body 1 to the point of collision
r2 - The vector from the center of mass of rigid body 2 to the point of collision
I1 - Inertia tensor for rigid body 1
I2 - Inertia tensor for rigid body 2

I apply the impulse to the objects as such.

v1 = v1 + (Jn)/m1
v2 = v2 + (-Jn)/m2
w1 = w1 + (r1 x Jn)/I1
w2 = w2 + (r2 x -Jn)/I2

v1 - Velocity of rigid body 1
v2 - Velocity of rigid body 2
w1 - Angular velocity of rigid body 1
w2 - Angular velocity of rigid body 2

I have tested this and everything works correctly. The thing that I'm clueless about is when I then decide to throw dynamic/kinetic friction into the equation. So far from the books I'm reading this is what I have been able to figure out.

This is the updated equations to update the rigid bodies.

v1 = v1 + [Jn + (uJ)t]/m1
v2 = v2 + [-Jn + (uJ)t]/m2
w1 = w1 + {r1 x [Jn + (uJ)t]}/I1
w2 = w2 + {r2 x [-Jn + (uJ)t]}/I2

u - The coefficient of friction
t - The tangent normal, this is a unit vector that is perpendicular to the contact normal.

The tangent normal is calculated with this equation.

t = [(n x Vr) x n]
t = t/|t|




I started out using the very same equations for friction, my source was the book "Physics for Game Developers - David M. Bourg". But I realized there is a problem with those equations; I was not able to simulate a tire(squashed sphere) rolling over a plane. Its easy to see why the equations would fail in this case:

J = -Vr.n(1+e) / {1/m1 + 1/m2 + n . [(r1 x n)/I1] x r1 + n . [(r2 x n)/I2] x r2}

It follows from the equation above that J is directly proportional to the normal velocity component. So if Vr is purely tangential to n - as in the case I described above, then J reduces to zero. And so does the change in linear and angular velocities.
Obviously this is not correct - if you gently place a rotating tire on the floor, you expect its linear velocity (forward) to increase and its angular velocity to decrease.

I've modified the impulse equation to handle this situation. Its best to treat the impulse as a vector rather than a scalar.

J = -Vr.n(1+e)n - Vr.t(u)t
ImpN = J; ImpN.Normalize();
J /= {1/m1 + 1/m2 + ImpN . [(r1 x ImpN)/I1] x r1 + ImpN . [(r2 x ImpN)/I2] x r2}

and ImpN = direction of impulse

Then,

dV = J/m
dL = rxJ

dL = change in angular momentum


These seem to work well enough..

#5 MrRowl   Crossbones+   -  Reputation: 1511

Like
0Likes
Like

Posted 22 September 2007 - 11:15 PM

"I was not able to simulate a tire(squashed sphere) rolling over a plane. Its easy to see why the equations would fail in this case..."

You can't expect equations derived for rigid bodies (i.e. a rigid sphere rolling on a rigid plane) to result in any rolling friction, since rolling friction comes from energy lost through compression/decompression of the objects.

The Guendelman paper describes a (presumably) more physically correct method for calculating rolling friction.



#6 chand81   Members   -  Reputation: 181

Like
0Likes
Like

Posted 23 September 2007 - 03:57 AM

Quote:
Original post by MrRowl
"I was not able to simulate a tire(squashed sphere) rolling over a plane. Its easy to see why the equations would fail in this case..."

You can't expect equations derived for rigid bodies (i.e. a rigid sphere rolling on a rigid plane) to result in any rolling friction, since rolling friction comes from energy lost through compression/decompression of the objects.

The Guendelman paper describes a (presumably) more physically correct method for calculating rolling friction.


Thanks MrRowl, can you please provide a link to necessary paper for simulating rolling friction?
I see your point, guess I need to spend more time reading papers than coding :)

KiroNeem, apologies if I put you off track a bit.

#7 BirdiePeep   Members   -  Reputation: 134

Like
0Likes
Like

Posted 23 September 2007 - 04:54 AM

xissburg - Thanks for the code, it has given me a little better of an idea of what I should be doing. Although seeing as angular motion is not taken into account It will still not completely work for what I'm trying to do.

DonDickieD - I started to print that PDF you linked to without even looking at how many pages it was. I came back to my computer with about two hundred pages on my desk and my printer beeping at me that it needed more paper. :p

This actually looks like a very good reference and I've begun to read over it, thank you very much. Once I find what I'm looking for I will post back my findings.

Chand81 - Actually, I might have some insight for you on what you are trying to do. From what I understand there are several variations of friction, the main ones in question are.

Static - This is a force that opposes motion of a body that is resting on another. This is related to rolling friction, but it is static friction that will cause a ball to roll across the ground when it has angular velocity.

Dynamic/Sliding/Kinetic - This is a force that will dampen motion of a body that is moving across another.

Rolling - This is a force that will slow a ball or tire from rolling due to deformation of the objects. This does not actually cause the linear motion as static friction does.

Here are some resources I have been reading over
http://en.wikipedia.org/wiki/Friction
http://en.wikipedia.org/wiki/Rolling_resistance
http://webphysics.davidson.edu/faculty/dmb/PY430/Friction/rolling.html

#8 Dirk Gregorius   Members   -  Reputation: 730

Like
0Likes
Like

Posted 23 September 2007 - 05:08 AM

Basically the friction impulse works like the collision impulse.

We seek a momentum conserving impulses...

v1' = v1 + M1^-1 * P1
w1' = w1 + I1^-1 * L1 = I1^-1 * ( r1 x P1 )

v2' = v2 + M2^-1 * P2
w2' = w2 + I2^-1 * L2 = I1^-1 * ( r2 x P2 )

...such that relative tangential velocity vanishes...

( v2 + w2 x r2 - (v1 + w1 x r1) ) * t = 0

...also knowing the direction of the impulses

P1 = -lambda * t
P2 = lambda * t = -P1

v := linear velocity
w := angualr velocity (read omega)
v' := linear velocity after impulse (read post-velocity)
w' := angular velocity after impulse
P := linear impulse applied at center of mass
L := angular impulse

r := vector from center of mass to contact point

w x r is the cross product of w and r ( read cross(w, r) )



If you have computed lambda this way you have found an impulse that eliminates the entire relative tantential velocity. Unfortunately there is a coupling to the normal impulse that means that you maybe cannot activate the full impulse, so you need to check this and clamp the impulse if required:

lambda < mu * lambda_n


The procedure for finding the impulses is usually the same. You blindly write down the impulse equation. Then you write the kinematic constraint that the post velocities must satisfy. Finally you usually also know the direction of impulse.


HTH,
-Dirk

#9 Dirk Gregorius   Members   -  Reputation: 730

Like
0Likes
Like

Posted 23 September 2007 - 05:12 AM

Chand:
"Thanks MrRowl, can you please provide a link to necessary paper for simulating rolling friction?"

Danny is refering to he second papers in my first post...


#10 BirdiePeep   Members   -  Reputation: 134

Like
0Likes
Like

Posted 29 September 2007 - 07:22 AM

"DonDickieD" - Thanks for your reply about how you went about handling friction. Based off of your implementation I was able to get something that is a lot closer to what I'm trying to do. Though, shouldn't time be taken into account when calculating how much friction acts on the tangent impulse? Or is that taken care of inherently by the fact that you are working with velocity instead of forces?

[Edited by - KiroNeem on September 29, 2007 3:22:42 PM]

#11 Dirk Gregorius   Members   -  Reputation: 730

Like
0Likes
Like

Posted 30 September 2007 - 12:24 AM

Sorry, I don't understand the question. Let me describe how you use my above formulas to deal with multiple simultaneous contacts. Maybe this helps:


for ( all contacts c )
Vec3 t1, t2;
PlaneSpace( c.Normal, t1, t2 ); // Find two mutally ortohgonal vectors

ApplyNormalImpulse( c );
ApplyFrictionImpulse( c, t1 );
ApplyFrictionImpulse( c, t2 );




One gotcha is the clamping of the normal and friction impulse. The straight forward way is to only to apply the normal impulse when then you have an approaching relative velocity at the contact point, but this naiive approach will result in jitter (see e.g. Erin Catto's Box2D slides for a more sorrow explanation). What actually is happening is that you are refining you global solution (the true contact impulse/force) by applying sequential delta impulses. So you clamp against the accumulated impulse and not the delta impulse. The same holds for the friction impulse, which is dependent on the normal impulse. So here again you accumulate the impulse and clamp this against the accumulated normal impulse.


Maybe you look at the mentioned papers by Erin Catto and let me know if you have further questions...


HTH,
-Dirk

#12 BirdiePeep   Members   -  Reputation: 134

Like
0Likes
Like

Posted 13 October 2007 - 02:06 PM

DonDickieD - I think there may have been a little misunderstanding about the last question. No problem though, after doing a bit more reading and testing I think the answer is what I had assumed earlier. However, I do again have to say thank you for all the help, honestly any piece of information can be helpful.

As far as the friction equation goes I believe I have finally figured it out, below is the algorithm. The main different from what I first posted is that the friction impulse is handled separately then the contact impulse. If anyone can verify what I'm doing is correct I would appreciate it.

...First we calculate the frictionless impulse

if(|Vrt| > u*J) //Dynamic friction
Jt = -u*J;
else //Static friction
Jt = -|Vrt|;

Jt = Jt/{1/m1 + 1/m2 + n . [(r1 x t)/I1] x r1 + n . [(r2 x t)/I2] x r2}

Jt - The scalar friction impulse
Vrt - The tangent component of the relative velocity
m1 - The mass of rigid body 1
m2 - The mass of rigid body 2
t - The tangent normal of the collision
r1 - The vector from the center of mass of rigid body 1 to the point of collision
r2 - The vector from the center of mass of rigid body 2 to the point of collision
I1 - Inertia tensor for rigid body 1
I2 - Inertia tensor for rigid body 2

Then you apply the friction impulse similar to the frictionless impulse.

v1 = v1 + (Jtt)/m1
v2 = v2 + (-Jtn)/m2
w1 = w1 + (r1 x Jtt)/I1
w2 = w2 + (r2 x -Jtt)/I2

v1 - Velocity of rigid body 1
v2 - Velocity of rigid body 2
w1 - Angular velocity of rigid body 1
w2 - Angular velocity of rigid body 2

This gives results that are really close to what I want, although there is one outstanding issue. The friction impulse will not apply to rotation around the relative contact point vector. This happens because of how the relative velocity is calculated. For example, if I have sphere spinning along the Y Axis on the ground then it will continue to spin forever.

Either this means that the above equation is wrong, or I just need to handle the friction for the rotation around the relative contact point separately (which is a simple task). As you can see below the cross product will basically cancel out any rotation along the axis created by the relative contact point.

Vr = (v1 + r1 x w1) - (v2 + r2 x w2)

Vr - The relative contact velocity

- Chase

#13 MrRowl   Crossbones+   -  Reputation: 1511

Like
0Likes
Like

Posted 14 October 2007 - 10:49 PM

Quote:
Original post by DonDickieD
ApplyNormalImpulse( c );
ApplyFrictionImpulse( c, t1 );
ApplyFrictionImpulse( c, t2 );


I'm pretty sure I managed to implement accumulated/clamped friction impulses with just one friction calculation... though I'm not exactly sure now! If so, it's in the jiglib source.

#14 chand81   Members   -  Reputation: 181

Like
0Likes
Like

Posted 15 October 2007 - 01:52 AM

Quote:
Original post by KiroNeem
This gives results that are really close to what I want, although there is one outstanding issue. The friction impulse will not apply to rotation around the relative contact point vector. This happens because of how the relative velocity is calculated. For example, if I have sphere spinning along the Y Axis on the ground then it will continue to spin forever.

The Guendelman paper describes a technique to handle Spinning & Rolling friction. The idea is very similar to the static & kinetic friction impulse, but this time a torque impulse is used - will only affect angular velocity and not linear velocity.


Quote:
Original post by MrRowl
Quote:
Original post by DonDickieD
ApplyNormalImpulse( c );
ApplyFrictionImpulse( c, t1 );
ApplyFrictionImpulse( c, t2 );


I'm pretty sure I managed to implement accumulated/clamped friction impulses with just one friction calculation... though I'm not exactly sure now! If so, it's in the jiglib source.


Thats true, I've seen the jiglib source as well. The idea is pretty simple, just project the velocity vector onto the contact surface to get the tangent component. The friction impulse is applied along the reverse tangent velocity direction.
I'm myself little confused why Dirk's implementation splits the friction impulse into 2 orthogonal components!!?


Here is my implementation, actually 2 implementations one based on the Guendelman paper and the other similar to JigLib. (both lot easier to read as currently only single dynamic object is considered)
Both approaches seem to give very similar results, I wonder if I have implemented the Guendelman approach incorrectly ??



//////////////////////////////////////////////////////////////////////////
//MrRowls Approach
//////////////////////////////////////////////////////////////////////////

gePt3Df Jn; float jn; //Normal Impulse Vector & Magnitude
float vn = Rvel.Dot(N); //Normal Speed
gePt3Df Vn = N * vn; //Normal Velocity

if(vn<0) //Approaching velocity
{
jn = vn*(-(1+elasticity));
jn /= massInv + ((cfg.iITwInv*(R.Cross(N))).Cross®).Dot(N);
Jn = N*jn;

dV += (Jn*massInv);
dL += R.Cross(Jn);
dw = cfg.iITwInv*dL;


//Rvel needs to be recomputed with the normal impulse having taken effect.
Rvel = (cfg.iAngVel+dw).Cross®;
Rvel += cfg.iLinVel + dV;
vn = Rvel.Dot(N);
Vn = N * vn; //Normal Velocity

gePt3Df Vt = Rvel - Vn; //Tangent Velocity
float vt = Vt.Length(); //Tangent speed
if( vt > 1e-7f)
{
gePt3Df Jt; float jt; //Tangent Impulse Vector & Magnitude
gePt3Df T = Vt/vt; //Tangent
jt = vt;
jt /= massInv + ((cfg.iITwInv*(R.Cross(T))).Cross®).Dot(T);

if (jt > jn*friction)
//Normal impulse is not pressing down hard enough.
//So the tangent impulse should allow slide.
jt = jn*friction;

Jt = T * (-jt);

dV += Jt*massInv;
dL += R.Cross(Jt);
dw = cfg.iITwInv * dL;
}

}






//////////////////////////////////////////////////////////////////////////
//Guendelman Approach
//////////////////////////////////////////////////////////////////////////
float vn = Rvel.Dot(N);

if(vn<0)
{
gePt3Df Vn = N * vn; //Normal Velocity

//Compute J such that Rvel's tangent component vanishes
gePt3Df J;
TMatrixf K(massInv,0,0,0, 0,massInv,0,0, 0,0,massInv,0, 0,0,0,1);
TMatrixf RSkew(0); RSkew.CreateSkew®;
K += RSkew.Transpose() * cfg.iITwInv * RSkew;
J = -Vn*elasticity - Rvel;
J = K.Inverse() * J;

float jn = J.Dot(N); //Impulse normal component
float jt = (J - N*jn).Length(); //Impulse tangent component

if (jt > friction*jn)
{
gePt3Df Vt = Rvel - Vn; //Tangent Velocity
float vt = Vt.Length(); //Tangent Velocity magnitude
gePt3Df T = Vt/vt; //Tangent

gePt3Df NUT = N - T*friction;
float j = -vn*(elasticity + 1);
j /= (K*NUT).Dot(N);
J = NUT*j;
}


dV += J*massInv;
dL += R.Cross(J);
dw = cfg.iITwInv*dL;

}




#15 Dirk Gregorius   Members   -  Reputation: 730

Like
0Likes
Like

Posted 15 October 2007 - 02:58 AM

Quote:
I'm pretty sure I managed to implement accumulated/clamped friction impulses with just one friction calculation... though I'm not exactly sure now! If so, it's in the jiglib source.


You did :-) You basically split the relative linear velocity at the contact into a normal and tangential component as in the Guendelman paper. I follow more the PGS approach implemented in terms of impulses and had some concerns with the warmstarting of the solver when not projecting onto a basis in the tangent plane. Actually I use central friction in practise so this is not a real issue anyway.

For those who are interested:
For central friction you only apply normal impulses at the contact points and the friction impulse is applied at the center of the contact patch/geometry. Since this doesn't account for twist (relative rotation around the contact normal) between the faces in contact you need to apply an additional (angular) twist impulse. For a contact patch of 4 points this gets you from 12 constraints (or impulses) down to seven.

#16 Dirk Gregorius   Members   -  Reputation: 730

Like
0Likes
Like

Posted 15 October 2007 - 03:08 AM

I remember now why I use two friction impulses. Basically I can't see how you can accumulate the friction impulse when the basis is changing during the iteration?
You see what I mean? In my opinion you need to project onto a staionary basis, while otherwise you are accumulating vectors (length) with changing direction which is not correct in my opinion.

Anyway it is good for the quality of the friction to align the first tangential axes with the relative velocity in contact plane (if it exists).

#17 erwincoumans   Members   -  Reputation: 272

Like
0Likes
Like

Posted 15 October 2007 - 07:25 AM

Quote:
Original post by DonDickieD
I remember now why I use two friction impulses. Basically I can't see how you can accumulate the friction impulse when the basis is changing during the iteration?
You see what I mean? In my opinion you need to project onto a staionary basis, while otherwise you are accumulating vectors (length) with changing direction which is not correct in my opinion.

Indeed. During the iterations, the velocity changes, and if you just one one friction direction, there will be no friction correction orthogonal to the normal and friction direction. This will result in very minor drift. This is noticeable in the Bullet CcdPhysicsDemo with a stack of cylinders (you can disable/comment out the second friction direction to see this).
Quote:

Anyway it is good for the quality of the friction to align the first tangential axes with the relative velocity in contact plane (if it exists).

Indeed, so in Bullet we use the projected normal onto the separating plane as first friction direction, and the orthogonal direction (normal.cross(firstFriction)) as second, unless the projected normal length is zero. In that case you still better use 2 orthogonal friction axis (during the iterations, the velocity changes and makes the projected normal non-zero).

Note that the straighforward accumulated impulse clipping results in 'box' shaped friction model, which in certain cases can cause objects to slide in unnatural directions: this happens when the friction gets clipped onto one of the box axis but not the other. This effect is reduced when using the projected normal as first friction direction.

Also note that we are making a lot of approximations. We can use the current normal impulse for the magnitude for the friction impulse. However, we have the 'chicken-and-egg' problem here, because friction and normal impulse have a two-way effect, ideally we solve a bigger system. However, the relaxation/iterative process usually converges fine and hides all those approximations.

To make a long story shorter: indeed, impulse based friction correction the same as normal collision correction, just using a different axis, and using an impulse magnitude dependent on its related normal collision impulse. That is why the friction has to happen at the end of all iterations. An alternative is to use the approximation of previous frames' normal impulse magnitude. Then the challenge becomes how to come up with an approximation the very first frame (otherwise a single-frame collision would not result in any friction).

To finish with the snippet the Bullet currently uses:

btVector3 frictionDir1 = vel - cp.m_normalWorldOnB * rel_vel;
btScalar lat_rel_vel = frictionDir1.length2();
if (lat_rel_vel > SIMD_EPSILON)//0.0f)
{
frictionDir1 /= btSqrt(lat_rel_vel);
addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);
btVector3 frictionDir2 = frictionDir1.cross(cp.m_normalWorldOnB);
frictionDir2.normalize();
addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);
} else
{
//re-calculate friction direction every frame, todo: check if this is really needed, because we use persistent contacts that keep the same normal
btVector3 frictionDir1,frictionDir2;
btPlaneSpace1(cp.m_normalWorldOnB,frictionDir1,frictionDir2);
addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);
addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);
}

Hope this helps,
Erwin
http://bulletphysics.com




Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.



PARTNERS