models/Eddy2D_H.hh

Go to the documentation of this file.

#ifndef Eddy2D_hh
#define Eddy2D_hh

#include "operator/sparseMatrix.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"

namespace hp2D {

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

  // forward declaration
  class InputEddy2D_H;

  // ***************************************************** Eddy2D_H_Interior **

  class Eddy2D_H_Interior : public concepts::OutputOperator {
  public:
    Eddy2D_H_Interior(const Real Omega_i, const Real H0) :
      Omega_i_(Omega_i), H0_(H0) {
      conceptsAssert(Omega_i > 0, concepts::Assertion());
    }
    Eddy2D_H_Interior(const Eddy2D_H_Interior& i) : 
      Omega_i_(i.Omega_i_), H0_(i.H0_) {}
    const Real Omega_i() const { return Omega_i_; }
    const Real H0() const { return H0_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << "Eddy2D_H_Interior(Omega_i = " << Omega_i_ 
    << ", H0_i = " << H0_ << ")";
    }
  private:
    const Real Omega_i_;
    const Real H0_;
  };

  // ************************************************************** Eddy2D_H **

  class Eddy2D_H : public AdaptiveModel<Cmplx>, public concepts::MaxwellModel {
  public:
    enum solverType { SUPERLU = 0, BICGSTAB = 1};
    Eddy2D_H(concepts::EddyGeometry2D& geom, const concepts::Formula<Real>& H0,
       const concepts::Formula<Real2d>& curlH0,
       const concepts::Formula<Real>* divgradH0 = 0,
       Eddy2D_H_Interior* interior = 0, const uint geomRefAttrib = 100, 
       const Real omega = OMEGA50,
       const Real mu = MU0, enum solverType type = SUPERLU);
    Eddy2D_H(concepts::EddyGeometry2D& geom, const Real H0,
       Eddy2D_H_Interior* interior = 0, const uint geomRefAttrib = 100, 
       const Real omega = OMEGA50,
       const Real mu = MU0, enum solverType type = SUPERLU);
    Eddy2D_H(concepts::EddyGeometry2D& geom, InputEddy2D_H& input,
       const uint geomRefAttrib = 100);
    virtual ~Eddy2D_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 const std::string mshAbbr_() { return geom_.meshAbbreviation(); }
  private:
    void constructor_(); 
    virtual hpFull& prebuild_() { return spc_->prebuild(); }
    virtual void solve_();
    void matrices_();
    void laplaceMatrix_();
    void identityMatrix_();
    bool connectedIdx_(uint& i);
    void systemMatrix_();
    void linearform_();

    std::auto_ptr<hpAdaptiveSpaceH1> spc_;
    concepts::EddyGeometry2D& geom_;
    std::auto_ptr<concepts::BoundaryConditions> bc_;
    std::auto_ptr<concepts::CellConditions> cc_;
    enum solverType type_;
    std::auto_ptr<concepts::Vector<Cmplx> > residual_;
    std::auto_ptr<Real> residualNorm_;
    std::auto_ptr<concepts::SparseMatrix<Real> > A_;
    std::auto_ptr<concepts::SparseMatrix<Real> > M_;
    std::auto_ptr<concepts::SparseMatrix<Cmplx> > S_;
    // load vector
    std::auto_ptr<concepts::Vector<Cmplx> > rhs_;
    concepts::PiecewiseFormulaFun<Real,Real> Sigma_Inv_;
    std::auto_ptr<const concepts::PiecewiseFormulaBase<Real> > H0_;
    std::auto_ptr<const concepts::PiecewiseFormulaBase<Real2d> > curlH0_;
    std::auto_ptr<const concepts::PiecewiseFormulaBase<Real> > divgradH0_;
     std::auto_ptr<Eddy2D_H_Interior> interior_;
    const Real omega_;
    const Real mu_;
    std::auto_ptr<Real> dissipation_;
    std::auto_ptr<Real> magnEnergy_;
    std::auto_ptr<const concepts::ElementFunction<Cmplx> > fun_;
    double solvetime_, matrixtime_, spacetime_;
  };

  // ****************************************************** InputEddy2D_H **

  class InputEddy2D_H : public concepts::InputParameter {
  public:
    InputEddy2D_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_; }
    const concepts::PiecewiseFormulaBase<Real>* H0() const { return H0_; }
    const concepts::PiecewiseFormulaBase<Real2d>* curlH0() const {
      return curlH0_; }
    const concepts::PiecewiseFormulaBase<Real>* divgradH0() const { 
      return divgradH0_; }
    bool interior() const { return false; }
    enum Eddy2D_H::solverType type() const { return type_; }
    bool solving() const { return solving_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    concepts::Sequence<Real> omega_;
    concepts::PiecewiseFormulaBase<Real>* H0_;
    concepts::PiecewiseFormulaBase<Real2d>* curlH0_;
    concepts::PiecewiseFormulaBase<Real>* divgradH0_;
    bool defaultOmega_;
    enum Eddy2D_H::solverType type_;
    bool solving_;
  };

} // namespace hp2D

