# EDIT:

This series ended up on Tuts+, and these are sort of deprecated. Please visit Tuts+ to see the finalized articles.

# Personal Update

Beyond C++ reflection I’ve taken a huge liking to physics programming for games. Reflection in C++ is a largely “solved” area of study. Though resources on the internet may be few and far between for creating custom introspection, people that know how to create such systems have explored much of what is to be accomplished.

However physics is another story entirely. Up until just a few years ago physics for games was in a highly primitive state; techniques just hadn’t been created due to the inability of hardware to perform necessary computation. There are a lot of new and exciting problems to solve.

So this next year while attending DigiPen I’ll be focusing my studies around: engine architecture, data oriented design and optimization, C++ introspection and game physics. This last topic being the topic I’m going to start an article series on.

# Game Physics High Level Detail

A game physics engine performs some functionality for ensuring that shapes behave in a particular way. To state this one could say: a physics engine removes degrees of freedom from moving bodies, and enforces these imposed rules.

From here on out, until specified otherwise I’ll be talking specifically about convex rigid bodies within a physics engine. A rigid body is just a shape that is implicitly rigid; the engine does not naturally support deformation of a set of points that make up a shape. Using this implicit definition the mathematics behind manipulating such shapes can be known intuitively and implemented efficiently.

Here’s a short feature list of some features a custom physics system may provide, assuming the engine uses rigid bodies:

- Collision detection
- Collision resolution
- Linear movement, angular rotation and integration
- Raycasting
- Islanding/sleeping/portals
- Friction simulation
- Spatial partitioning/bounding volumes
- Fracturing/splitting
- Multi-part bodies
- Various types of shapes
- Advanced mathematical constraints (joints, motors, contact constraints, etc.)

In this article series I will attempt to talk about all of the above topics in my own time. I’ve learned a lot about physics and physics engines in a short amount of time and by documenting what I have learned I hope to solidify my understanding, as well as help out others to achieve what I have. I would not have learned nearly as much as I currently know without help from others so it is only natural to want to do the same.

The best example of an open source physics engine that employs much of the feature list shown above would be Box2D by Erin Catto. The rest of this article series will be detailing the physics engine that I myself have written. There are of course other methods I choose not to talk about, or just don’t know about.

# Architecture

There are two main objects that make up a physics engine: shapes and bodies. A body is a small package of data that defines intrinsic properties of a physics object. Properties such as density, restitution (elasticity), friction coefficients, inertia tensors, along with any other flags or settings should be stored within the body. These bits of data are properties that can be isolated away from the shape definition. The shape itself is contained with the body through a pointer.

The shape definition defines what sort of shape a physics object represents. Some physics engines can attach multiple shapes to a body, which are referred to as fixtures. A shape stores any data required to define the shape itself, and provides various functions to operate on the shape, such as testing for line or point intersection, or generating a bounding volume.

Together the body and shape represent a physics object, which by the end of this article series will be able to perform a lot of interesting interactions with other physics objects.

All bodies should be contained within what is known as a scene or world. I refer to this object as a scene. The scene should contain a list of all live objects, as well as functionality inserting or removing bodies from the scene. The scene should also have callbacks for processing shape or ray queries. A query just checks to see if any bodies intersect with something like a point or ray.

The scene has one particular function called step, which steps the scene forward in time by a single delta time (dt). This step function steps all objects forward in time by integration. The integration step just moves the objects forward by using their velocity, position and acceleration to determine their next position.

During the step collisions are detected and then resolved. Often times a broadphase of some sort is used to detect possible collisions, and expensive collisions operations are only used when really needed.

All of this organization allows the user of your physics engine to worry about three main operations: creating and removing bodies, and stepping the scene. The rest is automated and taken care of within the physic system’s internals.

The last major isolated system would be the broadphase. There are two major phases in collision detection: the broad and narrow phases. The narrow phase is the one which determines if two shapes intersect or not. The broad phase determines if two shapes can possibly be intersecting or not. The broadphase should be constructed such that intersection queries are very very fast and simple. An axis-aligned bounding box (AABB) will suffice perfectly.

