models/Maxwell2D_H.hh

Go to the documentation of this file.

#ifndef Maxwell2D_H_hh
#define Maxwell2D_H_hh

#include "operator/sparseMatrix.hh"
#include "operator/domainDecomp.hh"
#include "hp2D/hpAdaptiveSpaceH1.hh"
#include "formula/boundary.hh"
#include "models/adaptiveModels.hh"
#include "models/maxwell.hh"
#include "models/maxwellConstants.hh"
#include "models/Eddy2D_geometries.hh"
#include "models/Maxwell2D_H_eField.hh"
#include "space/domainDecomp.hh"

namespace hp2D {

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

  // forward declaration
  class InputMaxwell2D_H;

  // ****************************************************** Maxwell2D_H_Base **

  class Maxwell2D_H_Base : public AdaptiveModel<Cmplx>,
         public concepts::MaxwellModel,
         public concepts::MaxwellBoundary {
    friend class concepts::ModelControl<Maxwell2D_H_Base>;
  public:
    enum solverType { SUPERLU = 0, SUPERLU2 = 1, BICGSTAB = 2, BICGSTAB2 = 3,
          BICGSTABSUPERLU = 4};
    Maxwell2D_H_Base(concepts::EddyGeometry2D& geom,
         enum boundaryType bType = PMC,
         enum solverType type = SUPERLU, bool diagPrecond = true,
         bool afterIter = false,
         const Real eps = EPS0, const Real omega = OMEGA50,
         const Real mu = MU0, const uint geomRefAttrib = 100);
    Maxwell2D_H_Base(concepts::EddyGeometry2D& geom, InputMaxwell2D_H& input,
         const uint geomRefAttrib = 100);
    virtual ~Maxwell2D_H_Base() {}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    virtual const std::string mshAbbr_() { return geom_.meshAbbreviation(); }
    concepts::SparseMatrix<Cmplx>*
    laplaceMatrix_(concepts::Space<Real>& spc,
       concepts::SparseMatrix<Cmplx>* S);
    concepts::SparseMatrix<Real>*
    identityMatrix_(concepts::Space<Real>& spc,
        concepts::SparseMatrix<Cmplx>* S);
    void linearform_();

    concepts::EddyGeometry2D& geom_;
    std::auto_ptr<concepts::BoundaryConditions> bc_;
    enum solverType type_;
    bool diagPrecond_;
    bool afterIter_, statusAfterIter_;
    std::auto_ptr<concepts::Vector<Cmplx> > residual_;
    std::auto_ptr<Real> residualNorm_;
    // load vector
    std::auto_ptr<concepts::Vector<Cmplx> > rhs_;
    concepts::PiecewiseFormulaFun<Cmplx, Real> iOmegaEps_plus_Sigma_Inv_;
    Real eps_;
    Real omega_;
    const Real mu_;
    std::auto_ptr<Real> dissipation_;
    std::auto_ptr<Real> magnEnergy_;
    double solvetime_, matrixtime_, spacetime_;
    uint iterations_;
  private:
    void constructor_(); 
  };

  // *********************************************************** Maxwell2D_H **
  
  class Maxwell2D_H : public Maxwell2D_H_Base {
    friend class concepts::ModelControl<Maxwell2D_H>;
  public:
    Maxwell2D_H(concepts::EddyGeometry2D& geom,
    enum concepts::MaxwellBoundary::boundaryType bType = PMC,
    enum solverType type = SUPERLU, bool diagPrecond = true,
    bool afterIter = false,
    const Real eps = EPS0, const Real omega = OMEGA50,
    const Real mu = MU0, const uint geomRefAttrib = 100);
    Maxwell2D_H(concepts::EddyGeometry2D& geom, InputMaxwell2D_H& input,
    const uint geomRefAttrib = 100);
    virtual ~Maxwell2D_H() {}
    virtual hpAdaptiveSpaceH1& space() const;
    virtual Real dissipation();
    virtual Real magnEnergy();
    concepts::ElementFormula<Cmplx>* hField();
    concepts::ElementFormula<concepts::Cmplx2d>* eField();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    virtual hpFull& prebuild_() { return spc_->prebuild(); }
  private:
    std::auto_ptr<hpAdaptiveSpaceH1> spc_;
    std::auto_ptr<concepts::SparseMatrix<Cmplx> > A_, S_;
    std::auto_ptr<concepts::SparseMatrix<Real> > M_;

    void constructSpace_();
    virtual void solve_();
    void matrices_();
  };

  // ******************************************************** Maxwell2D_H_DD **
  
  class Maxwell2D_H_DD : public Maxwell2D_H_Base {
    friend class concepts::ModelControl<Maxwell2D_H_DD>;
  public:
    Maxwell2D_H_DD(concepts::EddyGeometry2D& geom,
       enum concepts::MaxwellBoundary::boundaryType bType,
       enum solverType type = SUPERLU, bool diagPrecond = true,
       bool afterIter = false,
       const Real eps = EPS0, const Real omega = OMEGA50,
       const Real mu = MU0, const uint geomRefAttrib = 100,
       std::string domains = "(2)");
    Maxwell2D_H_DD(concepts::EddyGeometry2D& geom, InputMaxwell2D_H& input,
       const uint geomRefAttrib = 100,
       std::string domains = "(2)");
    virtual ~Maxwell2D_H_DD();

