vgrad.hpp

Go to the documentation of this file.

#ifndef VGRAD_HPP
#define VGRAD_HPP
 
#include <globaldofspace.hpp>
#include <vemcore.hpp>
#include <integralweight.hpp>
 
namespace HArDCore3D
{
 
  class VGrad : public GlobalDOFSpace
  {
  public:
    typedef std::function<double(const Eigen::Vector3d &)> FunctionType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> GradientType;
 
    struct LocalOperators
    {
      LocalOperators(
                     const Eigen::MatrixXd & _proj_gradient, 
                     const Eigen::MatrixXd & _dofs_gradient, 
                     const Eigen::MatrixXd & _proj_function 
                     )
        : proj_gradient(_proj_gradient),
          dofs_gradient(_dofs_gradient),
          proj_function(_proj_function)
      {
        // Do nothing
      }
      
      Eigen::MatrixXd proj_gradient;
      Eigen::MatrixXd dofs_gradient;
      Eigen::MatrixXd proj_function;
    };
    
    VGrad(const VEMCore & vem_core, bool use_threads = true, std::ostream & output = std::cout);
 
    const Mesh & mesh() const
    {
      return m_vem_core.mesh();
    }
    
    const size_t & degree() const
    {
      return m_vem_core.degree();
    }
    
    Eigen::VectorXd interpolate(
          const FunctionType & q, 
          const GradientType & grad_q, 
          const int doe_cell = -1, 
          const int doe_face = -1, 
          const int doe_edge = -1 
          ) const;
 
    inline const LocalOperators & edgeOperators(size_t iE) const
    {
      return *m_edge_operators[iE];
    }
 
    inline const LocalOperators & edgeOperators(const Edge & E) const
    {
      return *m_edge_operators[E.global_index()];
    }
    
    inline const LocalOperators & faceOperators(size_t iF) const
    {
      return *m_face_operators[iF];
    }
 
    inline const LocalOperators & faceOperators(const Face & F) const
    {
      return *m_face_operators[F.global_index()];
    }
 
    inline const LocalOperators & cellOperators(size_t iT) const
    {
      return *m_cell_operators[iT];
    }
 
    inline const LocalOperators & cellOperators(const Cell & T) const
    {
      return *m_cell_operators[T.global_index()];
    }
 
    inline const VEMCore::CellBases & cellBases(size_t iT) const
    {
      return m_vem_core.cellBases(iT);
    }
 
    inline const VEMCore::CellBases & cellBases(const Cell & T) const
    {
      return m_vem_core.cellBases(T.global_index());
    }
    
    inline const VEMCore::FaceBases & faceBases(size_t iF) const
    {
      return m_vem_core.faceBases(iF);
    }
 
    inline const VEMCore::FaceBases & faceBases(const Face & F) const
    {
      return m_vem_core.faceBases(F.global_index());
    }
 
    inline const VEMCore::EdgeBases & edgeBases(size_t iE) const
    {
      return m_vem_core.edgeBases(iE);
    }
 
    inline const VEMCore::EdgeBases & edgeBases(const Edge & E) const
    {
      return m_vem_core.edgeBases(E.global_index());
    }
 
    // The mass matrix of P^{k+1}(T) is the most expensive mass matrix in the calculation of this norm, which
    // is why there's the option of passing it as parameter if it's been already pre-computed when the norm is called.
    Eigen::MatrixXd computeL2Product(
                                     const size_t iT, 
                                     const double & penalty_factor = 1., 
                                     const Eigen::MatrixXd & mass_Pkmo_T = Eigen::MatrixXd::Zero(1,1), 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const;
 
    
  private:    
    LocalOperators _compute_edge_operators(size_t iE);
    LocalOperators _compute_face_operators(size_t iF);
    LocalOperators _compute_cell_operators(size_t iT);    
 
    const VEMCore & m_vem_core;
    bool m_use_threads;
    std::ostream & m_output;
 
    // Containers for local operators
    std::vector<std::unique_ptr<LocalOperators> > m_edge_operators;
    std::vector<std::unique_ptr<LocalOperators> > m_face_operators;
    std::vector<std::unique_ptr<LocalOperators> > m_cell_operators;
  };
 
} // end of namespace HArDCore3D
 
#endif