NewtonRaphsonSolver.h

#pragma once
 
#include <Caribou/config.h>
#include <sofa/core/behavior/OdeSolver.h>
#include <sofa/core/behavior/LinearSolver.h>
#include <sofa/core/objectmodel/Data.h>
#include <sofa/defaulttype/BaseMatrix.h>
#include <sofa/defaulttype/BaseVector.h>
#include <sofa/core/objectmodel/Link.h>
#include <SofaBaseLinearSolver/DefaultMultiMatrixAccessor.h>
#include <memory>
 
namespace SofaCaribou::ode {
 
class NewtonRaphsonSolver : public sofa::core::behavior::OdeSolver {
public:
    SOFA_CLASS(NewtonRaphsonSolver, sofa::core::behavior::OdeSolver);
 
    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>;
 
    NewtonRaphsonSolver();
 
    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> & { return p_times; }
 
    auto squared_residuals() const -> const std::vector<FLOATING_POINT_TYPE> & { return p_squared_residuals; }
 
    auto squared_initial_residual() const -> const FLOATING_POINT_TYPE & { return p_squared_initial_residual; }
 
private:
 
    virtual 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) = 0;
 
    virtual void assemble_system_matrix(const sofa::core::MechanicalParams & mechanical_parameters,
                                        sofa::component::linearsolver::DefaultMultiMatrixAccessor & matrix_accessor,
                                        sofa::defaulttype::BaseMatrix * A) = 0;
 
    virtual 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) = 0;
 
    bool has_valid_linear_solver () const;
 
    Data<unsigned> d_newton_iterations;
    Data<double> d_correction_tolerance_threshold;
    Data<double> d_residual_tolerance_threshold;
 
    Link<sofa::core::behavior::LinearSolver> l_linear_solver;
 
    Data<bool> d_converged;
 
 
    std::unique_ptr<sofa::defaulttype::BaseMatrix> p_A;
 
    std::unique_ptr<sofa::defaulttype::BaseVector> p_DX;
 
    std::unique_ptr<sofa::defaulttype::BaseVector> p_F;
 
    sofa::core::MultiVecDerivId p_U_id;
 
    std::vector<UNSIGNED_INTEGER_TYPE> p_times;
 
    std::vector<FLOATING_POINT_TYPE> p_squared_residuals;
 
    FLOATING_POINT_TYPE p_squared_initial_residual;
};
}