Archive for the 'General' Category

Block Solver

Box2D now implements a block solver for normal contact. If the contact manifold consists of two points, then the 2-by-2 MLCP is solved directly through enumeration. For small systems enumeration is easy to implement:

  1. Assume both constraints are active and solve the 2-by-2 equality problem. If both contact forces are positive then we are done.
  2. Assume contact 1 is active and contact 2 is inactive. Solve for contact force 1. If contact force 1 is positive and contact velocity 2 is non-negative then we are done.
  3. Same as step 2 except reversed.
  4. Set contact force 1 and 2 to zero.

Dirk Gregorius originally implemented this in Box2D_Lite in Jan. 07. I devised a scheme to make it work with accumlated impulses so that each block solve considers the total impulse applied over the step. Recently I ported this to Box2D 2.0.

The result is quite good and Box2D is now able to stack 24 boxes vertically without difficultly (1m cubes, 60Hz, 10 iterations, gravity = 10m/s^2). The behavior is similar to sphere stacking.

One additional contribution: you need to check the condition number of the block to avoid numerical problems. In the 2-by-2 case this is simple. If the condition number is too large then Box2D will ignore the second contact (because the constraints are somehow redundant).

GDC 2008

I once again participated in a physics tutorial at the GDC with the Essential Math crew.

You can download my slides here.

Just in case you missed it, Crayon Physics Deluxe won the IGF. Congratulations to Petri!

Position Correction

I tried the several algorithms for position correction of Box2D’s revolute joint.

I looked at these systems:

  • simple pendulum (1m diameter sphere on a massless 5m stick) with an initial angular velocity of 100 rad/s.
  • suspension bridge with 30 1m long planks.
  • multi-link chain with 30 1m long links.

Here are the algorithms:

Baumgarte - A fraction of the position error is added to the velocity error. There is no separate position solver.

Pseudo Velocities - After the velocity solver and position integration, the position error, Jacobian, and effective mass are recomputed. Then the velocity constraints are solved with pseudo velocities and a fraction of the position error is added to the pseudo velocity error. The pseudo velocities are initialized to zero and there is no warm-starting. After the position solver, the pseudo velocities are added to the positions.

This is also called the First Order World method or the Position LCP method. It is also very similar to the obsolete split impulses technique.

Modified Nonlinear Gauss-Seidel (NGS) - Like Pseudo Velocities except the position error is re-computed for each constraint and the positions are updated after the constraint is solved. The radius vectors (aka Jacobians) are recomputed too (otherwise the algorithm has horrible instability). The pseudo velocity states are not needed because they are effectively zero at the beginning of each iteration. Since we have the current position error, we allow the iterations to terminate early if the error becomes smaller than b2_linearSlop.

Full NGS or just NGS - Like Modified NGS except the effective mass is re-computed each time a constraint is solved.

I created four demos to test the different solvers. You can download them here.

Observations:

Baumgarte - this is the cheapest algorithm but it has some stability problems, especially with the bridge. The chain links separate easily close to the root and they jitter as they struggle to pull together. This is one of the most common methods in the field. The big drawback is that the position correction artificially affects the momentum, thus leading to instabilities and false bounce. I used a bias factor of 0.2. A larger bias factor makes the bridge less stable, a smaller factor makes joints and contacts more spongy.

Pseudo Velocities - this is more stable than the Baumgarte method. The bridge is stable. However, joints still separate with large angular velocities. Drag the simple pendulum in a circle quickly and the joint will separate. The chain separates easily and does not recover. I used a bias factor of 0.2. A larger value lead to the bridge collapsing when a heavy cube drops on it.

Modified NGS - this algorithm is better in some ways than Baumgarte and Pseudo Velocities, but in other ways it is worse. The bridge and chain are much more stable, but the simple pendulum goes unstable at high angular velocities. I used a bias factor 0f 1, so no tuning was necessary.

Full NGS - stable in all tests. The joints display good stiffness. The bridge still sags, but this is better than infinite forces. I used a bias factor of 1.

Recommendations:

Pseudo Velocities are not really worthwhile because the bridge and chain cannot recover from joint separation. In other cases the benefit over Baumgarte is small. Modified NGS is not a robust method for the revolute joint due to the violent instability seen in the simple pendulum. Perhaps it is viable with other constraint types, especially scalar constraints where the effective mass is a scalar.

This leaves Baumgarte and Full NGS. Baumgarte has small, but manageable instabilities and is very fast. I don’t think we can escape Baumgarte, especially in highly demanding cases where high constraint fidelity is not needed. Full NGS is robust and easy on the eyes. I recommend this as an option for higher fidelity simulation and certainly for suspension bridges and long chains.

Full NGS might be a good choice for ragdolls, especially motorized ragdolls where joint separation can be problematic. The number of NGS iterations can be reduced for better performance without harming robustness much.

