product_fe_space.h

#ifndef PROJECTS_DPG_PRODUCT_FE_SPACE_H
#define PROJECTS_DPG_PRODUCT_FE_SPACE_H
 
#include <lf/base/base.h>
#include <lf/fe/fe.h>
#include <lf/uscalfe/uscalfe.h>
 
#include <vector>
 
#include "dpg.h"
#include "product_dofhandler.h"
 
namespace projects::dpg {
 
template <typename SCALAR>
class ProductUniformFESpace {
 public:
  using scalar = SCALAR;
 
  ProductUniformFESpace() = default;
  ProductUniformFESpace(const ProductUniformFESpace &) = delete;
  ProductUniformFESpace(ProductUniformFESpace &&) noexcept = default;
  ProductUniformFESpace &operator=(const ProductUniformFESpace &) = delete;
  ProductUniformFESpace &operator=(ProductUniformFESpace &&) noexcept = default;
 
  ProductUniformFESpace(
      std::shared_ptr<const lf::mesh::Mesh> mesh_p,
      std::vector<
          std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
          rfs_tria_v,
      std::vector<
          std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
          rfs_quad_v,
      std::vector<
          std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
          rfs_edge_v,
      std::vector<
          std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
          rfs_point_v)
      : mesh_p_(std::move(mesh_p)),
        rfs_tria_v_(std::move(rfs_tria_v)),
        rfs_quad_v_(std::move(rfs_quad_v)),
        rfs_edge_v_(std::move(rfs_edge_v)),
        rfs_point_v_(std::move(rfs_point_v)) {
    numComponents_ = rfs_tria_v_.size();
    init();
  }
 
  [[nodiscard]] std::shared_ptr<const lf::mesh::Mesh> Mesh() const {
    LF_ASSERT_MSG(mesh_p_ != nullptr, "Invalid FE space, no mesh");
    return mesh_p_;
  }
 
  [[nodiscard]] const ProductUniformFEDofHandler &LocGlobMap() const {
    LF_ASSERT_MSG(mesh_p_ != nullptr, "Invalid FE space, no mesh");
    LF_ASSERT_MSG(dofh_p_ != nullptr, "Invalid FE space, no dofhandler");
    return *dofh_p_;
  }
 
  lf::fe::ScalarReferenceFiniteElement<SCALAR> const *ShapeFunctionLayout(
      lf::base::RefEl ref_el_type, size_type component) const;
 
  [[nodiscard]] size_type NumRefShapeFunctions(lf::base::RefEl ref_el_type,
                                               size_type component) const;
 
  [[nodiscard]] size_type NumComponents() const { return numComponents_; }
 
  std::shared_ptr<const lf::uscalfe::UniformScalarFESpace<SCALAR>>
  ComponentFESpace(size_type component) const {
    return std::make_shared<lf::uscalfe::UniformScalarFESpace<SCALAR>>(
        mesh_p_, rfs_tria_v_[component], rfs_quad_v_[component],
        rfs_edge_v_[component], rfs_point_v_[component]);
  }
 
  virtual ~ProductUniformFESpace() = default;
 
 private:
  std::shared_ptr<const lf::mesh::Mesh> mesh_p_;
 
  std::vector<
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
      rfs_tria_v_;
  std::vector<
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
      rfs_quad_v_;
  std::vector<
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
      rfs_edge_v_;
  std::vector<
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>>
      rfs_point_v_;
 
  std::vector<size_type> num_rsf_node_;
  std::vector<size_type> num_rsf_edge_;
  std::vector<size_type> num_rsf_tria_;
  std::vector<size_type> num_rsf_quad_;
 
  std::unique_ptr<ProductUniformFEDofHandler> dofh_p_;
 
  size_type numComponents_ = 0;
 
