Impulse based simulation - Jacobson, Baltman, Guendelman (Jiggle)

This topic is 4498 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

Recommended Posts

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

Share on other sites
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.

Share on other sites
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]

Share on other sites
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....

Share on other sites
Great! Thanks for taking the time...

Share on other sites
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.

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:

ra = Rat(p-xa)
rb = 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]

Share on other sites
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!!

Share on other sites
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

Share on other sites
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]

Share on other sites
Quote:
 Original post by DonDickieDyes 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.

Share on other sites
Hey Ury,

before I answer your questions, here is a new idea how to solve the problem. I thought about the problem and came up with the following idea. Here we go:

The velocity update for the two bodies is:

M * du /dt = f_ext + JT * lambda

Here u is the generalized velocity vector for body1 and body2, JT is the transposed jacobian and lambda is a vector of langrange multipliers.

Since we use linear dynamics we can use superposition and apply the external and constraint forces subsequently. So the velocity update because of the the constraint forces becomes:

M * du /dt = JT * lambda

The compatibility equation (velocity constraint) for the two bodies is:

J * u(t) = 0

Numerical integration of the Newton-Euler equation:

M * ( u(t+dt) - u(t) ) / dt = JT * lambda

<=> u(t+dt) = u(t) + M^-1*JT*lambda*dt

Semi-implicit integration, therefore

J * u(t+dt) = 0

Plug integrated Euler-Newton into compatibility equation:

J * [ u(t) + M^-1*JT*lambda*dt ] = 0
<=> J*M^-1*JT * lambda * dt + J*v(t) = 0

This is a linear system: A*x + b = 0

If we look now sharply we see that J*v(t) is the constraint error and lambda * dt is the wanted impulse. We define

A := J*M^-1*JT
x := lambda * dt
b := J*v(t)

=> lambda = A^-1 * constraint error (J*v(t))

So the wanted linear and angular impulses for body1 and body2 are simply:

P = JT * lambda

What do you think. Actually here we have your Langrange multipliers :-)

Edit: I just tried this formulation and it nicely blends over into your equation for the ball socket. Your K matrix is actually J * M^-1 * JT and and J*v(t) is the constraint error. I think I will give it a try. Hopefully it works for the the other joints as well...

[Edited by - DonDickieD on October 23, 2005 9:51:17 AM]

Share on other sites
Quote:
 Edit: I just tried this formulation and it nicely blends over into your equation for the ball socket. Your K matrix is actually J * M^-1 * JT and and J*v(t) is the constraint error. I think I will give it a try. Hopefully it works for the the other joints as well...

That's right. :)
Compared to my version: ur = J*v(t).

Quote:
 So the wanted linear and angular impulses for body1 and body2 are simply:P = JT * lambda

Careful! Using your formulation, the impulse is given by:
P = JT * lambda * dt

Anyways, the problem you are solving is actually static. There's no need for dt at all. Your original formula:
Quote:
 M * ( u(t+dt) - u(t) ) / dt = JT * lambda

Can be written:
M * ( u(t+dt) - u(t) ) = JT * gamma
where:
gamma = lambda * dt.

Now, after you find gamma, you get:
P = JT * gamma

See? No dynamics at all! Looks more elegant, right? Thats the problem. :)
Neither the system forces nor the time-step are considered.
If you follow this scheme, there are two steps that you have to undertake:
1. Integrate your velocities in time using the system forces.
2. Apply the constraints as a post-correction to your new velocity.

Now consider this:
( M U )( v1 ) = ( Fdt + Mv0 )
( L 0 )( lambda ) ( E )

where:
M - system mass matrix
L - J
U - -Jt
F - force
E - constraint error

Solving this, both steps are performed simultaneously!
Here we have "my" lagrange multipliers ;P

Oh.. in your formulation you don't handle distance drift. Don't forget it :)

A silly question...
f_ext includes coriolis forces as well, right?

[Edited by - ury on October 23, 2005 11:28:27 AM]

Share on other sites
Very cool! I have never thought of lagrange multiliers in this context. I always thought of them when thinking of the velocity constrained approach like in the ODE.

Regarding f_ext:

Yes, it includes the corioles force. Actually the system behaves better without it. Is seems as if long and thin objects with a high angular velocity gain energy when I consider this term. Any suggestions?

Little confusion - in your formula lambda is actually an impulse, right? Just to be sure.

In this form you solve a linear system. How would you solve limits in this context? I can think of two possible solutions:

