quadrature_lib.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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);
};



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);

  QGaussLogR(QGaussLogR<dim> &&) noexcept = default;

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);
};


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);

  QGaussRadauChebyshev(QGaussRadauChebyshev<dim> &&) noexcept = default;

private:
  const EndPoint ep;
};

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

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

  Quadrature<dim>
  compute_affine_transformation(
    const std::array<Point<dim>, dim + 1> &vertices) const;
};

class QTrianglePolar : public QSimplex<2>
{
public:
  QTrianglePolar(const Quadrature<1> &radial_quadrature,
                 const Quadrature<1> &angular_quadrature);

  QTrianglePolar(const unsigned int n);
};

class QDuffy : public QSimplex<2>
{
public:
  QDuffy(const Quadrature<1> &radial_quadrature,
         const Quadrature<1> &angular_quadrature,
         const double         beta = 1.0);

  QDuffy(const unsigned int n, const double beta);
};

template <int dim>
class QSplit : public Quadrature<dim>
{
public:
  QSplit(const QSimplex<dim> &base, const Point<dim> &split_point);
};

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

template <>
QGauss<1>::QGauss(const unsigned int n);
template <>
QGaussLobatto<1>::QGaussLobatto(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