From 8d394fafc54f93500bd864d1eeeab109926cd352 Mon Sep 17 00:00:00 2001 From: prill Date: Thu, 2 Nov 2006 14:40:04 +0000 Subject: [PATCH] Implemented Legendre-Gauss-Lobatto quadrature. git-svn-id: https://svn.dealii.org/trunk@14140 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature_lib.h | 73 ++++++++++- deal.II/base/source/quadrature_lib.cc | 137 ++++++++++++++++++++- deal.II/doc/news/changes.html | 5 + 3 files changed, 207 insertions(+), 8 deletions(-) diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 6aa72a8e12..c34825041d 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -49,19 +49,27 @@ class QGauss : public Quadrature * The Gauss-Lobatto quadrature rule. * * This modification of the Gauss quadrature uses the two interval end - * points as well. Being exact for polynomials of degree 2n-2, - * this formula is suboptimal by one degree. + * points as well. Being exact for polynomials of degree 2n-3, + * this formula is suboptimal by two degrees. * - * The quadrature points are interval end points plus the zeroes of + * The quadrature points are interval end points plus the roots of * the derivative of the Legendre polynomial Pn-1 of * degree n-1. The quadrature weights are * 2/(n(n-1)(Pn-1(xi)2). * - * @note The quadrature weights are not implemented yet. + * Note: This implementation has not yet been optimized concerning + * numerical stability and efficiency. It can be easily adapted + * to the general case of Gauss-Lobatto-Jacobi-Bouzitat quadrature + * with arbitrary parameters alpha, beta, of which + * the Gauss-Lobatto-Legendre quadrature (alpha = beta = 0) + * is a special case. * - * @sa http://en.wikipedia.org/wiki/Handbook_of_Mathematical_Functions + * @sa http://en.wikipedia.org/wiki/Handbook_of_Mathematical_Functions + * @sa Karniadakis, G.E. and Sherwin, S.J.: + * Spectral/hp element methods for computational fluid dynamics. + * Oxford: Oxford University Press, 2005 * - * @author Guido Kanschat, 2005, 2006 + * @author Guido Kanschat, 2005, 2006; F. Prill, 2006 */ template class QGaussLobatto : public Quadrature @@ -73,6 +81,47 @@ class QGaussLobatto : public Quadrature * (in each space direction). */ QGaussLobatto(const unsigned int n); + + protected: + /** + * Compute Legendre-Gauss-Lobatto quadrature + * points in the interval [-1, +1]. + * @param q number of points. + * @return vector containing nodes. + */ + std::vector compute_quadrature_points(const unsigned int q, + const int alpha, + const int beta) const; + + /** + * Compute Legendre-Gauss-Lobatto quadrature + * weights. + * @param x quadrature points. + * @return vector containing weights. + */ + std::vector compute_quadrature_weights(std::vector& x, + const int alpha, + const int beta) const; + + /** + * Evaluate a Jacobi polynomial + * \f$ P^{\alpha, \beta}_n(x) \f$. + * Note: The Jacobi polynomials are + * not orthonormal and defined on + * the interval [-1, +1]. + * @param x point of evaluation. + */ + double JacobiP(const double x, + const int alpha, + const int beta, + const unsigned int n) const; + + /** + * Evaluate the Gamma function + * \f[ \Gamma(n) = (n-1)! \f]. + * @param n point of evaluation (integer). + */ + unsigned int gamma(const unsigned int n) const; }; @@ -250,6 +299,18 @@ class QWeddle : public Quadrature template <> QGauss<1>::QGauss (const unsigned int n); template <> QGaussLobatto<1>::QGaussLobatto (const unsigned int n); +template <> +std::vector QGaussLobatto<1>:: +compute_quadrature_points(const unsigned int, const int, const int) const; +template <> +std::vector QGaussLobatto<1>:: +compute_quadrature_weights(std::vector&, const int, const int) const; +template <> +double QGaussLobatto<1>:: +JacobiP(const double, const int, const int, const unsigned int) const; +template <> +unsigned int QGaussLobatto<1>:: +QGaussLobatto<1>::gamma(const unsigned int n) const; template <> QGauss2<1>::QGauss2 (); template <> QGauss3<1>::QGauss3 (); template <> QGauss4<1>::QGauss4 (); diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 223765e532..49a41311dc 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -156,13 +156,139 @@ QGaussLobatto<1>::QGaussLobatto (unsigned int n) : Quadrature<1> (n) { - for (unsigned int k=0;k points = compute_quadrature_points(n, 1, 1); + std::vector w = compute_quadrature_weights(points, 0, 0); + + // scale points to the interval + // [0.0, 1.0]: + for (unsigned int i=0; iquadrature_points[i] = Point<1>(0.5 + 0.5*points[i]); + this->weights[i] = 0.5*w[i]; + } +} + + + +template <> +std::vector QGaussLobatto<1>:: +compute_quadrature_points(const unsigned int q, + const int alpha, + const int beta) const +{ + const unsigned int m = q-2; // no. of inner points + std::vector x(m); + + // compute quadrature points with + // a Newton algorithm. + const double epsilon = 1.e-10; // tolerance + + // we take the zeros of the Chebyshev + // polynomial (alpha=beta=-0.5) as + // initial values: + for (unsigned int i=0; i0) + r = (r + x[k-1])/2; + + do + { + s = 1.; + for (unsigned int i=0; i= epsilon); + + x[k] = r; + } // for + + // add boundary points: + x.insert(x.begin(), -1.); + x.push_back(+1.); + + return x; +} + + + +template <> +std::vector QGaussLobatto<1>:: +compute_quadrature_weights(std::vector& x, + const int alpha, + const int beta) const +{ + const unsigned int q = x.size(); + std::vector w(q); + double s = 0; + + 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; iquadrature_points[k] = Point<1>(.5-.5*std::cos(M_PI*k/(n-1))); + s = JacobiP(x[i], alpha, beta, q-1); + w[i] = factor/(s*s); } + w[0] *= (beta + 1); + w[q-1] *= (alpha + 1); + + return w; +} + + + +template <> +double QGaussLobatto<1>::JacobiP(const double x, + const int alpha, + const int beta, + const unsigned int n) const +{ + // the Jacobi polynomial is evaluated + // using a recursion formula. + std::vector p(n+1); + double v, a1, a2, a3, a4; + + // initial values P_0(x), P_1(x): + p[0] = 1.0; + if (n==0) return p[0]; + p[1] = ((alpha+beta+2)*x + (alpha-beta))/2; + if (n==1) return p[1]; + + for (unsigned int i=1; i<=(n-1); ++i) + { + v = 2*i + alpha + beta; + a1 = 2*(i+1)*(i + alpha + beta + 1)*v; + a2 = (v + 1)*(alpha*alpha - beta*beta); + a3 = v*(v + 1)*(v + 2); + a4 = 2*(i+alpha)*(i+beta)*(v + 2); + + p[i+1] = 1/a1*( (a2 + a3*x)*p[i] - a4*p[i-1]); + } // for + return p[n]; } + +template <> +unsigned int QGaussLobatto<1>::gamma(const unsigned int n) const +{ + unsigned int result = n - 1; + for (int i=n-2; i>1; --i) + result *= i; + return result; +} + + + template <> QGauss2<1>::QGauss2 () : @@ -473,6 +599,13 @@ QGauss::QGauss (const unsigned int n) +template +QGaussLobatto::QGaussLobatto (const unsigned int n) + : Quadrature (QGaussLobatto(n), QGaussLobatto<1>(n)) +{} + + + template QGauss2::QGauss2 () : diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index e041f93bc2..e44bf219de 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -390,6 +390,11 @@ inconvenience this causes.

base

    +
  1. Extended: The QGaussLobatto quadrature rule computes + Legendre-Gauss-Lobatto nodes and quadrature weights. +
    + (F. Prill 2006/11/02) +

  2. Extended: The contract function family contracts two tensors. There is now a new version of this function, that contracts a tensor of rank three with a second one of rank two over given indices -- 2.39.5