Making a 2D soft-body physics engine

Recently, I participated in Ludum Dare 53 and made a silly spaceship building game based on a soft-body physics engine.

I figured I'd share how this physics engine works, since it's actually remarkably simple :)

When making a physics engine, there is always a question of whether you simulate actual solid objects like boxes, disks, capsules, etc, or stick to having just point masses instead. Here are a few things to consider:

- The state of a solid object consists of its
**position**,**velocity**,**orientation**, and**angular velocity**. For point masses, you don't need the rotational stuff, just**position**and**velocity**.*(Rotation doesn't make sense for single points: the only point of this object is its center of mass, and the rotation of center of mass around itself is zero. The inertia tensor of a point mass is zero, etc.)*Rotation is tricky to handle, especially when your objects collide or you apply forces to specific points on an object (e.g. push by player, nearby explosion, etc). - Solid objects are more
**physically correct**. We can emulate e.g. a box with four point masses together with some constraints that force them to maintain a box-like shape, but constraints are also tricky to do right and will usually be voilated from time to time, meaning your box won't be a box. Also four connected vertices have mass concentrated in the vertices, therefore it will have an inertia tensor different to that of a box, and this will affect how our improvised box reacts to forces & rotation. Imagine a shallow wooden box compared to solid filled polymeric foam box – the difference is quite similar to our case. - Collisions between point masses are basically meaningless, unless we pretend that our "points" actually have a radius and treat them as disks just for collisions (they're still not
*rotating disks*, though, and you will notice that they behave a bit weird because of this). Collisions between point masses and static environment, on the other hand, are as easy as collisions go! Compare this to collisions between solids: they're easy if you only support disks (still need to correctly handle friction & rotation & etc), but if you want boxes/capsules/polygons you need to write some pretty non-trivial collision detection code, and in any case the collision response is quite involved.

So, neither option is better or worse, they're just different and suit different use cases. Most physics engines are built on solid objects because that's a more universal (albeit much more complicated) approach. In my case, however, I wanted a *soft-body* engine, meaning my so-called solid objects will deform freely and stop being solids. I'm sure there are ways to handle this inside a typical solid object-based engine, but I decided to use the point masses instead: they are simpler to code, and they already allow for arbitrary deformations, just need to figure out how to force the points to loosely maintain some shape.

Let's write a 2D point mass physics engine without any constraints or extra forces. It's really just a few lines of code:

```
struct point_mass
{
vec2 position;
vec2 velocity;
};
struct engine
{
std::vector<point_mass> points;
vec2 gravity;
void update(float dt)
{
for (auto & p : points)
p.velocity += gravity * dt;
for (auto & p : points)
p.position += p.velocity * dt;
}
}
```

First, we update the velocities using the external forces (called `gravity`

here, but it can be any external force, and it can differ for different objects / points in space). Then, we update the positions using the velocities, and that's it. This is the symplectic Euler integration method, known for good energy conservation properties. It is also remarkably simple; simpler than the universally beloved Vertet method.

Right now, we have something like this:

Since I'm planning to use these point masses only as building blocks for soft bodies, I don't care about collisions between two point masses. However, I still need collisions between points and the environment.

For the collision to work, we need a routine that computes the *collision normal* and *penetration depth* of the collision. For example, if my "environment" is just a floor below the Y=0 line, its collision normal points upwards, and the penetration depth for a point at (X, Y) is simply -Y:

```
struct collision
{
vec2 normal{0.f, 0.f};
float depth = -std::numeric_limits<float>::infinity();
};
collision find_collision(vec2 const & position)
{
return collision{vec2{0.f, 1.f}, -position.y};
}
```

Negative penetration depth means there is no collision.

My game actually only uses collisions with planets, i.e. with static disks, which looks like this:

```
struct planet
{
vec2 center;
float radius;
};
collision find_collision(vec2 const & position, planet const & planet)
{
vec2 delta = position - planet.center;
float distance = length(delta);
vec2 normal = delta / distance;
float depth = planet.radius - distance;
return collision{normal, depth};
}
```

If we have several colliding obstacles, we simply find the collision with the largest penetration depth – it is probably the most important one!

