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 <quadraturerule.hpp>
#include <basis.hpp>
#include <typeinfo>
 
namespace HArDCore3D {
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
// 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 + p(2) * b*b;
  }
};
 
typedef std::unordered_map<VectorZd, double, VecHash> MonomialCellIntegralsType;
    
std::vector<MonomialCellIntegralsType> IntegrateCellMonomials_onEdges
          (const Cell & T,      
          const size_t maxdeg           
          );
 
std::vector<MonomialCellIntegralsType> IntegrateCellMonomials_onFaces
          (const Cell & T,      
          const size_t maxdeg,           
          std::vector<MonomialCellIntegralsType> & integrals_edges  
          );
 
MonomialCellIntegralsType IntegrateCellMonomials
          (const Cell & T,      
          const int 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 = {} 
                    );
 
template<typename BasisType1, size_t N>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                                         
                    const GolyComplBasisCell & golycompl_basis,             
                    const TensorizedVectorFamily<BasisType1, N> & tens_family,  
                    MonomialCellIntegralsType mono_int_map = {}                 
                    )
  {
    size_t dim1 = golycompl_basis.dimension();
    size_t dim2 = tens_family.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = golycompl_basis.max_degree()+tens_family.max_degree();
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
 
    size_t dim_l1 = golycompl_basis.dimPkmo();
    size_t dim_l2 = tens_family.ancestor().dimension();
    gm.block(0, dim_l2, dim_l1, dim_l2) = GMGolyComplScalar(T, golycompl_basis, tens_family.ancestor(), 0, intmap, 2);
    gm.block(0, 2*dim_l2, dim_l1, dim_l2) = -GMGolyComplScalar(T, golycompl_basis, tens_family.ancestor(), 0, intmap, 1);
    gm.block(dim_l1, 0, dim_l1, dim_l2) = -GMGolyComplScalar(T, golycompl_basis, tens_family.ancestor(), 1, intmap, 2);
    gm.block(dim_l1, 2*dim_l2, dim_l1, dim_l2) = GMGolyComplScalar(T, golycompl_basis, tens_family.ancestor(), 1, intmap, 0);
    gm.block(2*dim_l1, 0, dim1-2*dim_l1, dim_l2) = GMGolyComplScalar(T, golycompl_basis, tens_family.ancestor(), 2, intmap, 1);
    gm.block(2*dim_l1, dim_l2, dim1-2*dim_l1, dim_l2) = -GMGolyComplScalar(T, golycompl_basis, tens_family.ancestor(), 2, intmap, 0);
    
    return gm;
  }
  
template<typename BasisType1, size_t N>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                                         
                    const TensorizedVectorFamily<BasisType1, N> & tens_family,  
                    const GolyComplBasisCell & golycompl_basis,             
                    MonomialCellIntegralsType mono_int_map = {}                 
                    )
  {
    return GramMatrix(T, golycompl_basis, tens_family, mono_int_map).transpose();
  }
  
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) );
  }
  
Eigen::MatrixXd GMGolyComplScalar(
                    const Cell & T, 
                    const GolyComplBasisCell & golycompl_basis, 
                    const MonomialScalarBasisCell & mono_basis, 
                    const size_t s, 
                    MonomialCellIntegralsType mono_int_map, 
                    const size_t m = 3, 
                    const size_t k1 = 3, 
                    const size_t k2 = 3 
                    );
  
template<typename BasisType>
Eigen::MatrixXd GMGolyComplScalar(
                     const Cell& T, 
                     const GolyComplBasisCell & basis1, 
                     const BasisType & basis2,  
                     const size_t s, 
                     MonomialCellIntegralsType mono_int_map, 
                     const size_t m = 3, 
                     const size_t k1 = 3, 
                     const size_t k2 = 3 
                     )
  {
    // 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', GMGolyComplScalar(T, basis1, basis2.ancestor(), s, mono_int_map, m, k1, k2) );
  }
  
