ElasticForm.hpp

Go to the documentation of this file.

#pragma once
 
#include "Form.hpp"
 
#include <polyfem/basis/ElementBases.hpp>
#include <polyfem/assembler/Assembler.hpp>
#include <polyfem/assembler/AssemblyValsCache.hpp>
 
#include <polyfem/utils/Jacobian.hpp>
#include <polyfem/utils/MatrixCache.hpp>
#include <polyfem/utils/Types.hpp>
 
#include <memory>
#include <optional>
#include <utility>
#include <vector>
 
namespace polyfem::solver
{
    enum class ElementInversionCheck
    {
        Discrete,
        Conservative,
    };
    NLOHMANN_JSON_SERIALIZE_ENUM(
        polyfem::solver::ElementInversionCheck,
        {{ElementInversionCheck::Discrete, "Discrete"},
         {ElementInversionCheck::Conservative, "Conservative"}})
 
    class ElasticForm : public Form
    {
        friend class ElasticForceDerivative;
 
    public:
        ElasticForm(const int n_bases,
                    std::vector<basis::ElementBases> &bases,
                    const std::vector<basis::ElementBases> &geom_bases,
                    const assembler::Assembler &assembler,
                    assembler::AssemblyValsCache &ass_vals_cache,
                    const double t, const double dt,
                    const bool is_volume,
                    const double jacobian_threshold = 0.,
                    const ElementInversionCheck check_inversion = ElementInversionCheck::Discrete,
                    const unsigned conservative_max_iter = 100000);
 
        std::string name() const override { return "elastic"; }
 
    protected:
        virtual double value_unweighted(const Eigen::VectorXd &x) const override;
 
        Eigen::VectorXd value_per_element_unweighted(const Eigen::VectorXd &x) const override;
 
        virtual void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override;
 
        void second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override;
 
    public:
        bool is_step_valid(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override;
 
        bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override;
 
        void update_quantities(const double t, const Eigen::VectorXd &x) override
        {
            t_ = t;
            x_prev_ = x;
        }
 
        double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override;
 
        void finish() override;
 
        void post_step(const polysolve::nonlinear::PostStepData &data) override;
 
    private:
        const int n_bases_;
        std::vector<basis::ElementBases> &bases_;
        const std::vector<basis::ElementBases> &geom_bases_;
 
        const assembler::Assembler &assembler_; 
        assembler::AssemblyValsCache &ass_vals_cache_;
        double t_;
        const double jacobian_threshold_;
        const ElementInversionCheck check_inversion_;
        const unsigned conservative_max_iter_;
        const double dt_;
        const bool is_volume_;
 
        StiffnessMatrix cached_stiffness_;                      
        mutable std::unique_ptr<utils::MatrixCache> mat_cache_; 
 
        void compute_cached_stiffness();
 
        void commit_refinement(const int id, utils::Tree &subdivision_tree) const;
 
        Eigen::VectorXd x_prev_;
 
        mutable std::vector<utils::Tree> quadrature_hierarchy_;
        int quadrature_order_;
 
        mutable std::optional<std::pair<int, utils::Tree>> pending_refinement_;
    };
} // namespace polyfem::solver