loc_comp_ellbvp.h

 
#ifndef LF_FE_LOCCOMPELLBVP
#define LF_FE_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 <map>
 
#include "scalar_fe_space.h"
#include "scalar_reference_finite_element.h"
 
namespace lf::fe {
template <typename SCALAR, typename DIFF_COEFF>
class DiffusionElementMatrixProvider final {
  static_assert(mesh::utils::isMeshFunction<DIFF_COEFF>);
 
 public:
  using Scalar =
      typename decltype(mesh::utils::MeshFunctionReturnType<DIFF_COEFF>() *
                        Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>())::Scalar;
  using ElemMat = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
 
  DiffusionElementMatrixProvider(const DiffusionElementMatrixProvider &) =
      delete;
  DiffusionElementMatrixProvider(DiffusionElementMatrixProvider &&) noexcept =
      default;
  DiffusionElementMatrixProvider &operator=(
      const DiffusionElementMatrixProvider &) = delete;
  DiffusionElementMatrixProvider &operator=(DiffusionElementMatrixProvider &&) =
      delete;
  DiffusionElementMatrixProvider(
      std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space, DIFF_COEFF alpha);
 
  bool isActive(const lf::mesh::Entity & /*cell*/) const { return true; }
  ElemMat Eval(const lf::mesh::Entity &cell) const;
 
  ~DiffusionElementMatrixProvider() = default;
 
 private:
  DIFF_COEFF alpha_;
  std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space_;
 
  quad::QuadRuleCache qr_cache_;
};
 
extern std::shared_ptr<spdlog::logger> diffusion_element_matrix_provider_logger;
 
template <class PTR, class DIFF_COEFF>
DiffusionElementMatrixProvider(PTR fe_space, DIFF_COEFF alpha)
    -> DiffusionElementMatrixProvider<typename PTR::element_type::Scalar,
                                      DIFF_COEFF>;
 
// First constructor (internal construction of quadrature rules
template <typename SCALAR, typename DIFF_COEFF>
DiffusionElementMatrixProvider<SCALAR, DIFF_COEFF>::
    DiffusionElementMatrixProvider(
        std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space, DIFF_COEFF alpha)
    : alpha_(std::move(alpha)), fe_space_(std::move(fe_space)) {}
 
// 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 <typename SCALAR, typename DIFF_COEFF>
typename lf::fe::DiffusionElementMatrixProvider<SCALAR, DIFF_COEFF>::ElemMat
DiffusionElementMatrixProvider<SCALAR, DIFF_COEFF>::Eval(
    const lf::mesh::Entity &cell) const {
  // 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(diffusion_element_matrix_provider_logger,
                      "{}, shape = \n{}", cell.RefEl(),
                      geo_ptr->Global(cell.RefEl().NodeCoords()));
  // Physical dimension of the cell
  const dim_t world_dim = geo_ptr->DimGlobal();
 
  // Get a quadrature rule of sufficiently high degree on the element
  const auto sfl = fe_space_->ShapeFunctionLayout(cell);
  const lf::quad::QuadRule qr = qr_cache_.Get(cell.RefEl(), 2 * sfl->Degree());
 
