whitney_zero_hodge_laplacian.h

#ifndef HLDO_SPHERE_OPERATORS_WHITNEY_ZERO_HODGE_LAPLACE_H
#define HLDO_SPHERE_OPERATORS_WHITNEY_ZERO_HODEE_LAPLACE_H
 
#include <laplace_matrix_provider.h>
#include <lf/assemble/assembler.h>
#include <lf/assemble/coomatrix.h>
#include <lf/assemble/dofhandler.h>
#include <lf/mesh/entity.h>
#include <lf/mesh/hybrid2d/mesh.h>
#include <lf/mesh/hybrid2d/mesh_factory.h>
#include <lf/mesh/mesh_interface.h>
#include <load_vector_provider.h>
#include <sphere_triag_mesh_builder.h>
 
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <iostream>
 
namespace projects::hldo_sphere {
 
namespace operators {
 
template <typename SCALAR>
class WhitneyZeroHodgeLaplace {
 public:
  WhitneyZeroHodgeLaplace()
      : coo_matrix_(lf::assemble::COOMatrix<SCALAR>(1, 1)) {
    // create mesh factory for basic mesh
    std::unique_ptr<lf::mesh::MeshFactory> factory =
        std::make_unique<lf::mesh::hybrid2d::MeshFactory>(3);
 
    std::shared_ptr<projects::hldo_sphere::mesh::SphereTriagMeshBuilder>
        sphere = std::make_shared<
            projects::hldo_sphere::mesh::SphereTriagMeshBuilder>(
            std::move(factory));
    sphere->setRefinementLevel(0);
    sphere->setRadius(1);
 
    mesh_p_ = sphere->Build();
 
    // create basic function
    auto f = [](Eigen::Matrix<double, 3, 1> x) -> SCALAR { return 0; };
    f_ = f;
  }
 
  void Compute() {
    // create element matrix provider
    projects::hldo_sphere::assemble::LaplaceMatrixProvider matrix_provider;
 
    const lf::assemble::DofHandler& dof_handler =
        lf::assemble::UniformFEDofHandler(mesh_p_,
                                          {{lf::base::RefEl::kPoint(), 1}});
 
    // create COO Matrix
    const lf::assemble::size_type n_dofs(dof_handler.NumDofs());
 
    lf::assemble::COOMatrix<SCALAR> coo_mat(n_dofs, n_dofs);
    coo_mat.setZero();
 
    lf::assemble::AssembleMatrixLocally<lf::assemble::COOMatrix<SCALAR>>(
        0, dof_handler, dof_handler, matrix_provider, coo_mat);
 
    // create element vector provider
    projects::hldo_sphere::assemble::LoadVectorProvider<SCALAR> vector_provider{
        f_};
 
    // create load vector
    Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> phi(n_dofs);
    phi.setZero();
 
    // assemble the global vector over entities with codim=0:
    AssembleVectorLocally(0, dof_handler, vector_provider, phi);
 
    // set the class attributes
    phi_ = phi;
    coo_matrix_ = coo_mat;
  }
 
  void SetMesh(std::shared_ptr<const lf::mesh::Mesh> mesh_p) {
    // check if cells are triagles
    for (const lf::mesh::Entity* tria : mesh_p->Entities(0)) {
      LF_ASSERT_MSG(
          tria->RefEl() == lf::base::RefEl::kTria(),
          "Mesh must be Triangular, unsupported cell " << tria->RefEl());
    }
 
    // check if dimension of the mesh is 3
    LF_ASSERT_MSG(mesh_p->DimWorld() == 3,
                  "World Dimension must be 3 but is" << mesh_p->DimWorld());
 
    // set mesh
    mesh_p_ = mesh_p;
  }
 
  void SetLoadFunction(
      std::function<SCALAR(const Eigen::Matrix<double, 3, 1>&)>& f) {
    f_ = f;
  }
 
  Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> GetLoadVector() { return phi_; }
 
  lf::assemble::COOMatrix<SCALAR> GetGalerkinMatrix() { return coo_matrix_; }
 
 private:
  std::shared_ptr<const lf::mesh::Mesh> mesh_p_;
  std::function<SCALAR(const Eigen::Matrix<double, 3, 1>&)> f_;
  lf::assemble::COOMatrix<SCALAR> coo_matrix_;
  Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> phi_;
};
 
}  // namespace operators
 
}  // namespace projects::hldo_sphere
 
#endif  // HLDO_SPHERE_OPERATORS_WHITNEY_ZERO_HODGE_LAPLACE_H