GMpoly_face.hpp

Go to the documentation of this file.

# ifndef _GMPOLY_FACE_HPP
# define _GMPOLY_FACE_HPP
 
/*
 
  Classes to implement the quadrature-free polynomial integration rules
  
*/
 
#include <cassert>
#include <cmath>
 
#include <algorithm>
#include <array>
#include <iostream>
#include <vector>
#include <unordered_map>
 
#include <mesh.hpp>
#include <quadraturerule.hpp>
#include <basis.hpp>
#include <typeinfo>
 
#include <GMpoly_cell.hpp>
 
namespace HArDCore3D {
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
// INTEGRALS OF MONOMIALS
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
struct VecFaceHash
{
  size_t operator()(const Eigen::Vector2i & p) const
  {
    constexpr size_t b = 100; // max degree must be less than the basis b
    return p(0) + p(1) * b;
  }
};
 
typedef std::unordered_map<Eigen::Vector2i, double, VecFaceHash> MonomialFaceIntegralsType;
    
std::vector<MonomialFaceIntegralsType> IntegrateFaceMonomials_onEdges
          (const Face & F,      
          const size_t maxdeg   
          );
 
MonomialFaceIntegralsType IntegrateFaceMonomials
          (const Face & F,      
          const int maxdeg   
          );
          
MonomialFaceIntegralsType CheckIntegralsDegree
         (const Face & F,              
           const size_t degree,         
           const MonomialFaceIntegralsType & mono_int_map = {}    
          );
 
 
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//  FUNCTIONS TO COMPUTE GRAM MATRICES
//----------------------------------------------------------------------
//----------------------------------------------------------------------
 
//
/* BASIC ONES: Monomial, tensorized, and generic template */
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(const Face& F,         
                     const BasisType1 & basis1,   
                     const BasisType2 & basis2,   
                     MonomialFaceIntegralsType mono_int_map = {} 
                     )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(useAncestor<BasisType1>() || useAncestor<BasisType2>(), "No method to compute this Gram matrix");
      
