quadrature_lib.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/polynomial.h>

#include <cmath>
#include <limits>
#include <algorithm>
#include <functional>


DEAL_II_NAMESPACE_OPEN


// please note: for a given dimension, we need the quadrature formulae
// for all lower dimensions as well. That is why in this file the check
// is for deal_II_dimension >= any_number and not for ==



template <>
QGauss<0>::QGauss (const unsigned int)
  :
  // there are n_q^dim == 1
  // points
  Quadrature<0> (1)
{
  // the single quadrature point gets unit
  // weight
  this->weights[0] = 1;
}



template <>
QGaussLobatto<0>::QGaussLobatto (const unsigned int)
  :
  // there are n_q^dim == 1
  // points
  Quadrature<0> (1)
{
  // the single quadrature point gets unit
  // weight
  this->weights[0] = 1;
}



template <>
QGauss<1>::QGauss (const unsigned int n)
  :
  Quadrature<1> (n)
{
  if (n == 0)
    return;

  std::vector<long double> points
    = Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);

  for (unsigned int i=0; i<(points.size()+1)/2; ++i)
    {
      this->quadrature_points[i][0] = points[i];
      this->quadrature_points[n-i-1][0] = 1.-points[i];

      // derivative of Jacobi polynomial
      const long double pp = 0.5*(n + 1)*Polynomials::jacobi_polynomial_value(n-1, 1, 1, points[i]);
      const long double x = -1. + 2.*points[i];
      const double w = 1./((1.-x*x)*pp*pp);
      this->weights[i] = w;
      this->weights[n-i-1] = w;
    }
}

namespace internal
{
  namespace QGaussLobatto
  {
    long double gamma(const unsigned int n)
    {
      long double result = n - 1;
      for (int i=n-2; i>1; --i)
        result *= i;
      return result;
    }



    std::vector<long double>
    compute_quadrature_weights(const std::vector<long double> &x,
                               const int alpha,
                               const int beta)
    {
      const unsigned int q = x.size();
      std::vector<long double> w(q);

      const long double factor = std::pow(2., alpha+beta+1) *
                                 gamma(alpha+q) *
                                 gamma(beta+q) /
                                 ((q-1)*gamma(q)*gamma(alpha+beta+q+1));
      for (unsigned int i=0; i<q; ++i)
        {
          const long double s = Polynomials::jacobi_polynomial_value(q-1, alpha, beta, x[i]);
          w[i] = factor/(s*s);
        }
      w[0]   *= (beta + 1);
      w[q-1] *= (alpha + 1);

      return w;
    }
  }
}



template <>
QGaussLobatto<1>::QGaussLobatto (const unsigned int n)
  :
  Quadrature<1> (n)
{
  Assert (n >= 2, ExcNotImplemented());

  std::vector<long double> points
    = Polynomials::jacobi_polynomial_roots<long double>(n-2, 1, 1);
  points.insert(points.begin(), 0);
  points.push_back(1.);
  std::vector<long double> w
    = internal::QGaussLobatto::compute_quadrature_weights(points, 0, 0);

  // scale weights to the interval [0.0, 1.0]:
  for (unsigned int i=0; i<points.size(); ++i)
    {
      this->quadrature_points[i][0] = points[i];
      this->weights[i]              = 0.5*w[i];
    }
}



template <>
QMidpoint<1>::QMidpoint ()
  :
  Quadrature<1>(1)
{
  this->quadrature_points[0] = Point<1>(0.5);
  this->weights[0] = 1.0;
}



template <>
QTrapez<1>::QTrapez ()
  :
  Quadrature<1> (2)
{
  static const double xpts[] = { 0.0, 1.0 };
  static const double wts[]  = { 0.5, 0.5 };

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(xpts[i]);
      this->weights[i] = wts[i];
    };
}



template <>
QSimpson<1>::QSimpson ()
  :
  Quadrature<1> (3)
{
  static const double xpts[] = { 0.0, 0.5, 1.0 };
  static const double wts[]  = { 1./6., 2./3., 1./6. };

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(xpts[i]);
      this->weights[i] = wts[i];
    };
}



template <>
QMilne<1>::QMilne ()
  :
  Quadrature<1> (5)
{
  static const double xpts[] = { 0.0, .25, .5, .75, 1.0 };
  static const double wts[]  = { 7./90., 32./90., 12./90., 32./90., 7./90. };

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(xpts[i]);
      this->weights[i] = wts[i];
    };
}



template <>
QWeddle<1>::QWeddle ()
  :
  Quadrature<1> (7)
{
  static const double xpts[] = { 0.0, 1./6., 1./3., .5, 2./3., 5./6., 1.0 };
  static const double wts[]  = { 41./840., 216./840., 27./840., 272./840.,
                                 27./840., 216./840., 41./840.
                               };

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(xpts[i]);
      this->weights[i] = wts[i];
    };
}


