uniform_scalar_fe_space.h

#ifndef LF_FESPACE_H
#define LF_FESPACE_H
/***************************************************************************
 * LehrFEM++ - A simple C++ finite element libray for teaching
 * Developed from 2018 at the Seminar of Applied Mathematics of ETH Zurich,
 * lead developers Dr. R. Casagrande and Prof. R. Hiptmair
 ***************************************************************************/
 
#include <lf/assemble/assemble.h>
#include <lf/fe/scalar_fe_space.h>
 
#include "lagr_fe.h"
 
namespace lf::uscalfe {
 
template <typename SCALAR>
class UniformScalarFESpace : public lf::fe::ScalarFESpace<SCALAR> {
 public:
  using Scalar = SCALAR;
 
  UniformScalarFESpace(const UniformScalarFESpace &) = delete;
  UniformScalarFESpace(UniformScalarFESpace &&) noexcept = default;
  UniformScalarFESpace &operator=(const UniformScalarFESpace &) = delete;
  UniformScalarFESpace &operator=(UniformScalarFESpace &&) noexcept = default;
  UniformScalarFESpace(
      std::shared_ptr<const lf::mesh::Mesh> mesh_p,
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
          rfs_tria_p,
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
          rfs_quad_p,
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
          rfs_edge_p = nullptr,
      std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
          rfs_point_p = nullptr)
      : lf::fe::ScalarFESpace<SCALAR>(),
        rfs_tria_p_(std::move(rfs_tria_p)),
        rfs_quad_p_(std::move(rfs_quad_p)),
        rfs_edge_p_(std::move(rfs_edge_p)),
        rfs_point_p_(std::move(rfs_point_p)),
        dofh_(InitDofHandler(std::move(mesh_p))) {}
 
  [[nodiscard]] std::shared_ptr<const lf::mesh::Mesh> Mesh() const override {
    return dofh_.Mesh();
  }
 
  [[nodiscard]] const lf::assemble::DofHandler &LocGlobMap() const override {
    LF_VERIFY_MSG(Mesh() != nullptr, "No valid FE space object: no mesh");
    return dofh_;
  }
 
  [[nodiscard]] lf::fe::ScalarReferenceFiniteElement<SCALAR> const *
  ShapeFunctionLayout(const lf::mesh::Entity &entity) const override;
 
  [[nodiscard]] lf::fe::ScalarReferenceFiniteElement<SCALAR> const *
  ShapeFunctionLayout(lf::base::RefEl ref_el_type) const;
 
  [[nodiscard]] size_type NumRefShapeFunctions(
      const lf::mesh::Entity &entity) const override;
 
  [[nodiscard]] size_type NumRefShapeFunctions(
      lf::base::RefEl ref_el_type) const;
 
  ~UniformScalarFESpace() override = default;
 
 private:
  std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
      rfs_tria_p_;
  std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
      rfs_quad_p_;
  std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
      rfs_edge_p_;
  std::shared_ptr<const lf::fe::ScalarReferenceFiniteElement<SCALAR>>
      rfs_point_p_;
  size_type num_rsf_node_{0}, num_rsf_edge_{0}, num_rsf_tria_{0},
      num_rsf_quad_{0};
 
  assemble::UniformFEDofHandler dofh_;
 
