All dynamic simulations assume to discretize the temporal evolution of the system through small time steps. This time step is usually noted dt. An integration scheme is the numerical method describing how to find the solution of ordinary differential equations (ODE) at the next time step. They are usually called ODESolver in SOFA.
Let’s write our ordinary differential equation as follows:
Two main categories of integration schemes exist:
explicit scheme: means that the new time step (t + dt) is computed based on information of the previous time step (t). For instance, in mechanics, internal or external forces would be computed on previous positions (x(t)). The ordinary differential equation looks like: Explicit schemes are usually known as being fast to solve (since the created linear system is lighter) but they require very small time steps, unless they may undergo stability issues.
implicit scheme: means that the new time step (t + dt) is computed based on information of this next time step (t + dt). For instance, in mechanics, internal or external forces would be computed on unknow positions at the next time step (x(t + dt)). The ordinary differential equation looks like: Implicit schemes are known as being slower to solve (the outcoming linear system is more complex) but they are more way more stable than explicit schemes.
Semi-implicit or semi-explicit are actually a combination of two explicit and implicit schemes together for different parts of the ordinary differential equation.
In the SOFA code
The integration scheme is described in the solve() function of the ODESolver. This solve() function will build the complete linear system The left hand side matrix is built using the function:
matrix = MechanicalMatrix(a,b,c);
with the matrix equals where is the mass matrix, is the damping matrix and is the stiffness matrix. The right hand side vector b is built through the function:
Previous positions, new positions, forces, etc. are MultiVec, i.e. vectors stored in the MechanicalState (MechanicalObject).
This matrix system is then sent to a LinearSolver in charge of finally solving the system according to the chosen scheme. Within the function ODESolver::solve(), the call to the LinearSolver will appear through the function call:
The Rayleigh damping is a numerical damping. This damping has therefore no physical meaning and must not be mixed up with physical damping (like DiagonalVelocityDampingForceField in SOFA). The Rayleigh damping corresponds to a damping matrix that is proportional to the mass or/and stiffness matrices using ceofficients, respectively Rayleigh stiffness factor or Rayleigh mass factor . This numerical damping is usually used to stabilize or ease convergence of the simulation. However, it has to be used carefully.
You can see the use of Rayleigh mass and stifness in the solve() function of the EulerImplicit class (see EulerImplicitSolver.cpp).
Last modified: 5 October 2018