loc_comp_ellbvp.h

 
#include <map>
 
#ifndef LF_LOCCOMPELLBVP
#define LF_LOCCOMPELLBVP
/***************************************************************************
 * 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/mesh/utils/utils.h>
#include <lf/quad/quad.h>
 
#include <iostream>
 
#include "precomputed_scalar_reference_finite_element.h"
#include "uscalfe.h"
 
namespace lf::uscalfe {
using quad_rule_collection_t = std::map<lf::base::RefEl, lf::quad::QuadRule>;
 
template <typename SCALAR, typename DIFF_COEFF, typename REACTION_COEFF>
class ReactionDiffusionElementMatrixProvider {
  static_assert(mesh::utils::isMeshFunction<DIFF_COEFF>);
  static_assert(mesh::utils::isMeshFunction<REACTION_COEFF>);
 
 public:
  using Scalar =
      typename decltype(mesh::utils::MeshFunctionReturnType<DIFF_COEFF>() *
                            Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>() +
                        mesh::utils::MeshFunctionReturnType<REACTION_COEFF>() *
                            Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>())::Scalar;
  using ElemMat = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
 
  ReactionDiffusionElementMatrixProvider(
      const ReactionDiffusionElementMatrixProvider &) = delete;
  ReactionDiffusionElementMatrixProvider(
      ReactionDiffusionElementMatrixProvider &&) noexcept = default;
  ReactionDiffusionElementMatrixProvider &operator=(
      const ReactionDiffusionElementMatrixProvider &) = delete;
  ReactionDiffusionElementMatrixProvider &operator=(
      ReactionDiffusionElementMatrixProvider &&) = delete;
  ReactionDiffusionElementMatrixProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
      DIFF_COEFF alpha, REACTION_COEFF gamma);
 
  ReactionDiffusionElementMatrixProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
      DIFF_COEFF alpha, REACTION_COEFF gamma,
      quad_rule_collection_t qr_collection);
 
  virtual bool isActive(const lf::mesh::Entity & /*cell*/) { return true; }
  ElemMat Eval(const lf::mesh::Entity &cell);
 
  virtual ~ReactionDiffusionElementMatrixProvider() = default;
 
 private:
  DIFF_COEFF alpha_;
  REACTION_COEFF gamma_;
  // fe_precomp_[i] contains precomputed reference finite element for ref_el i.
  std::array<PrecomputedScalarReferenceFiniteElement<SCALAR>, 5> fe_precomp_;
};
 
std::shared_ptr<spdlog::logger> &ReactionDiffusionElementMatrixProviderLogger();
 
template <class PTR, class DIFF_COEFF, class REACTION_COEFF>
ReactionDiffusionElementMatrixProvider(PTR fe_space, DIFF_COEFF alpha,
                                       REACTION_COEFF gamma)
    -> ReactionDiffusionElementMatrixProvider<
        typename PTR::element_type::Scalar, DIFF_COEFF, REACTION_COEFF>;
 
template <class PTR, class DIFF_COEFF, class REACTION_COEFF>
ReactionDiffusionElementMatrixProvider(
    PTR fe_space, DIFF_COEFF alpha, REACTION_COEFF gamma,
    std::map<lf::base::RefEl, lf::quad::QuadRule>)
    -> ReactionDiffusionElementMatrixProvider<
        typename PTR::element_type::Scalar, DIFF_COEFF, REACTION_COEFF>;
 
// First constructor (internal construction of quadrature rules
template <typename SCALAR, typename DIFF_COEFF, typename REACTION_COEFF>
ReactionDiffusionElementMatrixProvider<SCALAR, DIFF_COEFF, REACTION_COEFF>::
    ReactionDiffusionElementMatrixProvider(
        std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
        DIFF_COEFF alpha, REACTION_COEFF gamma)
    : alpha_(std::move(alpha)), gamma_(std::move(gamma)), fe_precomp_() {
  for (auto ref_el : {base::RefEl::kTria(), base::RefEl::kQuad()}) {
    auto fe = fe_space->ShapeFunctionLayout(ref_el);
    // Check whether shape functions for that entity type are available.
    // Note that the corresponding PrecomputedScalarReferenceFiniteElement local
    // object is not initialized if the associated description of local shape
    // functions is missing.
    if (fe != nullptr) {
      // Precompute cell-independent quantities based on quadrature rules
      // with twice the degree of exactness compared to the degree of the
      // finite element space.
      fe_precomp_[ref_el.Id()] = PrecomputedScalarReferenceFiniteElement(
          fe, quad::make_QuadRule(ref_el, 2 * fe->Degree()));
    }
  }
}
 
