pml_formula.hh

Go to the documentation of this file.

 
#ifndef waveprop_pmlformula_hh
#define waveprop_pmlformula_hh
 
 
#include <iostream>
#include <memory>
#include "stdio.h"
#include "basics.hh"
#include "formula/formula.hh"
#include "formula/matrixElementFormula.hh"
#include "geometry/formula.hh"
#include "hp2D.hh"
 
#include "pml_edge.hh"
 
#define pmlHolger_D 0
 
namespace concepts {
 
 
  // **************************************************** FormulaPMLPowerSigma **
  
  // Sigma for Cartesian PML, see Collino & Monk
  
  template<typename F=Real>
  class FormulaPMLPowerSigma : public Formula<F> {
  public:
    FormulaPMLPowerSigma(const Real offset, const int power=2
      , const F sigma0 = 5.0, const Real center=0) 
      : offset_(offset)
      , center_(center)
      , power_(power)
      , sigma0_(sigma0) 
    {}
    
    virtual FormulaPMLPowerSigma* clone() const {
      return new FormulaPMLPowerSigma(offset_, power_, sigma0_, center_);
    }
 
    bool inPMLregion(const Real p, const Real t = 0.0) 
    {
      Real arg = fabs(p - center_) - offset_;
      
      return arg >= 0;
    }
    
    virtual F operator() (const Real p, const Real t = 0.0) const {
      Real arg = fabs(p - center_) - offset_;
      if (arg < 0)
        return F(0);
      return sigma0_ * powi(arg , power_);
    }
    
    virtual F operator() (const Real2d& p, const Real t = 0.0) const {
      return (*this)(p[0], t);
    }
    
    virtual F operator() (const Real3d& p, const Real t = 0.0) const {
      return (*this)(p[0], t);
    }
    template<typename Real>
    static inline Real powi(Real x, int powercoeff) {
#if 1
      Real res(1);
      for(int i = 0; i < powercoeff; ++i) {
        res *= x;
      }
      return res;
#else
      return pow(x, powercoeff);
#endif
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"(" << sigma0_ << " * ( |x - (" << center_
                << ")| - " << offset_ << ")_+^" << power_ << ")))";
    }
  private:
    const Real offset_;
    const Real center_;
    const int  power_;
    const F    sigma0_; // PML strength
  };
  
  // ************************************************** FormulaPMLPowerSigma2D **
  
  template<typename F=Real>
  class FormulaPMLPowerSigma2D : public Formula<F> {
  public:
    FormulaPMLPowerSigma2D(const Real offset, const int power = 2, 
                           const F sigma0 = 5.0, 
                           const Real2d& center = Real2d(0,0))
      : offset_(offset), center_(center), power_(power), sigma0_(sigma0) {
    }
    
    virtual FormulaPMLPowerSigma2D* clone() const {
      return new FormulaPMLPowerSigma2D(offset_, power_, sigma0_, center_);
    }
    
    virtual F operator() (const Real p, const Real t = 0.0) const {
      conceptsThrowSimpleE("FormulaPMLPowerSigma2D currently only supports 2D!");
      assert(false);
    }
 
    inline bool inPMLregion(const Real2d& p, const Real t = 0.0) const
    {
 
      Real arg = (p-center_).l2() - offset_;
      return arg >= 0;
    }
    
    virtual F operator() (const Real2d& p, const Real t = 0.0) const {
      Real dist = (p-center_).l2() - offset_;
      if (dist < 0)
        return F(0);
      return sigma0_ * powi(dist, power_);
    }
    
    virtual F operator() (const Real3d& p, const Real t = 0.0) const {
      return (*this)( Real2d(p[0], p[1]), t);
    } 
    
    template<typename Real>
    static inline Real powi(Real x, int powercoeff) {
#if 1
      Real res(1);
      for (int i = 0; i < powercoeff; ++i) {
        res *= x;
      }
      return res;
#else
      return pow(x, powercoeff);
#endif
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << concepts::typeOf(*this)<<"(" << sigma0_ << " * ( |x - (" << center_[0]
                << "," << center_[1] << ")| - " << offset_ << ")_+^" << power_ << ")))";
    }
  private:
    const Real    offset_;
    const Real2d center_;
    // PML exponent
    const int     power_;
    // PML strength
    const F       sigma0_;
  };
  
