product_element_matrix_provider.h

#ifndef PROJECTS_DPG_PRODUCT_ELEMENT_MATRIX_PROVIDER
#define PROJECTS_DPG_PRODUCT_ELEMENT_MATRIX_PROVIDER
 
#include <lf/base/base.h>
#include <lf/quad/quad.h>
#include <lf/uscalfe/uscalfe.h>
 
#include <iostream>
#include <vector>
 
#include "dpg.h"
#include "product_fe_space.h"
#include "sub_element_matrix_provider.h"
 
namespace projects::dpg {
 
template <typename SCALAR>
class ProductElementMatrixProvider {
 public:
  using elem_mat_t = typename SubElementMatrixProvider<SCALAR>::elem_mat_t;
  using ElemMat = typename SubElementMatrixProvider<SCALAR>::ElemMat;
 
  ProductElementMatrixProvider(const ProductElementMatrixProvider&) = delete;
  ProductElementMatrixProvider(ProductElementMatrixProvider&&) noexcept =
      default;
  ProductElementMatrixProvider& operator=(const ProductElementMatrixProvider&) =
      delete;
  ProductElementMatrixProvider& operator=(ProductElementMatrixProvider&&) =
      delete;
 
  ProductElementMatrixProvider(
      std::shared_ptr<ProductUniformFESpace<SCALAR>> fe_space_trial,
      std::shared_ptr<ProductUniformFESpace<SCALAR>> fe_space_test,
      std::vector<std::shared_ptr<SubElementMatrixProvider<SCALAR>>>
          subproviders)
      : fe_space_trial_(std::move(fe_space_trial)),
        fe_space_test_(std::move(fe_space_test)),
        subproviders_(std::move(subproviders)) {}
 
  virtual bool isActive(const lf::mesh::Entity& /*cell*/) { return true; }
 
  ElemMat Eval(const lf::mesh::Entity& cell);
 
  virtual ~ProductElementMatrixProvider() = default;
 
 private:
  std::shared_ptr<ProductUniformFESpace<SCALAR>> fe_space_trial_;
  std::shared_ptr<ProductUniformFESpace<SCALAR>> fe_space_test_;
  std::vector<std::shared_ptr<SubElementMatrixProvider<SCALAR>>> subproviders_;
};
 
// template deduction hint
template <typename ptr, typename vector>
ProductElementMatrixProvider(ptr fe_space_tiral, ptr fe_space_test,
                             vector subproviders)
    -> ProductElementMatrixProvider<typename ptr::element_type::SCALAR>;
 
// evaluation method
template <typename SCALAR>
typename ProductElementMatrixProvider<SCALAR>::ElemMat
ProductElementMatrixProvider<SCALAR>::Eval(const lf::mesh::Entity& cell) {
  // retrive dofhandlers (for local assembly with the LocalStartIndex method
  LF_ASSERT_MSG(fe_space_trial_ != nullptr && fe_space_test_ != nullptr,
                "missing description of trial or test space");
  const ProductUniformFEDofHandler& dof_h_trial{fe_space_trial_->LocGlobMap()};
  const ProductUniformFEDofHandler& dof_h_test{fe_space_test_->LocGlobMap()};
 
  // initialize the element matrix:
  elem_mat_t mat(dof_h_test.NumLocalDofs(cell), dof_h_trial.NumLocalDofs(cell));
  mat.setZero();
 
  // build the element matrix from the "building block" sub element matrizes
  for (const auto& provider : subproviders_) {
    LF_ASSERT_MSG(provider != nullptr,
                  "invalid submatrix provider, is nullptr");
    if (provider->isActive(cell)) {
      // extract the involved components
      size_type trial_component = provider->TrialComponent();
      size_type test_component = provider->TestComponent();
      // evaluate the submatrix
      elem_mat_t sub_mat = provider->Eval(cell);
      // add it to the correct place in the element matrix.
      mat.block(dof_h_test.LocalStartIndex(cell, test_component),
                dof_h_trial.LocalStartIndex(cell, trial_component),
                sub_mat.rows(), sub_mat.cols()) += sub_mat;
    }
  }
  return mat;
}
 
}  // namespace projects::dpg
 
#endif  // PROJECTS_DPG_PRODUCT_ELEMENT_MATRIX_PROVIDER