Constraint Solver

The solver's first step is to mathematically define the physical limits of the system. This involves translating real-world restrictions into equality constraints (e.g., a weld joint must keep two bodies fixed together) and inequality constraints (e.g., contact points must not push objects into each other). Common formulations include:

• Linear Complementarity Problem (LCP): Models contact and friction where forces must be non-adhesive and non-penetrating.
• Jacobian matrices: Encode how a constraint violation relates to changes in body velocities.
• Baumgarte stabilization: A technique to add artificial damping to constraints, preventing drift and increasing numerical stability over many time steps.

Constraint solvers are broadly categorized by their solution approach. Iterative solvers, like the Projected Gauss-Seidel (PGS) method, solve constraints sequentially in loops, converging toward a solution. They are fast and memory-efficient for many constraints but may not find an exact solution. Direct solvers (e.g., solving a full LCP matrix) compute an exact mathematical solution in one step but are computationally expensive (O(n³)) and become impractical for large systems. Modern real-time engines like NVIDIA PhysX and Bullet primarily use iterative solvers for their performance, trading some accuracy for speed.

A primary function is resolving collisions. For each contact point generated by collision detection, the solver must calculate impulses (instantaneous forces) to:

• Enforce non-penetration: Apply forces to push intersecting objects apart.
• Model static and dynamic friction: Apply tangential forces that resist sliding motion.
• Handle stacking and resting contacts: Distribute forces correctly through piles of objects without introducing energy (preventing "jitter"). This is often the most computationally intensive part of the physics pipeline, as a single complex collision can generate dozens of contact constraints.

Beyond contact, solvers enforce connections between bodies. Joints (e.g., hinges, sliders, spherical) limit the degrees of freedom (DoF) between connected bodies. The solver calculates the torques needed to maintain these limits. For complex systems like robotic arms, efficient algorithms like the Articulated Body Algorithm (ABA) or Featherstone's algorithm are used to solve the forward dynamics of the entire kinematic chain in O(n) time, making real-time simulation of sophisticated robots possible.

The solver must be robust against numerical errors that accumulate over time. Key techniques include:

• Warm-starting: Using the solution from the previous time step as an initial guess, dramatically reducing the number of iterations needed for convergence.
• Sub-stepping: Advancing the constraint solution multiple times within a single render frame for fast-moving or small objects to prevent tunneling.
• Constraint prioritization: Solving critical constraints (like penetration) before less critical ones (like friction) to improve stability. Performance is measured in iterations per second and is highly dependent on the number of active constraints and the chosen solver tolerance.

The solver does not operate in isolation. It is a central stage in the physics pipeline, which follows a strict order:

1. Broadphase & Narrowphase Collision Detection: Identifies potential and exact contacts.
2. Contact Generation: Creates constraint data from collisions.
3. Constraint Solver Execution: Resolves all contacts and joint constraints.
4. Time Integration: Uses the solved impulses/forces to update body positions and velocities via methods like Symplectic Euler or Runge-Kutta. The solver's output directly drives the forward dynamics of every simulated object, making it the definitive authority on simulated motion.

A Linear Complementarity Problem (LCP) is the core mathematical formulation for modeling contact and friction in rigid body simulations. It provides a framework where the solver must find forces that satisfy two conditions simultaneously: non-penetration (bodies cannot push into each other) and non-adhesion (contact forces can only push, not pull). This is expressed as finding a force vector f such that:

• f ≥ 0 (forces are non-negative, i.e., pushing)
• Jv + b ≥ 0 (relative velocity at contact is separating or zero)
• fᵀ(Jv + b) = 0 (complementary condition: force is only applied if velocity is zero) Where J is the Jacobian matrix and b is a bias term. Solvers like Projected Gauss-Seidel are iterative methods for approximating LCP solutions.

The Projected Gauss-Seidel (PGS) solver is a prevalent iterative algorithm for solving the Linear Complementarity Problems arising from contact and joint constraints. It operates by:

• Sequentially processing each constraint equation.
• Calculating an impulse that would satisfy that single constraint in isolation.
• Projecting (clamping) that impulse to satisfy physical limits (e.g., non-negative for contact, bounded for friction).
• Immediately applying the projected impulse to the involved bodies' velocities, affecting subsequent constraint calculations in the same iteration. This sequential, immediate feedback often leads to faster convergence for contact-dominated scenes compared to batch methods. Its iterative nature means accuracy improves with more iterations, allowing a tunable trade-off between performance and physical correctness.

In physics simulation, a constraint is a mathematical rule that restricts the possible positions or velocities of simulated bodies. Constraints are the problems the solver must resolve. Major types include:

• Contact Constraints: Prevent inter-penetration between colliding bodies. Define that the relative velocity at the contact point must be non-negative (separating or zero).
• Friction Constraints: Model tangential forces at contacts, typically approximated by Coulomb friction cones.
• Joint Constraints: Connect bodies, removing degrees of freedom. Examples:
  • Revolute Joint: Allows only rotation about a single axis (like a door hinge).
  • Prismatic Joint: Allows only translation along a single axis (like a piston).
  • Spherical Joint: Allows three rotational degrees of freedom (like a ball-and-socket).
• Motor Constraints: Actively drive a system to a target position or velocity. The constraint solver's job is to compute the forces or impulses needed to satisfy all active constraints each simulation step.

The Jacobian matrix is a fundamental mathematical construct that linearly relates changes in system state to changes in constraint satisfaction. For a constraint function C(q)=0 (where q are body positions/orientations), the Jacobian J is defined as J = ∂C/∂q. Its critical roles are:

• Velocity Mapping: It projects body velocities into constraint space: the time derivative of C is Ḋ = Jv. The solver aims to make Jv = -b (for stabilization) or Jv ≥ 0 (for contact).
• Force Mapping: It transposes constraint forces from constraint space back into the force/torque space of each body: F = Jᵀλ, where λ is the Lagrange multiplier (the impulse/force calculated by the solver). A single constraint often affects multiple bodies; the Jacobian efficiently encodes these connections. The solver works primarily in the lower-dimensional constraint space defined by J to find λ, which is then applied as forces via Jᵀ.

The physics pipeline is the sequential computational workflow of a physics engine, where the constraint solver is one critical stage. A typical real-time pipeline for a single simulation step (frame) is:

1. Broadphase: Uses spatial structures (like BVHs or grids) to quickly find pairs of objects that are potentially colliding, culling obvious non-colliders.
2. Narrowphase & Contact Generation: Performs precise geometric tests on potential pairs to determine actual collisions, generating contact points, normals, and penetration depths.
3. Constraint Assembly: Converts contacts, joints, and other limits into formal constraint equations with associated Jacobian matrices and limits.
4. Constraint Solving (This Stage): The solver (e.g., PGS) processes all assembled constraints to compute the impulses/forces (λ) that resolve inter-penetration and satisfy joints.
5. Time Integration: Applies the solver's forces along with external forces (like gravity) to update body velocities and positions, using methods like semi-implicit Euler (Symplectic Euler). The solver sits at the heart of this pipeline, ensuring physical plausibility after collisions are detected.

The penalty method is an alternative to constraint-based solvers for enforcing contact and joint limits. Instead of solving for exact constraint forces, it models constraints as very stiff springs and dampers. When a violation occurs (e.g., penetration), the method:

• Applies a spring force proportional to the penetration depth, pushing bodies apart.
• Applies a damping force proportional to the relative velocity at the contact, reducing oscillation. Key Differences from Constraint Solvers:
• Pro: Simpler to implement, often computationally cheaper per contact.
• Con: Forces are approximate and depend on tuning stiffness/damping parameters.
• Con: Can cause "stiff" equations that require small time steps for stability, leading to "springy" or bouncy behavior.
• Con: Does not guarantee exact non-penetration; small overlaps may persist. Constraint solvers (implicit, solving LCPs) are generally preferred for rigid body simulation as they provide exact, parameter-free enforcement at the velocity level, though they are computationally more intensive per iteration.