Constraint Solver
The heart of the PBD algorithm is the constraint solver, which iteratively projects the predicted particle positions to satisfy the defined constraints. Since this projection must be done in a physically plausible manner, ideally conserving the system's total linear and angular momentum for internal constraints, we employ a nonlinear Gauss-Seidel Solver.
Constraint Projection: Momentum Conservation
To keep things as physically plausible as possible we make sure that for any internal constraint, the correction step should not introduce fictitious external forces, often referred to as "ghost forces." This is achieved by ensuring that the net change in linear and angular momentum is zero. Let be the position correction for particle . Linear momentum is conserved if the center of mass remains unchanged by the correction:
Angular momentum is conserved if the net torque produced by the corrections is zero with respect to a common center of rotation :
Remark 32.4.1 (Constraint Gradient). A key insight is that if the correction vector for the concatenated particle positions, , is chosen to be parallel to the constraint gradient , that is: then for many common internal constraints (which are independent of rigid body transformations), both momenta are automatically conserved when particle masses are equal.
As such the correction for each particle will be assumed to be parallel to the constraint gradient. This choice is supported by the fact that since the gradient points in the direction of the steepest ascent of the constraint function, and thus moving in the opposite direction is the most direct way to reduce the constraint error.
Position Correction (General Constraints)
Consider a single constraint involving particles. Given a set of predicted positions that violate this constraint (i.e., ), we seek a correction such that . To make this tractable, we linearize the constraint function using a first-order Taylor expansion around the current positions :
As established above, we restrict the correction to be in the direction of the constraint gradient. This allows us to define the correction for a single particle in terms of a single unknown scalar , which acts as a Lagrange multiplier: where is the inverse mass, ensuring that lighter particles are displaced more. The negative sign is a convention to align with the concept of a repulsive force for a positive constraint violation. In vector form for all involved particles, this is , where is the diagonal mass matrix.
To find the value of , we substitute the correction from Equation (32.4.2) back into the linearized constraint, Equation (32.4.1). For an equality constraint, this becomes: Solving for the Lagrange multiplier yields: With determined, the position correction for each particle is fully defined. We can also write the full correction in a single expression by substituting Equation (32.4.3) into Equation (32.4.2): This equation is the cornerstone of the PBD solver. It provides a direct, computationally efficient method to calculate the position corrections required to satisfy a single linearized constraint while respecting particle masses. This projection is solved multiple times for each constraint within a single time step. It is important to note that for constraints that are linear along their gradient, such as a simple distance constraint, this linearization is exact, and the constraint can be satisfied in a single step.
Example 32.4.1 (Distance Constraint). Let us consider one of the most fundamental constraints: the distance constraint, which enforces a fixed separation between two particles, and . This is a common building block used to model stretching resistance in springs, the edges of a cloth mesh, or rigid links between objects.
The constraint function is defined as the difference between the current distance and the desired rest distance : This function directly measures the error that needs to be corrected. While the general projection method from Equation (32.4.4) can be applied, for this specific and common case, the result simplifies to a very intuitive form. The final position corrections for each particle are given directly by: The term represents the total correction magnitude needed along the vector connecting the particles. This total correction is then distributed between the two particles based on their inverse mass ratio, . The correction is applied along the unit vector , moving the particles towards each other if they are too far apart () and away from each other if they are too close (). This formulation directly satisfies the conservation of linear momentum!
Hierarchical Solver
The standard Gauss-Seidel solver in PBD exhibits slow convergence for low-frequency, large-scale deformations because corrections propagate only locally. To accelerate this, Hierarchical Position Based Dynamics (HPBD) [Müller 2008] introduces a multi-resolution approach. The system is represented as a hierarchy of meshes, from coarse to fine. The solver operates first on the coarsest levels, allowing corrections for large-scale errors to propagate rapidly across the entire object. These corrections are then transferred and refined on successively finer levels in a single coarse-to-fine pass via interpolation (prolongation). This multigrid-inspired technique significantly improves convergence speed. For this to work correctly, constraints on the coarse levels must be unilateral (inequality constraints), resisting only stretching, to avoid artificially restricting large-scale bending and folding of the object.