Static system

This section details the formulation that is used to solve constraints in SOFA. During each step i of the simulation, a linearization of the internal forces is computed:

f(x_i) ≈ f(x_i−1) + K(x_i−1)dx              (1)

where f provides the volumetric internal stiffness forces at a given position x of the nodes, K(x) is the tangent stiffness matrix that depends on the actual position of the nodes and dx is the difference between consecutive positions in time dx = x_i − x_i−1. The lines and columns that correspond to fixed nodes are removed from the system to improve the condition number of the matrix K. Static equilibrium (the sum of external and internal forces equal to zero) is sought at each step:

−K(x_i−1)dx = p + f(x_i−1) + J^T λ        (2)

where p represents the external forces (e.g. gravity) that are known and J^T λ gathers the contributions of the Lagrange multipliers. The next step consists in the projection of the FEM model equations into the constraint space: the size of matrix K is often very large so an optimization in the motion space would be computationally very expensive. Instead, using the Schur complement of the constraint problem, we do a projection that drastically reduces the size of the search space. Thus three steps are performed. Step I: a free configuration x^free of the deformable model is found by solving equation (2) with λ = 0. We compute the constraint violation noted δ^free.Step II: This step is central in the method. It consists in projecting the mechanics into the constraint space. We obtain the smallest possible projection space:

δ = JK⁻¹J^Tλ  + δ^free         (3)

Step III: The final configuration, at the end of the time step, is corrected by using the value of the constraint response:

x = x^free + K⁻¹J^Tλ

Bibliography

• Nouvelles architectures parallèles pour simulations interactives médicales
• Real-time haptic simulation of medical procedures involving deformations and device-tissue interactions