A Simple Introduction to Position-Based Dynamics (PBD)

1. What is PBD?

The most popular methods for simulating dynamic systems are force-based, such as FEM, MPM, and SPH. In these methods, force is the primary quantity computed to update the system. More specifically, internal forces due to material deformation and external forces due to effects such as gravity and friction are evaluated. Then Newton's second law is applied to compute the acceleration caused by these forces. The velocity is updated using a time integration scheme, and finally the position is updated accordingly. Therefore, the general pathway is deformation ($\mathbf{x}^n$) / loading ($\mathbf{L}^n$) $\rightarrow$ force ($\mathbf{f}$) $\rightarrow$ acceleration ($\mathbf{a}$) $\rightarrow$ velocity ($\mathbf{v}$) $\rightarrow$ position ($\mathbf{P}^{n+1}$), where $n$ denotes the timestep.

The logic is rather natural. Nonetheless, an equally natural question is: could these intermediate quantities be omitted? Something like:

Yes, they can be. This is the idea of PBD, which updates positions directly. How does PBD do this? The key is constraints. This may sound awkward at first, but a constraint is just another word for a rule. Everything in the universe is dictated by rules. If we look more deeply at force-based methods and PBD, they actually share a similar philosophy. In a force-based method, $\mathbf{f}(\mathbf{x}^n,\mathbf{L}^n)$ specifies the force induced by material or element deformation. In PBD, $\mathbf{P}^{n+1}=\mathbf{P}^{n+1}(\mathbf{x}^n,\mathbf{L}^n)$ says that positions must satisfy some rules/constraints. For example, a spring tries to restore its original length despite the applied loading:

where $\mathbf{p}_1$ and $\mathbf{p}_2$ are the positions of the two ends of the spring, respectively. The scalar $d$ is the original length of the spring.

Therefore, the PBD update can be written as:

2. Algorithm

A general PBD algorithm is:

The algorithm above may look scary, but it is very easy to understand.

1. Lines 1--5: Initialize the system, where $m_i$ is the mass of vertex $i$. This is done only once at the beginning of the simulation.
2. Lines 6--25: Update the system at each timestep.
   1. Lines 7--9: Update vertex velocities by considering only the external force $\mathbf{f}_{\mathrm{ext}}$.
   2. Line 10: Add damping to the system. A perfectly elastic system would oscillate forever.
   3. Lines 11--13: Predict vertex positions using the updated velocities.
   4. Lines 14--16: Generate collision constraints if objects collide with each other.
   5. Lines 17--19: This is the core of the PBD system. Compute the updated vertex positions while considering all constraints. Note that multiple constraints can compete with each other, so multiple iterations are required to approximately satisfy them.
   6. Lines 20--24: Update the final vertex velocities and positions for this timestep.

3. Projecting the Constraints

We give a detailed explanation of constraint projection in this section.

Suppose we have $J$ constraints, $C_1,C_2,\dots,C_j,\dots,C_J$. Each constraint can be either an equality constraint, $C_i=0$, or an inequality constraint, $C_j\geq 0$. For an inequality constraint, we project it only when $C_j<0$. Now let us consider an internal equality constraint $C_j=0$, and assume that all vertex masses are equal. Let $\mathbf{p}$ be the concatenated vector of vertex positions. For an internal constraint that depends only on relative shape, not on global translation or rotation, $C$ is invariant under rigid-body modes. Therefore, the gradient of $C_j$ is perpendicular to those rigid-body modes, i.e., translation and rotation. We omit the subscript $j$ for simplicity in the following derivation. Given $\mathbf{p}$, we want to find a correction $\Delta \mathbf{p}$ such that the constraint is satisfied:

Using a first-order Taylor expansion:

If we choose $\Delta \mathbf{p}$ to be in the same direction as $\nabla_{\mathbf{p}} C(\mathbf{p})$, then the linear and angular momentum can be conserved for internal constraints with equal masses. That is:

where $\lambda$ is an unknown scalar.

Substituting this equation into the linearized constraint gives the final correction:

If the vertex masses are not equal, inverse-mass weighting can be applied. See paper [1] for details.

To apply a constraint-specific stiffness $k\in[0,1]$ to the system, we simply multiply it with the correction:

The stiffness parameter determines how much of the correction is applied per iteration.

4. Practical Considerations

4.1 Non-physical Stiffness

The biggest problem with PBD is its non-physical stiffness. Recall that we apply the stiffness $k$ by simply multiplying it with the position correction. If the initial constraint error is $r_0$, then after $n$ repeated projections, the remaining error is approximately

As the number of iterations increases, the system becomes more and more rigid. In other words, the stiffness is iteration-dependent rather than purely physical.

4.2 Gauss-Seidel Solver vs. Jacobi Solver

Given $J$ constraints, Gauss-Seidel processes constraints sequentially and immediately applies each correction. Jacobi computes all constraint corrections from the old positions first, accumulates them per particle, and then applies the combined corrections. Therefore, a Jacobi solver is generally slower to converge than a Gauss-Seidel solver. However, it is well suited for parallelization if conflicting updates are handled properly, for example by simple averaging.

5. Interactive Preview

Here is an interactive demo running a PBD simulation in the browser. You can adjust the stiffness, resolution, time step, and total simulation time, and use the view toggle to inspect the same motion from different perspectives. Press Run to compute the simulation, then use the progress slider to inspect the generated frames.

References

• Muller, M., Heidelberger, B., Hennix, M., and Ratcliff, J. (2007). "Position Based Dynamics." Journal of Visual Communication and Image Representation, 18(2), 109-118.
• Jakobsen, T. (2001). "Advanced Character Physics." GDC.