Rope Simulator in C++
For some time now, I’ve been missing quick and to-the-point tutorials on particle-based physics simulations. Hence, I decided to start sharing things I learned in a concise manner. For starters, we’re going to implement a simple rope simulator. However, before we do anything, a quick disclaimer.
This tutorial is aimed at achieving rapid, visual results. I do not follow perfect software engineering practices and do not explain the math in-depth.
Now that you’ve been warned, let’s get to it!
What is a rope?
What rope is, everybody sees. In order to create a digital representation of one, we divide it into multiple small segments, connected by points. Since rope and spring are two different things, those segments have equal and constant lengths — D.
This representation lets us focus on individual points (or particles) and their properties.
A bit of math
Now that we have the representation of a rope, how can we simulate it? First, we need to update the positions of each point based on the forces applied to the rope, such as gravity. Then, we need to make sure that our rope does not stretch or compress.
Verlet Integration
You’re probably familiar with the simplest of all — Euler integration. At a certain timestep, you update velocity based on acceleration and then update position based on velocity. We need to define a step size, for example, 0.01 seconds. In a more formal way, this method looks like this.
Everything looks great, but over time we lose a lot of precision, and the accuracy of the simulation is highly dependent on the chosen timestep. If you wanna see more details, please see this post. This issue brings us to Verlet integration, which is more stable and accurate. The basic idea is the following, we update our position based on the previous and instant velocity. I won’t derive the whole formula, cause it takes a bit of time to get right. In short, it looks like this.
After we repeat this procedure for every particle, we have successfully calculated the new state of the rope at the next time step. However, there is a problem. We might end up in a situation, where segments have different lengths. This can happen whenever two opposing forces are acting on neighboring particles. Let’s enforce this constraint!
Enforcing constraints
After the Verlet integration, segments may be slightly shorter or longer than before. In the beginning, we defined D to be our target segment length. This means, that if our rope consists of 10 particles, 9 segments, and D=1, then the rope should be 9 units long at all times. To make sure this is the case, we use the Jakobsen method. This method is used whenever a constraint of fixed distance between particles needs to be enforced, so feel free to use it in other projects, such as cloth simulations!
Let’s get to the algorithm. We take two particles (p,q) connected by a segment of length d=||p-q||. If d is smaller than the desired distance D, we push particles apart, each by 0.5(D-d). If they are too far, we do the opposite, we pull them closer, each by 0.5(d-D). Tricky things happen when a certain point is fixed. Think of it, as if the rope was pinned to the wall throughout the simulation. Then, we cannot move this particle around while enforcing constraint, hence we drop the 0.5 factor and move the neighboring particle by |d-D| in the appropriate direction. By repeating this procedure multiple times, our segments should have approximately appropriate lengths.
Code
Finally, we are ready to code it up. I will do it in C++, but feel free to use whatever you like. The translation should be rather easy.
First and foremost, let’s import the necessary libraries and put some global constants in place.
#include <math.h>
#include <stdlib.h>
#include <vector>
float GRAVITY = -9.81; // constant force applied to all particles
float TIMESTEP = 0.01; // Δt in our equations
unsigned JAKOBSEN_ITERATIONS = 50; // Number of times we will enforce the distance constraintFor convenience, let’s define a particle. Remember, we need to keep track of its previous position in order to perform Verlet integration.
struct Particle {
// Current position
float x;
float y;
// Previous position
float xPrevious;
float yPrevious;
// Whether or not particle is able to move
bool fixed;
};We will enclose the whole functionality in a class called Rope. I’ll define a simple constructor that creates a rope between two points. The first point will be fixed, while all others not. Depending on the application, we also want to control the total number of particles in a rope.
class Rope {
public:
Rope(float x1, float y1, float x2, float y2, unsigned nParticles) {
for (unsigned i = 0; i < nParticles; i++) {
// How close are we to the last point?
float w = (float)i / (nParticles - 1);
float x = w * x2 + (1 - w) * x1;
float y = w * y2 + (1 - w) * y1;
Particle p;
p.x = x;
p.y = y;
p.xPrevious = x;
p.yPrevious = y;
p.fixed = i == 0; // We fix only the first point
_particles.push_back(p);
}
// There is one less segment than there are particles
unsigned numberOfSegments = nParticles - 1;
float ropeLength = sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2));
_desiredDistance = ropeLength / numberOfSegments;
}
private:
// We need to store our particles somewhere
std::vector<Particle> _particles;
// Target distance of a single segment
float _desiredDistance;
};To advance the simulation by one step, we define a public function step. As we discussed, a single step consists of a Verlet integration and Jakobsen constraint enforcement.
void step() {
verletIntegration();
enforceConstraints();
}Let’s define a private method verletIntegration, that will take care of the fancy math.
void verletIntegration(){
for (Particle &p : _particles){
if(p.fixed)
continue; // We do not want to move fixed particles
float xCopy = p.x;
float yCopy = p.y;
// Updating particle's position
p.x = 2 * p.x - p.xPrevious + 0 * TIMESTEP * TIMESTEP;
p.y = 2 * p.y - p.yPrevious + GRAVITY* TIMESTEP * TIMESTEP;
p.xPrevious = xCopy;
p.yPrevious = yCopy;
}
}Perfect! Let’s enforce those constraints. We define another private method enforceConstraints.
void enforceConstraints(){
// We perform the enforcement multiple times
for (unsigned iteration = 0; iteration < JAKOBSEN_ITERATIONS; iteration++){
// We iterate over each pair of pa
for (size_t i = 1; i < _particles.size(); i++) {
Particle &p1 = _particles[i-1];
Particle &p2 = _particles[i];
// Calculating distance between the particles
float distance = sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2));
float distanceError = distance - _desiredDistance;
// The direction in which particles should be pulled or pushed
float xDifference = p2.x - p1.x;
float yDifference = p2.y - p1.y;
// We need to make it a unit vector
// This will allow us to easily scale the impact we have
// on each particle's position.
float xDirection = xDifference / sqrt(pow(xDifference, 2) + pow(yDifference, 2));
float yDirection = yDifference / sqrt(pow(xDifference, 2) + pow(yDifference, 2));
// Finally, we can update particles' positions
// We need to remember that fixed particles should stay in place
if (p1.fixed && !p2.fixed) {
p2.x -= xDirection * distanceError;
p2.y -= yDirection * distanceError;
} else if (p2.fixed && !p1.fixed) {
p1.x += xDirection * distanceError;
p1.y += yDirection * distanceError;
} else if (!p1.fixed && !p2.fixed){
p2.x -= 0.5 * xDirection * distanceError;
p2.y -= 0.5 * yDirection * distanceError;
p1.x += 0.5 * xDirection * distanceError;
p1.y += 0.5 * yDirection * distanceError;
}
}
}
}That’s it! Please see the GitHub Gist for the whole code. Our simulator is ready to be used. To show its functionality, I implemented a Python interface using pybind11 and created an animated plot using matplotlib and plot2vid.
Turns out that the technique I described is much more powerful than just ropes. You can easily transform this code to create more fun stuff, including cloths and solid bodies. For inspiration, check out this website.