  // ************************************************* FormulaPMLPowerSigmaB2D **
 
  template<typename F=Real>
  class FormulaPMLPowerSigmaB2D : public Formula<F> {
  public:
    FormulaPMLPowerSigmaB2D(const Real offset, const int power = 2, 
                            const F sigma0 = 5.0, 
                            const Real2d& center = Real2d(0,0))
      : offset_(offset), center_(center), power_(power), sigma0_(sigma0) {}
    
    virtual FormulaPMLPowerSigmaB2D* clone() const {
      return new FormulaPMLPowerSigmaB2D(offset_, power_, sigma0_, center_);
    }
    
    virtual F operator() (const Real p, const Real t = 0.0) const {
      conceptsThrowSimpleE
        ("FormulaPMLPowerSigmaB2D currently only supports 2D!"); assert(false);
    }
 
    inline bool inPMLregion(const Real2d& p, const Real t = 0.0) const
    {
      Real arg = (p-center_).l2() - offset_;
      return arg >= 0;
    }
    
    virtual F operator() (const Real2d& p, const Real t = 0.0) const {
      Real rho = (p-center_).l2();
      Real arg = rho - offset_;
      if (arg<0)
        return F(0);
      return sigma0_ * powi( arg , power_ + 1) / (rho) / (power_ + 1);
    }
    
    virtual F operator() (const Real3d& p, const Real t = 0.0) const {
      return (*this)( Real2d(p[0], p[1]), t);
    } 
    
    template<typename Real>
    static inline Real powi(Real x, int powercoeff) {
#if 1
      Real res(1);
      for (int i = 0; i < powercoeff; ++i) {
        res *= x;
      }
      return res;
#else
      return pow(x, powercoeff);
#endif
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"(" << sigma0_ << " * ( |x - ("
                << center_[0] << "," << center_[1] << ")| - " 
                << offset_ << ")_+^" << power_ << ")))";
    }
  private:
    const Real    offset_;  
    const Real2d  center_;
    // PML exponent
    const int     power_;
    // PML strength  
    const F       sigma0_;
  };
  
  // ********************************************************** FomulaPMLCart **
  
  class FormulaPMLCart : public ElementFormula<Cmplx> {
  public:
    enum PMLMode { DXDX, DYDY, IDENT, DX, DY };
    
    FormulaPMLCart(const ElementFormulaContainer<Cmplx> coeff_a
                 , const ElementFormulaContainer<Cmplx> coeff_b
                 , RCP<Formula<Real> > sigma_x
                 , RCP<Formula<Real> > sigma_y
                 , PMLMode mode
                 , double omega
                   ) 
      : coeff_a_(coeff_a)
      , coeff_b_(coeff_b)
      , sigma_x_(sigma_x)
      , sigma_y_(sigma_y)
      , mode_(mode)
      , omega_(omega)
    {}
    
    virtual FormulaPMLCart* clone() const {
      return new FormulaPMLCart(*this);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real3d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real2d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    Cmplx gammaX(Real x) const {
      return Cmplx(1, (*sigma_x_)(x) / omega_);
    }
    
    Cmplx gammaY(Real y) const {
      return Cmplx(1, (*sigma_y_)(y) / omega_);
    }
    
    template <class RealNd>
    inline Cmplx operator() (const ElementWithCell< Real > &elm, const RealNd &p, Real2d px, 
                             const Real t=0.0) const
    {
      switch(mode_) {
      case DXDX:  return gammaY(px[1]) / gammaX(px[0]) * coeff_a_(elm, p);
      case DX:    return gammaY(px[1]) * coeff_a_(elm, p);
      case DYDY:  return gammaX(px[0]) / gammaY(px[1]) * coeff_a_(elm, p);
      case DY:    return gammaX(px[0]) * coeff_a_(elm, p);
      case IDENT: return gammaX(px[0]) * gammaY(px[1]) * coeff_b_(elm, p);
      }
      conceptsThrowSimpleE("FormulaPMLCart has gone mad!"); assert(false);
    }
    
    
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"(coeffs: "
                << coeff_a_ << ", " << coeff_b_ << ", sigmas:" 
                << *sigma_x_ << ", " << *sigma_y_ 
                << ")";
    }
  private:
    const ElementFormulaContainer<Cmplx> coeff_a_;
    const ElementFormulaContainer<Cmplx> coeff_b_;
    RCP<Formula<Real> > sigma_x_;
    RCP<Formula<Real> > sigma_y_;
    PMLMode mode_;
    double omega_;
  };
  