1) Use an LCP and limit the lambda
2) Model the constraint error such that a limit is taken into account as well

Regarding contacts:

And now every physic programmers best friend - unilateral constraints :) I am currently implementing the Guendelman approach for modeling contacts. For integrating the constraints I worked with the sketch of R. Weinstein at the Siggraph this year. (If you are interessted I can send you the sketch plus some private correspondance where the idea is explained a little bit more).

1) Do see any problems to combine both apporaches?
2) How would you model contacts in a simulator where we deal with bilateral constraints as discussed in the previous posts?

Share on other sites
Quote:
 Yes, it includes the corioles force. Actually the system behaves better without it. Is seems as if long and thin objects with a high angular velocity gain energy when I consider this term. Any suggestions?

The system might be more stable but less realistic.
Try to remeber what happens if you spin a coin on a table. After loosing some of its speed, the coin falls on the table and starts "wobbling". This "wobbling" happens because of the coriolis force.

To make your simulation more stable you can:
1. Apply angular damping force (torque)
2. Clamp velocities.
3. Look for a bug :)

Check if smaller timesteps improve the simulation quality. If they don't, it must be a bug.

Quote:
 Regarding your formula:Little confusion - in your formula lambda is actually an impulse, right? Just to be sure.

No. Let me rewrite the equations:
1) Mv1 = Fdt + Mv0 + Jtlambda
2) Jv1 = E

You can clearly see that the correcting impulse is:
P = Jtlambda.

Let's assume that the valid momentum may lie in some hyperplane and let the rows of J be the normals that form this hypersurface. When we find the solution, we actually project the current momentum to this hyperplane.
The projection is done by using the normals from J which are scaled with the values from Lambda.

Quote:
 In this form you solve a linear system. How would you solve limits in this context? I can think of two possible solutions:1) Use an LCP and limit the lambda 2) Model the constraint error such that a limit is taken into account as well

Yes. LCP is used to formulate such constraints.

Quote:
 1) Do see any problems to combine both apporaches?

Quote:
 2) How would you model contacts in a simulator where we deal with bilateral constraints as discussed in the previous posts?

Usually, a MCP (Mixed Complementary Problem) is formulated.
Take a look at this thesis, it'll give you a general idea.

[Edited by - ury on October 24, 2005 4:31:23 AM]

Share on other sites
Here's a link to an excellent paper that describes an iterative method for solving lagrange multipliers:

http://www.gphysics.com/files/IterativeDynamics.pdf

It's the first paper on global solvers that I've been able to implement and from my experience, the results are more stable (no jiggle) and faster than the impulse methods I've previously tried.

-Steve.

Share on other sites
Ury:

Regarding Corioles:
Cool - now I have physical interpretation of this term. Thanks!

Regarding Lagrange Multipliers:
I once wrote down a small paper regarding this topic. I wrote it before reading Erin's paper (thanks sbroumley). You might find it interessting.

www.dirkgregorius.de/ConstrainedDynamics.doc

Regarding the sketch:

http://graphics.stanford.edu/~fedkiw/papers/stanford2005-05.pdf

Regarding physical simulators:
When using the velocity constrained approach like described in the PhD thesis by Kenny Erleben or in Erin's paper (and actually implemented in the ODE) you typically have the following simulation loop:

DetectCollisions();
SolveLCP();
DoTimeStep();
Render();

Basically you add the contact constraints and external forces and send it to your solver and ask him:

"Please find a velocity update regarding the current external forces and geometric constraints such that non of the geometric constraints are violated"

"Please lightning fast and with the most reasonable solution in a degenarated case"

The problem with this approach in my opinion is with that you might render invalid configurations all over the place. E.g. when the AI applies a big external force this will lead to a big displacement after the integration, which might end up in deep penetrations. We have this all over the place, besides convergence problems, etc...

What I do now is solve the system in the future (s.Guendelman) and then integrate the phase state such that no constraints are violated in the next frame in order to avoid the described artifacts. To do this I advance the system as if there were no constraints and then apply a post correction as you correctly stated above. Basically this is the Jacobsen projection, but on the velocity level and regarding rigid bodies. While a lot of people find solving the big system matrix more elegant in my opinion this modular approach will give me much more possibilities to correct something in a reasonable way when running into a singularity. Actually this is the same relaxation approach as a Jacoby or Gauss-Seidel iterative solver applies - at least for my understanding. Though on the mathematically level the big matrix approach solves a DAE while I solve a differential inclusion. Another advantage of this approach is that I can now easily use high order integrators, while this is quite difficult in the velocity constrained formulation. Since the constraints are not considered in the integration, there is no need to find the derivatives. Any comments, critics and ideas are appreciated...

