pml_formula.h

#pragma once
 
#include <iostream>
#include <memory>
#include "stdio.h"
#include "basics.hh"
#include "formula/elementFormulaRCP.hh"
#include "formula/formula.hh"
#include "geometry/formula.hh"
#include "hp2D.hh"
 
//TODO: properly expand this list
 
#define pmlHolger_D 0
 
namespace concepts {
 
// ******************************************************** FormulaExpImag1D **
 
class FormulaExpImag1D : public concepts::Formula<Cmplx> {
public:
  FormulaExpImag1D(const Real k, const Cmplx u = 1.0, Real x0 = 0.0) 
    : k_(k), x0_(x0), u_(u) {}
 
  virtual FormulaExpImag1D* clone() const {
    return new FormulaExpImag1D(k_, u_, x0_);
  }
 
  virtual Cmplx operator() (const Real p,    const Real t = 0.0) const {
    const Real arg = k_*(p-x0_);
    return Cmplx(cos(arg), sin(arg)) * u_;
  }
  virtual Cmplx operator() (const concepts::Real2d& p, const Real t = 0.0) const {
    return (*this)(p[0], t);
  }
  virtual Cmplx operator() (const concepts::Real3d& p, const Real t = 0.0) const {
    return (*this)(p[0], t);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "FormulaExpImag1D(" << u_ << " * exp(i " << k_ << "(x - " << x0_ << ")))";
  }
private:
  const Real k_, x0_;
  const Cmplx u_;
};
 
// ******************************************************** FormulaExpImag2D **
 
class FormulaExpImag2D : public concepts::Formula<Cmplx> {
public:
  FormulaExpImag2D(const Real2d k, const Cmplx u = 1.0, Real2d x0 = 0.0) 
    : k(k), x0(x0), u(u) {}
 
  virtual FormulaExpImag2D* clone() const {
    return new FormulaExpImag2D(k, u, x0);
  }
 
  virtual Cmplx operator() (const Real p,    const Real t = 0.0) const {
    conceptsThrowSimpleE("FormulaExpImag2D currently only supports 2D!"); assert(false);
  }
 
  virtual Cmplx operator() (const concepts::Real2d& p, const Real t = 0.0) const 
  {
    const Real arg = k[0]*(p[0]-x0[0]) 
      +              k[1]*(p[1]-x0[1]) ;
    //cout << "e^{i k x}: " << p << ", " << Cmplx(cos(arg), sin(arg)) * u << endl;
    return Cmplx(cos(arg), sin(arg)) * u;
  }
 
  virtual Cmplx operator() (const concepts::Real3d& p, const Real t = 0.0) const {
    return (*this)( Real2d(p[0], p[1]), t);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "FormulaExpImag2D(" << u << " * exp(i " << k << "(x - " << x0 << ")))";
  }
public:
  const Real2d k, x0;
  const Cmplx u;
};
 
// *********************************************** FormulaExpImag2DRadialDer **
 
class FormulaExpImag2DRadialDer : public concepts::Formula<Cmplx> {
public:
  //FormulaExpImag2DRadialDer(const Real2d k, const Cmplx u = 1.0, Real2d x0 = Real2d(0,0)) 
  FormulaExpImag2DRadialDer(const Real2d k, const Cmplx u = 1.0) 
    : k(k), x0(0,0), u(u) {}
 
  virtual FormulaExpImag2DRadialDer* clone() const {
    return new FormulaExpImag2DRadialDer(*this);
  }
 
  virtual Cmplx operator() (const Real p,    const Real t = 0.0) const {
    conceptsThrowSimpleE("FormulaExpImag2DRadialDer currently only supports 2D!"); assert(false);
  }
 
  virtual Cmplx operator() (const concepts::Real2d& p, const Real t = 0.0) const 
  {
    const Cmplx imag(0,1);
 
    const Real arg = k[0]*(p[0]-x0[0]) 
      +              k[1]*(p[1]-x0[1]) ;
    Cmplx val(cos(arg), sin(arg));
    val *= u;
 
    const Real len = sqrt( p[0]*p[0] + p[1]*p[1] );
 
    return Cmplx(0, (k[0]*p[0] + k[1]*p[1])/len) * val;
  }
 
