Simultaneous Colliding Contact

[SnprBoB86]

I am currently implementing constrained rigid body dynamics in my game engine. I am using David Baraff's paper on constrained rigid body dynamics (http://www-2.cs.cmu.edu/~baraff/sigcourse/notesd2.pdf) as my guide for handling colliding contact, but have run into a problem.

Suppose there is a box in the virtual world. This box has an infinite mass (1/Mass = 0). Another box, with a mass of 1, is placed directly above that box and dropped straight down due to the force of gravity. Upon detecting collision, 4 contact points are determined, once at each corner.

That all works fine, but if I understand the laws of physics correctly, the upper box should bounce straight up without any rotation. Unfortunately, my simulation causes the box to rotate. It appears that the force applied by the last contact causes significantly more torque than any of the other contacts and hence the box rotates. If I reverse the order of these contacts, the box spins the same amount in the opposite direction. Order of contact should be insignificant.

My code is essentially the same as the “collision” and “find_all_collisions” functions presented in the paper linked above. I will post code snippets if the pseudo-code presented is not sufficient.

Am I correct in assuming the box should bounce straight up because the order of the contacts shouldn’t matter? If so, what could I be doing wrong, and how can I go about fixing it?

Thanks in advanced,
Brandon Bloom

[hh10k]

Rather than performing collisions contact by contact, I just calculate the forces from all the contacts and apply the average of those forces to the objects. It works well enough for me, but you may want to consider weighting the contacts differently (eg by penetration depth) if you want something a little more accurate.

[SnprBoB86]

Well the way the algo presented in the paper works is that colliding contact (non resting or diverging contact points) are repelled by an instantanious impulse.

When using an elasticty of 0%, the box dropped on top correctly stops instantly (although shows a tiny hint of rotation). With an elasticity of 100%, the box correctly bounces back to its original height, but has incorrect angular velocity.

I have an ApplyImpulse method which accepts the impulse vector as well as a position in body space. It calls ApplyLinearImpulse as well as ApplyAngularImpulse. Since the box jumps to its original height with 100% elasticity, the ApplyLinearImpulse call is clearly correct. Either ApplyAngularImpulse is flawed, or the value passed to ApplyImpulse is flawed.

Here is some of the code:


/// <summary>
/// Instantaneously changed the linear and angular velocity of this body by applying an impulse.
/// </summary>
/// <param name="impulse">Impulse to be applied.</param>
/// <param name="at">Position in body space where this impulse is to be applied.</param>
public void ApplyImpulse(Vector3 impulse, Vector3 at)
{
this.ApplyLinearImpulse(impulse);
this.ApplyAngularImpulse(Vector3.Cross(at, impulse));
}

/// <summary>
/// Instantaneously changed the angular velocity of this body by applying an impulse.
/// </summary>
/// <param name="impulse">Impulse to be applied.</param>
public void ApplyAngularImpulse(Vector3 impulse)
{ // Does this function even make sense? Should it be this.angularVelocity += Vector3.TransformCoordinate(torque, this.InertiaTensorInverse) * duration;
// Is InertiaTensorInverse body space? or world space? or does that not make sense at all?
this.angularVelocity += Vector3.TransformCoordinate(impulse, this.InertiaTensorInverse);
}

/// <summary>
/// Instantaneously changed the linear velocity of this body by applying an impulse.
/// </summary>
/// <param name="impulse">Impulse to be applied. Measured in Newton * Seconds</param>
public void ApplyLinearImpulse(Vector3 impulse)
{
this.linearVelocity += impulse * this.MassInverse;
}


Maybe someone could double check this functions for my sanity?

[SnprBoB86]

Can anyone comment if those functions look correct??

[Reedbeta]

Aren't you supposed to split up the impulse into its parallel and perpendicular components (relative to the axis from the point of application to the body's center of mass), and then apply the parallel component as the linear impulse and the perpendicular component as the angular impulse?

I would try writing your ApplyImpulse function like the following:

Vector3 parallel = at * (Vector3.Dot(at, impulse) / Vector3.Dot(at, at));
Vector3 perp = impulse - parallel;
this.ApplyLinearImpulse(parallel);
this.ApplyAngularImpulse(Vector3.Cross(at, perp));

I believe that's the right projection formula...parallel and perp should be the components of impulse parallel and perpendicular to the at vector.

[SnprBoB86]

