explicitResidual.hh

Go to the documentation of this file.

 
#ifndef explicitResidual_hh
#define explicitResidual_hh
 
#include "estimator/estimator.hh"
 
//TODO: Precise the includes
#include "hp2D.hh"
#include "formula/elementFormulaContainer.hh"
#include "basics/exceptions.hh"
 
 
namespace concepts{
 
  template<class F>
  class ExplicitResidual: public LocalEstimator<F>{
 
  public:
 
    virtual ~ExplicitResidual(){};
 
    virtual ExplicitResidual<F>* clone() const = 0;
 
  void addINData(const Set<uint>& nSet, concepts::ElementFormulaContainer<F> frm){
    frms_.push_back(frm);
    iNAttrbs_.push_back(nSet);
  }
 
  void addHNData(const Set<uint>& nSet){
    hNAttrbs_ = hNAttrbs_ || nSet;
  }
  void addHNData(uint attrb) {
    hNAttrbs_.insert(attrb);
  }
 
 
  protected:
 
 
    ExplicitResidual(const SpaceOnCells<F>& spc,const Vector<F>& sol, const ElementFormulaContainer<F> res):
      LocalEstimator<F>(spc, sol),res_(res){};
 
    virtual std::ostream& info(std::ostream& os) const{
      return os << concepts::typeOf(*this);
    }
    //local reference to the residual
    const ElementFormulaContainer<F> res_;
 
    //Neumann formulas
    Sequence<ElementFormulaContainer<F> > frms_;
    //Neumann Attributes inhomogen
    Sequence<Set<uint> > iNAttrbs_;
    //Neumann attributes homogen
    Set<uint> hNAttrbs_;
 
  private:
 
    virtual void computeError_(const ElementFormulaContainer<F>& res, HashMap<Real>& jumpResidual) = 0 ;
    virtual void computeJumpPart_(HashMap<Real>& jumpResidual) const=0;
  };
 
 
}
 
 
 
 
 
namespace hp1D{
 
  //So far just a abstract class
  template<class F>
  class ExplicitResidual1D: concepts::ExplicitResidual<F>{
 
  public:
    ExplicitResidual1D(const concepts::SpaceOnCells<F>& spc,
                       const concepts::Vector<F>& sol,
                       const concepts::ElementFormulaContainer<F>& res);
 
    virtual ~ExplicitResidual1D();
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
 
  private:
    virtual void computeError_(const concepts::ElementFormulaContainer<F>& res, concepts::HashMap<Real>& jumpResidual);
    virtual void computeJumpPart_(concepts::HashMap<Real>& jumpResidual) const;
  };
 
}
 
 
namespace hp2D {
template<class F>
class ExplicitResidual2D: public concepts::ExplicitResidual<F> {
 
public:
 
  // Flag for the chose of the weight
    // DIAM :  default weight is just the diameter of cell :
    //            w_K = diam(K), w_E = diam(E)
    // DIAM_HP : residual type variant introduced by melenk and wohlmuth with the help of weighted norms
    //
    //
    // DIAM_P: polynomal degree weighted diameter :
    //            w_K = diam(K)/(1+p_K),  w_E = diam(E)/2*(1+p_E)
 
    // further possible weights :
    // DIAM_L  Diameter devided by eigenvalues lambda_K of the diffusion tensor elementwise.
    // where the diffussion tensor IK is assumed to be quasi-monoton.
    // for further details check the Diplomarbeit of Christian Merdon, section 2 and 3.
    //TODO: RETHINK NAME
    enum weight {DIAM, DIAM_HP, DIAM_P};
 
  ExplicitResidual2D(const concepts::SpaceOnCells<F>& spc,
             const concepts::Vector<F>& sol,
             const concepts::ElementFormulaContainer<F> res,
             const concepts::ElementFormulaContainer<F>* a = 0,
             bool square = false,
             enum weight w = DIAM,
             Real alpha = 0);
 
 
  virtual ~ExplicitResidual2D(){};
 
  /*
   * Starts the actual estimation of error.
   * You may start this routine after setting the Neumann boundary info data.
   */
  void compute();
 
  class Distance : public concepts::ElementFormula<Real>{
 
    public:
      //Constructor
      Distance(Real alpha=0) : alpha_(alpha){};
 
      virtual F operator() (const concepts::ElementWithCell<typename concepts::Realtype<F>::type>& elm, const concepts::Real p,const concepts::Real t = 0.0) const{
        throw concepts::MissingFeature("Operator not implemented.");
        return 0.0;
      }
 
