Given two bodies b1 and b2. In an impulse based simulator you integrate the bodies forward in time, find the constraint error and correct it. Mathematically I believe this are differential inclusions and you project the phase state of the bodies back onto the constraint manifold. This sounds difficult, but is extreme easy to implement. For the ball-socket joint you simple do something like this: b1.Step( dt ); b2.Step( dt ); // Let r1 and r2 be the offset vector from the CM the the joint anchor point v32 v1 = b1.GetVelocityAt( r1 ); v32 v2 = b2.GetVelocityAt( r2 ); // The relative velocity must be zero - so this is the constraints error v32 v_rel = v1 - v2; That's not too difficult. Actually this is the Jacoby entry, but you don't have to deal with large matrices and do all the bookkeeping. Now that we have found the constraint violation, we need to correct it - that is project back onto the constraint manifold. We don't want to simply divide the error between both bodies, but also consider the mass and moment of inertia somehow. We have to solve the following question: "Find corrections for the current linear velocities v1 and v2 plus corrections for the current angular velocities w1 and w1, such that the relative velocity at the joint anchor becomes zero. Preferable considering the masses and moments of inertia." In the literature (s. Papers) it is proposed to calculate the collision impulse and chose the direction of the relative velocity as collision normal and its length as relative normal velocity. E.g. in Jiggle the constraints are implemented this way. f32 j = -(1 + e)*vn / [ (1/m1) + (1/m2) + n*(J1^-1*(r1 x n) x r1) + n*(J2^-1*(r2 x n) x r2) Now the question: I assumed that after applying the impulse the relative velocity would be zero. This is sometimes correct, but in the majority of cases it seems as if the error is only reduced a bit. Isn't this bad for the convergence of the system? For two limits I don't have convergence anymore. Are there other methods to divide the error between the two bodies? I setup the impulse equations and avoided the calculation of the collison impulse - since this would be more numerically stable. But I think I have an error in the matrix, so I can't say how this works at the moment. Any suggestions? Regards, -Dirk

In an isolated system with a single joint I'm pretty sure that after applying the impulse my implementation gave zero velative velocity at the joint location.

Of course if the objects are part of other constraints (e.g. other joints, or resting on a surface under gravity etc) then applying impulses sequentially you only relax to the proper solution, but it should converge to a reasonable tolerance with a fairly small number of iterations, and any offset error (for normal conditions) can be accounted for by an error correction term. Without the error correction term you get some drift.

If I were you I'd check a single joint in isolation (no gravity) - if you don't almost get perfect behaviour there must be a problem in your implementation.

Of course you'll still get some drift when the objects can rotate, since even if the velocities at the joint position match perfectly after applying the impulse, when integrating the angular velocity over the timestep this will (generally) result in the joint position separating. Even if you solved all joints/constraints simultaneously you'd still get this source of drift.

You process all body pairs subsequently and apply 20 iterations. I believed that after you apply an impulse to the bodies, the constraint would be totally satisfied (at least for the moment). Of course subsequent constraints can again violate previous ones. This is why we do several iterations (relaxation).

Example:

Given two boxes b1 and b2.

Width1 = Width2 = 3
Height1 = Height2 = 2
Depth1 = Depth2 = 2

Mass1 = Mass = 12kg

From there I derived the moment of inertia:

Ixx = 8
Iyy = 13
Izz = 13

The initial orientations are identtity quaternion and the positions are

x1 = (0, 10, 0)
x2 = (0, 6, 0)

The ball joint is located exactly between the two bodies:

xj = (0, 8, 0)

Then the offset vectors from the centers of mass to the joint anchor are:

r1 = ( 0, -2, 0)
r2 = ( 0, 2, 0)

As initial phase state I chose:

v1( 0.0f, -13.0f, 0.0f )
v2( -2.0f, -3.0f, -2.0f )
w1( 0.0f, 0.0f, 0.0f )
w2( -1.0f, -1.0f, 2.0f )

The velocity constraint equation is:

v1 + w1 x r1 - ( v2 + w2 x r2 ) = 0

This is the relative velocity at the joint anchor. If I compute this from hand and in my system I get:

v_rel = (6, -10, 4)

Now I chose the direction of v_rel as collision normal and its length as the normal velocity:

v32 n = v_rel / |v_rel| = (0.48666424, -0.8111071, 0.32444283)
f32 vn = |v_rel| = 12.328828

Since we have identiy orientations the inverse global moment of inertia is:

Ixx^-1 = 1 / 8
Iyy^-1 = 1 / 13
Izz^-1 = 1 / 13

Finaly we apply a plastic (inealstic) collision, so we chose e = 0

