app-pkauf/PHC4.hh

Go to the documentation of this file.

#ifndef CONV_PHC4_HH_
#define CONV_PHC4_HH_

#include <fstream>
#include <memory>
#include <vector>
#include <complex>
#include "basics/outputMatlab.hh"
#include "basics/exceptions.hh"
#include "hp2D/spacePreBuilder.hh"
#include "geometry/meshImport2Dgeneral.hh"
#include "graphics/matlab.hh"
#include "basics/outputMatlab.hh"
#include "basics/exceptions.hh"
#include "hp2D/hpAdaptiveSpaceH1.hh"
#include "hp2D/bilinearForm.hh"
#include "operator/sparseMatrix.hh"
#include "operator/denseMatrix.hh"
#include "operator/superLU.hh"
#include "eigensolver/ARPACK.hh"
#include "space/formula.hh"
#include "toolbox/resourceMonitor.hh"
#include "Brillouin.hh"
#include "hp2D/aprioriRef2D.hh"
#include "toolbox/inputOutput.hh"

using concepts::Real;
using concepts::Cmplx;

namespace phc {

  // ****************************************************************** PHC4 **
  
  class PHC4 {
    public:
    /* Constructor: Imports the spaces, constructs the k-values (uint i)
       and the bilinearforms with dielectricity constants epsilon 
       (std::vector<Real> epsilon) on attributes 
       (std::vector<uint> attributes) and default epsilon (Real dflt).
       @param spc vector containing the adresses of the four spaces.
       @param i number of points on the boundary of the irreducible 
       Brillouin zone used in the computation.
       @param epsilon Must have the same length as attributes. Contains the 
       epsilon values belonging to the same index of the attribute vector.
       @param attributes see epsilon
       @param dflt The default value for the dielectric constant.*/
    PHC4(std::vector<hp2D::hpAdaptiveSpaceH1*> spc, Brillouin_Base& B, 
   uint i, std::vector<Real> epsilon, std::vector<uint> attributes,
   Real dflt = 1.0);
    PHC4(std::vector<hp2D::hpAdaptiveSpaceH1*> spc, Brillouin_Base& B, 
   concepts::Array<concepts::Real2d*> k_val, std::vector<Real> epsilon, 
   std::vector<uint> attributes, Real dflt = 1.0);
    ~PHC4();
    // rebuilding:
    void rebuild();
    hp2D::hpAdaptiveSpaceH1 get_spc(uint i);
    
    concepts::Array<concepts::Real2d*> get_k() { return k_; }
    // Constructing the bilinear forms with dielectric constants:
    void bilinearforms(std::vector<Real> epsilon, 
           std::vector<uint> attributes,
           Real dflt = 1.0);
    // Assembling the basic matrices:
    void assemble();
    
    /*get_16_matrices_k1_k2() saves the assembled matrices corresponding to
      the combinations of the 4 spaces as matlab files m_1, ... , m_16. 
      The k-vector array corresponding to the Brillouin zone iteration is saved
      to the file Brill_factors.m. Use the matlab routine Bandstructure.m to 
      compute the Bandstructure.*/
    
    void get_16_matrices_k1_k2_(std::string path = std::string(".")); 
    void get_16_matrices_k1_k2_TE(std::string path = std::string("."));
    void get_16_matrices_k1_k2_TM(std::string path = std::string("."));
    
    // Getting the multiplication factors out of k_x, k_y:
    void factor(Real k_x, Real k_y, concepts::Array<Real> v1,
    concepts::Array<Real> v2);