      virtual F operator() (const concepts::ElementWithCell<typename concepts::Realtype<F>::type>& elm, const concepts::Real2d& p,const concepts::Real t = 0.0) const{
 
        // __________
        // |  :     |
        // |--x p   |
        // |        |
        // |________|
        //dist = min(x,y,1-x,1-y)
        concepts::Real dist = std::min(std::min(p[0],p[1]), std::min(1-p[0],1-p[1]));
        return std::pow(dist, alpha_);
 
      }
      virtual F operator() (const concepts::ElementWithCell<typename concepts::Realtype<F>::type>& elm, const concepts::Real3d& p,const concepts::Real t = 0.0) const{
        //just take p0 and p1 since this class applicates for hp2D::IntegrableQuad
        concepts::Real dist = std::min(std::min(p[0],p[1]), std::min(1-p[0],1-p[1]));
        return std::pow(dist, alpha_);
      }
 
      virtual Distance* clone() const{return new Distance(*this);}
 
      virtual ~Distance(){};
 
    protected:
      virtual std::ostream& info(std::ostream& os) const{
        return os<<concepts::typeOf(*this) <<"( alpha = "<<alpha_<<")";
      }
 
    private:
      //representing alpha/2
      Real alpha_;
    };
 
 
 
    class EdgeWeight : public concepts::ElementFormula<Real>{
 
    public:
      //Constructor
      EdgeWeight(Real alpha=0) : alpha_(alpha){};
 
      virtual F operator() (const concepts::ElementWithCell<typename concepts::Realtype<F>::type>& elm, const concepts::Real p,const concepts::Real t = 0.0) const{
        // ((1-x)*x)^(alpha)
        return std::pow((1-p)*p, alpha_);
      }
 
      //for neumanntracelements
      virtual F operator() (const concepts::ElementWithCell<typename concepts::Realtype<F>::type>& elm, const concepts::Real2d& p,const concepts::Real t = 0.0) const{
        //take p[0] since this class applicates on neumanntraceElements which are hp1D::IntegrableElm
        // ((1-x)*x)^(alpha)
        return std::pow((1-p[0])*p[0], alpha_);
 
      }
      virtual F operator() (const concepts::ElementWithCell<typename concepts::Realtype<F>::type>& elm, const concepts::Real3d& p,const concepts::Real t = 0.0) const{
        //take p[0] since this class applicates on neumanntraceElements which are hp1D::IntegrableElm
        // ((1-x)*x)^(alpha)
        return std::pow((1-p[0])*p[0], alpha_);
      }
 
      virtual EdgeWeight* clone() const{return new EdgeWeight(*this);}
 
      virtual ~EdgeWeight(){};
 
    protected:
      virtual std::ostream& info(std::ostream& os) const{
        return os << concepts::typeOf(*this) << "( alpha = "<<alpha_<<")";
      }
 
    private:
      //representing alpha
      Real alpha_;
    };
 
    virtual const Real operator()() const {return 0.;}
 
    virtual ExplicitResidual2D<F>* clone() const{return new ExplicitResidual2D<F>(*this);}
 
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this);
    }
 
  private:
 
    //flag that marks all elements are squares, then diameter computation is faster.
    bool square_;
    //isotropic scalar coefficient function in higest derivative operator
    const concepts::ElementFormulaContainer<F>* a_;
    //weight type
    enum weight weight_;
 
    //weight function exponent
    Real alpha_;
 
    //faster data structure for accessing all neumann boundary attributes
    concepts::Set<uint> allNattrb_;
 
    virtual void computeError_(const concepts::ElementFormulaContainer<F>& res, concepts::HashMap<Real>& jumpResidual);
    virtual void computeJumpPart_(concepts::HashMap<Real>& jumpResidual) const;
 
  };
 
 
 
 
 
 
 
 
 
 
 
 
} //namespace 2D
 
 
//namespace hp3D{
//
//class ExplicitResidual3D: concepts::ExplicitResidual{
//public:
//
//  template<typename F>
//    ExplicitResidual3D(const concepts::SpaceOnCells<F>& spc, const concepts::Vector<F>& sol, const concepts::ElementFormula<F>& res);
//
//protected:
//  virtual std::ostream& info(std::ostream& os) const{
//    return os << "Explicit Residual Error Estimator()";
//  }
//
//
//private:
//  virtual void computeError_(const concepts::ElementFormula<F>& res, concepts::HashMap<Real>& jumpResidual);
//  virtual void computeJumpPart_(concepts::HashMap<Real>& jumpResidual) const;
//};
//
//}
 
#endif // explicitResidual_hh