mesh_function_whitney_two.h

#ifndef HLDO_SPHERE_MESH_FUNCTION_WHITNEY_TWO_H
#define HLDO_SPHERE_MESH_FUNCTION_WHITNEY_TWO_H
/*
 * @file mesh_function_whitney_two.h
 */
 
#include <lf/uscalfe/uscalfe.h>
 
namespace projects::hldo_sphere::post_processing {
 
template <typename SCALAR>
class MeshFunctionWhitneyTwo {
 public:
  MeshFunctionWhitneyTwo(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 same size as number of cells
    lf::base::size_type n = mesh->NumEntities(0);
    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<SCALAR> 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());
 
    // accoring to the lehrfem documentation the mesh index corresponds to the
    // dofhandler if all basis function are associated with entities off the
    // same codomain
    lf::base::size_type global_index = mesh_->Index(e);
 
    std::vector<SCALAR> vals(local.cols());
    for (int i = 0; i < local.cols(); i++) {
      vals[i] = mu_(global_index);
    }
 
    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_TWO_H