I think you are mistaken, Reedbeta. If you look at this paper (http://www-2.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf), (section 5.5 on page 29), you will read:

Now suppose that the same force F is applied off-center to the body as shown in figure 11. Since the force acting on the body is the same, the acceleration of the center of mass is the same.

That would lead me to beleive that this.ApplyLinearImpulse(impulse); would be correct. Also, as I stated earlier, the box does bounce back up to the correct height (according to its elasticity). So then it is the angular impulse that is primarily in question here. I believe my code is correct, but I am not certain.

Thanks for all the help so far, still trying to sort this out.

[Reedbeta]

Okay, I checked the Baraff paper and you are correct; the force is indeed "considered twice" for linear momentum and angular momentum.

The only other thing I can think of is to make sure that your InertiaTensorInverse is corrected for the body's current orientation - that is, InertiaTensorInverse = RotationMatrix * InertiaBodyTensorInverse * Transpose(RotationMatrix) where the InertiaBodyTensorInverse is as described in section 2.10, pgs 14-15 of the Baraff paper. The InertiaBodyTensorInverse is a constant, but the InertiaTensorInverse has to be recalculated each time the body's orientation changes.

[SnprBoB86]

Sorry for the delayed response. I corrected the rotated body inertia tensor issue. Thank you for finding that for me, but because I was using the identity matrix as my body intertia tensor it was clear it would not solve my original issue. Experimentation proved that to be true.

Thank you for your help.

Now I'm totally stumped...

[Baktash A. Shamshirsaz]

I think this will solve your problem:
http://3rdprogress.no-ip.com/archive/num7.html

In a nutshell, you need to take into account ALL contact points at the same time.

Cheers,
Baktash.

[SigKILL]

In a nutshell, you need to take into account ALL contact points at the same time.


Yup, thats correct. There is an explicit solution to this, which results in solving a LCP (Linear complementarity problem). The LCP is unfortunally a quadratic optimization problem, but there are specific algortihms (the fastest is probably an algorithm by dantzig, which is what ODE uses).

However, the most common method (IIRC Novodex does this) is probably something close to http://graphics.stanford.edu/papers/rigid_bodies-sig03/ , it is easier on the math (does not involve an LCP) and is fairly easy to implement. It is a nice read too, something that is kind of seldom in the field of computer science IMO.

-si

[Baktash A. Shamshirsaz]

Well, the tutorial that I wrote basically touched on the same deal. You COULD set constraints and use Linear Complementarity Problem but in most cases you can set up an augmentation matrix and solve it using Gauss/Jordan method. I've tried it in a lot of simple cases. It actually works this way as well :).

[SigKILL]

Well, the tutorial that I wrote basically touched on the same deal. You COULD set constraints and use Linear Complementarity Problem but in most cases you can set up an augmentation matrix and solve it using Gauss/Jordan method. I've tried it in a lot of simple cases. It actually works this way as well :).

hmm, your link doesn't seem to work here, so I'm not sure exactly what you're talking about. Are you basicly finding the amount of force to apply to avoid penetration ,using a linear system? I should really play with thison paper before I ask these questions ;). Or are you simply inverting the interia matrix (I usally store the inverse of the interia matrix with the body).
I was somehow mistaken this discussion of being about contacts and not collisions (collisions is actually quite simple when you know what to do ;) I think you'll find the solution in any dynamics book ). Do your method extend to include friction? (wich is a must IMO).

-si

[Baktash A. Shamshirsaz]

Ok, you mentioned a handful there :D.
First off, sorry for the link being down. The server was down for a few hours last night thanks to my no good ISP :S. You can check it out now.
Now let me clear up a few things about the analytical aproach and Impulses and Contact Forces.
----------The Analytical Method----------
For solving for contact forces, basically people like Baraff chose the analytical method (the same one I'm talking about) to calculate contact forces that exist at 'RESTING' objects. So first they go through the trouble of finding out whether they're resting or not and then for probably 'educational' purposes they go the extra mile to calculate the normal forces :). But the method that they use for the objects that are NOT resting and ARE colliding is prone to errors, because out of all collision or should I say 'contact' nodes, they just chose one or one at a time and calculate impulse(s). That is physically WRONG and does not match real life responses. As mentioned in several other papers (and even Baraff's) all simultaneous impulses/forces (or let's say responses) affect eachother at the same time, so they must be all solved together. So if a table hits the ground with all it's four legs all four impulses affect each other's magnitudes (as if there were 5 legs the magnitudes would be slightly less.) That is why instead of worrying about normal forces applied over long periods of time, I (and several other researchers like Brian Mirtich) take everything as impulses, even resting as multiple simultaneous micro collisions. Well, with this approach we can use the analytical method for solving for those impulses all at the same time and therefore create a result that is noticably close to real life situations (the smaller time steps, the smaller distance -- epsilon -- in between the rigid objects, the more realistic response.) Now a little bit on why I'm using Linear Systems to solve this. In most papers that I've read they first construct an equation for every node with a number of unknowns being equal to the number of collision nodes and then try to satisfy the equations using LCP in an iterative manner. The reason why is because they consider the final relative velocity with respect to contact normal to be equal or greater than coefficient of elasticity times the incoming relative velocity w.r.t. contact normal (because of a POSSIBLE simultaneous third collision behind the 2nd object.) For a situation like this your equations turn into inequalities (equal or greater than) and LCP would be a wise choice! Well instead of using this approach I just use only equations and settle for a maximum of 1 frame error possibility :) (in the next frame the problem will be automatically solved through the collision solver loop) and solve for the magnitudes right away! Although I do some pruning on the nodes first (I make sure I have the right normals and no collision node is repeated twice.) I must also stress on the fact that there are some cases where Gauss/Jordan will not be able to find a unique answer for the system and you HAVE to rely on LCP, but if you know your cases are simple enough you don't have to go for it. On all cases that I've tried this (mostly simple), I have gotten correct answers. Also I believe it is possible to incorporate friction here as well. I believe you can just think of it as an impulse applied on contact nodes in the opposite direction of the projection of their total velocities (linear + tangential) on the contact plane. You might ask 'what about it's magnitude? :)'. For the magnitude, we can use culomb's friction law! Fk = Mu * Fn! Here simply replace the forces by forces applied over collision time, or in other words impulses :)! Therefore we can have, Jk = Mu * Jn! Or normal impulse applied by contact surface times coefficient of friction. And I believe you can figure out how to deal with static/kinetic friction from here on! :)