// Second constructor (quadrature rules passed as arguments)
template <typename SCALAR, typename DIFF_COEFF, typename REACTION_COEFF>
ReactionDiffusionElementMatrixProvider<SCALAR, DIFF_COEFF, REACTION_COEFF>::
    ReactionDiffusionElementMatrixProvider(
        std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
        DIFF_COEFF alpha, REACTION_COEFF gamma,
        quad_rule_collection_t qr_collection)
    : alpha_(std::move(alpha)), gamma_(std::move(gamma)), fe_precomp_() {
  for (auto ref_el : {base::RefEl::kTria(), base::RefEl::kQuad()}) {
    // Obtain pointer to an object describing local shape functions
    auto fe = fe_space->ShapeFunctionLayout(ref_el);
    // Check whether shape functions for that entity type are available
    // Note that the corresponding PrecomputedScalarReferenceFiniteElement local
    // object is not initialized if the associated description of local shape
    // functions is missing.
    if (fe != nullptr) {
      // Obtain quadrature rule from user-supplied collection.
      auto qr_coll_ptr = qr_collection.find(ref_el);
      if (qr_coll_ptr != qr_collection.end()) {
        // A quadrature rule for the current entity type is available
        lf::quad::QuadRule qr = qr_coll_ptr->second;
        LF_ASSERT_MSG(qr.RefEl() == ref_el,
                      "qr.RefEl() = " << qr.RefEl() << " <-> " << ref_el);
        // Precomputations of cell-independent quantities
        fe_precomp_[ref_el.Id()] =
            PrecomputedScalarReferenceFiniteElement(fe, qr);
      }
    }
  }
}
 
// Main method for the computation of the element matrix
template <typename SCALAR, typename DIFF_COEFF, typename REACTION_COEFF>
typename lf::uscalfe::ReactionDiffusionElementMatrixProvider<
    SCALAR, DIFF_COEFF, REACTION_COEFF>::ElemMat
