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
GDC Tutorial '07
-
- Posts: 1
- Joined: Thu Apr 06, 2006 4:14 pm
- Location: Stockholm, Sweden
Re: GDC Tutorial '07
Thanks for posting this! You have a great way of explaining a complex subject in an accessible style.Erin Catto wrote:I revised my physics tutorial for this year's GDC.
On the page called "A Recipe for J", you write that one should steer clear of the partial derivatives. Could you explain this?
matti
-
- Posts: 316
- Joined: Fri Jul 01, 2005 5:29 am
- Location: Irvine
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.
So there are two methods to compute the Jacobian and I'm suggesting that one method is easier than another.
-
- Site Admin
- Posts: 4221
- Joined: Sun Jun 26, 2005 6:43 pm
- Location: California, USA
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?
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?
-
- Posts: 316
- Joined: Fri Jul 01, 2005 5:29 am
- Location: Irvine
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).
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).
-
- Site Admin
- Posts: 4221
- Joined: Sun Jun 26, 2005 6:43 pm
- Location: California, USA
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.
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.
-
- Posts: 7
- Joined: Mon Oct 02, 2006 12:37 am
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.
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.
-
- Posts: 861
- Joined: Sun Jul 03, 2005 4:06 pm
- Location: Kirkland, WA
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
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
-
- Posts: 7
- Joined: Mon Oct 02, 2006 12:37 am
-
- Posts: 44
- Joined: Sun Jan 22, 2006 4:31 am
-
- Posts: 7
- Joined: Mon Oct 02, 2006 12:37 am
-
- Posts: 861
- Joined: Sun Jul 03, 2005 4:06 pm
- Location: Kirkland, WA
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
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