Man that was a handful :D.

[Baktash A. Shamshirsaz]

Oh and here's one trick that I read in a paper by Mirtich that was published in 1998. It is called the SVD approach. You can't believe how easy it is and how funny it is that we have been overlooking it for so long.
In general, Af + b = 0 is the form of our equation for colliding contacts for two objects (if we ignore the greater than zero possibility) 'f' is the matrix or should I say the array of contact forces/impulses that we want to calculate. Here's a nice way:
1. Af = -b
2. f = A_inverse * (-b)
Therefore, f = A_inverse * (-b), and you're done! Make sure you have the correct nodes and more importantly the correct normals.
Researchers are right now working on finding out what types of contact nodes create an 'A' matrix that is not invertible so they can make this method bullet proof. Until then you're advised to use LCP.

[SigKILL]

I have a few 'issues' with you article, however, you really deal with the "easy" case (i.e. when contact is resolved instantly). The problem starts when you want to deal with what I (possibly wrongfully) call 'contacts' i.e. contacts that can't be resolved at the current timestep (like objects sliding or resting on eachother (i.e. stacking) and possibly other objects colliding with these). Well, two thing can happend, either there is enough force/impulse so the objects start moving away from each other (i.e. contact resolved), or there is just enough force/impulse so the objects don't penetrate. If you find a description of the LCP it will say what I said in the last sentence in a mathematical way.

Mirtich is talking about using singular value decomposition (SVD) to find
what he calls an "inverse" (not the real inverse). I'm not really sure what he means by this, but he is possibly reffering to an approximate inverse. The nice thing with the SVD is that you get two orthogonal matrices (i.e. non-singular) and one possibly singular diagonal matrix (which 'measures' the singularity of the matrix). My guess would be that he multiplies the right side with the orthogonal matrices, and solves with respect to the diagonal matrix (putting possibly zero where the value of the matrix is close to zero). I may be wrong, I just skimmed the article just now.

Anyways, if you really want to avoid solving an LCP, you should really read the paper I linked above. It provides a very nice solution, and simple to implement compared to alot of other techniques. From section 5 an onwards they really described very detailed how to create a robust rigid body simulator (for a realtime implementation you probably want to cut down on all the iterative stuff). They do not treat all collision points at once as you do though.

-si