template <>
QGaussLog<1>::QGaussLog(const unsigned int n,
                        const bool revert)
  :
  Quadrature<1> (n)
{
  std::vector<double> p = get_quadrature_points(n);
  std::vector<double> w = get_quadrature_weights(n);

  for (unsigned int i=0; i<this->size(); ++i)
    {
      // Using the change of variables x=1-t, it's possible to show
      // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
      // we can use this quadrature formula also with weight ln|1-x|.
      this->quadrature_points[i] = revert ? Point<1>(1-p[n-1-i]) : Point<1>(p[i]);
      this->weights[i]           = revert ? w[n-1-i] : w[i];
    }
}

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

  switch (n)
    {
    case 1:
      q_points[0] = 0.3333333333333333;
      break;

    case 2:
      q_points[0] = 0.1120088061669761;
      q_points[1] = 0.6022769081187381;
      break;

    case 3:
      q_points[0] = 0.06389079308732544;
      q_points[1] = 0.3689970637156184;
      q_points[2] = 0.766880303938942;
      break;

    case 4:
      q_points[0] = 0.04144848019938324;
      q_points[1] = 0.2452749143206022;
      q_points[2] = 0.5561654535602751;
      q_points[3] = 0.848982394532986;
      break;

    case 5:
      q_points[0] = 0.02913447215197205;
      q_points[1] = 0.1739772133208974;
      q_points[2] =  0.4117025202849029;
      q_points[3] = 0.6773141745828183;
      q_points[4] = 0.89477136103101;
      break;

    case 6:
      q_points[0] = 0.02163400584411693;
      q_points[1] = 0.1295833911549506;
      q_points[2] = 0.3140204499147661;
      q_points[3] = 0.5386572173517997;
      q_points[4] = 0.7569153373774084;
      q_points[5] = 0.922668851372116;
      break;


    case 7:
      q_points[0] = 0.0167193554082585;
      q_points[1] = 0.100185677915675;
      q_points[2] = 0.2462942462079286;
      q_points[3] = 0.4334634932570557;
      q_points[4] = 0.6323509880476823;
      q_points[5] = 0.81111862674023;
      q_points[6] = 0.940848166743287;
      break;

    case 8:
      q_points[0] = 0.01332024416089244;
      q_points[1] = 0.07975042901389491;
      q_points[2] = 0.1978710293261864;
      q_points[3] =   0.354153994351925;
      q_points[4] =   0.5294585752348643;
      q_points[5] = 0.7018145299391673;
      q_points[6] = 0.849379320441094;
      q_points[7] = 0.953326450056343;
      break;

    case 9:
      q_points[0] = 0.01086933608417545;
      q_points[1] = 0.06498366633800794;
      q_points[2] = 0.1622293980238825;
      q_points[3] = 0.2937499039716641;
      q_points[4] = 0.4466318819056009;
      q_points[5] = 0.6054816627755208;
      q_points[6] = 0.7541101371585467;
      q_points[7] = 0.877265828834263;
      q_points[8] = 0.96225055941096;
      break;

    case 10:
      q_points[0] = 0.00904263096219963;
      q_points[1] = 0.05397126622250072;
      q_points[2] =  0.1353118246392511;
      q_points[3] = 0.2470524162871565;
      q_points[4] = 0.3802125396092744;
      q_points[5] = 0.5237923179723384;
      q_points[6] = 0.6657752055148032;
      q_points[7] = 0.7941904160147613;
      q_points[8] = 0.898161091216429;
      q_points[9] = 0.9688479887196;
      break;


    case 11:
      q_points[0] = 0.007643941174637681;
      q_points[1] = 0.04554182825657903;
      q_points[2] = 0.1145222974551244;
      q_points[3] = 0.2103785812270227;
      q_points[4] = 0.3266955532217897;
      q_points[5] = 0.4554532469286375;
      q_points[6] = 0.5876483563573721;
      q_points[7] = 0.7139638500230458;
      q_points[8] = 0.825453217777127;
      q_points[9] = 0.914193921640008;
      q_points[10] = 0.973860256264123;
      break;

    case 12:
      q_points[0] = 0.006548722279080035;
      q_points[1] = 0.03894680956045022;
      q_points[2] = 0.0981502631060046;
      q_points[3] = 0.1811385815906331;
      q_points[4] = 0.2832200676673157;
      q_points[5] = 0.398434435164983;
      q_points[6] = 0.5199526267791299;
      q_points[7] = 0.6405109167754819;
      q_points[8] = 0.7528650118926111;
      q_points[9] = 0.850240024421055;
      q_points[10] = 0.926749682988251;
      q_points[11] = 0.977756129778486;
      break;

    default:
      Assert(false, ExcNotImplemented());
      break;
    }

  return q_points;
}


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

  switch (n)
    {
    case 1:
      quadrature_weights[0] = -1.0;
      break;
    case 2:
      quadrature_weights[0] = -0.7185393190303845;
      quadrature_weights[1] = -0.2814606809696154;
      break;

    case 3:
      quadrature_weights[0] = -0.5134045522323634;
      quadrature_weights[1] = -0.3919800412014877;
      quadrature_weights[2] = -0.0946154065661483;
      break;

    case 4:
      quadrature_weights[0] =-0.3834640681451353;
      quadrature_weights[1] =-0.3868753177747627;
      quadrature_weights[2] =-0.1904351269501432;
      quadrature_weights[3] =-0.03922548712995894;
      break;

    case 5:
      quadrature_weights[0] =-0.2978934717828955;
      quadrature_weights[1] =-0.3497762265132236;
      quadrature_weights[2] =-0.234488290044052;
      quadrature_weights[3] =-0.0989304595166356;
      quadrature_weights[4] =-0.01891155214319462;
      break;

    case 6:
      quadrature_weights[0] = -0.2387636625785478;
      quadrature_weights[1] = -0.3082865732739458;
      quadrature_weights[2] = -0.2453174265632108;
      quadrature_weights[3] = -0.1420087565664786;
      quadrature_weights[4] = -0.05545462232488041;
      quadrature_weights[5] = -0.01016895869293513;
      break;


    case 7:
      quadrature_weights[0] = -0.1961693894252476;
      quadrature_weights[1] = -0.2703026442472726;
      quadrature_weights[2] = -0.239681873007687;
      quadrature_weights[3] = -0.1657757748104267;
      quadrature_weights[4] = -0.0889432271377365;
      quadrature_weights[5] = -0.03319430435645653;
      quadrature_weights[6] = -0.005932787015162054;
      break;

    case 8:
      quadrature_weights[0] = -0.164416604728002;
      quadrature_weights[1] = -0.2375256100233057;
      quadrature_weights[2] = -0.2268419844319134;
      quadrature_weights[3] = -0.1757540790060772;
      quadrature_weights[4] = -0.1129240302467932;
      quadrature_weights[5] = -0.05787221071771947;
      quadrature_weights[6] = -0.02097907374214317;
      quadrature_weights[7] = -0.003686407104036044;
      break;

    case 9:
      quadrature_weights[0] = -0.1400684387481339;
      quadrature_weights[1] = -0.2097722052010308;
      quadrature_weights[2] = -0.211427149896601;
      quadrature_weights[3] = -0.1771562339380667;
      quadrature_weights[4] = -0.1277992280331758;
      quadrature_weights[5] = -0.07847890261203835;
      quadrature_weights[6] = -0.0390225049841783;
      quadrature_weights[7] = -0.01386729555074604;
      quadrature_weights[8] = -0.002408041036090773;
      break;

    case 10:
      quadrature_weights[0] = -0.12095513195457;
      quadrature_weights[1] = -0.1863635425640733;
      quadrature_weights[2] = -0.1956608732777627;
      quadrature_weights[3] = -0.1735771421828997;
      quadrature_weights[4] = -0.135695672995467;
      quadrature_weights[5] = -0.0936467585378491;
      quadrature_weights[6] = -0.05578772735275126;
      quadrature_weights[7] = -0.02715981089692378;
      quadrature_weights[8] = -0.00951518260454442;
      quadrature_weights[9] = -0.001638157633217673;
      break;


    case 11:
      quadrature_weights[0] = -0.1056522560990997;
      quadrature_weights[1] = -0.1665716806006314;
      quadrature_weights[2] = -0.1805632182877528;
      quadrature_weights[3] = -0.1672787367737502;
      quadrature_weights[4] = -0.1386970574017174;
      quadrature_weights[5] = -0.1038334333650771;
      quadrature_weights[6] = -0.06953669788988512;
      quadrature_weights[7] = -0.04054160079499477;
      quadrature_weights[8] = -0.01943540249522013;
      quadrature_weights[9] = -0.006737429326043388;
      quadrature_weights[10] = -0.001152486965101561;
      break;

    case 12:
      quadrature_weights[0] = -0.09319269144393;
      quadrature_weights[1] = -0.1497518275763289;
      quadrature_weights[2] = -0.166557454364573;
      quadrature_weights[3] = -0.1596335594369941;
      quadrature_weights[4] = -0.1384248318647479;
      quadrature_weights[5] = -0.1100165706360573;
      quadrature_weights[6] = -0.07996182177673273;
      quadrature_weights[7] = -0.0524069547809709;
      quadrature_weights[8] = -0.03007108900074863;
      quadrature_weights[9] = -0.01424924540252916;
      quadrature_weights[10] = -0.004899924710875609;
      quadrature_weights[11] = -0.000834029009809656;
      break;

    default:
      Assert(false, ExcNotImplemented());
      break;
    }

  return quadrature_weights;
}


