quadRule.hh

Go to the documentation of this file.

 
#ifndef quadRule_hh
#define quadRule_hh
 
#include <map>
 
#include "basics/typedefs.hh"
#include "basics/defines.hh"
#include "basics/outputOperator.hh"
#include "basics/vectorsMatrices.hh"
#include "quadrature.hh"
#include "geometry/cell2D.hh"
#include "toolbox/hashMap.hh"
 
//Precision libs small abcissas
//#include <gmp.h>
//#include <gmpxx.h>
 
namespace concepts {
 
  // ******************************************************** QuadratureRule **
 
  class QuadratureRule1d : public virtual OutputOperator {
  private:
    QuadratureRule1d(const QuadratureRule1d& other);
    QuadratureRule1d& operator=(const QuadratureRule1d& other);
  public:
    QuadratureRule1d() { }
    virtual ~QuadratureRule1d() { }
 
    virtual const Real* abscissas() const = 0;
 
    virtual const Real* weights() const = 0;
 
    virtual uint n() const = 0;
 
    virtual void printRule();
  };
 
  // ************************************************* QuadratureRule1dDynamic **
 
  class QuadratureRule1dDynamic : public QuadratureRule1d {
  public:
    QuadratureRule1dDynamic(const Real* abscissas = 0, const Real* weights = 0) :
      abscissas_(abscissas), weights_(weights) {}
    virtual const Real* abscissas() const { return abscissas_; }
    virtual const Real* weights() const { return weights_; }
  protected:
    const Real* abscissas_;
    const Real* weights_;
  };
 
  // ******************************************** QuadratureRule1dGaussLobatto **
 
