refinement.hh

#ifndef refinement_hh
#define refinement_hh
 
//#include "estimator.hh"
#include "estimator/marking.hh"
#include "space/hpMethod.hh"
#include "hp2D/hpAdaptiveSpace.hh"
#include "hp2D/quad.hh"
#include "hp2D/hpAdaptiveSpaceH1.hh"
//#include "lInfinityNorm.hh"
//#include "geometry/diameter.hh"
 
#if __GNUC__ == 2
#include <float.h>
#define EPS DBL_EPSILON
#define MAX_D DBL_MAX
#else
#include <limits>
#define MAX_D std::numeric_limits<double>::max()
#endif
 
using concepts::Real;
 
namespace estimator{
 
template<class F=Real>
class Refinement : public concepts::OutputOperator{
 
public:
  //apply routine that adjust the elements in the space,
  //e.g. elements that are marked for refinement
 
  virtual void apply(concepts::SpaceOnCells<F>& space) = 0;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const = 0;
 
};
 
 
//template<class F=Real>
//class H_Refinement: public concepts::Refinement<F>{
//public:
//  H_Refinement(const concepts::Marking& mark): mark_(mark){}
//
//private :
//    concepts::Marking& mark_;
//};
 
// Abstract Class for hp refinement strategies.
template<class F = Real, ushort dim = 2>
class HP_Refinement : public Refinement<F>{
public :
 
  virtual ~HP_Refinement(){};
 
  virtual void apply(concepts::SpaceOnCells<F>& space){
    hp2D::hpAdaptiveSpace<F>* spc = dynamic_cast<hp2D::hpAdaptiveSpace<F>*>(&space);
    if(spc){
      std::unique_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
      while (!sc->eos()) {
        concepts::ElementWithCell<Real>& element = (*sc)++;
        const hp2D::Quad<Real>* quad = dynamic_cast<const hp2D::Quad<Real>*> (&element);
        uint K = quad->support().key();
 
        typename concepts::HashMap<concepts::AdaptiveAdjustP<dim> >::const_iterator iter;
        //current element has a refinement rule
        if((iter = adjusts_.find(K)) != adjusts_.end())
          spc->adjust(*quad, iter->second);
      }
      spc->rebuild();
    }else{
      throw conceptsException(concepts::MissingFeature("Currently just hpAdaptiveSpace supported"));
    }
 
  }
 
protected :
  //map from cellKey to its adjust strategy
  concepts::HashMap<concepts::AdaptiveAdjustP<dim> > adjusts_;
 
};
 
template<class F = Real, ushort dim = 2>
class Prediction : public HP_Refinement<F,dim>{
public:
 
  //starting prediction assumption
  // if cell was not refined before :
  // H : h refinement is forced for first refinement
  // P : p refinement is forced for first refinement
  // after prediction computation is possible
  enum Mode {H = 0 , P = 1};
 
