dpg_element_vector_provider.h

#ifndef PROJECTS_DPG_DPG_ELEMENT_VECTOR_PROVIDER
#define PROJECTS_DPG_DPG_ELEMENT_VECTOR_PROVIDER
 
#include <lf/mesh/mesh.h>
 
#include "dpg.h"
#include "product_element_matrix_provider.h"
#include "product_element_vector_provider.h"
 
namespace projects::dpg {
template <typename SCALAR>
class DpgElementVectorProvider {
 public:
  using elem_vec_t = typename ProductElementVectorProvider<SCALAR>::elem_vec_t;
  using ElemVec = typename ProductElementVectorProvider<SCALAR>::ElemVec;
  using ElemMat = typename ProductElementMatrixProvider<SCALAR>::ElemMat;
 
  DpgElementVectorProvider(const DpgElementVectorProvider&) = delete;
  DpgElementVectorProvider(DpgElementVectorProvider&&) noexcept = default;
  DpgElementVectorProvider& operator=(const DpgElementVectorProvider&) = delete;
  DpgElementVectorProvider& operator=(DpgElementVectorProvider&&) = delete;
 
  DpgElementVectorProvider(
      std::shared_ptr<ProductElementVectorProvider<SCALAR>>
          extendedLoadVectorProvider,
      std::shared_ptr<ProductElementMatrixProvider<SCALAR>>
          extendedStiffnessMatrixProvider,
      std::shared_ptr<ProductElementMatrixProvider<SCALAR>> gramianProvider)
      : extendedLoadVectorProvider_(std::move(extendedLoadVectorProvider)),
        extendedStiffnessMatrixProvider_(
            std::move(extendedStiffnessMatrixProvider)),
        gramianProvider_(std::move(gramianProvider)) {}
 
  virtual bool isActive(const lf::mesh::Entity& /*cell*/) { return true; }
 
  ElemVec Eval(const lf::mesh::Entity& cell);
 
  virtual ~DpgElementVectorProvider() = default;
 
 private:
  std::shared_ptr<ProductElementVectorProvider<SCALAR>>
      extendedLoadVectorProvider_;
  std::shared_ptr<ProductElementMatrixProvider<SCALAR>>
      extendedStiffnessMatrixProvider_;
  std::shared_ptr<ProductElementMatrixProvider<SCALAR>> gramianProvider_;
};
 
// template deduction hint:
template <class PTR1, class PTR2>
DpgElementVectorProvider(PTR1 eLoadVectorProvider, PTR2 eStiffnessProvider,
                         PTR2 gramianProvider)
    -> DpgElementVectorProvider<typename PTR1::element_type::SCALAR>;
 
// evaluation method
template <typename SCALAR>
typename DpgElementVectorProvider<SCALAR>::ElemVec
DpgElementVectorProvider<SCALAR>::Eval(const lf::mesh::Entity& cell) {
  // check for nullptrs
  LF_ASSERT_MSG(extendedLoadVectorProvider_ != nullptr &&
                    extendedStiffnessMatrixProvider_ != nullptr &&
                    gramianProvider_ != nullptr,
                "nullptr error for some provider");
  LF_ASSERT_MSG(extendedLoadVectorProvider_->isActive(cell) &&
                    extendedStiffnessMatrixProvider_->isActive(cell) &&
                    gramianProvider_->isActive(cell),
                "Eval method called on inactive cell");
 
  // evaluate extended element vector l, extended stiffness matrix B and local
  // Gramian G.
  ElemVec extendedLoadVector = extendedLoadVectorProvider_->Eval(cell);
  ElemMat extendedStiffnessMatrix =
      extendedStiffnessMatrixProvider_->Eval(cell);
  ElemMat gramian = gramianProvider_->Eval(cell);
 
  // perform some size checks.
  LF_ASSERT_MSG(gramian.rows() == gramian.cols(),
                "non quadratic gramian of size ("
                    << gramian.rows() << ", " << gramian.cols() << ") on cell "
                    << cell << "\n");
 
  // evaluate element vector B^T G^-1 l
  return extendedStiffnessMatrix.transpose() *
         gramian.ldlt().solve(extendedLoadVector);
}
}  // namespace projects::dpg
#endif  // PROJECTS_DPG_DPG_ELEMENt_VECTOR_PROVIDER