  const Eigen::VectorXd determinants(geo_ptr->IntegrationElement(qr.Points()));
  LF_ASSERT_MSG(
      determinants.size() == qr.NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << qr.NumPoints());
  // Fetch the transformation matrices for the gradients
  const Eigen::MatrixXd JinvT(geo_ptr->JacobianInverseGramian(qr.Points()));
  LF_ASSERT_MSG(JinvT.cols() == 2 * qr.NumPoints(),
                "Mismatch " << JinvT.cols() << " <-> " << 2 * qr.NumPoints());
  LF_ASSERT_MSG(JinvT.rows() == world_dim,
                "Mismatch " << JinvT.rows() << " <-> " << world_dim);
  // Fetch values of diffusion coefficient alpha at quadrature points:
  auto alphaval = alpha_(cell, qr.Points());
  // The requested element matrix is square, size = number of local shape
  // functions
  ElemMat mat(sfl->NumRefShapeFunctions(), sfl->NumRefShapeFunctions());
  mat.setZero();
  // Compute the gradients of the reference shape functions in the quadrature
  // points
  const auto grsf = sfl->GradientsReferenceShapeFunctions(qr.Points());
  // Quadrature formula: loop over quadrature points
  for (base::size_type k = 0; k < qr.NumPoints(); ++k) {
    const double w = qr.Weights()[k] * determinants[k];
    // Transformed gradient in current quadrature node
    const auto trf_grad(JinvT.block(0, 2 * k, world_dim, 2) *
                        grsf.block(0, 2 * k, mat.rows(), 2).transpose());
    // Transformed gradients multiplied with coefficient
    mat += w * trf_grad.adjoint() * (alphaval[k] * trf_grad);
  }
  return mat;
}
 
template <typename SCALAR, typename REACTION_COEFF>
class MassElementMatrixProvider final {
  static_assert(mesh::utils::isMeshFunction<REACTION_COEFF>);
 
 public:
  using Scalar =
      typename decltype(mesh::utils::MeshFunctionReturnType<REACTION_COEFF>() *
                        Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>())::Scalar;
  using ElemMat = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
 
  MassElementMatrixProvider(const MassElementMatrixProvider &) = delete;
  MassElementMatrixProvider(MassElementMatrixProvider &&) noexcept = default;
  MassElementMatrixProvider &operator=(const MassElementMatrixProvider &) =
      delete;
  MassElementMatrixProvider &operator=(MassElementMatrixProvider &&) = delete;
  MassElementMatrixProvider(
      std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space,
      REACTION_COEFF gamma);
 
  bool isActive(const lf::mesh::Entity & /*cell*/) const { return true; }
  ElemMat Eval(const lf::mesh::Entity &cell) const;
 
  ~MassElementMatrixProvider() = default;
 
 private:
  REACTION_COEFF gamma_;
  std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space_;
 
  quad::QuadRuleCache qr_cache_;
};
 
extern std::shared_ptr<spdlog::logger> mass_element_matrix_provider_logger;
 
template <class PTR, class REACTION_COEFF>
MassElementMatrixProvider(PTR fe_space, REACTION_COEFF gamma)
    -> MassElementMatrixProvider<typename PTR::element_type::Scalar,
                                 REACTION_COEFF>;
 
// First constructor
template <typename SCALAR, typename REACTION_COEFF>
MassElementMatrixProvider<SCALAR, REACTION_COEFF>::MassElementMatrixProvider(
    std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space, REACTION_COEFF gamma)
    : gamma_(std::move(gamma)), fe_space_(std::move(fe_space)) {}
 
// 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 <typename SCALAR, typename REACTION_COEFF>
typename lf::fe::MassElementMatrixProvider<SCALAR, REACTION_COEFF>::ElemMat
MassElementMatrixProvider<SCALAR, REACTION_COEFF>::Eval(
    const lf::mesh::Entity &cell) const {
  // 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(mass_element_matrix_provider_logger, "{}, shape = \n{}",
                      cell.RefEl(), geo_ptr->Global(cell.RefEl().NodeCoords()));
  // Physical dimension of the cell
  const dim_t world_dim = geo_ptr->DimGlobal();
 
  // Get a quadrature rule of sufficiently high degree on the element
  const auto sfl = fe_space_->ShapeFunctionLayout(cell);
  const lf::quad::QuadRule qr = qr_cache_.Get(cell.RefEl(), 2 * sfl->Degree());
  // Metric factors in quadrature nodes
  const Eigen::VectorXd determinants(geo_ptr->IntegrationElement(qr.Points()));
  LF_ASSERT_MSG(
      determinants.size() == qr.NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << qr.NumPoints());
  // Fetch values of reaction coefficient gamma at quadrature points:
  auto gammaval = gamma_(cell, qr.Points());
  // Alocated element matrix
  ElemMat mat(sfl->NumRefShapeFunctions(), sfl->NumRefShapeFunctions());
  mat.setZero();
  // Compute the reference shape functions in the quadrature points
  const auto rsf = sfl->EvalReferenceShapeFunctions(qr.Points());
  // Loop over quadrature points to evaluate quadrature formula
  for (base::size_type k = 0; k < qr.NumPoints(); ++k) {
    const double w = qr.Weights()[k] * determinants[k];
    mat += w * ((gammaval[k] * rsf.col(k)) * (rsf.col(k).adjoint()));
  }
  return mat;
}
 
