EigenSparseSolver.h

#pragma once
 
#include <Caribou/config.h>
#include <SofaCaribou/Algebra/EigenMatrix.h>
#include <SofaCaribou/Algebra/EigenVector.h>
#include <SofaCaribou/Solver/LinearSolver.h>
 
#include <sofa/core/behavior/LinearSolver.h>
#include <SofaBaseLinearSolver/DefaultMultiMatrixAccessor.h>
#include <Eigen/Sparse>
 
namespace SofaCaribou::solver {
 
namespace internal {
template <typename T>
struct solver_traits {};
}
 
template <class EigenSolver_t>
class EigenSparseSolver : public sofa::core::behavior::LinearSolver, public SofaCaribou::solver::LinearSolver {
public:
    SOFA_CLASS(SOFA_TEMPLATE(EigenSparseSolver, EigenSolver_t), sofa::core::behavior::LinearSolver);
    using EigenSolver = EigenSolver_t;
    using SparseMatrix = std::remove_cv_t<typename EigenSolver::MatrixType>;
    using Vector = Eigen::Matrix<FLOATING_POINT_TYPE, Eigen::Dynamic, 1>;
 
    virtual void assemble (const sofa::core::MechanicalParams* mparams);
 
    void resetSystem() final;
 
    void setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) final;
 
    void setSystemRHVector(sofa::core::MultiVecDerivId b_id) final;
 
    void setSystemLHVector(sofa::core::MultiVecDerivId x_id) final;
 
    void solveSystem() override;
 
    inline virtual auto symmetric() const -> bool {return p_is_symmetric;}
 
    inline virtual void set_symmetric(bool is_symmetric) { p_is_symmetric = is_symmetric; }
 
    auto solver() const -> const EigenSolver & { return p_solver; }
 
    auto solver() -> EigenSolver & { return p_solver; }
 
    auto mechanical_params() const -> const sofa::core::MechanicalParams & { return p_mechanical_params; }
 
 
    auto matrix_accessor() const -> const sofa::component::linearsolver::DefaultMultiMatrixAccessor & { return p_accessor; }
 
    auto A() const -> const SparseMatrix & { return p_A; }
 
    auto b_id() const -> const sofa::core::MultiVecDerivId & { return p_b_id; }
 
    auto x_id() const -> const sofa::core::MultiVecDerivId & { return p_x_id; }
 
    auto A_is_factorized() const -> bool { return p_A_is_factorized; }
 
    // SOFA overrides
    static auto GetCustomTemplateName() -> std::string;
 
    template<class Derived>
    static auto canCreate(Derived* o, sofa::core::objectmodel::BaseContext* context, sofa::core::objectmodel::BaseObjectDescription* arg) -> bool;
 
private:
    sofa::defaulttype::BaseMatrix * create_new_matrix(sofa::Size rows, sofa::Size cols) const override {
        auto * matrix = new SofaCaribou::Algebra::EigenMatrix<SparseMatrix> (static_cast<Eigen::Index>(rows), static_cast<Eigen::Index>(cols));
        if (symmetric()) {
            matrix->set_symmetric(symmetric());
        }
        return matrix;
    }
 
    sofa::defaulttype::BaseVector * create_new_vector(sofa::Size n) const override {
        return new SofaCaribou::Algebra::EigenVector<Vector>(n);
    }
 
    bool solve(const sofa::defaulttype::BaseMatrix * A,
               const sofa::defaulttype::BaseVector * F,
               sofa::defaulttype::BaseVector * X) const override;
 
private:
 
    EigenSolver p_solver;
 
    sofa::core::MechanicalParams p_mechanical_params;
 
    sofa::component::linearsolver::DefaultMultiMatrixAccessor p_accessor;
 
    sofa::core::MultiVecDerivId p_b_id;
 
    sofa::core::MultiVecDerivId p_x_id;
 
    SparseMatrix p_A;
 
    Vector p_x;
 
    Vector p_b;
 
    bool p_A_is_factorized;
 
    bool p_is_symmetric = false;
};
 
} // namespace SofaCaribou::solver