Is my understanding of LCP solvers correct for rigid body physics

Started by
4 comments, last by Dirk Gregorius 16 years, 5 months ago
A little while back I was working through implementing a impulse based physics simulation. For the most part it worked, although there were a few inherent issues with the method that made me want to explore the realm of simultaneous contact resolution. Since then I have been researching this subject to learn how it works and how to implement it. So far I have learned that to solve simultaneous collisions that I will need to implement a LCP solver. From my understanding a LCP solver is a method for solving a matrix of linear equations such that they all satisfy similar constraints. These constants are listed below. Af + b >= 0 f >= 0 f(Af + b) = 0 Where for a rigid body simulation... A - A matrix of scalar values with the number of columns and rows equal to the number of bodies involved. The equation for calculating the values for this is long, so I will leave it out. From the looks of it the value represents the the energy distribution between the two objects along the contact normal. f - The scalar magnitude of the resulting impulse due to collision pair. This is what the LCP will be solving for. b - A vector that represents the desired delta velocity of a collision pair. It can be found with the formula below. b = Vrel(1 + e) Where Vrel - The relative velocity between a collision pair. e - The coefficient of restitution for the collision pair. Once A and b are found you can then implement the LCP solver to find f. From the looks of it there are several different kinds of LCP solvers, each which work slightly differently. Unfortunately, I've only been able to find one good explanation (called the Dantzig algorithm) in the "Fast Contact Force Computation for Nonpenetration Rigid Bodies" by David Baraff 94. All I'm really wanting to know at the moment is if what I understand so far is true. Additionally, if anyone has any advice I'd also love to hear. Thanks, - Kiro wingsout.com
Advertisement
Usually you use an interative method called Projected Gauss-Seidel which is similar to the Gauss-Seidel for systems of linear equations, but takes the constraints into account. An efficient implemenation of this solver can be found e.g. in the ODE or Open Tissue.

The ODE also has an implementation of the Dantzig solver mentioned by you. Another direct method for LCPs is the Lemke method. I think there is a GDC presentation of Chris Hecker somewhere about this. Note that you need to change the solver such that it can solve box constraints (e.g. for friction or limits and motors). In your notation this would be lo <= f <= hi.

Note that the method you described is patented by Ageia and I think the box Dantzing solver from the ODE is is patented somewhere as well.

Another method would be to use pairwise relaxation using impulses like e.g. found in Bullet and nicely presented by Erin Catto (www.gphyiscs.com)


Please also note that Baraff solves on the acceleration level and not on the velocity level as you mention in your post...
Hey DonDickieD, I actually want to say thanks. You were extremely helpful in my last thread about calculating impulses, I appreciate everything.

"Usually you use an interative method called Projected Gauss-Seidel which is similar to the Gauss-Seidel for systems of linear equations, but takes the constraints into account. An efficient implemenation of this solver can be found e.g. in the ODE or Open Tissue."
I actually have a paper about how to do this method which seems mostly understandable, "Iterative Dynamics with Temporal Coherence" Erin Catto 95. At first I skipped the paper since other methods appeared to have better results, although I may take a closer look at it.

"Note that the method you described is patented by Ageia and I think the box Dantzing solver from the ODE is is patented somewhere as well."
Wha?! Thats a bummer. I did a little research and I believe this is the patent you are speaking about. Although now I'm a little confused, looking through it they appear to be using a modified Gauss-Seidel method.
http://www.freepatentsonline.com/7079145.html

"Another method would be to use pairwise relaxation using impulses like e.g. found in Bullet and nicely presented by Erin Catto (www.gphyiscs.com)"
Hmm, I thought that I was currently implementing something very similar to this although sure now. I'll need to do a bit more research.

"Please also note that Baraff solves on the acceleration level and not on the velocity level as you mention in your post..."
I was piecing my information together from a few different sources, so I must have missed that. ^^

- Kiro
wingsout.com
If you want to look into PGS you might also look at the PhD of Kenny Erleben. This with the paper you mention of Erin is propably the best introduction into this topic you might get. Basically both describe the internal workings of the ODE.

ftp://ftp.diku.dk/diku/image/publications/erleben.050401.pdf


With respect to constraints make sure you understand that a constraint on the position usually imposes also constraints onto the velocity and the acceleration. For example for the ball socket joint you want a point on body A be coincident with a point on body B. In order that the points don't drift apart their relative velocity at this point must vanish. This is the called the velocity constraint. Mathematically the velocity constraint is the time derivative of the position constraint. So this is how they are related. So when you fix both positions and velocities errors (e.g. Baumgarte stabilization) make sure that you fix the position error caused by an error in the velocity constraint (e.g. the penetration because of an approaching relative velocity at a contact point). This is very often done wrong for angular constraints where people intuively choose some angle to correct the error without defining a proper position constraint.

Basically when you have a constraint dynamic systems you have your equations of motion (Newton-Euler) which are differential equations and algebraic equations which define the constraints. These together build a system of differential and algebraic equations (DAE). I think for constraint dynamics they build what is called an index 2 problem (might be wrong here since I do this from the back of my head now). So in order to solve such a system you basically derive the position constraint twice to find the acceleration constraint which you can then plug into the equations of motion and compute the Lagrange mutlipliers. This is what Baraff does. If you are interested in games you can forget about this approach, though I recommend it to get an overview over the system you are solving from a learning perspective. Dave Eberly has also a nice explanation about this in his "Game Physics" book.

An example of solving on the position level you can look at the papers of Jacobsen and Muller:


http://www.teknikus.dk/tj/gdc2001.htm
http://www.matthiasmueller.info/publications/posBasedDyn.pdf


Note that the Muller approach is also patented so it just describes the Jakobsen method which is actually the Rattle algorithm from molecular dynamics.


Finally a good read which introduces somewhat the concept of pairwise relaxation is the PhD of Barenbrug. It predates actually all patents so I guess it is prior art.

http://www.win.tue.nl/dynamo/publications/bartbthesis.pdf


HTH,
-Dirk

Quote:Original post by DonDickieD
... I think there is a GDC presentation of Chris Hecker somewhere about this. ...


http://chrishecker.com/How_to_Simulate_a_Ponytail

Specifically, http://chrishecker.com/images/e/e8/Gdmag200003-ponytail-1.pdf
and http://chrishecker.com/images/e/e8/Gdmag200003-ponytail-2.pdf
No, I actually meant this:

http://chrishecker.com/The_Mixed_Linear_Complementarity_Problem

This topic is closed to new replies.

Advertisement