    if constexpr (!useAncestor<BasisType1>() && useAncestor<BasisType2>()) {
      return transformGM(basis2, 'C', GramMatrix(F, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrix(F, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrix(F, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
 
  }
 
Eigen::MatrixXd GramMatrix(
                    const Face & F,                         
                    const MonomialScalarBasisFace & basis1, 
                    const MonomialScalarBasisFace & basis2, 
                    MonomialFaceIntegralsType mono_int_map = {} 
                    );
  
template<typename BasisType>
Eigen::MatrixXd GramMatrix(const Face& F, const BasisType & basis, MonomialFaceIntegralsType mono_int_map = {})
  {
    return GramMatrix(F, basis, basis, mono_int_map);
  }
 
template<typename BasisType1, typename BasisType2, size_t N>
Eigen::MatrixXd GramMatrix(
                     const Face& F, 
                     const TensorizedVectorFamily<BasisType1, N> & basis1, 
                     const TensorizedVectorFamily<BasisType2, N> & basis2,  
                     MonomialFaceIntegralsType mono_int_map = {} 
                     )
  {
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(basis1.dimension(), basis2.dimension());
 
    Eigen::MatrixXd anc_gm = GramMatrix(F, basis1.ancestor(), basis2.ancestor(), mono_int_map);
    size_t dim1 = anc_gm.rows();
    size_t dim2 = anc_gm.cols();
 
    for (size_t i=0; i<N; i++){
      gm.block(i*dim1, i*dim2, dim1, dim2) = anc_gm;    
    }
    return gm;
  }
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                             
                    const TangentFamily<BasisType1> & basis1,   
                    const TangentFamily<BasisType2> & basis2,   
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    Eigen::MatrixXd gm(basis1.dimension(), basis2.dimension());
    
    Eigen::MatrixXd anc_gm = GramMatrix(F, basis1.ancestor(), basis2.ancestor(), mono_int_map);
    size_t dim1 = anc_gm.rows();
    size_t dim2 = anc_gm.cols();
    
    Eigen::Matrix<double, 2, 3> generators1 = basis1.generators();
    Eigen::Matrix<double, 2, 3> generators2 = basis2.generators();
    
    for (size_t i=0; i<2; i++) {
      for (size_t j=0; j<2; j++) {
        gm.block(i*dim1, j*dim2, dim1, dim2) = generators1.row(i).dot(generators2.row(j)) * anc_gm;
      }
    }
    
    return gm;
  }
 
Eigen::MatrixXd GramMatrix(
                    const Face & F,                         
                    const RolyComplBasisFace & basis1,      
                    const RolyComplBasisFace & basis2,      
                    MonomialFaceIntegralsType mono_int_map = {} 
                    );
 
Eigen::MatrixXd GramMatrix(
                    const Face & F,                         
                    const GolyComplBasisFace & basis1,      
                    const GolyComplBasisFace & basis2,      
                    MonomialFaceIntegralsType mono_int_map = {} 
                    );
 
template<typename BasisType>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                             
                    const RolyComplBasisFace & basis1,          
                    const TangentFamily<BasisType> & basis2,    
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension()/2;
    size_t totaldegree = basis1.max_degree()+basis2.max_degree();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, 2*dim2);
    
    // Obtain integration data from IntegrateFaceMonomials
    MonomialFaceIntegralsType intmap = CheckIntegralsDegree(F, totaldegree, mono_int_map);
    
    std::vector<Eigen::MatrixXd> gmrs = {GMRolyComplScalar(F, basis1, basis2.ancestor(), 0, intmap), GMRolyComplScalar(F, basis1, basis2.ancestor(), 1, intmap)};
    Eigen::Matrix<double, 2, 3> jacobian1 = basis1.jacobian();
    Eigen::Matrix<double, 2, 3> generators2 = basis2.generators();
    Eigen::Matrix2d W = jacobian1 * generators2.transpose();
    
    for (size_t k=0; k<2; k++) {
      for (size_t j=0; j<2; j++) {
        gm.block(0, j*dim2, dim1, dim2) += W(k,j) * gmrs[k];
      }
    }
    
    return F.diam() * gm;
  }
 
template<typename BasisType>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                             
                    const TangentFamily<BasisType> & basis1,    
                    const RolyComplBasisFace & basis2,          
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    return GramMatrix(F, basis2, basis1, mono_int_map).transpose();
  }
 
Eigen::MatrixXd GMRolyComplScalar(
                    const Face & F, 
                    const RolyComplBasisFace & rolycompl_basis, 
                    const MonomialScalarBasisFace & mono_basis, 
                    const size_t m, 
                    MonomialFaceIntegralsType mono_int_map = {} 
                    );
 
template<typename BasisType>
Eigen::MatrixXd GMRolyComplScalar(
                     const Face & F, 
                     const RolyComplBasisFace & basis1, 
                     const BasisType & basis2,  
                     const size_t m, 
                     MonomialFaceIntegralsType mono_int_map = {} 
                     )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(BasisType::hasAncestor, "No method to compute this Gram matrix of derivatives");
    return transformGM(basis2, 'C', GMRolyComplScalar(F, basis1, basis2.ancestor(), m, mono_int_map) );
  }
  
Eigen::MatrixXd GMScalarDerivative(
                  const Face & F,  
                  const MonomialScalarBasisFace & basis1, 
                  const MonomialScalarBasisFace & basis2, 
                  const size_t m, 
                  MonomialFaceIntegralsType mono_int_map = {} 
                  );
  
