Verlet collision with impulse preservation

Nov. 29, 2010
comments

This article presents a system to integrate verlet physics collision with preservation of impulse. Normally verlet physics annihilates a lot of energy from a system, which makes it very stable but also quite unrealistic. Additionally simple methods of preserving impulse yield very unstable systems, a limitation which can be overcome by two steps of integration, one for at-rest acceleration canceling, and one with impulse preservation.

I'm building on the concept of verlet integration that I explained in a previous article. You might also be interested in different integration methods.

Contents

The solution sneak peak

The simulation below is the result of the tutorial to follow. Hover over it with the mouse to run the simulation. You can click into the canvas to create new spheres of random size.

You can put this on your own page if you want by pasting in the code blow.

<iframe
src         = "http://codeflow.org/frames/verlet.html#color=#f99e00,bg=none,start=20"
width       = "300"
height      = "200"
frameborder = "0"
scrolling   = "no"
/>

The arguments after the hash control the circle color, background color and start count of bodies, you can set any size you wish.

The Problem

Verlet physics is very good at resolving hard constraints in a stable fashion. However, in the process it destroys a lot of energy. The issue is visible with the code below.

Setup

The simulation consists of spherical bodies each of which is collided with every other body and the level boundary.

Body

Each body has a radius, position, previous position and accumulated acceleration.

var Body = function(x, y, radius){
    this.x = x;
    this.y = y;
    this.px = x;
    this.py = y;
    this.ax = 0;
    this.ay = 0;
    this.radius = radius;
}

Acceleration on a body is accumulated in the acceleration components and this method transfers the accumulated acceleration to a new current position.

accelerate: function(delta){
    this.x += this.ax * delta * delta;
    this.y += this.ay * delta * delta;
    this.ax = 0;
    this.ay = 0;
},

Transfer of existing inertia to a new position in a verlet integration step is done by the following method:

inertia: function(delta){
    var x = this.x*2 - this.px;
    var y = this.y*2 - this.py;
    this.px = this.x;
    this.py = this.y;
    this.x = x;
    this.y = y;
},

A body is drawn with the following code:

draw: function(ctx){
    ctx.beginPath();
    ctx.arc(this.x, this.y, this.radius, 0, Math.PI*2, false);
    ctx.fill();
},

Simulation

The simulation maintains a list of all bodies and runs this method for each simulation step:

var step = function(){
    var steps = 2;
    var delta = 1/steps;
    for(var i=0; i<steps; i++){
        // simulation code here
    }
    draw();
}

It has a few helper methods to apply the discrete steps of the simulation to the bodies (like acceleration, gravity, inertia and so on).

Pure verlet simulation steps

This is the content of the inner loop of each step for the energy annihilating pure verlet solution.

gravity();
accelerate(delta);
collide();
border_collide();
inertia(delta);

Gravity accumulates the gravitational acceleration on bodies, accelerate moves the bodies according to acceleration. Note that this results in a higher velocity (as it should) in the verlet integration scheme, since the current position moved away further from the previous position.

Collide resolves body vs. body constraints (two spheres may not overlap):

