GMpoly_cell.hpp

Go to the documentation of this file.

# ifndef _GMPOLY_CELL_HPP
# define _GMPOLY_CELL_HPP
 
/*
 
  Classes to implement the quadrature-free polynomial integration rules
  
*/
 
#include <cassert>
#include <cmath>
 
#include <algorithm>
#include <array>
#include <iostream>
#include <vector>
#include <unordered_map>
 
#include <mesh.hpp>
//#include <cell.hpp>
//#include <edge.hpp>
//#include <vertex.hpp>
#include <quadraturerule.hpp>
#include <basis.hpp>
#include <typeinfo>
 
namespace HArDCore2D {
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
// INTEGRALS OF MONOMIALS
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
 
struct VecHash
{
  size_t operator()(const VectorZd & p) const
  {
    constexpr size_t b = 100; // max degree must be less than the basis b
    return p(0) + p(1) * b;
  }
};
 
typedef std::unordered_map<VectorZd, double, VecHash> MonomialCellIntegralsType;
    
 
MonomialCellIntegralsType IntegrateCellMonomials
          (const Cell & T,      
          const size_t maxdeg           
          );
 
MonomialCellIntegralsType CheckIntegralsDegree
         (const Cell & T,              
           const size_t degree,         
           const MonomialCellIntegralsType & mono_int_map = {}    
          );
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//  TRANSFORMATIONS OF GRAM MATRICES FOR DERIVED BASES
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
template<typename BasisType>
inline Eigen::MatrixXd transformGM(const Family<BasisType> & family_basis,     
                               const char RC,                
                               const Eigen::MatrixXd & anc_GM          
                               )
  {
    if (RC=='R'){
      return family_basis.matrix() * anc_GM;
    }else{
      return anc_GM * (family_basis.matrix()).transpose();
    }
  }
 
template<typename BasisType>
inline Eigen::MatrixXd transformGM(const RestrictedBasis<BasisType> & restr_basis,     
                               const char RC,                
                               const Eigen::MatrixXd & anc_GM              
                               )
  {
    if (RC=='R'){
      return anc_GM.topRows(restr_basis.dimension());
    }else{
      return anc_GM.leftCols(restr_basis.dimension());
    }
  }
 
template<typename BasisType>
inline Eigen::MatrixXd transformGM(const ShiftedBasis<BasisType> & shifted_basis,     
                               const char RC,                
                               const Eigen::MatrixXd & anc_GM              
                               )
  {
    if (RC=='R'){
      return anc_GM.bottomRows(shifted_basis.dimension());
    }else{
      return anc_GM.rightCols(shifted_basis.dimension());
    }
  }
 
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//  FUNCTIONS TO COMPUTE GRAM MATRICES
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
 
/* BASIC ONES: Monomial, tensorized, and generic template */
 
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                         
                    const MonomialScalarBasisCell & basis1, 
                    const MonomialScalarBasisCell & basis2, 
                    MonomialCellIntegralsType mono_int_map = {} 
                    );
  
template<typename BasisType>
Eigen::MatrixXd GramMatrix(const Cell& T, const BasisType & basis, MonomialCellIntegralsType mono_int_map = {})
  {
    return GramMatrix(T, basis, basis, mono_int_map);
  };
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrix(
                     const Cell& T, 
                     const TensorizedVectorFamily<BasisType1, N> & basis1, 
                     const TensorizedVectorFamily<BasisType2, N> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
    Eigen::MatrixXd anc_gm = GramMatrix(T, 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;
  };
 
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                         
                    const RolyComplBasisCell & basis1,      
                    const RolyComplBasisCell & basis2,      
                    MonomialCellIntegralsType mono_int_map = {} 
                    );
 
template<typename BasisType1, size_t N>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                                         
                    const RolyComplBasisCell & rolycompl_basis,             
                    const TensorizedVectorFamily<BasisType1, N> & tens_family,  
                    MonomialCellIntegralsType mono_int_map = {}                 
                    )
  {
    size_t dim1 = rolycompl_basis.dimension();
    size_t dim2 = tens_family.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = rolycompl_basis.max_degree()+tens_family.max_degree();
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
    for (size_t m=0; m<N; m++){
      gm.block(0, m*dim2/N, dim1, dim2/N) = GMRolyComplScalar(T, rolycompl_basis, tens_family.ancestor(), m, intmap);
    }
    return gm;
  };
  
template<typename BasisType1, size_t N>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                                         
                    const TensorizedVectorFamily<BasisType1, N> & tens_family,  
                    const RolyComplBasisCell & rolycompl_basis,             
                    MonomialCellIntegralsType mono_int_map = {}                 
                    )
  {
    return GramMatrix(T, rolycompl_basis, tens_family, mono_int_map).transpose();
  };
  
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                         
                    const GolyComplBasisCell & basis1,      
                    const GolyComplBasisCell & basis2,      
                    MonomialCellIntegralsType mono_int_map = {} 
                    );
  
