numerics/quadrature.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_QUADRATURE_H
#define _MATHTL_QUADRATURE_H

#include <utils/array1d.h>
#include <utils/function.h>
#include <geometry/point.h>

namespace MathTL
{
  template <unsigned int DIM>
  class QuadratureRule
  {
  public:
    typedef QuadratureRule<DIM-1> SubQuadratureRule;

    QuadratureRule();

    QuadratureRule(const QuadratureRule<DIM>& Q);

    QuadratureRule(const SubQuadratureRule& Q,
                   const QuadratureRule<1>& Q1);
    
    virtual ~QuadratureRule();

    unsigned int get_N() const;

    void get_points(Array1D<Point<DIM> >& points) const;

    void get_weights(Array1D<double>& weights) const;

    double integrate(const Function<DIM, double>& f) const;

    double integrate(const Function<DIM, double>& f,
                     const Point<DIM>& a, const Point<DIM>& b) const;

  protected:
    Array1D<Point<DIM> > points_;

    Array1D<double> weights_;
  };

  class MidpointRule : public QuadratureRule<1>
  {
  public:
    MidpointRule();
  };

  class ClosedNewtonCotesRule : public QuadratureRule<1>
  {
  public:
    ClosedNewtonCotesRule(const unsigned int N);
  };

  class TrapezoidalRule: public QuadratureRule<1>
  {
  public:
    TrapezoidalRule();
  };

  class SimpsonRule: public QuadratureRule<1>
  {
  public:
    SimpsonRule();
  };

  template <unsigned int DIM>
  class CompositeRule
    : public QuadratureRule<DIM>
  {
  public:
    CompositeRule(const QuadratureRule<1>& Q,
                  const unsigned int N = 1);
    
  protected:
    const QuadratureRule<1> Q_;

    unsigned int N_;

  private:
    static bool uses_both_endpoints(const QuadratureRule<1>& Q);
  };
}


// include implementation of inline functions and specializations
#include <numerics/quadrature.cpp>

#endif