From f351817f66a77b2841792ea2520c4cdc4927c6ad Mon Sep 17 00:00:00 2001 From: prill Date: Wed, 22 Nov 2006 10:52:25 +0000 Subject: [PATCH] Internal changes in Lobatto quadrature: long double data type. git-svn-id: https://svn.dealii.org/trunk@14209 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature_lib.h | 24 ++++---- deal.II/base/source/quadrature_lib.cc | 66 +++++++++++++--------- 2 files changed, 50 insertions(+), 40 deletions(-) diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index cbf86eaaf2..fab71619ca 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -95,7 +95,7 @@ class QGaussLobatto : public Quadrature * * @return vector containing nodes. */ - std::vector + std::vector compute_quadrature_points (const unsigned int q, const int alpha, const int beta) const; @@ -109,8 +109,8 @@ class QGaussLobatto : public Quadrature * @p x denotes the quadrature points. * @return vector containing weights. */ - std::vector - compute_quadrature_weights (const std::vector &x, + std::vector + compute_quadrature_weights (const std::vector &x, const int alpha, const int beta) const; @@ -124,10 +124,10 @@ class QGaussLobatto : public Quadrature * the interval $[-1, +1]$. * @p x is the point of evaluation. */ - double JacobiP(const double x, - const int alpha, - const int beta, - const unsigned int n) const; + long double JacobiP(const long double x, + const int alpha, + const int beta, + const unsigned int n) const; /** * Evaluate the Gamma function @@ -313,14 +313,14 @@ 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>:: +std::vector QGaussLobatto<1>:: compute_quadrature_points(const unsigned int, const int, const int) const; template <> -std::vector QGaussLobatto<1>:: -compute_quadrature_weights(const std::vector&, const int, const int) const; +std::vector QGaussLobatto<1>:: +compute_quadrature_weights(const std::vector&, const int, const int) const; template <> -double QGaussLobatto<1>:: -JacobiP(const double, const int, const int, const unsigned int) const; +long double QGaussLobatto<1>:: +JacobiP(const long double, const int, const int, const unsigned int) const; template <> unsigned int QGaussLobatto<1>:: QGaussLobatto<1>::gamma(const unsigned int n) const; diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index b556166fc4..541af2ceed 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -156,14 +156,16 @@ QGaussLobatto<1>::QGaussLobatto (unsigned int n) : Quadrature<1> (n) { - std::vector points = compute_quadrature_points(n, 1, 1); - std::vector w = compute_quadrature_weights(points, 0, 0); + Assert (n >= 2, ExcNotImplemented()); + + std::vector 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->quadrature_points[i] = Point<1>(0.5 + 0.5*static_cast(points[i])); this->weights[i] = 0.5*w[i]; } } @@ -171,25 +173,33 @@ QGaussLobatto<1>::QGaussLobatto (unsigned int n) template <> -std::vector QGaussLobatto<1>:: +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); + std::vector x(m); // compute quadrature points with // a Newton algorithm. - const double epsilon = 1.e-10; // tolerance - + + // set tolerance +#ifdef HAVE_STD_NUMERIC_LIMITS + const long double + epsilon = static_cast(std::numeric_limits::epsilon()); +#else + const long double + epsilon = 1.e-19L; +#endif + // we take the zeros of the Chebyshev // polynomial (alpha=beta=-0.5) as // initial values: for (unsigned int i=0; i= epsilon); @@ -213,8 +223,8 @@ compute_quadrature_points(const unsigned int q, } // for // add boundary points: - x.insert(x.begin(), -1.); - x.push_back(+1.); + x.insert(x.begin(), -1.L); + x.push_back(+1.L); return x; } @@ -222,19 +232,19 @@ compute_quadrature_points(const unsigned int q, template <> -std::vector QGaussLobatto<1>:: -compute_quadrature_weights(const std::vector &x, +std::vector QGaussLobatto<1>:: +compute_quadrature_weights(const std::vector &x, const int alpha, const int beta) const { const unsigned int q = x.size(); - std::vector w(q); - double s = 0; + std::vector w(q); + long double s = 0.L; - const double factor = std::pow(2., alpha+beta+1) * - gamma(alpha+q) * - gamma(beta+q) / - ((q-1)*gamma(q)*gamma(alpha+beta+q+1)); + 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 &x, template <> -double QGaussLobatto<1>::JacobiP(const double x, - const int alpha, - const int beta, - const unsigned int n) const +long double QGaussLobatto<1>::JacobiP(const long 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; + std::vector p(n+1); + int v, a1, a2, a3, a4; // initial values P_0(x), P_1(x): - p[0] = 1.0; + p[0] = 1.0L; if (n==0) return p[0]; p[1] = ((alpha+beta+2)*x + (alpha-beta))/2; if (n==1) return p[1]; @@ -273,7 +283,7 @@ double QGaussLobatto<1>::JacobiP(const double x, 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]); + p[i+1] = static_cast( (a2 + a3*x)*p[i] - a4*p[i-1])/a1; } // for return p[n]; } -- 2.39.5