Once the broadphase chooses which objects could perhaps be colliding, it sends them off to the narrow phase. The narrow phase performs much more intensive operations to determine if two objects are colliding or not. The whole point of this is to reduce the amount of times the narrow phase has to be used, as it is expensive to operate with.

Once the narrow phase finds a pair of bodies to be colliding information about the collision is gathered and placed into a small container called the manifold. Do not ask why it is called a manifold, for I have no idea and neither does anyone else seem to! The manifold contains three important pieces of information:

- Collision penetration depth
- Direction to resolve in
- Points of contact

These pieces of information are used to fill in formulas that are used to resolve the penetration and solve for a single unknown: the magnitude of the resolution vector. Here’s a small diagram:

It is also useful to to store pointers or handles to the two objects that formed this collision info. This allows some useful functions to placed into the manifold object directly: solve and resolve. Solve would be the function to collect the three pieces of collision information. If no contact points are found, then there is no collision. If contact points are found, then resolving performs a resolution operation to force both objects to not be intersecting after integration.

# Velocity

Complex physics manipulation is performed on the velocities of objects. Trying to manipulate the positions of objects directly is very difficult, as it poses a problem that isn’t linear. By using derived position equations for velocity, the problem is thus simplified. Most of the time we will be only dealing with velocity manipulation.

# Impulse Resolution in 2D (No Rotation or Friction)

The act of resolving collisions is something that isn’t covered very well on the internet. There are some scattered documents and details, though the information isn’t always easy to find. Most of what I know I learned by talking with other people, but I know most people will not have such an opportunity. Hopefully this article series can provide a strong resource for understanding and constructing a simple physics engine.

The toughest place to code is in my opinion resolution of collision. There exists tons of information on collision detection and broadphase, and thus creating these portions of a physics engine is in my opinion not too difficult. Some resources for collision detection are: Real-Time Collision Detection by Christer Ericson, and Game Physics Engine Developement by Ian Millington. I have both of these books right next to me as I write this :)

Generating contact manifolds and resolving such manifolds are what most programmers get caught up in. So lets hit the ground running and tackle a portion of code that will bring your entire physics system to life: contact resolution.

The best type of contact resolution to start with is impulse resolution. The idea behind impulse resolution is to take your contact manifold and solve for a single velocity vector that can be used to add into an object’s current velocity, in order to make it not penetrating the other object during the next frame. The reason for starting with impulse resolution is that it’s quite simple and easy to get up and running, and more complicated and advanced techniques require you to understand impulse resolution anyway.

Now the contact manifold is assumed to contain the direction of our velocity vector we are solving for. I will cover how to generate the contact manifold in another article in this series. The only unknown left to solve for is the magnitude of this vector. It so happens that it’s possible to solve for this magnitude in relation to both objects. If this magnitude is known, you add velocity along the normal scaled by the magnitude for one object, and you do the same operation to the other object in the direction opposite to the manifold normal. Lets start from the top.

We have two objects moving along in our scene, and there exists a relative velocity between the two, with object A as the reference object at the origin:

1 2 |
Eq 1: VelocityRelative = VelocityA - VelocityB |

The relative velocity can be expressed in terms of the collision normal (from the collision manifold) with a dot product:

1 2 |
Eq 2: VelocityRelative dot ManifoldNormal = (VelocityA - VelocityB) dot ManifoldNormal |

This can be thought of as the relative normal velocity between the two objects. The next thing to include in this derivation is the coefficient of restitution (elasticity factor). Newton’s Law of Restitution states that:

1 2 |
Eq 3: VelocityBefore * restitution = VelocityAfter |

Often times the restitution identifier is specified by an e, or epsilon symbol. Knowing this it’s fairly simple to include it within our current equation:

1 2 |
Eq 4: VelocityRelativeAfter dot ManifoldNormal = -resitution(VelocityA - VelocityB) dot ManifoldNormal |

Now we need to go to another topic and model an impulse. I have said the term “impulse” quite a few times without defining it, so here is the definition I use: an impulse is an instantaneous change in velocity. We will use an impulse to change the velocity of an object. Here’s how you could use an impulse to modify the velocity of a single object:

1 2 |
Eq 5: VelocityNew = VelocityOld * Impulse |

Here Impulse would be a scalar floating point value. This isn’t too useful however, as it only scales an object’s current velocity, and thusly makes the object move slower or faster along the positive or negative direction of the vector.

What is needed is a way to do component-wise modification of the vector, so we can make it point slightly in one direction or another, allowing objects to make turns and change directions.

1 2 |
Eq 6: VelocityNew = VelocityOld + Impulse(Direction) |

In the above equation we can take a direction vector with a magnitude of 1, and scale it by our impulse. By adding this new vector to our velocity we can then modify the velocity in any way we wish. Our end goal is to solve for our Impulse scalar that will separate two objects from collision, so in the end we’ll need to distribute this scalar across two equations in terms of velocity.

Lets start with a simple momentum equation:

1 2 |
Eq 7: Momentum = Mass * Velocity |

An impulse is defined to be a change in momentum. Thus we get:

1 2 3 |
Eq 8: Impulse = MomentumAfter - MomentumBefore Impulse = Mass * VelocityAfter - Mass * VelocityBefore |

To isolate our velocity after we can rearrange into:

1 2 |
Eq 9: VelocityNew = VelocityOld + Impulse(Direction) / Mass |

Now lets change equation 4 into one that contains velocities under the influence of impulses. However we’ll want to express our VelocityNew as one that is acted upon by impulse, and substitute in equation 9:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
Eq 10: VelocityA + Impulse(Direction) / MassA - VelocityB + Impulse(Direction) / MassB = -Restitution(VelocityRelativeAtoB) * Direction Simplification: VelocityRelativeAtoB * Direction + Impulse( Direction / MassA + Direction / MassB ) + Restitution(VelocityRelativeAtoB) * Direction = 0 Isolate Impulse: Impulse = -(1 + Restitution)(VelocityRelativeAtoB dot Direction) ---------------------------------------------------- 1 + 1 ----- ----- MassA MassB |

Remember that the impulse is a scalar. Also note that all values on the right hand side of the equation are all known, including the Direction which was solved for by the collision detection.

All that is left here is to distribute this scalar impulse over each object’s velocity vector proportional to the total mass of the system (system being collision between both objects).

The total mass is MassA + MassB, so to get an even distribution you do: impulse * Mass / TotalMass. To simplify this one could use the following:

1 2 |
Eq 11 (applying an impulse proportional to mass of object): VelocityNew = VelocityOld + (1 / Mass) * Impulse(Direction) |

This can be done twice, once per object. The total impulse will be applied, except only a portion of the impulse will be applied to each object. This ensures smaller objects move more than larger ones during impact.

One thing you must ensure is that if the velocities are separating (objects moving away from one another) that you do nothing. Here’s a sample version of finalized code for impulse resolution in 2D without rotation or friction:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
// j = -(1 + epsilon) * N dot Vrel // ------------------------------- // invM_1 + invM_2 void Manifold::ApplyImpulse( void ) { Vec2 rv = B->m_velocity - A->m_velocity; real contactVel = Vec2::DotProduct( rv, normal ); // Do not resolve if velocities are separating if(contactVel > 0) return; // Calculate restitution real e = std::min(A->m_material.restitution, B->m_material.restitution); // Calculate impulse scalar real j = -(1.0f + e) * contactVel; j /= A->m_massdata.inv_mass + B->m_massdata.inv_mass; // Apply impulse Vec2 impulse = j * normal; A->m_velocity -= A->m_massdata.inv_mass * impulse; B->m_velocity += B->m_massdata.inv_mass * impulse; } |

# Rotation

Now that we have covered resolution without these two factors adding them in to our final equation (equation 10) will be quite a bit simpler. We will need to understand the concept of “mass” in terms of rotations. The inertia tensor of an object describes how “hard it is to turn” along some axis. In 2D there’s only a single axis of rotation (along the z) so a single tensor is all that’s needed. Here’s a new version of equation 11.