Eigen::MatrixXd GMRolyComplScalar(
                    const Cell & T, 
                    const RolyComplBasisCell & rolycompl_basis, 
                    const MonomialScalarBasisCell & mono_basis, 
                    const size_t m, 
                    MonomialCellIntegralsType mono_int_map 
                    );
 
template<typename BasisType>
Eigen::MatrixXd GMRolyComplScalar(
                     const Cell& T, 
                     const RolyComplBasisCell & basis1, 
                     const BasisType & basis2,  
                     const size_t m, 
                     MonomialCellIntegralsType mono_int_map 
                     )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(BasisType::hasAncestor, "No method to compute this Gram matrix of derivatives");
    return transformGM(basis2, 'C', GMRolyComplScalar(T, basis1, basis2.ancestor(), m, mono_int_map) );
  };  
  
template<typename BasisType>
constexpr bool useAncestor()
  {
    // For a given basis, we only use the ancestor if it has an ancestor and has the same rank as its
    //  ancestor. If it has a different rank than its ancestor, it must be dealt with a specific overload
 
    if constexpr(BasisType::hasAncestor){
      return (BasisType::tensorRank == BasisType::AncestorType::tensorRank);
    } else {
      return false;
    }
  };
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(const Cell& T,         
                     const BasisType1 & basis1,   
                     const BasisType2 & basis2,   
                     MonomialCellIntegralsType 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(T, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrix(T, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrix(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
 
  };
 
 
/* Gram matrix of scalar basis with one or two derivatives */
 
Eigen::MatrixXd GMScalarDerivative(
                  const Cell & T,  
                  const MonomialScalarBasisCell & basis1, 
                  const MonomialScalarBasisCell & basis2, 
                  const size_t m, 
                  MonomialCellIntegralsType mono_int_map = {} 
                  );
  
Eigen::MatrixXd GMScalarDerivative(
                  const Cell & T,  
                  const MonomialScalarBasisCell & basis1, 
                  const MonomialScalarBasisCell & basis2, 
                  const size_t m, 
                  const size_t l, 
                  MonomialCellIntegralsType mono_int_map = {} 
                  );
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GMScalarDerivative(
                     const Cell& T, 
                     const BasisType1 & basis1, 
                     const BasisType2 & basis2,  
                     const size_t m, 
                     MonomialCellIntegralsType mono_int_map = {} 
                             )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(BasisType1::hasAncestor || BasisType2::hasAncestor, "No method to compute this Gram matrix of derivatives");
      
    if constexpr (!BasisType1::hasAncestor && BasisType2::hasAncestor) {
      return transformGM(basis2, 'C', GMScalarDerivative(T, basis1, basis2.ancestor(), m, mono_int_map) );
    } else if constexpr (BasisType1::hasAncestor && !BasisType2::hasAncestor) {
      return transformGM(basis1, 'R', GMScalarDerivative(T, basis1.ancestor(), basis2, m, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), m, mono_int_map) ) );
    }
  };
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GMScalarDerivative(
                     const Cell& T, 
                     const BasisType1 & basis1, 
                     const BasisType2 & basis2,  
                     const size_t m, 
                     const size_t l, 
                     MonomialCellIntegralsType mono_int_map = {} 
                             )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(BasisType1::hasAncestor || BasisType2::hasAncestor, "No method to compute this Gram matrix of derivatives");
      
    if constexpr (!BasisType1::hasAncestor && BasisType2::hasAncestor) {
      return transformGM(basis2, 'C', GMScalarDerivative(T, basis1, basis2.ancestor(), m, l, mono_int_map) );
    } else if constexpr (BasisType1::hasAncestor && !BasisType2::hasAncestor) {
      return transformGM(basis1, 'R', GMScalarDerivative(T, basis1.ancestor(), basis2, m, l, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), m, l, mono_int_map) ) );
    }
    
  };
  
 
