BackwardEulerODESolver.h

#pragma once
 
#include <SofaCaribou/Ode/NewtonRaphsonSolver.h>
#include <SofaBaseLinearSolver/DefaultMultiMatrixAccessor.h>
#include <sofa/core/objectmodel/Data.h>
 
namespace SofaCaribou::ode {
 
class BackwardEulerODESolver : public NewtonRaphsonSolver {
public:
    SOFA_CLASS(BackwardEulerODESolver, NewtonRaphsonSolver);
 
    template <typename T>
    using Data = sofa::core::objectmodel::Data<T>;
 
    BackwardEulerODESolver();
 
    void solve (const sofa::core::ExecParams* params, SReal dt, sofa::core::MultiVecCoordId x_id, sofa::core::MultiVecDerivId v_id) override;
private:
 
    void assemble_rhs_vector(const sofa::core::MechanicalParams & mechanical_parameters,
                             const sofa::core::behavior::MultiMatrixAccessor & matrix_accessor,
                             sofa::core::MultiVecDerivId & f_id,
                             sofa::defaulttype::BaseVector * f) final;
 
    void assemble_system_matrix(const sofa::core::MechanicalParams & mechanical_parameters,
                                sofa::component::linearsolver::DefaultMultiMatrixAccessor & matrix_accessor,
                                sofa::defaulttype::BaseMatrix * A) final;
 
    void propagate_position_increment(const sofa::core::MechanicalParams & mechanical_parameters,
                                      const sofa::core::behavior::MultiMatrixAccessor & matrix_accessor,
                                      const sofa::defaulttype::BaseVector * dx,
                                      sofa::core::MultiVecCoordId & x_id,
                                      sofa::core::MultiVecDerivId & v_id,
                                      sofa::core::MultiVecDerivId & dx_id) final;
 
    Data<double> d_rayleigh_stiffness;
    Data<double> d_rayleigh_mass;
 
 
    sofa::core::MultiVecCoordId p_previous_x_id;
 
    sofa::core::MultiVecDerivId p_previous_v_id;
 
    sofa::core::MultiVecDerivId p_a_id;
};
 
} // namespace SofaCaribou::ode