EulerImplicitSolver

This component belongs to the category of integration schemes or ODE Solver. This scheme builds the system following an implicit scheme: forces are considered based on the state information at the next time step , unknown at the current time step.

Looking at continuum mechanics, the linear system arises from the dynamic equation. This dynamic is written as follows but other physics (like heat transfer) result in a similar equation:

where is the degrees of freedom, the mass matrix and a function of (and possibly its derivatives) acting on our system. In the case of the EulerImplicitSolver, this equation can be written:

by using a Taylor expansion, we get:

since we have: , then:

Finally, gathering the unknown (depending on ) in the left hand side, we have:

We can notice the appearance of the stiffness matrix : . The stiffness matrix is a symmetric matrix, can either be linear or non-linear regarding .

The computation of the right hand side is done by the ForceFields. Just like in the explicit case (see EulerExplicitSolver), the explicit contribution is implemented in the same function addForce(). The second part is computed by the function addDForce().

It is important to note that, depending on the choice of LinearSolver (direct or iterative), the API functions called to build the left hand side system matrix will not be the same:

• if a direct solver is used, the mass is computed in the addMToMatrix() and the stiffness part is computed in the function addKToMatrix() in ForceFields

• if an iterative solver is used, the mass is iteratively multiplied by the unknown within the addMDx(), as the stiffness part within the function addDForce() in ForceFields.

Considering viscosity

As you might have notice, the Taylor expansion detailed above does not take into account a possible dependency of the force on the velocity. By considering it, the effect of velocity will result in a viscosity effect through the damping matrix .

Let’s apply the Taylor expansion taking into account the velocity and we get:

Depending on the choice of LinearSolver (direct or iterative), the API functions called to build the damping matrix on the left hand side will not be the same:

• if a direct solver is used, the damping matrix is computed in the addBToMatrix() in ForceFields

• if an iterative solver is used, the damping is iteratively multiplied by the unknown within the addDForce() just as the stiffness part in the function addDForce() in ForceFields.

Dissipation

SOFA is a framework aiming at interactive simulations. For this purpose, dissipative schemes are very appropriate. The Euler scheme is an order 1 integration scheme (in time, since only using the current state and no older one like ). It is known to be a dissipative scheme. Moreover, only one Newton step is performed in the EulerImplicit, which might harm the energy conservation.

Activating the trapezoidalScheme option of the Euler implicit scheme will make the scheme less dissipative. This is due to the fact that the trapezoidal rule increases the order of the time integration. Moreover, higher order schemes are known to be less dissipative.

Sequence diagram

Data

The data trapezoidalScheme modifies the EulerImplicitSolver scheme and implements the trapezoidal rule:

This results in the following linear system:

The use of the trapezoidal rule is known to increase robustness and stability to the time integration due to the order 2 in time of this trapezoidal scheme.

The option is given to the user to add numerical Rayleigh damping using the data rayleighStiffness and rayleighMass. The description of the meaning and effect of these Rayleigh damping coefficients is given in ODESolver.

The data firstOrder enables to use the EulerImplicitSolver at the order 1, which means that only the first derivative of the DOFs (state) x appears in the equation. Higher derivatives are absent. This option is for instance well suited for heat diffusion equation using only the first derivative of the temperature field:

.

Usage

The EulerImplicitSolver requires:

• a LinearSolver to solve the linear system
• and a MechanicalObject to store the state vectors.

Example

This component is used as follows in XML format:

<EulerImplicitSolver name="ODEsolver" rayleighStiffness="0.1" rayleighMass="0.1" />

or using SofaPython3:

node.addObject('EulerImplicitSolver', name='ODEsolver', rayleighStiffness='0.1' rayleighMass='0.1')  

An example scene involving a StaticSolver is available in examples/Component/ODESolver/Backward/EulerImplicitSolver.scn