buergDoerfler.hh

#ifndef buerg_hh
#define buerg_hh
 
//#include "estimator.hh"
#include "meshGraph.hh"
 
#include "refinement.hh"
 
//#include "basics/vectorsMatrices.hh"
 
#include "graphics.hh"
 
#include "formula/frmE_component.hh"
 
 
 
 
namespace estimator{
 
 
 
//LocalProjection type can specified with this pattern :
  // H4    : bisection of core cell and anisotrop h refinement of neighbour cells
  // P1Iso : p_K+1 on core cell and isotrop P_{K'}+1 refinement on neighbours {K'}
  // P2Iso : p_K+2 on core cell and isotrop P_{K'}+2 refinement on neighbours {K'}
  enum pattern {NONE, H4, P1ISO, P2ISO, P1ANISO, P2ANISO, H2V, H2H};
 
 
//class representing a local for a given quad
class LocalMesh2 : public concepts::Mesh2{
 
  public:
  LocalMesh2(const concepts::MeshGraph2& meshGraph, uint K );
 
     ~LocalMesh2();
 
     //inline concepts::MutableMesh2& get(){ return mMesh_; }
     inline uint boundary() const { return bndAttrib_; }
     //returns the sources, ie. the original pointers of hp2D::Elements used
     inline const concepts::HashMap<const hp2D::Quad<Real>* >& sources() const {return elms_;}
    // inline const concepts::HashMap<concepts::Quad*>& connectors() const { return qMap_; }
     inline uint ncell() const { return cMap_.size(); }
     inline concepts::Scan2* scan() { return new S(cMap_); }
     //inline const concepts::Scan2* scan() const { return new S(cMap_); }
     virtual std::ostream& info(std::ostream& os) const{return os << "LocalMesh[]";}
     //Scanner of the LocalMesh
     private:
 
        class S : public concepts::Scan<concepts::Cell2> {
          concepts::HashMap<concepts::Quad2d*>::iterator iter_;
          concepts::HashMap<concepts::Quad2d*>& cells_;
        public:
          inline S(concepts::HashMap<concepts::Quad2d*>& cells):
          cells_(cells), iter_(cells.begin()){}
          inline S(const S &scan) : iter_(scan.iter_), cells_(scan.cells_) {};
 
          inline bool eos() const { return iter_ == cells_.end(); }
          inline concepts::Cell2& operator++(int) { return *((iter_++)->second); }
 
          inline concepts::Scan2* clone() const { return new S(*this); }
        };
 
    //map from vtx key to vertex
    concepts::HashMap<concepts::Vertex*> vMap_;
    concepts::HashMap<concepts::Edge*> eMap_;
    concepts::HashMap<concepts::Quad*> qMap_;
    concepts::HashMap<concepts::Quad2d*> cMap_;
 
    //new unique boundary (marked) attrbute for this local mesh
    uint bndAttrib_;
 
    //collection of the original hp2D Elements, here i.e Quads
    //those collection is needed outside, i.e. for the local projections
    concepts::HashMap<const hp2D::Quad<Real>* > elms_;
  }; //local Mesh Class
 
 
//********************************************* Class LOcalProjectionBase  *******
 
template<class F, class H = typename concepts::Realtype<F>::type >
class LocalProjectionBase : public concepts::OutputOperator {
 
public :
 
  LocalProjectionBase();
  LocalProjectionBase(enum pattern pat);
 
protected :
 
  //pre extraced coefficient due to coarse space element of original space
  concepts::HashMap<concepts::Array<F> > coeffs_;
 
  //the pattern type which specifies the projection, i.e. the map of PointProjection
  enum pattern pat_;
 
  void addH4_(const concepts::HashMap<const hp2D::Quad<H>* >& source,
         LocalMesh2& locMsh,
         const concepts::Vector<F>& sol);
 
  void addP1ISO_(const concepts::HashMap<const hp2D::Quad<H>* >& source,
           LocalMesh2& locMsh,
           const concepts::Vector<F>& sol);
 
  struct PointProjection{
 
  public:
    //Further refinements may possible
    enum subType {NONE, div4, divH2 , div2V};
 
    //default for hashmap construct, is never used on real info adding
    PointProjection(){ i = 5;};
 
    //i the number of the local quad child
    //K key of the coarse (original space) element
    PointProjection(enum subType qs, uint i,  const hp2D::Quad<H>* elm);
 
    PointProjection& operator=(const PointProjection& pP);
 
