Submit to StumbleUpon Share

Hello and welcome back to my blog!

This time i'm going to talk about the basic components that make up a physics engine and how to put them together; this tutorial is aimed at programmers who have a basic grasp of maths and geometry but would like to step into the world of simulation.

It is my hope that if, at the beginning of this article, you are able to code the physics behind the 1972 game Pong, by the end of the article you will be equally happy writing your own constraints to use in your own physics solver!

Can you code this?

I've always found the title of those books of from the '...for dummies' series reassuring; after all, if a dummy can learn this stuff you should stand a good chance, right?

...for dummies?

Hence the title of this article.

Plan of action

Ok, so i'm going to cover a few of the things you might want in a physics engine:

  • Rigid bodies
  • Collisions
  • Resting contact
  • Constraints (Joints)

Simulation

Bear with me, this is going to start out really basically, but hopefully later on it will become clear why.

Starting out with particles, which we assume will have a position P and a velocity V. Each frame we advance the position of P by adding the velocity V on to it. This is called integration.

You can use whatever units you like for you simulation; the usual choice is one unit/meter which is what i'm using here - the screen is two metres by two metres (in our world), so our velocities are specified in metres/second. To make this work, we must know how many seconds there are per frame, then we can do P += V*dt in our integration to advance the particles properly, where dt is the number of seconds per frame.

We can simulate multiple particles by storing an array of positions and velocities, and looping over integrating each one. To prevent our particles going off screen we would like to collide them against the screen's edges.

To bounce the particles off the edges we simply detect collisions against the edges and reverse the velocity in the offending axis, thus:

for all i
  if (P[i].x < -1 && V[i].x < 0)
  {
    V[i].x = -V[i].x;
  }
  if (P[i].y > 1 && V[i].y > 0)
  {
    V[i].y = -V[i].y;
  };
  if (P[i].x > 1 && V[i].x > 0)
  {
    V[i].x = -V[i].x;
  }
  if (P[i].y < -1 && V[i].y < 0)
  {
    V[i].y = -V[i].y;
  }

The first condition of each case checks for intersection of the particle in space and the second checks that the particle is actually heading towards the edge. The second condition is important otherwise we would do the collision response again next frame erronusly should the particles move more than one pixel per frame. This second condition carries though all of impulse based rigid body simulation and is what separates inequality contraints (e.g contacts) and equality constraints (most joints).

This is the kind of example code you will no doubt have written a long time ago when you first started getting interested in simulation. I'm including it here because I think there is a natural progression from this easy example to more complex physics code.

Whats going on in this simple simulation?

Well, without knowing it, we have assumed a physical matieral type for our particles including the coefficient of restituition and followed Newton's principle of conservation of momentum. We have also chosen that the world have infinite mass such that it doesn't feel a collision response when the particles impact on it. By choosing to reflect particles on impact we are also preseving their momentum (even though we have ignored their mass in our calculations), indicating that we have chosen a coefficient of restitution of 1, i.e. perfectly elastic - like a super-ball. Additionally, we have chosen an impulse/velocity model for the collision response, rather than a force/acceleration model - we did this simply by changing the velocity instantaneously, rather than over a period of time.

What we have here is actually a highly optimised special case physics simulation. Optimised how, you ask? Well, let me explain:

If we were to code the above "properly", not taking any short cuts, we would have to assume the following. Our environment is defined by its boundaries which are four axis aligned planes (actually, they're lines because we are in 2d). Each one has a normal (which points inwards) and a distance to the origin. They look like this:

Planes[4] =
{
   (1,  0, 1),
   (0, -1, 1),
   (-1, 0, 1),
   (0,  1, 1)
}

So, now we have to detect our collisions as we did before, but we can't take any shortcuts this time. In order to detect if our particles have penetrated the planes we must perform a dot product and add the plane's distance to origin:

for all particles i
{
  for all planes j
  {
    distance = P[i].x*Planes[j].a + P[i].y*Planes[j].b + Planes[j].c;
 
    if (distance < 0)
    {
      // collision responce
    }
  }
}

What this code is doing is finding, by projection, how much of the vector from the plane to the particle is in the direction of the plane normal. Because our plane's normals are unit length, this is a measure of the closest distance from the particle to the plane. Obviously, if this distance is less than 0, our particle has penetrated the plane and we need to take action to perform the collision response.

Now, if we take a closer look at that distance check above, including each plane's coefficents:

plane0dist = P[i].x*1  + P[i].y*0  + 1;
plane1dist = P[i].x*0  + P[i].y*-1 + 1;
plane2dist = P[i].x*-1 + P[i].y*0  + 1;
plane3dist = P[i].x*0  + P[i].y*1  + 1;

If we simplify that down a bit:

plane0dist =   P[i].x + 1;
plane1dist =  -P[i].y + 1;
plane2dist =  -P[i].x + 1;
plane3dist =   P[i].y + 1;

Re-arranging slightly we get these plane tests conditions:

if (P[i].x  < -1)
if (-P[i].y < -1) = if (P[i].y > 1)
if (-P[i].x < -1) = if (P[i].x > 1)
if (P[i].y  < -1)

Which are exactly the same as those which we used in our first simple example. The second condition for each plane must also be satisfied so we only do the collision response once per collision. This can be done by performing a dot product between the particle's velocity and the normal for each plane;

if (P[i].x < -1 && V[i]•N[i] < 0)

We call this the normal velocity since its the particle's velocity in the direction of the normal. If this results in a value less than zero then the particle is moving towards the plane and we allow the collision. Once the collision is resolved, the normal velocity will be >= 0 depending on the coefficient of restitution. I could perform the same analysis on the second condition as i did the first to prove that this general solution is the same as the first specific version, but i will leave that as an excercise for the reader.

This shows our two approaches are physically the same.

Ok, so what do we do now for the collision response?

We need a solution which achieves the same result as the original program but is more general. We can do this using the reflection vector from the law of reflection. The reflection vector is calculated thus:

R = V - 2*N*(VN)

Where V is the velocity vector and N the surface normal. We need to do a similar comparison with this operator as we did with the collision checks:

plane0vel x = V.x - 2* 1*(V.x* 1 + V.y* 0) = V.x -2*V.x   =   -V.x
plane0vel y = V.y - 2* 0*(V.x* 1 + V.y* 0) = V.y - 0      =    V.y
plane1vel x = V.x - 2* 0*(V.x* 0 + V.y*-1) = V.x - 0      =    V.x
plane1vel y = V.y - 2*-1*(V.x* 0 + V.y*-1) = V.y - 2*V.y  =   -V.y
plane2vel x = V.x - 2*-1*(V.x*-1 + V.y* 0) = V.x - 2*V.x  =   -V.x
plane2vel y = V.y - 2* 0*(V.x*-1 + V.y* 0)                =    V.y
plane3vel x = V.x - 2* 0*(V.x* 0 + V.y* 1)                =    V.x
plane3vel y = V.y - 2* 1*(V.x* 0 + V.y* 1) = V.y - 2*V.y  =   -V.y

Which, you should be able to see, are the exact same resulting velocities after collision with each plane in the first simple example.

You will have noticed that mass has been ignored in our calculations so far, this is because it just drops right out of the equations when you are dealing with collision against a body of infinite mass (i.e. our simple world).

So what does this new version look like so far in pseudo-code?

for all particles i
{
  for all planes j
  {
    N = {Planes.a, Planes.b};
    distance = P[i]•N + Planes[j].c;
 
    if (distance < 0 && V[i]•N < 0)
    {
      // collision response, reflect particle
      V[i] -= 2*N*N•V[i];
    }
  }
}

Great, so we've turned a nice simple example into one that's far more complicated than it needs to be. Why?

Firstly, to demonstrate that the roots of all our inital assumptions were grounded in physics and secondly, to demonstrate the advantages of the more complex implementation. Because we now have a way to deal with arbitrary 2d planes instead of axis aligned planes, it means we can rotate our world and still have the simulation work correctly:

Move the mouse over the demo to rotate the planes. Note that in this demo, i have done some correction to the points as the planes are rotated to make sure the points still lie within the planes.

Ok, so this is all very well but its not a very real world simulation. There is no gravity and the coefficient of restistution is far too high - in the real world, most everyday things have a near plastic coefficient of restitution, i.e near 0.

So, adding gravity is pretty simple, we just subtract gravity from our particle's velocities before we do any collision, remembering to account for the fact that gravity is an acceleration:

V[i] += G*dt

In this case G is the vector {0, -9.8} and dt is the number of seconds per frame as before.

We would also like to use a different coefficient of restitution as well because our particles might not be made of super-rubber. Recall the reflection vector equation that we used earlier:

R = V - 2*N*(VN)

We can actually re-write this equation to include the coefficient of restitution:

R = V - (1+e)*N*(VN)

Where e is the coefficient of restitution which varies from 0 (totally plastic) to 1 (totally elastic). You can see that these two equations are equivelent; if we set e to 1, we get the original equation. The particles in this example have coefficients of restitution varying from 0 to 1.

The only thing we can now do to this particle simulation to make it more realistic (without going to a full rigid body simulation with rotation) is to handle masses in our collisions. To accomplish this, we will have to use circles instead of particles, since its quite unlikely that two particles will ever collide what with them having no actual size!

Figure 1

As per Figure 1, two circles A and B intersect if the distance between them is less than the sum of their radii. The collision normal is simply the vector d between the two, but normalised.

In order to deal with the masses of our particles we need to do some thinking regarding the reflection equation that we already have. What we really need is some kind of mass ratio to apply a bias to the equation so that lighter particles get pushed out of the way by heavier particles, but we need it to conserve momentum as well.

ratioa = Mb / (Ma + Mb)
ratiob = Ma / (Ma + Mb)

The above two ratios will do the trick. Here is an example:

Ma = 1.0
Mb = 0.5

ratioa = 0.5 / (1.0 + 0.5) = 1/3
ratiob = 1.0 / (1.0 + 0.5) = 2/3

Obviously 1/3 + 2/3 = 1, so you can see that this conserves momentum correctly, and since body b is half the mass of body a, it should experience twice the reaction upon collision, which it does.

So we can plug these two ratios into our reflection vector equation to obtain two new equations for the reflected velocity of each body after collision:

Va = Va - (1+e)*N*((Vb-Va) N)*(Mb / (Ma+Mb))
Vb = Vb - (1+e)*-N*((Vb-Va)
-N)*(Ma / (Ma+Mb))

Our collision code gives us the normal pointing from body b towards body a, so we must reverse the normal to calculate body b's reflected velocity. You will have noticed that there are a lot of duplicated calculations in the two equations above. We would like to simplify this down a bit. What we can do is use the relative velocity of the two objects in question:

Vr = Va - Vb

Now that we have only one velocity to deal with the equation gets simpiler.

I = (1+e)*N*(Vr N)

Va = -I*(Mb / (Ma+Mb))
Vb = +I*(Ma / (Ma+Mb))

What we are doing here is to treat one object as being stationary relative to the other moving object. But there is still more we can do here to simplify these equations:

I = (1+e)*N*(Vr N)*(Ma*Mb)/( Ma+Mb)

Va - = I /  Ma
Vb + = I / Mb

By combining the ratios we can save some more calculation and reduce the final step to a simple divide though by the mass of the objects in question to produce the final velocities for those objects. You can see this works because

(Ma*Mb)/( Ma+Mb) / Ma Mb/( Ma+Mb)

and

(Ma*Mb)/( Ma+Mb) / Mb = Ma / (Ma+Mb)

Finally we can show that (from fractions)

(Ma+Mb)/( Ma*Mb) = 1/Ma + 1/Mb

Which means we can re-write the final linear impulse equation

I = (1+e)*N*(Vr N) / (1/Ma + 1/Mb)

This is handy because we can just store 1/Ma and 1/Mb inside the definition of our circle rigid bodies so we don't have to re-calculate. It also means that we now have a way to represent infinitly massive objects, such as our world (by storing the inverse mass as 0).

Va - = I * 1/Ma
Vb + = I * 1/Mb

In these demos you can manipulate the circles with the mouse.

Ok, so you probably already noticed we are starting to see problems due to penetration caused by the collision detection system detecting the collision after its already too late. This problem only gets more severe as we add more rigid-bodies, and is compounded by the fact that the system is currently doing nothing to correct the penetration:

To combat this problem, we need to be able to detect collisions before they actually happen and deal with penetration. Luckily i've covered this technique already in a previous article which i suggest you read.

Once this technique is applied the problem just goes away; and there are only a few lines of code required. The only caveat is that we have to assume that the coefficient of restitution is never required to be non-0; in the majority of real-world cases this isn't too much of a limitation, especially given the benefits.

Software engineering

In physics engines, the mathsy part is only half the story; the second half is how you organise your physics engine in terms of classes hierarchy and how you actually code it.

I'm going to present one possible solution here which has worked for me in the past; it may be over-simplified for your requirements but it should give you an idea. The language here is actionscript, but the principles apply to any OO language.

public class RigidBody
{
	protected var m_pos:Vector2;
	protected var m_vel:Vector2;
	protected var m_invMass:Number;
 
        ///
	///
	///
	public function RigidBody( pos:Vector2, vel:Vector2, invMass:Number )
	{
	        m_pos=pos;
		m_vel=vel;
		m_invMass=invMass;
	}
 
        ///
	///
	///
	public function GenerateContact( rb:RigidBody ):Contact
	{
		throw new Error("Not implemented on base class");
	}
 
	///
	///
	///
	public function Integrate( dt:Number ):void
	{
		m_pos=m_pos.Add( m_vel.MulScalar( dt ) );
	}
}

So, i've defined a RigidBody base class containing the three parameters that we need so far in our example. There is an integrate function to move foward in time and (what should be a virtual) function called GenerateContact() which will generate a collision normal and distance between any two RigidBodys.

public class Circle extends RigidBody
{
	private var m_radius:Number;
 
        ///
	///
	///
	public function Circle( pos:Vector2, radius:Number, invMass:Number )
	{
		m_radius = radius;
		super(pos, new Vector2(), invMass);
	}
 
        ///
	///
	///
	public override function GenerateContact( rb:RigidBody ) :Contact
	{
		if ( rb is Circle )
		{
			...
		}
                else if (rb is ...)
                {
                        ...
                }
                else
		{
		    throw new Error("unahandled case!");
                }
	}
}

And we have a Circle which derives from RigidBody and extends its functionality by having a radius parameter. Also, we have implemented the GenerateContact() function to return the required information when needed.

By separating functionality out into base class and concrete implementations we will be able to deal with multiple RigidBodys in the same lists, which greatly simplifies the code.

public class Plane extends RigidBody
{
	private var m_n:Vector2;
	private var m_d:Number
 
	public function Plane(n:Vector2, d:Number )
	{
		m_n=n;
		m_d=d;
 
		super(n.MulScalar(-d), new Vector2(), 0);
	}
 
        ///
	///
	///
	public override function GenerateContact( rb:RigidBody ):Contact
	{
		if ( rb is Particle )
		{
			...
		}
		else if ( rb is Circle )
		{
			...
		}
	        else
                {
			throw new Error("unhandled case!");
                }
	}

Above is another concrete implementation, this time representing an infinite plane (our screen edges).

public class Contact
{
	public var m_normal:Vector2;
	public var m_dist:Number;
 
	///
	///
	///
	public function Contact( n:Vector2, dist:Number )
	{
		m_normal = n;
		m_dist = dist;
	}
}

Above is the definition for the Contact class, which represents the information the system will need to resolve one collision.

private function Update( e:GameLoopEvent ):void
{
	const dt:Number = Math.min(e.m_elapsed, 1.0/15.0);
 
	// apply gravity
	for each ( var p:RigidBody in m_rigidBodies )
	{
		if ( p.m_InvMass>0 )
		{
			p.m_vel=p.m_vel.Add( Constants.kGravity.MulScalar( dt ) );
		}
	}
 
	// collide
	for ( var i:int=0; i<0||rbj.m_InvMass>0 )
	{
		const c:Contact=rbi.GenerateContact( rbj );
 
		...
 
		//resolve collision
	}
 
	// integrate
	for each ( p in m_rigidBodies )
	{
		p.Integrate( dt );
	}
}

Above is what our update loop currently looks like. The order here is fairly important;

  1. Firstly, external forces are applied - i.e. gravity.
  2. Then, all the collision detection is done to generate a contact which needs to be resolved for each pair of interacting rigid bodies.
  3. Finally, each RigidBody is integrated forward in time.

The time-step clamping at the beginning is so that as we debug the code, we don't end up with such a massive time-step the next frame that everything explodes!

Constraints

Ok, so we have a pile of spheres up and running and our physics engine is starting to take shape quite nicely, what about constraints?

Before i talk about this, i feel i should introduce some concepts needed to understand the problem.

Constraint

Loosely, this is something which enforces a condition between two rigid bodies. So, when there is a collision we create a constraint which enforces the rule that the colliding rigid body may not move through the object it collides with but also that it do this by only pushing. We already derrived the code for this above (not taking into account rotation yet):

I = (1+e)*N*(Vr N) / (1/Ma + 1/Mb)

All the constraints i'm going to work with are impulse velocity level constraints. Forces and accellerations will not come into it.

Inequality constraint

This is a constraint which only acts in one direction. So a collision constraint is an inequality constraint because it is only allowed to push and never pull. If it were allowed to pull, it would stick the rigid body to the object it collided with. The way we enforce this inequality is simply the check i mentioned at the beginning:

if (V•N < 0)
{
   // handle constaint
}

Little Big Planet contained a lot of constraints of this type; winch, piston, string are a few examples - they only allowed motion in one direction within the limits of the constraint. So a piston, for example was a length constraint with lower and upper bounds on the allowed length.

Equality constraint

This constraint type is allowed to pull and push equally. An example is a point to point constraint, where two rigid bodies are connected by a single point - one point on each body is forced to lie in the same place. With a hinge constraint the whole hinge axis is constrained to lie in the same place on both rigid bodies. A rod constraint would be a a length constraint where the length is never allowed to change.

Designing a constraint

The way to approach how to design a particular constraint is to think about what restricting effect it has on the two rigid-bodies its connected between.

For example, a distance constraint (a rod in LBP) simply prevents the two end points of the rod getting any closer than (or any further away from) the specified distance. Putting the end points exactly at the centre of mass of each rigid body, as we are about to describe keeps simple circular rigid bodies at a constant distance from each other, while still allowing both of them to continue to move around in 2d space.

Figure 2

In Figure 2, bodies A and B must not get any closer than distance d, but are still allowed to move around independently otherwise - this resolves as each body being able to orbit the other as they fly through space.

Impulse/velocity level

As mentioned previously, the simulator i'm describing is an impulse/velocity level simulator, rather than a force/acceleration level one. This means that all the constraints must only be solved by applying impulses to change velocities of bodies.

What this means when designing a constraint is that once you've worked out how it behaves on a position/distance level (i.e. on paper, with a diagram), you need to then think about how that translates into velocities.

So, to recap: our distance constraint stops the end points getting any closer or further away from each other. In velocity terms this is describing a one dimensional constraint; the relative velocity of the bodies is constrained to be zero in one particular axis. That axis is defined by vector between the end points.

Figure 3

In Figure 3, A and B have starting velocities which violate the length constraint that we've placed between them (from Figure 2). The first step in solving a constraint is to figure out where that velocity actually is. In this case you can see its the projection of Av and Bv on to the length constraint axis.

Once we have the velocities we want to remove, we convert into a single, relative velocity (as we did before with the contact constraint) and we already derived the maths required to solve one dimensional constraints (from the contact constraint):

I = (1+e)*N*(Vr • N) / (1/Ma + 1/Mb)

Lets re-write that so it becomes less of a mouth-full:

I = RelativeVelocityMagnitudeToRemove / (1/Ma + 1/Mb)

In addition to the velocity part, there will also be some positional fix-up to do, since with more than one constraint acting in the system simultaneously it will often fail to be resolved with one one solver iteration. In this simple example we just add an artificial correction to the velocity (current length - desired length) / time-step.

Above is the result - you can manipulate the circles with the mouse as before, although one is fixed in place to demonstrate the constraint chain.

Software engineering

So, lets cover the additions to the code.

public class Constraint
{
	protected var m_bodyA:RigidBody;
	protected var m_bodyB:RigidBody;
 
	/// 
	///
	/// 
	public function Constraint( bodyA:RigidBody, bodyB:RigidBody )
	{
		m_bodyA = bodyA;
		m_bodyB = bodyB;
 
		Assert( m_bodyA.m_InvMass>0||m_bodyB.m_InvMass>0, "Constraint between two infinite mass bodies not allowed" );
	}
 
	/// 
	///
	/// 
	public function ApplyImpulse( I:Vector2 ):void
	{
		m_bodyA.m_vel = m_bodyA.m_vel.Add( I.MulScalar(m_bodyA.m_InvMass) );
		m_bodyB.m_vel = m_bodyB.m_vel.Sub( I.MulScalar(m_bodyB.m_InvMass) );
	}
 
	/// 
	///
	/// 
	public function Solve( dt:Number ):void
	{
		throw new Error("base class doesn't implement!");
	}
}

Above you can see we've added another base class; this one defines the basic part of a constraint but doesn't actually define any concrete details - as before, we leave this to the derived classes. It does however provide a framework for us - each constraint must implement the Solve() function which does the meat of the work. Additionally, there is an ApplyImpulse() function defined to save us some work in the concrete implementations.

public class DistanceConstraint extends Constraint
{
	private var m_distance:Number;
 
	/// 
	///
	/// 
	public function DistanceConstraint( bodyA:RigidBody, bodyB:RigidBody, distance:Number )
	{
		super(bodyA, bodyB);
 
		m_distance = distance;
	}
 
	/// 
	///
	/// 
	public override function Solve( dt:Number ):void
	{
		// get some information that we need
		const axis:Vector2 = m_bodyB.m_Pos.Sub(m_bodyA.m_Pos);
		const currentDistance:Number = axis.m_Len;
		const unitAxis:Vector2 = axis.MulScalar(1/currentDistance);
 
		// calculate relative velocity in the axis, we want to remove this
		const relVel:Number = m_bodyB.m_vel.Sub(m_bodyA.m_vel).Dot(unitAxis);
 
		const relDist:Number = currentDistance-m_distance;
 
		// calculate impulse to solve
		const remove:Number = relVel+relDist/dt;
		const impulse:Number = remove / (m_bodyA.m_InvMass + m_bodyB.m_InvMass);
 
		// generate impulse vector
		const I:Vector2 = unitAxis.MulScalar(impulse);
 
		// apply
		ApplyImpulse(I);
	}
}

Above is the implementation for the actual distance constraint, it does the work as described above.

private var m_joints:Vector.<Constraint>;

In the solver we define a list of constraints - they are different to contacts in that they're more permanent, so they get their own list.

const dt:Number = Math.min(e.m_elapsed, 1.0/15.0);
 
// apply gravity
for each ( var p:RigidBody in m_rigidBodies )
{
	...
}
 
// collide
for ( var i:int=0; i<m_rigidBodies.length-1; i++ )
{
	const rbi:RigidBody=m_rigidBodies[i];
	for ( var j:int=i+1; j<m_rigidBodies.length; j++ )
	{
		const rbj:RigidBody=m_rigidBodies[j];
 
		if ( rbi.m_InvMass>0||rbj.m_InvMass>0 )
		{
			...
		}
	}
}
 
// solve constraints
for each(var jt:Constraint in m_joints)
{
	jt.Solve(dt);
} 
 
// integrate
for each ( p in m_rigidBodies )
{
	p.Integrate( dt );
}

The solver loop now looks like the above - we have gained a little loop for solving the constraints - i've chosen to do that after solving for contacts which makes the joints 'stronger' than the contacts (in that because they're resolved last, they will be more correct than the contacts which are done first), you can do it the other way around if you like.

Conclusion

This article has uncovered the physics lurking behind the most simple of 2d games (pong), shown that the principles therein are grounded in physics and expanded on the technique used in that game step by step all the way up to the point where designing a constraint in a physics engine can be seen as a simple extension of colliding two objects.

I hope this article has whet your appetite for physics simulation while also being easy to understand. I have run out of time with this article, but depending on its popularity i can write a lot more on the subject. Register your interest by posting in the comments! 🙂

Source code

As ever, if you enjoyed reading this article and would like to see more, please consider buying the source code which accompanies it; this will allow an indie developer like me to pay the rent and buy more food to eat! 🙂

You will receive all the framework code that i've laid out in this article, including the code behind the demos on this page which should serve as a good base to expand upon. They are written in actionscript 3.0 but the techniques are applicable to all OO languages.

After purchasing, you will be redirected to a page where you can download the source immediately.

USD 14.99

Subscribers can access the source here

Until next time, have fun!

Cheers, Paul.

Submit to StumbleUpon Share