hodge_laplacian_test.cc

 
#include <hodge_laplacian_experiment.h>
 
int main(int argc, char *argv[]) {
  if (argc != 4 && argc != 5) {
    std::cerr << "Usage: " << argv[0]
              << " max_refinement_level min_k max_k step=0.1 " << std::endl;
    exit(1);
  }
 
  double min_k;
  double max_k;
  double step = 0.1;
  const unsigned refinement_level = atoi(argv[1]);
  sscanf(argv[2], "%lf", &min_k);
  sscanf(argv[3], "%lf", &max_k);
  if (argc == 5) sscanf(argv[4], "%lf", &step);
  std::cout << "max_refinement_level : " << refinement_level << std::endl;
  std::cout << "min_k : " << min_k << ", max_k = " << max_k
            << " step = " << step << std::endl;
 
  // needs to be a reference such that it can be reassigned
  double &k = min_k;
 
  std::vector<unsigned> refinement_levels(refinement_level + 1);
  for (int i = 0; i < refinement_level + 1; i++) {
    refinement_levels[i] = i;
  }
 
  std::vector<double> ks;
  for (double i = min_k; i <= max_k; i += step) {
    ks.push_back(sqrt(i));
  }
 
  // mathematica function output requries the following helpers
  auto Power = [](double a, double b) -> double { return std::pow(a, b); };
  auto Sin = [](double a) -> double { return std::sin(a); };
  auto Cos = [](double a) -> double { return std::cos(a); };
  auto Sqrt = [](double a) -> double { return std::sqrt(a); };
 
  double r = 1.;
 
  // righthandside for the zero and two form
  auto f_zero = [&](const Eigen::Vector3d &x_vec) -> double {
    // first scale to the sphere
    Eigen::Vector3d x_ = x_vec;
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    // autogenerated by mathematica
    double ret =
        Power(k, 2) * (Cos(z) + Sin(x) + Sin(y)) +
        (2 * x * Cos(x) + 2 * y * Cos(y) + Power(x, 2) * Cos(z) +
         Power(y, 2) * Cos(z) + Power(y, 2) * Sin(x) + Power(z, 2) * Sin(x) +
         Power(x, 2) * Sin(y) + Power(z, 2) * Sin(y) - 2 * z * Sin(z)) /
            (Power(x, 2) + Power(y, 2) + Power(z, 2));
 
    return ret;
  };
 
  // righthandside for the one form
  auto f_one = [&](const Eigen::Vector3d &x_vec) -> Eigen::VectorXd {
    // first scale to the circle
    Eigen::Vector3d x_ = x_vec;
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    // autogenerated by mathematica
    Eigen::VectorXd ret(3);
    ret << (2 * z * (-Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(x) +
            2 * y * (-Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(y) -
            4 * x * y * z * Cos(z) - 2 * x * z * Sin(x) -
            x * Power(y, 2) * z * Sin(x) - x * Power(z, 3) * Sin(x) +
            2 * Power(y, 2) * Sin(y) + Power(x, 2) * Power(y, 2) * Sin(y) +
            2 * Power(z, 2) * Sin(y) + Power(x, 2) * Power(z, 2) * Sin(y) +
            Power(y, 2) * Power(z, 2) * Sin(y) + Power(z, 4) * Sin(y) -
            2 * x * y * Sin(z) - Power(x, 3) * y * Sin(z) -
            x * Power(y, 3) * Sin(z) +
            Power(k, 2) * (Power(x, 2) + Power(y, 2) + Power(z, 2)) *
                (-(x * z * Sin(x)) + (Power(y, 2) + Power(z, 2)) * Sin(y) -
                 x * y * Sin(z))) /
               Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 2),
        (-4 * x * y * z * Cos(x) +
         2 * x * (Power(x, 2) - Power(y, 2) + Power(z, 2)) * Cos(y) +
         2 * Power(x, 2) * z * Cos(z) - 2 * Power(y, 2) * z * Cos(z) +
         2 * Power(z, 3) * Cos(z) - 2 * y * z * Sin(x) -
         Power(y, 3) * z * Sin(x) - y * Power(z, 3) * Sin(x) -
         2 * x * y * Sin(y) - Power(x, 3) * y * Sin(y) -
         x * y * Power(z, 2) * Sin(y) + 2 * Power(x, 2) * Sin(z) +
         Power(x, 4) * Sin(z) + Power(x, 2) * Power(y, 2) * Sin(z) +
         2 * Power(z, 2) * Sin(z) + Power(x, 2) * Power(z, 2) * Sin(z) +
         Power(y, 2) * Power(z, 2) * Sin(z) +
         Power(k, 2) * (Power(x, 2) + Power(y, 2) + Power(z, 2)) *
             (-(y * z * Sin(x)) - x * y * Sin(y) +
              (Power(x, 2) + Power(z, 2)) * Sin(z))) /
            Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 2),
        (2 * x * (Power(x, 2) + Power(y, 2) - Power(z, 2)) * Cos(x) -
         4 * x * y * z * Cos(y) + 2 * Power(x, 2) * y * Cos(z) +
         2 * Power(y, 3) * Cos(z) - 2 * y * Power(z, 2) * Cos(z) +
         2 * Power(x, 2) * Sin(x) + 2 * Power(y, 2) * Sin(x) +
         Power(x, 2) * Power(y, 2) * Sin(x) + Power(y, 4) * Sin(x) +
         Power(x, 2) * Power(z, 2) * Sin(x) +
         Power(y, 2) * Power(z, 2) * Sin(x) - 2 * x * z * Sin(y) -
         Power(x, 3) * z * Sin(y) - x * Power(z, 3) * Sin(y) -
         2 * y * z * Sin(z) - Power(x, 2) * y * z * Sin(z) -
         Power(y, 3) * z * Sin(z) +
         Power(k, 2) * (Power(x, 2) + Power(y, 2) + Power(z, 2)) *
             ((Power(x, 2) + Power(y, 2)) * Sin(x) -
              z * (x * Sin(y) + y * Sin(z)))) /
            Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 2);
    return ret;
  };
 
  // Compute the analytic solution of the problem
  auto u_zero = [&](const Eigen::Vector3d x_vec) -> double {
    // first scale to the circle
    Eigen::Vector3d x_ = x_vec;
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    return Cos(z) + Sin(x) + Sin(y);
  };
 
  // Compute the analytic solution of the problem
  auto u_one = [&](const Eigen::Vector3d x_vec) -> Eigen::Vector3d {
    // first scale to the circle
    Eigen::Vector3d x_ = x_vec / x_vec.norm() * r;
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    Eigen::VectorXd ret(3);
 
    // mathematica autocompute
    ret << (-(x * z * Sin(x)) + (Power(y, 2) + Power(z, 2)) * Sin(y) -
            x * y * Sin(z)) /
               (Power(x, 2) + Power(y, 2) + Power(z, 2)),
        (-(y * z * Sin(x)) - x * y * Sin(y) +
         (Power(x, 2) + Power(z, 2)) * Sin(z)) /
            (Power(x, 2) + Power(y, 2) + Power(z, 2)),
        ((Power(x, 2) + Power(y, 2)) * Sin(x) - z * (x * Sin(y) + y * Sin(z))) /
            (Power(x, 2) + Power(y, 2) + Power(z, 2));
 
    return ret;
  };
 
  projects::hldo_sphere::experiments::HodgeLaplacianExperiment<double>
      experiment(u_zero, u_one, u_zero, f_zero, f_one, f_zero, k, "test");
 
  experiment.Compute(refinement_levels, ks);
 
  return 0;
}