    virtual concepts::DomainDecomp<hpAdaptiveSpaceH1>& space() const { 
      return *spc_;
    }
    Real dissipation();
    Real magnEnergy();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    virtual hpFull& prebuild_() { return prebuilds_; }
  private:
    hp2D::hpFull prebuilds_;
    std::auto_ptr<concepts::DomainDecomp<hp2D::hpAdaptiveSpaceH1> > spc_;
    concepts::Sequence<concepts::SparseMatrix<Cmplx>*> A_, S_;
    concepts::Sequence<concepts::SparseMatrix<Real>*> M_;

    void constructSpace_(const std::string& domainStr);

    virtual void solve_();
    void matrices_();
    void laplaceMatrix_();
    void identityMatrix_();
    

    template<class F>
    void clear_(concepts::Sequence<F*>& seq);
  };

  template<class F>
  void Maxwell2D_H_DD::clear_(concepts::Sequence<F*>& seq) {
    for (typename concepts::Sequence<F*>::iterator i = seq.begin(); 
   i != seq.end(); ++i)
      delete *i;
    seq.clear();
  }

  // ****************************************************** InputMaxwell2D_H **

  class InputMaxwell2D_H : public concepts::InputParameter {
  public:
    InputMaxwell2D_H(concepts::InOutParameters& input);
    virtual std::ostream& letters(std::ostream& os) const;
    virtual std::ostream& arguments(std::ostream& os) const;
    virtual std::ostream& description(std::ostream& os) const;
    virtual int input(int opt, const char* optarg);
    const concepts::Sequence<Real>& omega() const { return omega_; }
    enum Maxwell2D_H::solverType type() const { return type_; }
    enum concepts::MaxwellBoundary::boundaryType boundary() const { 
      return bType_;
    }
    bool solving() const { return solving_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    enum concepts::MaxwellBoundary::boundaryType bType_;
    concepts::Sequence<Real> omega_;
    bool defaultOmega_;
    enum Maxwell2D_H::solverType type_;
    bool solving_;
  };

} // namespace hp2D

namespace concepts {

  // ********************************** ModelControl<hp2D::Maxwell2D_H_Base> **

  template<>
  class ModelControl<hp2D::Maxwell2D_H_Base> : 
    public ModelControlBase<Model<Cmplx> > {
  public:
    ModelControl(hp2D::Maxwell2D_H_Base& model);
    virtual ~ModelControl() {}

    virtual hp2D::Maxwell2D_H_Base& model() {
      return static_cast<hp2D::Maxwell2D_H_Base&>(model_);
    }
    virtual const hp2D::Maxwell2D_H_Base& model() const {
      return static_cast<const hp2D::Maxwell2D_H_Base&>(model_);
    }
    virtual void matrices() = 0;
    void solverType(enum hp2D::Maxwell2D_H_Base::solverType type,
        bool diagPrecond);
    virtual Real solve() = 0;
    void storeMatricesToMatlab(const std::string matrixFile) const;
    void setOmega(const Real omega);
    void setEpsilon(const Real epsilon);
    void graphicsOut(const uint numPoints, const std::string path = "");
    hp2D::hpFull& prebuild() { return model().prebuild_(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // *************************************** ModelControl<hp2D::Maxwell2D_H> **

  template<>
  class ModelControl<hp2D::Maxwell2D_H> :
    public ModelControl<hp2D::Maxwell2D_H_Base> {
  public:
    ModelControl(hp2D::Maxwell2D_H& model);
    virtual ~ModelControl() {}
    virtual hp2D::Maxwell2D_H& model() {
      return static_cast<hp2D::Maxwell2D_H&>(model_);
    }
    virtual const hp2D::Maxwell2D_H& model() const {
      return static_cast<const hp2D::Maxwell2D_H&>(model_);
    }
    virtual void matrices();
    virtual Real solve();
    void storeMatricesToMatlab(const std::string matrixFile) const;
    void setOmega(const Real omega, bool assembleMatrices);
    void setEpsilon(const Real epsilon, bool assembleMatrices);
    void graphicsOut(const uint numPoints, const std::string path = "");
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // *********************************** ModelControl<hp2D::Maxwell2D_H_DD> **

  template<>
  class ModelControl<hp2D::Maxwell2D_H_DD> :
    public ModelControl<hp2D::Maxwell2D_H_Base> {
  public:
    ModelControl(hp2D::Maxwell2D_H_DD& model);
    virtual ~ModelControl() {}
    virtual hp2D::Maxwell2D_H_DD& model() {
      return static_cast<hp2D::Maxwell2D_H_DD&>(model_);
    }
    virtual const hp2D::Maxwell2D_H_DD& model() const {
      return static_cast<const hp2D::Maxwell2D_H_DD&>(model_);
    }
    virtual void matrices();
    virtual Real solve();
    void storeMatricesToMatlab(const std::string matrixFile) const;
    void setOmega(const Real omega, bool assembleMatrices);
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

} // namespace concepts



#endif // Maxwell2D_H_hh