/* Gram matrices of vector-valued basis with gradient, curl... */
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrix(
                     const Cell& T, 
                     const GradientBasis<BasisType1> & grad_basis, 
                     const TensorizedVectorFamily<BasisType2, N> & tens_family,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = grad_basis.dimension();
    size_t dim2 = tens_family.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = grad_basis.ancestor().max_degree()+tens_family.max_degree()-1;
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
    for (size_t m=0; m<N; m++){
      gm.block(0, m*dim2/N, dim1, dim2/N) = GMScalarDerivative(T, grad_basis.ancestor(), tens_family.ancestor(), m, intmap);
    }
    return gm/T.diam();
  };
 
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrix(
             const Cell& T, 
             const TensorizedVectorFamily<BasisType1, N> & tens_family, 
             const GradientBasis<BasisType2> & grad_basis,  
             MonomialCellIntegralsType mono_int_map = {} 
             )
  {
    return GramMatrix(T, grad_basis, tens_family, mono_int_map).transpose();
  };
 
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                     const Cell& T, 
                     const GradientBasis<BasisType1> & grad_basis1, 
                     const GradientBasis<BasisType2> & grad_basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = grad_basis1.dimension();
    size_t dim2 = grad_basis2.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = grad_basis1.ancestor().max_degree()+grad_basis2.ancestor().max_degree()-2;
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
    for (size_t m=0; m<2; m++){
      gm += GMScalarDerivative(T, grad_basis1.ancestor(), grad_basis2.ancestor(), m, m, intmap);
    }
    return gm/std::pow(T.diam(), 2);
  };
  
  
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                             
                    const CurlBasis<BasisType1> & basis1,       
                    const CurlBasis<BasisType2> & basis2,       
                    MonomialCellIntegralsType mono_int_map = {} 
                    )
  {
    // Dimension of the gram matrix
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    size_t totaldegree = basis1.ancestor().max_degree()+basis2.ancestor().max_degree()-2;
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1,dim2);
    
    // Obtain integration data from IntegrateCellMonomials
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
    
    for (size_t m=0; m<2; m++){
      gm += GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), m, m, intmap);
    }
    
    return gm/std::pow(T.diam(), 2);
  };
  
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                     const Cell& T, 
                     const CurlBasis<BasisType1> & curl_basis, 
                     const TensorizedVectorFamily<BasisType2, 2> & tens_family,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = curl_basis.dimension();
    size_t dim2 = tens_family.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = curl_basis.ancestor().max_degree()+tens_family.max_degree()-1;
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
    gm.block(0, 0, dim1, dim2/2) = GMScalarDerivative(T, curl_basis.ancestor(), tens_family.ancestor(), 1, intmap);
    gm.block(0, dim2/2, dim1, dim2/2) = -GMScalarDerivative(T, curl_basis.ancestor(), tens_family.ancestor(), 0, intmap);
 
    return gm/T.diam();
  };
 
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrix(
             const Cell& T, 
             const TensorizedVectorFamily<BasisType1, N> & tens_family, 
             const CurlBasis<BasisType2> & curl_basis,  
             MonomialCellIntegralsType mono_int_map = {} 
             )
  {
    return GramMatrix(T, curl_basis, tens_family, mono_int_map).transpose();
  };
 
 
template<typename BasisType1, typename BasisType2>
typename boost::disable_if<boost::is_same<BasisType2, MonomialCellIntegralsType>, Eigen::MatrixXd>::type GramMatrix(
                     const Cell& T,                          
                     const DivergenceBasis<BasisType1> & basis1,   
                     const BasisType2 & basis2,              
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    return GramMatrixDiv(T, basis1.ancestor(), basis2, mono_int_map);
  };
  
template<typename BasisType1, typename Basis2>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                         
                    const BasisType1 & basis1,                  
                    const DivergenceBasis<Basis2> & basis2,       
                    MonomialCellIntegralsType mono_int_map = {} 
                    )
  {
    return GramMatrixDiv(T, basis2.ancestor(), basis1, mono_int_map).transpose();
  };
  
template<typename BasisType1>
Eigen::MatrixXd GramMatrixDiv(
                    const Cell & T,                                         
                    const TensorizedVectorFamily<BasisType1, 2> & basis1,   
                    const MonomialScalarBasisCell & basis2,                 
                    MonomialCellIntegralsType mono_int_map = {}                 
                    )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree()-1;
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
    for (size_t i = 0; i < 2; i++) {
      gm.block(i*dim1/2, 0, dim1/2, dim2) = GMScalarDerivative(T, basis1.ancestor(), basis2, i, intmap);
    }
    
    return gm/T.diam();
  };
  
Eigen::MatrixXd GramMatrixDiv(
                    const Cell & T,                           
                    const RolyComplBasisCell & basis1,        
                    const MonomialScalarBasisCell & basis2,   
                    MonomialCellIntegralsType mono_int_map = {}   
                    );
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixDiv(
                     const Cell& T, 
                     const BasisType1 & basis1, 
                     const BasisType2 & basis2,  
                     MonomialCellIntegralsType 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', GramMatrixDiv(T, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrixDiv(T, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrixDiv(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
  };
 
