Adventures with 2D Mass-Spring Networks – Part I

Being faithful to Floating Sandbox’s roadmap, after releasing version 1.17.5 I’ve focused my attention again to the most important component of Floating Sandbox: the mass-spring network simulation. This is the sub-system that simulates 2D rigid bodies and makes sure that the set of particles that comprise a ship behaves as a solid whole, rather than as an incoherent soup of particles.

The mass-spring network simulation is responsible for more than 60% of the CPU time spent by Floating Sandbox when simulating a medium-sized ship, and it’s the primary cause of lagging when simulating larger ships. An improved implementation would allow Floating Sandbox to have a higher FPS rate with larger ships and, as we will see, would also improve the rigidity of structures, resulting in firmer ships that bend less for example when colliding with the ocean floor.

In this series of posts I’ll guide you through the experiments that I’ve conducted to better understand and improve the simulation, sharing what I’ve learnt and with the hope that my results will come handy to anyone else playing with the fantastic world of physics simulations. You should only need some basic high-school-grade physics and some programming knowledge in order to follow the series.

2D MASS-spring networks: the basics

A mass-spring network is a special case of a particle system – a set of particles, each being infinitely small (i.e. without a width nor a height) and each having its own mass. In a mass-spring network, these particles are connected to other particles by means of springs, which physically speaking behave exactly like springs in the real world: they basically spend their whole existence trying hard to maintain their original length, which we call rest length. The amount of “effort” with which they try to extend (when compressed) or reduce (when stretched) is the spring’s stiffness, which is measured as a simple (scalar) number.

In the 2D world of Floating Sandbox, the particles are laid out at the corners of a regular 2D grid, and are connected – via springs – to all of the particles that surround them:

By looking at this example, you should easily convince yourself of the following basic relationships, which will come handy when we delve into implementation details:

• Each particle may be connected to up to 8 other particles – in other words, each particle may be the endpoint of up to 8 springs.
  • To be frank, in Floating Sandbox there are also ropes, which may be connected to arbitrary particles, and thus in the game each particle may be connected to up to 9 other particles. In the rest of this series we’ll forego that however, so forget about it.
• The number of springs is roughly equal to 4 times the number of particles. You can see this by imagining adding particles to the example above, and realizing that most of the times, a new particle would generate 4 springs.

---

Physically-speaking, when a spring is extended or compressed, it exerts a force f on its two connected masses in order to move them back to the positions at which the length of the spring equals again its rest length. This force it exerts is linearly proportional to the amount of compression or extension, and as a matter of fact its formula (Hooke’s law) is really simple and states that:

Where d is the current length of the spring, and r is the spring’s rest length. Basically, the force is equal to the amount by which the spring is shorter or longer than its rest length, multiplied by the number that represents the spring stiffness, k. You can easily verify that a higher stiffness generates a stronger force, and that a spring at rest generates a zero force.

More in detail now, the network consists of N_p particles, each particle P_i having mass m_i, and connected to a number of other particles via springs S_j, each having stiffness k_j.

simulatiON

Simulating a mass-spring network consists of “evolving” the system of particles from a point in time to the next. Specifically, in Floating Sandbox the simulation operates as follows:

• Given the positions, velocities, and forces on all particles at a given moment in time (at this frame),
• Calculate the new positions (and velocities) of all particles at a subsequent moment in time (at the next frame),
• Taking into account the rigidity constraints exerted by the springs.

Now, the core of the question is how the simulation evolves the system from a moment to the next. To this end, there’s a number of algorithms in literature – ranging from simple physics-based simulators to more complex algorithms based on solving discrete optimization problems – each algorithm having its own pro’s and con’s and working better under different sets of circumstances and assumptions. In this series we’ll examine a few of these algorithms with an eye on their applicability to Floating Sandbox.

Before we embark on comparing an algorithm to the other, however, we need to decide how we intend to do so. Intuitively, the best algorithm for Floating Sandbox is an algorithm that satisfies the following two constraints:

1. The algorithm is fast enough to allow for a decent FPS rate even for large ships;
2. The algorithm provides a decent impression of rigidity – in other words, ships look firm and not jelly.

Ideally what we need now is an environment where we can benchmark algorithms against each other based on these two “goodness” measures. Enter SpringLab.

SpringLab

SpringLab is a C++ application (open-source, on GitHub) that I’ve written in order to evaluate mass-spring network simulation algorithms. It provides an environment for testing, fine-tuning, and measuring the effectiveness of mass-spring network simulation algorithms, allowing me to experiment with implementations and to compare implementations with each other.