  Prediction(const concepts::SpaceOnCells<F>& space, enum Mode mode = P,
         const Real hdec = sqrt(4.0), const Real pdec = sqrt(0.4), const Real nRdec = 1.0):
         hdecay_(hdec), pdecay_(pdec), noRefdecay_(nRdec)
  {
    const hp2D::hpAdaptiveSpace<F>* spc = dynamic_cast<const hp2D::hpAdaptiveSpace<F>* >(&space);
    if(spc){
      //information about the polynomial degree
      hp2D::hpFull& hpfull = spc->prebuild();
      std::unique_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
        while (*sc) {
          concepts::ElementWithCell<Real>& elm = (*sc)++;
          //key of current element
          uint K = elm.cell().connector().key().key();
          //set default prediction strategy according to chosen mode
          pred_[K] = mode ? MAX_D : 0;
          //add cell information
          cells_[K] = &elm.cell();
 
          const concepts::Connector2* cntr =
              dynamic_cast<const concepts::Connector2*>(&elm.cell().connector());
          //the cntr must be 2d since it is out of hpAdaptiveSpace
          conceptsAssert(cntr, concepts::Assertion());
          //set maximal polynomial degree (i.e inner), if no inner exists
          // default is 1 if there is a vertex dof on this cell
          const ushort* pMax = hpfull.innerDof(*cntr);
          //no inner dof
          if(pMax == 0){
            //check if there is at least one dof, i.e a vertex dof
            //this has 4 vertexList searches (costs?)
            bool hasVDof = (hpfull.vtxDof(*cntr, 0) || hpfull.vtxDof(*cntr, 1) ||
                  hpfull.vtxDof(*cntr, 2) || hpfull.vtxDof(*cntr, 3));
            if(hasVDof){
              //set default to one
              p_[K] = concepts::Array<ushort>(dim);
              for(ushort i = 0 ; i < dim; ++i)
                p_[K][i] = 1;
            }else{
              //no dofs on vertices, e.g. :
              //
              // Cell x for h=1 only has hanging nodes as vertices and is therefore not supported
              //  __________
              // |__________|
              // |   |x|    |
              // |   |_|
              // |___|_|____|
              throw conceptsException(concepts::MissingFeature("Cell has no vertex dofs, cannot be handled atm"));
            }
          }else{
            //has inner Dof and therefore a polynomial degree >= 2
            p_[K] = concepts::Array<ushort>(dim);
            for(ushort i = 0 ; i < dim; ++i)
               p_[K][i] = pMax[i];
          }
 
        }
    }else{
      throw conceptsException(concepts::MissingFeature("Currently just hpAdaptiveSpace supported"));
    }
  }
 
 
  virtual ~Prediction(){}
 
 
// builds the refinement strategys for given mark and estimator
void buildRefinement(const Marking& mark,  const concepts::LocalEstimator<F>& estimator){
 
  //all marked Cells (its keys)
  const concepts::Set<uint>& marked = mark.marks();
 
 
  //check all cells if they are marked
  concepts::HashMap<const concepts::Cell* >::const_iterator cIter = cells_.begin();
  for( ; cIter != cells_.end(); ++cIter){
    //current cell
    uint K = cIter->first;
    //cell is marked for refinement
    if(marked.exist(K)){
      //DEBUGL(1, "on K = "<< K << " error = "<< estimator(K) << "  pred = "<< pred_[K]);
      if (estimator(K) > pred_[K]){
        //DEBUGL(1, "on K = "<< K << " error = "<< estimator(K) << "  pred = "<< pred_[K] << " apply h ref");
        this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
        //add cell to as href marked, later needed to compute child prediction
        //those are not avaiable atm
        hRefCurCells_.push_back(cIter->second);
        //only isotropic p refinement is allowed
        conceptsAssert(p_[K][0]==p_[K][1], concepts::Assertion());
        //cell will have childs later, now we save child info for it here
        pred_[K] = 0.5 * hdecay_ * std::pow(0.5,p_[K][0] ) * estimator(K);
        //also polynomial degree is setted on kids etc
 
        //DEBUGL(1, "set prediction for error = "<< estimator(K) << " mult = "<< 0.5 * hdecay_ * std::pow(0.5,p_[K][0] )
        //        << " so pred = " << pred_[K]);
      }
      //set for p refinement
      else{
        //DEBUGL(1, "on K = "<< K << " error = "<< estimator(K) << "  pred = "<< pred_[K] << " apply p ref");
        //add isotropic p refinement plus 1
        this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(0, 1);
        //update info of polynomial degree on this element
        p_[K][0] += 1;  p_[K][1] += 1;
        //DEBUGL(1, "old pred_[K] = "<< pred_[K] << "  pdecay_ = "<< pdecay_ << " new  :");
        pred_[K] = pdecay_ * estimator(K);
      }
 
      //DEBUGL(1, "prediction( K  = "<< K << " ) "<< pred_[K]);
 
    }else{
      //not marked update with no refinement decay
      pred_[K] = noRefdecay_ * pred_[K];
    }
  }
}
 
 
virtual void apply(concepts::SpaceOnCells<F>& space){
  //refine space
  HP_Refinement<F,dim>::apply(space);
  // clear all adjusts as they needed to be evaluated new for same elements
  this->adjusts_.erase(this->adjusts_.begin(), this->adjusts_.end());
 
  //update informations for predictions, i.e. for h refined elements, that after apply
  //have new children
  concepts::Sequence<const concepts::Cell* >::const_iterator cIter;
  for(cIter = hRefCurCells_.begin(); cIter != hRefCurCells_.end(); ++cIter){
 
    //current adjacent cell key
    uint K = (*cIter)->connector().key().key();
    const concepts::Cell* child = 0;
    //current child it must be 4, else not supported
    //since prediction must be adapted for halfsplitting, e.g. in 2D
    uint i = 0 ;
    for( ; (child = (*cIter)->child(i)) !=0 ; ++i){
      //key of builded child cell
      uint cK = child->connector().key();
      //update polynomial degree to child cells
      p_[cK] = concepts::Array<ushort>(dim);
      for(ushort i = 0 ; i < dim; ++i)
        p_[cK][i] = p_[K][i];
 
      pred_[cK] = pred_[K];
      //add new cell to cell list
      cells_[cK] = child;
    }
    //erase old polynomial degree information
    p_.erase(p_.find(K));
    //erase old prediction info
    pred_.erase(pred_.find(K));
    //erase coarse cell
    cells_.erase(cells_.find(K));
 
 
  }
  //clear this map, since those cells no longer are finest cells in space
  hRefCurCells_.erase(hRefCurCells_.begin(), hRefCurCells_.end());
 
 
}
 
protected:
    virtual std::ostream& info(std::ostream& os) const{
      os << " Prediction[ type = hprefinement Strategy ]" << std::endl;
      os << "refinement chooses : "<< std::endl;
      return os << this->adjusts_;
      //further output : history
      // the parameters etc
    }
 
 
private:
  //error decay value for h refinement
  Real hdecay_;
  //error decay value for p refinement
  Real pdecay_;
  //error decay value for when no is applied refinement
  Real noRefdecay_;
 
