Jump to content

  • Log In with Google      Sign In   
  • Create Account

We need your help!

We need 1 more developer from Canada and 12 more from Australia to help us complete a research survey.

Support our site by taking a quick sponsored survey and win a chance at a $50 Amazon gift card. Click here to get started!

Dirk Gregorius

Member Since 24 Apr 2002
Offline Last Active Yesterday, 11:15 PM

#5250516 How can I implement warmstart?

Posted by Dirk Gregorius on 03 September 2015 - 08:27 PM

This has been exhaustively discussed in the Game Physics Tutorial at the GDC over the years. You can find all presentations here:



In particular "Fast and Simple Physics using Sequential Impulses" from 2006. It comes with a stripped down 2D implementation of Box2D (called Box2D Lite). I would starting studying this.  

#5250325 Odd issue with Separating Axis Theorem

Posted by Dirk Gregorius on 02 September 2015 - 03:36 PM

I gave two exhaustive talks about the SAT and contact creation. Maybe they are helpful:


Robust Contact Creation for Physics Simulations



The Separating Axis Test Between Convex Polyhedra


#5248635 Contact creation

Posted by Dirk Gregorius on 24 August 2015 - 04:36 PM

Oh, I might have misunderstood your question. In the case of the velocity iteration we still have the following:


dC/dt = ( v2 + w2 x r - v1 - w1 x r1 ) * n + ( x2 + r2 - x1 - r1 ) * dn/dt


But now either x2 + r2 - x1 - r1 = 0 if you setup the local contact points from global ones. Or in the case where the local points are disjoint (I think you asked about this earlier) you compute the normal from the difference vector anyway in which case the cross product  ( x2 + r2 - x1 - r1 ) x n would vanish since the vectors are parallel.


This might not matter in practice at all. You have a similar problem with prismatic joints. If you ignore the right term you get problems if the bodies are sliding apart quite a bit. Penetrations are usually small so it might not matter. Formally this would be the correct way though. 

#5248623 Contact creation

Posted by Dirk Gregorius on 24 August 2015 - 03:00 PM

No, velocity solving is a linear problem and during iterations you are projecting velocities onto the 'fixed' constraint manifold. Position projection is a non-linear problem (hence the N for Newton in NGS). The Jacobian is a function of the position. I guess you could try to anticipate the change in the position due to the velocity change, but I am not sure if this would makes sense.

#5247742 Cone/angle joints in impulse method

Posted by Dirk Gregorius on 19 August 2015 - 02:57 PM

You should be able to model a small ragdoll this way. Ragdolls are usually spherical joints with cone limits and revolute joints. Did you look at the angle limits? It should be fairly easy to convert this data. I have done this many times. 

#5247411 Cone/angle joints in impulse method

Posted by Dirk Gregorius on 18 August 2015 - 10:31 AM

This yields a 4x4 system when solved together with the point to point constraint. You have a couple of options now:

1) Solve each scalar constraint individually (this is essentially what the ODE dQuicksolve does)
2) You solve the point to point constraint as 3x3 block constraint and the limit as a scalar constraint. This is the common way these days from my experience
3) You solve the point to point constraint together with conical scalar constraint. This yields a 4x4 MLCP which can be solved with direct enumeration

Erin does the 3) in Box2D. This yields the best results as it boost converges since the point to point constraint and angular limit don't fight each other. In 3D I found this very expensive so I use 2) personally in Rubikon.

Since in an iterative solver the order in which we solve constraints matters the remaining question is if I solve the point to point constraint before the limit or vice versa. I solve the point to point constraint always last to emphasize it. I found that stretch at anchor points is more visually noticeable than a violated constraint.

Here is an example for Direct Enumaration from Murphy's webbook:

#5247409 Cone/angle joints in impulse method

Posted by Dirk Gregorius on 18 August 2015 - 10:24 AM

For each joint we need to define a body-local joint frame. Here is a sketch that hopefully illustrates the idea:


Of course you want the two local anchor point to be coincident using the point-to-point constraint. From what I understand you already figured that out, so I skip this for now, but we can talk about later if you like or need another example of what I will show now below.


Now we need to model the conical limit. A very simple way to do this would be using the following constraint:


C = u2 * u1  > cos( theta / 2 )


Here theta is the (full) cone angle and * denotes the dot product. From there we follow the basic procedure and build the velocity constraint:


dC/dt = du2/dt * u1 + u2 * du1/dt


The change of a fixed axis in a moving frame is simply omega x axis. Hence


dC/dt = ( omega2 x u2 ) * u1 + u2 * ( omega1 x u1 )


Using the scalar triple product identity and arranging terms should then yield:


dC/dt = ( u2 x u1 ) * ( omega2 - omega1 )


From the velocity constraint we can now identify the Jacobian by inspection:


J = ( 0^t -(u2 x u1)^T 0^T (u2 x u1)^T )


I used the ^T here to make clear that this is a 1 x 12 matrix. 


So in the solver we just check the position constraint. If it is satisfied we can ignore this constraint in this frame. Otherwise we need to solve:


J * W * JT * lambda = -J * v - beta * C / dt


This approach has a couple of advantages. Using the time derivative and utilizing that dC/dt = del_C/del_x * dx/dt + del_C/del_t is much easier than trying to build the gradients directly. Though a good geometrical understanding helps a lot to simplify constraints and Jacobians. You also have the stabilization term naturally handy and don't need to guess it like when you would start from the velocity constraint. BTW, the fancy formula is called the total derivative and a summary can be found here: https://en.wikipedia.org/wiki/Total_derivative