namespace concepts {

  // *************************************************** HField_CircularCoil **

  /* Formula for the magnetic field created by a circular coil

      @author  Kersten Schmidt, 2005
   */
  class HField_CircularCoil : public Formula<Real> {
  public:
    HField_CircularCoil(const Real R1, const Real R2, const Real h0 = 1.0) :
      R1_(R1), R2_(R2), h0_(h0) {
      conceptsAssert(R2 >= R1, Assertion());
      conceptsAssert(R1 >= 0, Assertion());
    }
    virtual Real operator() (const Real p, const Real t = 0.0) const {
      const Real r = std::abs(p);
      if (r >= R2_) return 0.0;
      if (r <= R1_) return h0_;
      return (R2_ - r)/(R2_ - R1_) * h0_;
    }
    virtual Real operator() (const Real2d& p, const Real t = 0.0) const {
      return (*this)(p.l2());
    }
    virtual Real operator() (const Real3d& p, const Real t = 0.0) const {
      return (*this)(p.l2());
    }
    virtual HField_CircularCoil* clone() const {
      return new HField_CircularCoil(R1_, R2_, h0_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      os << "HField_CircularCoil(r <= " << R1_ << " : " << h0_ << ", ";
      if (R2_ > R1_)
  os << R1_ << " < r <= " << R2_ << " : " << R2_/(R2_-R1_)*h0_ << "-("
     << h0_/(R2_-R1_) << ")*r, ";
      return os << "r > " << R2_ << " : 0.0)";
    }
  private:
    const Real R1_, R2_;
    const Real h0_;
  };

  // *********************************************** CurlHField_CircularCoil **

  /* Formula for the curl of the magnetic field created by a circular coil

      @author  Kersten Schmidt, 2006
   */
  class CurlHField_CircularCoil : public Formula<Real2d> {
  public:
    CurlHField_CircularCoil(const Real R1, const Real R2, const Real h0 = 1.0)
      : R1_(R1), R2_(R2), h0_(h0) {
      conceptsAssert(R2 >= R1, Assertion());
      conceptsAssert(R1 >= 0, Assertion());
    }
    virtual Real2d operator() (const Real p, const Real t = 0.0) const {
      return (*this)(Real2d(p,0));
    }
    virtual Real2d operator() (const Real2d& p, const Real t = 0.0) const {
      const Real r = p.l2();
      if (r >= R2_) return 0.0;
      if (r <= R1_) return 0.0;
      return Real2d(p[1], -p[0]) * (-h0_ /(R2_ - R1_) / r);
    }
    virtual Real2d operator() (const Real3d& p, const Real t = 0.0) const {
      return (*this)(Real2d(p));
    }
    virtual CurlHField_CircularCoil* clone() const {
      return new CurlHField_CircularCoil(R1_, R2_, h0_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      os << "CurlHField_CircularCoil(";
      if (R2_ > R1_)      
  os << "r <= " << R1_ << " : " << 0.0 << ", "
     << R1_ << " < r <= " << R2_ << " : (y,-x)^T/|r| * " 
     << h0_ / (R2_-R1_) << ", r > " << R2_ << " : 0.0";
      else os << "0.0";
      return os << ")";
    }
  private:
    const Real R1_, R2_;
    const Real h0_;
  };

  // ******************************************** DivGradHField_CircularCoil **

  /* Formula for the curl of the magnetic field created by a circular coil

      @author  Kersten Schmidt, 2006
   */
  class DivGradHField_CircularCoil : public Formula<Real> {
  public:
    DivGradHField_CircularCoil(const Real R1, const Real R2,
             const Real h0 = 1.0)
      : R1_(R1), R2_(R2), h0_(h0) {
      conceptsAssert(R2 >= R1, Assertion());
      conceptsAssert(R1 >= 0, Assertion());
    }
    virtual Real operator() (const Real p, const Real t = 0.0) const {
      return (*this)(Real2d(p,0));
    }
    virtual Real operator() (const Real2d& p, const Real t = 0.0) const {
      const Real r = p.l2();
      if (r >= R2_) return 0.0;
      if (r <= R1_) return 0.0;
      return -h0_/(R2_ - R1_)/r;
    }
    virtual Real operator() (const Real3d& p, const Real t = 0.0) const {
      return (*this)(Real2d(p));
    }
    virtual DivGradHField_CircularCoil* clone() const {
      return new DivGradHField_CircularCoil(R1_, R2_, h0_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      os << "DivGradHField_CircularCoil(";
      if (R2_ > R1_)      
  os << "r <= " << R1_ << " : " << 0.0 << ", "
     << R1_ << " < r <= " << R2_ << " : " << -h0_/(R2_ - R1_) << "/r"
     << ", r > " << R2_ << " : 0.0";
      else os << "0.0";
      return os << ")";
    }
  private:
    const Real R1_, R2_;
    const Real h0_;
  };


} // namespace concepts

#endif // Eddy2D_hh