template <>
QGaussLogR<1>::QGaussLogR(const unsigned int n,
                          const Point<1> origin,
                          const double alpha,
                          const bool factor_out_singularity) :
  Quadrature<1>( ( (origin[0] == 0) || (origin[0] == 1) ) ?
                 (alpha == 1 ? n : 2*n ) : 4*n ),
  fraction( ( (origin[0] == 0) || (origin[0] == 1.) ) ? 1. : origin[0] )
{
  // The three quadrature formulas that make this one up. There are
  // at most two when the origin is one of the extremes, and there is
  // only one if the origin is one of the extremes and alpha is
  // equal to one.
  //
  // If alpha is different from one, then we need a correction which
  // is performed with a standard Gauss quadrature rule on each
  // segment. This is not needed in the standard case where alpha is
  // equal to one and the origin is on one of the extremes. We
  // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
  // only need n quadrature points. In the most difficult one, we
  // need 2*n points for the first segment, and 2*n points for the
  // second segment.
  QGaussLog<1> quad1(n, origin[0] != 0);
  QGaussLog<1> quad2(n);
  QGauss<1> quad(n);

  // Check that the origin is inside 0,1
  Assert( (fraction >= 0) && (fraction <= 1),
          ExcMessage("Origin is outside [0,1]."));

  // Non singular offset. This is the start of non singular quad
  // points.
  unsigned int ns_offset = (fraction == 1) ? n : 2*n;

  for (unsigned int i=0, j=ns_offset; i<n; ++i, ++j)
    {
      // The first i quadrature points are the same as quad1, and
      // are by default singular.
      this->quadrature_points[i] = quad1.point(i)*fraction;
      this->weights[i] = quad1.weight(i)*fraction;

      // We need to scale with -log|fraction*alpha|
      if ( (alpha != 1) || (fraction != 1) )
        {
          this->quadrature_points[j] = quad.point(i)*fraction;
          this->weights[j] = -std::log(alpha/fraction)*quad.weight(i)*fraction;
        }
      // In case we need the second quadrature as well, do it now.
      if (fraction != 1)
        {
          this->quadrature_points[i+n] = quad2.point(i)*(1-fraction)+Point<1>(fraction);
          this->weights[i+n] = quad2.weight(i)*(1-fraction);

          // We need to scale with -log|fraction*alpha|
          this->quadrature_points[j+n] = quad.point(i)*(1-fraction)+Point<1>(fraction);
          this->weights[j+n] = -std::log(alpha/(1-fraction))*quad.weight(i)*(1-fraction);
        }
    }
  if (factor_out_singularity == true)
    for (unsigned int i=0; i<size(); ++i)
      {
        Assert( this->quadrature_points[i] != origin,
                ExcMessage("The singularity cannot be on a Gauss point of the same order!") );
        double denominator = std::log(std::abs( (this->quadrature_points[i]-origin)[0] )/alpha);
        Assert( denominator != 0.0,
                ExcMessage("The quadrature formula you are using does not allow to "
                           "factor out the singularity, which is zero at one point."));
        this->weights[i] /= denominator;
      }
}