NGS is tricky to implement for contacts because recomputing contact points can be expensive, however it is possible to estimate the new penetration depth based on the body movement and the contact normal. Another possibility would be to cache features and quickly redo the feature clipping to generate new points. Some points might disappear during iterations, so there are certainly some edge cases to consider.

Each joint in a can be handled differently in the position solver. So I recommend a system where the user can select the algorithm on a per joint basis. I would probably default to the slower Full NGS and let the user select the faster Baumgarte method in performance critical scenarios.

Buoyancy

I wrote a chapter in Game Programming Gems 6 about the buoyancy algorithm I used in Tomb Raider. The algorithm is suited for polyhedra on flat water surfaces. I posted the code on the downloads page.

Box2D Lives!

Box2D now has it own web site at box2d.org. The source is now hosted on SF and there are many new features and enhancements. There is even a forum at the new site! How modern.

I plan to continue supporting Box2D as I am addicted to physics programming. Please, somebody help me! :)

G++ Update and Games

I updated Box2D to support G++ and Linux. Please get it here and let me know if it compiles on Linux. This version is based on the GDC ‘07 version.

The fine folks over at Kloonigames have found Box2D useful and made a couple games with it. To mark the occasion, I have created a games page. Please let me know about any other games using Box2D.

Also, I have updated the theme once again to improve readability. Hopefully there will be no injuries this time.

GDC 2007

Hello all!

I again participated at the GDC, giving a physics tutorial with Jim Van Verth, Christer Ericson, Squirrel Eiserloh, Gino van den Bergen, and Marq Singer. I changed my presentation a bit this year because I wanted to coordinate better with the other speakers and I wanted to discuss Jacobians and more constraint types. I have added only a couple of enhancements to Box2D at this time.

Download my presentation here. I hope you find it useful.

Sequential Impulses

In 2005 I had a poster session at the GDC called “Iterative Dynamics with Temporal Coherence.” Obviously, such an obscure title is not likely to gain the headlines at the GDC. Further, trying to describe the PGS algorithm to game developers proved to be quite challenging. Nevertheless, a couple people out there managed to struggle through it.

I spent much of the next year talking with my co-worker and friend Gary Snethen about ways to simplify the PGS algorithm into something digestible for a wider audience. Gary had the brilliant idea that PGS is equivalent to applying impulses. I also read this and a paper from ETHZ 2005 (I can’t find it online now).

So I set out to prove that PGS is equivalent to sequential impulses. I scribbled some formulas down and it immediately became apparent that they are equivalent if you clamp accumulated impulses. So that was the birth of my version of SI. I presented all of this at the GDC 2006 tutorial with the Box2D mini-engine.

Now I don’t believe that SI is necessarily a new invention, however I had not seen the algorithm explicitly presented in any published work. I believe that accumulated clamping is the key to the algorithm, and I have not seen accumulated clamping presented explicitly in any other work. Also, SI with accumulated clamping is mathematically equivalent to PGS, so it is not a new mathematical device, on the other hand it is a new algorithm in the sense that a separate accumulator is not needed.

I hope this settles the origins of SI once and for all and properly credits Gary for his contribution.

The gPhysics Engine

Back in 1999 I wrote a physics engine using joint coordinates (Featherstone) and a direct LCP solver (Lemke). By no coincidence, the engine was called gPhysics.

The Lemke solver sucks. It is slow and has numerical problems. However, I am now again interested in joint coordinates. So I dusted off my old code and came across these demos.

The gPhysics engine proves that you can extend Featherstone’s algorithm to simulate trees, detached bodies, and contact. Now I just need to hook it up to a better contact solver.

Split Impulses and Post Projection

I’ve been tinkering with Box2D.

The original version uses Baumgarte stabilization of position errors. The problem with Baumgarte is that it affects the momentum, leading to bouncy or spongy constraints. I’ve attempted to address that problem with Split Impulses and Post Projection

Split Impulses works by making a separate set of pseudo velocities and impulses that serve to reduce position error and don’t influence the velocity solver. This is a fairly cheap solution because it is mixed in with the velocity solver. It works great for contact constraints. However, it does not work so well for joints. I don’t know why it doesn’t work well for joints, but on the other hand, I don’t know why it works so well for contact.

Post Projection is a more mathematically correct method of dealing with position errors. Post Projection (PP) waits until the velocity solver is done and the positions are integrated. PP then attempts to solve the position errors. PP introduces an outer iteration loop which recomputes the Jacobians and updates the positions. PP is the most expensive solution so far, but it is also the most effective.

See this post on Erwin’s site for more info: Bullet Forum Discussion

I also added correct handling of joint softness. This allow me to simulate a suspension bridge with soft joints and gets rid of the relaxation hack.

Box2D Super Split