9:03 PM
Continuation of this post
Previously, we created a basic particle-based physics engine. Each point was stored as simply a position and velocity, and iterated to produce natural-looking motion.
How can we extend this to more complex shapes?
Before we even try to simulate such shapes, we'll need a way to store them. A simple point can be represented fully with a position and velocity, as it has no defined rotation. To represent rotation and angular velocity, we just have to add two more fields to our collider struct:
struct collider { float px, py, vx, vy; // rotation = angular velocity float angle, rotation; };
It turns out that we'll be doing a lot of transformations on the position and velocity fields, and it'll be convenient to represent them as 2d vectors. We'll make a simple vec2 struct wrapping two floats.
struct vec2 { float x, y; };
struct collider { struct vec2 position, velocity; float angle, rotation; };
In order to use angle and rotation properly, we increment each collider's angle by its rotation in iterate(). Also, we can use some helper vector math functions instead of manually adding velocity to position.
void iterate (struct collider *c) { // vadd(a, b) returns the sum of two vectors c -> position = vadd(c -> position, c -> velocity); c -> angle += c -> rotation; }
But wait, this isn't everything! Now that we're working with arbitrarily shaped colliders instead of points, each collider can have a different shape, and that will affect its behavior. So, we need to represent the collider's shape in the struct.
I chose to define the shape of a collider by representing it as a collection of vertices. This isn't technically accurate for friction computations, but will do for now.
struct collider { size_t vertex_count; struct vec2 *vertices;
struct vec2 position, velocity; float angle, rotation; };
Ok, looks good now. Let's make a simple demo of a rotating square moving across the screen. We'll ignore collisions with the wall or any of that for now.
We don't actually have to change much from our previous point demo: iterate() is updated to work with rotation, so all we have to make is a helper function to display() our collider.
To make this simple, let's just connect all adjacent vertices with lines, and draw a point in the center to show the center of mass.
It'll also be convenient here to write a helper function screen_position(), which determines the screen-space position of any vertex. It will be useful for collision detection later.
struct vec2 screen_position (struct collider *c, size_t index) { // vrotate() does what it looks like - rotates a vector return vadd(c -> p, vrotate(c -> vertices[index], c -> angle)); }
void display (struct collider *c) { // assume collider to be properly centered around center of mass point((int) c -> position.x, (int) c -> position.y, 1);
// simply iterate through vertices and draw line to next one in sequence for (size_t vertex = 0; vertex < c -> vertex_count; ++vertex) { size_t next = vertex + 1; if (next == c -> vertex_count) next = 0;
struct vec2 p1 = screen_position(c, vertex), p2 = screen_position(c, next);
line((int) p1.x, (int) p1.y, (int) p2.x, (int) p2.y); } }
struct collider square;
square.vertex_count = 4; square.vertices = malloc(sizeof(struct vec2) * 4);
// nvec2() constructs a new vec2 square.vertices[0] = nvec2(-10.0, -10.0); square.vertices[1] = nvec2(-10.0, 10.0); square.vertices[2] = nvec2(10.0, 10.0); square.vertices[3] = nvec2(10.0, -10.0);
square.position = nvec2(50.0, height/2.0); square.velocity = nvec2(0.2, 0.0);
square.angle = 0.0; square.rotation = 0.03;
while (true) { iterate(&square);
clear(); display(&square); draw();
sleepms(20); }
Let's add collisions.
Conceptually the same logic applies as in the point simulation. When any vertex of the collider hits a wall, it exerts a force on the wall, and the same force is applied back to the collider.
To find the amount of force exerted, we multiply the collider's mass with the relative velocity between the colliding vertex and the wall. This is important if the collider is rotating. The location where force is applied matters, too, because it will change the rotation of the body.
Acceleration can be calculated with F = ma. Torque is calculated as the cross product of the force and offset, and rotational acceleration is torque divide moment of inertia.
Let's add two more fields to our collider struct to hold mass and radius (average distance from center of mass), and write a helper function apply_force() to keep things neat.
struct collider { size_t vertex_count; struct vec2 *vertices;
float mass, radius; struct vec2 position, velocity; float angle, rotation; };
// correct for units - velocity measured in pixels/frame, not pixels/second float damping = 60.0; // assume 60fps
void apply_force (struct collider *c, struct vec2 force, struct vec2 offset) { struct vec2 acceleration = vmult(force, 1.0 / c -> mass);
float torque = offset.x*force.y - offset.y*force.x, moment = c -> mass * c -> radius;
c -> v = vadd(c -> v, acceleration); c -> rotation += torque / moment / damping; }
Now, update iterate() to check and handle collisions. We essentially need to loop through all of the collider's vertices and check if they are contacting a wall. If so, compute the relative velocity, and call apply_force().
We already have a utility to find the screen-position of each vertex. Let's do the same, but for screen-velocity instead.
struct vec2 screen_velocity (struct collider *c, size_t index) { // velocity contributed by rotation is perpendicular // the vector from the center to vertex struct vec2 perpendicular = vrotate(c -> points[i], c -> angle + PI/2.0); return vadd(c -> velocity, vmult(perpendicular, c -> rotation)); }
We then add the following to iterate():
for (size_t vertex = 0; vertex < c -> vertex_count; ++vertex) { struct vec2 pos = screen_position(c, vertex), vel = screen_velocity(c, vertex);
// test for collision with walls if (p.x < 0.5 && v.x < 0.0) { struct vec2 force = nvec2(-v.x * c -> mass, 0.0), offset = vsub(pos, c -> position);
apply_force(c, force, offset); }
... }
Looks good! Let's try adding gravity. We can do the same as what we did in the particle engine - apply a constant force every frame.
Ah, that isn't good. Why does it happen?
Let's look at the code which calculates collisions with the bottom wall.
if (p.y > height - 0.5 && v.y > 0.0) { struct vec2 force = nvec2(0.0, -v.y * c -> mass), offset = vsub(pos, c -> position);
apply_force(c, force, offset); }
The restoring force is calculated as directly proportional to relative velocity. This causes zero restoring force to be applied if there is no relative motion. This works fine in a system with no external forces, but breaks when we add gravity.
In the real world, if an object is resting on a table or other surface, the surface exerts an upwards force on the object to counteract gravity. We can do the same here, by inserting a check into iterate():
if (p.y > height - 0.5) { struct vec2 force = nvec2(0.0, -0.07 * c -> mass), offset = vsub(pos, c -> position);
apply_force(c, force, offset); }
This fixes the clipping problem! However, we can see a few issues. Firstly, the box behaves as if there isn't any friction, and secondly, we see a bit of jitter.
We can add friction by modifying the collision handler - rather than only applying force directly perpendicular to the wall, also apply a small amount of force parallel to the wall, proportional to the relative horizontal velocity.
if (p.y > height - 0.5 && v.y > 0.0) { struct vec2 force = nvec2(-v.x * 0.3 * c -> mass, -v.y * c -> mass), offset = vsub(...
I'm using 0.3 as the coefficient of friction, but any value between 0 and 1 should work.
To increase stability we need to reduce numerical error. Generally, to mitigate the error of a newton iteration-based simulation, we reduce the step size. To do this, we can run iterate() multiple times per frame, and in each iteration, apply only a fraction of the collider's velocity to its position.
Let's add another parameter to iterate(): dt, the fraction of a full frame that each iteration simulates.
void iterate (struct collider *c, float dt) { c -> position = vadd(c -> position, vmult(c -> velocity, dt)); c -> angle += c -> rotation * dt; ...
while (true) { for (size_t i = 0; i < 20; ++i) { iterate(&square, 0.05); } ...
20 iterations is fine for our purposes, but more divisions will result in better stability.
Very nice.
tags: programming