My Rigid Body Simulator

As you may know, last year I completed AP Physics I and II. But that doesn't really mean much. The real question is, did I learn anything?

There's a couple ways to find out, but a rigid body simulator had been on my mind for a long time, and hey, why not?

The Objects

Each object is a collection of points connected by massless beams. The size of the dot is indicative of the mass of each point. Each object is assigned a velocity and angular velocity:

function createObject(points, velocity, angularVelocity, flags){
    var object = {
        points:points,
        velocity:velocity,
        angularVelocity:angularVelocity
    };

    setObjectMass(object);
    setCenterOfMass(object);
    setMomentofInertia(object, object.centerOfMass);

    object.fixed = (flags !== undefined && flags.indexOf('fixed') !== -1);
    objects.push(object);
    return object;
}

createObject([
    {x:-500,y:750, mass:9e99},
    {x:900,y:750, mass:9e99},
    {x:900,y:900, mass:9e99},
    {x:-100,y:900, mass:9e99}
    ], {x: 0, y: 0}, 0, 'fixed');
    

This snippet creates a block with arbitrarily large mass – used as a border.

Finding the center of mass and moment of inertia is easy enough if we treat an object as a collection of point masses:

//center of mass = weighted average of points
function setCenterOfMass(object){
    var sumX = 0, sumY = 0;
    for(var p = 0; p < object.points.length; p++){
        sumX += object.points[p].x * object.points[p].mass;
        sumY += object.points[p].y * object.points[p].mass;
    }
    object.centerOfMass = {x: sumX / object.mass, y: sumY / object.mass};
    }

//moment of inertia is equal to sum(mr^2). Note the axis will probably always be 
//the center of mass
//axisOfRotation is a point, not an actual axis. 
function setMomentofInertia(object, axisOfRotation){
    object.momentOfInertia = 0;
    object.axisOfRotation = axisOfRotation;
    for(var p = 0; p < object.points.length; p++){
        var r = distBetweenPoints(object.points[p], axisOfRotation);
        object.momentOfInertia += object.points[p].mass * r*r;
    }
}
    

From there we can simulate undisturbed motion as a combination of rotation and translation, with rotation occuring around the center of mass.

Why is the center of mass the axis of rotation? Take a look:

Each point, in order to continue in a circular path around its axis, must experience a force toward that axis with a magnitude \(F=\frac{mv^2}{r}\). If the object spins with with angular velocity \(\omega\), then \(F=\frac{m (\omega r)^2}{r}=mr \omega^2\). In vector form: \(\vec{F}=m \omega^2\, \vec{r}\) (where \(\vec{r}\) is a vector from the point to the axis of rotation). Because the object does not experience any outside force, the net force must be zero: \(\sum \vec{F} = (0,0)\). Substituting: \(\omega^2 \sum m_p\vec{r}=(0,0)\) and \(\sum m_p(\vec{c}-\vec{p})=(0,0)\) (where \(\vec{p}\) is the point location and \(\vec{c}\) is the axis of rotation). Continuing: \(\vec{c}\sum m_p-\sum m_p\vec{p}\) and finally \(\vec{c}=\frac{1}{m_{total}}\sum m_p\vec{p}\), which is the definition of the center of mass.

Cool! So now we have some rotating and moving objects:

Forces

Forces were one of the easier implementations, outside the collision ordeal, of course. Outside forces affect motion as per Newton's 2nd law: \[F=ma\] It is easy to merely say this, but in practice, in the discrete world of computers, it is a little more difficult. Even with window.requestAnimationFrame(), the time between frames is not garunteed.

Hence the use of the timeStep variable. Given an interval of time \(t\), the change in velocity can be approximated by \(\frac{d\,v}{d\,t}t\) or \(a t\). Likewise, as I neglected to mention in the previous section, \(\Delta p=\frac{d\,p}{dt}t=vt\) and \(\Delta \theta=\frac{d\,\omega}{d\,t}t=\omega t\).

var gravitationalField = {x: 0, y: 1000};
function computeForces(object, dt){
    var netForce = {x:0,y:0};
    netForce = add(netForce,scalar(object.mass, gravitationalField));

    var netAcceleration = scalar(1/object.mass, netForce);
    var deltaV = scalar(dt, netAcceleration);
    object.velocity = add(object.velocity, deltaV);
}
    

The only example currently implemented is gravity, though electric and magnetic influence would be cool next step. It does allow arbitrary directions of gravity though:

Collisions

Collisions were a little difficult. So difficult that I am not ashamed to say that I consulted Chris Hecker's awesome article for guidance.

But there were many other challenges besides collision response. How do you identify a collision? How do you prevent objects from passing through each other? How do you get a pile of squares to behave themselves, each with the square above pushing down on it and the square below pushing up?

My process for detecting collisions went like this:

1. For every pair of objects, determine if any point lies inside the other object. I used the magical algorithm found here.
2. If any points intersected, find the two line segments which intersect between the two objects.
3. Pick the line intersection pair farthest away from the intersecting point. This ensures that, should a point pass all the way through an object, the correct collision is chosen.
4. Project the intersection point onto the other object's line segment. Use this to find the amount that the object has been displaced.
5. Feed this information into Chris Hecker's formula.

The orange lines are all intersecting lines, the red point the intersecting point and the green point the collision point (assuming that the objects had not already passed through each other).

I then translated the offending object so that it was not intersecting the other polygon. Although this has the downside of creating a jerky look, it prevents objects that are smashed in a corner from passing through a wall.

In a continued effort to avoid objects passing through other objects, I tried something new. When an object has been involved in a collision in the current frame, the translation vector is projected onto a vector perpendicular to the vector from the center of mass to the collision:

//moves each object according to its linear/angular velocity
function moveObjects(timeStep){
    for(var p = 0; p < objects.length; p++){
        var object = objects[p];
        //delta x = v * t
        //experimental:
        //don't translate objects in the direction of a current collision
        var translate = {x: timeStep * object.velocity.x,y: timeStep * object.velocity.y};
        for(var k = 0; k < object.collisions.length; k++){
            var collVector = sub(object.collisions[k].collisionLoc, object.centerOfMass);
            //if translate has a component in the direction of collVector
            if(distBetweenPoints(collVector, translate) < magnitude(collVector)){
                //project onto a vector perpendicular to collVector
                translate = project(translate, rotateVec(collVector, Math.PI/2));
            }
        }
        translateObject(object, translate.x, translate.y);
        //delta theta = w * t
        rotateObject(object, object.centerOfMass, timeStep * object.angularVelocity);
    }
}
    

So there you have it. My own, homegrown rigid body simulation. Here's the completed (for now) project:

It definately has some major glitches. I think some predictive collision algorithms may be in order, but that it for another time.