  class QuadratureRule1dGaussLobatto : public QuadratureRule1dDynamic {
  public:
    QuadratureRule1dGaussLobatto(uint n) : rule_(n) {
      abscissas_ = rule_.abscissas();
      weights_ = rule_.weights();
    }
    virtual ~QuadratureRule1dGaussLobatto();
    virtual uint n() const { return rule_.n(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Quadrature<0> rule_;
  };
 
  // ********************************************* QuadratureRule1dGaussJacobi **
 
  class QuadratureRule1dGaussJacobi : public QuadratureRule1dDynamic {
  public:
    QuadratureRule1dGaussJacobi(uint n) : rule_(n) {
      abscissas_ = rule_.abscissas();
      weights_ = rule_.weights();
    }
    virtual ~QuadratureRule1dGaussJacobi();
    virtual uint n() const { return rule_.n(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Quadrature<4> rule_;
  };
 
  // ************************************************** QuadratureRule1dTrapeze **
 
  /* Trapeze quadrature rule including both endpoints -1 and 1.
     \f[ \int_{-1}^1 f(x) \, dx \approx \sum_{i=0}^(n-1) w_i f(x_i) \f]
     is only exact for \f$f \in P_{1}\f$.
     n must be greater or equal to 2.
     \n
     The abscissas \f$x_i = -1 + 2i/(n-1)\f$ are equidistant
     and the weights are
     \f[w_i = \frac{1}{n-1}\left\{\begin{array}{ll}1 & i \in \{0,n\}\\2 & \mbox{otherwise}\end{array}\right.\f]
 
     The computations and the storage of the values are done by the
     class Quadrature with template parameter 4. The difference
     between this class and Quadrature is that it is in a class
     hierarchy of quadrature rules. This has advantages when
     dynamically switching quadrature rules is needed. On the other
     hand, this class returns the values via a virtual function call
     \c abscissas() and \c weights() should therefore not be called
     to often (inside loops etc.).
 
     The abscissas are good for graphical output.
 
     @see Quadrature
     @author Kersten Schmidt, 2004
  */
 
  class QuadratureRule1dTrapeze : public QuadratureRule1dDynamic {
  public:
    QuadratureRule1dTrapeze(uint n)
      : rule_(n), shAbscissas_(0), shWeights_(0) {
      abscissas_ = rule_.abscissas();
      weights_   = rule_.weights();
    }
    QuadratureRule1dTrapeze(uint n, Real x0, Real xend);
    virtual ~QuadratureRule1dTrapeze();
 
    virtual const Real* abscissas() const { 
      return abscissas_ ? abscissas_ : shAbscissas_;
    }
    virtual const Real* weights() const { 
      return weights_ ? weights_ : shWeights_;
    }
 
    virtual uint n() const { return rule_.n(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Quadrature<5> rule_;
    Real *shAbscissas_, *shWeights_;
  };
 
  // ******************************************************* QuadRuleFactory **
 
  class QuadRuleFactory : public virtual OutputOperator {
  public:
    QuadRuleFactory(enum intRule type = GAUSS_JACOBI, uint constPoints = 10,
                    uint addPoints = 2, bool constant = false);
 
    QuadratureRule1d* operator()(const ushort p = 1) const;
 
    void set(enum concepts::intRule rule, bool constant, uint points);
 
    void reset();
 
    inline const uint count() const { return cnt_; }
 
    const std::string integrationRule() const;
 
    enum intRule type() const { return integrationType_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;    
  private:
    // TODO: adapt documentation: it should be set(), not setIntegrationRule()
    enum intRule integrationType_;
    uint constNumerOfPoints_;
    uint addNumberOfPoints_;
    bool useConstantNumberOfPoints_;
 
    uint cnt_;
  };
 
 
 
  // ********************************************************* QuadratureRule **
 
  class QuadratureRule : public virtual OutputOperator {
  private:
    QuadratureRule(const QuadratureRule1d& other);
    QuadratureRule& operator=(const QuadratureRule1d& other);
  public:
    QuadratureRule(){}
    virtual ~QuadratureRule() { }
 
    virtual const Real* abscissas(uint i) const = 0;
 
    virtual const Real* weights(uint i) const = 0;
 
    virtual const uint n(uint i=0) const = 0 ;
 
    //delivers ith quadraturepoint q on reference cell [-1,1]^n and ith weight, must be implemented in child classes
    virtual bool quadratureData(const uint i, Real3d& q, Real& w) const = 0;
   };
 
  // ******************************************************* QuadratureRule2d **
 
  class QuadratureRule2d : public QuadratureRule {
  private:
    QuadratureRule2d(const QuadratureRule2d& other);
    QuadratureRule2d& operator=(const QuadratureRule2d& other);
 
  public:
    QuadratureRule2d(): intX_(0), intY_(0) {};
    virtual ~QuadratureRule2d() { }
 
    virtual const Real* abscissas(uint i) const{
      conceptsAssert(i < 2, concepts::Assertion());
      return i ? intY_ : intX_;
    }
 
    virtual const Real* weights(uint i = 0) const = 0;
    virtual const uint n(uint i = 0) const = 0;
 
    virtual const bool domain() const = 0 ;
 
    virtual const bool tensor() const = 0 ;
 
  protected:
    const Real* intX_;
    const Real* intY_;
  };
 
  // ********************************************* QuadratureRule2dQuadTensor **
 
  class QuadratureRule2dQuadTensor : public QuadratureRule2d {
  public:
    QuadratureRule2dQuadTensor(const concepts::QuadratureRule1d* xRule,
                               const concepts::QuadratureRule1d* yRule);
 
    virtual const Real* weights(uint i) const ;
 
    virtual const uint n(uint i) const;
 
    virtual bool quadratureData(const uint i, Real3d& q, Real& w) const;
 
    //tensored rules are supposed to be computed on [-1,1]^2
    virtual const bool domain() const { return true; }
 
    // returns true as class is tensor structured
    virtual const bool tensor() const { return true; }
 
    //returns 1d quadrature rule in x direction
    //@pre quadrature rule exists by constructor
    const concepts::QuadratureRule1d* xRule() const{return xRule_.get();}
 
    //returns 1d quadrature rule in y direction
    //@pre quadrature rule exists by constructor
    const concepts::QuadratureRule1d* yRule() const{return yRule_.get();}
 
 
  protected:
      virtual std::ostream& info(std::ostream& os) const;
  private:
    std::unique_ptr<const concepts::QuadratureRule1d> xRule_;
    std::unique_ptr<const concepts::QuadratureRule1d> yRule_;
  };
 
  // ********************************************** QuadratureRule2dQuadDuffy **
 
  class QuadratureRule2dQuadDuffy : public QuadratureRule2d {
  public:
  QuadratureRule2dQuadDuffy(Real beta, uint nx, uint ny, uint noVtx = 0);
 
  virtual ~QuadratureRule2dQuadDuffy();
 
    virtual const Real* weights(uint i = 0) const{
      conceptsAssert(i < 1, concepts::Assertion());
      return weights_;
    }
 
   const uint n(uint i = 0) const{
      conceptsAssert(i < 2, concepts::Assertion());
      return n_;
   }
 
   // returns i-th quadrature point \c q = (x_i,y_i) and i-th weight \c w
   virtual bool quadratureData(const uint i, Real3d& q, Real& w) const;
 
   // abscissas and weights are computed in [0,1]^2
   virtual const bool domain() const { return false; }
 
   // non-tensor structure
   virtual const bool tensor() const { return false; }
 
   // prints the quadrature points
   void printData() const;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
 
  private :
      void compute_(const Real2d S,const Real2d B,const Real2d C,
            const Real beta, const uint nx, const uint ny,
            Real*& abscXP, Real*& abscYP, Real*& wP);
 
      // weights for each integration Point
      const Real* weights_;
      //duffy transformation parameter
      Real beta_;
      //number of quadrature points n_ = 2*(nx+ny)
      uint n_;
  };
 
  // ************************************************** QuadRuleFactoryBase2d **
 
  class QuadRuleFactoryBase2d: public virtual OutputOperator {
 
  private :
    QuadRuleFactoryBase2d(const QuadRuleFactoryBase2d& other);
    QuadRuleFactoryBase2d& operator=(const QuadRuleFactoryBase2d& other);
  public :
    QuadRuleFactoryBase2d() {}
    virtual QuadratureRule2d* operator()(const concepts::QuadNd& cell,
                                         const ushort* p) const = 0;
    
    // print of the quadrature rule
    virtual const std::string integrationRule() const = 0;
  };
 
 
  // ************************************************ QuadRuleFactoryTensor2d **
 
  //this class is a direct copy of the QuadRuleFactory
  class QuadRuleFactoryTensor2d : public QuadRuleFactoryBase2d {
  public:
    QuadRuleFactoryTensor2d(enum intRule type = GAUSS_JACOBI, uint constPoints = 10,
                            uint addPoints = 2, bool constant = false);
    
    virtual QuadratureRule2d* operator()(const concepts::QuadNd& cell, const ushort* p) const;
    
 
    virtual void setTensor(enum concepts::intRule rule, bool constant,
                           uint points);
 
    void reset();
    
    inline const uint count() const {return cnt_; }
    
    enum intRule type() const {return integrationType_; }
    
    virtual const std::string integrationRule() const;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
 
    // Tensor information datas
 
    // TODO: adapt documentation: it should be set(), not setIntegrationRule()
    enum intRule integrationType_;
    uint constNumerOfPoints_;
    uint addNumberOfPoints_;
    bool useConstantNumberOfPoints_;
 
    uint cnt_;
  };
 
  // ************************************************ QuadRuleFactoryTensorDuffy2d **
 
  class QuadRuleFactoryTensorDuffy2d : public QuadRuleFactoryTensor2d {
  private:
    struct DuffyData : public OutputOperator {
      // duffy parameter
      Real beta;
      // number of integration points for underlying tensor integration rule.
      uint npx;
      uint npy;
      
      DuffyData(const Real beta = 1 ,const uint npx = 20, const uint npy = 20) : beta(beta), npx(npx), npy(npy) {}
 
      DuffyData& operator=(const DuffyData& dd) {
        if (this != &dd){
          this->beta = dd.beta;
          this->npx = dd.npx;
          this->npy = dd.npy;
        }
        return *this;
      }
 
    protected:
      virtual std::ostream& info(std::ostream& os) const {
        return os << "DuffyData{beta = "<< beta <<" , (npx,npy) = ("<< npx << " , " << npy << ")}";
      }
    };
 
  public:
    QuadRuleFactoryTensorDuffy2d(enum intRule tensorType = GAUSS_JACOBI, uint constPointsT = 10,
                                 uint addPointsT = 2, bool constantT = false);
 
 
    inline void addSingularPoint(const uint vAttr,const Real beta, uint npx = 20 , uint npy = 20){
      vtxData_[vAttr]= DuffyData(beta,npx,npy);
    }
 
    virtual QuadratureRule2d* operator()(const concepts::QuadNd& cell, const ushort* p) const;
 
    virtual const std::string integrationRule() const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    
    // data from each attribute to its duffy quadrature points necessarc data, i.e. beta, npx, npy
    concepts::HashMap<DuffyData> vtxData_;
  };
 
 
} // namespace concepts
 
#endif // quadRule_hh