dirac_convergence_test.cc

 
#include <dirac_convergence_test.h>
 
using lf::uscalfe::operator-;
using complex = std::complex<double>;
 
int main(int argc, char *argv[]) {
  if (argc != 3) {
    std::cerr << "Usage: " << argv[0] << " max_refinement_level k" << std::endl;
    exit(1);
  }
 
  double kt = 0.5;
  double &k = kt;
  unsigned refinement_level = atoi(argv[1]);
  std::sscanf(argv[2], "%lf", &k);
  std::cout << "max_refinement_level : " << refinement_level << std::endl;
  std::cout << "k : " << k << std::endl;
 
  // mathematica function output requries the following helpers
  auto Power = [](complex a, complex b) -> complex { return std::pow(a, b); };
  auto Complex = [](double a, double b) -> complex { return complex(a, b); };
  auto Sin = [](complex a) -> complex { return std::sin(a); };
  auto Cos = [](complex a) -> complex { return std::cos(a); };
  auto Sqrt = [](complex a) -> complex { return std::sqrt(a); };
 
  // righthandside for the zero and two form
  auto f_zero = [&](const Eigen::Vector3d &x_vec) -> complex {
    // first scale to the circle
    Eigen::Vector3d x_ = x_vec / x_vec.norm();
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    // autogenerated by mathematica
    complex ret = Complex(0, 1) * k * (Cos(y) + Sin(x) + Sin(z)) +
                  (x * z * Cos(x) + x * y * Cos(y) + y * z * Cos(z) +
                   2 * z * Sin(x) + 2 * x * Sin(y) + 2 * y * 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::VectorXcd {
    // first scale to the circle
    Eigen::Vector3d x_ = x_vec / x_vec.norm();
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    // autogenerated by mathematica
    Eigen::VectorXcd ret(3);
    ret << ((Power(y, 2) + Power(z, 2)) * Cos(x) - x * z * Cos(z) +
            x * y * Sin(y) +
            Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
                (y * Cos(z) + z * Sin(y)) -
            Complex(0, 1) * k *
                (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)),
        -((x * y * Cos(x) + y * z * Cos(z) +
           Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
               (-(z * Cos(x)) + x * Cos(z)) +
           (Power(x, 2) + Power(z, 2)) * Sin(y) +
           Complex(0, 1) * k *
               (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))),
        (-(x * z * Cos(x)) + (Power(x, 2) + Power(y, 2)) * Cos(z) +
         y * z * Sin(y) -
         Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
             (y * Cos(x) + x * Sin(y)) +
         Complex(0, 1) * k *
             ((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;
  };
 
  // righthandside for the two form
  auto f_two = [&](const Eigen::Vector3d &x_vec) -> complex {
    // first scale to the circle
    Eigen::Vector3d x_ = x_vec / x_vec.norm();
    double x = x_(0);
    double y = x_(1);
    double z = x_(2);
 
    // autogenerated by mathematica
    complex ret = (y * Cos(x) + z * Cos(y) + x * Cos(z)) /
                      Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) +
                  Complex(0, 1) * k * (Cos(y) + Sin(x) + Sin(z));
    return ret;
  };
 
  projects::hldo_sphere::debugging::DiracConvergenceTest test(f_zero, f_one,
                                                              f_two, k);
 
  test.Compute(refinement_level);
 
  return 0;
}