One thing he mentions on his site which I found to work quite well is to form the MLCP of the full contact manifold per collision pair and solve it explicitly.

Basically, suppose you are in 2D and you have two boxes stacking. You form the full Jacobian of this system as:

[ -n1, -n1 x rAP1, n1, n1 x rBP1]
[ -n2, -n2 x rAP2, n2, n2 x rBP2]

(That is, it's a single large rectangular matrix).

Then you multiply out J * M-1 * JT to form a 2x2 matrix. Then you explicitly solve this 2x2 matrix to find your two contact impulses.

If they're both > 0, you're done. If one is < 0, clamp it to 0 and remove it from the system and solve the resulting 1x1 matrix for the other contact in isolation. If that one is less than 0, then clamp that one as well and you know that there are no impulses involved.

You can likewise expand this to include friction, but it's more complicated (still working through how it would work myself). Either way, it gives good stacking behavior and far less jitter.

Ignore for the moment that it's an LCP. You have a 2x2 system of equations. You can solve it using cramer's rule. You'll get two solution impulses: one for each contact point.

Now remember that it's an LCP, so the answers must all be >0. Check your two impulses. If they're both > 0, you're done. If one of them is <0, you need to feed your problem into an LCP solver. Since the system is so small, you can basically brute force it (direct enumeration method).

So you take the impulse that was <0, and you set it to 0 (using the complementarity condition, you're basically setting a slack variable to some non zero, positive value instead). Then you solve for the remaining impulse term as if it were the only contact. If it's <0, you clamp it to 0 as well.

...
Here's a web book on LCP. If you understand section 1.1 you'll be in good shape. It's in postscript format, so you'll need to download and install ghostscript and ghostview. If the download links are broken (they are sometimes), look through the download mirrors they have listed.

But here's the crib notes for it: In LCP format, basically you have Iw - Mz = q. M is your 2x2 matrix that you found from the Jacobian. q is the vector of initial conditions (ie: each element is (1+e)*(delta V * n)). I is the identity matrix, and w is a "slack" variable. If you remove Iw from the equation, you basically have your system of equations Mz = -q. Which you can solve normally using Cramer's rule. (The M matrix here is your J * M-1 * JT matrix).

The LCP has a non negativity condition (which is what you want in this case for collision impulses). So elements in z must be >0. If any are <0, they get clamped to 0.

The LCP also has what's called a complementarity condition. Which means that if an element in z is 0, the corresponding element in w is >= 0. So when you have to clamp an element in z to 0, you increase the corresponding element in w so the system will still solve correctly.

Solving the system is basically a combinatorial problem. You want to find what elements in z can be >0, and which need to be clamped to 0 and have the slack vector w handle it.

At the end of the process, you just discard w and the z vector are your impulses. For larger problems, you have to use more sophisticated solvers than direct enumeration, since it has an exponential order of complexity. But since this is a simple 2x2 system, you can do the algorithm I outlined at the beginning of the post: solve the system normally, and if any of the impulses are <0, clamp it to 0 and solve the remaining contact point in isolation.

Using this method in my little homebrew engine I was able to stack over 200 boxes in real time (~55-60 FPS) and have them be stable. You have to add friction to keep them from slipping out laterally, though (imagine stacking wet soap). And you need Baumgarte stabilization, too, or the stack will start to sag (imagine stacking blocks of depleted uranium on top of blocks of styrofoam). But that was it in my engine. No sequential impulses, etc.

I'm not sure if you need to add friction to the LCP, though, or if you can just sort of fake it. I tried adding it to my LCP, but I had to hack it to get it to work right, and even then it only works when the two bodies have a two point contact manifold. If it's vertex vs. edge it won't apply any friction for some reason.

When I first thought about friction, I was going to form a friction contact at each real contact (basically form the perpendicular to each contact normal for each contact point, and add that constraint). So for a two point contact manifold that would give 4 constraints. I would then use my matrix library based on SIMD and invert the matrix to solve the system.

However, if you try doing that you quickly run into a problem I didn't understand at the time: the matrix will be singular. The reason is actually sort of interesting.

So imagine you're in 2D. There are three degrees of freedom per body. So imagine now you have two bodies. There are six degrees of freedom for the system. Now add a constraint between the bodies. The constraint, no matter what it is, removes a single degree of freedom from both bodies. So each (binary) constraint removes two degrees of freedom from the system.

So in the 2D friction case I was working on, I was trying to add 4 constraints, which would remove a total of 8 degrees of freedom. But there's only 6 degrees of freedom total, so the system ends up being over constrained and no solution exists.

In my engine I just dropped one of the friction constraints arbitrarily. And if any of the contact impulses were <0, I clamped all of them to 0. A really terrible hack, but it worked well enough.

If you're building on my work, there are a couple of ways I can think to approach the problem:

1. First, solve the 2x2 LCP without friction, just from the contact impulses. Then solve the 2x2 friction linear system (without any constraints. Pretend you have infinite friction coefficients) for your friction impulses. Then clamp the friction impulses based on the values you got from the contact impulses. So if(contact impulse * friction coefficient < abs(friction impulse)) -> clamp friction impulse so its absolute value is equal to contact impulse * friction coefficient. It's not mathematically correct to handle it this way, but it might work just fine in practice.

2. Take the two constraints from the two point contact manifold, and determine a new constraint which removes the last two degrees of freedom. This is your implicit friction term. You can now solve the critically constrained (by which I mean you've removed exactly all degrees of freedom from the system) system. You then enforce the constraints that the two contact impulses must be >0, and the friction impulse must be bounded by +/- (coefficient friction * sum of contact impulses).

Unfortunately I have no idea how to solve this problem. It's not a LCP. I don't even think it's a MLCP (mixed linear complimentarity problem). It might be something individual solvers can handle (Lemke, Simplex, etc.), though.

If you graph the problem, the "feasible region" (the area in which a solution might exist given the constraints) is the intersection of cones and half planes, so it's very well defined. I'm thinking you might be able to do something like "find the closest point in the feasible region to what the solution would be if it were unconstrained", which is a least squares quadratic programming problem that's fairly well defined. You might even be able to use a modified GJK algorithm for it. I'm not really sure to be honest.

3. The ultimate in faking it: just take your sequential impulse solver engine exactly as is, and just replace solving the contact impulses in a manifold sequentially with the 2x2 LCP solver. Your friction terms are found sequentially just like normal. You'd probably be better off solving the friction term as if there's only one per contact manifold though (centroid (center) of the contact manifold).

...

Beyond that I have no idea.

[Edited by - Numsgil on July 27, 2009 3:15:46 PM]

Quote:Original post by AHashem
I'd also like to know your opinion in accumlating the impulses and using the features of collision to store collision points among frames...do you do this in your engine??

No, but my engine was pretty ghetto :) Basically I was applying impulses as soon as I found them instead of storing them across iterations. I wasn't doing any warm starting, either. The only caching I was doing was only performing collision detection once per frame.

You'll get better results (probably) with both methods.