BCLagrangianForm.hpp

Go to the documentation of this file.

#pragma once
 
#include "AugmentedLagrangianForm.hpp"
#include <polyfem/assembler/RhsAssembler.hpp>
 
#include <polyfem/mesh/Obstacle.hpp>
#include <polyfem/mesh/LocalBoundary.hpp>
#include <polyfem/utils/JSONUtils.hpp>
 
namespace polyfem::solver
{
    enum class BCLumpingMode
    {
        ROW_SUM, 
        HRZ      
    };
 
    NLOHMANN_JSON_SERIALIZE_ENUM(
        BCLumpingMode,
        {{BCLumpingMode::ROW_SUM, "row_sum"},
         {BCLumpingMode::HRZ, "hrz"}});
 
    class BCLagrangianForm : public AugmentedLagrangianForm
    {
    public:
        BCLagrangianForm(const int ndof,
                         const std::vector<int> &boundary_nodes,
                         const std::vector<mesh::LocalBoundary> &local_boundary,
                         const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
                         const QuadratureOrders &n_boundary_samples,
                         const StiffnessMatrix &mass,
                         const assembler::RhsAssembler &rhs_assembler,
                         const size_t obstacle_ndof,
                         const bool is_time_dependent,
                         const double t,
                         const BCLumpingMode lumping = BCLumpingMode::ROW_SUM);
 
        std::string name() const override
        {
            return "bc-alagrangian";
        }
 
        BCLagrangianForm(const int ndof,
                         const std::vector<int> &boundary_nodes,
                         const StiffnessMatrix &mass,
                         const size_t obstacle_ndof,
                         const Eigen::MatrixXd &target_x);
 
        double value_unweighted(const Eigen::VectorXd &x) const override;
 
        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:
        void update_quantities(const double t, const Eigen::VectorXd &x) override;
 
        void update_lagrangian(const Eigen::VectorXd &x, const double k_al) override;
 
        double compute_error(const Eigen::VectorXd &x) const override;
 
        virtual bool can_project() const override;
        virtual void project_gradient(Eigen::VectorXd &grad) const override;
        virtual void project_diag(Eigen::VectorXd &diag) const override;
        virtual void project_hessian(StiffnessMatrix &hessian) const override;
 
    private:
        const std::vector<int> &boundary_nodes_;
        const std::vector<mesh::LocalBoundary> *local_boundary_;
        const std::vector<mesh::LocalBoundary> *local_neumann_boundary_;
        const QuadratureOrders n_boundary_samples_;
        const int n_dofs_;
        Eigen::VectorXi constraints_;     
        Eigen::VectorXi not_constraints_; 
        std::vector<int> old_to_new_;
 
        const assembler::RhsAssembler *rhs_assembler_; 
        const bool is_time_dependent_;
 
        StiffnessMatrix masked_lumped_mass_sqrt_; 
        StiffnessMatrix masked_lumped_mass_;      
 
        BCLumpingMode lumping_ = BCLumpingMode::ROW_SUM; 
 
        void init_masked_lumped_mass(
            const StiffnessMatrix &mass,
            const size_t obstacle_ndof);
 
        void update_target(const double t);
    };
} // namespace polyfem::solver