Tag Archives: Physics

Implementing A Mass Spring System

Here’s an older post I wrote a few years ago when I was dabbling in OpenGL programming. The code is written in C++ and leaves a lot to be desired, but I figured I would post it in case anyone was wondering how to implement something like this.


This is a quick tutorial on how to implement a mass spring system in your game. You can use this to simulate objects like hair, string, cloth, and things of that nature. We will discuss the abstract implementation of a spring system, and then an example of this system simulating a string will be given (mainly because it is the most simple to create and imagine), but this can be extended to any other type of object fairly simply. Only the simulation will be defined for now, no collision detection will be discussed.

Internal spring forces:
First lets start off with how the springs will be defined. A spring has two nodes, a left and a right. These two nodes have an impact on the force created when a node is moved around relative to the other. Lets call them nodes “a” and “b” of the spring “A”. Imagine grabbing node “a” with your left hand and node “b” with your right hand and pulling them apart. The very nature of a spring says that it will experience a restoring force, and naturally want to pull back to its inertial state. If we calculate the force with respect to “a” we see that it pulls in the direction of node “b” and node “b” pulls in the opposite direction of that force (toward node “a”). This type of system is called a harmonic oscillator, and the force can be described by Hooke’s Law: F_{spring} = -k*\Delta{x}, where k is the elasticity coefficient (a.k.a the spring constant) and \Delta{x} is the change of the springs length from the inertial state. So when \Delta{x} is positive (i.e. the spring is being pulled apart) the force is negative (representing contraction), and the inverse is true when \Delta{x} is negative (the spring is being pushed together).

Mathematics of the net force:
To describe the NET force we use \vec{F}_{net} = \frac{d^2}{dt^2}(m\vec{u}) = \frac{d}{dt}(\vec{p}) = \frac{d}{dt}(m\vec{v}) = m\vec{a} (where \vec{u} is the spring nodes position, \vec{v} is its velocity, and \vec{p} = m\vec{v} = momentum). This formula should hopefully look familiar, its Newton’s famous second law of motion! Note that we assume mass is constant regardless of velocity, so we can factor it out when we differentiate. As you can see the formula states that the force on an object is directionally proportional to its mass and acceleration, and with some algebraic manipulation we can solve for acceleration and get \vec{a}=\vec{F}/m meaning that acceleration is directly proportional to force and indirectly proportional to mass.

External forces:
There are many forces that can act on a spring. Gravity is the most common, followed by air resistance. Gravity is a force vector causing an acceleration of approximately \begin{pmatrix}0 \\ 9.81 \\ 0 \end{pmatrix} m/s^2 (according to Newton’s acceleration due to gravity constant for earth little, dubbed “g” or little g). Wind resistance can be a vector of any direction. We also need to note that energy is conserved in a physical system, so we need a friction factor to dampen the force, or else the system would continue to operate forever. We calculate this factor by saying \vec{F}_{friction} = -k_{friction} * \vec{v}_{relative}, where k_{friction} is the friction coefficient and \vec{v}_{relative} = \vec{v}_{b} - \vec{v}_{a}. We apply this force in the opposite direction of the current velocity, hence the “-” in front of kFriction.

Anchor nodes:
The concept of an anchor node is simple. It is a node that does not let any force act on it. In the case of a string think of it as being a part of the string nailed to a wall. This node will not budge from its spot.

Putting it all together:
Heres what the algorithm for one iteration to animate a spring will look like in a nutshell:
– For each moving spring node -> ComputeExternalForces acting on it and add it to the node’s force vector
– For each moving spring node -> Compute the internal force neighbor nodes are exerting on this node and add the force to this node’s force vector
– For each moving spring node -> Compute the nodes acceleration (a=f/m) and add it to the spring’s velocity, and add the velocity to the spring’s position in space
– For each moving spring node -> Compute the energy loss (friction)


Here are some screen shots that show the end result:

An anchor node in the middle:

An anchor node on the left end:


And here’s the code:

struct SpringNode {

CVector3<double> position; // My current position
CVector3<double> velocity; // My velocity

CVector3<double> force;

double mass; //
std::vector<SpringNode> neighbors; // Neighbors
std::vector<double> neighborDistances; // How far am I away from my neighbors
};

struct AnchorNode {
SpringNode *node;
CVector3<double> position; // Anchor position
};

class ISpringSystem {

public:

/*!
*/
void ComputeSingleForce( SpringNode const *a, SpringNode const *b );

/*!
*/
void AddAnchorNode( int node, CVector3<double> const &anchorPosition );

/*!
*/
void Animate();

/*!
*/
void Draw();

/* The derived class will connect the springs for us.
*/
virtual void ConnectSprings( )=0;

protected:
double _springK; // Elasticity coefficient. The larger K is, the stiffer the spring.
double _inertialDistance; // Inertial or resting distance of the springs
CVector3<double> _force; // Temporary force vector for ComputeSingleForce
CVector3<double> _vel; // Temporary velocity vector for ComputeSingleForce
std::vector<SpringNode> _springNodes; // Vector of spring nodes
std::vector<AnchorNode> _anchorNodes; // Vector of anchor nodes
std::vector<SpringNode*> _movingNodes; // Vector of moving nodes
double _mass; // Mass of each spring node
double _energyLoss; // Energy loss to occur as the simulation progresses
CTimer *_timer; // Timer used in physics calculations

};

#define SPRING_TOLERANCE (0.0000000005)

