OO2Configs.h

#ifndef OO2CONFIGS_H
#define OO2CONFIGS_H

#include <complex>
#include <string>
#include <gmshfem/ScalarFunction.h>
#include <petsc.h>

#define im PetscScalar(0, 1)

struct OO2 {
  gmshfem::function::ScalarFunction<PetscScalar> S;
  gmshfem::function::ScalarFunction<PetscScalar> dS;
};

class AbstractOO2Config {
  public:
  virtual std::string name() const = 0;
  virtual OO2 coefs(gmshfem::function::ScalarFunction<PetscScalar> k) const = 0;
};

class TaylorOrderZero : public AbstractOO2Config
{
 public:
  double _alpha = 0;
  std::string _name;

  TaylorOrderZero(std::string name, double alpha) :
    _alpha(alpha), _name(name)
  {
  }

  std::string name() const override
  {
    return _name;
  }


  OO2 coefs(gmshfem::function::ScalarFunction<PetscScalar> k) const override
  {
    return {im * k * (PetscScalar(cos(_alpha/2), sin(_alpha/2))), 0.0};
  }
};

class TaylorOrderZeroReal : public AbstractOO2Config
{
 public:
  double _alpha = 0;
  std::string _name;

  TaylorOrderZeroReal(std::string name, double alpha) :
    _alpha(alpha), _name(name)
  {
  }

  std::string name() const override
  {
    return _name;
  }


  OO2 coefs(gmshfem::function::ScalarFunction<PetscScalar> k) const override
  {
    return {im * k * cos(_alpha/2), 0.0};
  }
};


class TaylorOrderTwo : public AbstractOO2Config
{
 public:
  double _alpha = 0;
  std::string _name;

  TaylorOrderTwo(std::string name, double alpha) :
    _alpha(alpha), _name(name)
  {
  }

  std::string name() const override
  {
    return _name;
  }


  OO2 coefs(gmshfem::function::ScalarFunction<PetscScalar> k) const override
  {
    // For operator a + b*laplacian, not b grad grad
    return {im * k * cos(_alpha / 2),
            im / (2.0 * k) * PetscScalar(cos(_alpha / 2), +sin(-_alpha / 2))};
  }
};

class TaylorOrderTwoFlipped : public AbstractOO2Config
{
 public:
  double _alpha = 0;
  std::string _name;

  TaylorOrderTwoFlipped(std::string name, double alpha) :
    _alpha(alpha), _name(name)
  {
  }

  std::string name() const override
  {
    return _name;
  }

  OO2 coefs(gmshfem::function::ScalarFunction<PetscScalar> k) const override
  {
    // Flipped in ORAS
    return {im * k * cos(_alpha / 2),
            -im / (2.0 * k) * PetscScalar(cos(_alpha / 2), +sin(-_alpha / 2))};
  }
};

class SommerfeldSecondOrderDamping : public AbstractOO2Config
{
 public:
  double _beta = 0;
  std::string _name;

  SommerfeldSecondOrderDamping(std::string name, double beta) :
    _beta(beta), _name(name)
  {
  }

  std::string name() const override
  {
    return _name;
  }


  OO2 coefs(gmshfem::function::ScalarFunction<PetscScalar> k) const override
  {
    return {im * k,
            _beta / k};
  }
};

#undef im // Not needed anymore

#endif // OO2CONFIGS_H