Eigen::MatrixXd GMScalarDerivative(
                  const Face & F,  
                  const MonomialScalarBasisFace & basis1, 
                  const MonomialScalarBasisFace & basis2, 
                  const size_t m, 
                  const size_t l, 
                  MonomialFaceIntegralsType mono_int_map = {} 
                  );
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GMScalarDerivative(
                     const Face& F, 
                     const BasisType1 & basis1, 
                     const BasisType2 & basis2,  
                     const size_t m, 
                     MonomialFaceIntegralsType mono_int_map = {} 
                             )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(BasisType1::hasAncestor || BasisType2::hasAncestor, "No method to compute this Gram matrix of derivatives");
      
    if constexpr (!BasisType1::hasAncestor && BasisType2::hasAncestor) {
      return transformGM(basis2, 'C', GMScalarDerivative(F, basis1, basis2.ancestor(), m, mono_int_map) );
    } else if constexpr (BasisType1::hasAncestor && !BasisType2::hasAncestor) {
      return transformGM(basis1, 'R', GMScalarDerivative(F, basis1.ancestor(), basis2, m, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GMScalarDerivative(F, basis1.ancestor(), basis2.ancestor(), m, mono_int_map) ) );
    }
  }
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GMScalarDerivative(
                     const Face& F, 
                     const BasisType1 & basis1, 
                     const BasisType2 & basis2,  
                     const size_t m, 
                     const size_t l, 
                     MonomialFaceIntegralsType mono_int_map = {} 
                             )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(BasisType1::hasAncestor || BasisType2::hasAncestor, "No method to compute this Gram matrix of derivatives");
      
    if constexpr (!BasisType1::hasAncestor && BasisType2::hasAncestor) {
      return transformGM(basis2, 'C', GMScalarDerivative(F, basis1, basis2.ancestor(), m, l, mono_int_map) );
    } else if constexpr (BasisType1::hasAncestor && !BasisType2::hasAncestor) {
      return transformGM(basis1, 'R', GMScalarDerivative(F, basis1.ancestor(), basis2, m, l, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GMScalarDerivative(F, basis1.ancestor(), basis2.ancestor(), m, l, mono_int_map) ) );
    }
    
  }
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                             
                    const CurlBasis<BasisType1> & basis1,       
                    const CurlBasis<BasisType2> & basis2,       
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    // Dimension of the gram matrix
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    size_t totaldegree = basis1.ancestor().max_degree()+basis2.ancestor().max_degree()-2;
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1,dim2);
    
    // Obtain integration data from IntegrateFaceMonomials
    MonomialFaceIntegralsType intmap = CheckIntegralsDegree(F, totaldegree, mono_int_map);
    
    for (size_t m=0; m<2; m++){
      gm += GMScalarDerivative(F, basis1.ancestor(), basis2.ancestor(), m, m, intmap);
    }
    
    return gm/std::pow(F.diam(), 2);
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                             
                    const CurlBasis<BasisType1> & basis1,       
                    const TangentFamily<BasisType2> & basis2,   
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    // Dimension of the gram matrix
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension()/2;
    size_t totaldegree = basis1.ancestor().max_degree()+basis2.max_degree()-1;
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1,dim2*2);
    
    // Obtain integration data from IntegrateFaceMonomials
    MonomialFaceIntegralsType intmap = CheckIntegralsDegree(F, totaldegree, mono_int_map);
    
    std::vector<Eigen::MatrixXd> gmd = {GMScalarDerivative(F, basis1.ancestor(), basis2.ancestor(), 0, intmap), GMScalarDerivative(F, basis1.ancestor(), basis2.ancestor(), 1, intmap)};
    // JF must be the jacobian (change of variable 3D->2D) of a scalar basis on the face (underlying basis1, perhaps after
    // several derived bases). We create one to make sure we capture the correct jacobian
    MonomialScalarBasisFace mono_P0_F(F, 0);
    Eigen::Matrix<double, 2, 3> JF = mono_P0_F.jacobian();
    VectorRd nF = F.normal();
    Eigen::Matrix<double, 2, 3> generators2 = basis2.generators();
    Eigen::Matrix2d W;
    W.col(0) = JF * nF.cross(generators2.row(0));
    W.col(1) = JF * nF.cross(generators2.row(1));
    
    for (size_t k=0; k<2; k++) {
      for (size_t j=0; j<2; j++) {
        gm.block(0, j*dim2, dim1, dim2) += W(k,j) * gmd[k];
      }
    }
    
    return gm;
  }
 
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                             
                    const TangentFamily<BasisType1> & basis1,   
                    const CurlBasis<BasisType2> & basis2,       
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    return GramMatrix(F, basis2, basis1, mono_int_map).transpose();
  }
 
 
