product_element_vector_provider.h

#ifndef PROJECTS_DPG_PRODUCT_ELEMENT_VECTOR_PROVIDER
#define PROJECTS_DPG_PRODUCT_ELEMENT_VECTOR_PROVIDER
 
#include <lf/mesh/mesh.h>
#include <lf/uscalfe/uscalfe.h>
 
#include <iostream>
#include <vector>
 
#include "dpg.h"
#include "product_fe_space.h"
#include "sub_element_vector_provider.h"
 
namespace projects::dpg {
 
using size_type = lf::uscalfe::size_type;
 
template <typename SCALAR>
class ProductElementVectorProvider {
 public:
  using elem_vec_t = typename SubElementVectorProvider<SCALAR>::elem_vec_t;
  using ElemVec = typename SubElementVectorProvider<SCALAR>::ElemVec;
 
  ProductElementVectorProvider(const ProductElementVectorProvider&) = delete;
  ProductElementVectorProvider(ProductElementVectorProvider&&) noexcept =
      default;
  ProductElementVectorProvider& operator=(const ProductElementVectorProvider&) =
      delete;
  ProductElementVectorProvider& operator=(ProductElementVectorProvider&&) =
      delete;
 
  ProductElementVectorProvider(
      std::shared_ptr<ProductUniformFESpace<SCALAR>> fe_space_test,
      std::vector<std::shared_ptr<SubElementVectorProvider<SCALAR>>>
          subproviders)
      : fe_space_test_(std::move(fe_space_test)),
        subproviders_(std::move(subproviders)) {}
 
  virtual bool isActive(const lf::mesh::Entity& /*cell*/) { return true; }
 
  ElemVec Eval(const lf::mesh::Entity& cell);
 
  virtual ~ProductElementVectorProvider() = default;
 
 private:
  std::shared_ptr<ProductUniformFESpace<SCALAR>> fe_space_test_;
  std::vector<std::shared_ptr<SubElementVectorProvider<SCALAR>>> subproviders_;
};
 
// template deduction hint
template <typename PTR, typename vector>
ProductElementVectorProvider(PTR fe_space_test, vector subproviders)
    -> ProductElementVectorProvider<typename PTR::element_type::SCALAR>;
 
// evaluation method
template <typename SCALAR>
typename ProductElementVectorProvider<SCALAR>::ElemVec
ProductElementVectorProvider<SCALAR>::Eval(const lf::mesh::Entity& cell) {
  // retrive dofhandler(for local assembly with LocalStartIndex method
  LF_ASSERT_MSG(fe_space_test_ != nullptr,
                "missing description of the trial space");
  const ProductUniformFEDofHandler& dof_h{fe_space_test_->LocGlobMap()};
 
  // initialize the element vector
  elem_vec_t vec(dof_h.NumLocalDofs(cell));
  vec.setZero();
 
  // build the building block from the "building block" subproviders
  for (const auto& provider : subproviders_) {
    LF_ASSERT_MSG(provider != nullptr,
                  "invalid subvector provider, is nullptr");
    if (provider->isActive(cell)) {
      // extract involved component
      size_type test_component = provider->TestComponent();
      // evaluate the subvector
      elem_vec_t sub_vec = provider->Eval(cell);
      // add it to the corresponding place in the element vector.
      vec.segment(dof_h.LocalStartIndex(cell, test_component),
                  sub_vec.size()) += sub_vec;
    }
  }
  return vec;
}
 
}  // namespace projects::dpg
 
#endif  // PROJECTS_DPG_PRODUCT_ELEMENT_VECTOR_PROVIDER.