  virtual Cmplx operator() (const concepts::Real3d& p, const Real t = 0.0) const {
    return (*this)( Real2d(p[0], p[1]), t);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "FormulaExpImag2DRadialDer(" << u << " * exp(i " << k << "(x - " << x0 << ")))";
  }
public:
  const Real2d k, x0;
  const Cmplx u;
};
 
// **************************************************** FormulaExpImag2DGrad **
 
class FormulaExpImag2DGrad : public concepts::Formula< concepts::Point<Cmplx, 2> > 
{
public:
  FormulaExpImag2DGrad(const Real2d k, const Cmplx u = 1.0, Real2d x0 = Real2d(0,0)) 
    : k(k), x0(x0), u(u) 
  { }
 
  virtual FormulaExpImag2DGrad* clone() const {
    return new FormulaExpImag2DGrad(k, u, x0);
  }
 
  virtual Cmplx2d operator() (const Real p,    const Real t = 0.0) const {
    conceptsThrowSimpleE("FormulaExpImag2DGrad currently only supports 2D!"); assert(false);
  }
 
  virtual Cmplx2d operator() (const concepts::Real2d& p, const Real t = 0.0) const 
  {
    const Cmplx imag(0,1);
 
    const Real arg = k[0]*(p[0]-x0[0]) 
      +              k[1]*(p[1]-x0[1]) ;
    Cmplx val(cos(arg), sin(arg));
    val *= u;
 
    //const Real len = sqrt( p[0]*p[0] + p[1]*p[1] );
 
    return Cmplx2d( Cmplx(0, k[0])*val, Cmplx(0, k[1])*val );
  }
 
  virtual Cmplx2d operator() (const concepts::Real3d& p, const Real t = 0.0) const {
    //conceptsThrowSimpleE("FormulaExpImag2DGrad currently only supports 2D!"); assert(false);
    return (*this)( Real2d(p[0], p[1]), t);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "FormulaExpImag2DGrad(" << u << " * exp(i " << k << "(x - " << x0 << ")))";
  }
public:
  const Real2d k, x0;
  const Cmplx u;
};
 
 
// ************************************************** FormulaNormalOuterSP2D **
 
template< class F >
class FormulaNormalOuterSP2D : public concepts::ElementFormula< F > 
{
public:
  static const bool OUTER = true;
  static const bool INNER = false;
 
  FormulaNormalOuterSP2D( 
      const ElementFormulaContainer< concepts::Point<F, 2> > vf, 
      Real2d center = Real2d(0,0), bool direction = OUTER) 
    : vf(vf)
    , center(center)
    , direction(direction)
  { }
 
  virtual FormulaNormalOuterSP2D* clone() const {
    return new FormulaNormalOuterSP2D(*this);
  }
 
  virtual F operator() (const concepts::ElementWithCell<Real>& elm,
      const Real p,    const Real t = 0.0) const 
  {
    const concepts::Cell& c = elm.cell();
 
    const Edge2d* cc = dynamic_cast<const concepts::Edge2d*> (&c);
    assert(cc);
 
    Real2d normal = cc->jacobian(p);
    normal.ortho();
    normal /= normal.l2();
 
    Real2d pp = elm.elemMap(p); //pp /= pp.l2();
    pp -= center;
 
    if ( ( normal.dot(pp) < 0 && direction == OUTER ) ||
         ( normal.dot(pp) > 0 && direction == INNER )
        ) 
    {  
      normal *= -1;
    }
 
    const concepts::Point<F, 2> u = vf.operator()(elm, p, t);
  
    //std::cout << "normal " << elm.elemMap(p) << ": " << normal 
    //  << "SP: " << normal[0] * u[0] + normal[1] * u[1] << "center: " << center << std::endl;
 
    return normal[0] * u[0] + normal[1] * u[1];
  }
 
  virtual F operator() (const concepts::ElementWithCell<Real>& elm, 
      const concepts::Real2d& p, const Real t = 0.0) const 
  {
    conceptsThrowSimpleE("FormulaNormalOuterSP2D currently only supports 1D!"); assert(false);
  }
 