template <>
unsigned int QGaussOneOverR<2>::quad_size(const Point<2> singularity,
                                          const unsigned int n)
{
  double eps=1e-8;
  bool on_edge=false;
  bool on_vertex=false;
  for (unsigned int i=0; i<2; ++i)
    if ( ( std::abs(singularity[i]  ) < eps ) ||
         ( std::abs(singularity[i]-1) < eps ) )
      on_edge = true;
  if (on_edge && (std::abs( (singularity-Point<2>(.5, .5)).norm_square()-.5)
                  < eps) )
    on_vertex = true;
  if (on_vertex) return (2*n*n);
  if (on_edge) return (4*n*n);
  return (8*n*n);
}

template <>
QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
                                  const Point<2> singularity,
                                  const bool factor_out_singularity) :
  Quadrature<2>(quad_size(singularity, n))
{
  // We treat all the cases in the
  // same way. Split the element in 4
  // pieces, measure the area, if
  // it's relevant, add the
  // quadrature connected to that
  // singularity.
  std::vector<QGaussOneOverR<2> > quads;
  std::vector<Point<2> > origins;
  // Id of the corner with a
  // singularity
  quads.emplace_back(n, 3, factor_out_singularity);
  quads.emplace_back(n, 2, factor_out_singularity);
  quads.emplace_back(n, 1, factor_out_singularity);
  quads.emplace_back(n, 0, factor_out_singularity);

  origins.emplace_back(0., 0.);
  origins.emplace_back(singularity[0], 0.);
  origins.emplace_back(0., singularity[1]);
  origins.push_back(singularity);

  // Lexicographical ordering.

  double eps = 1e-8;
  unsigned int q_id = 0; // Current quad point index.
  Tensor<1,2> dist;

  for (unsigned int box=0; box<4; ++box)
    {
      dist = (singularity-GeometryInfo<2>::unit_cell_vertex(box));
      dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
      double area = dist[0]*dist[1];
      if (area > eps)
        for (unsigned int q=0; q<quads[box].size(); ++q, ++q_id)
          {
            const Point<2> &qp = quads[box].point(q);
            this->quadrature_points[q_id] =
              origins[box]+
              Point<2>(dist[0]*qp[0], dist[1]*qp[1]);
            this->weights[q_id] = quads[box].weight(q)*area;
          }
    }
}


