Mickeal Verschoor and Andrei C. Jalba
ACM Transactions on Graphics, Volume 38, Issue 2, 2019. (pdf)
Abstract: Simulating (elastically) deformable models that can collide with each other and with the environment remains a challenging task. The resulting contact problems can be elegantly approached using Lagrange multipliers to represent the unknown magnitude of the response forces. Typical methods construct and solve a Linear Complementarity Problem (LCP) to obtain the response forces. This requires the inverse of the generalized mass matrix, which is generally hard to obtain for deformable-body problems. In this article, we tackle such contact problems by directly solving the Mixed Linear Complementarity Problem (MLCP) and omitting the construction of an LCP matrix. Since a convex quadratic program with linear constraints is equivalent to an MLCP,we propose to use a Conjugate Residual (CR) solver as the backbone of our collision response system. By dynamically updating the set of active constraints, the MLCP with inequality constraints can be solved efficiently. We also propose a simple yet efficient preconditioner that ensures faster convergence. Finally, our approach is faster than existing methods (at the same accuracy), and it allows accurate treatment of friction.
In the following months I will be working on cleaning the code and provide a number of examples. If you are interested in the code, please send me an email. Once the code is cleaned, I will upload it to my Bitbucket account and will push future updates through that repository.
Overview of the method: the method solves a Mixed Linear Complementarity Problem (MLCP) that is obtained by constraining the motion of the degrees of freedom of the simulation using contact and frictional constraints. Such problems are usually solved by creating and solving the associated Linear Complementarity Problem (LCP). The solution of this problem gives the unknown Lagrangian multipliers for each constraint that is associated with a contact point. Obtaining the LCP from the MLCP in case of deformable body simulations can be challenging, it often requires the inverse of the stiffness matrix of the deformable bodies. Computing this can be a performance bottleneck. In our approach we solve the MLCP directly using a modified version of the Conjugate Residual method. Doing so doesn’t require the inverse (or an approximation) of the stiffness matrix of the simulated system.
The advantage of using a Conjugate Residual method is that the method solves a Least Squares problem. This is potentially interesting in cases when the system becomes over-constrained. Where other numerical methods can fail in such cases, the CR method finds the solution that minimizes the error for all constraints. However, in case a Least Squares solution is found, not all of the constraints do satisfy their complementarity conditions. In our method, we allow the CR method to update the state of each constraint while solving the problem. When a constraint is found that violates its complementarity conditions, the constraint can be activated or deactivated, depending on the situation. This eventually leads to a situation in which a solution is found where none of the constraints are violated, i.e., the solution satisfies all constraints.
Another advantage of our approach is that we can solve the non-linear contact problem without recomputing the LCP matrix. When the method converges, the surface of the colliding entities may have deformed or changed. These changes are not reflected in the current configurations of the constraints. By updating the configuration of all constraints given the current state of the geometry, the method can compute a solution in which the constraints do match the obtained surface. This approach can be seen as a Newton method where each Newton iteration linearizes the contacts and solves the corresponding MLCP. Eventually, the method converges to a solution in which all collisions are perfectly resolved and the distances at contact points with active constraints are less than a specified tolerance. The number of non-linear updates depend on the complexity of the scene. One update of a constraint can potentially trigger other collisions.
Static friction is modeled by Lagrange multipliers. When the friction multipliers exceed the boundary given by the friction coefficient and the contact (normal) multiplier, kinetic friction is applied. Kinetic friction is modeled by continuously approximating the direction of the friction force and magnitude. The direction is kept aligned with the sliding direction, the magnitude is approximated given the current value of the normal multiplier. To improve stability and convergence, both the direction and magnitude approximation use a running average of the last N intermediate sliding directions and normal multiplier magnitudes. Eventually, a solution is obtained that satisfies the conditions for friction.
Preconditioning: The CR method solves two separate problems. 1) the velocity of the DoFs. 2) Lagrangian multipliers of the constraints. Since these problems have different scales, running the CR method on both of them will result in poor convergence. Depending on the problem, CR can prioritize to solve either the velocity or Lagrangian multipliers. Ideally, each step in the CR method should minimize both quantities equally. Using a standard diagonal preconditioner we can solve this problem. With more sophisticated preconditioners we can accelerate the convergence of the method further. The ideal preconditioner is the inverse of the matrix that is being solved, therefore we approximate the inverse of the system matrix by an approximation of its block-wise inverse. This approach accelerates the convergence further.
Additional optimizations: Every time a constraint changes its state (from active/inactive, or kinetic/static friction), the residual needs to be recomputed. This step requires a multiplication with the full system matrix of the simulation. This is usually an expensive operation that we ideally want to avoid. By splitting the residual vector in a pure unconstrained residual and a residual associated with the constraints, we can omit this multiplication. Only the residual of the constraints require to be recomputed. By adding the unconstrained residual and the constraint residual together, the actual residual vector is obtained.
Update May 2020: recently I noticed some very subtle bugs in the collision detection functions after cleaning the code. This resulted very sporadically in a missed collision. Due to this, some of the test cases failed to run properly. At the moment I’m running a large number of test cases to ensure that the collision detection runs without any issues. A missed collision can result in wrongly performed inside/outside tests. Based on this information constraints are initialized. If the inside/outside information is incorrect, also a constraint is incorrectly initialized, which could potentially result in counteracting constraints. With this now hopefully fixed, I’ll expect to upload the code soon.