  void init();
};
 
// Initialization method.
template <typename SCALAR>
void ProductUniformFESpace<SCALAR>::init() {
  // check validity and consistency of mesh pointer:
  LF_ASSERT_MSG(mesh_p_ != nullptr, "No mesh");
  LF_ASSERT_MSG(mesh_p_->DimMesh() == 2, "Only 2D cells");
 
  // we do not check, wheter a specification of reference finite elements
  // for cells (triangles and quadrilaterals) is given, since this may not
  // be desired necessary for some components in  a dpg method (fluxes)
 
  // check validity and consistency of the number of components:
  LF_ASSERT_MSG(numComponents_ != 0, "0 components FE space.");
  LF_ASSERT_MSG(numComponents_ == rfs_tria_v_.size() &&
                    numComponents_ == rfs_quad_v_.size() &&
                    numComponents_ == rfs_edge_v_.size() &&
                    numComponents_ == rfs_point_v_.size(),
                "Missmatch in number of shape functions passed for different "
                "entity tpyes.");
 
  // resize vectors saving the number of loal shape functions for each
  // component:
  num_rsf_tria_.resize(numComponents_);
  num_rsf_quad_.resize(numComponents_);
  num_rsf_edge_.resize(numComponents_);
  num_rsf_node_.resize(numComponents_);
 
  // compatibility checks and initialization of number of shape functions:
  for (size_type component = 0; component < numComponents_; component++) {
    num_rsf_tria_[component] = 0;
    num_rsf_quad_[component] = 0;
    num_rsf_edge_[component] = 0;
    num_rsf_node_[component] = 0;
 
    if (rfs_tria_v_[component] != nullptr) {
      // probe reference finite element
      LF_ASSERT_MSG(rfs_tria_v_[component]->RefEl() == lf::base::RefEl::kTria(),
                    "wrong type for triangle in component" << component);
      // initialize number of shape functions associated to entitites.
      num_rsf_node_[component] =
          rfs_tria_v_[component]->NumRefShapeFunctions(2);
      num_rsf_edge_[component] =
          rfs_tria_v_[component]->NumRefShapeFunctions(1);
      num_rsf_tria_[component] =
          rfs_tria_v_[component]->NumRefShapeFunctions(0);
    }
    if (rfs_quad_v_[component] != nullptr) {
      // probe reference finite element for quads.
      LF_ASSERT_MSG(rfs_quad_v_[component]->RefEl() == lf::base::RefEl::kQuad(),
                    "Wrong type for quad in component " << component);
      // initialize number of shape functions associated to entitites.
      num_rsf_node_[component] =
          rfs_quad_v_[component]->NumRefShapeFunctions(2);
      num_rsf_edge_[component] =
          rfs_quad_v_[component]->NumRefShapeFunctions(1);
      num_rsf_quad_[component] =
          rfs_quad_v_[component]->NumRefShapeFunctions(0);
    }
    if (rfs_edge_v_[component] != nullptr) {
      LF_ASSERT_MSG(
          rfs_edge_v_[component]->RefEl() == lf::base::RefEl::kSegment(),
          "Wrong type for edge in component " << component);
      num_rsf_node_[component] =
          rfs_edge_v_[component]->NumRefShapeFunctions(1);
      num_rsf_edge_[component] =
          rfs_edge_v_[component]->NumRefShapeFunctions(0);
    }
 
    // Compatibility check for numbers of local shape functions associated with
    // edges. Those must be the same for all reference shape function
    // descriptions in a component passed to the finite element space.
    if ((rfs_tria_v_[component] != nullptr) &&
        (rfs_quad_v_[component] != nullptr)) {
      LF_ASSERT_MSG((rfs_tria_v_[component]->NumRefShapeFunctions(2)) ==
                        (rfs_quad_v_[component]->NumRefShapeFunctions(2)),
                    "#RSF missmatch on nodes in component"
                        << component << ": "
                        << rfs_tria_v_[component]->NumRefShapeFunctions(2)
                        << "<->"
                        << rfs_quad_v_[component]->NumRefShapeFunctions(2));
      LF_ASSERT_MSG((rfs_tria_v_[component]->NumRefShapeFunctions(1)) ==
                        (rfs_quad_v_[component]->NumRefShapeFunctions(1)),
                    "#RSF missmatch on edges in component"
                        << component << ": "
                        << rfs_tria_v_[component]->NumRefShapeFunctions(1)
                        << "<->"
                        << rfs_quad_v_[component]->NumRefShapeFunctions(1));
    }
    if ((rfs_tria_v_[component] != nullptr) &&
        (rfs_edge_v_[component] != nullptr)) {
      LF_ASSERT_MSG((rfs_tria_v_[component]->NumRefShapeFunctions(2)) ==
                        (rfs_edge_v_[component]->NumRefShapeFunctions(1)),
                    "#RSF missmatch on nodes in component"
                        << component << ": "
                        << rfs_tria_v_[component]->NumRefShapeFunctions(2)
                        << "<->"
                        << rfs_edge_v_[component]->NumRefShapeFunctions(1));
      LF_ASSERT_MSG((rfs_tria_v_[component]->NumRefShapeFunctions(1)) ==
                        (rfs_edge_v_[component]->NumRefShapeFunctions(0)),
                    "#RSF missmatch on edges in component"
                        << component << ": "
                        << rfs_tria_v_[component]->NumRefShapeFunctions(1)
                        << "<->"
                        << rfs_edge_v_[component]->NumRefShapeFunctions(0));
    }
    if ((rfs_quad_v_[component] != nullptr) &&
        (rfs_edge_v_[component] != nullptr)) {
      LF_ASSERT_MSG((rfs_quad_v_[component]->NumRefShapeFunctions(2)) ==
                        (rfs_edge_v_[component]->NumRefShapeFunctions(1)),
                    "#RSF missmatch on nodes in component"
                        << component << ": "
                        << rfs_quad_v_[component]->NumRefShapeFunctions(2)
                        << "<->"
                        << rfs_edge_v_[component]->NumRefShapeFunctions(1));
      LF_ASSERT_MSG((rfs_quad_v_[component]->NumRefShapeFunctions(1)) ==
                        (rfs_edge_v_[component]->NumRefShapeFunctions(0)),
                    "#RSF missmatch on edges in component"
                        << component << ": "
                        << rfs_quad_v_[component]->NumRefShapeFunctions(1)
                        << "<->"
                        << rfs_edge_v_[component]->NumRefShapeFunctions(0));
    }
  }
 
  // initialization of dof handler
  // collect the number of interior reference shape functions for all
  // components
  std::vector<dof_map_t> rsf_layouts(numComponents_);
  for (size_type component = 0; component < numComponents_; component++) {
    rsf_layouts[component] = {
        {lf::base::RefEl::kPoint(), num_rsf_node_[component]},
        {lf::base::RefEl::kSegment(), num_rsf_edge_[component]},
        {lf::base::RefEl::kTria(), num_rsf_tria_[component]},
        {lf::base::RefEl::kQuad(), num_rsf_quad_[component]}};
  }
  dofh_p_ = std::make_unique<ProductUniformFEDofHandler>(mesh_p_, rsf_layouts);
}
 
template <typename SCALAR>
lf::fe::ScalarReferenceFiniteElement<SCALAR> const *
ProductUniformFESpace<SCALAR>::ShapeFunctionLayout(lf::base::RefEl ref_el_type,
                                                   size_type component) const {
  // Retrive specification of local shape functions for
  // a certain component.
  switch (ref_el_type) {
    case lf::base::RefEl::kPoint(): {
      return rfs_point_v_[component].get();
    }
    case lf::base::RefEl::kSegment(): {
      // Null pointer valid return value:
      // indicates that a shape function description for edges is missing
      return rfs_edge_v_[component].get();
    }
    case lf::base::RefEl::kTria(): {
      // Null pointer valid return value:
      // indicates that a shape function description for triangles is missing.
      return rfs_tria_v_[component].get();
    }
    case lf::base::RefEl::kQuad(): {
      // Null pointer valid return value:
      // indicates that a shape function description for quadrilaterals is
      // missing.
      return rfs_quad_v_[component].get();
    }
    default: {
      LF_ASSERT_MSG(false, "Illegal entity type");
    }
  }
  return nullptr;
}
 
template <typename SCALAR>
size_type ProductUniformFESpace<SCALAR>::NumRefShapeFunctions(
    lf::base::RefEl ref_el_type, size_type component) const {
  // Retrieve number of interior shape functions from rsf layouts
  // for a certain component.
  switch (ref_el_type) {
    case lf::base::RefEl::kPoint(): {
      return num_rsf_node_[component];
    }
    case lf::base::RefEl::kSegment(): {
      return num_rsf_edge_[component];
    }
    case lf::base::RefEl::kTria(): {
      return num_rsf_tria_[component];
    }
    case lf::base::RefEl::kQuad(): {
      return num_rsf_quad_[component];
    }
    dafault : { LF_ASSERT_MSG(false, "Illegal entity type"); }
  }
  return 0;
}
 
}  // namespace projects::dpg
 
#endif  // PROJECTS_DPG_PRODUCT_FE_SPACE_H