Regards,

-Dirk

Share on other sites
Quote:
 Original post by DonDickieDWhat I do now is solve the system in the future (s.Guendelman) and then integrate the phase state such that no constraints are violated in the next frame in order to avoid the described artifacts. To do this I advance the system as if there were no constraints and then apply a post correction as you correctly stated above. Basically this is the Jacobsen projection, but on the velocity level and regarding rigid bodies. While a lot of people find solving the big system matrix more elegant in my opinion this modular approach will give me much more possibilities to correct something in a reasonable way when running into a singularity. Actually this is the same relaxation approach as a Jacoby or Gauss-Seidel iterative solver applies - at least for my understanding. Though on the mathematically level the big matrix approach solves a DAE while I solve a differential inclusion. Another advantage of this approach is that I can now easily use high order integrators, while this is quite difficult in the velocity constrained formulation. Since the constraints are not considered in the integration, there is no need to find the derivatives. Any comments, critics and ideas are appreciated...

Sounds kind of familiar. Keep going.

Share on other sites
Quote:
 Regarding Lagrange Multipliers:I once wrote down a small paper regarding this topic. I wrote it before reading Erin's paper (thanks sbroumley). You might find it interessting. www.dirkgregorius.de/ConstrainedDynamics.doc

I'll read it as soon as I get the chance. :)

I have mixed feelings about Erin's paper. The paper is very clear and straight forward. The notations are very simple and the paper can be a VERY useful reference for a newbie since it describes the entire simulator from head to toe.
Unfortunately, the collision and friction models proposed in this paper are very poor. I have to credit Erin for not hiding this fact and openly stating this weakness. The most disappointing part in this paper is the iterative solver. He identifies his solver as Gauss-Seidel method. In reality, his solver is a close relative of the Jacobi's method with a bug. :) Jacobi's method on it's own has a very poor convergence rate. I am sure that the bug doesn't improve that either.

EDIT: Actually, there's no bug in Erin's solver. See Erin's explanation here.

For a good implementation of an iterative solver (SOR) and friction, please read Kenny Erleben's thesis. An implementation of this solver can be found in QuickStep which Kenny wrote for the ODE engine. BTW, there's even a more efficient friction method than the one proposed by Kenny. This method is formulated using box constrained LCP.
Tell me if you want to hear more about any of these subjects.

Quote:
 Regarding the sketch:http://graphics.stanford.edu/~fedkiw/papers/stanford2005-05.pdf

Seems like a nice idea. In fact, it's a very important thing to implement if you plan using Guendelman's approach.

Quote:
 Regarding physical simulators:.The problem with this approach in my opinion is with that you might render invalid configurations all over the place..What I do now is solve the system in the future (s.Guendelman).

For me, both approaches are almost the same.

The difference between:
1. Find new velocity.
2. Apply post-correction to confirm with a set of constraints.

and:
1. Find new velocity which confirms with a set of constraints.

is that the second approach can be more accurate.

If no external forces are applied, solving LCP is equivalent to performing a post-correction of velocities. In this case, you consider a much larger set of constraints, hence, the solution is more accurate. (Our last two posts are a good example for that). In fact, nobody's stopping you from using both aproaches at the same time. For instance, you could solve collisions using an impulse based approach and handle joints and contacts using LCP.

BTW, in my current engine, I do the following:

1. Collision detection using a continuous sweep.
2. Perform a full step.
3. Collision detection using a continuous sweep.
4. If I don't like what I see:
4.1 Step back in time.
4.2 Perform a full step taking pairs from 1 and 3 under consideration.

Look how closely it resembles Guendelman's approach:
a. Collision detection and modeling.
c. Contact resolution.

If collisions occur, then steps 1,2 can be treated as a, b and steps 3,4 as c, d.

Quote:
 Another advantage of this approach is that I can now easily use high order integrators, while this is quite difficult in the velocity constrained formulation. Since the constraints are not considered in the integration, there is no need to find the derivatives.

I don't follow you on this one :)
It doesn't matter which approach you choose, at the end of the day you have a black box that gives you a new velocity.