ReactionDiffusionElementMatrixProvider<
    SCALAR, DIFF_COEFF, REACTION_COEFF>::Eval(const lf::mesh::Entity &cell) {
  // Topological type of the cell
  const lf::base::RefEl ref_el{cell.RefEl()};
  // Obtain precomputed information about values of local shape functions
  // and their gradients at quadrature points.
  PrecomputedScalarReferenceFiniteElement<SCALAR> &pfe =
      fe_precomp_[ref_el.Id()];
  if (!pfe.isInitialized()) {
    // Accident: cell is of a type not covered by finite element
    // specifications or there is no quadrature rule available for this
    // reference element type
    std::stringstream temp;
    temp << "No local shape function information or no quadrature rule for "
            "reference element type "
         << ref_el;
    throw base::LfException(temp.str());
  }
 
  // Query the shape of the cell
  const lf::geometry::Geometry *geo_ptr = cell.Geometry();
  LF_ASSERT_MSG(geo_ptr != nullptr, "Invalid geometry!");
  LF_ASSERT_MSG((geo_ptr->DimLocal() == 2),
                "Only 2D implementation available!");
  SPDLOG_LOGGER_TRACE(ReactionDiffusionElementMatrixProviderLogger(),
                      "{}, shape = \n{}", ref_el,
                      geo_ptr->Global(ref_el.NodeCoords()));
 
  // Physical dimension of the cell
  const dim_t world_dim = geo_ptr->DimGlobal();
  // Gram determinant at quadrature points
  const Eigen::VectorXd determinants(
      geo_ptr->IntegrationElement(pfe.Qr().Points()));
  LF_ASSERT_MSG(
      determinants.size() == pfe.Qr().NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << pfe.Qr().NumPoints());
  // Fetch the transformation matrices for the gradients
  const Eigen::MatrixXd JinvT(
      geo_ptr->JacobianInverseGramian(pfe.Qr().Points()));
  LF_ASSERT_MSG(
      JinvT.cols() == 2 * pfe.Qr().NumPoints(),
      "Mismatch " << JinvT.cols() << " <-> " << 2 * pfe.Qr().NumPoints());
  LF_ASSERT_MSG(JinvT.rows() == world_dim,
                "Mismatch " << JinvT.rows() << " <-> " << world_dim);
 
  // compute values of coefficients alpha, gamma at quadrature points:
  auto alphaval = alpha_(cell, pfe.Qr().Points());
  auto gammaval = gamma_(cell, pfe.Qr().Points());
 
  // Element matrix
  ElemMat mat(pfe.NumRefShapeFunctions(), pfe.NumRefShapeFunctions());
  mat.setZero();
 
  // Loop over quadrature points
  for (base::size_type k = 0; k < pfe.Qr().NumPoints(); ++k) {
    const double w = pfe.Qr().Weights()[k] * determinants[k];
    // Transformed gradients
    const auto trf_grad(JinvT.block(0, 2 * k, world_dim, 2) *
                        pfe.PrecompGradientsReferenceShapeFunctions()
                            .block(0, 2 * k, mat.rows(), 2)
                            .transpose());
    // Transformed gradients multiplied with coefficient
    const auto alpha_trf_grad(alphaval[k] * trf_grad);
    mat += w * (trf_grad.adjoint() * alpha_trf_grad +
                (gammaval[k] * pfe.PrecompReferenceShapeFunctions().col(k)) *
                    (pfe.PrecompReferenceShapeFunctions().col(k).adjoint()));
  }
  return mat;
}
 
template <typename SCALAR, typename COEFF, typename EDGESELECTOR>
class MassEdgeMatrixProvider {
 public:
  using scalar_t =
      decltype(SCALAR(0) * mesh::utils::MeshFunctionReturnType<COEFF>(0));
  using ElemMat = Eigen::Matrix<scalar_t, Eigen::Dynamic, Eigen::Dynamic>;
 
  MassEdgeMatrixProvider(const MassEdgeMatrixProvider &) = delete;
  MassEdgeMatrixProvider(MassEdgeMatrixProvider &&) noexcept = default;
  MassEdgeMatrixProvider &operator=(const MassEdgeMatrixProvider &) = delete;
  MassEdgeMatrixProvider &operator=(MassEdgeMatrixProvider &&) = delete;
  MassEdgeMatrixProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space, COEFF gamma,
      EDGESELECTOR edge_selector = base::PredicateTrue{})
      : gamma_(std::move(gamma)),
        edge_sel_(std::move(edge_selector)),
        fe_precomp_() {
    auto fe = fe_space->ShapeFunctionLayout(base::RefEl::kSegment());
    LF_ASSERT_MSG(fe != nullptr, "No shape functions specified for edges");
    // Precompute entity-independent quantities based on a LehrFEM++ built-in
    // quadrature rule
    fe_precomp_ = PrecomputedScalarReferenceFiniteElement(
        fe, quad::make_QuadRule(base::RefEl::kSegment(), 2 * fe->Degree()));
  }
  MassEdgeMatrixProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space, COEFF gamma,
      lf::quad::QuadRule quadrule,
      EDGESELECTOR edge_selector = base::PredicateTrue{})
      : gamma_(std::move(gamma)),
        edge_sel_(std::move(edge_selector)),
        fe_precomp_() {
    auto fe = fe_space->ShapeFunctionLayout(base::RefEl::kSegment());
    LF_ASSERT_MSG(fe != nullptr, "No shape functions specified for edges");
    LF_ASSERT_MSG(quadrule.RefEl() == base::RefEl::kSegment(),
                  "Quadrature rule not meant for EDGE entities!");
    // Precompute entity-independent quantities
    fe_precomp_ = PrecomputedScalarReferenceFiniteElement(fe, quadrule);
  }
 
  bool isActive(const lf::mesh::Entity &edge) {
    LF_ASSERT_MSG(edge.RefEl() == lf::base::RefEl::kSegment(),
                  "Wrong type for an edge");
    return edge_sel_(edge);
  }
 
  ElemMat Eval(const lf::mesh::Entity &edge);
 
  virtual ~MassEdgeMatrixProvider() = default;
 
 private:
  COEFF gamma_;            // functor for coefficient
  EDGESELECTOR edge_sel_;  // Defines the active edges
  // Precomputed quantities at quadrature points
  PrecomputedScalarReferenceFiniteElement<SCALAR> fe_precomp_;
};
 
