quadrature_lib.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__quadrature_lib_h
#define dealii__quadrature_lib_h


#include <deal.II/base/config.h>
#include <deal.II/base/quadrature.h>

DEAL_II_NAMESPACE_OPEN


template <int dim>
class QGauss : public Quadrature<dim>
{
public:
  QGauss (const unsigned int n);
};


template <int dim>
class QGaussLobatto : public Quadrature<dim>
{
public:
  QGaussLobatto(const unsigned int n);

protected:
  std::vector<long double>
  compute_quadrature_points (const unsigned int q,
                             const int alpha,
                             const int beta) const;

  std::vector<long double>
  compute_quadrature_weights (const std::vector<long double> &x,
                              const int alpha,
                              const int beta) const;

  long double JacobiP(const long double x,
                      const int alpha,
                      const int beta,
                      const unsigned int n) const;

  static long double gamma(const unsigned int n);
};



template <int dim>
class QMidpoint : public Quadrature<dim>
{
public:
  QMidpoint ();
};


template <int dim>
class QSimpson : public Quadrature<dim>
{
public:
  QSimpson ();
};



template <int dim>
class QTrapez : public Quadrature<dim>
{
public:
  QTrapez ();
};



template <int dim>
class QMilne : public Quadrature<dim>
{
public:
  QMilne ();
};


template <int dim>
class QWeddle : public Quadrature<dim>
{
public:
  QWeddle ();
};



template <int dim>
class QGaussLog : public Quadrature<dim>
{
public:
  QGaussLog(const unsigned int n,
            const bool revert=false);

private:
  static
  std::vector<double>
  get_quadrature_points(const unsigned int n);

  static
  std::vector<double>
  get_quadrature_weights(const unsigned int n);
};




template <int dim>
class QGaussLogR : public Quadrature<dim>
{
public:
  QGaussLogR(const unsigned int n,
             const Point<dim> x0 = Point<dim>(),
             const double alpha = 1,
             const bool factor_out_singular_weight=false);

#ifdef DEAL_II_WITH_CXX11

  QGaussLogR(QGaussLogR<dim> &&) = default;
#endif

protected:
  const double fraction;
};


template <int dim>
class QGaussOneOverR : public Quadrature<dim>
{
public:
  QGaussOneOverR(const unsigned int n,
                 const Point<dim> singularity,
                 const bool factor_out_singular_weight=false);
  QGaussOneOverR(const unsigned int n,
                 const unsigned int vertex_index,
                 const bool factor_out_singular_weight=false);
private:
  static unsigned int quad_size(const Point<dim> singularity,
                                const unsigned int n);
};



template <int dim>
class QSorted : public Quadrature<dim>
{
public:
  QSorted (const Quadrature<dim> &quad);

private:
  bool compare_weights(const unsigned int a,
                       const unsigned int b) const;
};

template <int dim>
class QTelles: public Quadrature<dim>
{
public:
  QTelles (const Quadrature<1> &base_quad, const Point<dim> &singularity);
  QTelles (const unsigned int n, const Point<dim> &singularity);

};

template <int dim>
class QGaussChebyshev : public Quadrature<dim>
{
public:
  QGaussChebyshev(const unsigned int n);

private:
  static std::vector<double>
  get_quadrature_points(const unsigned int n);

  static std::vector<double>
  get_quadrature_weights(const unsigned int n);

};


template <int dim>
class QGaussRadauChebyshev : public Quadrature<dim>
{
public:
  /* EndPoint is used to specify which of the two endpoints of the unit interval
   * is used also as quadrature point
   */
  enum EndPoint
  {
    left,
    right
  };
  QGaussRadauChebyshev(const unsigned int n,
                       EndPoint ep=QGaussRadauChebyshev::left);

#ifdef DEAL_II_WITH_CXX11

  QGaussRadauChebyshev(QGaussRadauChebyshev<dim> &&) = default;
#endif

private:
  const EndPoint ep;
  static std::vector<double>
  get_quadrature_points(const unsigned int n, EndPoint ep);

  static std::vector<double>
  get_quadrature_weights(const unsigned int n, EndPoint ep);

};

template <int dim>
class QGaussLobattoChebyshev : public Quadrature<dim>
{
public:
  QGaussLobattoChebyshev(const unsigned int n);

private:
  static std::vector<double>
  get_quadrature_points(const unsigned int n);

  static std::vector<double>
  get_quadrature_weights(const unsigned int n);

};

/* -------------- declaration of explicit specializations ------------- */

template <> QGauss<1>::QGauss (const unsigned int n);
template <> QGaussLobatto<1>::QGaussLobatto (const unsigned int n);
template <>
std::vector<long double> QGaussLobatto<1>::
compute_quadrature_points(const unsigned int, const int, const int) const;
template <>
std::vector<long double> QGaussLobatto<1>::
compute_quadrature_weights(const std::vector<long double> &, const int, const int) const;
template <>
long double QGaussLobatto<1>::
JacobiP(const long double, const int, const int, const unsigned int) const;
template <>
long double
QGaussLobatto<1>::gamma(const unsigned int n);

template <> std::vector<double> QGaussLog<1>::get_quadrature_points(const unsigned int);
template <> std::vector<double> QGaussLog<1>::get_quadrature_weights(const unsigned int);

template <> QMidpoint<1>::QMidpoint ();
template <> QTrapez<1>::QTrapez ();
template <> QSimpson<1>::QSimpson ();
template <> QMilne<1>::QMilne ();
template <> QWeddle<1>::QWeddle ();
template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert);
template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag);
template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag);
template <> QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity);



DEAL_II_NAMESPACE_CLOSE

#endif