Baraff's Direct Solver

[avatar7]

Hi, all,

This is my first post here, I found this forum yesterday while googling, and was amazed to see so many famous developers and authors around.

For the last year or so I've been interested in constrained multi-body simulation, with accent over bi-lateral constraints rather than impact/contact. I believe I've gathered most of the free papers that are available, but while most of them use a velocity-based constraint approach and iterative solvers, Baraff's "Linear-Time Dynamics using Lagrange Multipliers (1996)" is aimed at the acceleration-based direct linear method of solving them.
I'm neither a mathematician nor a physics engineer, so most of the concepts are a bit complex to me, but most of that Baraff paper is really easy to follow, so in the last months I've been trying to write a small physics engine based on his algorithms and at first look it seemed it was working fine. All the multiplications and summations were correct, with proper-sized matrices etc., I even did them on paper myself, so unless there's a mistake in the paper itself, everything should've been working fine when I actually feed it real body and constraint data. It wasn't to be so though. The two "factor" and "solve" algorithms aside, they seem to be working perfectly, so either I'm not building the constraint Jacobians properly or I'm not using correctly the produced multipliers for getting the constraint forces.
My questions:

1) Has anybody actually written a code based on the algorithms from Baraff's paper and is it possible to take a look at some parts of it ?

2) Since in all papers I found they use velocity-based constraints, there's no explanation or an example of what exactly is the "c" vector (notion from Baraff's paper) that appears after differentiating a velocity constraint w.r.t. time. For example, in a ball-joint, I first had the impression I could leave the "c" vector at zero, since the constraint is holonomic, but I'm not sure at all if that's true. Atm, I'm just using the Jacobians produced on the velocity level, and I thought they should work for the acceleration level too. But they aren't and I'm not sure if it's something to do with deriving the "c" vector or it could be something wrong with either the algorithms or even something else.

I could try and make a small movie of my two bodies connected by a ball joint test, as it'd be better to see what's going on. It seems to me the torques are correct, the bodies rotate in a proper way in comparison to each other as if they're connected with the joint, but the linear velocity doesn't seem to follow the constraint at all. I tried manually entering static "c" vector values and it seemed to have an effect, so I believe I have to derive that each frame, but I simply have no idea what it's supposed to be.

As I said, I'm not much profficient in calculus, so please excuse my rather basic explanation. Thanks in advance for all answers!

[Erin Catto]

The term c arises from the chain rule of differentiation. Suppose you have the scalar _position_ constraint:

f(x) = x^2 + 2

The time derivative gives the velocity level constraint:

df/dt = x * dx/dt

So here the Jacobian is 'x'.

Differentiate once more to get the accelaration constraint:

d^2 f/dt^2 = x * d^2 x/dt^2 + (dx/dt)^2

The acceleration Jacobian is the same ('x'). The second term is the 'c' term.

[tasora]

Hallo,
I read the interesting paper of Baraff on the linear-time direct
solver (pay attention, however, that the linear-time complexity
happens only for acyclic topology - no closed kinematical chains)
This motivated the developement of the direct solver of my
multibody software CHRONO, which exploits the sparsity of
the system. It is a bit different from Baraff'one (still in case of acyclic
topology has linear complexity).
I use special sparse-matrix methods and an optimized LDLt
sparse decomposition (not as fast as Baraff method, but more
general, and if topology is not acyclic it becomes more and more
similar to a classic full LDLt decomposition. In most cases, however,
a fast linear-time behavior is obtained).
You can find it in:
"An optimized lagrangian-multiplier approach for interactive multibody simulation in kinematic and dynamical digital prototyping",
TASORA A.
VIII ISCSB, International Symposium on Computer Simulation in Biomechanics, 4-6 July 2001, Politecnico di Milano, Italy.
http://www.deltaknowledge.com/tasora/pu ... angian.pdf.
Also, I also recently adopted similar 'sparse' methods to obtain a
fast simplex solver for the LCP problem, but in general I am more
and more convinced that large-scale complementarity problems
(ex: with more than 100 bodies) should be attacked via an iterative
solver. I am testing a projected-SOR method, which I recently
implemented, and it performs very fast. Of course the precision
of an iterative solver is lower than the precision of a direct
solver...
Best regards

[avatar7]

Thanks for the replies!

After double and tripple checking today, and reading some other papers, I'm absolutely sure the Jacobian I use for my ball-joint is just as it should be, and updating it each frame is done properly too. Also, the 'c' term, as mentioned in a paper about how to make joints in ODE, is zero all the time. If I understood correctly, it's only non-zero when the constraint is non-holonomic or when it's used for correcting drift. So it must be a problem in the algorithms somewhere in Baraff's paper. Actually, I already found a mistake in it, but it was a small one. But perhaps a negative sign or something somewhere else is still causing the problems. Since I don't know much about sparse matrix and solvers theory, I can just follow the paper's pseudo code and hope it's correct when I resemble it in actual c++ code. That's why I was asking if someone has actually written those algorithms already, so they could point out if there are any mistakes in the actual paper. Would be very helpful !
I know an iterative solver is much better for games and real-time simulations when there are alot of bodies involved and when accuracy isn't so important, but in what I'm doing there are just a few bodies, but accuracy is of major importance. That's why I believe a direct solver would do much better job. Thanks, Tasora, for the paper, I'll check it out.

[Erin Catto]

tasora,

I enjoyed reading your paper. The sparse LDL is interesting. I was surprised you were able to deal with redundancies so easily. Have you tried simulating a Bricard mechanism with this algorithm? It is notorously difficult to simulate due to the constraint redundancy.

http://www.tu-chemnitz.de/ifm/english/e ... ricard.htm

I am curious about your sparse LCP solver. Do you plan to share your algorithm?

[kenny]

1) Has anybody actually written a code based on the algorithms from Baraff's paper and is it possible to take a look at some parts of it ?

I think Baraff and Witkin used to have a company, where they among other things made the initial rigid body simulator used in Maya. I think it was called Coriolis... I know from Baraff that Maya is still based on his method. Anyways if you play around with Maya, you'll have to agree it works pretty good, except for two things it scales badly and it can in no-way handle initially penetrating objects.

A few years ago we did have an implementation of Baraffs method in OpenTissue, but I have no idea of what happen with it. I know of many people who have tried to implement Baraff's 94 paper, but today they have all shifted to velocity-based formulations.

The major difference between Baraff 94 method and the typical state-of the art velocity based solution used by many today is:

1) The direct method used for the LCP
2) The back-tracking method needed to find TOI
3) Switching to impulse-momentum law
4) The need for the time derivatives of the jacobians

1) leads to O(n^3) complexity at best and O(n^4) at worst, 2) basically means you have no control over how many steps is needed to complete a single frame computation. You can even get into an infinite loop (this is what happens with initial penetrating objects). A more timely solution may be to replace the bisection method with something like CCD. 4) Is in my opinion really tedious, if you actually derive the equations they get widely, so one often end up with using approximations.