integration/quadRule.hh

Go to the documentation of this file.

/* 1D quadrature rules
 */

#ifndef quadRule_hh
#define quadRule_hh

#include <map>

#include "basics/typedefs.hh"
#include "basics/defines.hh"
#include "basics/outputOperator.hh"
#include "quadrature.hh"

namespace concepts {

  // ******************************************************** QuadratureRule **

  class QuadratureRule : public OutputOperator {
  private:
    QuadratureRule(const QuadratureRule& other);
    QuadratureRule& operator=(const QuadratureRule& other);
  public:
    QuadratureRule() { }
    virtual ~QuadratureRule() { }

    virtual const Real* abscissas() const = 0;

    virtual const Real* weights() const = 0;

    virtual uint n() const = 0;

    void printRule();
  };

  // ************************************************* QuadratureRuleDynamic **

  class QuadratureRuleDynamic : public QuadratureRule {
  public:
    QuadratureRuleDynamic(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_;
  };

  // ******************************************** QuadratureRuleGaussLobatto **

  class QuadratureRuleGaussLobatto : public QuadratureRuleDynamic {
  public:
    QuadratureRuleGaussLobatto(uint n) : rule_(n) {
      abscissas_ = rule_.abscissas();
      weights_ = rule_.weights();
    }
    virtual ~QuadratureRuleGaussLobatto();
    virtual uint n() const { return rule_.n(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Quadrature<0> rule_;
  };

  // ********************************************* QuadratureRuleGaussJacobi **

  class QuadratureRuleGaussJacobi : public QuadratureRuleDynamic {
  public:
    QuadratureRuleGaussJacobi(uint n) : rule_(n) {
      abscissas_ = rule_.abscissas();
      weights_ = rule_.weights();
    }
    virtual ~QuadratureRuleGaussJacobi();
    virtual uint n() const { return rule_.n(); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Quadrature<4> rule_;
  };

  // ************************************************** QuadratureRuleTrapeze **

  /* 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 QuadratureRuleTrapeze : public QuadratureRuleDynamic {
  public:
    QuadratureRuleTrapeze(uint n)
      : rule_(n), shAbscissas_(0), shWeights_(0) {
      abscissas_ = rule_.abscissas();
      weights_   = rule_.weights();
    }
    QuadratureRuleTrapeze(uint n, Real x0, Real xend);
    virtual ~QuadratureRuleTrapeze();

    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 OutputOperator {
  public:
    QuadRuleFactory(enum intRule type = GAUSS_JACOBI, uint constPoints = 10,
                    uint addPoints = 2, bool constant = false);

    QuadratureRule* 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:
    enum intRule integrationType_;
    uint constNumerOfPoints_;
    uint addNumberOfPoints_;
    bool useConstantNumberOfPoints_;

    uint cnt_;
  };

} // namespace concepts

#endif // quadRule_hh