/*------------- Templates for MatrixFamily  --------------------*/
 
    template<typename BasisType1, typename BasisType2, size_t N>
    Eigen::MatrixXd GramMatrix(
            const Cell& T, 
            const MatrixFamily<BasisType1, N> & basis1, 
            const MatrixFamily<BasisType2, N> & basis2,  
            MonomialCellIntegralsType mono_int_map = {} 
    )
    {
        Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
        Eigen::MatrixXd anc_gm = GramMatrix(T, 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*N; i++){
            gm.block(i*dim1, i*dim2, dim1, dim2) = anc_gm;
        }
        return gm;
    }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
            const Cell& T,                          
            const DivergenceBasis<MatrixFamily<BasisType1, dimspace>> & basis1,   
            const TensorizedVectorFamily<BasisType2, dimspace> & basis2,              
            MonomialCellIntegralsType mono_int_map = {} 
    )
    {
        Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
        // Grab the scalar bases underlying basis1 and basis2
        BasisType1 scalar_basis1 = basis1.ancestor().ancestor();
        BasisType2 scalar_basis2 = basis2.ancestor();
        size_t dim1 = scalar_basis1.dimension();
        size_t dim2 = scalar_basis2.dimension();
 
        // Integrals of monomials
        size_t totaldegree = scalar_basis1.max_degree() + scalar_basis2.max_degree();
        MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
        for (size_t j=0; j < dimspace; j++){
            Eigen::MatrixXd gm_anc = GMScalarDerivative(T, scalar_basis1, scalar_basis2, j, intmap);
            for (size_t i=0; i < dimspace; i++){
                gm.block(j*dimspace*dim1 + i*dim1, i*dim2, dim1, dim2) = gm_anc;
            }
        }
 
        return gm/T.diam();
    }
 
    template<typename BasisType1, typename BasisType2>
    Eigen::MatrixXd GramMatrix(
            const Cell& T,                          
            const GradientBasis<TensorizedVectorFamily<BasisType1, dimspace>> & basis1,   
            const MatrixFamily<BasisType2, dimspace> & basis2,              
            MonomialCellIntegralsType mono_int_map = {} 
    )
    {
        Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
        // Grab the scalar bases underlying basis1 and basis2
        BasisType1 scalar_basis1 = basis1.ancestor().ancestor();
        BasisType2 scalar_basis2 = basis2.ancestor();
        size_t dim1 = scalar_basis1.dimension();
        size_t dim2 = scalar_basis2.dimension();
 
        // Integrals of monomials
        size_t totaldegree = scalar_basis1.max_degree() + scalar_basis2.max_degree();
        MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
        for (size_t j=0; j < dimspace; j++){
            Eigen::MatrixXd gm_anc = GMScalarDerivative(T, scalar_basis1, scalar_basis2, j, intmap);
            for (size_t i=0; i < dimspace; i++){
                gm.block(i*dim1, j*dimspace*dim2 + i*dim2, dim1, dim2) = gm_anc;
            }
        }
 
        return gm/T.diam();
    }
 
    template<typename BasisType1, typename BasisType2>
    Eigen::MatrixXd GramMatrix(
            const Cell& T,                          
            const MatrixFamily<BasisType1, dimspace> & basis1,              
            const GradientBasis<TensorizedVectorFamily<BasisType2, dimspace>> & basis2,   
            MonomialCellIntegralsType mono_int_map = {} 
    )
    {
        return GramMatrix(T, basis2, basis1, mono_int_map).transpose();
    }
 
    template<typename BasisType1, typename BasisType2, size_t N>
    Eigen::MatrixXd GramMatrix(
            const Cell& T,                          
            const GradientBasis<TensorizedVectorFamily<BasisType1, N>> & basis1,   
            const GradientBasis<TensorizedVectorFamily<BasisType2, N>> & basis2,   
            MonomialCellIntegralsType mono_int_map = {} 
    )
    {
        Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
        // Grab the scalar bases underlying basis1 and basis2
        BasisType1 scalar_basis1 = basis1.ancestor().ancestor();
        BasisType2 scalar_basis2 = basis2.ancestor().ancestor();
        size_t dim1 = scalar_basis1.dimension();
        size_t dim2 = scalar_basis2.dimension();
 
        // Integrals of monomials
        size_t totaldegree = scalar_basis1.max_degree() + scalar_basis2.max_degree();
        MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
        // Gram matrices of sum_i \partial_i phi \partial_i psi for the two scalar bases phi, psi
        Eigen::MatrixXd gm_anc = Eigen::MatrixXd::Zero(dim1, dim2);
        for (size_t k=0; k < dimspace; k++){
            gm_anc += GMScalarDerivative(T, scalar_basis1, scalar_basis2, k, k, intmap);
        }
 
        for (size_t i=0; i < N; i++){
            gm.block(i*dim1, i*dim2, dim1, dim2) = gm_anc;
        }
 
        return gm/(std::pow(T.diam(),2));
    }
 
 
} // end of namespace HArDCore2D
 
#endif // end of _GMPOLY_CELL_HPP