std::shared_ptr<spdlog::logger> &MassEdgeMatrixProviderLogger();
 
// deduction guide:
template <class PTR, class COEFF, class EDGESELECTOR = base::PredicateTrue>
MassEdgeMatrixProvider(PTR, COEFF coeff,
                       EDGESELECTOR edge_predicate = base::PredicateTrue{})
    -> MassEdgeMatrixProvider<typename PTR::element_type::Scalar, COEFF,
                              EDGESELECTOR>;
 
// Eval() method
// TODO(craffael) remove const once
// https://
// developercommunity.visualstudio.com/content/problem/180948/vs2017-155-c-cv-qualifiers-lost-on-type-alias-used.html
// is resolved
template <class SCALAR, class COEFF, class EDGESELECTOR>
typename MassEdgeMatrixProvider<SCALAR, COEFF, EDGESELECTOR>::ElemMat
MassEdgeMatrixProvider<SCALAR, COEFF, EDGESELECTOR>::Eval(
    const lf::mesh::Entity &edge) {
  // Topological type of the cell
  const lf::base::RefEl ref_el{edge.RefEl()};
  LF_ASSERT_MSG(ref_el == lf::base::RefEl::kSegment(),
                "Edge must be of segment type");
  // Query the shape of the edge
  const lf::geometry::Geometry *geo_ptr = edge.Geometry();
  LF_ASSERT_MSG(geo_ptr != nullptr, "Invalid geometry!");
 
  // Obtain the metric factors for the quadrature points
  const Eigen::VectorXd determinants(
      geo_ptr->IntegrationElement(fe_precomp_.Qr().Points()));
  LF_ASSERT_MSG(determinants.size() == fe_precomp_.Qr().NumPoints(),
                "Mismatch " << determinants.size() << " <-> "
                            << fe_precomp_.Qr().NumPoints());
 
  // Element matrix
  ElemMat mat(fe_precomp_.NumRefShapeFunctions(),
              fe_precomp_.NumRefShapeFunctions());
  mat.setZero();
 
  auto gammaval = gamma_(edge, fe_precomp_.Qr().Points());
 
  // Loop over quadrature points
  for (long k = 0; k < determinants.size(); ++k) {
    // Build local matrix by summing rank-1 contributions
    // from quadrature points.
    const auto w =
        (fe_precomp_.Qr().Weights()[k] * determinants[k]) * gammaval[k];
    mat += ((fe_precomp_.PrecompReferenceShapeFunctions().col(k)) *
            (fe_precomp_.PrecompReferenceShapeFunctions().col(k).adjoint())) *
           w;
  }
  return mat;
}
 
template <typename SCALAR, typename MESH_FUNCTION>
class ScalarLoadElementVectorProvider {
  static_assert(mesh::utils::isMeshFunction<MESH_FUNCTION>);
 
 public:
  using scalar_t = decltype(
      SCALAR(0) * mesh::utils::MeshFunctionReturnType<MESH_FUNCTION>(0));
  using ElemVec = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
 