void ISpringSystem::ComputeSingleForce( SpringNode const *a, SpringNode const *b )
{
// Calculate vector between b and a (giving the direction and magnitude of the vector they make up)
_force = b->position - a->position;

// Calculate the relative velocity between a and b
_vel = b->velocity - a->velocity;

//
if ( _force.Length() - _inertialDistance < SPRING_TOLERANCE )
{
_force.Zero();
return;
}

// Calculate the intensity of the spring ( -k * deltaDistance )
// This is where the spring constant comes into play. A higher K will yield
// a stiffer spring.
double intensity = -_springK * (_force.Length() - _inertialDistance);

// Normalize the force direction
_force.Normalize();

// Apply the internal force and the friction force
//
_force += ((_force/_force.Length())*intensity) - (40000 * _vel);
}

void ISpringSystem::AddAnchorNode( int node, CVector3<double> const &anchorPosition )
{
// Setup the anchor nodes
_anchorNodes.resize(_anchorNodes.size()+1);
_anchorNodes[_anchorNodes.size()-1].node = &_springNodes[node];
_anchorNodes[_anchorNodes.size()-1].position = anchorPosition;
_movingNodes.erase( _movingNodes.begin() + node );
}

void ISpringSystem::Animate()
{
static bIsInit=false;

// Initialize the timer
if ( !bIsInit ) _timer->Init();

// Get the elapsed seconds
double elapsed = _timer->GetElapsedSeconds();

// Gravitational force
// To compute the gravitational force we use F = ma
// (m = mass, a = acceleration due to gravity (in units per second which is g=9.81m/s^2)
double fg = -9.81 * _mass;

// Compute all external forces on the springs
for( int i=0; i<_movingNodes.size(); ++i )
{
// Zero out the force vector
_movingNodes[i]->force.Zero();

// Factor in the gravitational force. Since it moves only in the y direction, we
// only modify the y component of the vector.
_movingNodes[i]->force[1] = fg;
}

// Compute internal forces on the springs
for( int i=0; i<_movingNodes.size(); ++i )
{
// Compute the force the neighbor nodes are exerting on each node
for( int j=0; j<_movingNodes[i]->neighbors.size(); ++j )
{
//
ComputeSingleForce( _movingNodes[i]->neighbors[j], _movingNodes[i] );
_movingNodes[i]->force += _force; // Add the force exerted to this nodes force vector
}
}

// Increment velocity and position of nodes
for( int i=0; i<_movingNodes.size(); ++i )
{
_movingNodes[i]->velocity += (_movingNodes[i]->force/_mass)*elapsed; // Add
_movingNodes[i]->position += _movingNodes[i]->velocity; // Add the velocity to pos
_movingNodes[i]->velocity *= _energyLoss;
}

}

void ISpringSystem::Draw()
{
// Push the current attribute list onto the OpenGL attribute stack
glPushAttrib(GL_CURRENT_BIT);
glLoadIdentity();

// Draw as lines
glBegin(GL_LINES);
glDisable(GL_CULL_FACE);

glColor3f(1.0f, 0.0f, 0.0f);

for( int i=0; i<_springNodes.size()-1; ++i )
{
glVertex3f( _springNodes[i].position[0], _springNodes[i].position[1], _springNodes[i].position[2] );
glVertex3f( _springNodes[i+1].position[0], _springNodes[i+1].position[1], _springNodes[i+1].position[2] );
}

glEnd();

// Pop the last element off OpenGL attribute stack
glPopAttrib();
}

// String spring system class
class CStringSpringSystem : public ISpringSystem {

public:
void ConnectSprings();

};

void CStringSpringSystem::ConnectSprings( )
{
_springNodes.resize( 80 );
_movingNodes.resize( 80 );

_timer = new CHiResTimer;
_timer->Init();

_springK = 8000; // Rigidity of springs
_inertialDistance = .05; // Springs are 5 cm apart
_mass = .05; // Each mass is .05 kg
_energyLoss = 1-.00001; // Amount of energy the system will lose as time elapses

// Set springNode distances along the positive x-axis
for ( int i=0; i<_springNodes.size(); ++i )
{
_springNodes[i].position = CVector3<double>( 0.0f + i*_inertialDistance, 1.0f, -5.0f );
_movingNodes[i] = &_springNodes[i];
_movingNodes[i]->velocity.Zero();
}

// Set up the first spring node
_springNodes[0].neighbors.resize(1);
_springNodes[0].neighborDistances.resize(1);
_springNodes[0].neighbors[0] = &_springNodes[1];
_springNodes[0].neighborDistances[0] = _inertialDistance;

// Set up the last spring node
_springNodes[_springNodes.size()-1].neighbors.resize(1);
_springNodes[_springNodes.size()-1].neighborDistances.resize(1);
_springNodes[_springNodes.size()-1].neighbors[0] = &_springNodes[_springNodes.size()-2];
_springNodes[_springNodes.size()-1].neighborDistances[0] = _inertialDistance;

// Set up all nodes in between
for( int i=1; i<_springNodes.size()-1; ++i )
{
_springNodes[i].neighbors.resize(2);
_springNodes[i].neighborDistances.resize(2);

//
_springNodes[i].neighbors[0] = &_springNodes[i-1];
_springNodes[i].neighborDistances[0] = _inertialDistance;

//
_springNodes[i].neighbors[1] = &_springNodes[i+1];
_springNodes[i].neighborDistances[1] = _inertialDistance;
}

// Make the middle node an anchor node (i.e. it does not move)
AddAnchorNode( (_springNodes.size()-1)/2, CVector3<double>( 0.0f, _springNodes[(_springNodes.size()-1)/2].position[1], -5.0f ) );
}