The merits of projecting twice

I started a couple of years ago writing a post called "The trouble with PBD contacts", but never finished. This was based on my practical experience of modeling contact. I hope I will finish that post one day, once I learn all the factors causing all that trouble. But long story short, PBD contacts are pretty crap and usually quite violent. As Stewart puts it: "the resulting discretization behaves as if it had a "random" coefficient of restitution when impacts occur". And the same applies to stabilized velocity contact constraints. In practice, this manifests as objects popping out sometimes at high speed. I call them "violent contacts". Bottom line: the velocity at impact should always be zero and this should be our constraint (unless we want to model non-zero restitution, but we can add that afterwards). If we do that, then all our contacts will be just fine and no explosions will appear. But what about the penetration error? I'm sure you will ask that as nobody likes objects on screen clipping each other. Stay tuned.

I will focus here mostly on what I call position-based methods, which is different from position-based dynamics (PBD), but it includes it. There are generally two types of such methods: explicit and implicit, pertaining of course to the type of integrator used. But, more importantly, is that it also makes the distinction between explicit and implicit constraint directions. Meaning, when is the direction of the constraint forces computed: at the beginning or the end of the time-step? In the first category I put SHAKE and RATTLE which are the darlings of molecular dynamics. They are explicit symplectic integrators of constrained systems. Hence, they can easily blow up or end up without a solution for large time steps, but people live with it as long as energy is conserved. On the other hand, implicit methods tend to be much more stable, as by construction they dissipate energy (and people live with it). This is why they are used so much in graphics and gaming. In this category I put the so-called "projection" methods, where you integrate up to a candidate position (without constraints) and then project the positions back onto the constraint manifold. They are also known as step-and-project (SAP). Here's where I put PBD and "fast projection" too, as they are approximate solvers for SAP. In fact, I will often say PBD instead of SAP, unless I make it clear that I'm referring to the nonlinear Gauss-Seidel projection popularized by Nvidia.

I will not go into the nitty gritty details of equations, but there is one important point I want to make. Neither SHAKE nor PBD preserve the hidden constraint that velocities along constraint directions should be zero (the normal component to the constraint manifold or the first derivative of the constraint equation). And I will prove it with a plot of the velocity constraint drift for the case of a simple pendulum.

RATTLE is the only scheme that calculates the correct velocity. And it does it at the price of an extra projection step at the end for the velocity. Given that the position constraint (explicit) projection at the beginning of RATTLE is nonlinear, and the velocity projection is always linear, it seems like a small price to pay. It's only one linear system to solve, compared to several if using Newton for the positions. So, here's what I did, I took RATTLE and I replaced the first part with PBD. And the results were great in terms of velocity. I called it PROJ2 for lack of a better name and because I project twice: first the positions and then the velocities.

But why should we bother? SHAKE and PBD work perfectly well the way they are. Indeed, if you are dealing with equality constraints. But if you are dealing with contacts (and friction) velocity starts to matter more. That high velocity error in the plot will become the impact reflected velocity. This is not an issue with equality constraints, as no matter what velocity they have along the constraint direction, the bodies will stay pinned to the constrained manifold (up to a solver tolerance). So, this is why I modified RATTLE like I did. It's basically the modification you need to do to PBD in order to make it explosion free for contacts.

I have not not tried PROJ2 on practical contacts yet, but I did something similar in the past. I projected the velocities before running PBD and everything was great. If you think about it, it's the same thing, if you move the last step of RATTLE at the beginning of the next frame. So, I know it works! It made contacts much smoother. The only reluctance I had, like everybody else, was that it was more expensive. And now I know the answer: it doesn't matter! It's worth the price if your contacts are fully inelastic and not exploding. And you also get rid of that stabilization term... 

Think about it: first project velocities and then positions. Does it sound familiar? The latter step is sometimes called "post-projection": where you correct the positions after you solve the problem at velocity level in order to remove any penetration. This technique is usually attributed to a paper by Cline and Pai, but ultimately it's just projection. Usually the problem is linearized, but when dealing with contact, the normals are usually considered constant throughout the solve, so it is linear. Of course, it would be much better if you used PBD instead, but the border between post-projection and PBD is quite gray to be honest anyway.

One more thing, what's wrong with stabilization? Nothing, if you are writing a rigid body simulator where you don't care so much about small penetrations. And then you can use a small stabilization factor and clamp the velocity to avoid explosions. But if you care about having zero penetrations and correct normal velocities, you should use the post-projection step. In fact, this has been a debate for years in the community, and for that reason Box2D uses post-projection, while most other engines uses stabilization or a mixture of the two. Now, I tend to agree with Box2D. By the way, if you care about that little analogy between implicit springs and the stabilization factor and the "constraint force mixing" factor, I would say forget it - I think it's pretty much wrong anyway (or at least not rigorous). 

In the past, I feared that post-projection dissipates artificially too much energy. But this phase plot for the pendulum made me feel better.

In fact, it preserves energy much better than PBD. So, it has two merits: it reduces dissipation for PBD and solves penetration problem for velocity level methods. From now on I think this is what I'm going to use as a contact solver (and constraints in general), no matter how expensive it is.