mesh_function_grad_fe.h

 
#ifndef INCGb6997524e2834b5b8e4bba019fb35cc6
#define INCGb6997524e2834b5b8e4bba019fb35cc6
#include "scalar_fe_space.h"
 
namespace lf::fe {
 
template <class SCALAR_FE, class SCALAR_COEFF>
class MeshFunctionGradFE {
 public:
  using Scalar = decltype(SCALAR_FE(0) * SCALAR_COEFF(0));  // NOLINT
 
  MeshFunctionGradFE(std::shared_ptr<const ScalarFESpace<SCALAR_FE>> fe_space,
                     Eigen::Matrix<SCALAR_COEFF, Eigen::Dynamic, 1> dof_vector)
      : fe_space_(std::move(fe_space)), dof_vector_(std::move(dof_vector)) {}
 
  [[nodiscard]] std::shared_ptr<const mesh::Mesh> getMesh() const {
    return fe_space_->Mesh();
  }
 
  [[nodiscard]] std::shared_ptr<const ScalarFESpace<SCALAR_FE>> getFESpace()
      const {
    return fe_space_;
  }
 
  std::vector<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>> operator()(
      const lf::mesh::Entity& e, const Eigen::MatrixXd& local) const {
    // Access information on local shape functions for the entity
    auto grad_sf_eval =
        fe_space_->ShapeFunctionLayout(e)->GradientsReferenceShapeFunctions(
            local);
    // Fetch the coefficients of the global shape functions associated with
    // the current mesh entity
    Eigen::Matrix<SCALAR_COEFF, 1, Eigen::Dynamic> local_dofs(
        1, grad_sf_eval.rows());
    auto global_dofs = fe_space_->LocGlobMap().GlobalDofIndices(e);
    for (Eigen::Index i = 0; i < grad_sf_eval.rows(); ++i) {
      local_dofs(i) = dof_vector_(global_dofs[i]);
    }
    // gradients w.r.t. reference element coordinates \hat{x}
    auto local_grads = (local_dofs * grad_sf_eval).eval();
    // Transform to Cartesian coordinates
    auto jac_t = e.Geometry()->JacobianInverseGramian(local);
    auto dim_local = base::narrow<Eigen::Index>(e.RefEl().Dimension());
    std::vector<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>> result(local.cols());
    // Transform all the local gradients in the evaluation points
    for (Eigen::Index i = 0; i < result.size(); ++i) {
      result[i] = jac_t.block(0, dim_local * i, jac_t.rows(), dim_local) *
                  local_grads.block(0, i * dim_local, 1, dim_local).transpose();
    }
    return result;
  }
 
 private:
  std::shared_ptr<const ScalarFESpace<SCALAR_FE>> fe_space_;
  Eigen::Matrix<SCALAR_COEFF, Eigen::Dynamic, 1> dof_vector_;
};
 
// deduction guide
template <class T, class SCALAR_COEFF>
MeshFunctionGradFE(std::shared_ptr<T>,
                   const Eigen::Matrix<SCALAR_COEFF, Eigen::Dynamic, 1>&)
    -> MeshFunctionGradFE<typename T::Scalar, SCALAR_COEFF>;
 
}  // namespace lf::fe
 
#endif  // INCGb6997524e2834b5b8e4bba019fb35cc6