move ODESolver to other folder to follow new design

30_Components/27_IntegrationScheme/10_Forward/10_EulerExplicitSolver.md

@@ -0,0 +1,41 @@
+EulerExplicitSolver  
+===================
+
+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)$.
+
+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: 
+
+$$
+\mathbf{M}\Delta v=dt\left(f(x(t))\right)
+$$
+
+since forces only depend on known state (at our current time step). These forces are computed by the ForceField in the `addForce()` function. The system matrix $\mathbf{A}$ is only equal to the mass matrix $\mathbf{M}$.
+
+Depending on whether the mass matrix is diagonal or not, SOFA supports two cases:
+
+1) The mass matrix is diagonal. It makes the resolution of the linear system trivial (best performances). In this case, the system matrix $\mathbf{A}$ equals a diagonal mass matrix $\mathbf{M}$ which is diagonal, and it can be stored as a vector $|m|$ . Moreover, its inverse can directly be obtained as: $\mathbf{M}^{-1}=|m|^{-1}=\frac{1}{|m|}$.
+   The solution $\Delta v_{sol}=dt\cdot \mathbf{M}^{-1}f(x(t))$ finally corresponds to a division operation of $f(x(t))$ by the mass. This computation is actually performed by the Mass component in the `accFromF()` function. Therefore, no LinearSolver is needed to compute directly or iteratively a solution.
+2) The mass matrix is not diagonal. Solving the system requires a linear solver. 
+
+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>
+
+
+
+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.

30_Components/27_IntegrationScheme/20_Backward/20_EulerImplicitSolver.md

