Physics Simulation: Rigid Bodies, Soft Bodies & Collision Detection | Game Development Hub

Introduction Physics simulation brings virtual worlds to life. This guide covers the algorithms and data structures behind production physics engines: rigid body dynamics with constraint solvers, hierarchical collision detection (broad/narrow phase), and soft body simulation techniques. Rigid Body Dynamics Rigid bodies are the foundation — objects that maintain fixed distances between internal points. State consists of position, orientation (quaternion), linear velocity, and angular velocity. Equations of Motion // Linear: F = ma → a = F/m v += a * dt x += v * dt // Angular: τ = Iα → α = I⁻¹τ ω += α * dt q += 0.5 * quat(0, ω) * q * dt q = normalize(q) Inertia Tensor The 3×3 inertia tensor (in body space) relates angular velocity to angular momentum: // Body-space inertia for common shapes (mass m) // Box (w×h×d): I = m/12 * diag(h²+d², w²+d², w²+h²) // Sphere (r): I = 2/5 * m * r² * Identity // Cylinder (r, h): I = m/12 * diag(3r²+h², 3r²+h², 6r²) // World-space inertia: I_world = R * I_body * R^T Integration Methods Method Order Stability Use Case Euler (explicit) 1 Poor Simple, non-stiff Semi-implicit Euler 1 Good Game standard (symplectic) Velocity Verlet 2 Excellent Energy conservation Runge-Kutta 4 4 Good Offline/simulation Semi-implicit Euler (symplectic) is the game industry standard — stable for stiff springs, conserves energy reasonably, and is trivial to implement. Collision Detection: Broad Phase Broad phase quickly eliminates non-colliding pairs using simple bounds. Target: O(n log n) or better for n bodies. Sweep and Prune (SAP) Sorts axis-aligned bounding box (AABB) intervals on each axis. Overlap on all three axes = potential collision. struct AABB { float3 min, max; }; // Sort endpoints on X axis struct Endpoint { float value; int bodyId; bool isMin; }; vector endpointsX; // ... fill and sort ... // Sweep int activeCount = 0; for (auto& ep : endpointsX) { if (ep.isMin) { // Check overlap with active bodies for (int other : activeBodies) { if (aabb[ep.bodyId].overlapsY(aabb[other]) && aabb[ep.bodyId].overlapsZ(aabb[other])) { addPair(ep.bodyId, other); } } activeBodies.add(ep.bodyId); } else { activeBodies.remove(ep.bodyId); } } Uniform Grid / Spatial Hash Divide world into uniform cells. Insert bodies into overlapping cells. Check pairs within each cell. // Spatial hash for dynamic worlds unordered_map > grid; uint64_t hashCell(int3 coord) { return ((uint64_t)(coord.x + 10000) << 42) | ((uint64_t)(coord.y + 10000) << 21) | (uint64_t)(coord.z + 10000); } void insertBody(int bodyId, AABB bounds) { int3 minCell = floor(bounds.min / cellSize); int3 maxCell = floor(bounds.max / cellSize); for (int x = minCell.x; x <= maxCell.x; x++) for (int y = minCell.y; y <= maxCell.y; y++) for (int z = minCell.z; z <= maxCell.z; z++) { grid[hashCell({x,y,z})].push_back(bodyId); } } Bounding Volume Hierarchy (BVH) Tree of bounding volumes (AABBs or spheres). Built bottom-up (SAH heuristic) or top-down. Excellent for static geometry, refit for dynamic. Collision Detection: Narrow Phase Narrow phase computes exact contact information (point, normal, penetration depth) for pairs from broad phase. GJK (Gilbert-Johnson-Keerthi) Iterative algorithm for distance between convex shapes using Minkowski difference and support mapping: // Support function for shape float3 support(Shape* shape, float3 dir) { // Return point furthest in direction dir return shape->getFurthestPoint(dir); } // GJK distance query bool gjk(Shape* A, Shape* B, float3& closestPoint) { Simplex simplex; float3 dir = float3(1, 0, 0); // Initial direction for (int iter = 0; iter < MAX_ITER; iter++) { float3 a = support(A, dir); float3 b = support(B, -dir); float3 p = a - b; // Minkowski difference point if (dot(p, dir) <= 0) return false; // Origin not in Minkowski sum simplex.add(p); if (simplex.containsOrigin(dir)) { closestPoint = simplex.closestPointToOrigin(); return true; } } return false; } EPA (Expanding Polytope Algorithm) After GJK detects intersection, EPA finds penetration depth and contact normal by expanding the simplex: // EPA expands simplex faces toward origin // Returns penetration depth and contact normal bool epa(Simplex& simplex, Shape* A, Shape* B, float& depth, float3& normal, float3& contactPoint) { // Build initial polytope from GJK simplex // Iteratively expand face closest to origin // until convergence } SAT (Separating Axis Theorem) For polyhedra: test separation on face normals of both shapes + cross products of edges. Simpler than GJK/EPA for box-box, triangle-mesh. Constraint Solvers Constraints enforce relationships: contact non-penetration, joints, motors. Solved as linear complementarity problem (LCP) or via iterative impulse methods. Sequential Impulse (PGS — Projected Gauss-Seidel) Iteratively solve each constraint, applying impulses to satisfy it. Order-independent convergence for well-conditioned systems. // Contact constraint: relative velocity at contact along normal ≥ 0 // J * v = -b (constraint equation) // Impulse: P = J^T * λ float solveContact(Body* A, Body* B, ContactManifold& manifold) { for (int iter = 0; iter < SOLVER_ITERATIONS; iter++) { for (Contact& c : manifold.contacts) { // Compute effective mass float3 rA = c.point - A->position; float3 rB = c.point - B->position; float3 Jn = c.normal; float3 JrA = cross(rA, c.normal); float3 JrB = cross(rB, c.normal); float invMass = A->invMass + B->invMass + dot(JrA, A->invInertiaWorld * JrA) + dot(JrB, B->invInertiaWorld * JrB); // Relative velocity at contact float3 vA = A->velocity + cross(A->angularVelocity, rA); float3 vB = B->velocity + cross(B->angularVelocity, rB); float vRel = dot(vB - vA, c.normal); // Bias (Baumgarte stabilization) float bias = -BAUMGARTE * c.penetration / dt; // Compute lambda float lambda = -(vRel + bias) / invMass; lambda = max(lambda + c.accumulatedImpulse, 0) - c.accumulatedImpulse; c.accumulatedImpulse += lambda; // Apply impulse float3 impulse = c.normal * lambda; A->applyImpulse(-impulse, rA); B->applyImpulse(impulse, rB); } } } XPBD (Extended Position-Based Dynamics) Modern alternative: solve constraints directly on positions, more stable for stiff materials. Used in Unity DOTS Physics, Bevy XPBD. // XPBD constraint projection // x_new = x_old + Δx // Δx = -α * ∇C(x) * (C(x) / |∇C(x)|²) // α = compliance = 1/stiffness void projectConstraint(Body* A, Body* B, Constraint& c, float dt) { float3 grad = c.gradient(A, B); // ∇C float C = c.value(A, B); // Constraint violation float alpha = c.compliance / (dt * dt); float denom = dot(grad, A->invMass * grad) + dot(grad, B->invMass * grad) + alpha; float deltaLambda = -C / denom; A->position += A->invMass * grad * deltaLambda; B->position -= B->invMass * grad * deltaLambda; } Soft Body Simulation Soft bodies deform — cloth, jelly, flesh. Represented as mass-spring systems, FEM tetrahedra, or position-based dynamics particles. Mass-Spring Systems Particles connected by springs (structural, shear, bend). Simple but suffers from stiffness instability. struct Spring { int a, b; float restLength; float stiffness; float damping; }; void solveSprings(vector & particles, vector & springs, float dt) { for (auto& s : springs) { Particle& pA = particles[s.a]; Particle& pB = particles[s.b]; float3 dir = pB.pos - pA.pos; float len = length(dir); float3 force = dir * (s.stiffness * (len - s.restLength) / len); float3 dampForce = (pB.vel - pA.vel) * s.damping; pA.applyForce(force + dampForce); pB.applyForce(-force - dampForce); } } FEM (Finite Element Method) Tetrahedral mesh with constitutive model (Neo-Hookean, StVK). Computes strain energy density, derives forces from energy gradient. Accurate but expensive. PBD/XPBD for Cloth & Soft Bodies Position-Based Dynamics: directly project positions to satisfy distance constraints. Unconditionally stable, artist-friendly parameters. // Distance constraint (cloth edge) float C = length(pB - pA) - restLength; float3 grad = normalize(pB - pA); float invMassSum = pA.invMass + pB.invMass; float deltaLambda = -C / (invMassSum + compliance/dt²); pA.pos += pA.invMass * grad * deltaLambda; pB.pos -= pB.invMass * grad * deltaLambda; // Volume constraint (tetrahedron for soft body) // C = volume(tet) - restVolume // Gradient w.r.t. 4 vertices via scalar triple product Conclusion Physics engines combine rigid body dynamics (constraint solvers), hierarchical collision detection (SAP/BVH broad phase + GJK/EPA narrow phase), and soft body methods (PBD/XPBD for stability). Modern engines like PhysX, Box2D, Jolt, and Bevy XPBD use iterative impulse solvers with position correction (XPBD) for robustness. Choose algorithms based on your simulation needs: rigid bodies for gameplay, PBD for cloth/soft bodies, FEM for offline accuracy.