  virtual F operator() (const concepts::ElementWithCell<Real>& elm,
      const concepts::Real3d& p, const Real t = 0.0) const {
    conceptsThrowSimpleE("FormulaNormalOuterSP2D currently only supports 1D!"); assert(false);
    //return (*this)( Real2d(p[0], p[1]), t);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "FormulaNormalOuterSP2D(" << vf << ", " << center << ", " << direction << ")";
  }
public:
  const ElementFormulaContainer< concepts::Point<F, 2> > vf;
  Real2d center;
  bool direction;
};
 
// **************************************************** ComposeFormulaMatVec **
 
template< class F, uint DIM, typename G = typename Realtype<F>::type>
class ComposeFormulaMatVec : public concepts::ElementFormula< concepts::Point<F, DIM>, G > 
{
public:
  ComposeFormulaMatVec( 
      RCP<const ElementFormula< concepts::Mapping<F, DIM>, G> > A,
      RCP<const ElementFormula< concepts::Point<F, DIM> , G> > vf
      ) 
    : A(A)
    , vf(vf)
  { }
 
  virtual ComposeFormulaMatVec* clone() const {
    return new ComposeFormulaMatVec(*this);
  }
 
  virtual concepts::Point<F, DIM> operator() (const concepts::ElementWithCell<G>& elm,
      const Real p,    const Real t = 0.0) const 
  {
    return A->operator()(elm, p, t) * vf->operator()(elm, p, t);
  }
 
  virtual concepts::Point<F, DIM> operator() (const concepts::ElementWithCell<G>& elm, 
      const concepts::Real2d& p, const Real t = 0.0) const 
  {
    //cout << "p: " << p << A->operator()(elm, p, t) << vf->operator()(elm, p, t) 
    //  << "  =  " << A->operator()(elm, p, t) * vf->operator()(elm, p, t)  << endl;
    return A->operator()(elm, p, t) * vf->operator()(elm, p, t);
  }
 
  virtual concepts::Point<F, DIM> operator() (const concepts::ElementWithCell<G>& elm,
      const concepts::Real3d& p, const Real t = 0.0) const 
  {
    return A->operator()(elm, p, t) * vf->operator()(elm, p, t);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "ComposeFormulaMatVec(A=" << *A << ", vf="<< *vf << ")";
  }
public:
  RCP<const ElementFormula< concepts::Mapping<F, DIM>, G> > A;
  RCP<const ElementFormula< concepts::Point<F, DIM>, G> > vf;
};
 
// *************************************************** ComposeFormulaVecEntry **
 
template< class F, uint DIM, typename G = typename Realtype<F>::type>
class ComposeFormulaVecEntry : public concepts::ElementFormula< F, G > 
{
public:
  ComposeFormulaVecEntry( 
      RCP<const ElementFormula< concepts::Point<F, DIM> , G> > vf,
      int index
      ) 
    : vf(vf)
    , index(index)
  { }
 
  virtual ComposeFormulaVecEntry* clone() const {
    return new ComposeFormulaVecEntry(*this);
  }
 
  virtual F operator() (const concepts::ElementWithCell<G>& elm,
      const Real p,    const Real t = 0.0) const 
  {
    return vf->operator()(elm, p, t)[index];
  }
 
  virtual F operator() (const concepts::ElementWithCell<G>& elm, 
      const concepts::Real2d& p, const Real t = 0.0) const 
  {
    return vf->operator()(elm, p, t)[index];
  }
 
  virtual F operator() (const concepts::ElementWithCell<G>& elm,
      const concepts::Real3d& p, const Real t = 0.0) const 
  {
    return vf->operator()(elm, p, t)[index];
  }
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "ComposeFormulaVecEntry(i=" << index << ", vf="<< *vf << ")";
  }
public:
  RCP<const ElementFormula< concepts::Point<F, DIM>, G> > vf;
  int index;
};
 
// *********************************************** FormulaIncPlaneWaveSource **
 
class FormulaIncPlaneWaveSource : public ElementFormula<Cmplx> {
public:
  typedef std::map<int, double> MID;
 
  FormulaIncPlaneWaveSource(ElementFormulaContainer<Cmplx> coeff_a
                            , ElementFormulaContainer<Cmplx> coeff_b
/* RCP< PiecewiseConstFormula<Cmplx> > coeff_a */
/*       , RCP< PiecewiseConstFormula<Cmplx> > coeff_b */
                            , RCP< concepts::FormulaExpImag2D > u_inc )
    : coeff_a(coeff_a) // CAN BE IMPROVED
    , coeff_b(coeff_b)
    , u_inc(u_inc)
  { }
    
