Hello, I’m Glenn Fiedler and welcome to Virtual Go, my project to simulate a Go board and stones.

In previous articles we mathematically defined the go stone, rendered it, determined how it moves and rotates, accelerated it by gravity and detected when it collides with the go board. We also discussed the inertia tensor and how it affects the movement of the go stone.

In this article we reach our first milestone: a go stone bouncing and coming to rest on the go board.

We do this using a technique called impulse-based collision response.

The concept is simple: apply an impulse at the point of collision to make the go stone bounce.

Exactly how to calculate and apply this impulse is rather complicated and interesting. So read on!

Linear Collision Response

We now pick up where we left off at the end of the collision detection article.

We have a contact point and a contact normal for the collision.

Lets start by calculating a collision response impulse without rotation.

First, take the dot product of the linear momentum of the go stone with the contact normal. If this value is less than zero, continue; otherwise, the stone is moving away from the board and you should not apply a collision impulse.

To calculate the impulse we need the concept of ‘elasticity’.

If the collision is perfectly elastic the go stone bounces off the board without losing any energy:

If the collision is perfectly inelastic, the go stone loses all its vertical motion post-collision and slides along the surface of the board:

What we want is something in between:

To support this we introduce a new concept called the ‘coefficient of restitution’. When this value is 1 the collision is perfectly elastic, when it is 0 the collision is inelastic. At 0.5, it’s halfway between.

This gives the following formula:

Where:

• j is the magnitude of the collision impulse
• e is the coefficient of restitution [0,1]
• p is the linear momentum of the go stone
• n in the contact normal for the collision

Note that the direction of the collision impulse is always along the contact normal so to apply the impulse just multiply the contact normal by j and add it to the linear momentum vector.

Here is the code:

void ApplyLinearCollisionImpulse( StaticContact & contact, float e )
{
    const float mass = contact.rigidBody->mass;
    const float j = max( - ( 1 + e ) * dot( contact.rigidBody->linearMomentum, contact.normal ), 0 );
    contact.rigidBody->linearMomentum += j * contact.normal;
}

And here is the result:

Collision Response With Rotation

Now lets calculate collision response with rotation.

To do this we will now calculate the velocity of the stone at the contact point, then take the dot product of this vs. the contact normal to check if the stone is moving towards the board. This is necessary because the stone is rotating so different points on the stone have different velocities.

Next we use the same basic strategy of applying an impulse along the contact normal with magnitude j except this impulse is now applied to the go stone at the contact point.

It’s actually quite tricky now to calculate the magnitude of this impulse such that the post-collision behavior matches the coefficient of restitution because the impulse at the contact point has both linear and angular effects, and the balance between them depends on the following factors:

• The point of application of the impulse
• The direction of the impulse relative to the center of the stone
• The inertia tensor of the go stone

Here is the general equation for calculating j of the collision response between two moving objects:

You can find an excellent derivation of this result on wikipedia.

Understandably this is quite complex, but in our case, the go board never moves, so we can simplify the equation by assigning zero velocity and infinite mass to the second body. This leads to the following, simpler, but still pretty complex equation:

Where:

• j is the magnitude of the collision impulse
• e is the coefficient of restitution [0,1]
• n in the contact normal for the collision
• v is the the go stone velocity at the contact point
• r is the contact point minus the center of the go stone
• I is the inertia tensor of the go stone
• m is the mass of the go stone

Which turns into the following code:

vec3f GetVelocityAtWorldPoint( const RigidBody & rigidBody, vec3f point ) const
{
    vec3f angularVelocity = transformVector( rigidBody.inverseInertiaTensor, 
                                             rigidBody.angularMomentum );

    return linearVelocity + cross( angularVelocity, point - rigidBody.position );
}

void ApplyCollisionImpulse( StaticContact & contact, float e )
{
    RigidBody & rigidBody = *contact.rigidBody;

    vec3f velocityAtPoint = GetVelocityAtWorldPoint( rigidBody, contact.point );

    const float vn = min( 0, dot( velocityAtPoint, contact.normal ) );

    // calculate inverse inertia tensor in world space

    mat4f rotation;
    rigidBody.orientation.toMatrix( rotation );
    mat4f transposeRotation = transpose( rotation );
    mat4f i = rotation * rigidBody.inverseInertiaTensor * transposeRotation;

    // apply collision impulse

    const vec3f r = contact.point - rigidBody.position;

    const float k = rigidBody.inverseMass + dot( cross( r, contact.normal ), transformVector( i, cross( r, contact.normal ) ) );

    const float j = - ( 1 + e ) * vn / k;

    rigidBody.linearMomentum += j * contact.normal;
    rigidBody.angularMomentum += j * cross( r, contact.normal );
}

And here is the result:

Coulomb Friction

