mesh_function_whitney_one.h

#ifndef HLDO_SPHERE_MESH_FUNCTION_WHITNEY_ONE_H
#define HLDO_SPHERE_MESH_FUNCTION_WHITNEY_ONE_H
/*
 * @file mesh_function_whitney_one.h
 */
 
#include <lf/uscalfe/uscalfe.h>
 
namespace projects::hldo_sphere::post_processing {
 
template <typename SCALAR>
class MeshFunctionWhitneyOne {
 public:
  MeshFunctionWhitneyOne(const Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>& mu,
                         const std::shared_ptr<const lf::mesh::Mesh> mesh)
      : mu_(mu), mesh_(mesh) {
    // make sure mu has the right size
    lf::base::size_type n = mesh->NumEntities(1);
    int mu_size = mu.size();
    LF_VERIFY_MSG(
        n == mu_size,
        "Not the right number of basis expansion coefficiants in mu expected: "
            << n << " given: " << mu_size);
  };
 
  std::vector<Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>> operator()(
      const lf::mesh::Entity& e, const Eigen::MatrixXd& local) const {
    // Only triangles are supported
    LF_VERIFY_MSG(e.RefEl() == lf::base::RefEl::kTria(),
                  "Unsupported cell type " << e.RefEl());
 
    // Get the geometry of the entity
    const auto* geom = e.Geometry();
 
    // Compute the global vertex coordinates
    Eigen::MatrixXd vertices = geom->Global(e.RefEl().NodeCoords());
 
    // get the lambdas
    const lf::uscalfe::FeLagrangeO1Tria<SCALAR> hat_func;
 
    // The gradients are constant on the triangle
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> ref_grads =
        hat_func.GradientsReferenceShapeFunctions(Eigen::VectorXd::Zero(2))
            .transpose();
 
    // The JacobianInverseGramian is constant on the triangle
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> J_inv_trans =
        geom->JacobianInverseGramian(Eigen::VectorXd::Zero(2));
 
    // get the gradients
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> grad =
        J_inv_trans * ref_grads;
 
    // get lambda values
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> lambda =
        hat_func.EvalReferenceShapeFunctions(local);
 
    // correct for orientation
    auto edgeOrientations = e.RelativeOrientations();
 
    // define return vector
    std::vector<Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>> vals(local.cols());
 
    // get global indices of the edges
    std::vector<lf::base::size_type> global_indices(3);
    auto edges = e.SubEntities(1);
    for (lf::base::size_type i = 0; i < 3; i++) {
      global_indices[i] = mesh_->Index(*edges[i]);
    }
 
    // for each point evaluate the basis functions and add the
    // values together
    for (int i = 0; i < local.cols(); i++) {
      Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> bs(grad.rows(), 3);
      bs.col(0) = lambda(0, i) * grad.col(1) - lambda(1, i) * grad.col(0);
      bs.col(1) = lambda(1, i) * grad.col(2) - lambda(2, i) * grad.col(1);
      bs.col(2) = lambda(2, i) * grad.col(0) - lambda(0, i) * grad.col(2);
 
      bs.col(0) *= lf::mesh::to_sign(edgeOrientations[0]);
      bs.col(1) *= lf::mesh::to_sign(edgeOrientations[1]);
      bs.col(2) *= lf::mesh::to_sign(edgeOrientations[2]);
 
      vals[i] = mu_(global_indices[0]) * bs.col(0) +
                mu_(global_indices[1]) * bs.col(1) +
                mu_(global_indices[2]) * bs.col(2);
    }
 
    return vals;
  }
 
 private:
  const Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>& mu_;
  const std::shared_ptr<const lf::mesh::Mesh> mesh_;
};
 
}  // namespace projects::hldo_sphere::post_processing
#endif  // HLDO_SPHERE_MESH_FUNCTION_WHITNEY_ONE_H