This post is part of my Game Physics Series.

Contact constraints are amongst the most important constraints for a physics engine, because you would want to resolve collision in most scenarios. And if not done properly, your rigid bodies can behave in a very visually jarring way to the player.

### Derivation

Remember that Sequential Impulse models constraints in the form of , where is the __Jacobian matrix__, is the __velocity vector__, and is the __bias term__. Erin Catto (author of Box2D) proposed a systematic way to derive velocity constraints: you first find the position constraint , and then differentiate it with respect to time to obtain the velocity constraint . Sometimes, it is more intuitive to directly derive the velocity constraint itself, and I think it is the case for contact constraints. However, for the purpose of demonstration, I will first find the equation for the position constraint, and then derive the velocity constraint.

Based on the contact format previously described, a possible interface for contact data may look as follows.

struct Contact { // contact point data Vec3 globalPositionA; Vec3 globalPositionB; Vec3 localPositionA; Vec3 localPositionB; // these 3 vectors form an orthonormal basis Vec3 normal; // points from colliderA to colliderB Vec3 tangent1, tangent2; // penetration depth float depth; // for clamping (more on this later) float normalImpulseSum; float tangentImpulseSum1; float tangentImpulseSUm2; Contact(void) : normalImpulseSum(0.0f) , tangentImpulseSum1(0.0f) , tangentImpulseSum2(0.0f) { } // there's typically more stuff // but omitted for simplicity's sake };

Let us look at the figure below before moving onto finding the position constraint.

and are centers of mass of collider A and collider B, respectively. and are two contact points (deepest penetrating points) on two colliders. and are vectors from centers of mass to contact points. is the contact normal (always pointing from collider A to collider B for consistency). and are contact tangents, and is pointing “into” the screen.

We basically want our penetration depth to be zero, so the position constraint states that: the vector from to must be in the same direction as (the dot products of both vectors, i.e. the penetration depth, is ).

We can see that:

So, the position constraint is:

By differentiating with respect to time, we get the velocity constraint :

By some algebraic manipulation and making use of the triple product identity , we get:

where:

Note that the form of our velocity constraint is instead of . However, don’t worry. This simply means that we impose the constraint when the two colliders are penetrating.

The geometric interpretation of the velocity constraint we just derived is as follows: __the projection of the relative velocity of the two contact points onto the contact normal is zero__. This means if we solve this constraint, the two contact points would not penetrate further.

Recall that the correction to the velocity vector is , where is the Lagrangian Multiplier:

and is the mass matrix:

### Clamping The Lagrangian Multiplier

We have found out how to solve the contact constraint. However, this is not the entire picture. Remember there are usually more than one contact constraints in the entire physics system. The Sequential Impulse method iterates through all contacts and solve each contact constraints one-by-one:

for (Manifold &manifold: manifolds) for (Contact &contact : manifold.contacts) SolveContactConstraints(contact);

The change to a collider’s velocity will affect the result of later resolution of the same collider. At one point, you may end up with a velocity change that actually pulls a collider pairs closer instead of pushing them apart. For this reason, we need to perform some clamping when applying velocity change using . This is the reason why in the interface for contact shown earlier contains the data members `normalImpulseSum`

.

We want to make sure that throughout all iterations within one time step, the following inequality is satisfied:

where is the Lagrangian Multiplier computed during the i-th iteration.

Each time we compute for a contact, we make a copy of `normalImpulseSum`

, add to `normalImpulseSum`

, clamp `normalImpulseSum`

between 0 and positive infinity, and then calculate the difference of the clamped `normalImpulseSum`

from the previously copied `normalImpulseSum`

. This difference is the actual Lagrangian Multiplier we want to use to solve the contact constraint.

### Baumgarte Stabilization

Solving the velocity constraint is not enough. Remember that by solving the velocity constraint, we only prevent the two colliders from further penetrating. It does nothing to actually push them apart. The accumulated error results in “positional drift”.

This is when Baumgarte Stabilization comes into play. Some call this method a hack, some view it as an elegant way to fix positional drifts. Either way, it works.

The basic idea is to feed the penetration depth (position error) back into the velocity constraint as a bias, so just a little bit of energy is introduced into the system to push the two colliders apart.

As previously shown, the bias term is zero. But it is not zero any more when Baumgarte Stabilization comes into the equation.

To compute the penetration depth , we use what we already got from the position constraint equation:

We then feed this back into the constraint equation after multiplying with a __Baumgarte term__:

where is typically a value between 0 and 1 (usually close to 0). Unfortunately, there is no one “correct” value for ; you will have to tweak this value based on the specific scenarios in your game. This is one of the minor quirks about constraint-based physics engines, lots of parameter tweaking. Erin Catto proposed that you start increasing the value from 0 until your physics system starts to become unstable, and then use half that value.

### Restitution

Wait, what about having the colliders “bouncing” apart? So far by solving the velocity contact constraints only make sure that the two colliders do not penetrate any further. There is another quantity we are going to add to the bias term in the constraint equation in addition to the Baumgarte term: the __restitution term__.

The __coefficient of restitution__ between two objects, denoted , is defined as __the ratio between the parting speed of two colliders after collision and the closing speed before collision__; for realistic result, you would want to set this value between zero and one.

Therefore, the contribution of the restitution term to the bias term is simply the coefficient of restitution times the projection of the closing velocity of two contact points onto the normal vector:

So our final bias term with contribution from both the Baumgarte term and restitution term is:

### Frictions

So far we have only covered the normal part of contact constraints. Now I’m going to talk about the tangential part of contact constraints, a.k.a. frictions.

Frictional forces are proportional to normal forces, so they are dependent on the result of contact resolution in the normal direction. Thus, we typically solve the tangential part of contact constraints after the normal part.

The normal part of contact constraints attempts to zero out the projection of the relative velocity of two contact points onto the contact normal (without considering the bias term , that is). The tangential part of contact constraints does something very similar: it attempts to zero out the projection of the relative velocity of two contact points onto the contact tangent plane (formed by the two contact tangents).

The tangential part and normal part of contact constraints are basically doing the same thing, except that the directions onto which the relative velocity of the two contact points are projected are different. So there is no surprise to see that the Jacobians for the two contact tangents look very similar to the Jacobian for the contact normal. Let’s relabel the Jacobian for the contact normal , and put it side-by-side with the Jacobians for the two contact tangents, denoted and .

Pretty similar, right?

What is also similar is that the Lagrangian Multipliers computed for the contact tangents also need clamping. However, instead of clamping between zero and positive infinity, the boundaries the Lagrangian Multipliers are clamped depend on the normal impulse sum, which is stored in the `normalImpulseSum`

data member. The clamping boundaries are negative and positive `normalImpulseSum`

multiplied by the __coefficient of friction__, denoted of the two colliders, which is typically between zero and one. Let us also denote `normalImpulseSum`

as .

The clamping process is extremely similar to the normal part: we make copies of `tagentImpulseSum1`

and `tangentImpulseSum2`

, add the Lagrangian Multipliers computed from both tangential Jacobians to `tangentImpulseSum1`

and `tangentImpulseSum2`

, clamp `tangentImpulseSum1`

and `tangentImpulseSum2`

between and , and then respectively calculate the difference of the clamped `tangentImpulseSum1`

from the previously copied `tangentImpulseSum1`

and the difference of the clamped `tangentImpulseSum2`

from the previously copied `tangentImpulseSum2`

. These differences are the actual Lagrangian Multipliers we want to use to solve the tangential part of contact constraints.

You might notice that this clamping method is not perfect. We can end up with clamping the magnitude of total tangential impulse to . However, this error is usually not noticeable, so I usually just let it be (Joshua Davis, a member of DigiPen Institute of Technology, informed me of this idea).

### End of Contact Constraints

This is the end of the derivation of contact constraints and demonstration on how to solve them. I will cover stability issues and how to improve stability using a technique called “warm starting” in later posts.

Helped me a lot to understand it! In 2D, what will the mass matrix be?

You just take out the rows and columns related to the Z dimension.

Hi!

How do you derive that differentiating rB is wB x rB?

Thx

The derivative of rB with respect to time is (vB + wB x rB), where vB is the linear part and (wB x rB) is the angular part of the velocity at rB.

Thanks, I think I get it now.

Hi, I’m working on a physics engine based on constraints. At the point “Clamping the Langrangian Multiplier” I dont understand what you are doing. What is the variable

`normalImpulseSum`

for? Is it just the sum of the impulses of the two bodies? And when I read in Erin Catto’s papers he made a large matrix containing all contraints at the same time. Is this a similiar method to the one you’re describing? By the way great website, took me a great step forward:)1. Throughout the iterations within one time step, sometimes you will get a negative normal impulse for one contact as a result of resolving other contacts., which will pull penetrating objects even closer; this is not what we want. The variable

