pml_edge.hh

#ifndef waveprop_pml_edge_hh
#define waveprop_pml_edge_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"
 
namespace concepts {
 
// ******************************************************** FormulaPMLEdgeRadia **
 
class FormulaPMLEdgeRadia: public ElementFormula<Cmplx> {
public:
  enum class Mode
    : int {
      A11,
    A22,
    A12,
    A21,
    M11,
    M22,
    M12,
    M21,
    N11,
    N22,
    N12,
    N21,
    IDENT,
    ROTROT
  };
 
  FormulaPMLEdgeRadia();
  FormulaPMLEdgeRadia(RCP<Formula<Real> > sigma, RCP<Formula<Real> > sigmaB,
      Mode mode, const Real2d& center = Real2d(0, 0));
 
  virtual FormulaPMLEdgeRadia* clone() const {
    return new FormulaPMLEdgeRadia(*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;
 
  RCP<Formula<Real> > sigma() const {
    return sigma_;
  }
  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_;
  Mode mode_;
  Real2d center_;
 
  Cmplx gamma_(Real2d &p) const {
    return Cmplx(1, (*sigma_)(p));
  }
  Cmplx gammaB_(Real2d &p) const {
    return Cmplx(1, (*sigmaB_)(p));
  }
};
 
// ******************************************************* RadialPML_2D **
 
class RadialPML_2D: public OutputOperator {
public:
  RadialPML_2D(const Real offset, const int power = 2,
      const Real sigma0 = 5.0, const Real2d& center = Real2d(0, 0));
 
  ElementFormulaContainer<Cmplx> getIdentityFormula() const {
    return *_u_v_formula_;
  }
 
  ElementFormulaContainer<Cmplx> getRotRotFormula() const {
    return *_rot_U_rot_V_formula_;
  }
 
  ElementFormulaContainer<MapCmplx2d> getVectorIdentityFormula() const {
    return *_U_V_formula_;
  }
 
  ElementFormulaContainer<MapCmplx2d> getLaplaceFormula() const {
    return *_grad_u_grad_v_formula_;
  }
 
  ElementFormulaContainer<MapCmplx2d> getGraduvFormula() const {
    return *_grad_u_V_formula_;
  }
 
protected:
  virtual std::ostream&
  info(std::ostream& os) const;
private:
  std::unique_ptr<FormulaPMLEdgeRadia> _u_v_formula_;
  std::unique_ptr<FormulaPMLEdgeRadia> _rot_U_rot_V_formula_;
  std::unique_ptr<MatrixElementFormula<Cmplx, 2> > _grad_u_grad_v_formula_;
  std::unique_ptr<MatrixElementFormula<Cmplx, 2> > _U_V_formula_;
  std::unique_ptr<MatrixElementFormula<Cmplx, 2> > _grad_u_V_formula_;
};
 
} // namespace concepts
 
#endif // waveprop_pml_edge_hh