basis.hpp

Go to the documentation of this file.

// Core data structures and methods for 2D schemes based on polynomial unknowns in the elements and on the edges
//
// Provides:
//  - Full and partial polynomial spaces on the element and edges
//  - Generic routines to create quadrature nodes over cells and edges of the mesh
//
// Authors: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr)
//          Jerome Droniou (jerome.droniou@monash.edu)
//
 
/*
 *
 *      This library was developed around HHO methods, although some parts of it have a more
 * general purpose. If you use this code or part of it in a scientific publication, 
 * please mention the following book as a reference for the underlying principles
 * of HHO schemes:
 *
 * The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications. 
 *  D. A. Di Pietro and J. Droniou. Modeling, Simulation and Applications, vol. 19. 
 *  Springer International Publishing, 2020, xxxi + 525p. doi: 10.1007/978-3-030-37203-3. 
 *  url: https://hal.archives-ouvertes.fr/hal-02151813.
 *
 */
 
#ifndef BASIS_HPP
#define BASIS_HPP
 
#include <boost/multi_array.hpp>
 
#include <mesh.hpp>
//#include <edge.hpp>
#include <iostream>
 
#include <polynomialspacedimension.hpp>
 
#include <quadraturerule.hpp>
 
namespace HArDCore2D
{
  
  constexpr int dimspace = 2;
  typedef Eigen::Matrix2d MatrixRd;
  typedef Eigen::Vector2d VectorRd;   
  typedef Eigen::Vector2i VectorZd;
  
  static const std::vector<VectorRd> basisRd = { VectorRd(1., 0.), VectorRd(0., 1.) };
  
  template <typename T>
  using BasisQuad = boost::multi_array<T, 2>; 
 
  template <typename T>
  using FType = std::function<T(const VectorRd &)>; 
 
  enum TensorRankE 
    {
      Scalar = 0,
      Vector = 1,
      Matrix = 2
    };
 
 
  //-----------------------------------------------------------------------------
  //                Powers of monomial basis
  //-----------------------------------------------------------------------------
 
  template<typename GeometricSupport>
  struct MonomialPowers
  {
    // Only specializations are relevant
  };
  
  template<>
  struct MonomialPowers<Cell>
  {
    static std::vector<VectorZd> complete(size_t degree)
    {
      std::vector<VectorZd> powers;
      powers.reserve( PolynomialSpaceDimension<Cell>::Poly(degree) );
      for (size_t l = 0; l <= degree; l++)
        {
          for (size_t i = 0; i <= l; i++)
            {
              powers.push_back(VectorZd(i, l - i));
            } // for i
        }   // for l
 
      return powers;
    }
  };
 
  //----------------------------------------------------------------------
  //----------------------------------------------------------------------
  //          STRUCTURE FOR GRADIENT OF MATRIX-VALUED FUNCTIONS
  //----------------------------------------------------------------------
  //----------------------------------------------------------------------
 
  template<size_t N>
  struct MatrixGradient
  {
    MatrixGradient(
      Eigen::Matrix<double, N, N> diff_x,
      Eigen::Matrix<double, N, N> diff_y
      ) : dx(diff_x),
          dy(diff_y)
      {
        // do nothing
      }
    
    MatrixGradient() : dx(Eigen::Matrix<double, N, N>::Zero()),
                       dy(Eigen::Matrix<double, N, N>::Zero())
      {
        // do nothing
      }
                     
    MatrixGradient operator+(const MatrixGradient& G){
      return MatrixGradient<N>(dx + G.dx, dy + G.dy);
    }
 
    MatrixGradient& operator+=(const MatrixGradient& G){
      this->dx += G.dx;
      this->dy += G.dy;
      return *this;
    }
 
    MatrixGradient operator-(const MatrixGradient& G){
      return MatrixGradient<N>(dx - G.dx, dy - G.dy);
    }
 
    // Directional derivatives
    Eigen::Matrix<double, N, N> dx;
    Eigen::Matrix<double, N, N> dy;
  
  };
 
  template <size_t N>
  MatrixGradient<N> operator*(double scalar, MatrixGradient<N> const & G){
    return MatrixGradient<N>(scalar * G.dx, scalar * G.dy);
  };
 
  //----------------------------------------------------------------------
  //----------------------------------------------------------------------
  //          SCALAR MONOMIAL BASIS
  //----------------------------------------------------------------------
  //----------------------------------------------------------------------
 
  class MonomialScalarBasisCell
  {
  public:
    typedef double FunctionValue;
    typedef VectorRd GradientValue;
    typedef VectorRd CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef MatrixRd HessianValue;
 
    typedef Cell GeometricSupport;    
 
    constexpr static const TensorRankE tensorRank = Scalar;
    constexpr static const bool hasAncestor = false;
    static const bool hasFunction = true;
    static const bool hasGradient = true;
    static const bool hasCurl = true;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = true;
    static const bool hasCurlCurl = false;
    
    MonomialScalarBasisCell(
                            const Cell &T, 
                            size_t degree  
                            );
 
    inline size_t dimension() const
    {
      return ( m_degree >= 0 ? (m_degree + 1) * (m_degree + 2) / 2 : 0);
    }
 
    FunctionValue function(size_t i, const VectorRd & x) const;
 
    GradientValue gradient(size_t i, const VectorRd & x) const;
 
    CurlValue curl(size_t i, const VectorRd & x) const;
 
    HessianValue hessian(size_t i, const VectorRd & x) const;
    
    inline size_t max_degree() const
    {
      return m_degree;
    }
 
    inline VectorZd powers(size_t i) const
    {
      return m_powers[i];
    }
 
  private:
    inline VectorRd _coordinate_transform(const VectorRd& x) const
    {
      return (x - m_xT) / m_hT;
    }
 
    size_t m_degree;
    VectorRd m_xT;
    double m_hT;
    std::vector<VectorZd> m_powers;
    Eigen::Matrix2d m_rot; // rotation pi/2
  };
 
  //------------------------------------------------------------------------------
  
  class MonomialScalarBasisEdge
  {
  public:
    typedef double FunctionValue;
    typedef VectorRd GradientValue;
    typedef VectorRd CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
    
    typedef Edge GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Scalar;
    constexpr static const bool hasAncestor = false;
    static const bool hasFunction = true;
    static const bool hasGradient = true;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
 
    MonomialScalarBasisEdge(
                            const Edge &E, 
                            size_t degree  
                            );
 
    inline size_t dimension() const
    {
      return m_degree + 1;
    }
 
    FunctionValue function(size_t i, const VectorRd &x) const;
 
    GradientValue gradient(size_t i, const VectorRd &x) const;
 
    inline size_t max_degree() const
    {
      return m_degree;
    }
 
  private:
    inline double _coordinate_transform(const VectorRd &x) const
    {
      return (x - m_xE).dot(m_tE) / m_hE;
    }
 
    size_t m_degree;
    VectorRd m_xE;
    double m_hE;
    VectorRd m_tE;
  };
 
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  //          DERIVED BASIS
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  //----------------------------FAMILY------------------------------------------
  //
  template <typename BasisType>
  class Family
  {
  public:
    typedef typename BasisType::FunctionValue FunctionValue;
    typedef typename BasisType::GradientValue GradientValue;
    typedef typename BasisType::CurlValue CurlValue;
    typedef typename BasisType::DivergenceValue DivergenceValue;
    typedef typename BasisType::RotorValue RotorValue;
    typedef typename BasisType::HessianValue HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = BasisType::tensorRank;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = BasisType::hasFunction;
    static const bool hasGradient = BasisType::hasGradient;
    static const bool hasCurl = BasisType::hasCurl;
    static const bool hasDivergence = BasisType::hasDivergence;
    static const bool hasRotor = BasisType::hasRotor;
    static const bool hasHessian = BasisType::hasHessian;
    static const bool hasCurlCurl = BasisType::hasCurlCurl;
 
    typedef BasisType AncestorType;
    
    Family(
           const BasisType &basis,       
           const Eigen::MatrixXd &matrix 
           )
      : m_basis(basis),
        m_matrix(matrix)
    {
      assert((size_t)matrix.cols() == basis.dimension() || "Inconsistent family initialization");
    }
 
    inline size_t dimension() const
    {
      return m_matrix.rows();
    }
 
