whitney_one_vector_provider.h

#ifndef HLDO_SPHERE_ASSEMBLE_WHITNEY_ONE_VECTOR_PROVIDER_H
#define HLDO_SPHERE_ASSEMBLE_WHITNEY_ONE_VECTOR_PROVIDER_H
 
#include <lf/base/base.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 WhitneyOneVectorProvider {
 public:
  WhitneyOneVectorProvider(
      const std::function<Eigen::Matrix<SCALAR, 3, 1>(const Eigen::Vector3d &)>
          &f)
      : f_(std::move(f)) {}
 
  Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> Eval(
      const lf::mesh::Entity &entity) const {
    LF_ASSERT_MSG(entity.RefEl() == lf::base::RefEl::kTria(),
                  "Only Triangles are supported no " << entity.RefEl());
 
    const auto *const geom = entity.Geometry();
 
    // Compute the global vertex coordinates
    Eigen::MatrixXd vertices = geom->Global(entity.RefEl().NodeCoords());
 
    // define quad rule with sufficiantly high degree
    lf::quad::QuadRule quadrule{lf::quad::make_TriaQR_EdgeMidpointRule()};
 
    // Construct the whitney one form basis functions
    const lf::uscalfe::FeLagrangeO1Tria<SCALAR> hat_func;
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> ref_grads =
        hat_func.GradientsReferenceShapeFunctions(Eigen::VectorXd::Zero(2))
            .transpose();
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> J_inv_trans =
        geom->JacobianInverseGramian(Eigen::VectorXd::Zero(2));
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> grads =
        J_inv_trans * ref_grads;
 
    // get edge orientations
    auto edgeOrientations = entity.RelativeOrientations();
    Eigen::Vector3d s;
    s(0) = lf::mesh::to_sign(edgeOrientations[0]);
    s(1) = lf::mesh::to_sign(edgeOrientations[1]);
    s(2) = lf::mesh::to_sign(edgeOrientations[2]);
 
    const auto b_hat = [&](Eigen::Vector2d x_hat)
        -> Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> {
      Eigen::Matrix<SCALAR, 3, 3> b_hat;
      Eigen::Matrix<SCALAR, 3, 1> lambda_hat =
          hat_func.EvalReferenceShapeFunctions(x_hat);
      for (int i = 0; i < 3; i++) {
        int ip1 = (i + 1) % 3;
 
        b_hat.col(i) = s(i) * (lambda_hat(i) * grads.col(ip1) -
                               lambda_hat(ip1) * grads.col(i));
      }
 
      return b_hat;
    };
 
    // returns evaluation of @f$\textbf{f} * \textbf{b}@f$ for at a given point
    // on the reference triangle f is evaluated at the global coordinates of the
    // triangle radially projected on the sphere
    const auto f_tilde_hat =
        [&](Eigen::Vector2d x_hat) -> Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> {
      Eigen::Vector3d x = geom->Global(x_hat);
      Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> result =
          b_hat(x_hat).transpose() * f_(x);
      return result;
    };
 
    // Compute the elements of the vector with given quadrature rule
    Eigen::Matrix<SCALAR, 3, 1> element_vector;
    element_vector.setZero();
    const Eigen::MatrixXd points = quadrule.Points();
    const Eigen::VectorXd weights =
        (geom->IntegrationElement(points).array() * quadrule.Weights().array())
            .matrix();
    for (lf::base::size_type n = 0; n < quadrule.NumPoints(); ++n) {
      const Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> f_eval =
          f_tilde_hat(points.col(n));
      element_vector += weights[n] * f_eval;
    }
    return element_vector;
  }
 
  bool isActive(const lf::mesh::Entity &entity) const { return true; }
 
 private:
  const std::function<Eigen::Matrix<SCALAR, 3, 1>(const Eigen::Vector3d &)> f_;
};
 
}  // end namespace assemble
 
}  // namespace projects::hldo_sphere
 
#endif  // HLDO_SPHERE_ASSEMBLE_WHITNEY_ONE_VECTOR_PROVIDER_H