template <>
QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
                                  const unsigned int vertex_index,
                                  const bool factor_out_singularity) :
  Quadrature<2>(2*n *n)
{
  // This version of the constructor
  // works only for the 4
  // vertices. If you need a more
  // general one, you should use the
  // one with the Point<2> in the
  // constructor.
  Assert(vertex_index <4, ExcIndexRange(vertex_index, 0, 4));

  // Start with the gauss quadrature formula on the (u,v) reference
  // element.
  QGauss<2> gauss(n);

  Assert(gauss.size() == n*n, ExcInternalError());

  // For the moment we only implemented this for the vertices of a
  // quadrilateral. We are planning to do this also for the support
  // points of arbitrary FE_Q elements, to allow the use of this
  // class in boundary element programs with higher order mappings.
  Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));

  // We create only the first one. All other pieces are rotation of
  // this one.
  // In this case the transformation is
  //
  // (x,y) = (u, u tan(pi/4 v))
  //
  // with Jacobian
  //
  // J = pi/4 R / cos(pi/4 v)
  //
  // And we get rid of R to take into account the singularity,
  // unless specified differently in the constructor.
  std::vector<Point<2> >      &ps = this->quadrature_points;
  std::vector<double>         &ws = this->weights;
  double pi4 = numbers::PI/4;

  for (unsigned int q=0; q<gauss.size(); ++q)
    {
      const Point<2> &gp = gauss.point(q);
      ps[q][0] = gp[0];
      ps[q][1] = gp[0] * std::tan(pi4*gp[1]);
      ws[q]    = gauss.weight(q)*pi4/std::cos(pi4 *gp[1]);
      if (factor_out_singularity)
        ws[q] *= (ps[q]-GeometryInfo<2>::unit_cell_vertex(0)).norm();
      // The other half of the quadrilateral is symmetric with
      // respect to xy plane.
      ws[gauss.size()+q]    = ws[q];
      ps[gauss.size()+q][0] = ps[q][1];
      ps[gauss.size()+q][1] = ps[q][0];
    }

  // Now we distribute these vertices in the correct manner
  double theta = 0;
  switch (vertex_index)
    {
    case 0:
      theta = 0;
      break;
    case 1:
      //
      theta = numbers::PI/2;
      break;
    case 2:
      theta = -numbers::PI/2;
      break;
    case 3:
      theta = numbers::PI;
      break;
    }

  double R00 =  std::cos(theta), R01 = -std::sin(theta);
  double R10 =  std::sin(theta), R11 =  std::cos(theta);

  if (vertex_index != 0)
    for (unsigned int q=0; q<size(); ++q)
      {
        double x = ps[q][0]-.5,  y = ps[q][1]-.5;

        ps[q][0] = R00*x + R01*y + .5;
        ps[q][1] = R10*x + R11*y + .5;
      }
}


template <int dim>
QSorted<dim>::QSorted(const Quadrature<dim> &quad) :
  Quadrature<dim>(quad)
{
  std::vector<unsigned int> permutation(quad.size());
  for (unsigned int i=0; i<quad.size(); ++i)
    permutation[i] = i;

  std::sort(permutation.begin(),
            permutation.end(),
            std::bind(&QSorted<dim>::compare_weights,
                      std::ref(*this),
                      std::placeholders::_1,
                      std::placeholders::_2));

  // At this point, the variable is_tensor_product_flag is set
  // to the respective value of the given Quadrature in the base
  // class copy constructor.
  // We only call a quadrature formula 'tensor product'
  // if the quadrature points are also sorted lexicographically.
  // In particular, any reordering destroys that property
  // and we might need to modify the variable accordingly.
  for (unsigned int i=0; i<quad.size(); ++i)
    {
      this->weights[i]           = quad.weight(permutation[i]);
      this->quadrature_points[i] = quad.point(permutation[i]);
      if (permutation[i] != i)
        this->is_tensor_product_flag = false;
    }
}


template <int dim>
bool QSorted<dim>::compare_weights(const unsigned int a,
                                   const unsigned int b) const
{
  return (this->weights[a] < this->weights[b]);
}


// construct the quadrature formulae in higher dimensions by
// tensor product of lower dimensions

template <int dim>
QGauss<dim>::QGauss (const unsigned int n)
  :  Quadrature<dim> (QGauss<dim-1>(n), QGauss<1>(n))
{}



template <int dim>
QGaussLobatto<dim>::QGaussLobatto (const unsigned int n)
  :  Quadrature<dim> (QGaussLobatto<dim-1>(n), QGaussLobatto<1>(n))
{}



template <int dim>
QMidpoint<dim>::QMidpoint ()
  :
  Quadrature<dim> (QMidpoint<dim-1>(), QMidpoint<1>())
{}



template <int dim>
QTrapez<dim>::QTrapez ()
  :
  Quadrature<dim> (QTrapez<dim-1>(), QTrapez<1>())
{}



template <int dim>
QSimpson<dim>::QSimpson ()
  :
  Quadrature<dim> (QSimpson<dim-1>(), QSimpson<1>())
{}



template <int dim>
QMilne<dim>::QMilne ()
  :
  Quadrature<dim> (QMilne<dim-1>(), QMilne<1>())
{}


template <int dim>
QWeddle<dim>::QWeddle ()
  :
  Quadrature<dim> (QWeddle<dim-1>(), QWeddle<1>())
{}

