RBFWithQuadraticLagrange.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/quadrature/Quadrature.hpp>
#include <polyfem/assembler/Assembler.hpp>
#include <Eigen/Dense>
 
namespace polyfem
{
    namespace basis
    {
        // See `RBFWithQuadratic.cpp` for a detail commented version of the code.
        // This class implements the same consistency constraints, but uses Lagrange
        // multipliers to solve the constrained least-square system, instead of
        // eliminating the constraints explicitly from the system.
        class RBFWithQuadraticLagrange
        {
        public:
            RBFWithQuadraticLagrange(const assembler::LinearAssembler &assembler, const Eigen::MatrixXd &centers, const Eigen::MatrixXd &collocation_points,
                                     const Eigen::MatrixXd &local_basis_integral, const quadrature::Quadrature &quadr,
                                     Eigen::MatrixXd &rhs, bool with_constraints = true);
 
            void basis(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) const;
 
            void grad(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) const;
 
            void bases_values(const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const;
 
            void bases_grads(const int axis, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const;
 
        private:
            bool is_volume() const { return centers_.cols() == 3; }
 
            // Computes the matrix that evaluates the kernels + polynomial terms on the given sample points
            void compute_kernels_matrix(const Eigen::MatrixXd &samples, Eigen::MatrixXd &A) const;
 
            // Computes the constraint matrix C that we want to impose (C w = d)
            void compute_constraints_matrix_2d_old(const int num_bases, const quadrature::Quadrature &quadr, Eigen::MatrixXd &C) const;
            void compute_constraints_matrix_2d(const assembler::LinearAssembler &assembler, const int num_bases, const quadrature::Quadrature &quadr, Eigen::MatrixXd &C) const;
 
            // Computes the constraint matrix C that we want to impose (C w = d)
            void compute_constraints_matrix_3d(const int num_bases, const quadrature::Quadrature &quadr, Eigen::MatrixXd &C) const;
 
            // Computes the weights by solving a (possibly constrained) linear least square
            void compute_weights(const assembler::LinearAssembler &assembler, const Eigen::MatrixXd &collocation_points,
                                 const Eigen::MatrixXd &local_basis_integral, const quadrature::Quadrature &quadr,
                                 Eigen::MatrixXd &rhs, bool with_constraints);
 
        private:
            // #C x dim matrix of kernel center positions
            Eigen::MatrixXd centers_;
 
            // (#C + dim + 1) x #B matrix of weights extending the #B bases that are non-vanishing on the polytope
            Eigen::MatrixXd weights_;
        };
    } // namespace basis
} // namespace polyfem