New integration schemes - new doc !

20_Simulation_Principles/20_MechanicalObject.md

@@ -122,7 +122,7 @@ Symbolic ids
 The MultiVec entries are not directly accessible by the solvers. The MultiVec are represented by identificators. The operations on the vectors are implemented using visitors which contain the identificators of the relevant vectors. The MultiVec identificators (MultiVecId) have different types, depending on the data they contain
 (positions or their derivatives) and the access mode.
 
-The use of symbolic identificators (MultiVecId) prevent other components (like solvers) from handling state vectors directly and allow to easily work with abstract MultiVec by using their ids. These symbolic ids are widely used by specialized visitors, like the ones used in [ODESolver](./../system-resolution/integration-scheme/).
+The use of symbolic identificators (MultiVecId) prevent other components (like solvers) from handling state vectors directly and allow to easily work with abstract MultiVec by using their ids. These symbolic ids are widely used by specialized visitors, like the ones used in [IntegrationScheme](./../system-resolution/integration-scheme/).
 
 ``` cpp
 typedef TMultiVecId<V_COORD, V_READ>  ConstMultiVecCoordId;

30_Components/40_ODESolver/10_Forward/10_EulerExplicitSolver.md → 30_Components/40_IntegrationScheme/10_Forward/10_EulerExplicitIntegrationScheme.md

@@ -1,15 +1,15 @@
-EulerExplicitSolver  
+EulerExplicitIntegrationScheme  
 ===================
 
-The EulerExplicitSolver component belongs to the category of [integration schemes or ODE Solver](../../../../simulation-principles/system-resolution/integration-scheme/). This scheme allows to solve dynamic systems explicitly: all forces will be computed based on the state information at the current time step $x(t)$.
+The EulerExplicitIntegrationScheme component belongs to the category of [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/). This scheme allows to solve dynamic systems explicitly: all forces will be computed based on the state information at the current time step $x(t)$.
 
 Looking at continuum mechanics, the linear system $\mathbf{A}x=b$ arises from the dynamic equation. This dynamic is written as follows but other physics (like heat transfer) result in a similar equation:
 
 $$
 \mathbf{M}\Delta v=dt\left(f(x,t)\right)
 $$
 
-where $x$ is the degrees of freedom, $\mathbf{M}$ the mass matrix and $f(x,t)$ a function of $x$ (and possibly its derivatives) acting on our system. In the case of the EulerExplicitSolver, this equation can be written: 
+where $x$ is the degrees of freedom, $\mathbf{M}$ the mass matrix and $f(x,t)$ a function of $x$ (and possibly its derivatives) acting on our system. In the case of the EulerExplicitIntegrationScheme, this equation can be written: 
 
 $$
 \mathbf{M}\Delta v=dt\left(f(x(t))\right)
@@ -26,16 +26,10 @@ Depending on whether the mass matrix is diagonal or not, SOFA supports two cases
 Note that the **symplectic** data allows to modify the scheme to make it [symplectic](https://en.wikipedia.org/wiki/Semi-implicit_Euler_method), i.e. velocities are updated before the positions.
 It allows to update the positions from the newly computed velocities, instead of velocities from the previous time step.
 This option makes the scheme more robust in time.
-EulerExplicitSolver is symplectic by default.
-
-Sequence diagram
-----------------
-
-<a href="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/EulerExplicitSolver.png?raw=true"><img src="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/EulerExplicitSolver.png?raw=true" title="Flow diagram for a EulerExplicitSolver"/></a>
-
+EulerExplicitIntegrationScheme is symplectic by default.
 
 
 Usage  
 -----  
 
-The EulerExplicitSolver **requires** a MechanicalObject to store the state vectors. However, as explained above, no LinearSolver is needed and the EulerExplicitSolver is **only working using a [UniformMass](../../../mass/uniformmass/) or [DiagonalMass](../../../mass/diagonalmass/)**, which ensures to have a diagonal system matrix.
+The EulerExplicitIntegrationScheme **requires** a MechanicalObject to store the state vectors. However, as explained above, no LinearSolver is needed and the EulerExplicitIntegrationScheme is **only working using a [UniformMass](../../../mass/uniformmass/) or [DiagonalMass](../../../mass/diagonalmass/)**, which ensures to have a diagonal system matrix.

30_Components/40_IntegrationScheme/20_Backward/20_EulerImplicitIntegrationScheme.md

@@ -0,0 +1,67 @@
+EulerImplicitIntegrationScheme  
+===================
+
+This component belongs to the category of [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/). This scheme builds the system following an implicit scheme: forces are considered based on the state information at the next time step $x(t+dt)$, unknown at the current time step.
+
+This is the most broadly used integration scheme thank's to its unconditional stability and simplicity.  
+
+It is recommended to read the theoretical part of the [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/) page to understand the following. 
+
+The EulerImplicitIntegrationScheme inherits from VelocityBaseIntegrationScheme because its integration scheme in velocity is invertible. Its equations are the following : 
+
+$$
+\begin{aligned}
+g_{\boldsymbol{x}}^{(t,h)} &: \boldsymbol{v}, \boldsymbol{a} \mapsto \boldsymbol{x}_t + h \boldsymbol{v} \\
+g_{\boldsymbol{v}}^{(t,h)} &: \boldsymbol{a} \mapsto \boldsymbol{v}_t + h \boldsymbol{a}
+\end{aligned}
+$$
+
+
+#### API specialization
+
+As explained in the [IntegrationScheme](../../../../simulation-principles/system-resolution/integration-scheme/) documentation, the specilization consist in implementing 4 methods that requires the knowledge of four terms/expressions. for Euler implicit they are the following : 
+$$
+\begin{aligned}
+&\tilde{g}_{\boldsymbol{x}}^{(t,h)} \equiv g_{\boldsymbol{x}}^{(t,h)} \\
+&g_{\boldsymbol{v}}^{(t,h)-1} : \boldsymbol{v} \mapsto \frac{1}{h}(\boldsymbol{v} - \boldsymbol{v}_t)
+\end{aligned}
+\qquad \qquad \qquad 
+\begin{aligned}
+&\frac{\mathrm{d} \tilde{g}_{\boldsymbol{x}}^{(t,h)}}{\mathrm{d} \boldsymbol{v}} = h \\
+&\frac{\mathrm{d} g_{\boldsymbol{v}}^{(t,h)-1}}{\mathrm{d} \boldsymbol{v}} = \frac{1}{h}
+\end{aligned}
+$$
+
+
+#### Trapezoidal rule
+
+Activating the trapezoidalScheme option of the Euler implicit scheme will make the scheme less dissipative. This will apply the [trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule) to the velocity integration, and thus will increase the order of the time integration, which is known to be less dissipative.
+It is also known to increase robustness and stability to the time integration due to the order 2 in time of this trapezoidal scheme. The modified scheme is the following:
+
+$$
+\begin{aligned}
+g_{\boldsymbol{x}}^{(t,h)} &: \boldsymbol{v}, \boldsymbol{a} \mapsto \boldsymbol{x}_t + \frac{h}{2} \boldsymbol{v_t} + \frac{h}{2} \boldsymbol{v} \\
+g_{\boldsymbol{v}}^{(t,h)} &: \boldsymbol{a} \mapsto \boldsymbol{v}_t + h \boldsymbol{a}
+\end{aligned}
+$$
+
+This results in the updated term: 
+
+$$
+
+\frac{\mathrm{d} \tilde{g}_{\boldsymbol{x}}^{(t,h)}}{\mathrm{d} \boldsymbol{v}} = h/2 
+$$
+
+To activate this trapezoidal rule, you need to use the data `trapezoidalScheme=true`.
+
+
+
+
+Usage  
+-----  
+
+The EulerImplicitIntegrationScheme **requires**:
+
+- a [LinearSolver](../../../../simulation-principles/system-resolution/linear-solver/) to solve the linear system
+- and a MechanicalObject to store the state vectors.
+

30_Components/40_ODESolver/20_Backward/50_StaticSolver.md → 30_Components/40_IntegrationScheme/20_Backward/50_StaticEquilibriumIntegrationScheme.md

@@ -1,17 +1,17 @@
-StaticSolver  
+StaticEquilibriumIntegrationScheme  
 ============  
 
-This component belongs to the category of [integration schemes or ODE Solver](../../../../simulation-principles/system-resolution/integration-scheme/).  
+This component belongs to the category of [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/).  
 
 In the field of mechanics, statics consists in finding the equilibrium taking into account the loads (internal forces, external forces and torques) acting on the physical system, that do not experience an acceleration ( $a=0$ ). Finding a static equilibrium means finding a solution to: $\textstyle \sum F=0$ where $F$ is the sum of all loads, one of which might be unknown.  
 
 In a static analysis, the inertia and damping effects are ignored, i.e. the dynamic effect of the mass is ignored. It can thus be written: $M=I \alpha=0$. In the same way, when running a static simulation, time is not elapsing and time steps should rather be considered as convergence steps.  
 
 In a static simulation involving elasticity, the linear system that we solve corresponds to $K \Delta u=f$ where $K$ is the stiffness matrix (derivative of elastic forces), $\Delta u$ is a vector describing the total increment of displacement and $f$ are all explicit forces. We realize here that the static solver is in fact an implicit scheme, since the $K$ matrix is present in the left-hand side of the equation. The solution $\Delta u$ is obtained iteratively. At each iteration _i_, the displacement is incremented $\Delta u_{i+1}=\Delta u_{i}+\delta u_i$, thus resulting in the following system to solve: $K_i \delta u_i=f$.  
 
-In case of non-linear elasticity, $K_i$ is a linearization which must be updated with regards to the increment of displacement $\delta u_i$. In such cases, several iterations of Newton Raphson are required to find an appropriate approximate solution. In one step of the StaticSolver, the number of Newton Raphson iterations is ruled by the data field **newton_iterations**.
+In case of non-linear elasticity, $K_i$ is a linearization which must be updated with regards to the increment of displacement $\delta u_i$. In such cases, several iterations of Newton Raphson are required to find an appropriate approximate solution. In one step of the StaticEquilibriumIntegrationScheme, the number of Newton Raphson iterations is ruled by the data field **newton_iterations**.
 
-_Reminder_: the [Newton Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method) is an iterative algorithm aiming at finding the solution of the system $f(x)=0$ where $f(x)$ is non-linear. At each iteration of Newton Raphson algorithm, we find a new approximate solution:
+_Reminder_: the Newton Raphson method is an iterative algorithm aiming at finding the solution of the system $f(x)=0$ where $f(x)$ is non-linear. At each iteration of Newton Raphson algorithm, we find a new approximate solution:
 
 $x^{n+1}=x^n-\frac{f(x^n)}{f'(x^n)}$ where $f'(x^n) = \frac{df}{dx}(x^n)$
 
@@ -22,20 +22,14 @@ $$
 $$
 
 
-Sequence diagram
-----------------
-
-<a href="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/StaticSolver.png?raw=true"><img src="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/StaticSolver.png?raw=true" title="Flow diagram for the StaticSolver"/></a>
-
-
 Usage  
 -----  
 
-At each simulation step and each Newton Raphson iteration, the StaticSolver **requires**:
+At each simulation step and each Newton Raphson iteration, the StaticEquilibriumIntegrationScheme **requires**:
 
 - a [LinearSolver](../../../../simulation-principles/system-resolution/linear-solver/) to solve the linear system
 - and a MechanicalObject to store the state vectors.
 
-A StaticSolver must be used in simulations where the dynamics has no or a negligible effect on the system. A StaticSolver would also be relevant for systems with low mass. In such case, we fall into the quasi-static analysis.
+A StaticEquilibriumIntegrationScheme must be used in simulations where the dynamics has no or a negligible effect on the system. A StaticEquilibriumIntegrationScheme would also be relevant for systems with low mass. In such case, we fall into the quasi-static analysis.
 
 In some loading configuration, applying the full forces and torques might not lead to any converging simulation. It is then relevant to go for an incremental loading, i.e. loads are applied incrementally at each simulation step  $i$. This incremental loading has to be done in the associated ForceField. If you want to use this solver with Newton Raphson iterations, it is in the user's hand to make sure the external forces used in the scene (pressure, traction, etc.) only get incremented at each time step, and not at each calls to addForce (which is currently the case for most force fields).

30_Components/40_IntegrationScheme/20_Backward/60_NewmarkImplicitIntegrationScheme.md

@@ -0,0 +1,43 @@
+NewmarkIntegrationScheme  
+=====================
+
+This component belongs to the category of [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/). This scheme builds the system following an implicit scheme: forces are considered based on the state information at the next time step $x(t+dt)$, unknown at the current time step.
+
+It is recommended to read the theoretical part of the [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/) page to understand the following. 
+
+The Newmark-$\beta$ integration scheme is a parametrized integration scheme. It presents a set of two parameters $\zeta = (\beta, \gamma)$, balancing the 'implicitness' of the solver:
+1. With $\zeta =(0,0.5)$, the scheme is an explicit central difference scheme
+2. With $\zeta =(0.25,0.5)$ it behaves like an average constant acceleration scheme
+3. With $\zeta =(1/6,0.5)$ it behaves like a linear accelerations scheme
+4. If $2 \beta \geq \gamma \geq 1/2$, then the Newmark-$\beta$ method is stable regardless of the size of the time-step
+
+It inherits from AccelerationBAsedIntegrationScheme as only the acceleration integration is implicit. Its equations are the following : 
+
+$$
+\begin{aligned}
+g_{\boldsymbol{x}}^{(t,h)} &: \boldsymbol{v}, \boldsymbol{a} \mapsto \boldsymbol{x}_t + h \boldsymbol{v}_t + h^2\left[\left(\frac{1}{2} - \beta\right)\boldsymbol{a}_t + \beta \boldsymbol{a}\right] \\
+g_{\boldsymbol{v}}^{(t,h)} &: \boldsymbol{a} \mapsto \boldsymbol{v}_t + h \left[(1-\gamma)\boldsymbol{a}_t + \gamma \boldsymbol{a}\right]
+\end{aligned}
+$$ 
+
+
+#### API specialization
+
+As explained in the [IntegrationScheme](../../../../simulation-principles/system-resolution/integration-scheme/) documentation, the specialization consist in implementing 5 methods that requires the knowledge of four terms/expressions. For Newmark-$\beta$ they are the following : 
+$$
+\frac{\mathrm{d} g_{\boldsymbol{x}}^{(t,h)}}{\mathrm{d} \boldsymbol{v}} = 0 
+\qquad \qquad \qquad 
+\frac{\mathrm{d} g_{\boldsymbol{x}}^{(t,h)}}{\mathrm{d} \boldsymbol{a}} = h^2\beta 
+\qquad \qquad \qquad 
+\frac{\mathrm{d} g_{\boldsymbol{v}}^{(t,h)}}{\mathrm{d} \boldsymbol{a}} = h \gamma 
+$$
+
+
+
+Usage  
+-----  
+
+At each simulation step and each Newton Raphson iteration, the NewmarkIntegrationScheme **requires**:
+
+- a [LinearSolver](../../../../simulation-principles/system-resolution/linear-solver/) to solve the linear system
+- and a MechanicalObject to store the state vectors.

30_Components/40_IntegrationScheme/20_Backward/BDFIntegrationScheme.md

@@ -0,0 +1,69 @@
+BDFIntegrationScheme
+============
+
+This component belongs to the category of [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/).
+It is an implicit method for the numerical integration of the ODE resulting from Newton's second law of motion.
+
+It is recommended to read the theoretical part of the [integration schemes](../../../../simulation-principles/system-resolution/integration-scheme/) page to understand the following. 
+
+The method relies on [Backward Differentiation Formula](https://en.wikipedia.org/wiki/Backward_differentiation_formula) (BDF). This is a special case of [linear multistep method](https://en.wikipedia.org/wiki/Linear_multistep_method).  
+To integrate the ODE in time, it uses information from the previous time steps to compute the next step.
+It establishes a linear combination of the unknown states, the previous states and the values of the ODE function when applied on those states.
+
+It the specific case of BDF, the coefficients of the linear combination come from the approximation of the function by a Lagrange interpolation polynomial.
+The order of the BDF is the number of previous time steps required to approximate the interpolation polynomial.
+The first-order BDF requires a single time step in the past to compute the next.
+It corresponds to the [backward Euler method](EulerImplicitSolver.md).
+The coefficients are unique for a given order, but can be influenced by a change of time step size.
+The SOFA component supports any order, and any change of time step size.
+
+In the following we are going to present the generic linear multistep method equations from which BFD derives. It has to be noted that we made the choice to only apply the linear multi step method to the position integration, not on the velocity. Using it to both decreases the stability margin of the integration scheme, leading to big instabilities in the scenes. 
+
+The LinearMultistepIntegrationScheme inherits from VelocityBaseIntegrationScheme because its integration scheme in velocity is invertible. Its equations are the following : 
+
+$$
+\begin{aligned}
+g_{\boldsymbol{x}}^{(t,h)} &: \boldsymbol{v}, \boldsymbol{a} \mapsto - \sum_{j=1}^{n} \frac{\alpha_j}{\alpha_{n+1}} \boldsymbol{x}_{t-n+j} + h \frac{\beta_{n+1}}{\alpha_{n+1}} \boldsymbol{v} + h \sum_{j=1}^{n} \frac{\beta_j}{\alpha_{n+1}} \boldsymbol{v}_{t-n+j} \\
+g_{\boldsymbol{v}}^{(t,h)} &: \boldsymbol{a} \mapsto \boldsymbol{v}_{t} + h \boldsymbol{a}
+\end{aligned}
+$$
+
+with $n$ the order of the integration scheme
+
+#### API specialization
+
+As explained in the [IntegrationScheme](../../../../simulation-principles/system-resolution/integration-scheme/) documentation, the specilization consist in implementing 4 methods that requires the knowledge of four terms/expressions. for Euler implicit they are the following : 
+$$
+\begin{aligned}
+&\tilde{g}_{\boldsymbol{x}}^{(t,h)} \equiv g_{\boldsymbol{x}}^{(t,h)} \\
+&g_{\boldsymbol{v}}^{(t,h)-1} : \boldsymbol{v} \mapsto \frac{\boldsymbol{v} - \boldsymbol{v}_t}{h}
+\end{aligned}
+\qquad \qquad \qquad 
+\begin{aligned}
+&\frac{\mathrm{d} \tilde{g}_{\boldsymbol{x}}^{(t,h)}}{\mathrm{d} \boldsymbol{v}} = h \frac{\beta_{n+1}}{\alpha_{n+1}} \\
+&\frac{\mathrm{d} g_{\boldsymbol{v}}^{(t,h)-1}}{\mathrm{d} \boldsymbol{v}} = \frac{1}{h}
+\end{aligned}
+$$
+
+
+#### BDF implementation
+
+The generic implementation of LinearMultistepIntegrationScheme requires to overide only one virtual funciton being :
+
+```cpp
+// Method that will compute the $\alpha$ and $\beta factors$
+virtual void computeFactors() = 0;
+```
+
+The the implementation of BDFIntegrationScheme only computes those terms using the previously cited Lagrange polynomial interpolation. 
+
+
+-----  
+
+The BDFIntegrationScheme **requires**:
+
+- a [LinearSolver](../../../../simulation-principles/system-resolution/linear-solver/) to solve the linear system
+- and a MechanicalObject to store the state vectors.
+
+
+