ImpulseBased Dynamic Simulation


IBDS algorithm:This example shows the impulsebased dynamic simulation (IBDS) method [BFS05, BS06]. IBDS is a predictorcorrector method to perform a simulation step for a constrained system. For each position constraint $C(\mathbf x) = 0$ a corresponding velocity constraint is defined as the time derivative $\dot C(\mathbf x) = 0$.
1. Position constraint solver loopCorrection impulses are computed and applied in a loop using a fixed number of iterations or until the error is within a userdefined tolerance.2. Time integrationIn the first step of the position constraint solver 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.3. Compute Lagrange multipliersThe general form to compute the Lagrange multipliers in IBDS is $$\mathbf J \mathbf M^{1} \mathbf J^T \boldsymbol \lambda = \frac{1}{\Delta t} \mathbf C(\mathbf x^*),$$ 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)}{\Delta t \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. To compute the Lagrange multipliers, the constraint gradients of the distance constraints 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. Correction impulse:The correction impulse is determined as: $$\mathbf p = \mathbf J^T \lambda$$5. Velocity update:Finally, the velocity of a constrained particle is updated: $$\Delta \mathbf v = \frac{1}{m} \mathbf p$$Velocity constraint solver loopThe velocity constraints $\dot C(\mathbf x) = 0$ are solved analogously to the position constraints. The Lagrange multipliers are determined by solving: $$\mathbf J \mathbf M^{1} \mathbf J^T \boldsymbol \lambda = \mathbf {\dot C}(\mathbf x).$$ However, note that no position prediction is required since the solve is performed after the time integration step.Time integrationThe particles are advected by numerical time integration. In our case we use a symplectic Euler method: $$\begin{align*} \mathbf v(t + \Delta t) &= \mathbf v(t) + \frac{\Delta t}{m} \mathbf f^{\text{ext}}(t) \\ \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t), \end{align*}$$ where $\mathbf f^{\text{ext}}$ are the external forces.References