  virtual FormulaIncPlaneWaveSource* clone() const {
    return new FormulaIncPlaneWaveSource(*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>
  Cmplx operator() (const ElementWithCell< Real > &elm, const RealNd &p, Real2d px, 
      const Real t=0.0) const
  {
    Real2d k = u_inc->k;
    Cmplx k2 = k[0]*k[0] + k[1]*k[1]; // (can be complex)
    //Attribute atr = elm.cell().connector().attrib();
    Cmplx a = coeff_a(elm, p);
    Cmplx b = coeff_b(elm, p);
 
    if (a.real() != 1. || b.real() != 1.) {
      DEBUGL(pmlHolger_D, "a , b= " << a << ", " << b);
      DEBUGL(pmlHolger_D, "elemFrm value = " <<
             (b - a) * k2 * u_inc->operator()(elm.elemMap(p), t));
    } else {
      //DEBUGL(1, "a,b,b-a = " << a << b << b-a);
    }
    // NOTE: the term k2=omega^2 is not contained in coeff_b
    // (difference to the documentation)
    return (b - a) * k2 * u_inc->operator()(elm.elemMap(p), t);
    //return (b - a * k2) * u_inc->operator()(elm.elemMap(p), t);
  }
 
  static PiecewiseConstFormula<Cmplx> genTMCoeff(
          const MID& attToEps);
 
  static PiecewiseConstFormula<Cmplx> genTECoeff(
          const MID& attToEps);
 
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "FormulaIncPlaneWaveSource(" << coeff_a << ", " << coeff_b << ")";
  }
private:
  concepts::ElementFormulaContainer<Cmplx> coeff_a, coeff_b;
/*   RCP< concepts::ElementFormula<Cmplx> > coeff_a; */
/*   RCP< concepts::ElementFormula<Cmplx> > coeff_b; */
  RCP< concepts::FormulaExpImag2D > u_inc;
};
 
PiecewiseConstFormula<Cmplx> FormulaIncPlaneWaveSource::genTMCoeff(
        const MID& attrToEps)
{
    PiecewiseConstFormula<Cmplx> coeffs(1.);
    for(MID::const_iterator it = attrToEps.begin();
            it != attrToEps.end(); ++it) 
    {
        coeffs[it->first] =      it->second;
 
        if(it->second < 0) {
          double abs_fact = 1e-2;
          printf("Adding absorbing material (%f %f) for attr %d\n", 
              -it->second, abs_fact, it->first);
          coeffs[it->first] =      -it->second;
#if __cplusplus >= 201103L
          coeffs[it->first].imag(abs_fact);
#else
          //          imag(coeffs[it->first]) = abs_fact;
          coeffs[it->first].imag() = abs_fact;
#endif
 
        } else {
          printf("Adding non-absorbing material (%f %f) for attr %d\n", 
              it->second, 0., it->first);
        }
    }
    return coeffs;
}
 
PiecewiseConstFormula<Cmplx> FormulaIncPlaneWaveSource::genTECoeff(
        const MID& attrToEps)
{
    PiecewiseConstFormula<Cmplx> coeffs(1.);
    for(MID::const_iterator it = attrToEps.begin();
            it != attrToEps.end(); ++it) 
    {
        coeffs[it->first] =  1. / it->second;
 
        if(it->second < 0) {
          coeffs[it->first] =      -it->second;
#if __cplusplus >= 201103L
          coeffs[it->first].imag(1e-3);
#else
          //          imag(coeffs[it->first]) = 1e-3;
          coeffs[it->first].imag() = 1e-3;
#endif
        }
    }
    return coeffs;
}
 
// **************************************************** FormulaPMLPowerSigma **
 
template<typename F=Real>
class FormulaPMLPowerSigma : public concepts::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 concepts::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;
      //printf("sigma(%f), arg = %f, cen %f, off %f\n", p, arg, center, offset);
      if (arg < 0)
        return F(0);
      return sigma0 * powi(arg , power);
    }
 
    virtual F operator() (const concepts::Real2d& p, const Real t = 0.0) const {
      return (*this)(p[0], t);
    }
 
    virtual F operator() (const concepts::Real3d& p, const Real t = 0.0) const {
      return (*this)(p[0], t);
    }
    template<typename Real>
      static inline Real powi(Real x, int power) {
#if 1
        Real res(1);
        for(int i = 0; i < power; ++i) {
          res *= x;
        }
        return res;
#else
        return pow(x, power);
#endif
      }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << "FormulaPMLPowerSigma(" << sigma0 << " * ( |x - (" << center 
        << ")| - " << offset << ")_+^" << power << ")))";
    }
  private:
    const Real offset;
    const Real center;
    const int power;
    const F sigma0;
};
 