    //transforms the point due to quadchild number (i , located from left to right anticlockwise
    void operator()(concepts::Real2d& p) const;
 
    //quad child number position
    uint i;
    //the subdevision type
    enum subType sub;
    uint K;
    const hp2D::Quad<H>* elm;
  };
 
 
   //maps from finer cell key to pointprojection data, i.e quad child number, father quad of original space
   // which ie. has T matrix, and subdevision info, used for the application operators for integration
   // this is used to avoid output operator
   __gnu_cxx::hash_map<uint, PointProjection> data_;
 
 
};
 
 
//********************************************* Class LocalProjections_S  *******
 
// Scalar returner point
 
template<class F, class H = typename concepts::Realtype<F>::type >
class LocalProjections_S : public LocalProjectionBase<F,H>, public concepts::ElementFormula<F, H >{
  //concepts::ElementFormula<concepts::Point<F,dim>,H>{ //for vectorvalued
 
public:
 
  LocalProjections_S(const enum pattern pat,
               const concepts::HashMap<const hp2D::Quad<H>* >& source,
             LocalMesh2& locMsh,
             const concepts::Vector<F>& sol,
             const concepts::ElementFunction<F,H>& f);
 
  LocalProjections_S(const LocalProjections_S<F,H>& frm);
 
  //TODO: This must be done correctly deletion wont work atm ... cannot double delete stuff
  virtual ~LocalProjections_S();
 
  virtual F operator() (const concepts::ElementWithCell<H>& elm, const Real p, const Real t = 0.0) const;
 
  virtual F operator() (const concepts::ElementWithCell<H>& elm, const concepts::Real2d& p,const Real t = 0.0) const;
 
  virtual F operator() (const concepts::ElementWithCell<H>& elm, const concepts::Real3d& p, const Real t = 0.0) const;
 
  virtual LocalProjections_S<F,H>* clone() const{ return new LocalProjections_S<F,H>(*this);}
 
 
 
 
protected:
  virtual std::ostream& info(std::ostream& os) const;
 
private:
  //type of projection, e.g. Value, Grad, Laplacian, ...
  std::auto_ptr<concepts::ElementFunction<F,H> > f_;
};
 
 
 
 
//********************************************* Class LocalProjections_P  *******
// localprojections returning points
 
//template<uint dim, typename F, typename G>
template<uint dim, class F, class H = typename concepts::Realtype<F>::type >
class LocalProjections_P : public LocalProjectionBase<F,H>, public concepts::ElementFormula<concepts::Point<F,dim>, H >{    //concepts::ElementFormula<concepts::Point<F,dim>,H>{ //for vectorvalued
 
public:
  LocalProjections_P(const enum pattern pat,
               const concepts::HashMap<const hp2D::Quad<H>* >& source,
             LocalMesh2& locMsh, //TODO: Evtl den space nehmen nicht local Mesh, damit sichergestellt wird dass nicht falsches mesh evaluiert wird da kinder dann gar nicht vorhanden
             //local mesh ist notwendig da sonst die Vater zu kind beziehung flöten geht
             const concepts::Vector<F>& sol,
             const concepts::ElementFunction<F,H>& f);
 
  LocalProjections_P(const LocalProjections_P<dim,F,H>& frm);
 
 
  //TODO: This must be done correctly deletion wont work atm ... cannot double delete stuff
  virtual ~LocalProjections_P();
 
 
  virtual concepts::Point<F,dim> operator() (const concepts::ElementWithCell<H>& elm, const Real p, const Real t = 0.0) const;
 
  virtual concepts::Point<F,dim> operator() (const concepts::ElementWithCell<H>& elm, const concepts::Real2d& p,const Real t = 0.0) const;
 
  virtual concepts::Point<F,dim> operator() (const concepts::ElementWithCell<H>& elm, const concepts::Real3d& p, const Real t = 0.0) const;
  virtual LocalProjections_P<dim,F,H>* clone() const{ return new LocalProjections_P<dim,F,H>(*this);}
 
protected:
  virtual std::ostream& info(std::ostream& os) const;
 
private:
  //type of projection, e.g. Value, Grad, Laplacian, ...
  std::auto_ptr<concepts::ElementFunction<F,H> > f_;
};
 
 
 
//********************************************* Class BuerGDoerfler  *******
 
template<class F = Real>
class BuergDoerfler : public Marking<F>{
 
public :
 
  BuergDoerfler(const concepts::SpaceOnCells<Real>& spc, const LocalEstimator<F>& estimator, const concepts::Vector<F>& sol);
 