We don’t often get to see frictionless collisions in the real world so the previous result looks a bit strange.

In order to get realistic behavior out of the go stone, we need to add friction.

We’ll model sliding friction between the go stone and the board using the Coulomb friction model.

In this model the friction impulse is proportional the magnitude of the normal impulse ‘j’ and limited by a friction cone defined by the coefficient of friction ‘u’:

Lower friction coefficient values mean less friction, higher values mean more friction. Typical values for the coefficient of friction are in the range [0,1].

Calculation of the Coulomb friction impulse is performed much like the calculation of the normal impulse except this time the impulse is in the tangent direction against the direction of sliding.

Here is the formula for calculating the magnitude of the friction impulse:

Where:

• j_t is the magnitude of the friction impulse (pre-cone limit)
• u is the coefficient of friction [0,1]
• t in the tangent vector in the direction of sliding
• v is the the go stone velocity at the contact point
• r is the contact point minus the center of the go stone
• I is the inertia tensor of the go stone
• m is the mass of the go stone

Which turns into this code:

void ApplyCollisionImpulseWithFriction( StaticContact & contact, float e, float u )
{
    RigidBody & rigidBody = *contact.rigidBody;

    vec3f velocityAtPoint = GetVelocityAtWorldPoint( rigidBody, contact.point );

    const float vn = min( 0, velocityAtPoint, contact.normal );

    // calculate inverse inertia tensor in world space

    mat4f rotation;
    rigidBody.orientation.toMatrix( rotation );
    mat4f transposeRotation = transpose( rotation );
    mat4f i = rotation * rigidBody.inverseInertiaTensor * transposeRotation;

    // apply collision impulse

    const vec3f r = contact.point - rigidBody.position;

    const float k = rigidBody.inverseMass + dot( cross( r, contact.normal ), transformVector( i, cross( r, contact.normal ) ) );

    const float j = - ( 1 + e ) * vn / k;

    rigidBody.linearMomentum += j * contact.normal;
    rigidBody.angularMomentum += j * cross( r, contact.normal );

    // apply friction impulse

    velocityAtPoint = rigidBody.GetVelocityAtWorldPoint( contact.point );

    vec3f tangentVelocity = velocityAtPoint - contact.normal * dot( velocityAtPoint, contact.normal );

    if ( length_squared( tangentVelocity ) > 0.001f * 0.001f )
    {
        vec3f tangent = normalize( tangentVelocity );

        const float vt = dot( velocityAtPoint, tangent );

        const float kt = rigidBody.inverseMass + dot( cross( r, tangent ), transformVector( i, cross( r, tangent ) ) );

        const float jt = clamp( -vt / kt, -u * j, u * j );

        rigidBody.linearMomentum += jt * tangent;
        rigidBody.angularMomentum += jt * cross( r, tangent );
    }
}

And gives the following result:

Rolling Friction

This looks really good but there is one small problem remaining.

Due to its shape, the go stone really prefers to rotate about axes on the xz plane. This means that if you you attempt to spin up around the y axis it stands up and spins on its side like a coin:

This is kinda cool. The problem is that it spins like this forever.

Why? Shouldn’t coulomb friction handle this for us?

No. Coulomb friction only handles friction when the two surfaces are sliding relative to each other.

Sliding friction is just one type of friction and there are many others.

What we have here specifically is combination of rolling and spinning friction.

At this point I have very little patience so I came up with my own hack approximation of spinning and rolling friction that gives me the result that I want: vibrant motion at high energies but slightly damped so the stone slows down, collapses from spinning, wobbles a bit and come to rest.

Here is my horrible code that does this:

    // this is a *massive* hack to approximate rolling/spinning
    // friction and it is completely made up and not accurate at all!

    if ( collided )
    {
        float momentum = length( stone.rigidBody.angularMomentum );
        const float factor_a = DecayFactor( 0.9925f, dt );
        const float factor_b = DecayFactor( 0.9995f, dt );
        const float a = 0.0f;
        const float b = 1.0f;
        if ( momentum >= b || appliedSpin )
        {
            stone.rigidBody.angularMomentum *= factor_b;
        }
        else if ( momentum <= a )
        {
            stone.rigidBody.angularMomentum *= factor_a;
        }
        else
        {
            const float alpha = ( momentum - a ) / ( b - a );
            const float factor = factor_a * ( 1 - alpha ) + factor_b * alpha;
            stone.rigidBody.angularMomentum *= factor;
        }
    }

    // apply damping

    const float factor = DecayFactor( 0.99999f, dt );

    stone.rigidBody.linearMomentum *= factor;
    stone.rigidBody.angularMomentum *= factor;

And here is the end result:

Coming Soon: Digitizing The Go Stone

If you enjoyed this article please donate.

Donations offset hosting costs and encourage me to write more articles!