// ************************************************** FormulaPMLPowerSigma2D **
 
template<typename F=Real>
class FormulaPMLPowerSigma2D : public concepts::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);
    }
 
    bool inPMLregion(const concepts::Real2d& p, const Real t = 0.0) 
    {
      Real arg = sqrt((p[0]-center[0])*(p[0]-center[0]) + 
          (p[1]-center[1])*(p[1]-center[1]))  -offset;
 
      return arg >= 0;
    }
 
    virtual F operator() (const concepts::Real2d& p, const Real t = 0.0) const {
      Real arg = sqrt((p[0]-center[0])*(p[0]-center[0]) + 
          (p[1]-center[1])*(p[1]-center[1])) - offset;
      if (arg<0)
        return F(0);
      return sigma0 * powi(arg,power);
    }
 
    virtual F operator() (const concepts::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 power) {
#if 1
        Real res(1);
        for (int i = 0; i < power; ++i) {
          res *= x;
        }
        return res;
#else
        return pow(x, power);
#endif
      }
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << "FormulaPMLPowerSigma2D(" << sigma0 << " * ( |x - (" << center[0] 
        << "," << center[1] << ")| - " << offset << ")_+^" << power << ")))";
    }
  private:
    const Real offset;
    const Real2d& center;
    const int power;
    const F sigma0;
};
 
// ************************************************* FormulaPMLPowerSigmaB2D **
 
// Sigma bar 2D for Ridia case, see Collino & Monk
 
template<typename F=Real>
class FormulaPMLPowerSigmaB2D : public concepts::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("FormulaPMLPowerSigma Bar 2D currently only supports 2D!"); assert(false);
    }
 
    bool inPMLregion(const concepts::Real2d& p, const Real t = 0.0) 
    {
      Real arg = sqrt((p[0]-center[0])*(p[0]-center[0]) + 
          (p[1]-center[1])*(p[1]-center[1]))  -offset;
 
      return arg >= 0;
    }
 
    virtual F operator() (const concepts::Real2d& p, const Real t = 0.0) const {
      Real rho = sqrt((p[0]-center[0])*(p[0]-center[0]) + (p[1]-center[1])*(p[1]-center[1]));
      Real arg = rho - offset;
      if (arg<0)
        return F(0);
      return sigma0 * powi( arg , power + 1) / (rho) / (power + 1);
    }
 
    virtual F operator() (const concepts::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 power) {
#if 1
        Real res(1);
        for (int i = 0; i < power; ++i) {
          res *= x;
        }
        return res;
#else
        return pow(x, power);
#endif
      }
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << "FormulaPMLPowerSigma2D(" << sigma0 << " * ( |x - (" << center[0] 
        << "," << center[1] << ")| - " << offset << ")_+^" << power << ")))";
    }
  private:
    const Real offset;
    const Real2d& center;
    const int power;
    const F sigma0;
};
 
// ********************************************************** FormulaPMLCart **
 
class FormulaPMLCart : public ElementFormula<Cmplx> {
  public:
    enum PMLMode { DXDX, DYDY, IDENT, DX, DY };
 
