eXtended Position-Based Dynamics (XPBD)
|
|||||||||||||||||||||||||||||
XPBD algorithm:This example shows the eXtended Position-Based Dynamics (XPBD) method introduced by Macklin et al. [MMC16,BMM17]. When simulating deformable solids using the original PBD method, we have the problem that the constraint stiffness depends on the time step size and the iteration count. XPBD solves this problem by using a compliance formulation. The method performs the following steps:
1. Time integrationIn 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.2. Solver loopPosition corrections are computed and applied in a loop using a fixed number of iterations.3. Compute Lagrange multipliersPBDThe 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. XPBDWhen using XPBD, we do not directly compute a Lagrange multiplier in each iteration. Instead we start with $\lambda_i = 0$ and then in each iteration we compute a change of the Lagrange multiplier as $$\Delta \lambda_i = \frac{-C_i(\mathbf x) - \tilde \alpha_i \lambda_i}{\sum_j \frac{1}{m_j} \| \partial C_i(\mathbf x) / \partial \mathbf x_j\|^2 + \tilde \alpha_i},$$ where $\tilde \alpha = \frac{\alpha}{\Delta t^2}$ and $\alpha$ is the compliance value which is in our example $\alpha = \frac 1 k$. Note, in this example we set $\alpha = 0$ if $k$ is set to 0 by the user. In this special case the algorithm is equivalent to PBD. Constraint gradients:To compute the Lagrange multipliers, the constraint gradients are required which are computed as: $$\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 change of the Lagrange multiplier: $$\begin{align*} \Delta \mathbf x_{i_1} &= \frac{1}{m_{i_1}} \frac{\partial C_i}{\partial \mathbf x_{i_1}} \Delta \lambda_i \\ \Delta \mathbf x_{i_2} &= \frac{1}{m_{i_2}} \frac{\partial C_i}{\partial \mathbf x_{i_2}} \Delta \lambda_i \end{align*}$$ and the correction is applied to update $\mathbf x^*$ as $$\mathbf x^* := \mathbf x^* + \Delta \mathbf x.$$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*} $$References
|