template <int dim>
QTelles<dim>::QTelles (
  const Quadrature<1> &base_quad, const Point<dim> &singularity)
  :
  // We need the explicit implementation if dim == 1. If dim > 1 we use the
  // former implementation and apply a tensorial product to obtain the higher
  // dimensions.
  Quadrature<dim>(
    dim == 2 ?
    QAnisotropic<dim>(
      QTelles<1>(base_quad, Point<1>(singularity[0])),
      QTelles<1>(base_quad, Point<1>(singularity[1]))) :
    dim == 3 ?
    QAnisotropic<dim>(
      QTelles<1>(base_quad, Point<1>(singularity[0])),
      QTelles<1>(base_quad, Point<1>(singularity[1])),
      QTelles<1>(base_quad, Point<1>(singularity[2]))) :
    Quadrature<dim>())
{
}

template <int dim>
QTelles<dim>::QTelles (
  const unsigned int n, const Point<dim> &singularity)
  :
  // In this case we map the standard Gauss Legendre formula using the given
  // singularity point coordinates.
  Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
{}



template <>
QTelles<1>::QTelles (
  const Quadrature<1> &base_quad, const Point<1> &singularity)
  :
  // We explicitly implement the Telles' variable change if dim == 1.
  Quadrature<1>(base_quad)
{
  // We define all the constants to be used in the implementation of
  // Telles' rule
  const double eta_bar = singularity[0] * 2. - 1.;
  const double eta_star = eta_bar * eta_bar - 1.;
  double gamma_bar;

  std::vector<Point<1> > quadrature_points_dummy(quadrature_points.size());
  std::vector<double> weights_dummy(weights.size());
  unsigned int cont = 0;
  const double tol = 1e-10;
  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
    {
      if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
        {
          quadrature_points_dummy[d-cont] = quadrature_points[d];
          weights_dummy[d-cont] = weights[d];
        }
      else
        {
          // We need to remove the singularity point from the quadrature point
          // list. To do so we use the variable cont.
          cont = 1;
        }

    }
  if (cont == 1)
    {
      quadrature_points.resize(quadrature_points_dummy.size()-1);
      weights.resize(weights_dummy.size()-1);
      for (unsigned int d = 0; d < quadrature_points.size(); ++d)
        {
          quadrature_points[d] = quadrature_points_dummy[d];
          weights[d] = weights_dummy[d];
        }
    }
  // We need to check if the singularity is at the boundary of the interval.
  if (std::abs(eta_star) <= tol)
    {
      gamma_bar = std::pow((eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0)
                  + std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0)
                  + eta_bar;
    }
  else
    {
      gamma_bar = (eta_bar * eta_star + std::abs(eta_star))/std::abs(eta_bar * eta_star + std::abs(eta_star))*
                  std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0)
                  + (eta_bar * eta_star - std::abs(eta_star))/std::abs(eta_bar * eta_star - std::abs(eta_star))*
                  std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0)
                  + eta_bar;
    }
  for (unsigned int q = 0; q < quadrature_points.size(); ++q)
    {
      double gamma = quadrature_points[q][0] * 2 - 1;
      double eta = (std::pow(gamma - gamma_bar, 3.0)
                    + gamma_bar * (gamma_bar * gamma_bar + 3))
                   / (1 + 3 * gamma_bar * gamma_bar);

      double J = 3 * ((gamma - gamma_bar) *(gamma - gamma_bar))
                 / (1 + 3 * gamma_bar * gamma_bar);

      quadrature_points[q][0] = (eta + 1) / 2.0;
      weights[q] = J * weights[q];

    }
}

namespace internal
{
  namespace QGaussChebyshev
  {
    std::vector<double>
    get_quadrature_points(const unsigned int n)
    {

      std::vector<double> points(n);
      // n point quadrature: index from 0 to n-1
      for (unsigned short i=0; i<n; ++i)
        // would be cos((2i+1)Pi/(2N+2))
        // put + Pi so we start from the smallest point
        // then map from [-1,1] to [0,1]
        points[i] = 1./2.*(1.+std::cos(numbers::PI*(1.+double(2*i+1)/double(2*(n-1)+2))));

      return points;
    }



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

      std::vector<double> weights(n);

      for (unsigned short i=0; i<n; ++i)
        {
          // same weights as on [-1,1]
          weights[i] = numbers::PI/double(n);
        }

      return weights;

    }
  }
}


template <>
QGaussChebyshev<1>::QGaussChebyshev(const unsigned int n)
  :
  Quadrature<1> (n)
{

  Assert(n>0,ExcMessage("Need at least one point for the quadrature rule"));
  std::vector<double> p=internal::QGaussChebyshev::get_quadrature_points(n);
  std::vector<double> w=internal::QGaussChebyshev::get_quadrature_weights(n);

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(p[i]);
      this->weights[i]           = w[i];
    }
}


template <int dim>
QGaussChebyshev<dim>::QGaussChebyshev (const unsigned int n)
  :
  Quadrature<dim> (QGaussChebyshev<1>(n))
{}


