whitney_two_vector_provider.h

#ifndef HLDO_SPHERE_ASSEMBLE_WHITNEY_TWO_VECTOR_PROVIDER_H
#define HLDO_SPHERE_ASSEMBLE_WHITNEY_TWO_VECTOR_PROVIDER_H
 
#include <lf/mesh/entity.h>
#include <lf/mesh/utils/mesh_data_set.h>
#include <lf/quad/quad.h>
#include <lf/uscalfe/lagr_fe.h>
 
#include <Eigen/Dense>
#include <functional>
#include <utility>
 
namespace projects::hldo_sphere {
 
namespace assemble {
 
template <typename SCALAR>
class WhitneyTwoVectorProvider {
 public:
  WhitneyTwoVectorProvider(
      const std::function<SCALAR(const Eigen::Vector3d &)> &f)
      : f_(std::move(f)) {}
 
  Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> Eval(
      const lf::mesh::Entity &entity) const {
    const auto *const geom = entity.Geometry();
 
    // Only triangles are supported
    LF_VERIFY_MSG(entity.RefEl() == lf::base::RefEl::kTria(),
                  "Unsupported cell type " << entity.RefEl());
 
    // Compute the global vertex coordinates
    Eigen::MatrixXd vertices = geom->Global(entity.RefEl().NodeCoords());
 
    // define quad rule with sufficiantly high degree since the
    // baricentric coordinate functions have degree 1
    lf::quad::QuadRule quadrule{lf::quad::make_TriaQR_EdgeMidpointRule()};
 
    // Compute the elements of the vector with given quadrature rule
    const Eigen::MatrixXd points = geom->Global(quadrule.Points());
    const Eigen::VectorXd weights =
        (geom->IntegrationElement(quadrule.Points()).array() *
         quadrule.Weights().array())
            .matrix();
 
    SCALAR sum = 0;
    for (lf::base::size_type n = 0; n < quadrule.NumPoints(); ++n) {
      Eigen::VectorXd x = points.col(n);
      sum += weights[n] * f_(x);
    }
 
    Eigen::Matrix<SCALAR, 1, 1> element_vector;
    element_vector(0) = sum;
 
    return element_vector;
  }
 
  bool isActive(const lf::mesh::Entity &entity) const { return true; }
 
 private:
  const std::function<SCALAR(const Eigen::Vector3d &)> f_;
};
 
}  // end namespace assemble
 
}  // namespace projects::hldo_sphere
 
#endif  // HLDO_SPHERE_ASSEMBLE_WHITNEY_TWO_VECTOR_PROVIDER_H