dirac_operator_experiment.h

#ifndef HLDO_SPHERE_EXPERIMENTS_DIRAC_OPERATOR_EXAMPLE_H
#define HLDO_SPHERE_EXPERIMENTS_DIRAC_OPERATOR_EXAMPLE_H
 
#include <dirac_operator_source_problem.h>
#include <lf/io/vtk_writer.h>
#include <lf/mesh/hybrid2d/mesh_factory.h>
#include <lf/mesh/utils/tp_triag_mesh_builder.h>
#include <lf/mesh/utils/utils.h>
#include <lf/quad/quad.h>
#include <lf/refinement/refinement.h>
#include <lf/uscalfe/uscalfe.h>
#include <mesh_function_whitney_one.h>
#include <mesh_function_whitney_two.h>
#include <mesh_function_whitney_zero.h>
#include <norms.h>
#include <results_processing.h>
#include <sphere_triag_mesh_builder.h>
 
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
 
namespace projects::hldo_sphere::experiments {
 
using complex = std::complex<double>;
 
class DiracOperatorExperiment {
 public:
  DiracOperatorExperiment(
      std::function<complex(const Eigen::Matrix<double, 3, 1> &)> u_zero,
      std::function<
          Eigen::Matrix<complex, 3, 1>(const Eigen::Matrix<double, 3, 1> &)>
          u_one,
      std::function<complex(const Eigen::Matrix<double, 3, 1> &)> u_two,
      std::function<complex(const Eigen::Matrix<double, 3, 1> &)> f_zero,
      std::function<
          Eigen::Matrix<complex, 3, 1>(const Eigen::Matrix<double, 3, 1> &)>
          f_one,
      std::function<complex(const Eigen::Matrix<double, 3, 1> &)> f_two,
      double &k, std::string name)
      : u_zero_(u_zero),
        u_one_(u_one),
        u_two_(u_two),
        f_zero_(f_zero),
        f_one_(f_one),
        f_two_(f_two),
        k_(k),
        name_(name) {}
 
  void Compute(std::vector<unsigned> refinement_levels,
               std::vector<double> ks) {
    int nl = refinement_levels.size();
    int nk = ks.size();
 
    // Initialize solution wrapper
    projects::hldo_sphere::post_processing::ProblemSolutionWrapper<complex>
        solutions;
    solutions.k = ks;
    solutions.levels = refinement_levels;
    solutions.mesh = std::vector<std::shared_ptr<const lf::mesh::Mesh>>(nl);
    solutions.solutions = std::vector<
        projects::hldo_sphere::post_processing::ProblemSolution<complex>>(nl);
 
    projects::hldo_sphere::operators::DiracOperatorSourceProblem lse_builder;
 
    // functions are passed by reference hence changing the k still influences
    // the functions
    lse_builder.SetLoadFunctions(f_zero_, f_one_, f_two_);
 
    // Read the mesh from the gmsh file
    std::unique_ptr<lf::mesh::MeshFactory> factory =
        std::make_unique<lf::mesh::hybrid2d::MeshFactory>(3);
 
    projects::hldo_sphere::mesh::SphereTriagMeshBuilder sphere =
        projects::hldo_sphere::mesh::SphereTriagMeshBuilder(std::move(factory));
 
    // we always take the same radius
    sphere.setRadius(1.);
 
    // start timer for total time
    std::chrono::steady_clock::time_point start_time_total =
        std::chrono::steady_clock::now();
 
    for (unsigned il = 0; il < nl; ++il) {
      unsigned lvl = refinement_levels[il];
      // Initialize solution vectors
      solutions.solutions[il].mu_zero.resize(nk);
      solutions.solutions[il].mu_one.resize(nk);
      solutions.solutions[il].mu_two.resize(nk);
 
      for (int ik = 0; ik < nk; ik++) {
        k_ = ks[ik];
        std::cout << "\nStart computation of refinement_level " << lvl
                  << " and k = " << k_ << std::flush;
 
        // start timer
        std::chrono::steady_clock::time_point start_time =
            std::chrono::steady_clock::now();
 
        // get mesh
        sphere.setRefinementLevel(lvl);
        const std::shared_ptr<lf::mesh::Mesh> mesh = sphere.Build();
 
        // end timer
        std::chrono::steady_clock::time_point end_time =
            std::chrono::steady_clock::now();
        double elapsed_sec =
            std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
                                                                  start_time)
                .count() /
            1000.;
 
        std::cout << " -> Built Mesh " << elapsed_sec << " [s] " << std::flush;
 
        solutions.mesh[lvl] = mesh;
 
        // setup the problem
        lse_builder.SetMesh(mesh);
        lse_builder.SetK(k_);
 
        // start timer
        start_time = std::chrono::steady_clock::now();
 
        // compute the system
        lse_builder.Compute();
 
        // end timer
        end_time = std::chrono::steady_clock::now();
        elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(
                          end_time - start_time)
                          .count() /
                      1000.;
 
        std::cout << " -> Computed LSE " << elapsed_sec << " [s] "
                  << std::flush;
 
        // start timer
        start_time = std::chrono::steady_clock::now();
 
        // solve the system
        lse_builder.Solve();
 
        // end timer
        end_time = std::chrono::steady_clock::now();
 
        elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(
                          end_time - start_time)
                          .count() /
                      1000.;
 
        std::cout << " -> Solved System " << elapsed_sec << " [s] "
                  << std::flush;
 
        // store solutions
        solutions.solutions[lvl].mu_zero[ik] = lse_builder.GetMu(0);
        solutions.solutions[lvl].mu_one[ik] = lse_builder.GetMu(1);
        solutions.solutions[lvl].mu_two[ik] = lse_builder.GetMu(2);
      }  // end loop k
    }    // end loop level
 
    // end timer total
    std::chrono::steady_clock::time_point end_time_total =
        std::chrono::steady_clock::now();
    double elapsed_sec_total =
        std::chrono::duration_cast<std::chrono::milliseconds>(end_time_total -
                                                              start_time_total)
            .count() /
        1000.;
 
    std::cout << "\nTotal computation time for all levels " << elapsed_sec_total
              << " [s]\n";
 
    projects::hldo_sphere::post_processing::process_results<
        decltype(u_zero_), decltype(u_one_), decltype(u_two_), complex>(
        name_, solutions, u_zero_, u_one_, u_two_, k_);
  }
 
 private:
  std::function<complex(const Eigen::Matrix<double, 3, 1> &)> u_zero_;
  std::function<Eigen::Matrix<complex, 3, 1>(
      const Eigen::Matrix<double, 3, 1> &)>
      u_one_;
  std::function<complex(const Eigen::Matrix<double, 3, 1> &)> u_two_;
  std::function<complex(const Eigen::Matrix<double, 3, 1> &)> f_zero_;
  std::function<Eigen::Matrix<complex, 3, 1>(
      const Eigen::Matrix<double, 3, 1> &)>
      f_one_;
  std::function<complex(const Eigen::Matrix<double, 3, 1> &)> f_two_;
  double &k_;
  std::string name_;
};
 
}  // namespace projects::hldo_sphere::experiments
 
#endif  // HLDO_SPHERE_EXPERIMENTS_DIRAC_OPERATOR_EXAMPLE_H