    /* The routine handling the Brillouin iterations:
       @param kmax Number of eigenvalues to compute (should not be too small)
       @param maxiter maximal number of Arnoldi iterations to be performed
       @param eps tolerance up to which Arnoldi should converge
       @param n saving the eigenfunctions to the nth eigenvalue (counted from 
       below) if EF_graph == true.
       @EF_graph if true the eigenfunctions are plotted
       The output of this routine are matlab files containing an eigenvalue 
       matrix (Brillouin Point x Eigenvalue), Eigenvalue in increasing order 
       at first Brillouin point. If EF_graph == true you get another set of 
       matlab files containing the basis functions, indexed with the 
       Brillouin point index..
    */
    void compute(uint kmax, uint maxiter, Real eps, uint n,
     bool EF_graph);
    /* If you want to compute the generalized eigenvalues with matlab, use 
       this routine. 
       Specify the Brillouin Point at which you want to look and indicate 
       the refinement step for the name of the output matrices. With the 
       boolian TE you can choose whether you want TE modes (TE_TM = 0) or 
       TM modes (TE_TM = 1).
       The output m-files with the matrices will be called 
       'Laplace_T*_"Brillouin_Point"_"Ref_Step"_.m' and
       'Identity_T*_"Brillouin_Point"_"Ref_Step"_.m', with * being 'E' or 'M'
       corresponding to TE_TM. Also a file with the times needed to assemble
       the matrices is written.
    */ 
    void give_matrices(uint Brillouin_Point, uint Ref_Step,  
                       std::string tag);
    void rcount() { refine_++; }
    /* If you want to compute the eigenvalues of a refined problem 
       use this routine:
       Same output as compute, but with additional parameter in the
       matlab file for number
       of recomputation step (this parameter is zero in compute(...).
    */
    void recompute(uint kmax, uint maxiter, Real eps, uint n, bool EF_graph);
    uint nDoF() { return nDoF_; }
    // returns pointer to  the eigenvalue matrix for TE / TM modes:
    concepts::SparseMatrix<Real>* get_ev_TE() { return ev_TE_matrix_; }
    concepts::SparseMatrix<Real>* get_ev_TM() { return ev_TM_matrix_; }
    void set_tag(std::string tag){ tag_ = tag; set_tag_ = true; }
    double get_time(){ return t_; }
    
  private:
    // The four spaces:
    hp2D::hpAdaptiveSpaceH1& spc1_;
    hp2D::hpAdaptiveSpaceH1& spc2_;
    hp2D::hpAdaptiveSpaceH1& spc3_;
    hp2D::hpAdaptiveSpaceH1& spc4_;
    std::vector<hp2D::hpAdaptiveSpaceH1*> spc_;
    hp2D::hpFull& prebuild_;
    uint refine_;
    uint nDoF_;
    Brillouin_Base& B_;
    concepts::Array<concepts::Real2d*> k_;
    // Bilinearforms for TE and TM modes:
    hp2D::Laplace<> *la_TE_, *la_TM_;
    hp2D::Identity<> *id_TE_, *id_TM_;
    // for the assembly routine:
    bool assembled_;
    bool rebuilt_;
    bool computed_;
    bool EF_graph_;
    // The different matrices without factors:
    std::vector<concepts::SparseMatrix<Real>* > mat_la_TE_;
    std::vector<concepts::SparseMatrix<Real>* > mat_la_TM_;
    std::vector<concepts::SparseMatrix<Real>* > mat_id_TE_;
    std::vector<concepts::SparseMatrix<Real>* > mat_id_TM_;
    // The array of multiplication factors:
    concepts::Array<Cmplx> factor_;
    // The single space factors:
    concepts::Array<Cmplx> single_factor_;
    // matrices of eigenvalues:
    concepts::SparseMatrix<Real>* ev_TE_matrix_;
    concepts::SparseMatrix<Real>* ev_TM_matrix_;
    // The residual starting vectors:
    concepts::Array<Cmplx>* resid_TE_;
    concepts::Array<Cmplx>* resid_TM_;
    // intermediate eigenvalue arrays:
    concepts::Array<concepts::Array<Cmplx>*>* ev_TE_;
    concepts::Array<concepts::Array<Cmplx>*>* ev_TM_;
    concepts::Vector<Cmplx>* EF_TE_;
    concepts::Vector<Cmplx>* EF_TM_;
    // a tag that is appended to e.g. ev_TE_"tag"_1.m
    std::string tag_;
    bool set_tag_;
    double t_;

    // Sorting the eigenvalues / eigenvectors:
    concepts::Array<uint> sort_(concepts::Array<Cmplx>* unsorted);

    void brillouinIteration_(uint kmax, uint maxiter, Real eps, uint n, 
           uint j, concepts::LiCoI<Real>* Id_OP);
  };

} // namespace phc

#endif /*CONV_PHC4_HH_*/