@@ -0,0 +1,117 @@
+EulerImplicitSolver  
+===================
+
+This component belongs to the category of [integration schemes or ODE Solver](../../../../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.
+
+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 EulerImplicitSolver, this equation can be written: 
+
+$$
+\mathbf{M} \Delta v=dt \cdot f(x(t+dt))
+$$
+
+by using a Taylor expansion, we get:
+
+$$
+\mathbf{M} \Delta v=dt \cdot \left( f(x(t))+\cdot \frac{\partial f}{\partial x} \Delta x \right)
+$$
+
+since we have: $\Delta x=dt(v(t)+\Delta v)$, then:
+
+$$
+\mathbf{M}\Delta v=dt\cdot \left( f(x(t))+dt\cdot \frac{\partial f}{\partial x}v(t)+dt\cdot \frac{\partial f}{\partial x} \Delta v \right)
+$$
+
+Finally, gathering the unknown (depending on $\Delta v$) in the left-hand side, we have:
+
+$$
+\left( \mathbf{M}-dt^2 \cdot \frac{\partial f}{\partial x} \right) \Delta v=dt\cdot f(x(t))+dt^2\cdot \frac{\partial f}{\partial x}v(t)
+$$
+
+We can notice the appearance of the stiffness matrix : $\mathbf{K}_{ij}=\textstyle\frac{\partial f_i}{\partial x_j}$. The stiffness matrix $\mathbf{K}$ is a symmetric matrix, can either be linear or non-linear regarding $x$.
+
+$$
+\left( \mathbf{M}-dt^2 \cdot \mathbf{K} \right) \Delta v=dt\cdot f(x(t))+dt^2\cdot \mathbf{K}v(t)
+$$
+
+The computation of the **right hand side** is done by the ForceFields. Just like in the explicit case (see [EulerExplicitSolver](../../forward/eulerexplicitsolver/)), the explicit contribution $dt\left(f(x(t))\right)$ is implemented in the same function `addForce()`. The second part $dt^2\cdot \frac{\partial f}{\partial x}v(t)$ 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 $\mathbf{A}=\left( M-dt^2 \cdot \frac{\partial f}{\partial x} \right)$ will not be the same:
+
+  - if a direct solver is used, the mass $\mathbf{M}$ is computed in the `addMToMatrix()` and the stiffness part $-dt^2 \cdot \frac{\partial f}{\partial x}$ is computed in the function `addKToMatrix()` in ForceFields
+
+  - if an iterative solver is used, the mass is iteratively multiplied by the unknown $\mathbf{M} \Delta v$ within the `addMDx()`, as the stiffness part $-dt^2 \cdot \frac{\partial f}{\partial x} \Delta v$ within the function `addDForce()` in ForceFields.
+
+
+#### Considering viscosity
+
+
+As you might have noticed, the Taylor expansion detailed above does not take into account a possible dependency of the force  $f(x,t)$ on the velocity. By considering it, the effect of velocity will result in a viscosity effect through the damping matrix $\mathbf{B}$.
+
+Let's apply the Taylor expansion taking into account the velocity and we get:
+
+$$
+\mathbf{M} \Delta v=dt \cdot \left( f(x(t), v(t))+\cdot \frac{\partial f}{\partial x} \Delta x+\cdot \frac{\partial f}{\partial v} \Delta v \right)
+$$
+
+$$
+\left( \mathbf{M}-dt \cdot \frac{\partial f}{\partial v}-dt^2 \cdot \frac{\partial f}{\partial x} \right) \Delta v=dt\cdot f(x(t),v(t))+dt^2\cdot \frac{\partial f}{\partial x}v(t)
+$$
+
+$$
+\left( \mathbf{M}-dt \cdot \mathbf{B}-dt^2 \cdot \mathbf{K} \right) \Delta v=dt\cdot f(x(t),v(t))+dt^2\cdot \mathbf{K}v(t)
+$$
+
+Depending on the choice of LinearSolver (direct or iterative), the API functions called to build the $\mathbf{B}$ damping matrix on the left hand side will not be the same:
+
+  - if a direct solver is used, the damping matrix $\mathbf{B}$ is computed in the `addBToMatrix()` in ForceFields
+  - if an iterative solver is used, the damping is iteratively multiplied by the unknown $\mathbf{B} \Delta v$ 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 time integration scheme (since only using the current state $x(t)$ and no older one like $x(t-dt)$). It is known to be a dissipative scheme. Moreover, only one Newton step is performed in the EulerImplicit, which might harm the energy conservation.
+
+#### Numerical damping
+
+With Rayleigh damping, the option is given to the user to add numerical damping. The description of the meaning and effect of these Rayleigh damping coefficients is given in [ODESolver](../../../../simulation-principles/system-resolution/integration-scheme/#rayleigh-damping).
+
+
+#### Trapezoidal rule
+
+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](https://en.wikipedia.org/wiki/Trapezoidal_rule) increases the order of the time integration. Moreover, higher order schemes are 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:
+
+$$
+y_{n+1}-y_n=\frac{dt}{2}(f(y_{n+1})+f(y_n))
+$$
+
+This results in the following linear system:
+
+$$
+\left( \mathbf{M}-\frac{dt^2}{4} \frac{\partial f}{\partial x}\right) \Delta v=dt\cdot f(x(t))+\frac{dt^2}{2}\cdot \frac{\partial f}{\partial x}v(t)
+$$
+
+
+Sequence diagram
+----------------
+
+<a href="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/EulerImplicitSolver.png?raw=true"><img src="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/EulerImplicitSolver.png?raw=true" title="Flow diagram for the EulerImplicitSolver"/></a>
+
+
+
+
+Usage  
+-----  
+
+The EulerImplicitSolver **requires**:
+
+- a [LinearSolver](../../../../simulation-principles/system-resolution/linear-solver/) to solve the linear system
+- and a MechanicalObject to store the state vectors.
+

30_Components/27_IntegrationScheme/20_Backward/50_StaticSolver.md

@@ -0,0 +1,41 @@
+StaticSolver  
+============  
+
+This component belongs to the category of [integration schemes or ODE Solver](../../../../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**.
+
+_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:
+
+$x^{n+1}=x^n-\frac{f(x^n)}{f'(x^n)}$ where $f'(x^n) = \frac{df}{dx}(x^n)$
+
+In our elasticity case, the system to solve is $K_i \delta u_i-f=0$. At each iteration of Newton Raphson algorithm $n$ at simulation step $i$, we therefore find:
+
+$$
+\delta u_i^{n+1}=\delta u_i^{n}-\frac{(K_i^n \delta u_i^n-f)}{K_i^n}
+$$
+
+
+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**:
+
+- 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.
+
+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/27_IntegrationScheme/20_Backward/60_NewmarkImplicitSolver.md

@@ -0,0 +1,54 @@
+NewmarkImplicitSolver  
+=====================
+
+This component belongs to the category of [integration schemes or ODE Solver](../../../../simulation-principles/system-resolution/integration-scheme/).  
+
+This scheme is an implicit time integrator for dynamic system using the Newmark scheme. To compute the new position or new velocity, the NewmarkImplicitSolver is based on the following equations:
+
+$$
+x_{t+h}=x_t+h v_t+\frac{h^2}{2}((1-2\beta)a_t+2\beta a_{t+h})
+$$
+
+$$
+v_{t+h}=v_t+h((1-\gamma)a_t+\gamma a_{t+h})
+$$
+
+
+Applied to a mechanical system where $\small Ma_t+(r_MM+r_KK)v_t+Kx_t=f_{ext}$, we need to solve the following system:
+
+
+$$
+\tiny Ma_{t+h}+(r_MM+r_KK)v_{t+h}+Kx_{t+h}=f_{ext}
+$$
+
+$$
+\tiny Ma_{t+h}+(r_MM+r_KK)(v_t+h((1-\gamma)a_t+\gamma a_{t+h}))+K(x_t+hv_t+\frac{h^2}{2}((1-2\beta)a_t+2\beta a_{t+h}))=f_{ext}
+$$
+
+$$
+\tiny (M+h\gamma(r_MM+r_KK)+h^2\beta K)a_{t+h}=f_{ext}-(r_MM+r_KK)(v_t+h(1-\gamma)a_t)-K(x_t+hv_t+\frac{h^2(1-2\beta)}{2}a_t)
+$$
+
+$$
+\tiny ((1+h\gamma r_M)M+(h^2\beta +h\gamma r_K)K)a_{t+h}=f_{ext}-(r_MM+r_KK)v_t-Kx_t-(r_MM+r_KK)(h(1-\gamma)a_t)-K(hv_t+\frac{h^2(1-2\beta)}{2}a_t)
+$$
+
+$$
+\tiny ((1+h\gamma r_M)M+(h^2\beta+h\gamma r_K)K)a_{t+h}=a_t-(r_MM+r_KK)(h(1-\gamma)a_t)-K(hv_t+\frac{h^2(1-2\beta)}{2}a_t)
+$$
+
+
+
+Sequence diagram
+----------------
+
+<a href="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/NewmarkImplicitSolver.png?raw=true"><img src="https://github.com/sofa-framework/doc/blob/master/images/integrationscheme/NewmarkImplicitSolver.png?raw=true" title="Flow diagram for the NewmarkImplicitSolver"/></a>
+
+
+Usage  
+-----  
+
+At each simulation step and each Newton Raphson iteration, the NewmarkImplicitSolver **requires**:
+
+- a [LinearSolver](../../../../simulation-principles/system-resolution/linear-solver/) to solve the linear system
+- and a MechanicalObject to store the state vectors.

30_Components/27_IntegrationScheme/20_Backward/BDFOdeSolver.md

@@ -0,0 +1,70 @@
+BDFOdeSolver
+============
+
+This component belongs to the category of [integration schemes or ODE Solver](../../../../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.
+
+The method relies on [Backward Differentiation Formula](https://en.wikipedia.org/wiki/Backward_differentiation_formula) (BDF).
+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.
+This equation leads to a nonlinear function to solve.
+That is why this component requires a [NewtonRaphsonSolver](NewtonRaphsonSolver.md), which the purpose is to solve nonlinear functions.
+
+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.
+
+Details
+-------
+
+The ODE resulting from Newton's second law of motion is:
+
+$$
+\begin{bmatrix}
+\frac{d q}{d t}\\
+M \frac{d \dot{q}}{d t}
+\end{bmatrix}
+=
+\begin{bmatrix}
+\dot{q}\\
+F(q, \dot{q})
+\end{bmatrix}
+$$
+
+where $q$ and $\dot{q}$ are respectively the position and the velocity, $M$ is the mass matrix, and $F$ is the sum of forces.
+
+We define $y(t)=\begin{bmatrix} q \\ \dot{q} \end{bmatrix}$, and $f(t,y)=\begin{bmatrix} \dot{q} \\ M^{-1} F(q,\dot{q}) \end{bmatrix}$, such that the ODE is $y'=f(t,y)$.
+We are interested in computing $y(t_{n+s})$ and we know the values of $y(t_{n+j})$ for $j \lt s$.
+The Lagrange interpolation polynomial is the linear combination $L(t) = \sum_{j=0}^s y(t_n+j) l_j(t)$, where $l_j$ is the basis polynomials defined as $l_j(t)=\prod_{0 \leq m \leq s, m \neq j} \frac{t-t_{n+m}}{t_{n+j}-t{n+m}}$.
+Then, we approximate $y'$ by $L'$, leading to the equation $\sum_{j=0}^s y_{n+j} l'_j(t_{n+s}) = f(t_{n+s},y_{n+s})$.
+We can now define a nonlinear function $r(q, \dot{q}) = \left[\sum_{j=0}^s y_{n+j} l'_j(t_{n+s})\right] - f(t_{n+s},y_{n+s})$. 
+The roots of this function corresponds to $y_{tn+s}$, i.e. the next values of the state.
+
+To find the root of this function, we use a Newton-Raphson algorithm.
+It means the derivative of the function is necessary.
+At each iteration $i$ of the algorithm, we solve the following linear system:
+
+$$
+\begin{bmatrix}
+l'_s(t_{n+s}) I & -dt I \\
+-dt \frac{\partial F}{\partial q} & l'_s(t_{n+s}) M - dt \frac{\partial F}{\partial \dot{q}}
+\end{bmatrix}
+\begin{bmatrix}
+q^{i+1} - q^i \\
+\dot{q}^{i+1} - \dot{q}^{i}
+\end{bmatrix}
+= -r(q^i, \dot{q}^i)
+$$
+
+This system is solved by a block Gaussian elimination leading to the velocity-based linear system:
+
+$$
+\left(l'_s(t_{n+s}) M - dt \frac{\partial F}{\partial \dot{q}} - \frac{dt^2}{l'_s(t_{n+s})} \frac{\partial F}{\partial q} \right)
+(\dot{q}^{i+1} - \dot{q}^i)
+= -r_{\dot{q}}(q^i, \dot{q}^i) - \frac{dt}{l'_s(t_{n+s})} \frac{\partial F}{\partial q} r_{q}(q^i, \dot{q}^i)
+$$
+
+