  // ************************************************ FormulaPMLBoxRestriction **
  
  template <class F, typename G = typename Realtype<F>::type>
  class FormulaPMLBoxRestriction : public ElementFormula<F, G> {
  public:
    FormulaPMLBoxRestriction(
                             RCP<FormulaPMLPowerSigma<Real> > sigma_x
                           , RCP<FormulaPMLPowerSigma<Real> > sigma_y
                           , RCP<ElementFormula<F, G> > formula
                             )
      : sigma_x(sigma_x)
      , sigma_y(sigma_y)
      , formula(formula)
    { }
    
    virtual F operator() (const ElementWithCell< Real > &elm, const Real3d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    virtual F operator() (const ElementWithCell< Real > &elm, const Real2d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    virtual F operator() (const ElementWithCell< Real > &elm, const Real p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    template <class RealNd>
    F operator() (const ElementWithCell< Real > &elm, const RealNd &p, Real2d px, 
                  const Real t=0.0) const
    {
      if ( sigma_x->inPMLregion(px[0], t) || sigma_y->inPMLregion(px[1], t)) {
        return formula->operator()(elm, p, t);
      } else {
        //cout << "out: " << px << p << formula->operator()(elm, p, t) << endl;
        
        return F(0);
        //return formula->operator()(elm, p, t);
      }
    }
    
    virtual FormulaPMLBoxRestriction* clone() const {
      return new FormulaPMLBoxRestriction(*this);
    }
    
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"(" << formula << " :" << *sigma_x << ", " << sigma_y << ")))";
    }
  private:
    RCP<FormulaPMLPowerSigma<Real> > sigma_x;
    RCP<FormulaPMLPowerSigma<Real> > sigma_y;
    RCP<ElementFormula<F> > formula;
  };
  
  // ******************************************************** FormulaPMLRadia **
  
  class FormulaPMLRadia : public ElementFormula<Cmplx> {
  public:
    enum PMLMode { AD1, AD2, AS1, AS2, IDENT};
    
    FormulaPMLRadia(RCP<Formula<Real> > sigma, RCP<Formula<Real> > sigmaB,
                    PMLMode mode, const Real2d& center = Real2d(0,0));
    
    virtual FormulaPMLRadia* clone() const {
      return new FormulaPMLRadia(*this);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm,
                              const Real p, const Real t = 0.0) const;
    virtual Cmplx operator() (const ElementWithCell< Real > &elm,
                              const Real2d &p, const Real t = 0.0) const;
    virtual Cmplx operator() (const ElementWithCell< Real > &elm,
                              const Real3d &p, const Real t = 0.0) const;
 
    inline RCP<Formula<Real> > sigma()  const { return sigma_;  }
    inline RCP<Formula<Real> > sigmaB() const { return sigmaB_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    RCP<Formula<Real> > sigma_;
    RCP<Formula<Real> > sigmaB_;
    PMLMode mode_;
    Real2d center_;
    
 
    inline Cmplx gamma_(Real2d &p) const{
      return Cmplx(1, (*sigma_) (p));
    }
    inline Cmplx gammaB_(Real2d &p) const{
      return Cmplx(1, (*sigmaB_)(p));
    }    
  };
 
  // ******************************************************* RadialPMLFormulas **
 
  class RadialPMLFormulas : public OutputOperator {
  public:
    RadialPMLFormulas(const Real offset, const int power = 2, 
                      const Real sigma0 = 5.0, 
                      const Real2d& center = Real2d(0,0));
 
    inline ElementFormulaContainer<MapCmplx2d> A() const {
      return A_;
    }
    inline ElementFormulaContainer<Cmplx> M() const {
      return M_;
    }
 
    inline RCP<Formula<Real> > sigma()  const { return M_.sigma();  }
    inline RCP<Formula<Real> > sigmaB() const { return M_.sigmaB(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    MatrixElementFormula<Cmplx, 2> A_;
    FormulaPMLRadia M_;
  };
 
 
 
 
  // ******************************************************* new cartesian PML **
  
  class FormulaPMLCartNew : public ElementFormula<Cmplx> {
  public:
    enum PMLMode { DXDX, DYDY, IDENT, DX, DY , ZERO};
    
    FormulaPMLCartNew(  RCP<Formula<Real> > sigma_x, 
                        RCP<Formula<Real> > sigma_y, 
                        PMLMode mode ) 
      : sigma_x_(sigma_x), 
        sigma_y_(sigma_y),
        mode_(mode)
    {}
    
    virtual FormulaPMLCartNew* clone() const {
      return new FormulaPMLCartNew(*this);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real3d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real2d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    inline RCP<Formula<Real> > sigmax() const { return sigma_x_; }
    inline RCP<Formula<Real> > sigmay() const { return sigma_y_; }
 
    Cmplx gammaX(Real x) const {
      return Cmplx(1., (*sigma_x_)(x) );
    }
    
    Cmplx gammaY(Real y) const {
      return Cmplx(1., (*sigma_y_)(y) );
    }
    
    template <class RealNd>
    inline Cmplx operator() (const ElementWithCell< Real > &elm, const RealNd &p, Real2d px, const Real t=0.0) const
    {
      switch(mode_) {
      case DXDX:  return gammaY(px[1]) / gammaX(px[0]);
      case DX:    return gammaY(px[1]) ;
      case DYDY:  return gammaX(px[0]) / gammaY(px[1]) ;
      case DY:    return gammaX(px[0]) ;
      case ZERO : return 0.0;
      case IDENT: return gammaX(px[0]) * gammaY(px[1]) ;
      }
      conceptsThrowSimpleE("FormulaPMLCart has gone mad!"); assert(false);
    }
    
    
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"( sigma x ="
                << *sigma_x_ << ", sigma y = " << *sigma_y_ 
                << ")";
    }
  private:
    RCP<Formula<Real> > sigma_x_;
    RCP<Formula<Real> > sigma_y_;
    PMLMode mode_;
 
  };
 
  class CartesianPMLFormulas : public OutputOperator {
    
  private:
    MatrixElementFormula<Cmplx, 2> A_;
    FormulaPMLCartNew M_;
    
  public:
    CartesianPMLFormulas(const Real offset_x, const Real offset_y, 
       const int power = 2,                       
       const Real sigma0 = 5.0, 
       const Real center_x =0,
       const Real center_y = 0);
 
    inline ElementFormulaContainer<MapCmplx2d> A() const {
      return A_;
    }
    inline ElementFormulaContainer<Cmplx> M() const {
      return M_;
    }
 
    inline RCP<Formula<Real> > sigmax() const { return M_.sigmax(); }
    inline RCP<Formula<Real> > sigmay() const { return M_.sigmay(); }
  
    friend std::ostream& operator<<(std::ostream& out, const CartesianPMLFormulas &ca)
    {
      return(ca.info(out));
    } 
protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
 
 
  // *********************************************************** hamburger PML **
  class FormulaPMLHamburger : public ElementFormula<Cmplx>{
  public:
    enum PMLMode {M1, M2, M3, M4, IDENT};
  private:
    Real R;
    Real d;
    int power;
    Real sigma0;
    Real center_x, center_y;
    Real2d center_point_up, center_point_down;
    PMLMode mode;
    FormulaPMLCartNew Cart;
    FormulaPMLRadia RadiaUp;
    FormulaPMLRadia RadiaDown;
 
    FormulaPMLCartNew::PMLMode convert_mode_to_cart(FormulaPMLHamburger::PMLMode mode)
    {
      // modeCart 
      switch (mode){
      case M1: return(FormulaPMLCartNew::DXDX);
      case M2: return(FormulaPMLCartNew::ZERO);
      case M3: return(FormulaPMLCartNew::ZERO);
      case M4: return(FormulaPMLCartNew::DYDY);
      case IDENT: return(FormulaPMLCartNew::IDENT);
      }
      return(FormulaPMLCartNew::IDENT);
    }
    FormulaPMLRadia::PMLMode convert_mode_to_radia(FormulaPMLHamburger::PMLMode mode)
    {
      // modeRadia 
      switch (mode){
      case M1: return(FormulaPMLRadia::AD1);
      case M2: return(FormulaPMLRadia::AS1);
      case M3: return(FormulaPMLRadia::AS2);
      case M4: return(FormulaPMLRadia::AD2);
      case IDENT: return(FormulaPMLRadia::IDENT);
      }
      return(FormulaPMLRadia::IDENT);
    }
  public:
    FormulaPMLHamburger(const Real R_, 
                        const Real d_,
                        const int power_, 
                        const Real sigma0_,
                        const Real center_x_, 
                        const Real center_y_,
                        PMLMode mode_)
      : R(R_), 
        d(d_), 
        power(power_), 
        sigma0(sigma0_), 
        center_x(center_x_), 
        center_y(center_y_), 
        center_point_up(center_x,center_y+d/2.),
        center_point_down(center_x,center_y-d/2.),
        mode(mode_),
        Cart(
            makeRCP(new FormulaPMLPowerSigma<Real>(R, power, sigma0, center_x)), 
            makeRCP(new ConstFormula<Real>(0)),
            convert_mode_to_cart(mode)),
        RadiaUp(
          makeRCP(new FormulaPMLPowerSigma2D<Real>(R, power, sigma0)), 
          makeRCP(new FormulaPMLPowerSigmaB2D<Real>(R, power,sigma0)), 
          convert_mode_to_radia(mode),
          center_point_up),
        RadiaDown(
          makeRCP(new FormulaPMLPowerSigma2D<Real>(R, power, sigma0)), 
          makeRCP(new FormulaPMLPowerSigmaB2D<Real>(R, power,sigma0)), 
          convert_mode_to_radia(mode),
          center_point_down)
 
    {
 
    }
   
    virtual FormulaPMLHamburger* clone() const {
      return new FormulaPMLHamburger(*this);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real3d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real2d &p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real p, const Real t=0.0) const
    {
      Real2d px = elm.elemMap(p);
      return operator()(elm, p, px, t);
    }
    
    
    template <class RealNd>
    inline Cmplx operator() (const ElementWithCell< Real > &elm, const RealNd &p, Real2d px, const Real t=0.0) const
    {
      
      // px is the point in global coordinates
      if (px[1]>d/2.){
        // Upper "Bread" Part
        return RadiaUp(elm, p, t);
      }
      else if (px[1]< -d/2.){
        // Lower "Bread" Part
        return(RadiaDown(elm, p, t));
      }
      else {
        // "Steak" Part
        return(Cart(elm, p, px, t));
      }
    }
    
 
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"( sigma x ="
                << *Cart.sigmax() << ", sigma y = " << *Cart.sigmay()
                << ", sigma up = "<< *RadiaUp.sigma() 
                << ", sigma down = " << *RadiaDown.sigma() << ")";
    }
  };
 
  class HamburgerPMLFormulas : public OutputOperator {
    
  private:
    Real R_, d_;
    int power_;
    Real sigma0_, center_x_, center_y_;
    MatrixElementFormula<Cmplx, 2> A_;
    FormulaPMLHamburger M_;
    
  public:
    HamburgerPMLFormulas(const Real R, 
                         const Real d, 
                         const int power = 2,                       
                         const Real sigma0 = 5.0, 
                         const Real center_x =0,
                         const Real center_y = 0);
 
    inline ElementFormulaContainer<MapCmplx2d> A() const {
      return A_;
    }
    inline ElementFormulaContainer<Cmplx> M() const {
      return M_;
    }
 
 
  
    friend std::ostream& operator<<(std::ostream& out, const HamburgerPMLFormulas &ca)
    {
      return(ca.info(out));
    } 
protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
 
} // namespace concepts
 
#endif // waveprop_pmlformula_hh