Building a 2D physics engine, Part 1 – Basic motion


Let’s kick off this series with some basics, and walk through the design & approach used. I’m assuming that everyone is familiar with the math concepts I’ve covered in the past, so I won’t dwell too heavily on the details of the math or solving the equations, except where it’s new. I will, however, introduce most concepts in their mathematical form first so we understand the relationships between various quantities, and then discuss how they appear in code (which is often much simpler than the math version). But first, let’s start with an overview of the physics engine we’re building.

We’ll be modeling rigid bodies, which are the individual objects that interact in the world. Some of them move about, while others are stationary (think the walls or floors of a level). Each has some shape, mass, and other physical properties, and they cannot pass through each other. They are called rigid since their individual shape & composition doesn’t change over time (versus soft bodies, such as cloth, fluids, gel, etc… which are deformable). Together, these bodies make up a simulation system. The simulation system is just the grouping of these objects along with some additional properties such as gravity and time. The simulation is the component that advances time and updates the bodies.

As each object moves about, it may encounter other objects. These encounters may result in a collision, which will need to be resolved in a realistic manner. This may mean the objects bounce off of each other, affecting their velocities. It might mean an object slides down a slope, affected by friction. It could be many things, and we’ll model a basic system that can be used to simulate several interesting effects. We will have code to detect such collisions, and return some details in the form of contact data. The contact data will then be used to resolve the collision in some way.

For this simple tutorial, we’ll be focused on basic concepts, and so I’m going to keep the code and the algorithms very naïve, avoiding obfuscating optimizations. When we move on to the 3D physics engine, since the basics will have already been covered, we’ll explore and tackle more advanced topics like broad phase pair detection, more exotic collision detection methods & contact point generation, continuous collision detection, time of impact, joints, and so on.

Linear Newtonian Dynamics

The first thing we need to understand, before any of the cool stuff we just talked about can happen, is how to model and move our objects. We will use Newton’s basic laws of motion, which may sound familiar from high school physics courses, to model our simulation.

First, let’s discuss the individual objects themselves. If we assume our object is just a particle with no shape or orientation, we can model it simply as a position with a mass. This position, since we’re doing a 2D physics engine, is just a 2D point. Let’s write it as p. It’s written in bold since it’s a vector quantity with 2 components, x and y. The mass is a simple scalar, let’s call it m. In code, this might look like:

Vector2 position;
float mass;

We’ll cover shapes in the next post, and rotation & angular acceleration in a later post. Therefore, for the remainder of this discussion we’ll assume we’re talking about a simple particle with a mass.

We know that moving objects have velocity. What is velocity? Well, it describes the change in position over time. The rate of change of one variable as a result of a change of another variable is the core principle of Calculus, and is called a derivative. Stated using Calculus terminology, we’d say that velocity is the first derivative of position with respect to time. This means that it is a measure of how much the position changes as a result of a change in time. This is often written as:

\boldsymbol{v}=\boldsymbol{p'}

The prime mark (‘) on the p represents first derivative with respect to time. This says the velocity, v, is the first derivative of position p with respect to time. So, if we have a velocity and a change in time, how do we determine how much it changes the position? We need to take the reverse of the derivative in order to solve it. The reverse of a derivative is called an anti-derivative, or integral. Finding and using the integral to solve for change in position in this way is called integration.

The integral for the velocity above would be written as:

\int \boldsymbol(v)dt

Which means the integral of the velocity, for a change in time of dt. dt is sometimes written as ∆t to mean “change in time” or “delta time”. Solving the above equation yields:

\boldsymbol{p}=\boldsymbol{v}\Delta t + \boldsymbol{c}

That is, the new position is equal to the velocity multiplied by the change in time, added to some constant c. The constant, in this case, is the starting position before the change in time occurred. So:

\boldsymbol{p}=\boldsymbol{v}\Delta t + {\boldsymbol{p}}_{0}

Armed with this, we can now write code to update our position, given the velocity and a change in time:

float dt; // Change in time this frame, in seconds
Vector2 velocity;

position = velocity * dt + position;

Or, more simply:

position += velocity * dt;

Great, so we know how to update our position for a given change in time, if we have a velocity. How do we update our velocity? Well, again we can state this as the first derivative of velocity with respect to time. We call this quantity acceleration.

\boldsymbol{a}=\boldsymbol{v'}=\boldsymbol{p''}

That is, acceleration is equal to the first derivative of velocity (with respect to time), or the second derivative of position. We integrate the same way as above to update the velocity using acceleration:

\boldsymbol{v}=\boldsymbol{a}\Delta t + {\boldsymbol{v}}_{0}

And in code form:

Vector2 acceleration;

velocity += acceleration * dt;

So, given an acceleration, we can update our velocity. And given that updated velocity, we can update our position. The last piece remaining then, is to find a way to generate accelerations. Of course, we could assign them directly, but it’s hard to know exactly what acceleration should be applied to something. One exception to this is gravity, since we know that gravity is an equal acceleration felt by all objects, so we can just apply that one directly. In the real world, accelerations are the result of another concept called forces, and by modeling forces for the rest of our interactions, we get more realistic accelerations as a result.

The 3 laws of motion from Newton state:

  1. An object remains at rest, or continues to move at constant velocity, until acted on by an external force.
  2. The amount of acceleration affecting an object is proportional to the force applied divided by the object’s mass.
  3. When an object exerts a force on another object, the other object exerts and equivalent force in the opposite direction back.

We can use the second law above to come up with this equation:

\boldsymbol{F}=m\boldsymbol{a}

Rearranged to:

\boldsymbol{a}=\frac{\boldsymbol{F}}{m}

This means that given a force on an object, we can determine exactly what acceleration is generated based on the object’s mass. We end up encountering a division by mass often enough, that it’s usually best to just store the inverse of the mass directly so we can just multiply later. We now have all the basic equations we need in order to do simple simulations of our objects.

The code

As a reminder, all of the source code for this tutorial is maintained in this github repository. The code is not yet complete at the time of this writing, but it will always be at least up to or ahead of the current blog post being published. Meaning, it has everything in this blog post (and more) in it, but it may not be the final code with all of the features until we actually get through those posts. Therefore, if you run it and it doesn’t look complete or contains stability issues, this is the reason. Please feel free to report any bugs you find in the issue tracker on Github, or via comments to the posts. If it’s not something I’m already aware of and fixing, I’ll be sure to get it fixed soon.

The sample is broken up into the following components:

    1. Main.cpp – Boilerplate Windows application level stuff, like implementing the main entry point and starting up a sample PhysicsWorld instance, throwing some objects into it, and updating it each frame.
    2. DebugRenderer*.* – A simple debug rendering facility to render shapes with. This is used to render some simple shapes representing the objects in the simulation, and some debug features like contact points, etc…
    3. Vector2.h, Matrix2.h – Minimal math code we need to implement this sample.
    4. PhysicsWorld.h/cpp – Implements a single simulation. This allows configuring and inserting objects into it, setting global attributes like gravity and iteration counts, and updating the simulation each frame with an update call.
    5. RigidBody.h/cpp – The rigid body objects we’re simulating. This is the object that contains the properties we discussed in this blog post, for instance.
    6. RigidBodyPair.h/cpp – Represents a pair of 2 interacting bodies. As we’ll see in later posts, this will track the contact (if any) between two objects.
    7. Shape.h – Shape contains the implementations of the various shapes we can simulate in this sample.
      Collision.cpp – Implements the simple collision detection routines we use for this sample.

In this post, we encountered the following pieces:

From RigidBody.h:

// A rigid body is a dynamic object which can be affected by forces and constraints
class RigidBody
{
public:
    // The rigid body takes over the shape's lifetime.
    // When the rigid body is destroyed, it will delete the shape
    RigidBody(Shape* shape, float mass);
    ~RigidBody();

    const Shape* GetShape() const { return _shape; }

    const Vector2& Position() const { return _position; }
    Vector2& Position() { return _position; }

    const Vector2& LinearVelocity() const { return _linearVelocity; }
    Vector2& LinearVelocity() { return _linearVelocity; }

    const Vector2& Force() const { return _force; }
    Vector2& Force() { return _force; }

    const float Mass() const { return _mass; }
    const float InvMass() const { return _invMass; }

    const float Rotation() const { return _rotation; }
    float& Rotation() { return _rotation; }

    const float AngularVelocity() const { return _angularVelocity; }
    float& AngularVelocity() { return _angularVelocity; }

    const float Torque() const { return _torque; }
    float& Torque() { return _torque; }

    const float I() const { return _I; }
    const float InvI() const { return _invI; }

private:
    Shape* _shape;

    // Linear
    Vector2 _position;
    Vector2 _linearVelocity;
    Vector2 _force;
    float _mass, _invMass;

    // Angular
    float _rotation;
    float _angularVelocity;
    float _torque;
    float _I, _invI;
};

From PhysicsWorld.cpp:

void PhysicsWorld::Update(float dt)
{
    float invDt = dt > 0.0f ? 1.0f / dt : 0.0f;

    // Determine overlapping bodies and update contact points
    UpdatePairs();

    // Integrate forces to obtain updated velocities
    for (auto& body : _bodies)
    {
        if (body->InvMass() == 0.0f)
            continue;

        body->LinearVelocity() += dt * (_gravity + body->InvMass() * body->Force());
        body->AngularVelocity() += dt * body->InvI() * body->Torque();

        // Dampen the velocities to simulate friction (we'll add friction simulation later)
        static const float DampeningTerm = 0.001f;
        body->LinearVelocity() += -body->LinearVelocity().Normalized() * DampeningTerm;
        body->AngularVelocity() += (body->AngularVelocity() > 0) ? -DampeningTerm : DampeningTerm;

        // Clamp to 0 if the value becomes too low
        static const float ClampThreshold = 0.01f;
        if (body->LinearVelocity().LengthSq() < ClampThreshold)
        {
            body->LinearVelocity() = Vector2(0, 0);
        }
        if (fabsf(body->AngularVelocity()) < ClampThreshold)
        {
            body->AngularVelocity() = 0.0f;
        }
    }

    // Do all one time init for the pairs
    for (auto& pair : _pairs)
    {
        pair.second.PreSolve(invDt);
    }

    // Sequential Impulse (SI) loop. See Erin Catto's GDC slides for SI info
    for (int i = 0; i < _maxIterations; ++i)
    {
        for (auto& pair : _pairs)
        {
            pair.second.Solve();
        }
    }

    // Integrate new velocities to obtain final state vector (position, rotation).
    // Also clear out any forces in preparation for the next frame
    for (auto& body : _bodies)
    {
        if (body->InvMass() == 0.0f)
            continue;

        body->Position() += dt * body->LinearVelocity();
        body->Rotation() += dt * body->AngularVelocity();

        body->Force() = Vector2(0, 0);
        body->Torque() = 0.0f;
    }
}

For now, we can ignore the non-highlighted portions, which implement angular motion, collision detection, and collision response. Other than some artificial placeholder drag (we’ll implement friction later), the code implements exactly what we discussed in the post.

That’s it for this post. I want to keep each post scoped to make them easy to follow. If you have questions or want me to drill into something deeper, please just let me know in the comments. In the next post, we’ll start to discuss shapes and basic collision detection!

Advertisements
This entry was posted in Math, Physics, Physics & Collision Detection, Tutorial and tagged , , , , , , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s