Integrator.h

Go to the documentation of this file.

 
#ifndef PBAT_SIM_ALGORITHM_NEWTON_INTEGRATOR_H
#define PBAT_SIM_ALGORITHM_NEWTON_INTEGRATOR_H
 
#include "Config.h"
#include "pbat/Aliases.h"
#include "pbat/fem/Tetrahedron.h"
#include "pbat/io/Archive.h"
#include "pbat/math/optimization/LineSearch.h"
#include "pbat/math/optimization/Newton.h"
#include "pbat/physics/StableNeoHookeanEnergy.h"
#include "pbat/sim/contact/MultibodyTetrahedralMeshSystem.h"
#include "pbat/sim/dynamics/FemElastoDynamics.h"
 
#ifdef PBAT_USE_SUITESPARSE
    #include <Eigen/CholmodSupport>
#else
    #include <Eigen/SparseCholesky>
#endif // PBAT_USE_SUITESPARSE
#include <memory>
#include <optional>
#include <vector>
 
namespace pbat::sim::algorithm::newton {

class Integrator
{
  public:
    static auto constexpr kDims  = 3; 
    static auto constexpr kOrder = 1; 
 
    using ScalarType     = Scalar; 
    using IndexType      = Index;  
    using MeshSystemType = contact::MultibodyTetrahedralMeshSystem<IndexType>; 
    using ElastoDynamicsType = dynamics::FemElastoDynamics<
        fem::Tetrahedron<1>,
        kDims,
        physics::StableNeoHookeanEnergy<kDims>,
        ScalarType,
        IndexType>; 
    using HessianMatrixType =
        Eigen::SparseMatrix<ScalarType, Eigen::ColMajor, IndexType>; 
#ifdef PBAT_USE_SUITESPARSE
    using DecompositionType =
        Eigen::CholmodDecomposition<HessianMatrixType, Eigen::Upper>; 
#else
    using DecompositionType =
        Eigen::SimplicialLDLT<HessianMatrixType, Eigen::Upper>; 
#endif // PBAT_USE_SUITESPARSE
    Integrator(Integrator const& other);
    Integrator(Integrator&& other) noexcept = default;
    Integrator& operator=(Integrator const& other);
    Integrator& operator=(Integrator&& other) noexcept = default;
    Integrator(Config config, MeshSystemType meshSystem, ElastoDynamicsType elastoDynamics);
    void Step(std::optional<io::Archive> archive = std::nullopt);
    [[maybe_unused]] auto GetElastoDynamics() const -> ElastoDynamicsType const&
    {
        return mElastoDynamics;
    }

    [[maybe_unused]] auto GetElastoDynamics() -> ElastoDynamicsType& { return mElastoDynamics; }
    [[maybe_unused]] auto GetMeshSystem() const -> MeshSystemType const& { return mMeshes; }
    [[maybe_unused]] auto GetMeshSystem() -> MeshSystemType& { return mMeshes; }
    [[maybe_unused]] auto GetConfig() const -> Config const& { return mConfig; }
    [[maybe_unused]] auto GetConfig() -> Config& { return mConfig; }
 
  protected:
    void AssembleHessian(ScalarType bt2);
    void ApplyConfig();
 
  private:
    Config mConfig; 
    contact::MultibodyTetrahedralMeshSystem<IndexType>
        mMeshes;                                    
    ElastoDynamicsType mElastoDynamics;             
    math::optimization::Newton<ScalarType> mNewton; 
    math::optimization::BackTrackingLineSearch<ScalarType> mLineSearch; 
    std::vector<Eigen::Triplet<ScalarType, IndexType>>
        mTriplets;                                      
    std::unique_ptr<DecompositionType> mInverseHessian; 
    HessianMatrixType mHessian;                         
    Eigen::Vector<ScalarType, Eigen::Dynamic> mGrad;    
};
 
} // namespace pbat::sim::algorithm::newton
 
#endif // PBAT_SIM_ALGORITHM_NEWTON_INTEGRATOR_H