Within SpringLab you create a mass-spring network by loading an image (created e.g. via Microsoft Paint) whose pixels have colors that correspond to a small set of “materials” (if you’re familiar with Floating Sandbox, you’ll find this technique somewhat familiar). Once a mass-spring network is loaded, you can tweak the environment and run any of the simulation algorithms implemented within the application, taking measurements that tell how well an algorithm fits our needs.

Now that we have the means of evaluating the goodness of a simulation algorithm, it’s time to get our hands dirty and start delving into the intricacies of writing mass-spring network simulations. We are going to begin with the simplest algorithm one may think of.

Algorithm 1: Classic Simulator

One simple way of going about simulating a mass-spring network is to use the simple laws of Newtonian mechanics to predict where a particle shall be in the next frame given its current state.

Sir Isaac Newton

Newton’s laws state that:

• A particle’s position changes from a moment to the next because of the particle’s velocity;
• A particle’s velocity changes from a moment to the next because of the particle’s acceleration, which in turn is imparted to the particle by means of a force.

More specifically:

• After a period of time , a particle with a given velocity has moved by
• A force applied to a particle of mass causes a change in velocity (acceleration) equal to
• After a period of time , an acceleration has changed the particle’s velocity by

Assuming that in a given time interval both the particle’s velocity and the force applied to it are constant (and this is a big assumption, as we’ll see later), these laws can be applied to evolve the system from time t to time as follows:

Great, so now we know how to calculate the next state of the system given current positions, velocities, and forces. Our next task is to calculate the forces that our particles are subject to during the simulation’s frame. In our mass-spring network, these forces consist of two components:

• External forces: these are “environmental” forces that we want to apply to the particles in our simulation. One trivial such force is gravity, but in our simulation we can go as crazy as we want and add, for example, buoyancy, blast forces from explosions, and forces generating from collisions with immovable objects (e.g. the ocean floor in Floating Sandbox).
• Spring forces: these are the forces exerted by each spring in our network to the two particles at its endpoints, which we have previously seen as being .

We now have all the building blocks that we need in order to simulate a mass-spring network using basic Newtonian mechanics. This algorithm is implemented in SpringLab as the Classic Simulator – a very naive implementation of the formulas given above.

Let’s see how the algorithm behaves when we move one endpoint of a simple spring after having pinned the other endpoint:

Nice, the spring seems to behave realistically enough. What happens with a slightly more complex structure?

Uhm, indeed, you can tell the square of particles is trying hard to maintain its structure, moving particles up and down in an effort to bring all springs’ lengths back to their rest lengths, but the structure is too “jelly”. You can see this also in the previous animation, the spring quickly recoils to its original length, but the particle oscillates around its rest position for a bit too long.

What we see happening in the first animation is the result of the particle gaining velocity while the extended spring recovers its rest length (remember, the spring force increases the velocity of the particle towards the center of the spring), and this velocity causes the particle to overshoot the rest position and land on the “other side” of it. Once on the other side, the spring is now compressed, and thus a new spring force causes the particle to go back towards the center of the spring, causing it to overshoot again – albeit by a bit less – and the whole cycle begins again. Luckily enough, the “overshoot” velocities generated by each of these oscillations get smaller and smaller, and thus the oscillations die out quite quickly.

One possible way of reducing the oscillations is to transform our spring-mass network into a spring-mass-damper network. Such a network is a mass-spring network with added “shock absorbers” which damp the velocities that are gained by particles when they are subject to spring forces. Conceptually, our basic device in the mass-spring-damper network would look like this:

The force exerted by a damper onto a moving particle is linearly proportional to the velocity component of that movement which is parallel to the damper, via a (scalar) coefficient which we’ll call .

Without adding many more details, after introducing damping the two simulations above look as follows:

Slightly better, huh? The oscillations are still visible, but they have a shorter duration now and we can eventually damp them further by increasing the coefficient.

The complete Classic Simulator algorithm looks as follows:

In the algorithm:

• The procedure takes as input the positions and velocities of all particles, together with all of the springs
• Lines 11, 12: the Hooke force is calculated by measuring the displacement between the two particles, subtracting from it the spring’s rest length , and multiplying this difference in length (compression or extension) with the spring’s stiffness ; this force is directed along the spring direction – i.e. along the unit vector
• Lines 14, 15: the damper force is calculated next by projecting the relative velocity of the two particles onto the spring direction (by means of a dot-product), scaled by the spring’s damping coefficient
• Lines 17, 18: the forces just calculated are added to the total spring forces of each of the two particles

How does the Classic Simulator fare with respect to rigidity? Does it simulate rigid bodies convincingly enough? In order to answer these questions, SpringLab implements a bending measurements that we can use to compare rigidity across different simulation algorithms. The bending measurement basically evaluates the displacement of the following ideal model of a bar of particles, subject to gravity and with one side anchored to a hypothetical “wall”:

Ideal bending model

Let’s see how our Classic Simulator fares with bending:

Ugh – it’s not even bending like a bar, it bends like a stripe of soft rubber. Fine, we know what’s going on here: the Hooke forces are too weak for the structure to maintain its shape. We know that we can control the strength of Hooke forces by increasing the springs’ stiffness. Let’s do that and see what happens when we increase the Classic Simulator’s spring stiffness parameter to 139,000 (slow motion):

What has just happened there? The network has exploded. Welcome to problem #1 with mass-spring network simulations: integration.

problem #1: Euler’s Integration and DISCRETE time

Newton’s laws are valid in our continuous world, in which time moves from an instant to the next by infinitesimal time intervals. Computer simulations and mathematical models, however, deal with finite time intervals; the system being simulated evolves from a point in time to another point in time . We can choose our to be as short as we want – at the cost of slowing down the system evolution in the process, like a movie in slow-motion – but this time interval will always be finite nonetheless.

As we’ve noted earlier, when we calculate the position and velocity of the particle at the next time point , we assume that velocities and forces are constant during that time interval. This assumption has a nasty side-effect, because it causes a particle swept up in Hooke’s forces during its spring’s relaxation to overshoot its target position – the position at which the spring is at its rest length. The following diagram should give you an idea of how this overshoot happens:

In this diagram we see that after the spring is extended by 0.5 in frame B, the resulting Hooke force, assumed to be constant during the time interval, causes the particle to overshoot its rest position to the “other side” (frame C), at which point the spring is compressed by 0.6 – more compressed than the amount by which the spring was extended in the previous frame. This compression causes now a new Hooke force – opposite to the previous one, and larger in magnitude – which in turn pushes the particle to the opposite side during the next interval, after which the spring is extended by 0.7 (frame D). You can tell that this is a vicious circle, and after enough cycles, the particles will be oscillating so wide that the system drowns in floating point overflows.

In the real world, the spring force applied to a particle that is attached to a recoiling spring would not be constant, but it would smoothly weaken instead as the particle travels towards its target position during infinitesimal time steps and the spring’s extension becomes smaller and smaller. The force would eventually decay to zero as the particle reaches its destination.

There are at least two ways in which we can avoid this pitfall.

For starters, this method of predicting the particle’s position and velocity in the next time step by assuming constant forces during the time interval – called Explicit Euler Integration – is not the only method available to us. There are other methods, among which Implicit Euler Integration, that do not have this downside, but which are in turn difficult to implement and time-consuming to run with every frame of the simulation.

Another cheaper to implement and more intuitive way of counteracting this issue is to choose for a smaller . Referring to the diagram above, say that shortly after frame B we were to re-sample the spring extension and recalculate a new Hooke force. This force would certainly be smaller, given that the spring extension would now be shorter than 0.5. This new and weaker Hooke force would then be less likely to overshoot the rest position, and the explosion we have witnessed earlier would not happen. Indeed, you can try this out yourself by changing the time step interval to 0.006198 and the spring stiffness to 229376.0; the bar becomes stiffer, at the expense of the simulation becoming annoyingly slow.

The problem is that the value that we choose for our time step is closely interconnected to how real-time we want our simulation to be. Say that we’re targeting 64 frames per second (like in Floating Sandbox); if we wanted a 1:1 correspondence between the simulation and the real world (i.e. one second in the simulation corresponding to one second in the real world), we’d have to choose 1s/64=15.625ms as our . Anything smaller than that, when played at 64 frames per second, would result in a slower-than-real simulation.

An alternative would be to allow for more than 64 frames per second; for example, we could choose and play 1000 frames per second. Surely that would be fantastic! But you probably already know where the gotcha is – in order to play 1000 frames per second, the CPU time required to simulate one frame can’t be more than 1ms. If your implementation of the mass-spring network simulation takes more than 1ms per frame, you’ll produce less than 1000 frames per second, causing the simulation to “lag” and look like a movie in slow-motion.

blankets too short

In conclusion, it appears that we are dealing with a blanket that is too short: to counteract Euler’s integration weaknesses and have stable simulations with a reasonably rigidity, we need a small ; however, the smaller the , the more the number of frames we need to be able to simulate in one second, and there’s a limit to how many frames we can simulate in one second, which is dictated by our CPU and by the number of calculations we need to do in a frame (which is also related to how large the mass-spring network is).

In the following posts we’ll see how I’ve tackled these issues in Floating Sandbox, and then we’ll look at a couple of alternative mass-spring network simulation algorithms which promise more stiffness for less CPU time. Stay tuned!