dpg_element_matrix_provider.h

#ifndef PROJECTS_DPG_DPG_ELEMENT_MATRIX_PROVIDER
#define PROJECTS_DPG_DPG_ELEMENT_MATRIX_PROVIDER
 
#include <lf/mesh/mesh.h>
 
#include "dpg.h"
#include "product_element_matrix_provider.h"
 
namespace projects::dpg {
 
template <typename SCALAR>
class DpgElementMatrixProvider {
 public:
  using elem_mat_t = typename ProductElementMatrixProvider<SCALAR>::elem_mat_t;
  using ElemMat = typename ProductElementMatrixProvider<SCALAR>::ElemMat;
 
  DpgElementMatrixProvider(const DpgElementMatrixProvider&) = delete;
  DpgElementMatrixProvider(DpgElementMatrixProvider&&) noexcept = default;
  DpgElementMatrixProvider& operator=(const DpgElementMatrixProvider&) = delete;
  DpgElementMatrixProvider& operator=(DpgElementMatrixProvider&&) = delete;
 
  DpgElementMatrixProvider(
      std::shared_ptr<ProductElementMatrixProvider<SCALAR>>
          extendedStiffnessMatrixProvider,
      std::shared_ptr<ProductElementMatrixProvider<SCALAR>> gramianProvider)
      : extendedStiffnessMatrixProvider_(
            std::move(extendedStiffnessMatrixProvider)),
        gramianProvider_(std::move(gramianProvider)) {}
 
  virtual bool isActive(const lf::mesh::Entity& /*cell*/) { return true; }
 
  ElemMat Eval(const lf::mesh::Entity& cell);
 
  virtual ~DpgElementMatrixProvider() = default;
 
 private:
  std::shared_ptr<ProductElementMatrixProvider<SCALAR>>
      extendedStiffnessMatrixProvider_;
  std::shared_ptr<ProductElementMatrixProvider<SCALAR>> gramianProvider_;
};
 
// template deduction hint
template <class PTR>
DpgElementMatrixProvider(PTR stiffness_provider, PTR gramianProvider)
    -> DpgElementMatrixProvider<typename PTR::element_type::SCALAR>;
 
// evaluation method.
template <typename SCALAR>
typename DpgElementMatrixProvider<SCALAR>::ElemMat
DpgElementMatrixProvider<SCALAR>::Eval(const lf::mesh::Entity& cell) {
  // check for nullptrs
  LF_ASSERT_MSG(extendedStiffnessMatrixProvider_ != nullptr &&
                    gramianProvider_ != nullptr,
                "nullptr error for some provider");
  // check that the method is called on an active cell
  LF_ASSERT_MSG(gramianProvider_->isActive(cell) &&
                    extendedStiffnessMatrixProvider_->isActive(cell),
                "Eval method called on inactive cell. " << cell);
 
  // evaluate extended stiffness matrix B and local Gramian G
  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");
  LF_ASSERT_MSG(
      extendedStiffnessMatrix.rows() == gramian.cols(),
      "size missmatch between gramian & extended Stiffness matrix on cell "
          << cell);
 
  // evaluate DPG element matrix A = B^T G^-1 B
  return extendedStiffnessMatrix.transpose() *
         gramian.ldlt().solve(extendedStiffnessMatrix);
}
 
}  // namespace projects::dpg
 
#endif  // PROJECTS_DPG_ELEMENT_MATRIX_PROVIDER