  virtual ~BuergDoerfler();
 
  //MeshGraph2(const concepts::SpaceOnCells<Real>& spc, bool regular = true)
 
  //FIX ME: rethink template here, for sol this is ok, seems ok
  void addRhs_f(const concepts::ElementFormula<F>& f, const Real c = 1);
 
  //TODO: Implement this method to possible generalize burg doerfler
  //add information for neumann/robin boundary condition
  void addRhs_g(const concepts::ElementFormula<F>& g, const concepts::Set<uint>& set);
 
  //TODO: Implement this method to possible generalize burg doerfler
  //add information to build weak formulation i.e. the billinearforms
  void addLhs(const concepts::BilinearForm<F>&, const Real c = 1);
 
  void compute();
 
 
  //FIXME This is global just for testing ...
 
  //convergence indicator for h4 pattern, maps from Key to indicator
  concepts::HashMap<Real> indicator_H4_;
  //convergence indicator for p1iso pattern, maps from Key to indicator
  concepts::HashMap<Real> indicator_P1Iso_;
 
 
 
  concepts::HashMap<Real> sup_H4_;
  concepts::HashMap<Real> sup_P1Iso_;
 
  //dofs of corresponding refinement pattern spaces :
  concepts::HashMap<uint> dofs_H4_;
  concepts::HashMap<uint> dofs_P1Iso_;
 
protected:
  virtual std::ostream& info(std::ostream& os) const;
 
 
private :
  //all rhs f data added
  concepts::Sequence<const concepts::ElementFormula<F>* > rhs_f_;
  //and their building coefficients
  concepts::Sequence<Real> c_f_;
  //bool controls if rhs_f data is added
  bool rhs_f_b_;
 
 
 
 
  const concepts::SpaceOnCells<Real>& spc_;
  const concepts::Vector<F>& sol_;
 
 
  concepts::MeshGraph2 meshgraph_;
  const LocalEstimator<F>& estimator_;
 
 
  //convergence indicators, maps from cell key to its refinement patterns
  //corresponding solution of the local problem
 
//  //convergence indicators from h bisection refinement
//  concepts::HashMap<Real> ind_hbisect;
//
//  //convergence indicators from p increase plus 1
//  concepts::HashMap<Real> ind_p1;
 
 
  void compute_H4Iso_Ind_K(LocalMesh2& localMesh,
             uint K,
             Real n_K,
             const hp2D::hpFull& globPrebuild,
             concepts::BoundaryConditions& bc);
 
 
 
  void compute_P1Iso_Ind_K(LocalMesh2& localMesh,
               uint K,
               Real n_K,
               const hp2D::hpFull& globPrebuild,
               concepts::BoundaryConditions& bc);
 
  void compute_();
 
  hp2D::hpFull* hbisectPattern_(concepts::Mesh2& msh, uint K,const hp2D::hpFull& globPrebuild);
 
 
  //p(cell)  + 1 pattern on each  cell
  hp2D::hpFull* p1Pattern_(concepts::Mesh2& msh, uint K,const hp2D::hpFull& globPrebuild);
 
};
 
 
template<class F = Real>
class BuergRefinement : public estimator::HP_Refinement<F,2> {
 
 
 
 
    struct weightedIndictor : public concepts::OutputOperator{
 
      weightedIndictor(uint K, Real sup, enum pattern pat);
 
 
      bool operator<(const weightedIndictor& wI) const;
 
        bool operator==(const weightedIndictor& wI) const;
 
        //returns the value
        Real get() const{ return sup_2;}
 
        //key of corresponding
        uint K;
        //the supremum squared = indicator(pat) * n_K
        Real sup_2;
       // Real dof;
      //  Real ind;
        //pattern type
        enum pattern type_;
 
    protected:
        virtual std::ostream& info(std::ostream& os) const;
 
    };
 
    public:
    //friend class BuerGDoerfler<F>;
 
    BuergRefinement(Real theta = 0.35);
 
    virtual ~BuergRefinement();
 
    // builds the refinement strategys for given mark and estimator
    void buildRefinement(const Marking<F>& mark,  const LocalEstimator<F>& estimator);
 
  virtual std::ostream& info(std::ostream& os) const;
 
    private:
      //current cells of given space
      //concepts::HashMap<const concepts::Cell* > cells_;
 
      //possible adaptive bulk criteria parameter
      Real theta_;
 
 
 
};
 
 
 
 
 
 
 
 
}
 
 
 
 
 
#endif // buerg_hh