Now we can enter the impulse formula:

f32 j = -(1 + e)*vn / [ (1/m1) + (1/m2) + n*(J1^-1*(r1 x n) x r1) + n*(J2^-1*(r2 x n) x r2)

I get j = -29.517485

Now I apply the ( correction ) impuls.

v1' = v1 + j * n / m1 = ( -1.1970921, -11.004847, -0.79806137 )
w1' = w1 + I1^-1 * ( r1 x j * n ) = ( 2.3941841, 0.00000000, -2.2100160 )
v2' = v2 - j * n / m2 = ( -0.80290794, -4.9951534, -1.2019386 )
w2' = w2 - I2^-1 * ( r2 x j * n ) = ( 1.3941841, -1.0000000, -0.21001601 )

If I recalculate the relative velocity here I get:

V_rel' = ( 5.2342482, -6.0096931, -7.1728592 )

This is the isolated example. The problem is that this works! The ball joint, the hinge with limits, the slider with limits and bounce, the universal, all work. But when I add a second axis limit to the universal joint there is no more convergence. If I run a constraint verification after each applied impulse the constraint error is nearly never repaired. It seems as if the collison impulses only repair a bit of the error, which is way too slow in my opinion.

BTW: My fist guess was that some sign was wrong, but I tried all four permutations for +/-vn and +/-n. None worked. Any help is greatly appreciated.

-Dirk

[Edited by - DonDickieD on October 22, 2005 4:07:30 PM]

OK - I'll try to find time to set up the same experiment in my code, but am not sure when I'll get a chance to do this....

Great! Thanks for taking the time...

Dirk, there's absolutely no reason why your approach should work at all.

Although the impulse that you generate at the end makes sure that your new relative velocity is orthogonal to your normal, it's hardly zero.
Even if that wasn't enough, you have to compensate for any distance drift that may occur between the anchors of the joint. Because of that, making sure that the relative velocity is zero is simply not enough.

Here's something that might help you:

First some terminology:
x - position of the center of mass
v - velocity of the center of mass
R - rotation matrix
w - angular velocity
IL - inertia matrix in body's local coordinates
I - inertia matrix in world coordinates.
M - mass matrix of a body (a diagonal matrix)

- cross matrix.

u is equivalent to p x u.

I = R * IL * Rt

Let's set up a ball socket joint between two bodies, A and B, at some point p.
First we record the anchor position in body's local coordinates:

r`a = Rat(p-xa)
r`b = Rbt(p-xb)

These vectors should stay constant for the entire simulation.

For every iteration, we define for each body:
r = R * r`,
K = M-1 - [r] * I-1 * [r]

Now we can define the distance between anchors in world space as:
dist = xa + ra - (xb + rb)

We want to make sure that dist stays zero all the time.

Velocity of an anchor point in world space is:
u = v + w x r

And so, the relative velocity of the anchor points is:
ur = ua - ub.

Let's say the we have some impulse J. If we apply J on body A and -J on body B, the new relative velocity of the anchor points is given by:
Kt = Ka + Kb
ur` = ur + KtJ

If we wanted our new relative velocity to be zero, then we could find the appropriate J by solving:
-ur = KtJ

Since we want to compensate for any distance drift, we can try and solve:
alpha = -e/(2*dt)
alpha * dist - ur = KtJ

where:
dt - the size of your time step.
e - some empirical coefficient. (A good value to start with is 1.0)

Good luck, tell me if it works.

PS. If I were you, I'd really consider using lagrange multipliers :)

[Edited by - ury on October 22, 2005 1:05:31 PM]

ury - yes of course you're right, it's pretty obvious now you've pointed it out :) Rats!

DDD - I did set up the same experiment, but I actually get slightly different numbers to you - perhaps because your moment of inertia calculation is wrong. If you have a block that has dimensions (large, small, small) the moments of inertia around the axes should be (small, large, large) - you have two smalls and one large(!), which is wrong whatever your choice of axes.

So, before the impulse the joint velocity error magnitude is 12.3 (same as you get).

My impulse is 32.69, and afterwards the joint velocity error is 9.5 and pretty well perfectly orthogonal to the impulse direction.

Guess I need to work on this a bit :)

Big apologies to anyone that looked at my code and assumed it was good!!

Ury:

Thanks for the suggestion. I will try it out this afternoon and post the result here. I haven't looked to deep into the solution, but I am quite sure that when the corrected velocities are calculated we have to consider the conservation of momentum. So I came up with this solution. I scanned my calculations and uploaded it here:

www.dirkgregorius.de/Kugelgelenk.jpg