  ScalarLoadElementVectorProvider(const ScalarLoadElementVectorProvider &) =
      delete;
  ScalarLoadElementVectorProvider(ScalarLoadElementVectorProvider &&) noexcept =
      default;
  ScalarLoadElementVectorProvider &operator=(
      const ScalarLoadElementVectorProvider &) = delete;
  ScalarLoadElementVectorProvider &operator=(
      ScalarLoadElementVectorProvider &&) = delete;
  ScalarLoadElementVectorProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
      MESH_FUNCTION f);
  ScalarLoadElementVectorProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
      MESH_FUNCTION f, quad_rule_collection_t qr_collection);
  virtual bool isActive(const lf::mesh::Entity & /*cell*/) { return true; }
  /*
   * @brief Main method for computing the element vector
   *
   * @param cell current cell for which the element vector is desired
   * @return local load vector as column vector
   *
   */
  ElemVec Eval(const lf::mesh::Entity &cell);
 
  virtual ~ScalarLoadElementVectorProvider() = default;
 
 private:
  MESH_FUNCTION f_;
 
  std::array<PrecomputedScalarReferenceFiniteElement<SCALAR>, 5> fe_precomp_;
};
 
std::shared_ptr<spdlog::logger> &ScalarLoadElementVectorProviderLogger();
 
// Deduction guide
template <class PTR, class MESH_FUNCTION>
ScalarLoadElementVectorProvider(PTR fe_space, MESH_FUNCTION mf)
    -> ScalarLoadElementVectorProvider<typename PTR::element_type::Scalar,
                                       MESH_FUNCTION>;
 
// Constructors
template <typename SCALAR, typename MESH_FUNCTION>
ScalarLoadElementVectorProvider<SCALAR, MESH_FUNCTION>::
    ScalarLoadElementVectorProvider(
        std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
        MESH_FUNCTION f)
    : f_(std::move(f)) {
  for (auto ref_el : {base::RefEl::kTria(), base::RefEl::kQuad()}) {
    auto fe = fe_space->ShapeFunctionLayout(ref_el);
    // Check whether shape functions for that entity type are available
    if (fe != nullptr) {
      // Precompute cell-independent quantities based on quadrature rules
      // with twice the degree of exactness compared to the degree of the
      // finite element space.
      fe_precomp_[ref_el.Id()] =
          PrecomputedScalarReferenceFiniteElement<SCALAR>(
              fe, quad::make_QuadRule(ref_el, 2 * fe->Degree()));
    }
  }
}
 
template <typename SCALAR, typename MESH_FUNCTION>
ScalarLoadElementVectorProvider<SCALAR, MESH_FUNCTION>::
    ScalarLoadElementVectorProvider(
        std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space,
        MESH_FUNCTION f, quad_rule_collection_t qr_collection)
    : f_(std::move(f)) {
  for (auto ref_el : {base::RefEl::kTria(), base::RefEl::kQuad()}) {
    auto fe = fe_space->ShapeFunctionLayout(ref_el);
    // Check whether shape functions for that entity type are available
    if (fe != nullptr) {
      // Obtain quadrature rule from user-supplied collection.
      auto qr_coll_ptr = qr_collection.find(ref_el);
      if (qr_coll_ptr != qr_collection.end()) {
        // A quadrature rule for the current entity type is available
        lf::quad::QuadRule qr = qr_coll_ptr->second;
        LF_ASSERT_MSG(qr.RefEl() == ref_el,
                      "qr.RefEl() = " << qr.RefEl() << " <-> " << ref_el);
        // Precompute cell-independent quantities using the user-supplied
        // quadrature rules
        fe_precomp_[ref_el.Id()] =
            PrecomputedScalarReferenceFiniteElement<SCALAR>(fe, qr);
      } else {
        // Quadrature rule is missing for an entity type for which
        // local shape functions are available
        LF_ASSERT_MSG(false, "Quadrature rule missing for " << ref_el);
      }
    }
  }
}
 