    FunctionValue function(size_t i, const VectorRd &x) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      FunctionValue f = m_matrix(i, 0) * m_basis.function(0, x);
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          f += m_matrix(i, j) * m_basis.function(j, x);
        } // for j
      return f;
    }
 
    FunctionValue function(size_t i, size_t iqn, const boost::multi_array<FunctionValue, 2> &ancestor_value_quad) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      FunctionValue f = m_matrix(i, 0) * ancestor_value_quad[0][iqn];
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          f += m_matrix(i, j) * ancestor_value_quad[j][iqn];
        } // for j
      return f;
    }
   
    
    GradientValue gradient(size_t i, const VectorRd &x) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      GradientValue G = m_matrix(i, 0) * m_basis.gradient(0, x);
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          G += m_matrix(i, j) * m_basis.gradient(j, x);
        } // for j
      return G;
    }
 
    GradientValue gradient(size_t i, size_t iqn, const boost::multi_array<GradientValue, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      GradientValue G = m_matrix(i, 0) * ancestor_gradient_quad[0][iqn];
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          G += m_matrix(i, j) * ancestor_gradient_quad[j][iqn];
        } // for j
      return G;
    }
 
    CurlValue curl(size_t i, const VectorRd &x) const
    {
      static_assert(hasCurl, "Call to curl() not available");
 
      CurlValue C = m_matrix(i, 0) * m_basis.curl(0, x);
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          C += m_matrix(i, j) * m_basis.curl(j, x);
        } // for j
      return C;
    }
 
    CurlValue curl(size_t i, size_t iqn, const boost::multi_array<CurlValue, 2> &ancestor_curl_quad) const
    {
      static_assert(hasCurl, "Call to curl() not available");
 
      CurlValue C = m_matrix(i, 0) * ancestor_curl_quad[0][iqn];
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          C += m_matrix(i, j) * ancestor_curl_quad[j][iqn];
        } // for j
      return C;
    }
 
    DivergenceValue divergence(size_t i, const VectorRd &x) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      DivergenceValue D = m_matrix(i, 0) * m_basis.divergence(0, x);
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          D += m_matrix(i, j) * m_basis.divergence(j, x);
        } // for j
      return D;
    }
 
    DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array<DivergenceValue, 2> &ancestor_divergence_quad) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      DivergenceValue D = m_matrix(i, 0) * ancestor_divergence_quad[0][iqn];
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          D += m_matrix(i, j) * ancestor_divergence_quad[j][iqn];
        } // for j
      return D;
    }
 
    DivergenceValue rotor(size_t i, const VectorRd &x) const
    {
      static_assert(hasRotor, "Call to rotor() not available");
 
      RotorValue R = m_matrix(i, 0) * m_basis.rotor(0, x);
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          R += m_matrix(i, j) * m_basis.rotor(j, x);
        } // for j
      return R;
    }
 
    RotorValue rotor(size_t i, size_t iqn, const boost::multi_array<RotorValue, 2> &ancestor_rotor_quad) const
    {
      static_assert(hasRotor, "Call to rotor() not available");
 
      RotorValue R = m_matrix(i, 0) * ancestor_rotor_quad[0][iqn];
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          R += m_matrix(i, j) * ancestor_rotor_quad[j][iqn];
        } // for j
      return R;
    }
 
    inline HessianValue hessian(size_t i, const VectorRd & x) const
    {
      static_assert(hasHessian, "Call to hessian() not available");
 
      HessianValue H = m_matrix(i, 0) * m_basis.hessian(0, x);
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          H += m_matrix(i, j) * m_basis.hessian(j, x);
        } // for j
      return H;
    }
 
    HessianValue hessian(size_t i, size_t iqn, const boost::multi_array<HessianValue, 2> &ancestor_hessian_quad) const
    {
      static_assert(hasHessian, "Call to hessian() not available");
 
      HessianValue H = m_matrix(i, 0) * ancestor_hessian_quad[0][iqn];
      for (auto j = 1; j < m_matrix.cols(); j++)
        {
          H += m_matrix(i, j) * ancestor_hessian_quad[j][iqn];
        } // for j
      return H;
    }
 
    inline const Eigen::MatrixXd & matrix() const
    {
      return m_matrix;
    }
 
    constexpr inline const BasisType &ancestor() const
    {
      return m_basis;
    }
 
    inline size_t max_degree() const
    {
      return m_basis.max_degree();
    }
 
  private:
    BasisType m_basis;
    Eigen::MatrixXd m_matrix;
  };
 
  //----------------------TENSORIZED--------------------------------------------------------
 
 
  template <typename ScalarFamilyType, size_t N>
  class TensorizedVectorFamily
  {
  public:
    typedef typename Eigen::Matrix<double, N, 1> FunctionValue;
    typedef typename Eigen::Matrix<double, N, dimspace> GradientValue;
    typedef GradientValue CurlValue;
    typedef VectorRd CurlCurlValue;
    typedef double DivergenceValue;
    typedef double RotorValue;
    typedef void HessianValue;    
 
    typedef typename ScalarFamilyType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Vector;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = ScalarFamilyType::hasFunction;
    static const bool hasGradient = ScalarFamilyType::hasGradient;
    // We know how to compute the curl, scalar rotor, and divergence if gradient is available
    static const bool hasDivergence = ( ScalarFamilyType::hasGradient && N==dimspace );
    static const bool hasRotor = ( ScalarFamilyType::hasGradient && N==dimspace );
    static const bool hasCurl = ( ScalarFamilyType::hasGradient && N==dimspace );
    // In future versions, one could include the computation of second derivatives
    // checking that BasisType has information on the Hessian
    static const bool hasHessian = false;
    static const bool hasCurlCurl = (ScalarFamilyType::hasHessian && N == dimspace);
   
    typedef ScalarFamilyType AncestorType;
       
    TensorizedVectorFamily(const ScalarFamilyType & scalar_family)
      : m_scalar_family(scalar_family)
    {
      static_assert(ScalarFamilyType::tensorRank == Scalar,
                    "Vector family can only be constructed from scalar families");
      m_rot.row(0) << 0., 1.;
      m_rot.row(1) << -1., 0.;
 
    }
 
    inline size_t dimension() const
    {
      return m_scalar_family.dimension() * N;
    }
 
    FunctionValue function(size_t i, const VectorRd &x) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      FunctionValue ek = Eigen::Matrix<double, N, 1>::Zero();
      ek(i / m_scalar_family.dimension()) = 1.;
      return ek * m_scalar_family.function(i % m_scalar_family.dimension(), x);
    }
 
    FunctionValue function(size_t i, size_t iqn, const boost::multi_array<double, 2> &ancestor_value_quad) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      FunctionValue ek = Eigen::Matrix<double, N, 1>::Zero();
      ek(i / m_scalar_family.dimension()) = 1.;
      return ek * ancestor_value_quad[i % m_scalar_family.dimension()][iqn];
    }
 
    GradientValue gradient(size_t i, const VectorRd &x) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      GradientValue G = Eigen::Matrix<double, N, dimspace>::Zero();
      G.row(i / m_scalar_family.dimension()) = m_scalar_family.gradient(i % m_scalar_family.dimension(), x);
      return G;
    }
 
    GradientValue gradient(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      GradientValue G = Eigen::Matrix<double, N, dimspace>::Zero();
      G.row(i / m_scalar_family.dimension()) = ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn];
      return G;
    }
 
    CurlValue curl(size_t i, const VectorRd &x) const
    {
      static_assert(hasCurl, "Call to curl() not available");
      
      return m_rot * gradient(i, x);
    }
 
    CurlValue curl(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasCurl, "Call to curl() not available");
 
      return m_rot * ancestor_gradient_quad[i][iqn];
    }
 
    DivergenceValue divergence(size_t i, const VectorRd &x) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      return m_scalar_family.gradient(i % m_scalar_family.dimension(), x)(i / m_scalar_family.dimension());
    }
 
    DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      return ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn](i / m_scalar_family.dimension());      
    }
 
    RotorValue rotor(size_t i, const VectorRd &x) const
    {
      static_assert(hasRotor && hasCurl, "Call to rotor() not available");
 
      return curl(i, x).trace();
    }
    
    RotorValue rotor(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasRotor, "Call to rotor() not available");
 
      return curl(i, iqn, ancestor_gradient_quad).trace();
    }
    
    CurlCurlValue curlcurl(size_t i, const VectorRd &x) const
    {
      static_assert(hasCurlCurl, "Call to curlcurl() not available");
 
      FunctionValue ek = Eigen::Matrix<double, N, 1>::Zero();
      int p = i / m_scalar_family.dimension();
      ek(p) = 1.;
      typename ScalarFamilyType::HessianValue hess = m_scalar_family.hessian(i % m_scalar_family.dimension(), x);
      return - hess.trace() * ek + hess.col(p);
    }
 
    CurlCurlValue curlcurl(size_t i, size_t iqn, const boost::multi_array<typename ScalarFamilyType::HessianValue, 2> &ancestor_hessian_quad) const
    {
      static_assert(hasCurlCurl, "Call to curlcurl() not available");
 
      FunctionValue ek = Eigen::Matrix<double, N, 1>::Zero();
      int p = i / m_scalar_family.dimension();
      ek(p) = 1.;
      typename ScalarFamilyType::HessianValue hess = ancestor_hessian_quad[i % m_scalar_family.dimension()][iqn];
      return - hess.trace() * ek + hess.col(p);
    }
    
    constexpr inline const ScalarFamilyType &ancestor() const
    {
      return m_scalar_family;
    }
 
    inline size_t max_degree() const
    {
      return m_scalar_family.max_degree();
    }
    
  private:
    ScalarFamilyType m_scalar_family;
    Eigen::Matrix2d m_rot;  // Rotation of pi/2
  };
  
  //----------------------MATRIX FAMILY--------------------------------------------------------
 
 
  // Useful formulas to manipulate this basis (all indices starting from 0):
  //      the i-th element in the basis is f_{i%r} E_{i/r}
  //      the matrix E_m has 1 at position (m%N, m/N)
  //      the element f_aE_b is the i-th in the family with i = b*r + a
  template<typename ScalarFamilyType, size_t N>
  class MatrixFamily
  {
  public:
    typedef typename Eigen::Matrix<double, N, N> FunctionValue;
    typedef MatrixGradient<N> GradientValue;
    typedef VectorRd CurlValue;
    typedef typename Eigen::Matrix<double, N, 1> DivergenceValue;
    typedef void RotorValue;    
    typedef void HessianValue;    
 
    typedef typename ScalarFamilyType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Matrix;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = ScalarFamilyType::hasFunction;
    static const bool hasGradient = true;
    static const bool hasDivergence = ( ScalarFamilyType::hasGradient && N==dimspace );
    static const bool hasRotor = false;
    static const bool hasCurl = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
   
    typedef ScalarFamilyType AncestorType;
 
    MatrixFamily(const ScalarFamilyType & scalar_family)
      : m_scalar_family(scalar_family),
        m_E(N*N),
        m_transposeOperator(Eigen::MatrixXd::Zero(scalar_family.dimension()*N*N, scalar_family.dimension()*N*N))
    {
      static_assert(ScalarFamilyType::tensorRank == Scalar,
                    "Vector family can only be constructed from scalar families");
                    
      // Construct the basis for NxN matrices
      for (size_t j = 0; j < N; j++){
        for (size_t i = 0; i < N; i++){
          m_E[j*N + i] = Eigen::Matrix<double, N, N>::Zero();
          m_E[j*N + i](i,j) = 1.;
        }        
      }
      
      // Construct the transpose operator
      for (size_t i = 0; i < dimension(); i++){
        // The i-th basis element is sent to the j-th basis element s.t. j%r=i%r, (j/r)%N=(i/r)/N, (j/r)/N=(i/r)%N
        size_t r = m_scalar_family.dimension();
        size_t j = ( ( (i/r)%N )*N + (i/r)/N ) * r + (i%r);
        m_transposeOperator(i,j) = 1.;
      }
    }
 
    inline size_t dimension() const
    {
      return m_scalar_family.dimension() * N * N;
    }
 
    FunctionValue function(size_t i, const VectorRd & x) const
    {
      static_assert(hasFunction, "Call to function() not available");
            
      return m_scalar_family.function(i % m_scalar_family.dimension(), x) * m_E[i / m_scalar_family.dimension()];
    }
    
    FunctionValue function(size_t i, size_t iqn, const boost::multi_array<double, 2> &ancestor_value_quad) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      return ancestor_value_quad[i % m_scalar_family.dimension()][iqn] * m_E[i / m_scalar_family.dimension()];
    }
 
    DivergenceValue divergence(size_t i, const VectorRd & x) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      Eigen::Matrix<double, N, 1> V = m_E[i/m_scalar_family.dimension()].col( (i/m_scalar_family.dimension()) /N);
      return m_scalar_family.gradient(i % m_scalar_family.dimension(), x)( (i / m_scalar_family.dimension()) / N) * V;
    }
    
    DivergenceValue divergence(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      Eigen::Matrix<double, N, 1> V = m_E[i/m_scalar_family.dimension()].col( (i/m_scalar_family.dimension()) /N);
      return ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn]( (i / m_scalar_family.dimension()) / N) * V;      
    }
 
    GradientValue gradient(size_t i, const VectorRd & x) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      const VectorRd gradient_scalar_function = m_scalar_family.gradient(i % m_scalar_family.dimension(), x);
      const Eigen::Matrix<double, N, N> matrix_basis = m_E[i / m_scalar_family.dimension()];
      return GradientValue(gradient_scalar_function.x() * matrix_basis, gradient_scalar_function.y() * matrix_basis); 
    }
 
    GradientValue gradient(size_t i, size_t iqn, const boost::multi_array<VectorRd, 2> &ancestor_gradient_quad) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      const VectorRd gradient_scalar_function = ancestor_gradient_quad[i % m_scalar_family.dimension()][iqn];
      const Eigen::Matrix<double, N, N> matrix_basis = m_E[i / m_scalar_family.dimension()];
      return GradientValue(gradient_scalar_function.x() * matrix_basis, gradient_scalar_function.y() * matrix_basis); 
    }
 
    inline const ScalarFamilyType &ancestor() const
    {
      return m_scalar_family;
    }
    
    inline const size_t matrixSize() const
    {
      return N;
    }
    
    inline const Eigen::MatrixXd transposeOperator() const
    {
      return m_transposeOperator;
    }
    
    inline const Eigen::MatrixXd symmetriseOperator() const
    {
      return (Eigen::MatrixXd::Identity(dimension(), dimension()) + m_transposeOperator)/2;
    }
 
    const Eigen::MatrixXd symmetricBasis() const
    {
      size_t dim_scalar =  m_scalar_family.dimension();
      Eigen::MatrixXd symOpe = symmetriseOperator();
      
      Eigen::MatrixXd SB = Eigen::MatrixXd::Zero(dim_scalar * N*(N+1)/2, dimension());
      
      // The matrix consists in selecting certain rows of symmetriseOperator, corresponding to the lower triangular part
      // To select these rows, we loop over the columns of the underlying matrices, and get the rows in symmetriseOperator()
      // corresponding to the lower triangular part of each column
      size_t position = 0;
      for (size_t i = 0; i<N; i++){
        // Skip the first i*dim_scalar elements in the column (upper triangle) in symOpe, take the remaining (N-i)*dim_scalar
        SB.middleRows(position, (N-i)*dim_scalar) = symOpe.middleRows(i*(N+1)*dim_scalar, (N-i)*dim_scalar);
        // update current position
        position += (N-i)*dim_scalar;
      }
      
      return SB;
    }
 
 
    inline const Eigen::MatrixXd traceOperator() const
    {
      size_t r =  m_scalar_family.dimension();
      
      Eigen::MatrixXd Tr = Eigen::MatrixXd::Zero(r, dimension());
      
      /* According to the formulas described at the start of this class, the l-th element of the
       matrix family is f_{l%r} E_{l/r}.
       Given the construction of E_m=E_{l/r}, the trace is either 0, or 1 if m%N=m/N, that is,
       (l/r)%N = (l/r)/N.
       So Tr_{kl}=1 only if k=l%r and (l/r)%N = (l/r)/N
      */
      for (size_t l=0; l < dimension(); l++){
        if ( size_t(l/r)%N == size_t(size_t(l/r)/N) ){
          size_t k = l%r;
          Tr(k,l) = 1.;
        }
      }
             
      return Tr;
    }
    
  private:
    ScalarFamilyType m_scalar_family;
    std::vector<Eigen::Matrix<double, N, N>> m_E;
    Eigen::MatrixXd m_transposeOperator;
  };
 
  
  //----------------------TANGENT FAMILY--------------------------------------------------------
 
  template <typename ScalarFamilyType>
  class TangentFamily
  {
  public:
    typedef VectorRd FunctionValue;
    typedef VectorRd GradientValue;
    typedef VectorRd CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;    
 
    typedef Edge GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Vector;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = ScalarFamilyType::hasHessian;
 
    typedef ScalarFamilyType AncestorType;
    
    TangentFamily(
                  const ScalarFamilyType &scalar_family,    
                  const VectorRd &generator                   
                  )
      : m_scalar_family(scalar_family),
        m_generator(generator)
    {
      static_assert(ScalarFamilyType::hasFunction, "Call to function() not available");
      static_assert(std::is_same<typename ScalarFamilyType::GeometricSupport, Edge>::value,
                    "Tangent families can only be defined on edges");
    }
 
    inline size_t dimension() const
    {
      return m_scalar_family.dimension();
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      return m_scalar_family.function(i, x) * m_generator;
    }
 
    constexpr inline const ScalarFamilyType &ancestor() const
    {
      return m_scalar_family;
    }
 
    inline size_t max_degree() const
    {
      return m_scalar_family.max_degree();
    }
 
    inline VectorRd generator() const
    {
      return m_generator;
    }
 
  private:
    ScalarFamilyType m_scalar_family;
    VectorRd m_generator;
  };
 
 
  //--------------------SHIFTED BASIS----------------------------------------------------------
 
 
  template <typename BasisType>
  class ShiftedBasis
  {
  public:
    typedef typename BasisType::FunctionValue FunctionValue;
    typedef typename BasisType::GradientValue GradientValue;
    typedef typename BasisType::CurlValue CurlValue;
    typedef typename BasisType::DivergenceValue DivergenceValue;
    typedef typename BasisType::RotorValue RotorValue;
    typedef typename BasisType::HessianValue HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = BasisType::tensorRank;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = BasisType::hasFunction;
    static const bool hasGradient = BasisType::hasGradient;
    static const bool hasCurl = BasisType::hasCurl;
    static const bool hasDivergence = BasisType::hasDivergence;
    static const bool hasRotor = BasisType::hasRotor;
    static const bool hasHessian = BasisType::hasHessian;
    static const bool hasCurlCurl = BasisType::hasCurlCurl;
 
    typedef BasisType AncestorType;
    
    ShiftedBasis(
                 const BasisType &basis, 
                 const int shift         
                 )
      : m_basis(basis),
        m_shift(shift)
    {
      // Do nothing
    }
 
    inline size_t dimension() const
    {
      return m_basis.dimension() - m_shift;
    }
    
    constexpr inline const BasisType &ancestor() const
    {
      return m_basis;
    }
 
    inline size_t max_degree() const
    {
      return m_basis.max_degree();
    }
    
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      return m_basis.function(i + m_shift, x);
    }
 
    inline GradientValue gradient(size_t i, const VectorRd &x) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      return m_basis.gradient(i + m_shift, x);
    }
 
    inline CurlValue curl(size_t i, const VectorRd &x) const
    {
      static_assert(hasCurl, "Call to curl() not available");
 
      return m_basis.curl(i + m_shift, x);
    }
 
    inline DivergenceValue divergence(size_t i, const VectorRd &x) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      return m_basis.divergence(i + m_shift, x);
    }
 
    inline RotorValue rotor(size_t i, const VectorRd &x) const
    {
      static_assert(hasRotor, "Call to rotor() not available");
 
      return m_basis.rotor(i + m_shift, x);
    }
 
    inline HessianValue hessian(size_t i, const VectorRd & x) const
    {
      static_assert(hasHessian, "Call to hessian() not available");
 
      return m_basis.hessian(i + m_shift, x);
    }
    
  private:
    BasisType m_basis;
    int m_shift;
  };
 
  //-----------------------RESTRICTED BASIS-------------------------------------------------------
 
 
  template <typename BasisType>
  class RestrictedBasis
  {
  public:
    typedef typename BasisType::FunctionValue FunctionValue;
    typedef typename BasisType::GradientValue GradientValue;
    typedef typename BasisType::CurlValue CurlValue;
    typedef typename BasisType::DivergenceValue DivergenceValue;
    typedef typename BasisType::RotorValue RotorValue;
    typedef typename BasisType::HessianValue HessianValue;
    
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = BasisType::tensorRank;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = BasisType::hasFunction;
    static const bool hasGradient = BasisType::hasGradient;
    static const bool hasCurl = BasisType::hasCurl;
    static const bool hasDivergence = BasisType::hasDivergence;
    static const bool hasRotor = BasisType::hasRotor;
    static const bool hasHessian = BasisType::hasHessian;
    static const bool hasCurlCurl = BasisType::hasCurlCurl;
    
    typedef BasisType AncestorType;
    
    RestrictedBasis(
                    const BasisType &basis, 
                    const size_t &dimension 
                    )
      : m_basis(basis),
        m_dimension(dimension)
    {
      // Make sure that the provided dimension is smaller than the one
      // of the provided basis
      assert(dimension <= basis.dimension());
    }
 
    inline size_t dimension() const
    {
      return m_dimension;
    }
    
    constexpr inline const BasisType &ancestor() const
    {
      return m_basis;
    }
    
    inline size_t max_degree() const
    {
      // We need to find the degree based on the dimension, assumes the basis is hierarchical
      size_t deg = 0;
      while (PolynomialSpaceDimension<GeometricSupport>::Poly(deg) < m_dimension){
        deg++;
      }
      return deg;
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      static_assert(hasFunction, "Call to function() not available");
 
      return m_basis.function(i, x);
    }
 
    GradientValue gradient(size_t i, const VectorRd &x) const
    {
      static_assert(hasGradient, "Call to gradient() not available");
 
      return m_basis.gradient(i, x);
    }
 
    CurlValue curl(size_t i, const VectorRd &x) const
    {
      static_assert(hasCurl, "Call to curl() not available");
 
      return m_basis.curl(i, x);
    }
 
    DivergenceValue divergence(size_t i, const VectorRd &x) const
    {
      static_assert(hasDivergence, "Call to divergence() not available");
 
      return m_basis.divergence(i, x);
    }
 
    RotorValue rotor(size_t i, const VectorRd &x) const
    {
      static_assert(hasRotor, "Call to rotor() not available");
 
      return m_basis.rotor(i, x);
    }
 
  private:
    BasisType m_basis;
    size_t m_dimension;
  };
 
  //--------------------GRADIENT BASIS----------------------------------------------------------
 
 
  template <typename BasisType>
  class GradientBasis
  {
  public:
    typedef typename BasisType::GradientValue FunctionValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
    typedef VectorRd CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Vector;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = true;
    // In future versions, one could include the computation of derivatives
    // checking that BasisType has information on the Hessian
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
    
    typedef BasisType AncestorType;
    
    GradientBasis(const BasisType &basis)
      : m_scalar_basis(basis)
    {
      static_assert(BasisType::hasGradient,
                    "Gradient basis requires gradient() for the original basis to be available");
      // Do nothing
    }
 
    inline size_t dimension() const
    {
      return m_scalar_basis.dimension();
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      return m_scalar_basis.gradient(i, x);
    }
 
    constexpr inline const BasisType &ancestor() const
    {
      return m_scalar_basis;
    }
 
 
  private:
    BasisType m_scalar_basis;
  };
 
  //-------------------CURL BASIS-----------------------------------------------------------
 
 
  template <typename BasisType>
  class CurlBasis
  {
  public:
    typedef VectorRd FunctionValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Vector;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
 
    typedef BasisType AncestorType;
    
    CurlBasis(const BasisType &basis)
      : m_basis(basis)
    {
      static_assert(BasisType::tensorRank == Scalar && std::is_same<typename BasisType::GeometricSupport, Cell>::value,
                    "Curl basis can only be constructed starting from scalar bases on elements");
      static_assert(BasisType::hasCurl,
                    "Curl basis requires curl() for the original basis to be available");
    }
 
    inline size_t dimension() const
    {
      return m_basis.dimension();
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      return m_basis.curl(i, x);
    }
    
    constexpr inline const BasisType &ancestor() const
    {
      return m_basis;
    }
 
 
  private:
    size_t m_degree;
    BasisType m_basis;
  };
 
 
  //--------------------DIVERGENCE BASIS----------------------------------------------------------
 
 
  template <typename BasisType>
  class DivergenceBasis
  {
  public:
    typedef double FunctionValue;
    typedef VectorRd GradientValue;
    typedef VectorRd CurlValue;
    typedef double DivergenceValue;
    typedef void HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Scalar;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
    
    typedef BasisType AncestorType;
    
    DivergenceBasis(const BasisType &basis)
      : m_vector_basis(basis)
    {
      static_assert(BasisType::tensorRank == Vector || BasisType::tensorRank == Matrix,
                    "Divergence basis can only be constructed starting from vector bases");
      static_assert(BasisType::hasDivergence,
                    "Divergence basis requires divergence() for the original basis to be available");
      // Do nothing
    }
 
    inline size_t dimension() const
    {
      return m_vector_basis.dimension();
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      return m_vector_basis.divergence(i, x);
    }
 
    constexpr inline const BasisType &ancestor() const
    {
      return m_vector_basis;
    }
 
  private:
    BasisType m_vector_basis;
  };
 
 
  //--------------------ROTOR BASIS----------------------------------------------------------
 
 
  template <typename BasisType>
  class RotorBasis
  {
  public:
    typedef double FunctionValue;
    typedef VectorRd GradientValue;
    typedef VectorRd CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Scalar;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
    
    typedef BasisType AncestorType;
    
    RotorBasis(const BasisType &basis)
      : m_vector_basis(basis)
    {
      static_assert(BasisType::tensorRank == Vector,
                    "Rotor basis can only be constructed starting from vector bases");
      static_assert(BasisType::hasRotor,
                    "Rotor basis requires rotor() for the original basis to be available");
      // Do nothing
    }
 
    inline size_t dimension() const
    {
      return m_vector_basis.dimension();
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      return m_vector_basis.rotor(i, x);
    }
 
    constexpr inline const BasisType &ancestor() const
    {
      return m_vector_basis;
    }
 
  private:
    BasisType m_vector_basis;
  };
 
 
  //--------------------HESSIAN BASIS----------------------------------------------------------
 
 
  template <typename BasisType>
  class HessianBasis
  {
  public:
    typedef MatrixRd FunctionValue;
    // The following types do not make sense in this context and are set to void
    typedef void GradientValue;
    typedef void CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
 
    typedef typename BasisType::GeometricSupport GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Matrix;
    constexpr static const bool hasAncestor = true;
    static const bool hasFunction = true;
    // In future versions, one could include the computation of derivatives
    // checking that BasisType has information on the Hessian
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
 
    typedef BasisType AncestorType;
    
    HessianBasis(const BasisType &basis)
      : m_scalar_basis(basis)
    {
      static_assert(BasisType::tensorRank == Scalar,
                    "Hessian basis can only be constructed starting from scalar bases");
      static_assert(BasisType::hasHessian,
                    "Hessian basis requires hessian() for the original basis to be available");
      // Do nothing
    }
 
    inline size_t dimension() const
    {
      return m_scalar_basis.dimension();
    }
 
    inline FunctionValue function(size_t i, const VectorRd &x) const
    {
      return m_scalar_basis.hessian(i, x);
    }
 
    constexpr inline const BasisType &ancestor() const
    {
      return m_scalar_basis;
    }
 
 
  private:
    BasisType m_scalar_basis;
  };
 
  //---------------------------------------------------------------------
  //      BASES FOR THE KOSZUL COMPLEMENTS OF G^k, R^k, H^k
  //---------------------------------------------------------------------
  class RolyComplBasisCell
  {
  public:
    typedef VectorRd FunctionValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> CurlValue;
    typedef double DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
 
    typedef Cell GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Vector;
    constexpr static const bool hasAncestor = false;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = true;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
 
    RolyComplBasisCell(
                       const Cell &T, 
                       size_t degree  
                       );
 
    inline size_t dimension() const
    {
      return PolynomialSpaceDimension<Cell>::RolyCompl(m_degree);
    }
 
    FunctionValue function(size_t i, const VectorRd &x) const;
 
    DivergenceValue divergence(size_t i, const VectorRd &x) const;
    
    inline size_t max_degree() const
    {
      return m_degree;
    }
    
    inline VectorZd powers(size_t i) const
    {
      return m_powers[i];
    }
    
  private:
    inline VectorRd _coordinate_transform(const VectorRd &x) const
    {
      return (x - m_xT) / m_hT;
    }
 
    size_t m_degree;
    VectorRd m_xT;
    double m_hT;
    std::vector<VectorZd> m_powers;
  };
 
 
  class GolyComplBasisCell
  {
  public:
    typedef VectorRd FunctionValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> GradientValue;
    typedef Eigen::Matrix<double, dimspace, dimspace> CurlValue;
    typedef void DivergenceValue;
    typedef double RotorValue;
    typedef void HessianValue;
 
    typedef Cell GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Vector;
    constexpr static const bool hasAncestor = false;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = false;
    static const bool hasRotor = true;
    static const bool hasHessian = false;
    static const bool hasCurlCurl = false;
 
 
    GolyComplBasisCell(
                       const Cell &T, 
                       size_t degree  
                       );
 
    inline size_t dimension() const
    {
      return PolynomialSpaceDimension<Cell>::GolyCompl(m_degree);
    }
 
    FunctionValue function(size_t i, const VectorRd &x) const;
 
    RotorValue rotor(size_t i, const VectorRd &x) const;
 
    inline const std::shared_ptr<RolyComplBasisCell> &rck() const
    {
      return m_Rck_basis;
    }
 
  private:
    size_t m_degree;
    Eigen::Matrix2d m_rot;  // Rotation of pi/2: the basis of Gck is the rotated basis of Rck
    std::shared_ptr<RolyComplBasisCell> m_Rck_basis;  // A basis of Rck, whose rotation gives a basis of Gck
  };
 
  class HolyComplBasisCell
  {
  public:
    typedef MatrixRd FunctionValue;
    // The following types do not make sense in this context and are set to void
    typedef void GradientValue;
    typedef void CurlValue;
    typedef void DivergenceValue;
    typedef void RotorValue;
    typedef void HessianValue;
 
    typedef Cell GeometricSupport;
 
    constexpr static const TensorRankE tensorRank = Matrix;
    constexpr static const bool hasAncestor = false;
    static const bool hasFunction = true;
    static const bool hasGradient = false;
    static const bool hasCurl = false;
    static const bool hasDivergence = true;
    static const bool hasRotor = false;
    static const bool hasHessian = false;
 
    HolyComplBasisCell(
                       const Cell &T, 
                       size_t degree  
                       );
 
    inline size_t dimension() const
    {
      return PolynomialSpaceDimension<Cell>::HolyCompl(m_degree);
    }
 
    FunctionValue function(size_t i, const VectorRd &x) const;
 
    inline size_t max_degree() const
    {
      return m_degree;
    }    
    
  private:
    inline VectorRd _coordinate_transform(const VectorRd &x) const
    {
      return (x - m_xT) / m_hT;
    }
 
    size_t m_degree;
    VectorRd m_xT;
    double m_hT;
    Eigen::Matrix2d m_rot;
    TensorizedVectorFamily<MonomialScalarBasisCell, dimspace> m_Pkmo;
  };
 
  //------------------------------------------------------------------------------
  //            Free functions
  //------------------------------------------------------------------------------
 
  enum BasisFunctionE
    {
      Function,
      Gradient,
      SymmetricGradient, //AC
      SkewsymmetricGradient, //AC
      Curl,
      Divergence,
      Rotor,
      Hessian,
      CurlCurl,
    };
  
  template <typename outValue, typename inValue, typename FunctionType>
  boost::multi_array<outValue, 2> transform_values_quad(
                                                        const boost::multi_array<inValue, 2> &B_quad,   
                                                        const FunctionType &F                          
                                                        )
  {
    boost::multi_array<outValue, 2> transformed_B_quad( boost::extents[B_quad.shape()[0]][B_quad.shape()[1]] );
    
    std::transform( B_quad.origin(), B_quad.origin() + B_quad.num_elements(), transformed_B_quad.origin(), F);
 
    return transformed_B_quad;
  };
 
 
  template <typename ScalarBasisType, size_t N>
  Family<TensorizedVectorFamily<ScalarBasisType, N>> GenericTensorization(
          const ScalarBasisType & B,  
          const std::vector<Eigen::VectorXd> & v 
  )
  {
      size_t r = B.dimension();
      size_t k = v.size();
 
      Eigen::MatrixXd M = Eigen::MatrixXd::Zero(r*k, r*N);
 
      for (size_t i=0; i < k; i++){
          // Check the dimension of vector v[i]
          assert(v[i].size() == N);
 
          // Fill in M
          for (size_t j=0; j < N; j++){
              M.block(i*r, j*r, r, r) = v[i](j) * Eigen::MatrixXd::Identity(r,r);
          }
      }
 
      TensorizedVectorFamily<ScalarBasisType, N> tensbasis(B);
      return Family<TensorizedVectorFamily<ScalarBasisType, N>>(tensbasis, M);
  }
 
 
  Eigen::MatrixXd PermuteTensorization( const size_t a,  
                                        const size_t b   
                                      );
 
  inline static std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> 
  symmetrise_matrix = [](const Eigen::MatrixXd & x)->Eigen::MatrixXd { return 0.5*(x+x.transpose());};
 
  inline static std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> 
  skew_symmetrise_matrix = [](const Eigen::MatrixXd & x)->Eigen::MatrixXd { return 0.5*(x-x.transpose());};
 
 
  template <typename ScalarBasisType, size_t N>
  Family<MatrixFamily<ScalarBasisType, N>> IsotropicMatrixFamily(
                                                          const ScalarBasisType & B  
                                                          )
  {
    size_t r = B.dimension();
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(r, r*N*N);
 
    for (size_t k=0; k < N; k++){
        M.block(0, k*r*(N+1), r, r) = Eigen::MatrixXd::Identity(r, r);
    }
 
    MatrixFamily<ScalarBasisType, N> matbasis(B);
    return Family<MatrixFamily<ScalarBasisType, N>>(matbasis, M);
  }
  //------------------------------------------------------------------------------
  //        BASIS EVALUATION
  //------------------------------------------------------------------------------
  
  namespace detail {
 
    template<typename BasisType, BasisFunctionE BasisFunction>
    struct basis_evaluation_traits {};
 
    //---------- GENERIC EVALUATION TRAITS -----------------------------------------
    // Evaluate the function value at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, Function>
    {
      static_assert(BasisType::hasFunction, "Call to function not available");
      typedef typename BasisType::FunctionValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
      {
        return basis.function(i, x);
      }
           
      // Computes function value at quadrature node iqn, knowing values of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const BasisType &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<ReturnValue, 2> &ancestor_value_quad
                                         )
      {
        return basis.function(i, iqn, ancestor_value_quad);
      }
 
    };
 
    // Evaluate the gradient value at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, Gradient>
    {
      static_assert(BasisType::hasGradient, "Call to gradient not available");
      typedef typename BasisType::GradientValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
      {
        return basis.gradient(i, x);
      }
      
      // Computes gradient value at quadrature node iqn, knowing gradients of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const BasisType &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<ReturnValue, 2> &ancestor_gradient_quad
                                         )
      {
        return basis.gradient(i, iqn, ancestor_gradient_quad);
      }
    };
 
    // AC - it's supposed to be used only in case of vector variables (vhhospace)
    // Evaluate the symmetric part of gradient at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, SymmetricGradient>
    {
        static_assert(BasisType::hasGradient, "Call to gradient not available");
        typedef typename BasisType::GradientValue ReturnValue;
        static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
        {
            return 0.5 * (basis.gradient(i, x) + basis.gradient(i, x).transpose());
        }
 
        // Computes gradient value at quadrature node iqn, knowing gradients of ancestor basis at quadrature nodes
        static inline ReturnValue evaluate(
                const BasisType &basis,
                size_t i,
                size_t iqn,
                const boost::multi_array<ReturnValue, 2> &ancestor_gradient_quad
        )
        {
            return 0.5 * (basis.gradient(i, iqn, ancestor_gradient_quad) + basis.gradient(i, iqn, ancestor_gradient_quad).transpose());
        }
    };
 
      // AC - it's supposed to be used only in case of vector variables (vhhospace)
      // Evaluate the symmetric part of gradient at x
      template<typename BasisType>
      struct basis_evaluation_traits<BasisType, SkewsymmetricGradient>
      {
          static_assert(BasisType::hasGradient, "Call to gradient not available");
          typedef typename BasisType::GradientValue ReturnValue;
          static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
          {
              return 0.5 * (basis.gradient(i, x) - basis.gradient(i, x).transpose());
          }
 
          // Computes gradient value at quadrature node iqn, knowing gradients of ancestor basis at quadrature nodes
          static inline ReturnValue evaluate(
                  const BasisType &basis,
                  size_t i,
                  size_t iqn,
                  const boost::multi_array<ReturnValue, 2> &ancestor_gradient_quad
          )
          {
              return 0.5 * (basis.gradient(i, iqn, ancestor_gradient_quad) - basis.gradient(i, iqn, ancestor_gradient_quad).transpose());
          }
      };
 
    // Evaluate the curl value at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, Curl>
    {
      static_assert(BasisType::hasCurl, "Call to curl not available");
      typedef typename BasisType::CurlValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
      {
        return basis.curl(i, x);
      }
      
      // Computes curl value at quadrature node iqn, knowing curls of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const BasisType &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<ReturnValue, 2> &ancestor_curl_quad
                                         )
      {
        return basis.curl(i, iqn, ancestor_curl_quad);
      }
    };
 
    // Evaluate the divergence value at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, Divergence>
    {
      static_assert(BasisType::hasDivergence, "Call to divergence not available");
      typedef typename BasisType::DivergenceValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
      {
        return basis.divergence(i, x);
      }
      
      // Computes divergence value at quadrature node iqn, knowing divergences of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const BasisType &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<ReturnValue, 2> &ancestor_divergence_quad
                                         )
      {
        return basis.divergence(i, iqn, ancestor_divergence_quad);
      }
    };
 
    // Evaluate the rotor value at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, Rotor>
    {
      static_assert(BasisType::hasRotor, "Call to rotor not available");
      typedef typename BasisType::RotorValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
      {
        return basis.rotor(i, x);
      }
      
      // Computes rotor value at quadrature node iqn, knowing rotors of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const BasisType &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<ReturnValue, 2> &ancestor_rotor_quad
                                         )
      {
        return basis.rotor(i, iqn, ancestor_rotor_quad);
      }
    };
    
    // Evaluate the Hessian value at x
    template<typename BasisType>
    struct basis_evaluation_traits<BasisType, Hessian>
    {
      static_assert(BasisType::hasHessian, "Call to Hessian not available");
      typedef typename BasisType::HessianValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType & basis, size_t i, const VectorRd & x)
      {
        return basis.hessian(i, x);
      }
      
      // Computes Hessian value at quadrature node iqn, knowing Hessians of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const BasisType &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<ReturnValue, 2> &ancestor_hessian_quad
                                         )
      {
        return basis.hessian(i, iqn, ancestor_hessian_quad);
      }
    };
    
        // Evaluate the curl curl value at x
    template <typename BasisType>
    struct basis_evaluation_traits<BasisType, CurlCurl>
    {
      static_assert(BasisType::hasCurlCurl, "Call to curl curl not available");
      typedef typename BasisType::CurlCurlValue ReturnValue;
      static inline ReturnValue evaluate(const BasisType &basis, size_t i, const VectorRd &x)
      {
        return basis.curlcurl(i, x);
      }
            
      // Computes hessian value at quadrature node iqn, knowing hessian of ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                    const BasisType &basis, 
                    size_t i, 
                    size_t iqn,
                    const boost::multi_array<ReturnValue, 2> &ancestor_curlcurl_quad
                    )
      {
        return basis.curlcurl(i, iqn, ancestor_curlcurl_quad);
      }
 
    };
 
    //---------- TENSORIZED FAMILY EVALUATION TRAITS -----------------------------------------
    // Evaluate the value at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, Function>
    {
      static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasFunction, "Call to function not available");
      
      typedef typename TensorizedVectorFamily<ScalarBasisType, N>::FunctionValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Function; // Type of values needed from ancestor basis
      typedef double AncestorBasisFunctionValue; 
      
      // Computes function values at x
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.function(i, x);
      }
      
      // Computes function value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<double, 2> &ancestor_basis_quad
                                         )
      {
        return basis.function(i, iqn, ancestor_basis_quad);
      }
    };
 
    // Evaluate the gradient at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, Gradient>
    {
      static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
      
      typedef typename TensorizedVectorFamily<ScalarBasisType, N>::GradientValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
      typedef VectorRd AncestorBasisFunctionValue;
          
      // Computes gradient value at x
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.gradient(i, x);
      }
 
      // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
                                         )
      {
        return basis.gradient(i, iqn, ancestor_basis_quad);
      }
    };
 
    // AC - supposed to be used only in case of vector variables (vhhospace)
    // Evaluate the symmetric gradient at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, SymmetricGradient>
    {
        static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
 
        typedef typename TensorizedVectorFamily<ScalarBasisType, N>::GradientValue ReturnValue;
        static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
        typedef VectorRd AncestorBasisFunctionValue;
 
        // Computes gradient value at x
        static inline ReturnValue evaluate(
                const TensorizedVectorFamily<ScalarBasisType, N> &basis,
                size_t i,
                const VectorRd &x
        )
        {
            return 0.5 * (basis.gradient(i, x) + basis.gradient(i, x).transpose()) ;
        }
 
        // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
        static inline ReturnValue evaluate(
                const TensorizedVectorFamily<ScalarBasisType, N> &basis,
                size_t i,
                size_t iqn,
                const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
        )
        {
            return 0.5 * (basis.gradient(i, iqn, ancestor_basis_quad) + basis.gradient(i, iqn, ancestor_basis_quad).transpose());
        }
    };
 
      // AC - supposed to be used only in case of vector variables (so gradient is a matrix)
      // Evaluate the symmetric gradient at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
      template <typename ScalarBasisType, size_t N>
      struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, SkewsymmetricGradient>
      {
          static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
 
          typedef typename TensorizedVectorFamily<ScalarBasisType, N>::GradientValue ReturnValue;
          static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
          typedef VectorRd AncestorBasisFunctionValue;
 
          // Computes gradient value at x
          static inline ReturnValue evaluate(
                  const TensorizedVectorFamily<ScalarBasisType, N> &basis,
                  size_t i,
                  const VectorRd &x
          )
          {
              return 0.5 * (basis.gradient(i, x) - basis.gradient(i, x).transpose()) ;
          }
 
          // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
          static inline ReturnValue evaluate(
                  const TensorizedVectorFamily<ScalarBasisType, N> &basis,
                  size_t i,
                  size_t iqn,
                  const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
          )
          {
              return 0.5 * (basis.gradient(i, iqn, ancestor_basis_quad) - basis.gradient(i, iqn, ancestor_basis_quad).transpose());
          }
      };
 
    // Evaluate the curl at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, Curl>
    {
      static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasCurl, "Call to curl not available");
      
      typedef typename TensorizedVectorFamily<ScalarBasisType, N>::CurlValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
      typedef VectorRd AncestorBasisFunctionValue;
 
      // Computes curl values at x
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.curl(i, x);
      }
 
      // Computes curl values at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
                                         )
      {
        return basis.curl(i, iqn, ancestor_basis_quad);
      }
      
    };
 
    // Evaluate the divergence at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, Divergence>
    {
      static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasDivergence, "Call to divergence not available");
      
      typedef typename TensorizedVectorFamily<ScalarBasisType, N>::DivergenceValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
      typedef VectorRd AncestorBasisFunctionValue;
 
      // Computes divergence value at x
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.divergence(i, x);
      }
 
      // Computes divergence value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
                                         )
      {
        return basis.divergence(i, iqn, ancestor_basis_quad);
      }
      
    };
 
    // Evaluate the rotor at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, Rotor>
    {
      static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasRotor, "Call to rotor not available");
      
      typedef typename TensorizedVectorFamily<ScalarBasisType, N>::RotorValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
      typedef VectorRd AncestorBasisFunctionValue;
 
      // Computes divergence value at x
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.rotor(i, x);
      }
 
      // Computes rotor value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
                                         )
      {
        return basis.rotor(i, iqn, ancestor_basis_quad);
      }
      
    };
 
        // Evaluate the curl curl at x of a TensorizedVectorFamily (includes information on the ancestor basis, to optimise eval_quad for tensorized bases)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<TensorizedVectorFamily<ScalarBasisType, N>, CurlCurl>
    {
      static_assert(TensorizedVectorFamily<ScalarBasisType, N>::hasCurlCurl, "Call to curl curl not available");
      
      typedef typename TensorizedVectorFamily<ScalarBasisType, N>::CurlCurlValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Hessian; // Type of values needed from ancestor basis
      typedef MatrixRd AncestorBasisFunctionValue;
 
      // Computes curl curl value at x
      static inline ReturnValue evaluate(
                    const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                    size_t i, 
                    const VectorRd &x
                    )
      {
        return basis.curlcurl(i, x);
      }
      
     
 
      // Computes curl curl value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                    const TensorizedVectorFamily<ScalarBasisType, N> &basis, 
                    size_t i, 
                    size_t iqn,
                    const boost::multi_array<MatrixRd, 2> &ancestor_basis_quad
                    )
      {
        return basis.curlcurl(i, iqn, ancestor_basis_quad);
      }
      
    };
 
 
    //---------- MATRIX FAMILY EVALUATION TRAITS -----------------------------------------
    // Evaluate the value at x of a MatrixFamily (includes information on the ancestor basis, to optimise eval_quad for matrix families)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<MatrixFamily<ScalarBasisType, N>, Function>
    {
      static_assert(MatrixFamily<ScalarBasisType, N>::hasFunction, "Call to function not available");
      
      typedef typename MatrixFamily<ScalarBasisType, N>::FunctionValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Function; // Type of values needed from ancestor basis
      typedef double AncestorBasisFunctionValue; 
      
      // Computes function values at x
      static inline ReturnValue evaluate(
                                         const MatrixFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.function(i, x);
      }
      
      // Computes function value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const MatrixFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<double, 2> &ancestor_basis_quad
                                         )
      {
        return basis.function(i, iqn, ancestor_basis_quad);
      }
    };
 
    // Evaluate the divergence at x of a MatrixFamily (includes information on the ancestor basis, to optimise eval_quad for matrix families)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<MatrixFamily<ScalarBasisType, N>, Divergence>
    {
      static_assert(MatrixFamily<ScalarBasisType, N>::hasDivergence, "Call to divergence not available");
      
      typedef typename MatrixFamily<ScalarBasisType, N>::DivergenceValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
      typedef VectorRd AncestorBasisFunctionValue;
 
      // Computes divergence value at x
      static inline ReturnValue evaluate(
                                         const MatrixFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.divergence(i, x);
      }
 
      // Computes divergence value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const MatrixFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
                                         )
      {
        return basis.divergence(i, iqn, ancestor_basis_quad);
      }
      
    };
 
 
    // Evaluate the gradient at x of a MatrixFamily (includes information on the ancestor basis, to optimise eval_quad for matrix families)
    template <typename ScalarBasisType, size_t N>
    struct basis_evaluation_traits<MatrixFamily<ScalarBasisType, N>, Gradient>
    {
      static_assert(MatrixFamily<ScalarBasisType, N>::hasGradient, "Call to gradient not available");
      
      typedef typename MatrixFamily<ScalarBasisType, N>::GradientValue ReturnValue;
      static const BasisFunctionE AncestorBasisFunction = Gradient; // Type of values needed from ancestor basis
      typedef VectorRd AncestorBasisFunctionValue;
 
      // Computes gradient value at x
      static inline ReturnValue evaluate(
                                         const MatrixFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         const VectorRd &x
                                         )
      {
        return basis.gradient(i, x);
      }
 
      // Computes gradient value at quadrature node iqn, knowing ancestor basis at quadrature nodes
      static inline ReturnValue evaluate(
                                         const MatrixFamily<ScalarBasisType, N> &basis, 
                                         size_t i, 
                                         size_t iqn,
                                         const boost::multi_array<VectorRd, 2> &ancestor_basis_quad
                                         )
      {
        return basis.gradient(i, iqn, ancestor_basis_quad);
      }
      
    };
 
  } // end of namespace detail
 
 
  //-----------------------------------EVALUATE_QUAD--------------------------------------
  
  template<BasisFunctionE BasisFunction>
  struct evaluate_quad {
    template<typename BasisType>
    static boost::multi_array<typename detail::basis_evaluation_traits<BasisType, BasisFunction>::ReturnValue, 2>
    compute(
            const BasisType & basis,    
            const QuadratureRule & quad 
            )
    {
      typedef detail::basis_evaluation_traits<BasisType, BasisFunction> traits;
 
      boost::multi_array<typename traits::ReturnValue, 2>
        basis_quad( boost::extents[basis.dimension()][quad.size()] );
      
      for (size_t i = 0; i < basis.dimension(); i++) {
        for (size_t iqn = 0; iqn < quad.size(); iqn++) {
          basis_quad[i][iqn] = traits::evaluate(basis, i, quad[iqn].vector());
        } // for iqn
      } // for i
   
      return basis_quad;
    }
 
    template<typename BasisType>
    static boost::multi_array<typename detail::basis_evaluation_traits<Family<BasisType>, BasisFunction>::ReturnValue, 2>
    compute(
            const Family<BasisType> & basis, 
            const QuadratureRule & quad      
            )
    {
      typedef detail::basis_evaluation_traits<Family<BasisType>, BasisFunction> traits;
 
      boost::multi_array<typename traits::ReturnValue, 2>
        basis_quad( boost::extents[basis.dimension()][quad.size()] );
 
      // Compute values at quadrature note on ancestor basis, and then do a transformation
      const auto & ancestor_basis = basis.ancestor();
      boost::multi_array<typename traits::ReturnValue, 2>
        ancestor_basis_quad = evaluate_quad<BasisFunction>::compute(ancestor_basis, quad);
  
      Eigen::MatrixXd M = basis.matrix();
      for (size_t i = 0; i < size_t(M.rows()); i++) {
        for (size_t iqn = 0; iqn < quad.size(); iqn++) {
          basis_quad[i][iqn] = M(i,0) * ancestor_basis_quad[0][iqn];
        }
        for (size_t j = 1; j < size_t(M.cols()); j++) {
          for (size_t iqn = 0; iqn < quad.size(); iqn++) {
            basis_quad[i][iqn] += M(i,j) * ancestor_basis_quad[j][iqn];
          }
        } // for j
      } // for i
   
      return basis_quad;
    }
    
    template <typename BasisType, size_t N>
    static boost::multi_array<typename detail::basis_evaluation_traits<TensorizedVectorFamily<BasisType, N>, BasisFunction>::ReturnValue, 2>
    compute(
            const TensorizedVectorFamily<BasisType, N> &basis, 
            const QuadratureRule &quad      
            )
    {
      typedef detail::basis_evaluation_traits<TensorizedVectorFamily<BasisType, N>, BasisFunction> traits;
 
      boost::multi_array<typename traits::ReturnValue, 2>
        basis_quad(boost::extents[basis.dimension()][quad.size()]);
 
      const auto &ancestor_basis = basis.ancestor();
      const boost::multi_array<typename traits::AncestorBasisFunctionValue, 2> ancestor_basis_quad 
        = evaluate_quad<traits::AncestorBasisFunction>::compute(ancestor_basis, quad);
 
      for (size_t i = 0; i < basis.dimension(); i++){
        for (size_t iqn = 0; iqn < quad.size(); iqn++){
          basis_quad[i][iqn] = traits::evaluate(basis, i, iqn, ancestor_basis_quad);
        }
      }
      return basis_quad;
    }
  
    template <typename BasisType, size_t N>
    static boost::multi_array<typename detail::basis_evaluation_traits<MatrixFamily<BasisType, N>, BasisFunction>::ReturnValue, 2>
    compute(
            const MatrixFamily<BasisType, N> &basis, 
            const QuadratureRule &quad      
            )
    {
      typedef detail::basis_evaluation_traits<MatrixFamily<BasisType, N>, BasisFunction> traits;
 
      boost::multi_array<typename traits::ReturnValue, 2>
        basis_quad(boost::extents[basis.dimension()][quad.size()]);
 
      const auto &ancestor_basis = basis.ancestor();
      const boost::multi_array<typename traits::AncestorBasisFunctionValue, 2> ancestor_basis_quad 
        = evaluate_quad<traits::AncestorBasisFunction>::compute(ancestor_basis, quad);
 
      for (size_t i = 0; i < basis.dimension(); i++){
        for (size_t iqn = 0; iqn < quad.size(); iqn++){
          basis_quad[i][iqn] = traits::evaluate(basis, i, iqn, ancestor_basis_quad);
        }
      }
      return basis_quad;
    }
  };
 
  //------------------------------------------------------------------------------
  //          ORTHONORMALISATION
  //------------------------------------------------------------------------------
 
  template<typename T>
  Eigen::MatrixXd gram_schmidt(
                               boost::multi_array<T, 2> & basis_eval, 
                               const std::function<double(size_t, size_t)> & inner_product 
                               )
  {
    auto induced_norm = [inner_product](size_t i)->double
    {
      return std::sqrt( inner_product(i, i) );
    };
 
    // Number of basis functions
    size_t Nb = basis_eval.shape()[0];
    // Number of quadrature nodes
    size_t Nqn = basis_eval.shape()[1];
 
    Eigen::MatrixXd B = Eigen::MatrixXd::Zero(Nb, Nb);
    
    // Normalise the first element
    double norm = induced_norm(0);
    for (size_t iqn = 0; iqn < Nqn; iqn++) {
      basis_eval[0][iqn] /= norm;
    } // for iqn
    B(0,0) = 1./norm;
    
    for (size_t ib = 1; ib < Nb; ib++) {
      // 'coeffs' represents the coefficients of the ib-th orthogonal function on the ON basis functions 0 to ib-1
      Eigen::RowVectorXd coeffs = Eigen::RowVectorXd::Zero(ib);
      for (size_t pb = 0; pb < ib; pb++) {
        coeffs(pb) = - inner_product(ib, pb);
      } // for pb
      
      // store the values of the orthogonalised ib-th basis function
      for (size_t pb = 0; pb < ib; pb++) {
        for (size_t iqn = 0; iqn < Nqn; iqn++) {
          basis_eval[ib][iqn] += coeffs(pb) * basis_eval[pb][iqn];
        } // for iqn
      } // for pb
      
      // normalise ib-th basis function
      double norm = induced_norm(ib);
      for (size_t iqn = 0; iqn < Nqn; iqn++) {
        basis_eval[ib][iqn] /= norm;
      } // for iqn
      coeffs /= norm;
      // Compute ib-th row of B.
      // B.topLeftCorner(ib, ib) contains the rows that represent the ON basis functions 0 to ib-1 on
      // the original basis. Multiplying on the left by coeffs gives the coefficients of the ON basis
      // functions on the original basis functions from 0 to ib-1.
      B.block(ib, 0, 1, ib) = coeffs * B.topLeftCorner(ib, ib);
      B(ib, ib) = 1./norm;
    }
    return B;
  }
 
  //------------------------------------------------------------------------------
 
  double scalar_product(const double & x, const double & y);
 
  double scalar_product(const VectorRd & x, const VectorRd & y);
 
  template<int N>
  double scalar_product(const Eigen::Matrix<double, N, N> & x, const Eigen::Matrix<double, N, N> & y)
  {
    return (x.transpose() * y).trace();
  }
 
  template <typename Value>
  boost::multi_array<double, 2> scalar_product(
                                               const boost::multi_array<Value, 2> &basis_quad, 
                                               const Value &v                                  
                                               )
  {
    boost::multi_array<double, 2> basis_dot_v_quad(boost::extents[basis_quad.shape()[0]][basis_quad.shape()[1]]);
    std::transform(basis_quad.origin(), basis_quad.origin() + basis_quad.num_elements(),
                   basis_dot_v_quad.origin(), [&v](const Value &x) -> double { return scalar_product(x,v); });
    return basis_dot_v_quad;
  }
 
  template <typename Value>
  boost::multi_array<double, 2> scalar_product(
                                               const boost::multi_array<Value, 2> & basis_quad, 
                                               const std::vector<Value> & v 
                                               )
  {
    boost::multi_array<double, 2> basis_dot_v_quad(boost::extents[basis_quad.shape()[0]][basis_quad.shape()[1]]);
    for (std::size_t i = 0; i < basis_quad.shape()[0]; i++) {
      for (std::size_t iqn = 0; iqn < basis_quad.shape()[1]; iqn++) {
        basis_dot_v_quad[i][iqn] = scalar_product(basis_quad[i][iqn], v[iqn]);
      } // for iqn
    } // for i
    return basis_dot_v_quad;
  }
 
  boost::multi_array<VectorRd, 2> matrix_vector_product(
                                                        const boost::multi_array<MatrixRd, 2> & basis_quad, 
                                                        const std::vector<VectorRd> & v_quad 
                                                        );
 
  boost::multi_array<VectorRd, 2> matrix_vector_product(
                                                        const std::vector<MatrixRd> & m_quad, 
                                                        const boost::multi_array<VectorRd, 2> & basis_quad 
                                                        );
 
  boost::multi_array<VectorRd, 2> vector_matrix_product(
                                                        const std::vector<VectorRd> & v_quad, 
                                                        const boost::multi_array<MatrixRd, 2> & basis_quad 
                                                        );
  
  template<typename BasisType>
  Family<BasisType> l2_orthonormalize(
                                      const BasisType & basis, 
                                      const QuadratureRule & qr,  
                                      boost::multi_array<typename BasisType::FunctionValue, 2> & basis_quad 
                                      )
  {
    // Check that the basis evaluation and quadrature rule are coherent
    assert ( basis.dimension() == basis_quad.shape()[0] && qr.size() == basis_quad.shape()[1] );
 
    // The inner product between the i-th and j-th basis vectors
    std::function<double(size_t, size_t)> inner_product
      = [&basis_quad, &qr](size_t i, size_t j)->double
      {
        double r = 0.;
        for(size_t iqn = 0; iqn < qr.size(); iqn++) {
          r += qr[iqn].w * scalar_product(basis_quad[i][iqn], basis_quad[j][iqn]);
        } // for iqn
        return r;
      };
 
    Eigen::MatrixXd B = gram_schmidt(basis_quad, inner_product);
 
    return Family<BasisType>(basis, B);    
  }
 
  /* This method is more robust that the previous one based on values at quadrature nodes */
  template <typename BasisType>
  Family<BasisType> l2_orthonormalize(
                                      const BasisType &basis,       
                                      const Eigen::MatrixXd & GM    
                                      )
  {
    // Check that the basis and Gram matrix are coherent
    assert(basis.dimension() == size_t(GM.rows()) && GM.rows() == GM.cols());
 
    Eigen::MatrixXd L = GM.llt().matrixL();
 
    return Family<BasisType>(basis, L.inverse());
  }
 
 
  //------------------------------------------------------------------------------
  //      GRAM MATRICES
  //------------------------------------------------------------------------------
 
  template<typename FunctionValue>
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<FunctionValue, 2> & B1, 
                                      const boost::multi_array<FunctionValue, 2> & B2, 
                                      const QuadratureRule & qr,        
                                      const size_t nrows,   
                                      const size_t ncols,   
                                      const std::string sym = "nonsym"     
                                      )
  {
    // Check that the basis evaluation and quadrature rule are coherent
    assert ( qr.size() == B1.shape()[1] && qr.size() == B2.shape()[1] );
    // Check that we don't ask for more members of family than available
    assert ( nrows <= B1.shape()[0] && ncols <= B2.shape()[0] );
 
    // Re-cast quadrature weights into ArrayXd to make computations faster
    Eigen::ArrayXd qr_weights = Eigen::ArrayXd::Zero(qr.size());
    for (size_t iqn = 0; iqn < qr.size(); iqn++){
      qr_weights(iqn) = qr[iqn].w;
    }
 
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(nrows, ncols);
    for (size_t i = 0; i < nrows; i++) {
      size_t jcut = 0;
      if ( sym=="sym" ) jcut = i;
      for (size_t j = 0; j < jcut; j++){
        M(i, j) = M(j, i);
      }
      for(size_t j = jcut; j < ncols; j++) {
        std::vector<double> tmp(B1.shape()[1]);
        // Extract values at quadrature nodes for elements i of B1 and j of B2
        auto B1i = B1[ boost::indices[i][boost::multi_array_types::index_range(0, B1.shape()[1])] ];
        auto B2j = B2[ boost::indices[j][boost::multi_array_types::index_range(0, B1.shape()[1])] ];
        // Compute scalar product of B1i and B2j, and recast it as ArrayXd
        std::transform(B1i.begin(), B1i.end(), B2j.begin(), tmp.begin(), [](FunctionValue a, FunctionValue b)->double { return scalar_product(a,b);});
        Eigen::ArrayXd tmp_array = Eigen::Map<Eigen::ArrayXd, Eigen::Unaligned>(tmp.data(), tmp.size());
        // Multiply by quadrature weights and sum (using .sum() of ArrayXd makes this step faster than a loop)
        M(i,j) = (qr_weights * tmp_array).sum();
 
 
        // Simple version with loop (replaces everything above after the for on j)
        /*
          for (size_t iqn = 0; iqn < qr.size(); iqn++) {
          M(i,j) += qr[iqn].w * scalar_product(B1[i][iqn], B2[j][iqn]);
          } // for iqn
        */
 
      } // for j
    } // for i
    return M;
  }
 
  template<typename FunctionValue>
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<FunctionValue, 2> & B1, 
                                      const boost::multi_array<FunctionValue, 2> & B2, 
                                      const QuadratureRule & qr,        
                                      const std::string sym = "nonsym"  
                                      )
  {
    return compute_gram_matrix<FunctionValue>(B1, B2, qr, B1.shape()[0], B2.shape()[0], sym);
  }
 
  template<typename FunctionValue>
  inline Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<FunctionValue, 2> & B, 
                                             const QuadratureRule & qr        
                                             )
  {
    return compute_gram_matrix<FunctionValue>(B, B, qr, "sym");
  }
 
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<VectorRd, 2> & B1,  
                                      const boost::multi_array<double, 2> & B2,         
                                      const QuadratureRule & qr                  
                                      );
 
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<double, 2> & B1, 
                                      const boost::multi_array<double, 2> & B2, 
                                      const QuadratureRule & qr,        
                                      const size_t nrows,   
                                      const size_t ncols,   
                                      const std::string sym = "nonsym"     
                                      );
 
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<double, 2> & B1, 
                                      const boost::multi_array<double, 2> & B2, 
                                      const QuadratureRule & qr,        
                                      const std::string sym = "nonsym"     
                                      );
 
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<VectorRd, 2> & B1, 
                                      const boost::multi_array<VectorRd, 2> & B2, 
                                      const QuadratureRule & qr,        
                                      const size_t nrows,   
                                      const size_t ncols,   
                                      const std::string sym = "nonsym"     
                                      );
 
  Eigen::MatrixXd compute_gram_matrix(const boost::multi_array<VectorRd, 2> & B1, 
                                      const boost::multi_array<VectorRd, 2> & B2, 
                                      const QuadratureRule & qr,        
                                      const std::string sym = "nonsym"     
                                      );
                                      
  template<typename ScalarFamilyType, size_t N>
  Eigen::MatrixXd compute_gram_matrix(const MatrixFamily<ScalarFamilyType, N> & MatFam, 
                                      const boost::multi_array<double, 2> & scalar_family_quad, 
                                      const QuadratureRule & qr        
                                      )
  {
    Eigen::MatrixXd Gram = Eigen::MatrixXd::Zero(MatFam.dimension(), MatFam.dimension());
        
    // Gram matrix of the scalar family
    Eigen::MatrixXd scalarGram = compute_gram_matrix<double>(scalar_family_quad, qr);
    for (size_t i = 0; i < MatFam.dimension(); i = i+scalarGram.rows()){
      Gram.block(i, i, scalarGram.rows(), scalarGram.cols()) = scalarGram;
    }   
    
    return Gram;
  }
                                      
 
 
  template <typename T>
  Eigen::VectorXd integrate(
                            const FType<T> &f,        
                            const BasisQuad<T> &B,    
                            const QuadratureRule &qr, 
                            size_t n_rows = 0         
                            )
  {
    // If default, set n_rows to size of family
    if (n_rows == 0)
      {
        n_rows = B.shape()[0];
      }
 
    // Number of quadrature nodes
    const size_t num_quads = qr.size();
 
    // Check number of quadrature nodes is compatible with B
    assert(num_quads == B.shape()[1]);
 
    // Check that we don't ask for more members of family than available
    assert(n_rows <= B.shape()[0]);
 
    Eigen::VectorXd V = Eigen::VectorXd::Zero(n_rows);
 
    for (size_t iqn = 0; iqn < num_quads; iqn++)
      {
        double qr_weight = qr[iqn].w;
        T f_on_qr = f(qr[iqn].vector());
        for (size_t i = 0; i < n_rows; i++)
          {
            V(i) += qr_weight * scalar_product(B[i][iqn], f_on_qr);
          }
      }
 
    return V;
  }
 
  template <typename T, typename U>
  Eigen::MatrixXd compute_weighted_gram_matrix(
                                               const FType<U> &f,               
                                               const BasisQuad<T> &B1,          
                                               const BasisQuad<T> &B2,          
                                               const QuadratureRule &qr,        
                                               size_t n_rows = 0,               
                                               size_t n_cols = 0,               
                                               const std::string sym = "nonsym" 
                                               )
  {
    // If default, set n_rows and n_cols to size of families
    if (n_rows == 0 && n_cols == 0)
      {
        n_rows = B1.shape()[0];
        n_cols = B2.shape()[0];
      }
 
    // Number of quadrature nodes
    const size_t num_quads = qr.size();
    // Check number of quadrature nodes is compatible with B1 and B2
    assert(num_quads == B1.shape()[1] && num_quads == B2.shape()[1]);
    // Check that we don't ask for more members of family than available
    assert(n_rows <= B1.shape()[0] && n_cols <= B2.shape()[0]);
 
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(n_rows, n_cols);
    for (size_t iqn = 0; iqn < num_quads; iqn++)
      {
        double qr_weight = qr[iqn].w;
        U f_on_qr = f(qr[iqn].vector());
        for (size_t i = 0; i < n_rows; i++)
          {
            T f_B1 = f_on_qr * B1[i][iqn];
            size_t jcut = 0;
            if (sym == "sym")
              jcut = i;
            for (size_t j = 0; j < jcut; j++)
              {
                M(i, j) = M(j, i);
              }
            for (size_t j = jcut; j < n_cols; j++)
              {
                M(i, j) += qr_weight * scalar_product(f_B1, B2[j][iqn]);
              }
          }
      }
    return M;
  }
 
  template <typename T, typename U>
  Eigen::MatrixXd compute_weighted_gram_matrix(
                                               const FType<U> &f,        
                                               const BasisQuad<T> &B1,   
                                               const BasisQuad<T> &B2,   
                                               const QuadratureRule &qr, 
                                               const std::string sym     
                                               )
  {
    return compute_weighted_gram_matrix(f, B1, B2, qr, B1.shape()[0], B2.shape()[0], sym);
  }
 
  Eigen::MatrixXd compute_weighted_gram_matrix(
                                               const FType<VectorRd> &f,      
                                               const BasisQuad<VectorRd> &B1, 
                                               const BasisQuad<double> &B2,   
                                               const QuadratureRule &qr,      
                                               size_t n_rows = 0,             
                                               size_t n_cols = 0              
                                               );
  
 
  Eigen::MatrixXd compute_weighted_gram_matrix(
                                               const FType<VectorRd> &f,      
                                               const BasisQuad<double> &B1,   
                                               const BasisQuad<VectorRd> &B2, 
                                               const QuadratureRule &qr,      
                                               size_t n_rows = 0,             
                                               size_t n_cols = 0              
                                               );
  
  template <typename Value>
  Eigen::MatrixXd compute_closure_matrix(
                                         const boost::multi_array<Value, 2> & f_quad, 
                                         const boost::multi_array<Value, 2> & g_quad, 
                                         const QuadratureRule & qr_f,                 
                                         const QuadratureRule & qr_g                  
                                         )
  {
    const size_t dimf = f_quad.shape()[0];
    const size_t dimg = g_quad.shape()[0];
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(dimf, dimg);
  
    std::vector<Value> int_f(dimf);
    for (size_t i = 0; i < dimf; i++){
      int_f[i] = qr_f[0].w * f_quad[i][0];
      for (size_t iqn = 1; iqn < qr_f.size(); iqn++){
        int_f[i] += qr_f[iqn].w * f_quad[i][iqn];
      }      
    }
    std::vector<Value> int_g(dimg);
    for (size_t j = 0; j < dimg; j++){
      int_g[j] = qr_g[0].w * g_quad[j][0];
      for (size_t iqn = 1; iqn < qr_g.size(); iqn++){
        int_g[j] += qr_g[iqn].w * g_quad[j][iqn];
      }      
    }
    
    // Compute (int w_i) and create dot product with (int w_j) for j<=i, then fill in M by symmetry
    for (size_t i = 0; i < dimf; i++){
      for (size_t j = 0; j < dimg; j++){
        M(i,j) = scalar_product(int_f[i], int_g[j]);
      }
    }
  
    return M;  
  }
 
  template <typename Value>
  Eigen::MatrixXd compute_closure_matrix(
                                         const boost::multi_array<Value, 2> & f_quad, 
                                         const boost::multi_array<Value, 2> & g_quad, 
                                         const QuadratureRule & qr                    
                                         )
  {
    return compute_closure_matrix(f_quad, g_quad, qr, qr);
  }
 
  template <typename Value>
  Eigen::MatrixXd compute_closure_matrix(
                                         const boost::multi_array<Value, 2> & f_quad, 
                                         const QuadratureRule & qr                    
                                         )
  {
    return compute_closure_matrix(f_quad, f_quad, qr, qr);
  }
 
  //------------------------------------------------------------------------------
  //        L2 projection of a function
  //------------------------------------------------------------------------------
 
  template<typename BasisType>
  Eigen::VectorXd l2_projection(
                                const std::function<typename BasisType::FunctionValue(const VectorRd &)> & f, 
                                const BasisType & basis, 
                                QuadratureRule & quad, 
                                const boost::multi_array<typename BasisType::FunctionValue, 2> & basis_quad, 
                                const Eigen::MatrixXd &mass_basis = Eigen::MatrixXd::Zero(1,1)               
                              )
  {
    Eigen::MatrixXd Mass = mass_basis;
    if (Mass.norm() < 1e-13){
      Mass = compute_gram_matrix(basis_quad, quad);
    }
 
    Eigen::LDLT<Eigen::MatrixXd> cholesky_mass(Mass);
    Eigen::VectorXd b = Eigen::VectorXd::Zero(basis.dimension());
    for (size_t i = 0; i < basis.dimension(); i++) {
      for (size_t iqn = 0; iqn < quad.size(); iqn++) {
        b(i) += quad[iqn].w * scalar_product(f(quad[iqn].vector()), basis_quad[i][iqn]);
      } // for iqn
    } // for i
 
    return cholesky_mass.solve(b);
  }
 
  //------------------------------------------------------------------------------
  //        Decomposition on a polynomial basis
  //------------------------------------------------------------------------------
 
 
  template <typename BasisType>
  struct DecomposePoly
  {
 
    DecomposePoly(
            const Edge &E,    
            const BasisType &basis  
            ):
          m_dim(basis.dimension()),
          m_basis(basis),
          m_on_basis(basis, Eigen::MatrixXd::Identity(m_dim,m_dim))  // Need to initialise before we compute the real ON basis
      {
        /* The decomposition will be performed using an orthonormalised version of m_basis, for stability.
         The scalar product for which we orthonormalise is a minimal one (much less expensive than integrating
         over the edge): \sum_i phi(x_i) psi(x_i) where (x_i)_i are nodes that are sufficient to determine
         entirely polynomials up to the considered degree.
         The nodes we build here correspond to equidistant P^k nodes on the edge */
 
        // Nodes 
        m_nb_nodes = basis.max_degree() + 1;
        m_nodes.reserve(m_nb_nodes);
        double h = 1./std::max(1., double(basis.max_degree()));
        for (size_t i=0; i<m_nb_nodes; i++){
          VectorRd xi = E.vertex(0)->coords() + double(i) * h * (E.vertex(1)->coords() - E.vertex(0)->coords());
          m_nodes.emplace_back(xi.x(), xi.y(), 1.); 
        }
 
        // We then orthonormalise Orthonormalise basis for simple scalar product
        m_on_basis_nodes.resize(boost::extents[m_dim][m_nb_nodes]);
        m_on_basis_nodes = evaluate_quad<Function>::compute(m_basis, m_nodes);
        m_on_basis = Family<BasisType>(l2_orthonormalize(m_basis, m_nodes, m_on_basis_nodes));
      };
      
    DecomposePoly(
            const Cell &T,    
            const BasisType &basis  
            ):
          m_dim(basis.dimension()),
          m_basis(basis),
          m_on_basis(basis, Eigen::MatrixXd::Identity(m_dim,m_dim))  // Need to initialise before we compute the real ON basis
      {
        /* The decomposition will be performed using an orthonormalised version of m_basis, for stability.
         The scalar product for which we orthonormalise is a minimal one (much less expensive than integrating
         over the cell): \sum_i phi(x_i) psi(x_i) where (x_i)_i are nodes that are sufficient to determine
         entirely polynomials up to the considered degree.
         The nodes we build here correspond to P^k nodes on a tetrahedron */
        
        // Create nodes
        std::vector<Eigen::Vector2i> indices = MonomialPowers<Cell>::complete(basis.max_degree());
        VectorRd e1 = VectorRd(1., 0.);
        VectorRd e2 = VectorRd(0., 1.);
        VectorRd x0 = T.center_mass() - T.diam()*VectorRd(1., 1.)/2.0;
        
        // Nodes 
        m_nb_nodes = indices.size();
        m_nodes.reserve(m_nb_nodes);
        for (size_t i=0; i<m_nb_nodes; i++){
          VectorRd xi = x0 + T.diam()*(double(indices[i](0)) * e1 + double(indices[i](1)) * e2)/std::max(1.,double(basis.max_degree()));
          m_nodes.emplace_back(xi.x(), xi.y(), 1.); 
        }
 
        // We then orthonormalise Orthonormalise basis for simple scalar product
        m_on_basis_nodes.resize(boost::extents[m_dim][m_nb_nodes]);
        m_on_basis_nodes = evaluate_quad<Function>::compute(m_basis, m_nodes);
        m_on_basis = Family<BasisType>(l2_orthonormalize(m_basis, m_nodes, m_on_basis_nodes));
      };
 
    inline QuadratureRule get_nodes() const 
    {
      return m_nodes;    
    };
 
    Family<BasisType> family(
            boost::multi_array<typename BasisType::FunctionValue, 2> &values  
            )
    {
      // Gram matrix of the polynomials and the ON basis
      Eigen::MatrixXd Gram_P_onbasis = compute_gram_matrix(values, m_on_basis_nodes, m_nodes);
      // L=matrix of P as a family of the ON basis. In theory, Gram_onbasis is identity, but due to rounding errors it 
      // can be a bit different. Computing L as if it's not the case improves stability 
      Eigen::MatrixXd Gram_onbasis = compute_gram_matrix(m_on_basis_nodes, m_nodes);
      Eigen::MatrixXd L = Gram_onbasis.ldlt().solve(Gram_P_onbasis.transpose());
      return Family<BasisType>(m_basis, L.transpose() * m_on_basis.matrix());
    };
 
 
    // Member variables
    size_t m_dim;   // dimension of basis
    BasisType m_basis; // Basis on which we decompose
    Family<BasisType> m_on_basis; // Orthonormalised basis (for simple scalar product)
    boost::multi_array<typename BasisType::FunctionValue, 2> m_on_basis_nodes; // Values of ON basis at the nodes
    size_t m_nb_nodes;    // nb of nodes
    QuadratureRule m_nodes; 
    
  };
 
  
 
} // end of namespace HArDCore2D
 
#endif