app-bholger/gfem/gfemDiffraction.h

Go to the documentation of this file.

#pragma once

#include <string>

#include "app-bholger/pml/pml_formula.h"
#include "basics.hh"
#include "toolbox.hh"
#include "formula.hh"
#include "geometry.hh"
#include "space.hh"
#include "graphics.hh"
#include "function.hh"
#include "operator.hh"
#include "linearFEM.hh"
#include "linDG1D.hh"
#include "hp2D.hh"
#include "gfemIdentity.h"
#include "gfemLaplace.h"
#include "hp2D/formula.hh"
#include <string>
#include "gfemAdvectionIgnore.h" 
#include "gfemIdentityIgnore.h" 
#include "hp1D/bilinearForm.hh"

namespace concepts {
namespace gfem {

typedef std::map<int, double> MID;

class GfemDiffraction {
public:
  enum EquationType {TRANSVERSAL_ELECTRIC, TRANSVERSAL_MAGNETIC};

  GfemDiffraction(const SpaceOnCells<Real>& space, 
          EquationType type, 
          const MID& attrToEpsilon, 
          double omega,
          const BoundaryConditions& bc,
          Set<uint> nonRefl_bdAttr
          );

  static PiecewiseConstFormula<Cmplx> genTMCoeff(
          const MID& attToEps);

  static PiecewiseConstFormula<Cmplx> genTECoeff(
          const MID& attToEps);

  void setRhs(const std::string& rhs_r, const std::string& rhs_i);
        
  void adjustFloquetT(int dir, Cmplx k, int pml_cellAttrib 
      , RCP<const ElementFormula<Cmplx> >    uinc_F
      , RCP<const ElementFormula<Cmplx2d> > Duinc_F
      , concepts::Set<uint> pml_innerAttrib
      , Real2d center
      );
  
  void setNonRefl(Set<uint> nonRefl_bdAttr, 
      std::string u_inc_r, std::string u_inc_i,  
      std::string Du_inc_r, std::string Du_inc_i  
      );
 
  static const bool CONTRIB_RHS = true;
  static const bool NOCONTRIB_RHS = false;

  void set0thSommerfeld(Set<uint> nonRefl_bdAttr
      , RCP< concepts::ElementFormula< Cmplx > > rhs_F
      , RCP< concepts::ElementFormula< Cmplx > > dn_rhs_F
      , bool contribRHS = CONTRIB_RHS 
      );
 
  void setPMLBC(Real2d center, Real2d dist, Real strength, Real powercoeff
                , Real omega
                , RCP<const ElementFormula<Cmplx> > uinc_F
                , RCP<const ElementFormula<Cmplx2d> > Duinc_F
                , int pml_cellAttrib
                , concepts::Set<uint> pml_innerAttrib
                , concepts::Set<uint> pml_outerAttrib
                , RCP<const ElementFormula<Cmplx> > b2_0_F
      );

  // TODO: extend for the case where a(x) != const
  void setPMLBC_sc(
      Real2d center, Real2d dist, Real strength, Real powercoeff
      , Real omega
      , RCP<const ElementFormula<Cmplx> > uinc_F
      , RCP<const ElementFormula<Cmplx2d> > Duinc_F
      , RCP<const ElementFormula<Cmplx> > DaDuinc_F
      , concepts::Set<uint> innerAttrib
      , concepts::Set<uint> pml_outerAttrib
      , RCP<const ElementFormula<Cmplx> > b2_0_F
      );

  void setDtN( Set<uint> nonRefl_bdAttrTop
      , Set<uint> nonRefl_bdAttrBottom
      , Set<uint> nonRefl_bdAttr
      , std::string u_inc_r, std::string u_inc_i
      , std::string Du_inc_r, std::string Du_inc_i
      , int approx_oder
      , double min_y = 0
      );

  void solve(Vector<Cmplx>& sol );

  void storeBFeig(std::string path, bool verbose=false);


  // not yet supported
  //void setKPer(Cmplx k0, Cmplx k1);

  void setKPer(Cmplx k0);
  
  void compressSystemM() {
    system_M->compress();
  }

 #if 0

  void solveQR(Vector<Cmplx>& sol );
#endif
  
  void storeSystemM(const std::string& filename, Vector<Cmplx>* sol = NULL);

  void setOmega(double omega);

  void setB0(double b0) {
    this->b0 = b0;
  }

  double getb0() {
    return b0;
  }

  const hp2D::TraceSpace& getTraceSpace() { 
    return trace_space;
  }

  void assemble();
  
  void regularizeL2Penalty(double eps);

private:
  const SpaceOnCells<Real>& space;
  hp2D::TraceSpace trace_space;
  BoundaryConditions bc;

  double omega;
  double b0; 
  Cmplx kPer; 

public:
  static PiecewiseConstFormula<Cmplx> zero_PW;
  static PiecewiseConstFormula<Cmplx> one_PW;

  PiecewiseConstFormula<Cmplx> coeffId;
  PiecewiseConstFormula<Cmplx> coeffLap;
  RCP< const ElementFormula<Cmplx> > pmlcoeff_ident;
  RCP< const ElementFormula<Cmplx> > pmlLapCoeff_ident;
  RCP< const ElementFormula<Cmplx> > pmlLapCoeff_advectx;
  RCP< const ElementFormula<Cmplx> > pmlLapCoeff_advecty;
  RCP< concepts::MatrixElementFormula<Cmplx, 2> > coeffLapMat;

  RCP< SparseMatrix<Cmplx> > id_trace_M;

private:
  GfemLaplace lap_BF;
  GfemIdentity id_mass_BF;
  //GfemIdentity id_trace_BF;
  hp1D::Identity< Cmplx > id_trace_BF;
  RCP< GfemAdvectionIgnore > advect; // i * k0 * int(a u \conj v)
  RCP< GfemIdentityIgnore > id_mass_k2;

  Vector<Cmplx> rhs;
  RCP< SparseMatrix<Cmplx> > lap_M;
  RCP< SparseMatrix<Cmplx> > lapPML_M;
  RCP< SparseMatrix<Cmplx> > id_mass_M;
  RCP< SparseMatrix<Cmplx> > id_mass_k2_M;
  RCP< SparseMatrix<Cmplx> > advect_M;

  RCP< SparseMatrix<Cmplx> > system_M;

  RCP< FormulaPMLPowerSigma<Real> > sigma_x;
  RCP< FormulaPMLPowerSigma<Real> > sigma_y;
};

} // namespace
} // namespace