// TODO(craffael) remove const once
// http://developercommunity.visualstudio.com/content/problem/180948/vs2017-155-c-cv-qualifiers-lost-on-type-alias-used.html
// is resolved
template <typename SCALAR, typename MESH_FUNCTION>
typename ScalarLoadElementVectorProvider<SCALAR, MESH_FUNCTION>::ElemVec
ScalarLoadElementVectorProvider<SCALAR, MESH_FUNCTION>::Eval(
    const lf::mesh::Entity &cell) {
  // Type for source function
  using source_fn_t = mesh::utils::MeshFunctionReturnType<MESH_FUNCTION>;
  // Topological type of the cell
  const lf::base::RefEl ref_el{cell.RefEl()};
  // Obtain precomputed information about values of local shape functions
  // and their gradients at quadrature points.
  auto &pfe = fe_precomp_[ref_el.Id()];
  // Accident: cell is of a type not coverence by finite element specifications
  LF_ASSERT_MSG(
      pfe.isInitialized(),
      "No local shape function information for entity type " << ref_el);
 
  // Query the shape of the cell
  const lf::geometry::Geometry *geo_ptr = cell.Geometry();
  LF_ASSERT_MSG(geo_ptr != nullptr, "Invalid geometry!");
  LF_ASSERT_MSG((geo_ptr->DimLocal() == 2),
                "Only 2D implementation available!");
  SPDLOG_LOGGER_TRACE(ScalarLoadElementVectorProviderLogger(),
                      "{}, shape = \n{}", ref_el,
                      geo_ptr->Global(ref_el.NodeCoords()));
 
  // Obtain the metric factors for the quadrature points
  const Eigen::VectorXd determinants(
      geo_ptr->IntegrationElement(pfe.Qr().Points()));
  LF_ASSERT_MSG(
      determinants.size() == pfe.Qr().NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << pfe.Qr().NumPoints());
  SPDLOG_LOGGER_TRACE(ScalarLoadElementVectorProviderLogger(),
                      "LOCVEC({}): Metric factors :\n{}", ref_el,
                      determinants.transpose());
 
  // Element vector
  ElemVec vec(pfe.NumRefShapeFunctions());
  vec.setZero();
 
  auto fval = f_(cell, pfe.Qr().Points());
 
  // Loop over quadrature points
  for (long k = 0; k < determinants.size(); ++k) {
    SPDLOG_LOGGER_TRACE(ScalarLoadElementVectorProviderLogger(),
                        "LOCVEC: [{}] -> [weight = {}]",
                        pfe.Qr().Points().transpose(), pfe.Qr().Weights()[k]);
    // Contribution of current quadrature point
    vec += (pfe.Qr().Weights()[k] * determinants[k] * fval[k]) *
           pfe.PrecompReferenceShapeFunctions().col(k).conjugate();
  }
 
  SPDLOG_LOGGER_TRACE(ScalarLoadElementVectorProviderLogger(), "LOCVEC = \n{}",
                      vec.transpose());
 
  return vec;
}
 
template <class SCALAR, class FUNCTOR, class EDGESELECTOR = base::PredicateTrue>
class ScalarLoadEdgeVectorProvider {
 public:
  static_assert(mesh::utils::isMeshFunction<FUNCTOR>,
                "FUNCTOR does not fulfill the concept of a mesh function.");
  using Scalar =
      decltype(SCALAR(0) * mesh::utils::MeshFunctionReturnType<FUNCTOR>(0));
  using ElemVec = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
 
  ScalarLoadEdgeVectorProvider(const ScalarLoadEdgeVectorProvider &) = delete;
  ScalarLoadEdgeVectorProvider(ScalarLoadEdgeVectorProvider &&) noexcept =
      default;
  ScalarLoadEdgeVectorProvider &operator=(
      const ScalarLoadEdgeVectorProvider &) = delete;
  ScalarLoadEdgeVectorProvider &operator=(ScalarLoadEdgeVectorProvider &&) =
      delete;
  ScalarLoadEdgeVectorProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space, FUNCTOR g,
      EDGESELECTOR edge_sel = base::PredicateTrue{})
      : g_(std::move(g)), edge_sel_(std::move(edge_sel)), pfe_() {
    auto fe = fe_space->ShapeFunctionLayout(base::RefEl::kSegment());
    LF_ASSERT_MSG(fe != nullptr, "No shape functions specified for edges");
    // Precompute entity-independent quantities based on a LehrFEM++ built-in
    // quadrature rule
    pfe_ = PrecomputedScalarReferenceFiniteElement(
        fe, quad::make_QuadRule(base::RefEl::kSegment(), 2 * fe->Degree()));
  }
 
  ScalarLoadEdgeVectorProvider(
      std::shared_ptr<const UniformScalarFESpace<SCALAR>> fe_space, FUNCTOR g,
      lf::quad::QuadRule quadrule,
      EDGESELECTOR edge_sel = base::PredicateTrue{})
      : g_(std::move(g)), edge_sel_(std::move(edge_sel)), pfe_() {
    auto fe = fe_space->ShapeFunctionLayout(base::RefEl::kSegment());
    LF_ASSERT_MSG(fe != nullptr, "No shape functions specified for edges");
    LF_ASSERT_MSG(quadrule.RefEl() == base::RefEl::kSegment(),
                  "Quadrature rule not meant for EDGE entities!");
    // Precompute entity-independent quantities
    pfe_ = PrecomputedScalarReferenceFiniteElement(fe, quadrule);
  }
 
  virtual bool isActive(const lf::mesh::Entity &cell) {
    return edge_sel_(cell);
  }
  /*
   * @brief Main method for computing the element vector
   *
   * @param cell current cell for which the element vector is desired
   * @return local load vector as column vector
   *
   */
  ElemVec Eval(const lf::mesh::Entity &edge);
 
  virtual ~ScalarLoadEdgeVectorProvider() = default;
 
 private:
  FUNCTOR g_;              // source function
  EDGESELECTOR edge_sel_;  // selects edges
  PrecomputedScalarReferenceFiniteElement<SCALAR> pfe_;
};
 