Eigen::MatrixXd GMGolyCompl(
                    const Cell & T, 
                    const GolyComplBasisCell & basis1, 
                    const GolyComplBasisCell & basis2, 
                    const size_t s1, 
                    const size_t s2, 
                    MonomialCellIntegralsType mono_int_map, 
                    const size_t m1 = 3, 
                    const size_t m2 = 3, 
                    const size_t k1 = 3, 
                    const size_t k2 = 3 
                    );
  
  
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<3; 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 = {} 
                     )
  {
    return GramMatrixCurlCurl(T, basis1.ancestor(), basis2.ancestor(), mono_int_map);
  }
  
template<typename BasisType1, typename BasisType2>
typename boost::disable_if<boost::is_same<BasisType2, MonomialCellIntegralsType>, Eigen::MatrixXd>::type GramMatrix(
                     const Cell& T,                          
                     const CurlBasis<BasisType1> & basis1,   
                     const BasisType2 & basis2,              
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    return GramMatrixCurl(T, basis1.ancestor(), basis2, mono_int_map);
  }
  
template<typename BasisType1, typename Basis2>
Eigen::MatrixXd GramMatrix(
                    const Cell & T,                         
                    const BasisType1 & basis1,                  
                    const CurlBasis<Basis2> & basis2,       
                    MonomialCellIntegralsType mono_int_map = {} 
                    )
  {
    return GramMatrixCurl(T, basis2.ancestor(), basis1, mono_int_map).transpose();
  }
  
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrixCurlCurl(
                     const Cell& T, 
                     const TensorizedVectorFamily<BasisType1, N> & basis1, 
                     const TensorizedVectorFamily<BasisType2, N> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree()-2;
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
    
    Eigen::MatrixXd GMxx = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 0, 0, intmap);
    Eigen::MatrixXd GMyy = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 1, 1, intmap);
    Eigen::MatrixXd GMzz = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 2, 2, intmap);
    
    gm.block(0, 0, dim1/N, dim2/N) = GMzz + GMyy;
    gm.block(0, dim2/N, dim1/N, dim2/N) = -GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 1, 0, intmap);
    gm.block(0, 2*dim2/N, dim1/N, dim2/N) = -GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 2, 0, intmap);
    gm.block(dim1/N, 0, dim1/N, dim2/N) = -GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 0, 1, intmap);
    gm.block(dim1/N, dim2/N, dim1/N, dim2/N) = GMzz + GMxx;
    gm.block(dim1/N, 2*dim2/N, dim1/N, dim2/N) = -GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 2, 1, intmap);
    gm.block(2*dim1/N, 0, dim1/N, dim2/N) = -GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 0, 2, intmap);
    gm.block(2*dim1/N, dim2/N, dim1/N, dim2/N) = -GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 1, 2, intmap);
    gm.block(2*dim1/N, 2*dim2/N, dim1/N, dim2/N) = GMyy + GMxx;
    
    return gm/std::pow(T.diam(), 2);
  }
 
Eigen::MatrixXd GramMatrixCurlCurl(
                     const Cell& T, 
                     const GolyComplBasisCell & basis1, 
                     const GolyComplBasisCell & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     );
 
template<typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrixCurlCurl(
                     const Cell& T, 
                     const GolyComplBasisCell & basis1, 
                     const TensorizedVectorFamily<BasisType2, N> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm(dim1, dim2);
    
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree();
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
    
    size_t dim_l1 = basis1.dimPkmo();
    gm.block(0, 0, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 1, 0, 2) - GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 2, 0, 1);
    gm.block(0, dim2/N, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 2, 0, 0) + GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 1, 1, 2)
                                        + GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 2, 2, 2) + 2/T.diam()*GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 3, 3, 2);
    gm.block(0, 2*dim2/N, dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 1, 0, 0) - GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 1, 1, 1)
                                            -GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 2, 2, 1) - 2/T.diam()*GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 3, 3, 1);
    gm.block(dim_l1, 0, dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 2, 1, 1) - GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 0, 0, 2)
                                          -GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 2, 2, 2) - 2/T.diam()*GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 3, 3, 2);
    gm.block(dim_l1, dim2/N, dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 0, 1, 2) + GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 2, 1, 0);
    gm.block(dim_l1, 2*dim2/N, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 0, 1, 1) + GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 0, 0, 0)
                                                +GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 2, 2, 0) + 2/T.diam()*GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 3, 3, 0);
    gm.block(2*dim_l1, 0, dim1-2*dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 1, 2, 2) + GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 0, 0, 1)
                                                   +GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 1, 1, 1) + 2/T.diam()*GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 3, 3, 1);
    gm.block(2*dim_l1, dim2/N, dim1-2*dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 0, 2, 2) - GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 0, 0, 0)
                                                        -GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 1, 1, 0) - 2/T.diam()*GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 3, 3, 0);
    gm.block(2*dim_l1, 2*dim2/N, dim1-2*dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 0, 2, 1) - GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 1, 2, 0);
    
    return gm;
  }
 
