In Update loop a square hits the edge and turns around. An if is enough for that. As soon as two moving bodies meet, it isn't: who gets how much speed, which way does the crate tip, and why does a stack stay up in the first place?
The demo has circles and boxes, a floor, two side walls and a ramp. At the start a pyramid and a tower are already standing there, the rest is still falling in.
A click on empty space spawns a body, dragging a body pulls it through the scene on a spring. The chips switch restitution, gravity and the contact points (orange), R clears the scene. The pyramid stands until the small circle lands on its tip. It arrives at around 10 m/s, after that the pyramid is rubble.
Underneath runs the loop from Fixed timestep: the world takes steps of exactly 1/120 second out of an accumulator, which at 60 fps means two steps per frame. The count is shown in the bar at the top right.
Where, how deep, which direction
A contact consists of a normal pointing from body A to body B, a penetration depth and a point. For two circles that is cheap:
const Vec2 d = b.pos - a.pos;
const float dist = length(d);
const float radii = a.radius + b.radius;
if (dist >= radii) return m;
m.normal = dist > 1e-6F ? d * (1.0F / dist) : Vec2{0.0F, 1.0F};
m.contacts[0].separation = dist - radii;
separation is negative as long as the bodies overlap. Circle against box works in the box's coordinate system: clamp the circle center to the box, the distance to that point is the distance to the box. If the center is inside the box, the circle leaves through the nearest face.
Box against box is the big one. By the separating axis theorem two convex shapes only overlap if they overlap on every edge normal, for two boxes that is four axes. The axis with the smallest overlap becomes the contact normal. For the points you take the edge of one box that points most against the normal and clip it against the side planes of the other box. What lies behind the reference face becomes a contact, so at most two points. This is the method of Box2D Lite, I followed it when building this. The boxes are real rotated rectangles, no simplified variant.
The impulse
At the contact point only one thing counts: how fast the two bodies approach each other along the normal. The relative velocity at the point includes the rotation: vel + angVel x r, with r as the offset from the center of mass to the contact. The impulse that brings it to the wanted value:
j = -(1 + e) * dot(vRel, n) / (1/mA + 1/mB + cross(rA, n)^2 / IA + cross(rB, n)^2 / IB)
On top is the relative velocity, e is the restitution (0: the hit swallows everything, 1: fully elastic). On the bottom is the effective mass. Because the inertia I is part of it, a hit far from the center of mass pushes the crate away less than it spins it. In code that is the reciprocal, computed once per contact:
c.massNormal = 1.0F / (a.invMass + b.invMass + a.invInertia * rn1 * rn1 + b.invInertia * rn2 * rn2);
Two circles of equal mass with e = 1 swap their velocities this way, against a static wall (invMass 0) the velocity is mirrored. Both are tests in test.cpp, calculated by hand.
Two things I had to add. Restitution only applies above 1 m/s impact speed, otherwise every crate at rest would bounce: gravity alone adds 0.08 m/s per step. And the solver does not remove the overlap through the impulse alone. On top comes a small target speed, 20 percent of the depth beyond 1 cm of tolerance (slop) per step, at most 4 m/s. The slop lets bodies sink in a tiny bit so a contact doesn't flip between touching and apart in every step.
Friction
The same calculation once more, along the tangent (the normal turned by 90 degrees). It is capped by the normal impulse times the friction coefficient:
const float maxFriction = friction * c.pn; c.pt = std::clamp(oldPt + c.massTangent * -vt, -maxFriction, maxFriction);
The coefficient is fixed at 0.5. A box sliding over the floor at 3 m/s brakes with μ·g, about 4.9 m/s², and stops after a little over half a second.
Why stacks shake
The contacts depend on each other: the bottom crate carries everything above it, and every impulse at one contact disturbs its neighbors. The solver corrects the contacts one after another (sequential impulses) and repeats the round, eight times per step here. That approaches the right answer without ever reaching it. Warm starting helps a lot: the impulse a contact needed in the last step is applied again right at the start of the next one. A contact is recognized by its pair of bodies and a point within 10 cm.
I measured it with the demo's code, 10 seconds of simulation, boxes of 1 x 1 m. The tower of five sinks 1 cm at the top with eight iterations and warm starting (that is the slop) and doesn't move after that. Without warm starting it is 5 cm and it still twitches at 3 cm/s. If I also drop the position correction, it falls apart. For the tower of ten, a single iteration brings it down, with two it stands but twitches at up to 6 cm/s, from four on it is calm. Without warm starting the tower of ten topples even at eight iterations.
What's missing
Nothing is interpolated, the last step is what gets drawn. Bodies never sleep, and skipping resting bodies saves a lot in big scenes. There is no continuous collision detection, so very fast bodies can punch through thin walls. The collision search checks every pair with a bounding circle, which is enough for the demo's 90 bodies at most.
The code is in the repo under rigid-bodies. The project Game Dev Snippets collects all demos of the series.