Pavan K Rao

Custom Engine Devlog 4 | Umbra | Physics Part 3 Contact Resolution


game-engine umbra physics

Part 3 — Contact Resolution: Sequential Impulse and Constraint Solving

We now know where bodies collide and how deep they overlap. The solver’s job: compute impulses that prevent penetration, simulate friction, and do it all stably enough for boxes to stack.

The Sequential Impulse Method

UMBRA uses Sequential Impulse (SI), the same algorithm powering Box2D and most real-time physics engines. The idea:

  1. Precompute constraint data once per frame (effective masses, bias terms)
  2. Iterate over all contacts multiple times, applying small corrective impulses
  3. Each iteration improves the solution — with enough iterations, the system converges

This is essentially Gauss-Seidel iteration applied to the complementarity problem of contact constraints.

Precomputation: PrecomputeContactConstraints()

Before the solver iterates, we compute everything that doesn’t change between iterations:

Lever Arms

contact.rA = contactPoint - bodyA.Position;
contact.rB = contactPoint - bodyB.Position;

These vectors from body centers to the contact point determine how impulses create torque.

Effective Mass (Normal)

The “effective mass” along the contact normal accounts for both bodies’ masses and the rotational coupling through the lever arms:

1/m_eff = 1/mA + 1/mB + (rA × n)²/IA + (rB × n)²/IB

float rACrossN = Cross2D(rA, normal);
float rBCrossN = Cross2D(rB, normal);
float normalDenom = invMassSum
    + rACrossN * rACrossN * invInertiaA
    + rBCrossN * rBCrossN * invInertiaB;
contact.normalMass = 1.0f / normalDenom;

The (rA × n)² terms represent how much of a normal impulse at the contact point goes into rotation rather than translation. A contact far from the center of mass is “softer” because the impulse is partly absorbed by angular acceleration.

Effective Mass (Tangent)

The exact same formula but with the tangent direction (perpendicular to the normal):

Vector2f tangent(-normal.y, normal.x);
float rACrossT = Cross2D(rA, tangent);
// ... same pattern
contact.tangentMass = 1.0f / tangentDenom;

Restitution Velocity Bias

For bouncy collisions, we need to add energy. The bias is computed from the closing speed and the coefficient of restitution:

Vector2f velA = bodyA.Velocity + Vector2f(-rA.y, rA.x) * bodyA.AngularVelocity;
Vector2f velB = bodyB.Velocity + Vector2f(-rB.y, rB.x) * bodyB.AngularVelocity;
float closingSpeed = Dot(velB - velA, normal);
contact.velocityBias = closingSpeed < -1.0f ? -e * closingSpeed : 0.0f;

The -1.0f threshold prevents micro-bounces at rest — if the closing speed is tiny, we just absorb it completely. This is a critical stability feature; without it, resting contacts jitter endlessly.

The Velocity Solver: ResolveContacts()

This runs VelocityIterations times (default: 8). Each iteration loops over every contact and does two things:

Normal Impulse with Accumulated Clamping

// Relative velocity at contact point
Vector2f relVel = velB - velA;
float velAlongNormal = Dot(relVel, normal);

// Delta impulse this iteration
float dj = (-velAlongNormal + velocityBias) * normalMass;

// Accumulated clamping (the key to SI stability)
float oldAccum = contact.normalImpulseAccum;
contact.normalImpulseAccum = max(oldAccum + dj, 0.0f);
dj = contact.normalImpulseAccum - oldAccum;  // actual delta to apply

// Apply
bodyA.Velocity -= normal * dj * invMassA;
bodyA.AngularVelocity -= Cross2D(rA, normal * dj) * invInertiaA;
bodyB.Velocity += normal * dj * invMassB;
bodyB.AngularVelocity += Cross2D(rB, normal * dj) * invInertiaB;

Why accumulated clamping instead of per-iteration clamping?

Per-iteration clamping (dj = max(dj, 0)) can only add normal impulse. But when multiple contacts share bodies, correcting one contact can over-correct another. With accumulated clamping, the total impulse is clamped (max(total, 0)), allowing the delta to be negative — pulling back impulse that was applied in a previous iteration. This converges to the correct solution much faster.

Friction Impulse with Coulomb Clamping

After the normal impulse updates velocities, we recompute relative velocity and apply friction along the fixed tangent:

float velAlongTangent = Dot(relVel, tangent);
float djt = -velAlongTangent * tangentMass;

// Coulomb friction: |friction| <= mu * normalForce
float maxFriction = friction * contact.normalImpulseAccum;
float oldTangentAccum = contact.tangentImpulseAccum;
contact.tangentImpulseAccum = clamp(oldTangentAccum + djt, -maxFriction, maxFriction);
djt = contact.tangentImpulseAccum - oldTangentAccum;

The Coulomb friction model limits the tangential impulse to mu * normal_impulse. If the relative tangential velocity can be zeroed without exceeding this limit, it’s static friction (the objects stick). If the limit is hit, it’s dynamic friction (the objects slide). The accumulated clamping ensures both limits cooperate correctly across iterations.

The friction coefficient is the geometric mean of both bodies’ friction values:

float friction = sqrt(bodyA.DynamicFriction * bodyB.DynamicFriction);

Geometric mean (rather than arithmetic mean or min) means that a zero-friction surface is truly frictionless, while two high-friction surfaces are grippier than either alone.

Position Correction: PositionContraction()

The velocity solver prevents future penetration, but bodies may already be overlapping from the current frame. Position correction directly nudges bodies apart:

const float percent = 0.4f;  // correct 40% per iteration
const float slop = 0.01f;    // allow 1cm of overlap before correcting

float correctionMag = max(penetration - slop, 0.0f);
Vector2f correction = normal * (correctionMag / invMassSum) * percent;

bodyA.Position -= correction * invMassA;
bodyB.Position += correction * invMassB;

Two important stability parameters:

  • Slop (penetration allowance): without it, the solver fights numerical drift every frame, causing visible jitter. The 1cm slop means contacts are allowed to slightly overlap, which is invisible at normal scale.

  • Correction percentage (0.4): correcting 100% in one shot causes oscillation. Correcting a fraction and iterating (PositionIterations = 3) converges smoothly. This is essentially Baumgarte stabilization applied as direct position correction rather than velocity bias.

The correction is split between bodies proportional to their inverse masses — lighter bodies move more, and infinite-mass (static) bodies don’t move at all.

Sleeping Body Treatment

Throughout the solver, sleeping bodies are treated as having InverseMass = 0 and InverseInertia = 0:

float invMassA = bodyA.bIsSleeping ? 0.0f : bodyA.InverseMass;

This means a sleeping body acts like a static body in the constraint solver — it participates in collisions (pushing dynamic bodies away) but doesn’t move itself. This is cheaper than a full wake-up and prevents cascading wake-ups from minor contacts.

Why Multiple Iterations?

A single pass through all contacts produces a valid solution for each contact in isolation, but contacts share bodies. Fixing contact 1 might violate contact 2. Each additional iteration reduces the remaining error exponentially.

The default 8 velocity iterations is a practical sweet spot: stable stacking of 5-10 objects, responsive collisions, and no visible jitter. More iterations improve stability at the cost of CPU time — it’s a tunable quality knob.