template<typename BasisType1, typename BasisType2>
typename boost::disable_if<boost::is_same<BasisType2, MonomialFaceIntegralsType>, Eigen::MatrixXd>::type GramMatrix(
                     const Face& F,                          
                     const DivergenceBasis<BasisType1> & basis1,   
                     const BasisType2 & basis2,              
                     MonomialFaceIntegralsType mono_int_map = {} 
                     )
  {
    return GramMatrixDiv(F, basis1.ancestor(), basis2, mono_int_map);
  }
  
template<typename BasisType1, typename Basis2>
Eigen::MatrixXd GramMatrix(
                    const Face & F,                         
                    const BasisType1 & basis1,                  
                    const DivergenceBasis<Basis2> & basis2,       
                    MonomialFaceIntegralsType mono_int_map = {} 
                    )
  {
    return GramMatrixDiv(F, basis2.ancestor(), basis1, mono_int_map).transpose();
  }
  
/* This will work only if the generators of the TangentFamily are an orthonormal version of the Jacobian (change of variable 3D->2D) on the face */
template<typename BasisType1>
Eigen::MatrixXd GramMatrixDiv(
                    const Face & F,                                         
                    const TangentFamily<BasisType1> & basis1,   
                    const MonomialScalarBasisFace & basis2,                 
                    MonomialFaceIntegralsType mono_int_map = {}                 
                    )
  {
    size_t dim1 = basis1.dimension();
    size_t dim2 = basis2.dimension();
    Eigen::MatrixXd gm = Eigen::MatrixXd::Zero(dim1, dim2);
 
    // Integrals of monomials
    size_t totaldegree = basis1.max_degree()+basis2.max_degree()-1;
    MonomialFaceIntegralsType intmap = CheckIntegralsDegree(F, totaldegree, mono_int_map);
 
    for (size_t i = 0; i < 2; i++) {
      gm.block(i*dim1/2, 0, dim1/2, dim2) = GMScalarDerivative(F, basis1.ancestor(), basis2, i, intmap);
    }
    
    return gm/F.diam();
  }
  
Eigen::MatrixXd GramMatrixDiv(
                    const Face & F,                           
                    const RolyComplBasisFace & basis1,        
                    const MonomialScalarBasisFace & basis2,   
                    MonomialFaceIntegralsType mono_int_map = {}   
                    );
  
template<typename BasisType1, typename BasisType2>
Eigen::MatrixXd GramMatrixDiv(
                     const Face& F, 
                     const BasisType1 & basis1, 
                     const BasisType2 & basis2,  
                     MonomialFaceIntegralsType mono_int_map = {} 
                     )
  {
    // If no ancestor is to be used, we shouldn't be in this overload
    static_assert(useAncestor<BasisType1>() || useAncestor<BasisType2>(), "No method to compute this Gram matrix");
      
    if constexpr (!useAncestor<BasisType1>() && useAncestor<BasisType2>()) {
      return transformGM(basis2, 'C', GramMatrixDiv(F, basis1, basis2.ancestor(), mono_int_map) );
    } else if constexpr (useAncestor<BasisType1>() && !useAncestor<BasisType2>()) {
      return transformGM(basis1, 'R', GramMatrixDiv(F, basis1.ancestor(), basis2, mono_int_map) );
    } else {
      return transformGM(basis1, 'R', transformGM(basis2, 'C', GramMatrixDiv(F, basis1.ancestor(), basis2.ancestor(), mono_int_map) ) );
    }
  }
 
 
 
 
} // end of namespace HArDCore3D
 
#endif // end of _GMPOLY_FACE_HPP