Real-Time Physics Simulation Forum

I revised my physics tutorial for this year's GDC. You can download it here:

http://www.gphysics.com/downloads

Eventually, you should be able to find the other sections of the tutorial at:

http://www.essentialmath.com

Erin Catto wrote:I revised my physics tutorial for this year's GDC.

Thanks for posting this! You have a great way of explaining a complex subject in an accessible style.

On the page called "A Recipe for J", you write that one should steer clear of the partial derivatives. Could you explain this?

matti

Normally in textbooks they ask you to compute Jacobians by computing partial derivatives. If you sit down and try to do this, you'll find that it is quite tedious. On the other hand, you get Jacobian in a simpler fashion by computing the time derivative of the position constraint and isolating the velocity terms.

So there are two methods to compute the Jacobian and I'm suggesting that one method is easier than another.

Erin, apart from those 2 methods I think there is a 3rd popular method to compute Jacobian:

Empirical: measuring the effect of applying a unit-impulse. It is mentioned in Bart Barenbrug PhD thesis and Dynamo library.

I've seen some people using this in combination with Featurestone-style constraint solving.

Have you tried that?

Erwin, in my tutorial I show that the Jacobian maps constraint impulses from constraint space into Cartesian space. In order to apply the test impulse (lambda) you must convert it to Cartesian space (P). This implies that you already know the Cartesian space Jacobian.

Now when we consider a joint coordinate formulation (Featherstone), we need the joint space Jacobian. Actually, what we are after is the effective mass:

effective_mass = inv(J * invM * JT)

Given the current velocity error (Cdot), the effective mass tells us how to compute lambda to send the velocity error to zero:

effective_mass * delta_Cdot = lambda

Compare the above formula to:

mass * acceleration = force

The two formulas are completely analogous, the just operate in different spaces.

To compute the effective mass we need to apply lambda in constraint space, then compute the resulting acceleration in Cartesian space, then project the result back into constraint space.

In detail, we provide a test lambda (lambda_test), then use our Cartesian space Jacobian to compute P. When then apply P to the body tree and use Featherstone's algorithm to compute the change in the joint velocities. From the change in the joint velocities we can get the change in the body velocities. We then apply the Cartesian space Jacobian to the body velocity delta to get the constraint space velocity delta (delta_Cdot_test). Then we can determine the effective mass:

effective_mass = lambda_test / delta_Cdot_test

This is a fairly involved process, but it only takes order N time. I have implemented this in a joint coordinate version of Box2D. I will release the source for this soon.

I think Barenbrug's dissertation is pushing the idea of the effective mass, but he operations with position constraints, so he computes:

effective_mass = lambda_test / delta_C_test

I'm not yet willing to operate at the position level for rigid bodies because impulses are the natural language for dealing with collisions and friction, and you still can get a reaction force (joint breaking, etc).

It is a terminology issue, on page 38 of his PhD, Bart Barenbrug defines the Jacobian as:

Jacobian dC/dR as 'dependency between constraint error and restriction value (impulse/force etc).

Indeed, it is a unit-impulse method to measure the effective mass, given you know the axis/direction in the right space.

Hi,

I have found the slides and code very useful it understanding rigid body systems. I have been looking at the joint implementation in the Box2D demo and am unsure how it relates to the formulas from the slides. I cannot work what method the joint is using to calculate the constraint impulse and was wondering if someone could point me in the right direction?

Cheers,
Mark.

The formula goes like this:

1) Position constraint (x := position, R := Orientation matrix, r' vector from COM to anchor in body space)
C(x1, R1, x2, R2) = x2 + R2 * r2' - x1 - R1 * r1' = x2 + r2 - x1 - r1

2) Velocity constraint
dC/dt = J * v = v2 + w2 x r2 - v1 - w1 x r1

3) Jacobian (through inspection)
J = [ -I3 | skew(r1) | I3 | -skew(r2) ]



Now the algorithm:

a) We seek an impulse P...
v' = v + W * P

b) ...such that the post-velocities satisfy the velocity constraint...
J * v' = 0

c) ...also knowing the impulse direction
P = JT * lambda


Finally you get:
J * W * JT * lambda = - J * v -> lambda = -J * v / ( J * W * JT )

Not that the effective mass matrix J * W * JT can vary between 1 x 1 and 6 x 6 for binary constraints depending on the removed number of degrees of freedom. Finally plugging in the Baumgarte term (stabilization) yields:

lambda = -( 0.1* C / dt + J * v ) / ( J * W * JT )


HTH,
-Dirk

Thanks for that Dirk that clears up a lot. I was wondering what I3 is? I know it is probably going to be a diagonal matrix related to the mass of the object but I'm not sure on the exact components?

Cheers,
Mark.

If I am reading it correctly - I believe I3 is the 3x3 identity matrix.

--EK

Of course, I was thinking of W when I replied , speaking of which leads me to another stupid question, in W do we use the local or world inverse inertia matrix for each of the bodies?

Cheers,
Mark.

W contains the global inertia tensor (a 3x3 matrix). So the W matrix is block diagonal. It contains four 3x3 matrices on its diagonal. Note that you usually don't build this matrix in memory. If you multiply J * W * JT by hand or e.g. Matlab you end up with an easy formula for your effective mass. (I do it by hand since it is really good practise).

One hint:
Often the result is something like this (usually four additative terms):
K = n1 * w1 * n1 + n2 * w2 * n2 + ...

Since length(n1) == 1 you end up with K = w1 + w2 + ...

A good exercise is to build the effective mass for a non-penetration constraint and then end up with the formula in Baraff's Siggraph course. After this I would do the same for a spherical joint (ball-socket). The result must be the same as the collision matrix K in Mirtich's PhD. Try it and come back if you have problems and I try to help.

Note:
w1 == inv(mass1)
n1 is some axis in world space, e.g. the collision normal or a hinge axis


HTH,
-Dirk