Caribou

Implementation of an implicit backward euler solver compatible with nonlinear materials.
We are trying to solve to following
Where is the mass matrix, is the damping matrix, is the (possibly nonlinear) internal elastic force residual and is the external force vector (for example, gravitation force or surface traction).
We first transform this secondorder differential equation to a first one by introducing two independant variables:
Using the backward Euler scheme, we pose the following approximations:
where is the delta time between the steps and .
Substituting inside gives
And the problem becomes:
where is the vector of unknown accelerations.
Finally, assuming is nonlinear in , we iteratively solve for using the NewtonRaphson method and using the previous approximations to back propagate it inside the velocity and position vectors.
Let be the Newton iteration number for a given time step , we pose
where and are the position and velocity vectors at the beginning of the time step, respectively, and remains constant throughout the Newton iterations.
We then solve for with
And propagate back the new acceleration using and .
In addition, this component implicitly adds a Rayleigh's damping matrix . We therefore have
#include <BackwardEulerODESolver.h>
Public Types  
template<typename T >  
using  Data = sofa::core::objectmodel::Data< T > 
Public Types inherited from SofaCaribou::ode::NewtonRaphsonSolver  
template<typename T >  
using  Data = sofa::core::objectmodel::Data< T > 
template<typename T >  
using  Link = sofa::core::objectmodel::SingleLink< NewtonRaphsonSolver, T, sofa::core::objectmodel::BaseLink::FLAG_STRONGLINK > 
Public Member Functions  
SOFA_CLASS (BackwardEulerODESolver, NewtonRaphsonSolver)  
void  solve (const sofa::core::ExecParams *params, SReal dt, sofa::core::MultiVecCoordId x_id, sofa::core::MultiVecDerivId v_id) override 
Public Member Functions inherited from SofaCaribou::ode::NewtonRaphsonSolver  
SOFA_CLASS (NewtonRaphsonSolver, sofa::core::behavior::OdeSolver)  
void  init () override 
void  solve (const sofa::core::ExecParams *params, SReal dt, sofa::core::MultiVecCoordId x_id, sofa::core::MultiVecDerivId v_id) override 
auto  iteration_times () const > const std::vector< UNSIGNED_INTEGER_TYPE > & 
List of times (in nanoseconds) that each NewtonRaphson iteration took to compute in the last call to Solve().  
auto  squared_residuals () const > const std::vector< FLOATING_POINT_TYPE > & 
The list of squared residual norms (r^2) of every newton iterations of the last solve call.  
auto  squared_initial_residual () const > const FLOATING_POINT_TYPE & 
The initial squared residual (r0^2) of the last solve call.  