  //map from cell Key to prediction
  concepts::HashMap<Real> pred_;
 
  //map from cellKey to its highest polynomial degree, atm just isotropic, therefore only one ushort
  concepts::HashMap<concepts::Array<ushort> > p_;
 
  //current cells of given space
  concepts::HashMap<const concepts::Cell* > cells_;
 
  //map from (static) key to belonging that is hrefined in the space and therefore has kid elements
  //that need updated prediction values in the moment when the space is rebuilded
  //current and updated with build / apply routine
  concepts::Sequence<const concepts::Cell* > hRefCurCells_;
 
 
};
 
 
 
 
 
 
 
//template<class F = Real, ushort dim = 2>
//class Embedding : public HP_Refinement<F,dim>{
//public:
//
//  // H1 just H1 norm
//  // H2Partial2 is  H2 norm with just dxx and dyy part of H2 seminorm
//  enum indicator {H1=0 , H2Partial2};
//
// *
// * @param reg regularity barrier parameter that decides which refinement should be applied
// *        i.e for indicator smaller as reg we apply h-refinement
// *            else p refinement
// */
//  Embedding(const concepts::SpaceOnCells<F>& space, const concepts::Vector<F>& sol,const Real reg, enum indicator ind = H2Partial2)
//    :spc_(space), sol_(sol), reg_(reg), ind_(ind)
//  {
//
//    const hp2D::hpAdaptiveSpace<F>* spc = dynamic_cast<const hp2D::hpAdaptiveSpace<F>* >(&space);
//        if(spc){
//          //information about the polynomial degree
//          hp2D::hpFull& hpfull = spc->prebuild();
//          std::auto_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
//            while (*sc) {
//              concepts::ElementWithCell<Real>& elm = (*sc)++;
//
//              const hp2D::Quad<Real>* quad = dynamic_cast<const hp2D::Quad<Real>*>(&elm);
//              conceptsAssert(quad, concepts::Assertion());
//              //key of current element
//              uint K = elm.cell().connector().key().key();
//              //add cell information
//              quads_[K] = quad;
//            }
//        }else{
//        throw conceptsException(concepts::MissingFeature("Currently just hpAdaptiveSpace supported"));
//        }
//  }
//
//
//
//void buildRefinement(const Marking<F>& mark,  const concepts::LocalEstimator<F>& estimator){
//
//
//    //all marked Cells (its keys)
//    const concepts::Set<uint>& marked = mark.marks();
//
//    //diameter map from each Element /TODO: Just take part via Set constructor
//    concepts::Diameter h_K(spc_,concepts::Diameter::PLANE,true);
//
//    //potential used function evaluations
//    concepts::ElementFormulaVector<1> apprxSol(spc_, sol_, hp2D::Value<Real>());
//    concepts::ElementFormulaVector<2> apprxGrad(spc_, sol_, hp2D::Grad<Real>());
//    concepts::ElementFormulaVector<1> partial_xx(spc_, sol_, hp2D::Partial_xx<Real,hp2D::Quad<Real> >());
//    concepts::ElementFormulaVector<1> partial_yy(spc_, sol_, hp2D::Partial_yy<Real,hp2D::Quad<Real> >());
//
//    //check all cells if they are marked
//    concepts::Set<uint>::const_iterator mIter = marked.begin();
//    for( ; mIter != marked.end(); ++mIter){
//      //current marked cell key
//      uint K = *mIter;
//      //cell is marked for refinement
//
//      const hp2D::Quad<Real>* quad = quads_[K];
//
//      Real indicator = 0;
//      //TODO: build static change of point rate and method
//      estimator::LInfinity maximum(*quad, sol_,1,2);
//      // L^2(K) norm squared
//      Real L_2_K = concepts::L2product(*quad, apprxSol);
//      // H1 seminorm squared
//      Real H_1_K = concepts::L2product(*quad, apprxGrad);
//
//
//      if(ind_ == H1){
//      indicator = std::pow(maximum.get(),2) / ( std::pow( 0.5 * h_K(K), -2) * L_2_K + H_1_K );
//        //DEBUGL(1, "ind H (K = " << K <<")  = "<< indicator);
//      }
//      else if (ind_ == H2Partial2){
//        // H2 partial xx seminorm
//        Real H_2xx = concepts::L2product(*quad, partial_xx);
//        // H2 partial yy seminorm
//        Real H_2yy = concepts::L2product(*quad, partial_yy);
//          indicator = std::pow(maximum.get(), 2) / (std::pow(
//            0.5 * h_K(K), -2) * L_2_K + H_1_K + std::pow(0.5 * h_K(
//            K), 2) * (H_2xx + H_2yy));
//        //DEBUGL(1, "ind G (K = " << K <<")  = "<< indicator);
//      }else{
//        throw conceptsException(concepts::MissingFeature("Other indicator types not implemented"));
//      }
//
//
//      // set adaptive strategy
//      if( indicator < reg_ ){
//        this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
//        DEBUGL(1, "ind (K = " << K <<")  = "<< indicator << "  -->  h-refine");
//      }
//      else{
//        this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(0, 1);
//        DEBUGL(1, "ind (K = " << K <<")  = "<< indicator << "  -->  p-refine");
//      }
//  }//loop over all marked cells
//}//build refinement
//
//protected:
//  virtual std::ostream& info(std::ostream& os) const{
//    os << " Embedding[ type = hp - Refinement Strategy :" << std::endl;
//    os << "Indicator type :"<< ( (ind_) ? "H2Partial2" : "H1" ) << std::endl;
//    return os << "regularity barrier = "<< reg_ << "]";
//      //further output : history
//      // the parameters etc
//    }
//
//
//
//private:
//  const concepts::SpaceOnCells<F>& spc_;
//
//  const concepts::Vector<F>& sol_;
//
//  concepts::HashMap<const hp2D::Quad<Real>* > quads_;
//  //regularity barrier
//  Real reg_;
//
//  //indicator type
//  enum indicator ind_;
//
//};
 
 
 
 
 
 
 
 
 
 
 
template<class F = Real, ushort dim = 2>
class AprioriVertex : public HP_Refinement<F,dim>{
public:
 
