Position-Based Dynamics

Your browser does not support the HTML5 canvas tag.

Controls

0.00 s 0.00 ms 0 0

PBD algorithm:

Instead of computing forces and integrating accelerations, as in classical force-based simulation, position-based dynamics (PBD) [MHHR06,BMM17] formulates physical relationships as geometric constraints and enforces them by directly displacing particle positions. This example simulates a piece of cloth: a regular grid of particles connected by distance constraints and hanging under gravity from one or more fixed corners. You can grab and drag any of the free particles with the mouse.

Every time step the following steps are performed:

  1. time integration to predict particle positions
  2. loop
  3. compute Lagrange multipliers
  4. determine position correction
  5. update velocities
In this example we use distance constraints $$C_i(\mathbf{x}_{i_1}, \mathbf{x}_{i_2}) = \| \mathbf x_{i_1} -\mathbf x_{i_2} \|-d,$$ where $d$ is the rest length between particles $\mathbf{x}_{i_1}$ and $\mathbf{x}_{i_2}$; a constraint is satisfied once $C_i = 0$. Besides the horizontal and vertical structural constraints along the grid, diagonal constraints are added within each grid cell to resist shearing and keep the cloth from collapsing sideways.

1. Time integration

In the first step the particle positions are advected using a symplectic Euler method to obtain predicted positions $\mathbf x^*$: $$\begin{align*} \mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a_\text{ext}(t) \\ \mathbf x^* &= \mathbf x(t) + \Delta t \mathbf v^*, \end{align*}$$ where $\mathbf a_\text{ext}$ are accelerations due to external forces (here: gravity). Note that this prediction ignores the distance constraints entirely, so $\mathbf x^*$ will generally violate them — correcting this is exactly the job of the solver loop below.

2. Solver loop

The position corrections of all constraints are computed and applied one after another in a single sweep, and this sweep is repeated for a fixed number of iterations (Gauss-Seidel-style relaxation). Because each correction immediately updates the positions seen by the next constraint in the sweep, more iterations propagate the corrections further through the cloth and bring the particles closer to satisfying all constraints simultaneously. Unlike a spring stiffness constant in a force-based model, here it is the number of iterations (relative to the time step) that determines how stretchy or stiff the cloth appears — try reducing it in the controls panel and observe the effect.

3. Compute Lagrange multipliers

The general form to compute the Lagrange multipliers in position-based dynamics is $$\mathbf J \mathbf M^{-1} \mathbf J^T \boldsymbol \lambda = -\mathbf C(\mathbf p),$$ where $\mathbf J = (\partial C / \partial \mathbf x)^T$ is the Jacobian of the constraint function and $\mathbf M$ denotes the mass matrix.

Since the distance constraint is a scalar function the Lagrange multiplier is determined as $$\lambda_i = - \frac{C_i(\mathbf x)}{\sum_j \frac{1}{m_j} \| \partial C_i(\mathbf x) / \partial \mathbf x_j\|^2},$$ where $j$ denotes the indices of the particles which are influenced by the constraint. Fixed particles are excluded from this sum (their inverse mass is treated as zero), which is why they do not move under the correction below.

To compute the Lagrange multipliers, the constraint gradients are required which are computed as:

Constraint gradients:

$$\begin{align*} \frac{\partial C_i}{\partial \mathbf x_{i_1}} &= \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \\ \frac{\partial C_i}{\partial \mathbf x_{i_2}} &= - \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \end{align*}$$

4. Position correction:

The position correction for a particle is determined using the Lagrange multiplier: $$\begin{align*} \Delta \mathbf x_{i_1} &= \frac{1}{m_{i_1}} \frac{\partial C_i}{\partial \mathbf x_{i_1}} \lambda_i = -\frac{\frac{1}{m_{i_1}}}{\frac{1}{m_{i_1}} + \frac{1}{m_{i_2}}}(\|\mathbf{x}_{i_1} - \mathbf{x}_{i_2}\|-d)\frac{\mathbf{x}_{i_1}- \mathbf{x}_{i_2}}{\|\mathbf{x}_{i_1} - \mathbf{x}_{i_2}\|} \\ \Delta \mathbf x_{i_2} &= \frac{1}{m_{i_2}} \frac{\partial C_i}{\partial \mathbf x_{i_2}} \lambda_i = +\frac{\frac{1}{m_{i_1}}}{\frac{1}{m_{i_1}} + \frac{1}{m_{i_2}}}(\|\mathbf{x}_{i_1} - \mathbf{x}_{i_2}\|-d)\frac{\mathbf{x}_{i_1}- \mathbf{x}_{i_2}}{\|\mathbf{x}_{i_1} - \mathbf{x}_{i_2}\|} \end{align*}$$ and the correction is applied to update $\mathbf x^*$ as $$\mathbf x^* := \mathbf x^* + \Delta \mathbf x.$$ Note that heavier particles (larger mass $m$) receive a proportionally smaller correction — this is what makes the constraint respect the mass ratio between the two particles instead of always splitting the correction evenly.

5. Velocity update:

The final positions and velocities are computed as: $$\begin{align*} \mathbf x_i(t+\Delta t) &= \mathbf x_i^* \\ \mathbf v_i(t+\Delta t) &= \frac{1}{\Delta t} (\mathbf x_i(t+\Delta t) - \mathbf x_i(t)) \end{align*} $$ Rather than integrating forces, the velocity is simply reconstructed from how far the position moved during the step. This automatically accounts for all constraint corrections as if they had been caused by impulses, and it introduces a mild, built-in numerical damping (energy is removed whenever a correction shortens a particle's displacement) that contributes to the method's robustness.

Properties

The panel on the right also reports the measured time per simulation step, which grows with the number of particles, constraints and iterations — try increasing the grid resolution or the iteration count to see this trade-off between accuracy and performance directly.

References