namespace internal
{
  namespace QGaussRadauChebyshev
  {
    // Computes the points of the quadrature formula.
    std::vector<double>
    get_quadrature_points(const unsigned int n,
                          ::QGaussRadauChebyshev<1>::EndPoint ep)
    {

      std::vector<double> points(n);
      // n point quadrature: index from 0 to n-1
      for (unsigned short i=0; i<n; ++i)
        // would be -cos(2i Pi/(2N+1))
        // put + Pi so we start from the smallest point
        // then map from [-1,1] to [0,1]
        switch (ep)
          {
          case ::QGaussRadauChebyshev<1>::left:
          {
            points[i] = 1./2.*(1.-std::cos(numbers::PI*(1+2*double(i)/(2*double(n-1)+1.))));
            break;
          }

          case ::QGaussRadauChebyshev<1>::right:
          {
            points[i] = 1./2.*(1.-std::cos(numbers::PI*(2*double(n-1-i)/(2*double(n-1)+1.))));
            break;
          }

          default:
            Assert (false, ExcMessage ("This constructor can only be called with either "
                                       "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
                                       "second argument."));
          }

      return points;
    }



    // Computes the weights of the quadrature formula.
    std::vector<double>
    get_quadrature_weights(const unsigned int n,
                           ::QGaussRadauChebyshev<1>::EndPoint ep)
    {

      std::vector<double> weights(n);

      for (unsigned short i=0; i<n; ++i)
        {
          // same weights as on [-1,1]
          weights[i] = 2.*numbers::PI/double(2*(n-1)+1.);
          if (ep==::QGaussRadauChebyshev<1>::left && i==0)
            weights[i] /= 2.;
          else if (ep==::QGaussRadauChebyshev<1>::right && i==(n-1))
            weights[i] /= 2.;
        }

      return weights;
    }
  }
}


template <>
QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n,
                                              EndPoint ep)
  :
  Quadrature<1> (n),
  ep (ep)
{

  Assert(n>0,ExcMessage("Need at least one point for quadrature rules"));
  std::vector<double> p=internal::QGaussRadauChebyshev::get_quadrature_points(n,ep);
  std::vector<double> w=internal::QGaussRadauChebyshev::get_quadrature_weights(n,ep);

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(p[i]);
      this->weights[i]           = w[i];
    }
}


template <int dim>
QGaussRadauChebyshev<dim>::QGaussRadauChebyshev (const unsigned int n,
                                                 EndPoint ep)
  :
  Quadrature<dim> (QGaussRadauChebyshev<1>(n,static_cast<QGaussRadauChebyshev<1>::EndPoint>(ep))),
  ep (ep)
{}



namespace internal
{
  namespace QGaussLobattoChebyshev
  {
    // Computes the points of the quadrature formula.
    std::vector<double>
    get_quadrature_points(const unsigned int n)
    {

      std::vector<double> points(n);
      // n point quadrature: index from 0 to n-1
      for (unsigned short i=0; i<n; ++i)
        // would be cos(i Pi/N)
        // put + Pi so we start from the smallest point
        // then map from [-1,1] to [0,1]
        points[i] = 1./2.*(1.+std::cos(numbers::PI*(1+double(i)/double(n-1))));

      return points;
    }

    // Computes the weights of the quadrature formula.
    std::vector<double>
    get_quadrature_weights(const unsigned int n)
    {

      std::vector<double> weights(n);

      for (unsigned short i=0; i<n; ++i)
        {
          // same weights as on [-1,1]
          weights[i] = numbers::PI/double((n-1));
          if (i==0 || i==(n-1))
            weights[i] /= 2.;
        }

      return weights;
    }
  }
}



template <>
QGaussLobattoChebyshev<1>::QGaussLobattoChebyshev(const unsigned int n)
  :
  Quadrature<1> (n)
{

  Assert(n>1,ExcMessage("Need at least two points for Gauss-Lobatto quadrature rule"));
  std::vector<double> p=internal::QGaussLobattoChebyshev::get_quadrature_points(n);
  std::vector<double> w=internal::QGaussLobattoChebyshev::get_quadrature_weights(n);

  for (unsigned int i=0; i<this->size(); ++i)
    {
      this->quadrature_points[i] = Point<1>(p[i]);
      this->weights[i]           = w[i];
    }
}


template <int dim>
QGaussLobattoChebyshev<dim>::QGaussLobattoChebyshev (const unsigned int n)
  :
  Quadrature<dim> (QGaussLobattoChebyshev<1>(n))
{}



template<int dim>
QSimplex<dim>::QSimplex(const Quadrature<dim> &quad)
{
  std::vector<Point<dim> > qpoints;
  std::vector<double > weights;

  for (unsigned int i=0; i<quad.size(); ++i)
    {
      double r=0;
      for (unsigned int d=0; d<dim; ++d)
        r += quad.point(i)[d];
      if (r <= 1+1e-10)
        {
          this->quadrature_points.push_back(quad.point(i));
          this->weights.push_back(quad.weight(i));
        }
    }
}