MrRowl:

Thanks for taking the time to setup the experiment. Of course you are correct with the moment of inertia. This was only a copy error. In my test case this is correct, so strange that we don't get the same impulse.

Please don't excuse. There is no reason for doing this!


Note: I edited my previous post and corrected the inertia

Ury:

yes your solution works! Thanks for taking the time and post the solution. For the protocol - my solution (s. previous thread) works as well, though yours is more elegant and more efficient.

Regarding drift:

Currently I am handling the drift through projection - that is I simply translate and rotate the bodies into a valid position. I tried your method as well. Actually I am flipping back and forth between both of them. BTW: Is this the same as Baumgarte Stabilization, but in a slightly different context? What I would like to have is your solution and a post-stabilization like proposed by Baraff.

Regarding Lagrange Multipliers:

Why do you think that a constraint-based simulator is superior to an impulse based simulator? The constraint based systems handles all these scientific examples quite well, but for games I am not so sure. And what about scaleability on the next generation hardware. I started implementing a velocity constraint formulation, but I find the impulse based method much more preferable.

Here are some examples of my simulator. Please feel free to have a look. Simply drag and drop the *.phx files onto the framework.exe. Note that the system uses all over the wrong correction impulses, since I haven't changed it today. Under these circumstances the system behaves already quite well. I haven't implemented any auto-disabling nor shock-stepping, so the stacking doesn't behave so good right now. Since I don't have stacked objects at the moment this has lower priority. Constraints with limits and bounds, a car controller with suspension and ragdolls need be solved first.

LMB: Rotate
MW: Zoom
RMB: Apply force to body

www.dirkgregorius.de/Primo.rar

[Edited by - DonDickieD on October 23, 2005 11:42:29 AM]

Quote:Original post by DonDickieD
yes your solution works! Thanks for taking the time and post the solution. For the protocol - my solution (s. previous thread) works as well, though yours is more elegant and more efficient.

Glad to hear and thank you. :)

Quote:
BTW: Is this the same as Baumgarte Stabilization, but in a slightly different context?

Hardly... Although the formulae might look the same, the idea behind them is different.
My approach is a simple geometric observation. If we have some distance between our anchor points, defined by the vector 'dist', give the bodies enough velocity at the anchor points, so that they would each travel dist/2 towards each other during the next frame and meet (with gods help :)) at the same point.

Baumgarte enforced a constraint by introducing an "augmented" ODE for the constraint. Baumgarte stabilization behaves very much like a restrained harmonic oscillator and can be viewed as a spring connected between the two anchors. Tell me if you want to hear more details about it.

As you can see the ideas are extremely different but I guess that in the end, we get pretty much the same result. :)

Quote:
What I would like to have is your solution and a post-stabilization like proposed by Baraff.

What sort of post-stabilization are you talking about?

Quote:
Why do you think that a constraint-based simulator is superior to an impulse based simulator? The constraint based systems handles all these scientific examples quite well, but for games I am not so sure.

What do you call scientific examples? If you mean overly complex scenes with lots of interactions and constraints then you are looking at the future of gaming physics.

If you use velocity constraint formulation with lagrange multipliers (LM), both techniques find correcting impulses. It can even be argued that when solving a LM problem, an iterative solver is performing roughly the same operations as you would do in an impulse based (IB) simulator (although, the convergence of an iterative solver is usually better). I guess that if you look at it this way, both techniques are very close relatives.

But!
The solution for a LM problem is "optimal" because it is found considering all of your constraints at once (including joints, contacts, friction, etc)
You can expect to have less drift in your joints. This allows you to use more complex articulated bodies.
The technique can be generalized to support other types of bodies (such as deformable bodies).
I personally find LM more elegant and modular.

A downfall of LM is that it's linear but this problem can be remedied since you can always fallback to performing your corrections using "handmade" impulses just like in IB.

Quote:
And what about scaleability on the next generation hardware. I read that Novodex had/has quite some problems on the PS3 and I *assume* that they use either a velocity constraint formulation (Stewart/Trinkle and Anitescu) or an accelaration constraint formulation (Baraff/Wittkin).

Problems? That's strange. PS3 is supposed to have 7 dedicates processors just to perform SIMD floating-point operations. I think that you can simulate an atomic explosion in real-time with that monster. :)
As far as I can recall, Novodex really do use a velocity based LM system and solve it using some sort of an iterative solver.

Quote:
I started implementing a velocity constraint formulation, but I find the impulse based method much more preferable.

Why? What happened? :)

Quote:
Here are some examples of my simulator.

Good work, man! Keep us in loop as you improve your simulator.