Physics-Based Simulation

Penalty Method

At the beginning of each time step towards time $n+1$, we evaluate nodal position ${\hat{\mathbf{x}}}_k^{n+1}$ for each BC node $k$ based on their prescribed motions. During each Newton iteration $i$, for the iterate $x^i$, we define a velocity residual to assess how close each BC node is to meeting its target: $$r_{\text{BC},k}^i=\frac{1}{h}\Vert {\mathbf{x}}_k^i-{\hat{\mathbf{x}}}_k^{n+1}\Vert .$$ When $r_{\text{BC},k}^i$ falls below a specific tolerance $ϵ$ for any BC node $k$, we can fix the node at its current location ${\mathbf{x}}_k^i\approx {\hat{\mathbf{x}}}_k^{n+1}$ and apply the DOF elimination method in the subsequent iterations. This is particularly straightforward in scenes with only static BCs, where the DOF elimination method is directly applied.

For other BC nodes $k$ that are far from their target locations, we introduce new penalty terms to the Incremental Potential for each of these nodes: $$\frac{{\kappa }_{\text{M}}}{2}{m}_k\Vert {\mathbf{x}}_k-{\hat{\mathbf{x}}}_k^{n+1}{\Vert }^2.\text{\hspace{1em}(11.1.1)}$$ Here, $m_k$ represents the nodal mass, allowing for intuitive setting of the penalty stiffness ${\kappa }_{\text{M}}$, as the Hessian of the penalty term with respect to BC nodes is simply ${\kappa }_{\text{M}}$ times that of the inertia term.

Remark 11.1.1. For collision obstacles (CO), precisely calculating node masses is challenging due to unknown factors like density. A practical approach is to assume a density similar to that of the simulated solids in the scene. This assumption makes the diagonal entries on the Hessian of the penalty terms roughly ${\kappa }_{\text{M}}$ times that of the inertia term.

For codimensional COs such as shells, rods, and particles, the key is to consider a reasonably large thickness when calculating their volumes. This helps in ensuring that their physical properties align more closely with those of the main simulation bodies.

Setting the penalty stiffness ${\kappa }_{\text{M}}$ appropriately can be challenging. If ${\kappa }_{\text{M}}$ is set too low, it may not effectively move the BC node towards its target. Conversely, a too high ${\kappa }_{\text{M}}$ can lead to numerical issues. Thus, we initially set ${\kappa }_{\text{M}}$ to a reasonably large value and adaptively increase it as necessary.

During the Newton solve, if there are BC nodes $k$ where $r_{\text{BC},k}^i\ge ϵ$ at the point of Newton convergence, we double the penalty stiffness ${\kappa }_{\text{M}}$ to $2\times$ its current value and continue the Newton solve. This process is repeated until all BCs are satisfactorily met at convergence.

Remark 11.1.2. In practice, with double precision floating-point numbers, initializing ${\kappa }_{\text{M}}$ below $10^6$ is typically sufficient, given that the Hessian of the stiff penalty terms is purely diagonal. However, if certain BCs remain unsatisfied even when ${\kappa }_{\text{M}}$ is increased to above $10^{10}$, the optimization process may stall due to severe numerical errors. This stalling occurs because extremely stiff penalty terms are in conflict with the contact barriers. However, such a scenario would likely only occur under a rare CO/BC setting in a manner far more extreme than what is tested in Figure 2.3.1.