  assemble::UniformFEDofHandler InitDofHandler(
      std::shared_ptr<const lf::mesh::Mesh> mesh_p);
  [[nodiscard]] bool check_ptr() const {
    LF_VERIFY_MSG(Mesh() != nullptr, "No valid FE space object: no mesh");
    LF_VERIFY_MSG((rfs_quad_p_ != nullptr) && (rfs_quad_p_ != nullptr),
                  "No valid FE space object: no rsfs for cells");
    return true;
  }
 
};  // end class definition UniformScalarFESpace
 
const unsigned int kUniformScalarFESpaceOutMesh = 1;
 
const unsigned int kUniformScalarFESpaceOutDofh = 2;
 
const unsigned int kUniformScalarFESpaceOutRsfs = 4;
 
template <class SCALAR>
void PrintInfo(std::ostream &o, const UniformScalarFESpace<SCALAR> &fes,
               unsigned int ctrl = 0);
 
template <typename SCALAR>
std::ostream &operator<<(std::ostream &o,
                         const UniformScalarFESpace<SCALAR> &fes);
 
// Initialization methods
template <typename SCALAR>
assemble::UniformFEDofHandler UniformScalarFESpace<SCALAR>::InitDofHandler(
    std::shared_ptr<const lf::mesh::Mesh> mesh_p) {
  // Check validity and consistency of mesh pointer
  LF_VERIFY_MSG(mesh_p != nullptr, "Missing mesh!");
  LF_VERIFY_MSG((rfs_quad_p_ != nullptr) || (rfs_tria_p_ != nullptr),
                "Missing FE specification for cells");
  LF_VERIFY_MSG((mesh_p->DimMesh() == 2), "Only for 2D meshes");
 
  // Check whether all required finite element specifications are provided
  LF_VERIFY_MSG(
      (mesh_p->NumEntities(lf::base::RefEl::kTria()) == 0) ||
          (rfs_tria_p_ != nullptr),
      "Missing FE specification for triangles though mesh contains some");
  LF_VERIFY_MSG((mesh_p->NumEntities(lf::base::RefEl::kQuad()) == 0) ||
                    (rfs_quad_p_ != nullptr),
                "Missing FE specification for quads though mesh contains some");
 
  // Compatibility checks and initialization of numbers of shape functions
  // In particular only a single shape function may be associated to a node
  // in the case of a SCALAR finite element space
  if (rfs_tria_p_ != nullptr) {
    // Probe local shape functions on a triangle
    LF_VERIFY_MSG((*rfs_tria_p_).RefEl() == lf::base::RefEl::kTria(),
                  "Wrong type for triangle!");
    LF_VERIFY_MSG((*rfs_tria_p_).NumRefShapeFunctions(2) <= 1,
                  "At most one shape function can be assigned to each vertex");
    // Initialize numbers of shape functions associated to entities
    num_rsf_node_ = (*rfs_tria_p_).NumRefShapeFunctions(2);
    num_rsf_edge_ = (*rfs_tria_p_).NumRefShapeFunctions(1);
    num_rsf_tria_ = (*rfs_tria_p_).NumRefShapeFunctions(0);
  }
  if (rfs_quad_p_ != nullptr) {
    // Probe local shape functions for QUADs
    LF_VERIFY_MSG((*rfs_quad_p_).RefEl() == lf::base::RefEl::kQuad(),
                  "Wrong type for quad!");
    LF_VERIFY_MSG((*rfs_quad_p_).NumRefShapeFunctions(2) <= 1,
                  "At most one shape function can be assigned to each vertex");
    // Initialize numbers of shape functions associated to entities
    num_rsf_node_ = (*rfs_quad_p_).NumRefShapeFunctions(2);
    num_rsf_edge_ = (*rfs_quad_p_).NumRefShapeFunctions(1);
    num_rsf_quad_ = (*rfs_quad_p_).NumRefShapeFunctions(0);
  }
  if (rfs_edge_p_ != nullptr) {
    LF_VERIFY_MSG((*rfs_edge_p_).RefEl() == lf::base::RefEl::kSegment(),
                  "Wrong type for edge!");
    LF_VERIFY_MSG((*rfs_edge_p_).NumRefShapeFunctions(1) <= 1,
                  "At most one shape function can be assigned to each vertex");
    num_rsf_node_ = (*rfs_edge_p_).NumRefShapeFunctions(1);
    num_rsf_edge_ = (*rfs_edge_p_).NumRefShapeFunctions(0);
  }
 
  // Compatibility check for numbers of local shape functions associated with
  // edges. Those must be the same for all reference shape function descriptions
  // passed to the finite element space.
  if ((rfs_tria_p_ != nullptr) && (rfs_quad_p_ != nullptr)) {
    LF_ASSERT_MSG(((*rfs_tria_p_).NumRefShapeFunctions(2) ==
                   (*rfs_quad_p_).NumRefShapeFunctions(2)),
                  "#RSF mismatch on nodes "
                      << (*rfs_tria_p_).NumRefShapeFunctions(2) << " <-> "
                      << (*rfs_quad_p_).NumRefShapeFunctions(2));
    LF_ASSERT_MSG(((*rfs_tria_p_).NumRefShapeFunctions(1) ==
                   (*rfs_quad_p_).NumRefShapeFunctions(1)),
                  "#RSF mismatch on edges "
                      << (*rfs_tria_p_).NumRefShapeFunctions(1) << " <-> "
                      << (*rfs_quad_p_).NumRefShapeFunctions(1));
  }
  if ((rfs_tria_p_ != nullptr) && (rfs_edge_p_ != nullptr)) {
    LF_ASSERT_MSG(((*rfs_tria_p_).NumRefShapeFunctions(2) ==
                   (*rfs_edge_p_).NumRefShapeFunctions(1)),
                  "#RSF mismatch on nodes "
                      << (*rfs_tria_p_).NumRefShapeFunctions(2) << " <-> "
                      << (*rfs_edge_p_).NumRefShapeFunctions(1));
    LF_ASSERT_MSG(((*rfs_tria_p_).NumRefShapeFunctions(1) ==
                   (*rfs_edge_p_).NumRefShapeFunctions(0)),
                  "#RSF mismatch on edges "
                      << (*rfs_tria_p_).NumRefShapeFunctions(1) << " <-> "
                      << (*rfs_edge_p_).NumRefShapeFunctions(0));
  }
  if ((rfs_quad_p_ != nullptr) && (rfs_edge_p_ != nullptr)) {
    LF_ASSERT_MSG(((*rfs_quad_p_).NumRefShapeFunctions(2) ==
                   (*rfs_edge_p_).NumRefShapeFunctions(1)),
                  "#RSF mismatch on edges "
                      << (*rfs_quad_p_).NumRefShapeFunctions(2) << " <-> "
                      << (*rfs_edge_p_).NumRefShapeFunctions(1));
    LF_ASSERT_MSG(((*rfs_quad_p_).NumRefShapeFunctions(1) ==
                   (*rfs_edge_p_).NumRefShapeFunctions(0)),
                  "#RSF mismatch on edges "
                      << (*rfs_quad_p_).NumRefShapeFunctions(1) << " <-> "
                      << (*rfs_edge_p_).NumRefShapeFunctions(0));
  }
 
  // Initialization of dof handler starting with collecting the number of
  // interior reference shape functions
  lf::assemble::UniformFEDofHandler::dof_map_t rsf_layout{
      {lf::base::RefEl::kPoint(), num_rsf_node_},
      {lf::base::RefEl::kSegment(), num_rsf_edge_},
      {lf::base::RefEl::kTria(), num_rsf_tria_},
      {lf::base::RefEl::kQuad(), num_rsf_quad_}};
  return lf::assemble::UniformFEDofHandler(std::move(mesh_p), rsf_layout);
}
 
template <typename SCALAR>
lf::fe::ScalarReferenceFiniteElement<SCALAR> const *
UniformScalarFESpace<SCALAR>::ShapeFunctionLayout(
    const lf::mesh::Entity &entity) const {
  return ShapeFunctionLayout(entity.RefEl());
}
 
template <typename SCALAR>
lf::fe::ScalarReferenceFiniteElement<SCALAR> const *
UniformScalarFESpace<SCALAR>::ShapeFunctionLayout(
    lf::base::RefEl ref_el_type) const {
  // Retrieve specification of local shape functions
  switch (ref_el_type) {
    case lf::base::RefEl::kPoint(): {
      return rfs_point_p_.get();
    }
    case lf::base::RefEl::kSegment(): {
      // Null pointer valid return value: indicates that a shape function
      // description for edges is missing.
      // LF_ASSERT_MSG(rfs_edge_p_ != nullptr, "No RSF for edges!");
      return rfs_edge_p_.get();
    }
    case lf::base::RefEl::kTria(): {
      // Null pointer valid return value: indicates that a shape function
      // description for triangular cells is missing.
      // LF_ASSERT_MSG(rfs_tria_p_ != nullptr, "No RSF for triangles!");
      return rfs_tria_p_.get();
    }
    case lf::base::RefEl::kQuad(): {
      // Null pointer valid return value: indicates that a shape function
      // description for quadrilaterals is missing.
      // LF_ASSERT_MSG(rfs_quad_p_ != nullptr, "No RSF for quads!");
      return rfs_quad_p_.get();
    }
    default: {
      LF_VERIFY_MSG(false, "Illegal entity type");
    }
  }
  return nullptr;
}
 
template <typename SCALAR>
size_type UniformScalarFESpace<SCALAR>::NumRefShapeFunctions(
    const lf::mesh::Entity &entity) const {
  return NumRefShapeFunctions(entity.RefEl());
}
 
/* number of _interior_ shape functions associated to entities of various types
 */
template <typename SCALAR>
size_type UniformScalarFESpace<SCALAR>::NumRefShapeFunctions(
    lf::base::RefEl ref_el_type) const {
  LF_ASSERT_MSG((rfs_quad_p_ != nullptr) && (rfs_quad_p_ != nullptr),
                "No valid FE space object: no rsfs");
  // Retrieve number of interior shape functions from rsf layouts
  switch (ref_el_type) {
    case lf::base::RefEl::kPoint(): {
      return num_rsf_node_;
    }
    case lf::base::RefEl::kSegment(): {
      return num_rsf_edge_;
    }
    case lf::base::RefEl::kTria(): {
      return num_rsf_tria_;
    }
    case lf::base::RefEl::kQuad(): {
      return num_rsf_quad_;
    }
    default: {
      LF_VERIFY_MSG(false, "Illegal entity type");
    }
  }
  return 0;
}
 
template <typename SCALAR>
void PrintInfo(std::ostream &o, const UniformScalarFESpace<SCALAR> &fes,
               unsigned ctrl) {
  o << "Uniform scalar FE space, dim = " << fes.LocGlobMap().NumDofs()
    << std::endl;
 
  if ((ctrl & kUniformScalarFESpaceOutMesh) != 0) {
    o << *fes.Mesh() << std::endl;
  }
  if ((ctrl & kUniformScalarFESpaceOutDofh) != 0) {
    o << fes.LocGlobMap() << std::endl;
  }
  if ((ctrl & kUniformScalarFESpaceOutRsfs) != 0) {
    o << fes.NumRefShapeFunctions(lf::base::RefEl::kPoint()) << " rsfs @ nodes"
      << std::endl;
    o << fes.NumRefShapeFunctions(lf::base::RefEl::kSegment())
      << " rsfs @ edges" << std::endl;
    o << fes.NumRefShapeFunctions(lf::base::RefEl::kTria())
      << " rsfs @ triangles" << std::endl;
    o << fes.NumRefShapeFunctions(lf::base::RefEl::kQuad()) << " rsfs @ quads"
      << std::endl;
  }
}
 
template <typename SCALAR>
std::ostream &operator<<(std::ostream &o,
                         const UniformScalarFESpace<SCALAR> &fes) {
  fes.PrintInfo(o, 0);
  return o;
}
 
}  // namespace lf::uscalfe
 
#endif