    FormulaPMLCart(RCP<const ElementFormula<Cmplx> > coeff_a
        , RCP<const ElementFormula<Cmplx> > coeff_b
        , RCP<concepts::Formula<Real> > sigma_x
        , RCP<concepts::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 {
      std::cout << "FormulaPMLCart::clone() called " << *this << std::endl;
      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 {
      //std::cout << "gamma_y(" << y << ") " << Cmplx(1, (*sigma_y)(y) / omega) << std::endl;
      
      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 << "FormulaPMLCart(coeffs: " 
        << *coeff_a << ", " << *coeff_b << ", sigmas:" 
        << *sigma_x << ", " << *sigma_y 
        << ")";
    }
  private:
    RCP<const concepts::ElementFormula<Cmplx> > coeff_a;
    RCP<const concepts::ElementFormula<Cmplx> > coeff_b;
    RCP<concepts::Formula<Real> > sigma_x;
    RCP<concepts::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<concepts::FormulaPMLPowerSigma<Real> > sigma_x
      , RCP<concepts::FormulaPMLPowerSigma<Real> > sigma_y
      , RCP<concepts::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 << "FormulaPMLBoxRestriction(" << formula << " :" << *sigma_x << ", " << sigma_y << ")))";
  }
private:
  RCP<concepts::FormulaPMLPowerSigma<Real> > sigma_x;
  RCP<concepts::FormulaPMLPowerSigma<Real> > sigma_y;
  RCP<concepts::ElementFormula<F> > formula;
};
 
// ******************************************************* FormulaPMLRadia **
 
class FormulaPMLRadia : public ElementFormula<Cmplx> {
 public:
  enum PMLMode { AD1, AD2, AS, IDENT};
  
  FormulaPMLRadia(const ElementFormulaContainer<Cmplx> coeff_a
                  , const ElementFormulaContainer<Cmplx> coeff_b
                  , RCP<concepts::Formula<Real> > sigma
                  , RCP<concepts::Formula<Real> > sigmaB
                  , PMLMode mode
                  , double omega
                  )
    : coeff_a(coeff_a)
    , coeff_b(coeff_b)
    , sigma(sigma)
    , sigmaB(sigmaB)
    , mode(mode)
    , omega(omega)
    {}
    
    virtual FormulaPMLRadia* clone() const {
      return new FormulaPMLRadia(*this);
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real3d &p
                              , const Real t=0.0) const
    {
      throw std::string("FormulaPMLRadia currently only supports 2D!");
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real p
                              , const Real t=0.0) const
    {
      throw std::string("FormulaPMLRadia currently only supports 2D!");
    }
    Cmplx gamma(Real2d &p) const{
      return Cmplx(1, (*sigma )(p) / omega);
    }
    Cmplx gammaB(Real2d &p) const{
      return Cmplx(1, (*sigmaB)(p) / omega);
    }
    
    Real c_c(Real2d &p) const{
      return (p[0]*p[0])/(p[0]*p[0]+p[1]*p[1]); //\f$\cos^2\theta\f$
    }
    
    Real s_s(Real2d &p) const{
      return (p[1]*p[1])/(p[0]*p[0]+p[1]*p[1]); //\f$\sin^2\theta\f$
    }
    
    Real s_c(Real2d &p) const{
      return (p[0]*p[1])/(p[0]*p[0]+p[1]*p[1]); //\f$\cos\theta*\sin\theta\f$
    }
    
    virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real2d &p, const Real t = 0.0) const
    {
      Real2d px = elm.elemMap(p);
      
      switch(mode) {
      case AD1: return (gammaB(px)/gamma(px)*c_c(px) + gamma(px)/gammaB(px)*s_s(px)) * (coeff_a)(elm, p);
      case AD2: return (gammaB(px)/gamma(px)*s_s(px) + gamma(px)/gammaB(px)*c_c(px)) * (coeff_a)(elm, p);
      case AS : return ((gammaB(px)/gamma(px) - gamma(px)/gammaB(px))*s_c(px)) * (coeff_a)(elm, p);
      case IDENT: return gamma(px) * gammaB(px) * (coeff_b)(elm, p);
      }
      throw std::string("FormulaPMLRadia has gone mad!!");
    }
    
 protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << "FormulaPMLRadia(coeffs: " << coeff_a << ", " << coeff_b << 
        ", sigma" << *sigma << ")";
    }
 private:
    const ElementFormulaContainer<Cmplx> coeff_a;
    const ElementFormulaContainer<Cmplx> coeff_b;
    RCP<concepts::Formula<Real> > sigma;
    RCP<concepts::Formula<Real> > sigmaB;
    PMLMode mode;
    double omega;
};
 
} // namespace concepts