```
collision find_collision(vec2 const & position, std::vector<planet> const & planets)
{
collision result;
for (auto const & planet : planets)
if (auto c = collision(position, planet); c.depth > result.depth)
result = c;
return result;
}
```

Notice that we've initialized the depth to `-inf`

in the `collision`

constructor, so that a default-constructed `collision`

acts as no collision, and its `depth`

value is the neutral element with respect to `max`

.

I'm very annoyed that most resources on physics engine collisions only mention collision detection. Like, it's usually the easy part! Especially with rotating solids. Especially in 3D!

Anyway, we got our collision normal & penetration depth. There are three things we need to do:

- Push the point mass along the collision normal so that no collision is present. Or, in fancy language, prevent collision constraint violation. This is important to do because even if we do the velocity update described next, something (external forces, springs, soft body constraints, etc) might still push our object inside the collider. This is known as
*sinking*, or more generally as*constraint drifting*: we have fixed the*time derivative*of the constraint function, but we didn't fix the constraint itself. This pushing along collision normal can be problematic if you have some interaction forces, e.g. gravitational attraction: changing the position of an object changes the potential energy of this interaction, meaning our collisions will randomly add or remove energy to/from the system. Yikes!*(see my another post where I explore the specific problem of collision resolution messing up with gravity simulation)*. There is an alternative option to only resolve collision if the object is moving towards the collider, not away from it (i.e. if the dot product of the object's velocity and the collision normal is negative), but it only helps in very simple situations. - Update the normal component of the velocity. If you throw a ball towards a wall, it will bounce back. This is what things typically do when they collide something: they bounce back, maybe losing some energy in the process. This is usually captured as the
*elasticity coefficient*, or*bounciness*: elasticity = 0 means the object will stop upon collision, elasticity = 1 means the object will retain its kinetic energy in full. - Update the tangentinal component of the velocity, or, in simple terms, apply friction. Our point masses don't rotate, so this will be relatively easy: we multiply the tangential component by a factor of the form
`(1 - friction * dt)`

, or a much more stable & precise alternative`exp(- friction * dt)`

*(see another post of mine for an explanation of this formula)*. This means that while our point mass is in contact with the collider, it's speed*along*the collider will get smaller over time. Notice that while elasticity is a unitless coefficient in the range [0..1], friction has units`1/TIME`

and can be any positive number (at least if you use the stable formula, otherwise it cannot be larger than`1/dt`

) and 1/friction is roughly*the time it takes for the object to drop 63% of it's speed*, something like that.*(0.63 is the value of*.`1 - exp(-1)`

; again see this post to learn where this comes from)

Enough talking, here's the code:

```
void engine::update(float dt)
{
// ... velocity integration ...
for (auto & p : points)
{
collision c = find_collision(p.position, the_world);
// check if collision took place
if (c.depth < 0.f) continue;
// resolve the constraint
p.position += c.normal * c.depth;
// compute the normal & tangential velocity
auto vn = c.normal * dot(c.normal, p.velocity);
auto vt = p.velocity - vn;
// apply bouncing
vn = - elasticity * vn;
// apply friction
vt *= std::exp(- friction * dt);
// add up the new velocity
p.velocity = vn + vt;
}
}
```

Just a few formulas here and there and we get this:

Collisions are handled as if the points have a nonzero radius, but otherwise the logic is the same.

Say, I want to connect two objects in my physics engine with an invisible stick of a certain length. Or maybe I want to *constrain* the movement of an object to only happen along some predefined curve. This is what *constraints*, also known as *joints*, are for. Usually we're interested in *equality constraints*, meaning we want some function of our system to be equal to zero: \( f(\text{system}) = 0 \). For example, if I want a point mass to only move along the X=Y line, I'd do \( f(\text{system}) = p_X - p_Y \). If I want the distance between two points to be fixed to some value R, I'd do \( f(\text{system}) = |p_0 - p_1| - R \), etc.

Collisions are an example of an *inequality constraint* \( f(\text{system}) \leq 0 \). We resolve the collision if the penetration depth is positive, but once it is negative, we don't care. Some engines (like Box2D) treat collisions together with other constraints in a single constraint resolution framework; others (like Bullet, iirc) treat them as a separate thing entirely.

