Jump to content

  • Log In with Google      Sign In   
  • Create Account

Dirk Gregorius

Member Since 24 Apr 2002
Offline Last Active Today, 12:08 PM

#5252537 Realtime collision with dynamic concave meshes

Posted by Dirk Gregorius on 16 September 2015 - 10:57 AM

This is a very tough problem to solve and I am not aware of any solution for this in realtime. If I needed to solve this I would start from this paper using SDF:



The SDF are very intensive in memory (e.g. PhysX had an implementation of this called PMeshes for a while, but dropped it again because of the high memory usage) so you would need to come up with some kind of compression. I heard there has been some good experience with low band SDF in the FX field recently, but I have no practical experience. In games you usually take much simpler approaches. Dennis gave a presentation on this at GDC here:



There are some good SDF libraries out there which might be a good start. The authors of the libraries might be good initial contacts as well:



HTH and good luck!


#5251957 Problem implementing sequential impulse mass over 100Kg (solved, it was a mis...

Posted by Dirk Gregorius on 12 September 2015 - 06:55 PM

Did you disable position correction? I can be everything from wrong contact points, to wrong Jacobians, to wrong inertia. Whatever this is it is unrelated to the mass!


Is z your up axis?

And use temporaries! How do you expect that anyone can read this: vr.set(vr2.set(b2.rotation).cross(r2).add(b2.velocity)).sub(vr1.set(b1.rotation).cross(r1).add(b1.velocity));

#5251588 How to determine closest point on plane vs aabb

Posted by Dirk Gregorius on 10 September 2015 - 01:08 PM

I don't understand why you want to return the projected center in case the plane is axis aligned. The vertices should be equally good. Here is how I would break this up. Maybe it helps.

Vector3 ClosestPointOnPlane( Bounds3 aabb, Plane plane )
    Vector3 support = aabb->Support( -plane.m_Normal );
    return plane.Project( support );

Vector3 Bounds3::Support( Vector3 v )
    Vector support;
    support.x = v.x < 0 ? m_Min.x : m_Max.x;
    support.x = v.y < 0 ? m_Min.y : m_Max.y;
    support.x = v.z < 0 ? m_Min.z : m_Max.z;

    return support;

Vector Plane::Project( Vector3 p )
    float distance = dot( m_Normal, p ) - m_Offset;
    return p - distance * m_Normal;

#5251425 How to determine closest point on plane vs aabb

Posted by Dirk Gregorius on 09 September 2015 - 04:47 PM

Just get the support point on the box in the negative direction of the plane normal. This will give you both the closest point and penetration

#5250656 How can I implement warmstart?

Posted by Dirk Gregorius on 04 September 2015 - 11:07 PM

Base your engine on Box2d Lite and study all of the GDC presentations and you will be fine! :)

#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