template <typename SCALAR, typename COEFF, typename EDGESELECTOR>
class MassEdgeMatrixProvider final {
 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 ScalarFESpace<SCALAR>> fe_space,
                         COEFF gamma,
                         EDGESELECTOR edge_selector = base::PredicateTrue{})
      : gamma_(std::move(gamma)),
        edge_sel_(std::move(edge_selector)),
        fe_space_(fe_space) {}
 
  bool isActive(const lf::mesh::Entity &edge) const {
    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) const;
 
  ~MassEdgeMatrixProvider() = default;
 
 private:
  COEFF gamma_;            // functor for coefficient
  EDGESELECTOR edge_sel_;  // Defines the active edges
  // FE Space
  std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space_;
 
  quad::QuadRuleCache qr_cache_;
};
 
extern std::shared_ptr<spdlog::logger> mass_edge_matrix_provider_logger;
 
// 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) const {
  const lf::geometry::Geometry *geo_ptr = edge.Geometry();
  // Get the shape function layout on the edge
  const auto sfl = fe_space_->ShapeFunctionLayout(edge);
  // Compute a quadrature rule on the given entity
  const lf::quad::QuadRule qr = qr_cache_.Get(edge.RefEl(), 2 * sfl->Degree());
  // Obtain the metric factors for the quadrature points
  const Eigen::VectorXd determinants(geo_ptr->IntegrationElement(qr.Points()));
  LF_ASSERT_MSG(
      determinants.size() == qr.NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << qr.NumPoints());
 
  // Element matrix
  ElemMat mat(sfl->NumRefShapeFunctions(), sfl->NumRefShapeFunctions());
  mat.setZero();
 
  auto gammaval = gamma_(edge, qr.Points());
 
  // Compute the reference shape functions
  const auto rsf = sfl->EvalReferenceShapeFunctions(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 = (qr.Weights()[k] * determinants[k]) * gammaval[k];
    mat += ((rsf.col(k)) * (rsf.col(k).adjoint())) * w;
  }
  return mat;
}
 
template <typename SCALAR, typename MESH_FUNCTION>
class ScalarLoadElementVectorProvider final {
  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 ScalarFESpace<SCALAR>> fe_space, MESH_FUNCTION f);
  bool isActive(const lf::mesh::Entity & /*cell*/) const { 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) const;
 
  ~ScalarLoadElementVectorProvider() = default;
 
 private:
  MESH_FUNCTION f_;
 
  std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space_;
 
  quad::QuadRuleCache qr_cache_;
};
 
extern std::shared_ptr<spdlog::logger>
    scalar_load_element_vector_provider_logger;
 
// 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 ScalarFESpace<SCALAR>> fe_space, MESH_FUNCTION f)
    : f_(std::move(f)), fe_space_(std::move(fe_space)) {}
 