std::shared_ptr<spdlog::logger> &ScalarLoadEdgeVectorProviderLogger();
 
// deduction guide
template <class PTR, class FUNCTOR, class EDGESELECTOR = base::PredicateTrue>
ScalarLoadEdgeVectorProvider(PTR, FUNCTOR, EDGESELECTOR = base::PredicateTrue{})
    -> ScalarLoadEdgeVectorProvider<typename PTR::element_type::Scalar, FUNCTOR,
                                    EDGESELECTOR>;
 
// Eval() method
// TODO(craffael) remove const once
// https://developercommunity.visualstudio.com/content/problem/180948/vs2017-155-c-cv-qualifiers-lost-on-type-alias-used.html
// is resolved
template <class SCALAR, class FUNCTOR, class EDGESELECTOR>
typename ScalarLoadEdgeVectorProvider<SCALAR, FUNCTOR, EDGESELECTOR>::ElemVec
ScalarLoadEdgeVectorProvider<SCALAR, FUNCTOR, EDGESELECTOR>::Eval(
    const lf::mesh::Entity &edge) {
  // Topological type of the cell
  const lf::base::RefEl ref_el{edge.RefEl()};
  LF_ASSERT_MSG(ref_el == lf::base::RefEl::kSegment(),
                "Edge must be of segment type");
  // Query the shape of the edge
  const lf::geometry::Geometry *geo_ptr = edge.Geometry();
  LF_ASSERT_MSG(geo_ptr != nullptr, "Invalid geometry!");
 
  // Quadrature points on physical edge
  const Eigen::MatrixXd mapped_qpts(geo_ptr->Global(pfe_.Qr().Points()));
  LF_ASSERT_MSG(
      mapped_qpts.cols() == pfe_.Qr().NumPoints(),
      "Mismatch " << mapped_qpts.cols() << " <-> " << pfe_.Qr().NumPoints());
 
  // Obtain the metric factors for the quadrature points
  const Eigen::VectorXd determinants(
      geo_ptr->IntegrationElement(pfe_.Qr().Points()));
  LF_ASSERT_MSG(
      determinants.size() == pfe_.Qr().NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << pfe_.Qr().NumPoints());
 
  // Element vector
  ElemVec vec(pfe_.NumRefShapeFunctions());
  vec.setZero();
 
  auto g_vals = g_(edge, pfe_.Qr().Points());
 
  // Loop over quadrature points
  for (base::size_type k = 0; k < pfe_.Qr().NumPoints(); ++k) {
    // Add contribution of quadrature point to local vector
    const auto w = (pfe_.Qr().Weights()[k] * determinants[k]) * g_vals[k];
    vec += pfe_.PrecompReferenceShapeFunctions().col(k).conjugate() * w;
  }
  return vec;
}
 
}  // namespace lf::uscalfe
 
#endif