Just pick your favorite integrator and plug this velocity in. :)
In fact, you can even get even kinkier than that. In my engine, since I have a very stiff equation set, I am forced to use an implicit integration scheme. I do that, by finding a root for the equations of motion:
x = x0 + dt*v( x )

[Edited by - ury on October 26, 2005 5:55:03 AM]

Share on other sites
Quote:
 Tell me if you want to hear more about any of these subjects

Yes please! Especially, I like to know where the bug in the projected Gauss-Seidel is. The friction method is very interssting as well. Since you mention boxed LCP - I am still looking how to transform the classic LCP into a boxed LCP. Any reference?

Novodex uses the Gauss-Seidel as well (s. Lecture 6)

Technically you are absolutely right, there is very few difference. What might be an issue is that you need to call the collision engine twice (or even more), while I get away with calling it only ones. But I am not so sure here :-/

BTW: How do you generate contact points? Like in the ODE and OpenTissue or liked proposed in Guendelman using signed distance maps. If you have any good references they would be highly appreciated...

Regarding higher integrators:
The Newton-Euler equation in a constrained systems has the form:

M * du/dt = f_ext + f_c

When you want to use a higher integrator you need the forces at some sampling points, e.g. for the MidPoint Method at t+dt/2. I haven't looked into this problem so far, but I can't directly see how to compute f_c(t+dt/2)? Since I don't consider them in the first place, I have no need to compute them.

Share on other sites
ury. I'm sorry.. But you are a math punk. The Guendelman paper is the 3rd most important paper in physics simulation. Sure his collision stuff is craziness, but he was also handling mesh-mesh. That paper is important because it finally showed the academic snobs that PLAUSIBLE is the key.

Others have hinted at it in random papers. Even Baraff himself. But Guendelman came out and said.. "Damn it. It's just not that hard." If that paper was out before I had "figured things out" I would of saved myself a lot of time. The best thing about that paper is it makes physics APPROACHABLE. That in itself allows others to bring great inovation to the problem in the future. Unforunately most academics litter their papers with unecessary complexity just so that they can look smart or impress their famous advisors.

btw. Just because an algorithm has a faster convergence rate on paper does not make the algorithm take less wall clock time to execute on a computer.

Share on other sites
Quote:
 Unforunately most academics litter their papers with unecessary complexity just so that they can look smart or impress their famous advisors

I think I have read every paper and master or PhD thesis regarding physically-based animation. It is sometimes amazing how even the most simple problems are turned into really sophisticated ones. So I can't agree more :-)

Anyway, Ury has a very sorrow understanding of both mathematics and dynamics. And I am glad that he shares his konowledge and has the patience to answer all my silly questions. There is never no such thing like one way to go. I think there are some descent simulators out there that apply his technique, while there are others that use maybe another approach. This was a very good discussion so far and it would be sad if we end up in some kind of flame war here.

Regards,

-Dirk

Share on other sites
Quote:
 Original post by billy_zelsnackury. I'm sorry.. But you are a math punk. The Guendelman paper is the 3rd most important paper in physics simulation. Sure his collision stuff is craziness, but he was also handling mesh-mesh. That paper is important because it finally showed the academic snobs that PLAUSIBLE is the key.

I cannot disagree more with this statement. It would not call the paper “the more important…”, I would not even classify it as a paper worth reading, the only thing the Guendelman demonstrates is that it is a good example of what not to do.

Share on other sites
Quote:
 Original post by uryThe most disappointing part in this paper is the iterative solver. He identifies his solver as Gauss-Seidel method. In reality, his solver is a close relative of the Jacobi's method with a bug. :) Jacobi's method on it's own has a very poor convergence rate. I am sure that the bug doesn't improve that either.I guess that this is one of the reasons why he needed contact caching. :)

Ury, so I'm really curious now - what is the bug?!

thanks
-Steve.

Share on other sites
No flame war at all intended. I agree. This has been a very nice thread and one I sure others will reference in the future. Ury has presented some very good information. I just wanted to chime up on the stacking paper.

Share on other sites
Quote:
 I cannot disagree more with this statement. It would not call the paper “the more important…”, I would not even classify it as a paper worth reading, the only thing the Guendelman demonstrates is that it is a good example of what not to do

Could you give reasons for this? Without justification such a posting makes very few sense. I mean what shall one think after such a post. You could be either a very experienced physic programmer in which case it would be interessting to hear the justification for this statement or you are simply not able to implement the approach described in the paper...