GMpoly_edge.hpp

Go to the documentation of this file.

# ifndef _GMPOLY_EDGE_HPP
# define _GMPOLY_EDGE_HPP
 
/*
 
  Classes to implement the quadrature-free polynomial integration rules
  
*/
 
#include <cassert>
#include <cmath>
 
#include <mesh.hpp>
#include <quadraturerule.hpp>
#include <basis.hpp>
#include <GMpoly_cell.hpp>
 
namespace HArDCore2D {
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
// INTEGRALS OF MONOMIALS
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
 
typedef std::vector<double> MonomialEdgeIntegralsType;
 
MonomialEdgeIntegralsType IntegrateEdgeMonomials
          (const Edge & E,      
          const size_t maxdeg   
          );
 
 
MonomialEdgeIntegralsType CheckIntegralsDegree
          (const Edge & E,              
           const size_t degree,         
           const MonomialEdgeIntegralsType & mono_int_map = {}    
          );
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//  FUNCTIONS TO COMPUTE GRAM MATRICES
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
//
/* SCALAR bases */
 
Eigen::MatrixXd GramMatrix(
                    const Edge & E,                         
                    const MonomialScalarBasisEdge & basis1, 
                    const MonomialScalarBasisEdge & basis2, 
                    MonomialEdgeIntegralsType mono_int_map = {} 
                    );
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(const Edge& E,         
                     const BasisType1 & basis1,   
                     const BasisType2 & basis2,   
                     MonomialEdgeIntegralsType mono_int_map = {} 
                     )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(useAncestor<BasisType1>() || useAncestor<BasisType2>(), "No method to compute this Gram matrix");
      
    if constexpr (!useAncestor<BasisType1>() && useAncestor<BasisType2>()) {
      return transformGM(basis2, 'C', GramMatrix(E, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrix(E, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrix(E, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
 
  };
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrix(
            const Edge& E, 
            const TensorizedVectorFamily<BasisType1, N> & basis1, 
            const TensorizedVectorFamily<BasisType2, N> & basis2,  
            MonomialEdgeIntegralsType mono_int_map = {} 
    )
    {
        Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
        Eigen::MatrixXd anc_gm = GramMatrix(E, basis1.ancestor(), basis2.ancestor(), mono_int_map);
        size_t dim1 = anc_gm.rows();
        size_t dim2 = anc_gm.cols();
 
        for (size_t i=0; i<N; i++){
            gm.block(i*dim1, i*dim2, dim1, dim2) = anc_gm;
        }
        return gm;
    }
 
template<typename BasisType>
Eigen::MatrixXd GramMatrix(const Edge& E, const BasisType & basis, MonomialEdgeIntegralsType mono_int_map = {})
  {
    return GramMatrix(E, basis, basis, mono_int_map);
  };
 
/* GRADIENT-SCALAR bases */ 
 
Eigen::MatrixXd GMDer(
                  const Edge & E,  
                  const MonomialScalarBasisEdge & basis1, 
                  const MonomialScalarBasisEdge & basis2, 
                  MonomialEdgeIntegralsType mono_int_map = {} 
                  );
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GMDer(const Edge& E,         
                     const BasisType1 & basis1,   
                     const BasisType2 & basis2,   
                     MonomialEdgeIntegralsType mono_int_map = {} 
                     )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(useAncestor<BasisType1>() || useAncestor<BasisType2>(), "No method to compute this Gram matrix");
      
    if constexpr (!useAncestor<BasisType1>() && useAncestor<BasisType2>()) {
      return transformGM(basis2, 'C', GMDer(E, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GMDer(E, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GMDer(E, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
 
  };
 
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Edge & E,                         
                    const GradientBasis<BasisType1> & basis1, 
                    const BasisType2 & basis2, 
                    MonomialEdgeIntegralsType mono_int_map = {} 
                    )
  {
    return GMDer(E, basis1.ancestor(), basis2, mono_int_map)/E.diam();
  };
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Edge & E,                         
                    const BasisType1 & basis1, 
                    const GradientBasis<BasisType2> & basis2, 
                    MonomialEdgeIntegralsType mono_int_map = {} 
                    )
  {
    return (GMDer(E, basis2.ancestor(), basis1, mono_int_map).transpose())/E.diam();
  };
 
 
 
 
} // end of namespace HArDCore2D
 
#endif // end of _GMPOLY_EDGE_HPP