Constraint resolution is a pretty complicated topic in general, but in some cases we can hack them a bit. For example, to resolve the distance constraint between a pair of point masses, we can simply move the two points along the line between them so that the distance is as required:

```
struct distance_constraint
{
uint32_t index0, index1;
float distance;
};
struct engine
{
// ...
std::vector<distance_constraint> constraints;
};
void engine::update(float dt)
{
// ... velocity integration ...
// ... collision resolution ...
for (auto const & c : constraints)
{
auto & p0 = points[c.index0].position;
auto & p1 = points[c.index1].position;
auto delta = p1 - p0;
auto distance = length(delta);
auto required_delta = delta * (c.distance / distance);
auto offset = required_delta - delta;
p0 -= offset / 2.0;
p1 += offset / 2.0;
}
}
```

This is an example of a *hard constraint*: we force the constraint to be always satisfied, to behave as a *rigid* thing.

We could instead move the objects bound by a constraint a bit every frame, so that the constraint gets satisfied if we wait long enough:

```
void engine::update(float dt)
{
// ... velocity integration ...
// ... collision resolution ...
for (auto const & c : constraints)
{
auto & p0 = points[c.index0].position;
auto & p1 = points[c.index1].position;
auto delta = p1 - p0;
auto distance = length(delta);
auto required_delta = delta * (c.distance / distance);
float damping_factor = 1.f - std::exp(- constraint_damping * dt);
auto offset = (required_delta - delta) * damping_factor;
p0 -= offset / 2.0;
p1 += offset / 2.0;
}
}
```

Or we could use *force-based* constraints instead: apply forces to the constraint-bound objects that would move the system towards a state when the constraint is satisfied. For the distance constraint, the corresponding force-based thing is called a *spring*:

```
void engine::update(float dt)
{
// ... velocity integration ...
// ... collision resolution ...
for (auto const & c : constraints)
{
auto p0 = points[c.index0].position;
auto p1 = points[c.index1].position;
auto delta = p1 - p0;
auto distance = length(delta);
auto required_delta = delta * (c.distance / distance);
auto force = spring_force * (required_delta - delta);
points[c.index0].velocity -= force * dt;
points[c.index1].velocity += force * dt;
}
}
```

And this is what we get:

Much better for a *soft-body* engine! However, we probably don't want it to oscillate indefinitely, so we need some damping. This is a bit tricky: we only want to dampen the wiggling along the spring, but we don't want to touch the movement of these two point masses as a whole, and we want to preserve its rotation.

To achieve this, we compute the relative velocity along the spring's direction, compute the dampened velocity, and split the velocity correction between the two point masses:

```
void engine::update(float dt)
{
// ... velocity integration ...
// ... collision resolution ...
for (auto const & c : constraints)
{
auto p0 = points[c.index0].position;
auto p1 = points[c.index1].position;
auto & v0 = points[c.index0].velocity;
auto & v1 = points[c.index1].velocity;
auto delta = p1 - p0;
auto distance = length(delta);
auto direction = delta / distance;
auto required_delta = direction * c.distance;
auto force = spring_force * (required_delta - delta);
v0 -= force * dt;
v1 += force * dt;
auto vrel = dot(v1 - v0, direction);
auto damping_factor = exp(- spring_damping * dt);
auto new_vrel = vrel * damping_factor;
auto vrel_delta = new_vrel - vrel;
v0 -= vrel_delta / 2.0;
v1 += vrel_delta / 2.0;
}
}
```

With this damped spring, we get a nice soft distance constraint:

We can crank up the constants a bit and get something that looks like a rigid, hard constraint, but is still a soft constraint:

Of course, soft constraints have their downsides. They are *force-based*, so if there is a stronger force in the system (gravity, or another constraints), the constraint will be violated a lot. There are certain bounds on the constraint force (something like `spring_force * dt * dt < 1`

), otherwise the system will be unstable. But for a soft-body physics engine, they work pretty great!

So, how do we use all that information to make a soft-body physics engine? Like, how do we make a soft box? The first idea that comes to mind is to connect all vertices of the box with springs, and it indeed works:

but there is a problem with this approach: it does exactly what it says, meaning it only maintains pairwise distances between the points and doesn't care about their overall shape. In particular, this box made of 6 springs can be easily turned inside-out to form a different box:

and this spring-based "wheel" can be messed up so bad that is starts to levitate on it's own:

We could add extra constraints on the *angles* between consecutive points, but this would only hide and overcomplicate the problem. *(However, this is very common in molecular dynamics simulations; they even use the angle between the two planes formed by four consecutive points!)*

It turns out that there is a better way of simulating soft bodies.

I don't know a well-established name for this, and a quick google search failed to reveal anything of releavance, so I will call this method *shape matching*. If you know some resources on this, I would love to know them, since I had to derive all the equations myself :) It doesn't only apply to soft bodies or physics simulations, though: it is a much more general technique.

So, we want our points to *preserve their shape*. By *shape* we mean some solid object – a box or a polygon, for example. We want our points to form the vertices of this shape. However, in general they won't do that due to various forces acting on them.

Instead, let's do this: given the *current positions* of our points, let's find the *ideal position* of our solid shape so that it matches the *current shape* as much as possible. Since the *ideal shape* is solid, it can only move and rotate as a whole. *If I'm not mistaken, JellyCar physics engine uses the same concept.*

*If you don't care about the formula derivations, you can safely skip this part!*

Let's say that we have N points, and our *ideal shape* is specified by the vectors \( q_i \) from the shape's center of mass towards the shape's vertices. Our *current points* are \( p_i \). First, let's compute the *current* center of mass:

Then, compute the *current* points relative to the center of mass:

It would be wrong to expect \( r_i = q_i \) because our *ideal shape* can be *rotated*! Let \( R_\theta \) be the operator of rotation by an angle \( \theta \):

What we want is \( r_i = R_\theta\cdot q_i \), i.e. that the rotated *ideal shape* coincides with the *current shape*. This will, in general, not be true (again because of various other forces acting on our point masses), so we need to find a way to *force* it be true.

We can employ a classic trick of turning an equation into a minimization problem: let's try to minimize the sum of squared distances between the *current positions* and the *ideal positions*; the *ideal shape* that minimizes this value will be the shape we are looking for!

I've called it E because it resembles some energy formula, like for a harmonic oscillator.

The only unknown in the above equation is \( \theta \), so we need to differentiate the above formula with respect to this angle and set the derivative be equal to zero:

\[ \frac{dE}{d\theta} = \sum 2(r_i-R_\theta\cdot q_i)\cdot\left(-\frac{dR_\theta}{d\theta}\cdot q_i\right) \]We're using the fact that the derivative of a squared length of a vector can be expressed using a dot product of the vector and it's derivative:

\[ \frac{dA^2}{d\theta} = 2A\cdot \frac{dA}{d\theta} \]The expression for \( \frac{dE}{d\theta} \) above looks a bit intimidating, but don't fret, we'll deal with it step by step.

First, ignore the factor of 2, it can be simply divided away (i.e. dividing both sides by 2).

Next, each summand is a dot product of two vectors: \( r_i - R_\theta\cdot q_i \) and \( -\frac{dR_\theta}{d\theta}\cdot q_i \), the latter is a product of a matrix \( -\frac{dR_\theta}{d\theta} \) and a vector \( q_i \).

It so happens that \( \frac{dR_\theta}{d\theta} = R_{\theta+\frac{\pi}{2}} \), meaning this matrix rotates by an angle of \( \theta \) plus 90 degrees. This also means that for any vector \( v \)

\[ (R_\theta\cdot v)\cdot \left(\frac{dR_\theta}{d\theta} \cdot v\right) = 0 \]because these two vectors are rotated by a straight angle, i.e. by 90 degrees, meaning they are orthogonal, and their dot product is zero.

What this means for us is that in each summand of our big formula, the \( (R\theta\cdot q_i)\cdot (\frac{dR\theta}{d\theta} \cdot q_i) \) part is zero, and we are only left with \( r_i \cdot \left(\frac{dR_\theta}{d\theta}\cdot q_i\right) \) (I've removed the minus sign which also doesn't add anything now).

So, we're left with this equation:

\[ \sum r_i \cdot \left(\frac{dR_\theta}{d\theta}\cdot q_i\right) = 0 \]We probably could do some more tricks to derive a concise equation for \( \theta \) in matrix form, but I'll just expand everything in coordinates instead:

\[ \frac{dR_\theta}{d\theta} = \begin{pmatrix} -\sin\theta & -\cos\theta \\ \cos\theta & -\sin\theta \end{pmatrix} \] \[ \frac{dR_\theta}{d\theta}\cdot q_i = \begin{pmatrix} -q_i^X\sin\theta -q_i^Y\cos\theta \\ q_i^X\cos\theta -q_i^Y\sin\theta \end{pmatrix} \] \[ r_i \cdot \left(\frac{dR_\theta}{d\theta}\cdot q_i \right) =r_i^X \cdot (-q_i^X\sin\theta -q_i^Y\cos\theta) + r_i^Y \cdot (q_i^X\cos\theta -q_i^Y\sin\theta) = \] \[ = \sin\theta \cdot (-r_i^Xq_i^x - r_i^Yq_i^Y) + \cos\theta\cdot (-r_i^Xq_i^Y + r_i^Yq_i^X) = -\sin\theta\cdot(r_i \cdot q_i) - \cos\theta \cdot (r_i \times q_i) \]Here, a cross \( r_i\times q_i \) means the exterior product, also known as the *skew product* or the *pseudoscalar product*. It is the 2D analogue of a cross product and equals the signed area of a parallelogram spanned by two vectors.

Our final equation now looks like this:

\[ -\sin\theta\cdot(\sum r_i \cdot q_i) - \cos\theta \cdot (\sum r_i \times q_i) = 0 \]Which can be solved easily:

\[ \tan\theta = -\frac{\sum r_i\times q_i}{\sum r_i\cdot q_i} \] \[ \theta = - \text{atan2}\left(\sum r_i\times q_i, \sum r_i\cdot q_i\right) \]That's an extremely neat formula, if you ask me!

So we have a formula for the angle \( \theta \), meaning we can compute the *ideal positions* of our points that are as close to the *current positions* as possible. What do we do with them, though? We could simply move the points to these ideal positions, but that would be a *hard constraint* (although it would work much better than a distance constraint-based thing).

To make it into a *soft body*, we can add a spring-like force for each point mass towards it's *ideal position*:

```
float cross(vec2 a, vec2 b)
{
return a.x * b.y - a.y * b.x;
}
struct soft_body
{
struct vertex
{
uint32_t index;
vec2 position;
};
std::vector<vertex> vertices;
};
struct engine
{
// ...
std::vector<soft_body> soft_bodies;
};
void engine::update(float dt)
{
// ... velocity integration ...
// ... collision resolution ...
for (auto const & b : soft_bodies)
{
// compute the center of mass
vec2 center = vec2(0.f, 0.f);
for (auto const & v : b.vertices)
center += points[v.index].position;
center /= (float)b.indices.size();
// compute the shape rotation angle
float A = 0.f, B = 0.f;
for (auto const & v : b.vertices)
{
auto r = points[v.index].position - center;
A += dot(r, v.position);
B += cross(r, v.position);
}
float angle = -atan2(B, A);
// apply spring forces
for (auto const & v : b.vertices)
{
auto target = center + rotate(v.position, angle);
auto delta = target - points[v.index].position;
points[v.index].velocity += spring_force * delta * dt;
}
}
}
```

This works great. Here's how it looks if every `soft_body`

is just a square containing 4 points (note that different bodies can share points, allowing for complex structures made of several soft bodies):

However, we see that it wiggles forever, because our springs aren't damped. In my implementation, I added damping between any pair of consecutive vertices of a soft body, but I think it is better to compute the velocity of the center of mass (i.e. just the average velocity of the points of a soft body), compute the average angular velocity of rotation around the center of mass, compute the *target velocity* of each point mass treating the whole soft body as a solid body, and then add damping by slowly interpolating the velocities towards these target velocities.

These two methods of adding damping seem to produce quite similar results:

As usual, I wanted to explain something small but ended up writing a whole tutorial :) Hope you liked it. Check out my other article generalizing this technique to 3D. And, as always, thanks for reading!