yet another physics engine tentative

[zoso]

Hello everyone,

I've been working for some time at a physics engine as a hobby. Recently I got more involved as I decided to make it my degree paper. Unfortunately I kind of got into a dead end, nothing I tried seemed to work so I finally decided to get a holiday. Now that I got back after two weeks, I want to get back to work as I have a presentation in two weeks, so I decided to ask for help.

I've read part of this forum, but it's quite large to view it all. My problem is stability too and then performance. I'll describe you briefly my engine and afterwards you can tell me what could be wrong or what could be improved.

I use a time-stepping scheme/velocity based formulation as described in Catto and Erleben: symplectic Euler integrator, a PGS solver with SOR and a Guendelman approach to contact and collision resolution. Although these should be some of the best approaches that I've read about, I get a lot of instability especially in my bilateral constraints. I'll give you an outline of the engine to be more precise:

1. relax joints
2. detect collision
3. update velocities
4. solve contact/constraints
5. update positions

1 means something similar to Jakobsen method in order to avoid large Baumgarte factors, then in 2 I detect pentration mostly between OBBs (althoug I have spheres and capped cylinders too) and get contact info in pairs of closest points. I add them to a list and then for each of them (one by one) I apply collision impulses like in the Baraff-Witkin tutorial (for pairs of objects) and do a lot of projection to decrease inter-penetration. After 3 I solve the mixed LCP that I've built using the list of contacts and the joints. I have only dynamic friction, as described in some Baraff or Anitescu papers, so I don't use any friction cone approximation or tangential directions or constraint forces. After the PGS I clear the contact list and go to 5 and back to 1 looping of course , recalculating the contact list everytime, as well as the jacobian, constraint function and the A matrix.

I use a time step of 0.001, otherwise it blows up, but I do about 20 updates per frame, although I'm still far from 60hz. Gravity is -9, the restitution factor in the impulse-based part is 0.02, beta for joints and contact is 0.1, 15 iterations for PGS and w=1.2 for SOR, and I also relax the lambdas by 0.8. For now I only have ball-socket joints and I built an articulated rag-doll with OBBs and a sphere for the head (4 limbs each with 2 parts). It gets very unstable at violent impacts and when falling on other boxes. It gets a lot of extra energy and angular velocity (could be because I don't use hinges?) and the system gets into rest only in some particular situations.