[Baktash A. Shamshirsaz]

If you think about it, Baraff is doing the whole thing over the long run. He is applying 'forces' to keep them from interpenetrating. Here we're giving the object little velocity and by the very little impulse applied we're keeping it from moving into the other object. If you zoom in let's say in ODE or even in Half-Life 2 (which uses Havok) you will see that objects go below/over epsilon distance very fast but for very small amounts of time (usually you can find a jittering pixel for some time.) Here (I believe) it is because of the micro collisions and little impulses applied to keep objects from interpenetrating. Even if they're not using the impulse-based approach (it doesn't really matter to me which approach one choses), if they did it would produce the jittering pixel as I've seen it in my own tests :).
Mirtich is talking about using singular value decomposition (SVD) to find
what he calls an "inverse" (not the real inverse). I'm not really sure what he means by this, but he is possibly reffering to an approximate inverse. The nice thing with the SVD is that you get two orthogonal matrices (i.e. non-singular) and one possibly singular diagonal matrix (which 'measures' the singularity of the matrix). My guess would be that he multiplies the right side with the orthogonal matrices, and solves with respect to the diagonal matrix (putting possibly zero where the value of the matrix is close to zero). I may be wrong, I just skimmed the article just now.
I'm afraid you did not get the main idea. Just for the note, the 'b' vector in the impulse-based approach is (1+e)incoming_relative_velocity_wrt_normal and the 'A' matrix is the matrix of relative changes of total velocities wrt normals and 'f' is the array of impulses (I've clearly explained this in my tutorial.) f = A_inv * (-b) would clearly give you the impulses. Just as if you were to make an augmentation matrix of Af=-b and solve it using Gauss/Jordan.Whether the objects just 'stick' or move away depends on how much the coefficient of elasticity is (or in other words how much moment is transffered through a collision.) I actually worked a bit on a special Gauss/Jordan eliminator that calculates correct impulses for all contact points and even does a little hack around non-invertible 'A' matricies to avoid leaving us with no answer :). It has passed a few tests for now. But I'm still putting more complex tests on it. I'm sure it will work for those as well.

For objects sliding on each other, again, as I said, solving friction using culomb's friction law is a breeze and the 'resting' part is solved using micro collision as I also mentioned before. You should try reading some of Mirtich's papers on Impulse-based method for dynamic simulation or even his Ph.D thesis which is on the subject. You will see that contact 'forces' and LCP are not the only approaches that we can take. Aside from all this there are also energy based approaches that I unfortunately never got to look through. Anyway, thinking of EVERY thing as impulses works as well, if we do it over small time steps.

They do not treat all collision points at once as you do though.

If you read the 1st post in this thread, you will realize how WRONG NOT doing that is with more than 1 simultaneous collision point. That is mostly what my tutorial is trying to teach the reader. You can refer to the 1st reference mentioned in my tutorial for more information. Although you can solve impulses separately inbetween different 'pairs' of objects and apply impulses.
-------------------
Also from what I've observed, I know where our fundamental differences of thoughts arise from :). You're pretty much following Baraff's contact force calculation using LCP or basically QP where results (contact 'forces') are applied as a 'forces' directly in the force resolution loop. I'm generally taking a different twist on the subject by considering everything as impulses (even resting as millions of simultaneous impulses) and applying results as impulses on the object (changes of momentum no matter how large or small) but again we're solving the problem through an almost identical mathematical way (solving contact node equations.) Although the way we're solving this is quite different. You're taking a non-linear approach and I'm taking a linear approach. Regardsless of how we're doing this, we're both keeping multiple objects from interpenetrating in a physically realistic fashion.

[SigKILL]

About a year ago a math professor at my uni told me about an annoying kid from high school trying to teach him how to differentiate (the kid didn't really understand the subject properly). Now I think I know how the professor felt...

-si

[Baktash A. Shamshirsaz]

I'm really not going to argue with you over this as I've got more important things to attend to. Unfortunately after times and times of trying to make you understand a new approach you're still trying to push all the tired details in your paper into my mind. Instead of doing this go read the papers that I mentioned in my tutorial's reference and Mirtich's introduction to the impulse based approach. And just for the note, my simulation results are matching others and the details that I've given in my posts are all physically and logically correct (go ask anyone who's done this before.) If you have a problem understanding them, that's not my problem, that's yours. This will be my last post. If you get it, good for you! If you don't, well, I couldn't care less.