`normalImpulseSum`

is for keeping track of the total normal impulse for a contact throughout the iterations to ensure that it never becomes negative.2. Theoretically, we will get the optimal results by solving all constraints at once in a big matrix. But this is not possible due to limits on computational power. So we use the Sequential Impulse approach, where we solve one to a few constraints at a time, and the we iterate through all constraints multiple times in a single time step.

Hi Allen, i have tried to implement it, but does not work, can i see your implementation please?

The only implementation I have is the one I did for my school project at DigiPen, which I am not allowed to disclose because it belongs to DigiPen. However, you can check out Erin Catto’s implementation in Box2D. I used it as a reference.

Or, you can show your implementation using online paste tools like Pastebin so I can take a look at what’s wrong.

Hi, according to this method, every thing will be solved by velocity changes, but i’m thinking of where is the position constraint solver?

don’t we need to apply Baumgarte directly on our bodies position?

By applying Baumgarte stabilization to solve velocity constraints, we are feeding back positional errors into velocity changes. This has the effect of correcting the positional errors over a few time steps, instead of instantaneously. To apply instant positional correction, you basically integrate the position using the velocity change computed from Baumgarte stabilization and the delta time. This approach can potentially give you faster positional correction, but might introduce more jitter if you have too many conflicting positional constraints. I personally prefer just using Baumgarte stabilization to solve for velocity constraints. It takes multiple frames to correct positional errors, but it’s visually smoother and more stable.