Please help me! Hints, pieces of advice, tell me what do I wrong, possible mistakes, future improvements, anything. I can post more info i fyou need, and even demos or source code (although I'm nor very proud of them ).

Thank you.

[Erin Catto]

I wouldn't bother with SOR. At most it might save you an iteration and you will have added energy. Try setting the restitution to 0. You shouldn't need to relax the lambdas. The only time I needed to relax lambdas was for a suspension bridge (now I just use softness).

Are you limbs colliding at the joints? You would have contacts fighting with the joints.

Try skipping your step 1. Just use Baumgarte. Make a variable length suspended chain and throw objects at it. Is it stable?

[zoso]

Thank you for your fast reply.

I've tried doing all the things you mentioned, and indeed I had collisions everywhere, including the joints. I'm not very glad about it, as the limbs could pass through one another in certain situations, especially at high speeds. Now the rag-doll seems to get much often and quicker into a stable position. But still when falling violently on a box the balls jump out of their sockets and get back quite slow. When trying to increase beta I get more energy and angular velocities, even though it slowly reaches stability.

I'm currently working on the chain test, although I'm not yet sure I've build it right. Otherwise it seems much stable, but starts to suffer the same problems as above for higher betas and speeds.

I'll get back to you after I do some more tests.

[zoso]

Hello again,

It seems that the chain is very unstable, especially when it catches angular effect after being hit somewhere in the middle. I don't understand anymore. It doesn't work even if I put contact and friction away and solve the system with LU decomposition. If I increase Baumgarte then I get quickly a lot of fast spinning. Decreasing gravity or the time step doesn't work as the system will still blow up, but after a longer time.

I must be doing something terribly wrong!

[Erin Catto]

It seems like you might have a problem with your joint Jacobian or a bug in your joint solver. Try to simulate a pendulum and compare the result to the numerical integration of:

x' = v
v' = -(g/L) sin(x)

You should see better agreement as the time step decreases and/or the iterations increase.

[zoso]

Thank you very much. I was so dumb. You were right, there was a problem in the jacobian - one index. I'm so ashamed I didn't see it from the start.

Now it is so much stable. It is stable even without contact and friction constraints, just impulses. Even more stable it seems. And the joints remain glued even if I put back the collision between the limbs or restrict it to a certain distance from the joint.

Curiously, it still becomes rather unstable in some situations (less often), especially when the rag doll falls when 2 stacked boxes. Or when I try to increase the time step. It takes much longer time than before, but eventually it gets that odd spining and energy gain. And there is the penetration too. I no longer use projection as it seems to be worse.

Could you explain me why there is the need to solve the LCP for contact and friction when we have the impulse handling? Why is it better? Should there be both? And why did you advise me in the first place to set the restitution to 0?

I use friction both in collisions and contact, but I can't use a too large coefficient for the dynamic friction because it will blow up - the LCP hasn't got a solution I think. I tried using only the impulsive and it was quite OK. How should I approach this? Should I try to implement the friction like in your paper?

Thanks again.

[Erin Catto]

It sounds like you have a separate solver for collision and contact. I prefer to work with a just a contact solver that adds a velocity bias for bounce. Solving impulses sequentially with accumulated clamping is equivalent to solving an LCP. See if this helps.

[zoso]

I've rewritten the code, mainly after your paper "Iterative Dynamics with Temporal Coherence", except for the contact caching final part and it seems to run quite well at a time step of 0.01 (100Hz). As for the contact and collision part I decided to make two variants: one with impulses and one with contact forces. And they both work, with the first one appearing to be a bit better: both fast and stable.

Now, I've re-read the Guendelman paper and there is a clear distinction between collision (impulses) and contact (forces). As I understood, Baraff used a threshold velocity for separating colliding contact from resting contact, and here is where the Guendelman approach comes handy.

I've also read your latest slides and it seems you've change your method. These sequential impulses are another way of doing the PGS? Are they computed for pairs of two objects? (as I do with my collision impulses) In that case isn't this equivalent to zero-restitution micro-collisions applied in an iterative fashion?

So what do you think? Should I stick to impulses or contact? Should I bother to use the Guendelman scheme when I have just one contactat solver? (I noticed it helped decreasing Baumgarte for contact, but I think it turns the time-stepping method into an implicit one) Or is this new "sequential impulses" method better?

[Erin Catto]

SI can be mathematically equivalent to PGS. The key is to clamp accumulated impulses. This makes the method convergent. SI is naturally expressed on pairs of bodies, but it can easily be extended to handle n-way interactions.

The reason I prefer SI over PGS is that it is easier to understand and more flexible. However, if they are setup correctly, they will give mathematically identical results.

SI is not equivalent to Mirtich's micro-collisions because he did not clamp accumulated impulses. I have not tried Guendelman's algorithm, but if you combine contact+restitution into one solver you will get better performance.

[Dirk Gregorius]

The key difference between Erin's SI and Guendelman/Mirtich is that Erin clamps against the accumulated impulse while Guendelman and Mirtich require each delta impulse to be greater than zero. In Erin's presentation you see that this is bad for the convergence. Erin's method is equivalent to PGS so you can see it the other way around. While PGS/SI solve mathematically founded the mixed LCP there is no such guarentee for Mirtich/Guendelman. Personally I would say those methods are inferior to PGS/SI. In the particular case of Guendelman you see that he needs to play with the coefficients of restitution (< 0) during the contact phase (to achieve a better sampling) and also a shockstep. The Guendelman approach works nice for exactly one box on the ground. Through the shockstep he proceeds than through the layers reducing the problem to something he can handle - one box on the ground. When you extend the Guendelman approach to deal with joints (e.g. a ragdoll lying on the ground) you run into severe problems e.g. jittering.

In shouldn't make a difference whether you use impulses or forces since they are linearly related (F = P/dt) in your examples.

[zoso]

Hello again,

I'm sorry I haven't had the chance until now to thank you for your replies, but I've been busy. They were quite useful hints, as I finally got to the conclusion that the Guendelman scheme has no relevance with this simulator. Also, I finally mixed impulsive restitution with constraints using a threshold both for normal relative velocity and penetration depth. I had some problems with objects passing through one another at high speed and I had to use a bit of projection too. But at low speeds it's only PGS. Of course it could be solved through continuous collision detection, but I need to do a lot of improvement in this part.

Here is a demo with what I've done. I've implemented some kind of mouse manipulation, but it's still quite buggy. The engine is still slow: I have to use a step of 0.005 so 3 stacked boxes won't jitter. Most parameters can be modified from the input xml file. But I've noticed that it runs quite slow anyway on some computers. So it still needs a lot of optimization.

The friction is both dynamic and the one in Erin's Paper with the constant friction cone. I'll try to do the polyhedral approximation in the future. I've also tried to implement the hinge joint, but again it seems I'm doing something wrong - I don't understand if the rotation axis should be kept in world space, in one body's space or both.

The rag-doll is still unstable sometimes. I know it's because of the contacts at the joints too, fighting the bilateral constraints, but until I implement joint limits I wouldn't want to eliminate them. Maybe there are other reasons too. Another odd behavior is when I hang the rag-doll from a fixed point without touching anything and it preserves its momentum and keeps on spinning, even faster it seems, when it should have reached an inert state. Is this normal? Or it needs friction at the joint or with the air?

That's about all for now. Thanks again and I hope we can continue this discussion.