template<typename BasisType1, size_t N>
Eigen::MatrixXd GramMatrixCurlCurl(
                     const Cell& T, 
                     const TensorizedVectorFamily<BasisType1, N> & basis1, 
                     const GolyComplBasisCell & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    return GramMatrixCurlCurl(T, basis2, basis1, mono_int_map).transpose();
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixCurlCurl(
                     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', GramMatrixCurlCurl(T, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrixCurlCurl(T, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrixCurlCurl(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
  }
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrixCurl(
                    const Cell & T,                                         
                    const TensorizedVectorFamily<BasisType1, N> & basis1,   
                    const TensorizedVectorFamily<BasisType2, N> & 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);
 
    Eigen::MatrixXd GMx = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 0, intmap);
    Eigen::MatrixXd GMy = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 1, intmap);
    Eigen::MatrixXd GMz = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), 2, intmap);
    
    gm.block(0, dim2/N, dim1/N, dim2/N) = GMz;
    gm.block(0, 2*dim2/N, dim1/N, dim2/N) = -GMy;
    gm.block(dim1/N, 0, dim1/N, dim2/N) = -GMz;
    gm.block(dim1/N, 2*dim2/N, dim1/N, dim2/N) = GMx;
    gm.block(2*dim1/N, 0, dim1/N, dim2/N) = GMy;
    gm.block(2*dim1/N, dim2/N, dim1/N, dim2/N) = -GMx;
    
    return gm/T.diam();
  }
  