var collide = function(){
    for(var i=0, l=bodies.length; i<l; i++){
        var body1 = bodies[i];
        for(var j=i+1; j<l; j++){
            var body2 = bodies[j];
            var x = body1.x - body2.x;
            var y = body1.y - body2.y;
            var slength = x*x+y*y;
            var length = Math.sqrt(slength);
            var target = body1.radius + body2.radius;

            // if the spheres are closer
            // then their radii combined
            if(length < target){ 
                var factor = (length-target)/length;
                // move the spheres away from each other
                // by half the conflicting length
                body1.x -= x*factor*0.5;
                body1.y -= y*factor*0.5;
                body2.x += x*factor*0.5;
                body2.y += y*factor*0.5;
            }
        }
    }

Inertia processes a bodies velocity derived from its previous position.

It is easy to see why this relatively simple scheme does not preserve impulse energy. The verlet constraint resolution just moves the bodies around in order to be conflict free, and the velocity is encoded in the previous position, which is not adjusted for impulse preservation.

Preserving impulse

We could move the previous position around at each collision according to the laws of deflection and impulse transfer as to reflect an energetically more correct system state.

Simple model for deflection impulse

The impulse that two spherical bodies can exchange when colliding acts co-linearly to the normal on their contact surface, leading trough their centers of mass. Friction and rotational impulse are ignored for this tutorial.

The velocities of the bodies have to be decomposed into two components that represent the size of the velocity in the direction of the possible impulse force. This can be done by vector projection. Here the velocity is the red vectors and the impulse direction component is green.

Now the impulses need to be swapped between the bodies, which means subtracting each from one body and adding it to the other which results in a new velocity (blue):

The code in the collision loop to do this is below.

if(length < target){
    // record previous velocity
    var v1x = body1.x - body1.px;
    var v1y = body1.y - body1.py;
    var v2x = body2.x - body2.px;
    var v2y = body2.y - body2.py;

    // resolve the body overlap conflict
    var factor = (length-target)/length;
    body1.x -= x*factor*0.5;
    body1.y -= y*factor*0.5;
    body2.x += x*factor*0.5;
    body2.y += y*factor*0.5;

    // compute the projected component factors
    var f1 = (damping*(x*v1x+y*v1y))/slength;
    var f2 = (damping*(x*v2x+y*v2y))/slength;

    // swap the projected components
    v1x += f2*x-f1*x;
    v2x += f1*x-f2*x;
    v1y += f2*y-f1*y;
    v2y += f1*y-f2*y;

    // the previous position is adjusted
    // to represent the new velocity
    body1.px = body1.x - v1x;
    body1.py = body1.y - v1y;
    body2.px = body2.x - v2x;
    body2.py = body2.y - v2y;
}

This is better in the preservation of impulse, however it is not stable. Ghost impulse is being created by the application of acceleration and the subsequent preservation. See the example below

If you look closely you see a jittering effect. It would be possible to minimize the jittering by increasing the step count, but this would be very expensive method, and the jittering would remain (just at smaller scales).

Acceleration at rest

If you think about the problem it's obvious why the ghost impulse is introduced. If a sphere lies at rest in equilibrium on top of another sphere or the floor, then the acceleration creates velocity where none would have existed in the first place. The characteristic of acceleration at rest is that the pushback from the constraint exactly cancels the acceleration out.

The idea of correcting the simulation for this is to go it in two steps.

  1. Accumulate forces and move object, resolve constraints but without preserving impulse. This gives an object the chance to cancel acceleration impulse out when in equilibrium.
  2. The resulting velocity remaining after step 1 must have been unconstrained, therefore it is now real impulse, so we can move the object by its inertia and preserve its impulse when constrained.

For this the core step loop is adjusted:

gravity();
accelerate(delta);
collide(false);
border_collide();
inertia(delta);
collide(true);
border_collide_preserve_impulse();

The collide method now takes an argument (true or false) which indicates if impulse preservation should be performed (the boundary collision is also split).

if(length < target){
    var v1x = body1.x - body1.px;
    var v1y = body1.y - body1.py;
    var v2x = body2.x - body2.px;
    var v2y = body2.y - body2.py;

    var factor = (length-target)/length;
    body1.x -= x*factor*0.5;
    body1.y -= y*factor*0.5;
    body2.x += x*factor*0.5;
    body2.y += y*factor*0.5;

    if(preserve_impulse){
        var f1 = (damping*(x*v1x+y*v1y))/slength;
        var f2 = (damping*(x*v2x+y*v2y))/slength;

        v1x += f2*x-f1*x;
        v2x += f1*x-f2*x;
        v1y += f2*y-f1*y;
        v2y += f1*y-f2*y;

        body1.px = body1.x - v1x;
        body1.py = body1.y - v1y;
        body2.px = body2.x - v2x;
        body2.py = body2.y - v2y;
    }
}

You can see that the obnoxious jittering is gone. To achieve similarly stable results without the at-rest optimization, 20 or more internal steps would have to be performed.

Wrap up

I've demonstrated a method to implement very hard collision contacts and an optimization that makes it possible to have very nice at-rest state while also preserving impulse. I hope you enjoyed the tutorial.