template<int dim>
Quadrature<dim>
QSimplex<dim>::compute_affine_transformation(const std::array<Point<dim>, dim+1>& vertices) const
{
  Tensor<2,dim> B;
  for (unsigned int d=0; d<dim; ++d)
    B[d] = vertices[d+1]-vertices[0];

  B = transpose(B);
  const double J = std::abs(determinant(B));

  // if the determinant is zero, we return an empty quadrature
  if (J < 1e-12)
    return Quadrature<dim>();

  std::vector<Point<dim> > qp(this->size());
  std::vector<double> w(this->size());

  for (unsigned int i=0; i<this->size(); ++i)
    {
      qp[i] = Point<dim>(vertices[0]+B*this->point(i));
      w[i] = J*this->weight(i);
    }

  return Quadrature<dim>(qp, w);
}



QTrianglePolar::QTrianglePolar(const Quadrature<1> &radial_quadrature,
                               const Quadrature<1> &angular_quadrature) :
  QSimplex<2>(Quadrature<2>())
{
  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
  this->quadrature_points.resize(base.size());
  this->weights.resize(base.size());
  for (unsigned int i=0; i<base.size(); ++i)
    {
      const auto q = base.point(i);
      const auto w = base.weight(i);

      const auto xhat = q[0];
      const auto yhat = q[1];

      const double t = numbers::PI_2*yhat;
      const double pi = numbers::PI;
      const double st = std::sin(t);
      const double ct = std::cos(t);
      const double r = xhat/(st+ct);

      const double J = pi*xhat/(2*(std::sin(pi*yhat) + 1));

      this->quadrature_points[i] = Point<2>(r*ct, r*st);
      this->weights[i] = w*J;
    }
}



QTrianglePolar::QTrianglePolar(const unsigned int &n)
  :QTrianglePolar(QGauss<1>(n),QGauss<1>(n))
{}



QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
               const Quadrature<1> &angular_quadrature,
               const double beta) :
  QSimplex<2>(Quadrature<2>())
{
  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
  this->quadrature_points.resize(base.size());
  this->weights.resize(base.size());
  for (unsigned int i=0; i<base.size(); ++i)
    {
      const auto q = base.point(i);
      const auto w = base.weight(i);

      const auto xhat = q[0];
      const auto yhat = q[1];

      const double x = std::pow(xhat, beta)*(1-yhat);
      const double y = std::pow(xhat, beta)*yhat;

      const double J = beta * std::pow(xhat, 2.*beta-1.);

      this->quadrature_points[i] = Point<2>(x,y);
      this->weights[i] = w*J;
    }
}



QDuffy::QDuffy(const unsigned int n, const double beta)
  :QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
{}



template<int dim>
QSplit<dim>::QSplit(const QSimplex<dim> &base,
                    const Point<dim> &split_point)
{
  AssertThrow(GeometryInfo<dim>::is_inside_unit_cell(split_point, 1e-12),
              ExcMessage("The split point should be inside the unit reference cell."));

  std::array<Point<dim>, dim+1> vertices;
  vertices[0] = split_point;

  // Make a simplex from the split_point and the first dim vertices of each
  // face. In dimension three, we need to split the face in two triangles, so
  // we use once the first dim vertices of each face, and the second time the
  // the dim vertices of each face starting from 1.
  for (unsigned int f=0; f< GeometryInfo<dim>::faces_per_cell; ++f)
    for (unsigned int start=0; start < (dim >2 ? 2 : 1); ++start)
      {
        for (unsigned int i=0; i<dim; ++i)
          vertices[i+1] =
            GeometryInfo<dim>::unit_cell_vertex(
              GeometryInfo<dim>::face_to_cell_vertices(f, i+start)
            );
        const auto quad = base.compute_affine_transformation(vertices);
        if (quad.size())
          {
            this->quadrature_points.insert(
              this->quadrature_points.end(),
              quad.get_points().begin(),
              quad.get_points().end());
            this->weights.insert(
              this->weights.end(),
              quad.get_weights().begin(),
              quad.get_weights().end());
          }
      }
}



// explicit specialization
// note that 1d formulae are specialized by implementation above
template class QGauss<2>;
template class QGaussLobatto<2>;
template class QMidpoint<2>;
template class QTrapez<2>;
template class QSimpson<2>;
template class QMilne<2>;
template class QWeddle<2>;

template class QGauss<3>;
template class QGaussLobatto<3>;
template class QMidpoint<3>;
template class QTrapez<3>;
template class QSimpson<3>;
template class QMilne<3>;
template class QWeddle<3>;

template class QSorted<1>;
template class QSorted<2>;
template class QSorted<3>;

template class QTelles<1> ;
template class QTelles<2> ;
template class QTelles<3> ;

template class QGaussChebyshev<1>;
template class QGaussChebyshev<2>;
template class QGaussChebyshev<3>;

template class QGaussRadauChebyshev<1>;
template class QGaussRadauChebyshev<2>;
template class QGaussRadauChebyshev<3>;

template class QGaussLobattoChebyshev<1>;
template class QGaussLobattoChebyshev<2>;
template class QGaussLobattoChebyshev<3>;

template class QSimplex<1>;
template class QSimplex<2>;
template class QSimplex<3>;

template class QSplit<1>;
template class QSplit<2>;
template class QSplit<3>;

DEAL_II_NAMESPACE_CLOSE