  AprioriVertex(const concepts::SpaceOnCells<F>& space, const uint attrb){
    attrbs_.insert(attrb);
    //this easily extends to non hp2D::hpAdaptiveSpace
    std::unique_ptr<typename concepts::SpaceOnCells<F>::Scanner> sc(space.scan());
        while (*sc) {
          concepts::ElementWithCell<F>& elm = (*sc)++;
          //key of current element
          uint K = elm.cell().connector().key().key();
          //add cell information
          cells_[K] = &elm.cell();
        }
  }
 
 
  virtual ~AprioriVertex(){}
 
  void addVertex(const uint attrb){
    attrbs_.insert(attrb);
  }
 
  //to generalize, here estimator is not needed
  void buildRefinement(const Marking& mark,  const concepts::LocalEstimator<F>& estimator){
    const concepts::Set<uint>& marks = mark.marks();
    for(concepts::Set<uint>::const_iterator mIter = marks.begin(); mIter != marks.end(); ++mIter ){
      //key of marked element
      uint K = *mIter;
 
//      this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
 
      const concepts::Connector2* cntr =
           dynamic_cast<const concepts::Connector2*>(&(cells_[K]->connector()));
      if(cntr){
        bool pRef = true;
        //current vertex
        const concepts::Connector0* vertex = 0;
        for( uint i = 0; (vertex = cntr->vertex(i) ) != 0; ++i){
          //quad has a singular marked vertex
          if(attrbs_.exist(vertex->attrib().attrib())){
            this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
            pRef = false;
            break;
          }
        }
        //if no h refinement is applied
        if(pRef)
          this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(0, 1);
 
      }else{
        throw conceptsException(concepts::MissingFeature("Currently Only 2d elements supported"));
      }
    }//loop over marked elements
  }
 
 
protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << " AprioriVertex[ type = hp - Refinement Strategy ]";
      //further output : history
      // the parameters etc
    }
 
 
private:
 
  //set of all  vertice attributes at which h-refinement should be applied
  concepts::Set<uint> attrbs_;
 
  //current cells of given space
  concepts::HashMap<const concepts::Cell* > cells_;
 
};
 
 
}
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
#endif // refinement_hh