template<typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrixCurl(
                    const Cell & T,                                         
                    const GolyComplBasisCell & basis1,                      
                    const TensorizedVectorFamily<BasisType2, N> & basis2,   
                    MonomialCellIntegralsType mono_int_map = {}                 
                    )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm(dim1, dim2);
    
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree();
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
    
    size_t dim_l1 = basis1.dimPkmo();
    gm.block(0, 0, dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 1, 1)
                                     -GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 2, 2)
                                     -2*GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap)/T.diam();
    gm.block(0, dim2/N, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 1, 0);
    gm.block(0, 2*dim2/N, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 0, intmap, 2, 0);
    gm.block(dim_l1, 0, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 0, 1);
    gm.block(dim_l1, dim2/N, dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 0, 0)
                                               -GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 2, 2)
                                               -2*GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap)/T.diam();
    gm.block(dim_l1, 2*dim2/N, dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 1, intmap, 2, 1);
    gm.block(2*dim_l1, 0, dim1-2*dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 0, 2);
    gm.block(2*dim_l1, dim2/N, dim1-2*dim_l1, dim2/N) = GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 1, 2);
    gm.block(2*dim_l1, 2*dim2/N, dim1-2*dim_l1, dim2/N) = -GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 0, 0)
                                                          -GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap, 1, 1)
                                                          -2*GMGolyComplScalar(T, basis1, basis2.ancestor(), 2, intmap)/T.diam();
    
    return gm;
  }
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixCurl(
                     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', GramMatrixCurl(T, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrixCurl(T, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrixCurl(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixCurl(
                     const Cell& T, 
                     const BasisType1 & basis1, 
                     const CurlBasis<BasisType2> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    if constexpr (!useAncestor<BasisType1>()) {
      return GramMatrixCurlCurl(T, basis1, basis2.ancestor(), mono_int_map);
    } else {
      return transformGM(basis1, 'R', GramMatrixCurlCurl(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) );
    }
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(const Cell& T,         
                     const DivergenceBasis<BasisType1> & basis1,   
                     const DivergenceBasis<BasisType2> & basis2,   
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    return GramMatrixDivDiv(T, basis1.ancestor(), basis2.ancestor(), mono_int_map);
  }
 
Eigen::MatrixXd GMRolyComplScalarDiv(
                    const Cell & T, 
                    const MonomialScalarBasisCell & mono_basis, 
                    const RolyComplBasisCell & rolycompl_basis, 
                    const size_t k, 
                    MonomialCellIntegralsType mono_int_map 
                    );
  
template<typename BasisType>
Eigen::MatrixXd GMRolyComplScalarDiv(
                     const Cell& T, 
                     const BasisType & basis1, 
                     const RolyComplBasisCell & basis2, 
                     const size_t k, 
                     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(basis1, 'R', GMRolyComplScalarDiv(T, basis1.ancestor(), basis2, k, mono_int_map) );
  }
  
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrixDivDiv(
                     const Cell& T, 
                     const TensorizedVectorFamily<BasisType1, N> & basis1, 
                     const TensorizedVectorFamily<BasisType2, N> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree()-2;
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
    
    for (size_t i = 0; i < N; i++) {
      for (size_t j = 0; j < N; j++) {
        gm.block(i*dim1/N, j*dim2/N, dim1/N, dim2/N) = GMScalarDerivative(T, basis1.ancestor(), basis2.ancestor(), i, j, intmap);
      }
    }
    
    return gm/std::pow(T.diam(), 2);
  }
 
Eigen::MatrixXd GramMatrixDivDiv(
                     const Cell& T, 
                     const RolyComplBasisCell & basis1, 
                     const RolyComplBasisCell & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     );
 
template<typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrixDivDiv(
                     const Cell& T, 
                     const RolyComplBasisCell & basis1, 
                     const TensorizedVectorFamily<BasisType2, N> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    return GramMatrixDivDiv(T, basis2, basis1, mono_int_map).transpose();
  }
 
template<typename BasisType1, size_t N>
Eigen::MatrixXd GramMatrixDivDiv(
                     const Cell& T, 
                     const TensorizedVectorFamily<BasisType1, N> & basis1, 
                     const RolyComplBasisCell & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm(dim1, dim2);
    
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree();
    MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
    
    for (size_t i = 0; i < N; i++) {
      gm.block(i*dim1/N, 0, dim1/N, dim2) = GMRolyComplScalarDiv(T, basis1.ancestor(), basis2, i, mono_int_map);
    }
    
    return gm;
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixDivDiv(
                     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', GramMatrixDivDiv(T, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrixDivDiv(T, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrixDivDiv(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
  }
 
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, size_t N>
Eigen::MatrixXd GramMatrixDiv(
                    const Cell & T,                                         
                    const TensorizedVectorFamily<BasisType1, N> & 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 < N; i++) {
      gm.block(i*dim1/N, 0, dim1/N, 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) ) );
    }
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixDiv(
                     const Cell& T, 
                     const BasisType1 & basis1, 
                     const DivergenceBasis<BasisType2> & basis2,  
                     MonomialCellIntegralsType mono_int_map = {} 
                     )
  {
    if constexpr (!useAncestor<BasisType1>()) {
      return GramMatrixDivDiv(T, basis1, basis2.ancestor(), mono_int_map);
    } else {
      return transformGM(basis1, 'R', GramMatrixDivDiv(T, basis1.ancestor(), basis2.ancestor(), mono_int_map) );
    }
  }
 
/*------------- Templates for MatrixFamily (including Divergence thereof, and vs Gradient<Tensorized> --------------------*/
 
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 TensorizedVectorFamily<BasisType1, dimspace> & basis1,              
                  const DivergenceBasis<MatrixFamily<BasisType2, dimspace>> & basis2,   
                  MonomialCellIntegralsType mono_int_map = {} 
                  )
  {
    return GramMatrix(T, basis2, basis1, mono_int_map).transpose();
  }
 
 
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 HArDCore3D
 
#endif // end of _GMPOLY_CELL_HPP