RBFWithQuadratic.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/quadrature/Quadrature.hpp>
#include <polyfem/assembler/ElementAssemblyValues.hpp>
#include <polyfem/assembler/Assembler.hpp>
 
#include <Eigen/Dense>
 
namespace polyfem
{
    namespace basis
    {
        class RBFWithQuadratic
        {
        public:
            inline static int index_mapping(const int alpha, const int beta, const int d, const int ass_dim)
            {
                return ass_dim * ass_dim * d + ass_dim * beta + alpha;
            }
 
            static void setup_monomials_vals_2d(const int star_index, const Eigen::MatrixXd &pts, assembler::ElementAssemblyValues &vals);
            static void setup_monomials_strong_2d(const int dim, const assembler::LinearAssembler &assembler, const Eigen::MatrixXd &pts, const QuadratureVector &da, std::array<Eigen::MatrixXd, 5> &strong);
 
            RBFWithQuadratic(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 relationship w = L v + t between the unknowns (v) and the weights w
            void compute_constraints_matrix_2d_old(const int num_bases, const quadrature::Quadrature &quadr,
                                                   const Eigen::MatrixXd &local_basis_integral, Eigen::MatrixXd &L, Eigen::MatrixXd &t) const;
 
            // Computes the relationship w = L v + t between the unknowns (v) and the weights w
            void compute_constraints_matrix_2d(const assembler::LinearAssembler &assembler, const int num_bases, const quadrature::Quadrature &quadr,
                                               const Eigen::MatrixXd &local_basis_integral, Eigen::MatrixXd &L, Eigen::MatrixXd &t) const;
 
            // Computes the relationship w = L v + t between the unknowns (v) and the weights w
            void compute_constraints_matrix_3d(const assembler::LinearAssembler &assembler, const int num_bases, const quadrature::Quadrature &quadr,
                                               const Eigen::MatrixXd &local_basis_integral, Eigen::MatrixXd &L, Eigen::MatrixXd &t) 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