// 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) const {
  // Type for source function
  using source_fn_t = mesh::utils::MeshFunctionReturnType<MESH_FUNCTION>;
 
  // Get the shape function layout for the given cell
  const auto sfl = fe_space_->ShapeFunctionLayout(cell);
 
  // Initialize a quadrature rule of sufficiently high degree
  const lf::quad::QuadRule qr = qr_cache_.Get(cell.RefEl(), 2 * sfl->Degree());
 
  // 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(scalar_load_element_vector_provider_logger,
                      "{}, shape = \n{}", cell.RefEl(),
                      geo_ptr->Global(cell.RefEl().NodeCoords()));
 
  // Obtain the metric factors for the quadrature points
  const Eigen::VectorXd determinants(geo_ptr->IntegrationElement(qr.Points()));
  LF_ASSERT_MSG(
      determinants.size() == qr.NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << qr.NumPoints());
  SPDLOG_LOGGER_TRACE(scalar_load_element_vector_provider_logger,
                      "LOCVEC({}): Metric factors :\n{}", cell.RefEl(),
                      determinants.transpose());
  // Element vector
  ElemVec vec(sfl->NumRefShapeFunctions());
  vec.setZero();
 
  auto fval = f_(cell, qr.Points());
 
  // Compute the reference shape functions
  const auto rsf = sfl->EvalReferenceShapeFunctions(qr.Points());
 
  // Loop over quadrature points
  for (long k = 0; k < determinants.size(); ++k) {
    SPDLOG_LOGGER_TRACE(scalar_load_element_vector_provider_logger,
                        "LOCVEC: [{}] -> [weight = {}]",
                        qr.Points().transpose(), qr.Weights()[k]);
    // Contribution of current quadrature point
    vec +=
        (qr.Weights()[k] * determinants[k] * fval[k]) * rsf.col(k).conjugate();
  }
 
  SPDLOG_LOGGER_TRACE(scalar_load_element_vector_provider_logger,
                      "LOCVEC = \n{}", vec.transpose());
  return vec;
}
 
template <class SCALAR, class FUNCTOR, class EDGESELECTOR = base::PredicateTrue>
class ScalarLoadEdgeVectorProvider final {
 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 ScalarFESpace<SCALAR>> fe_space, FUNCTOR g,
      EDGESELECTOR edge_sel = base::PredicateTrue{})
      : g_(std::move(g)), edge_sel_(std::move(edge_sel)), fe_space_(fe_space) {}
 
  bool isActive(const lf::mesh::Entity &cell) const { 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) const;
 
  ~ScalarLoadEdgeVectorProvider() = default;
 
 private:
  FUNCTOR g_;              // source function
  EDGESELECTOR edge_sel_;  // selects edges
  std::shared_ptr<const ScalarFESpace<SCALAR>> fe_space_;
 
  quad::QuadRuleCache qr_cache_;
};
 
extern std::shared_ptr<spdlog::logger> scalar_load_edge_vector_provider_logger;
 
// 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) const {
  // Query the shape of the edge
  const lf::geometry::Geometry *geo_ptr = edge.Geometry();
  LF_ASSERT_MSG(geo_ptr != nullptr, "Invalid geometry!");
  LF_ASSERT_MSG(geo_ptr->DimLocal() == 1, "The passed entity is not an edge!");
 
  // Get the shape function layout of the given edge
  const auto sfl = fe_space_->ShapeFunctionLayout(edge);
 
  // Quadrature points on physical edge
  const lf::quad::QuadRule qr = qr_cache_.Get(edge.RefEl(), 2 * sfl->Degree());
 
  // Obtain the metric factors for the quadrature points
  const Eigen::VectorXd determinants(geo_ptr->IntegrationElement(qr.Points()));
  LF_ASSERT_MSG(
      determinants.size() == qr.NumPoints(),
      "Mismatch " << determinants.size() << " <-> " << qr.NumPoints());
 
  // Element vector
  ElemVec vec(sfl->NumRefShapeFunctions());
  vec.setZero();
 
  auto g_vals = g_(edge, qr.Points());
 
  // Compute the reference shape functions
  const auto rsf = sfl->EvalReferenceShapeFunctions(qr.Points());
 
  // Loop over quadrature points
  for (base::size_type k = 0; k < qr.NumPoints(); ++k) {
    // Add contribution of quadrature point to local vector
    const auto w = (qr.Weights()[k] * determinants[k]) * g_vals[k];
    vec += rsf.col(k).conjugate() * w;
  }
  return vec;
}
 
}  // namespace lf::fe
 
#endif