Hi there Erin, I imagine this is because the constraint space for contact constraints is relative to the collision normal itself. Would this be why the orientation of the collision normal in world space is irrelevant?

Hi Allen. I like your site and you have a nice blog! Also congratulations on your position at Naughty Dog!

When you differentiate C, did you consider that the normal vector is not constant? It can rotate with time. But somehow it doesn’t matter that you skipped it. Why?

Wow, it’s Erin Catto himself! I am a big fan of your work, and I love your presentations at GDC and DigiPen.

Actually, I didn’t consider the fact that the normal vector is not constant when I derived the constraint equation, and everything works quite fine for me. Now that you mentioned it, I should probably go back and check my derivation…

[Edit]

On a second thought, a new normal is computed upon every new contact data generation. Normals are treated as constant only when caching (warm starting) is at work, which happens due to frame coherency, so the normals shouldn’t change too much. Do you think it still matters?

Hi Allen, thanks a lot for writing this blog. It’s awesome and I’ve learned a lot. I was wondering about this issue too… The neglected term is basically the penetration depth vector dotted with the derivative of the normal vector. If you look at Erin Catto’s paper “iterative dynamics with temporal coherence”, he mentions that this term is usually neglected as an approximation because the penetration depth is usually small. I have two questions about that… firstly, since we usually use discrete collision detection, why is this assumption valid? Can’t the penetration depth be arbitrary? And also, since the chain rule guarantees C dot=Jv, doesn’t that mean we should still be able to simplify the entire expression into this form, including the neglected term?

Best articles on the subject I have come across! I will need to implement constraints in a project I’m working on and thanks to these articles I’m starting to understand how it all works. Thank you very much. Keep writing them!

You are very welcome. I’m glad this article helped.