Friction and LCP

[tkcastro]

Hello,
I'm developing a physics engine at my university (UFJF, Brazil) and I've followed Eberly's book, Game Physics.
Now, i can deal with multiple contacts and resting contacts with restitution using LCP (Lemke-Howson algorithm), but I have some questions and I'm asking for your help:
-The first LCP Solver receives the A matrix and a vector called "Vector B" and gives me the post impulse velocities and the impulses that i have to apply. (This Vector B has elements bi = 0 if (preImpulseVelocity >=0 ), else it is bi= 2 * preImpulseVelocity. I didn't understand something: Why i have to multiply the preImpulseVelocity by 2?)


-Now I need to add friction, but i don't know exactly what to do.. I've found a lot of papers, but i'm a little confused yet. For example, in http://www.dh.aist.go.jp/~kawachi/ca98/ca98kawachi.pdf he says:

• i have to decompose my velocity into normal and tangent variables,
  that's ok, all i have to do is to project my preImpulseVelocity into the normal and the tangential direction, right?
• well, i have doubled the size of my preImpulseVelocity vector, so i need to double the size of A matrix, but how? what's the new formula for it?

Thank you very much ,
TK

[jiangwei]

Hi
Do not use Lemke method to solving LCP,its slow and can not be control.
So PGS or SI (developed by Erin catto)is better choice to solving velocity-level constraints
But i think the position-level theory is the future

[Dirk Gregorius]

But i think the position-level theory is the future

Can you give arguments for this? How will you handle friction and motors (or non-holonomic constraints in general)?

[KenB]

Yes, besides the obvious modeling limitations, I am also curios what the principles would be that support such a conclusion. Position based methods simply neglect most of the right hand side of the first principle equations, and thus can be formally shown to be an approximate limit of the full dynamic theory. They do not contribute with any new fundamental insight or new principles of simulation - they just choose one of the simplest approximations - and they also give obvious artefacts.
However, they are simple to implement and grasp, and they do a reasonable job in certain a context, so I tend to like them for this (but whether I like them or not, they don't solve my problems...).

In addition (in a comment to a previous post in this thread), PGS (or SI if you wish) has been around for many years, and Catto's formulation is a slightly simplified version of Claude Lacoursiére's paper (it's using a simpler stepping equation) from 3-4 years ago (that doesn't claim to be first with PGS either). Catto had a reference to this paper on his web page for a while, but I wasn't able to find it now. Catto has done a great job making this available to the community!

[Erwin Coumans]

PGS (or SI if you wish) has been around for many years, and Catto's formulation is a slightly simplified version of Claude Lacoursiére's paper (it's using a simpler stepping equation) from 3-4 years ago (that doesn't claim to be first with PGS either).

Moreau is the first publication I'm aware of that mentioned use of Projected Gauss Seidel in the context of rigid body dynamics, or is there any earlier?

Several people realized that PGS and sequential impulse methods are essentially the same. Erin's contribution is mainly distributing existing knowledge to a wider audience through his excellent game developers conference slides/papers.

As mentioned in the recommended resources, the work from Kenny Erleben and the good introduction book/thesis by Helmut Garstenauer are highly recommended:

• Helmut Garstenauer, A Unified Framework for Rigid Body Dynamics
• Kenny Erleben, "Stable, Robust, and Versatile Multibody Dynamics Animation", ftp://ftp.diku.dk/diku/image/publicatio ... 050401.pdf

Given that you already have the contact normal, and applied normal impulse, this simplified friction model (in Bullet) works pretty well in practice:

Code: Select all

{
	btVector3 frictionDir1 = vel - cp.m_normalWorldOnB * rel_vel;
	btScalar lat_rel_vel = frictionDir1.length2();
	if (lat_rel_vel > SIMD_EPSILON)//0.0f)
	{
		//use projected normal direction as first friction direction and cross product of normal and frictionDir1 as second friction direction
		frictionDir1 /= btSqrt(lat_rel_vel);
		addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
		btVector3 frictionDir2 = frictionDir1.cross(cp.m_normalWorldOnB);
		frictionDir2.normalize();
		addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
	} else
	{
		//calculate two orthogonal vectors to normal direction
		//re-calculate friction direction every frame, todo: check if this is really needed
		btVector3	frictionDir1,frictionDir2;
		btPlaneSpace1(cp.m_normalWorldOnB,frictionDir1,frictionDir2);
		addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
		addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
	}

}

Thanks,
Erwin

[Erin Catto]

Correct, I did not invent PGS or SI. I just kind of merged the two and tried to present the hybrid algorithm in a digestible form along with some simple code to demonstrate the theory. Also, I did not come up with the idea of the equivalence of PGS and SI, I just made it work in some publicly available code and documented the steps to make it work. Indeed, much of my inspiration for this came from Claude's paper on splitting methods and from suggestions made by my friend Gary Snethen.

I know this has caused a lot of controversy about proper credit. But as with most research, improvements are incremental and this is where my work on Box2D+SI lies.

I guess if people want to refer to my algorithms, they should refer to Box2D rather than SI/PGS in general.

[Erin Catto]

Several people realized that PGS and sequential impulse methods are essentially the same. Erin's contribution is mainly distributing existing knowledge to a wider audience through his excellent game developers conference slides/papers.

Thanks Erwin. I agree with this mostly, except that you make it seem as if I simply copied someone else's algorithm and published it as my own. I would rather say that I took existing ideas from PGS and made them work in the context of an SI algorithm. A crucial ingredient to this is clamping of accumulated impulses, which I had not seen published anywhere. This is not simply a bug fix, but it is a fundamental ingredient for making SI convergent.

Further, many SI algorithms did not show how to handle generic joints and position error (i.e. Mirtich, Guendelman, etc). The original version of Box2D and my 2006 slides showed how to handle these cases. Again, these were not new ideas, but AFAIK at the time no one had published the details of how to make these work in an SI setting.

Much of this comes down to semantics. What is an algorithm? Well Bullet and Box2D both use a version of SI and both are somewhat different in the details. Casey Muratori presented a very nice form of GJK, but he did not invent GJK. In some ways the form of an algorithm is just as relevant as the core algorithm. For example, some forms of SI are not convergent.

So I don't take credit for PGS/SI at all, but I am responsible for the form of SI in Box2D and in my published works. Can you agree with this?

[Antonio Martini]

A crucial ingredient to this is clamping of accumulated impulses, which I had not seen published anywhere.

yes this is what i thought that it was new when i read your paper/presentation for the first time.

I dont remember the details of the Claude's paper, can anybody point out in which paper and where in the paper the method published by Erin was described?

what i dont understand is: if we have different SI methods let's say methods A and B and they are equivalent to PGS then it must mean that A = B. How can we have A != B but both equivalent to PGS? so (PGS = A) AND (PGS = B) AND A != B ?????

for example we are saying that clamping or not clamping they are both equivalent to PGS? so why do we clamp? clamping should have not effect given the previous argument. This definition of equivalence looks a bit too broad for my taste.

cheers,
Antonio

[Dirk Gregorius]

Erin's implementation of SI is equivalent to PGS. SI like described by Mirtich or Guendelman are not. Actually the differences only holds for inequality constraints when you need to clamp. All SI implementations are equivalent to PGS for equality constraints. Erin was indeed the first to publish a correct and formal description of it in his papers, but Bullet used SI for joints already at this point to my knowledge. Still there are differences between Bullet's implementation and Erin's approach (e.g. in the hinge constraints). I think R. Weinstein also published a paper on this as a follow up on the Guendelman paper, but this was after Erin if I remember correctly. We should also mention that the Doom3 SDK has an impulse solver, but I don't remember if it uses accumulated impulses, nor when it was published.

Anything I did forget?

[Antonio Martini]

Erin's implementation of SI is equivalent to PGS. SI like described by Mirtich or Guendelman are not. Actually the differences only holds for inequality constraints when you need to clamp. All SI implementations are equivalent to PGS for equality constraints. Erin was indeed the first to publish a correct and formal description of it in his papers, but Bullet used SI for joints already at this point to my knowledge. Still there are differences between Bullet's implementation and Erin's approach (e.g. in the hinge constraints). I think R. Weinstein also published a paper on this as a follow up on the Guendelman paper, but this was after Erin if I remember correctly. We should also mention that the Doom3 SDK has an impulse solver, but I don't remember if it uses accumulated impulses, nor when it was published.

Anything I did forget?

the definition of PGS includes inequality constraints, if we state that an algorithm is equivalent to PGS when only equality constraints are present, we are saying that it is equivalent to GS not PGS.
so you are saying that the Erin's version is equivalent to PGS while the other versions you mention are equivalent to GS for equality constraints + something else that is not PGS for inequality constraints. Right?

cheers,
Antonio

[Dirk Gregorius]

I think you can say it like this

[Erwin Coumans]

Main argument here is wether the PGS algorithm differs from the SI algorithm, or not. To me the algorithms are the same. If we are to talk about credit, as far as I can track down, PGS was first introduced within rigid body dynamics context by Moreau in 1999.

I agree with this mostly, except that you make it seem as if I simply copied someone else's algorithm and published it as my own. I would rather say that I took existing ideas from PGS and made them work in the context of an SI algorithm

Let's call this algorithm PGS. As far as I know, your publication described one route starting from relaxation scheme (sequentially applying impulses) that reaches PGS. You gave this sequential impulse route its own name, but it is still PGS in my opinion. I didn't intent to make it look like you published PGS as 'your own'. Your main contribution as far as I see it was showing the importance of clamping the accumulated impulse. This makes sequential impulse based methods, as implemented now in Bullet and Box2D, equivalent to PGS.

This is not simply a bug fix, but it is a fundamental ingredient for making SI convergent.

This is a fundamential ingredient for making PGS convergent. Some SI implementations lacked this accumulated clamping, which was the only missing step to reach PGS. In Bullet's PGS/SI implemenation I failed to recognized this clamping difference with quickstep, the ODE PGS implementation. With your help, this issue was fixed at the end of 2005.

Further, many SI algorithms did not show how to handle generic joints and position error (i.e. Mirtich, Guendelman, etc). The original version of Box2D and my 2006 slides showed how to handle these cases. Again, these were not new ideas, but AFAIK at the time no one had published the details of how to make these work in an SI setting.

Open Dynamics Engine quickstep, a PGS implementation had joints around 2004, I think. Bullet PGS/SI solver had generic joints in the PGS/SI context, a point to point constraint and hinge joint was implemented and published as open sourced in July 2005.

There are many other parts in the constraint setup and solving process that affect the quality more then just this accumulated impulse. Friction constraint formulation is one of them, and that was the original topic of this thread.

[Dirk Gregorius]

Did Moreau really introduce PGS? I read this paper a long time ago, but I think I remember that it deals more with pairwise relaxation similar to a GS style solver?

[Erwin Coumans]

Did Moreau really introduce PGS? I read this paper a long time ago, but I think I remember that it deals more with pairwise relaxation
similar to a GS style solver?

Well, Moreau discusses Gauss Seidel within the context of unilateral contacts.

9.2: Here is an iteration technique a la Gauss-Seidel which consists of treating a succession of single-contact problems.

Wouldn't contacts become a bit sticky without the projection part? I assumed the 'a la' bit includes projection. In fact, the contact law that Moreau mentions, see (40) in Numerical aspects of the sweeping process includes the complementarity condition. To quote his paper: if f(t,x)<0 this imposes no restriction on velocity U, and if f(t,q)>=0 this imposes U.n>=0. The first f(t,x)<0 should probably read f(t,q), with t=time and q=position, perhaps a typo?

Thanks,
Erwin

By the way: note that once we have bilateral and unilateral, we should talk about a mixed lcp (MLCP). Also, note that quickstep is actually implementing successive overrelaxation (SOR), very similar to PGS.

[tkcastro]

Thanks everyone for helping me

Since i'm a graduate student just in the beginning of this physics journey sketching a simple physics engine for study purposes i hope you understand my difficulties

Given that you already have the contact normal, and applied normal impulse, this simplified friction model (in Bullet) works pretty well in practice

I've read those papers and i already have the contact normal, and applied normal impulse. I have this LCP Solver and i just don't know how to calculate q, w, M and z. There are some papers explaining this, but i've tried and i didn't get some things, e.g in Erleben's work i need do compute D or E matrix, and they are diagonal?

Do not use Lemke method to solving LCP,its slow and can not be control.
So PGS or SI (developed by Erin catto)is better choice to solving velocity-level constraints

I was following Eberly's Game Physics until now.I just want a fast, stable and simple way to handle friction and all my constraints. Is the PGS method stable and fast, capable of handling all the constraints i have?


Thank you again (someday i hope will help other people with this too )