A little implementation gotcha. In your implementation you would store the constraint axis local in each body frame and then u1 = R1 * u1_local and u2 = R2 * u2_local to identify the current joint frames. Then you basically simply inspect the relative orientation.

#5247277 Contact creation

Posted by Dirk Gregorius on 17 August 2015 - 05:10 PM

Erin did a good summary on position projection here:



The Jacobians are a function of the position. What makes position correction expensive is that you need to rebuild the Jacobians to do proper position correction. You can keep them constant. This is referred to as split-impulse in Bullet. This works kind of ok with contacts, but not so much with joints from my experience.


I recommend looking at Box2D instead. It has what Erin calls a full NGS implementation.

#5247003 Contact creation

Posted by Dirk Gregorius on 16 August 2015 - 07:17 PM

Interesting! I would need to check how my code handles this case. 


I didn't add restitution to the bias (actually I don't use a bias a at all these days since I use a separate penetration recovery sweep - sometimes called position projection or NGS). I think I took the max of the two. You also don't apply restitution if the relative velocity falls below a threshold. I made good experience with 1 m/s as upper bound, iirc the default in PhysX is 2 m/s (YMMV). And another thing two be careful about is to compute the relative from the body velocities before you applied the warm-starting impulse.


Here is how this looks in my code:

if ( Restitution > 0.0f )
	// IMPORTANT: Compute the relative velocity at the contact point from the bodies and
	// not from the velocity buffer (which can already contain the warm starting impulses). 
	// Ideally the velocity should not contain integrated gravity as well! 
	rnVector3 V1 = Body1->GetLinearVelocity();
	rnVector3 W1 = Body1->GetAngularVelocity();
	rnVector3 V2 = Body2->GetLinearVelocity();
	rnVector3 W2 = Body2->GetAngularVelocity();

	float RelV = rnDot( V2 + rnCross( W2, Offset2 ) - V1 - rnCross( W1, Offset1 ), Normal );
		Out.Bias = Restitution * RelV;

#5246806 Contact creation

Posted by Dirk Gregorius on 15 August 2015 - 07:29 PM

Oh, also make sure to cache a separating axis. Basically in the case penetration you re-clip or rebuild closest points. If you found a separating axis test this first the next frame before you run the SAT as well.


You should see the idea here. We only run SAT if we are forced with guns and knives to do so smile.png

#5246804 Contact creation

Posted by Dirk Gregorius on 15 August 2015 - 07:14 PM

In the edge/edge case I simply recompute the closest points. You need to detect if the edges slide off each other. Luckily this is simple. I just compute the closest points between the lines through both edge vertices just check that they are actually on the edge. If t < 0 or t > 1 re-run the SAT as well.


Yes, you can use the smart SAT as well. The test is even easier. Randy wrote a summary here: 



I didn't check if it is correct though, but it looks good. Maybe a bit more complicated than it needs to be.

#5246489 Rigid Body Physics - Inertia Tensor Space

Posted by Dirk Gregorius on 14 August 2015 - 10:06 AM

If you multiply omega by dt the length of the resulting vector is Phi. It is kind of a simplification since you assume a constant direction of omega over the integration interval. 


f(x) = sin(x) / x in the limit x -> 0 is of the form 0/0. Interestingly the actual limit is f( x -> 0 ) = 1 which can be shown with the l'hopital rule.



I don't think that overloading operator+ is a good idea because if someone else reads your code it might not be immediately apparent to him what your are trying to do. Writing a simple function IntegrateOrientation() would be much more meaningful in this case (e.g. look for self-documenting code). I also like the term of 'navigation depth' in that context. It means if I read a function in some code how deep do I need to navigate into sub functions to understand the code. Ideally I can understand the overall function by just looking at the top level code. It makes you a lot of friends writing readable code.


Regarding the original problem I think we have been through all details and you should be able to figure out the problem yourself. Just take your Unity example and then step through both programs side by side until you find the problem. This shouldn't really be too hard! 

#5246327 Rigid Body Physics - Inertia Tensor Space

Posted by Dirk Gregorius on 13 August 2015 - 05:21 PM

Yes, I found this better for "small" angular velocities. What you describe is sometimes referred to as "exponential map" or "Grassia's" method. You need to take some extra care if Theta -> 0 if I remember correctly.  


Here is the reference if you are interested:



Section 3.0 is the part of practical relevance. A good read in general though. 

#5246315 Rigid Body Physics - Inertia Tensor Space

Posted by Dirk Gregorius on 13 August 2015 - 04:18 PM

The quaternion multiplication and using body and not the world space angular velocity in the quaternion integrator are two things that come to mind here which might be causing your problem.  The first is difficult to tell without seeing your quaternion multiplication. 


The quaternion derivative is defined as q' = 1/2 * w * q. Here I am assuming that we use a quaternion multiplication such that q = q2 * q1 (e,g. Shoemake) and w is the pure quaternion holding the world space angular velocity. If you want to use the body space angular velocity the formula conveniently changes to q' = 1/2 * q * w. Here w now holds the body space angular velocity. This is very easy to prove. Just plug w_global = q * w_local * conjugate( q ) into the formula. 


Also note that you need to normalize the quaternion after integration. And overwriting operator+ for the integration is a poor choice and totally misleading in my opinion.

#5246062 Contact creation

Posted by Dirk Gregorius on 12 August 2015 - 02:47 PM

Thinking about it, putting the contact points always onto the surface of one shape might be a good choice for games as well. E.g. the normal points from A to B and the contact points are always on the surface of A is consistent and should help gameplay programmers.


The contact stuff has been solved long time ago and I didn't really think about it recently. I agree a more consistent choice might be better!