1 2 3 4 5 6 |
Eq 12: RotationalVelocityNew = RotationalVelocityOld + (1 / InertiaTensor) * PointRelativeToCenterofMass cross Direction * Impulse Shortened: V' += I^-1 * r x (n * j) |

For the sake of brevity here’s the final equation, where r is the vector from center of mass to a point on the body (contact point). The velocity of a point on the body relative to it’s center is given by:

1 2 |
Eq 13: Vpoint = VCenterMass + AngularVelocity cross r |

Final equation (Direction substituted for n):

1 2 3 4 5 6 |
Eq 14: Impulse = -(1 + Restitution) * (VelocityRelativeAtoB dot n) ------------------------------------------------- 1 1 (rA cross n)^2 (rB cross n)^2 ----- + ----- + -------------- + -------------- MassA MassB InertiaTensorA InertiaTensorB |

And there we have it! This will solve for a separating impulse given a specific contact point. Now you might have noticed there are a couple cross products that are kinda strange in 2D. I won’t cover what they are here, but they still hold true. I’ll just give them to you:

1 2 3 4 5 |
Cross vector and scalar in 2D: V' = x * V.y, -x * V.x Cross scalar and vector in 2D: V' = -x * V.y, x * V.x |

1 2 |
Cross two vectors in 2D: A.x * B.y - A.y * B.x |

# Friction

Friction is the simplest thing to do in this entire resolution article, assuming we’re dealing with 2D! I actually haven’t derived this portion myself, so again I’m going to just throw the equations at you. Assuming we have our resolution direction in the manifold (the collision normal), we can calculate a tangent vector, which is a vector perpendicular to the collision normal. Just replace all instances of n in equation 14 with this tangent vector.

To solve for the tangent vector you do:

1 2 |
TangentVelocity = VRelAtoB - (VRelAtoB dot n) * n TangentVector = Normalize( TangentVelocity ) |

Again, just take the above TangentVector and replace it in equation 14 for n.

There is something I have missed. When solving for the force of friction there’s a max, which is the manifold normal * the coefficient of friction (from the body definition). This is due to the Coloumb friction law. Treat the coefficient of friction * normal as your “friction cap” when solving for your tangential impulse. This is just a simple capping of the final impulse vector.

# Penetration Resolution

Due to floating point error energy will always be lost from the system over time. As such, when objects start coming to a rest they will sink into each other due to gravity. How to solve? Well the solution is to actually just handle the leftover penetration manually by just moving the objects apart directly. This will prevent objects from sinking into one another, though doesn’t add any energy into the system.

To resolve penetration, simply move each object along the manifold normal in opposite directions. However there’s an art to doing so; you need to be gentle. If you just resolve 100% of the penetration each frame objects underneath other objects will jitter violently due to the naive penetration resolution scheme I am presenting. To avoid this, try only resolving a percentage of the leftover penetration, perhaps 20 to 80%.

Additionally you only want to resolve penetration if the penetration depth after the impulse is applied is above some constant arbitrary (but small) threshold, (try 0.01 to 0.1). If the penetration is below this, then don’t move either object.

This method of penetration resolution is called linear projection. Here’s a snippet of C++ code demonstrating this:

1 2 3 4 5 6 7 8 |
void Manifold::PositionalCorrection( void ) { const real k_slop = 0.01f; // Penetration allowance const real percent = 0.2f; // Penetration percentage to correct Vec2 correction = (std::max( penetration - k_slop, 0.0f ) / (A->m_massdata.inv_mass + B->m_massdata.inv_mass)) * percent * normal; A->m_tx.position -= A->m_massdata.inv_mass * correction; B->m_tx.position += B->m_massdata.inv_mass * correction; } |

# Iterative Solving

There is one one additional tweak you can do to increase the believability of your physics simulation: iterate over all contacts and solve + resolve the impulses many times (perhaps 5 to 20 iterations). Since the large equation 14 has relative velocity within it, each iteration will feed in the previous result to come up with a new one.

This will allow the step to propagate energy throughout multiple objects contacting one another within a single timestep. This is essential for allowing the “billiards balls” effect to ensue.

This simple iteration is a very easy way to vastly improve the results of resolution.

# Resources/Links

- http://chrishecker.com/images/e/e7/Gdmphys3.pdf
- http://www.cs.cmu.edu/~baraff/sigcourse/notesd2.pdf

Vittorio RomeoHey, great article. Gamedevtuts brought me here.

I implemented my version of it and it mostly works, except the “Iterative Solving” part. Can you post your implementation of it?

Thanks.

Randy GaulPost authorThanks! Here, I hope this helps:

Vittorio RomeoThanks for the quick reply – unfortunately I’m still confused, as many of the methods’ implementations in your snippet are unknown to me. My code structure is very different and I can’t really correlate them.

Hope I’m not being too pretentious, but I’d really, really appreciate if you could post the entire engine source code (if that’s ok with you), either as an archive to download or a git (GitHub?) repository.

I’ve been working on a physics engine for quite a while and I would love too see how some things have been implemented in your own.

Randy GaulPost authorI sent you an email :)

JoshNot a bad intro to the topic! :) I feel I should point out, however, that there are some issues with your mathematical language. For example, your Eq. 12 shows division by InertiaTensor; you can’t do that. Furthermore, just above that you say that in 2D there is only one tensor since you can only rotate around the z axis. This is slightly sloppy language, since even in 3D you only ever deal with one inertia tensor per object. What you want to say instead is that because you only deal with rotation along one axis, the inertia tensor reduces to just a scalar value.

namuolThanks for the guide; lots of practical advice!

I think I may be missing a few pieces, though. I apply a position correction, but objects still sink (a better word would be

squeeze) into each other. They do notcontinueto sink — rather, the amount they “squeeze” into the object beneath them is proportional to the number of objects stacked on their top.Could it have something to do with the order of things? Here is my approach:

This

mostlyworks. It looks like I’ve got my impulse calculations right with zero gravity, and I get the same results with friction disabled.I’m wondering if the only way to achieve that nice rigid stacking behavior is to take an iterative approach, which I’m trying to avoid. Thoughts?

Thanks again!

Randy GaulPost authorThat approach will do decently well, but it’s best to iterate over solutions. The best professional physics engines are iteratively solving these sorts of physics problems over the course of many frames. It gets so important to iteratively solve, that previous calculation can be used to jumpstart future calculations. This is called “warm starting”.

For your physics engine, if you want better stack you’ll need to do some tricks in order to get good stacking. This is because impulse resolution is a pretty naive resolution approach, and linear projection (the penetration correction I talk about in this post) is even more naive.

Feel free to email me with more specific questions!

AkshayHey, thanks for this article!

I am only a third way into it and it has already helped me clear many confusions regarding architecture. This and qu3e are a good resource!

ErlendHi!

I have a problem understanding how you deduced Eq. 10 given Eq. 4 and 9. I suppose ‘Direction’ here is what was formerly ‘ManifoldNormal’?

Is ‘VelocityRelativeAtoB’ just a substitute for ‘VelocityA – VelocityB’, which would make the right side equal to earlier, or is it implied that also these velocities undergo an impulse?

I’m also wondering what happened to the ‘ManifoldNormal’ on the left side, and why VelocityB is negative, but not its respective impulse vector.

Thanks!

ErlendI found the answer to my confusions, I think. However, I still don’t know why VelocityB is negative. For anyone else being confused, here is a correction:

LawrenceHi Randy,

thanks for the cool blog.

In eq.10:

VelocityA + Impulse(Direction) / MassA – VelocityB + Impulse(Direction) / MassB = -Restitution(VelocityRelativeAtoB) * Direction

does VelocityA and VelocityB refer to the “old” velocity (before collision) and VelocityRelativeAtoB the “new” velocity (after collision)?

Cheers

Randy GaulPost authorAlmost. VelocityA is the velocity of body A. VelocityB is the velocity of body B. AtoB implies vector subtraction. To find a vector that goes from A to B we do: endpoint – startpoint. A to B means we start at A and end at B, so we do: B – A. So, VelocityRelativeAtoB is VelocityB – VelocityA.