From cffc510c9b385828ed4e4b3e0880aad78b1c5bf0 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 19 Sep 2003 15:55:37 +0000 Subject: [PATCH] Fix computation of Gauss points. git-svn-id: https://svn.dealii.org/trunk@7989 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/source/quadrature_lib.cc | 29 ++++++++++++++++++++------- deal.II/doc/news/c-4-0.html | 9 +++++++++ 2 files changed, 31 insertions(+), 7 deletions(-) diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 86d92a06e2..f4709f74b4 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -14,6 +14,7 @@ #include #include +#include // please note: for a given dimension, we need the quadrature formulae @@ -34,15 +35,29 @@ QGauss<1>::QGauss (const unsigned int n) Quadrature<1> (n) { const unsigned int m = (n+1)/2; - long double z; - long double pp; - long double p1, p2, p3; - for (unsigned int i=1;i<=m;++i) + // tolerance for the Newton + // iteration below. we need to make + // it adaptive since on some + // machines (for example PowerPC) + // long double is the same as + // double -- in that case we can + // only get to a certain multiple + // of the accuracy of double there, + // while on other machines we'd + // like to go further down + const long double tolerance + = std::max (static_cast(std::numeric_limits::epsilon() / 100), + static_cast(std::numeric_limits::epsilon() * 5)); + + for (unsigned int i=1; i<=m; ++i) { - z = std::cos(deal_II_numbers::PI * (i-.25)/(n+.5)); + long double z = std::cos(deal_II_numbers::PI * (i-.25)/(n+.5)); - // Newton-iteration + long double pp; + long double p1, p2, p3; + + // Newton iteration do { // compute L_n (z) @@ -57,7 +72,7 @@ QGauss<1>::QGauss (const unsigned int n) pp = n*(z*p1-p2)/(z*z-1); z = z-p1/pp; } - while (abs(p1/pp) > 1.e-19); + while (abs(p1/pp) > tolerance); double x = .5*z; this->quadrature_points[i-1] = Point<1>(.5-x); diff --git a/deal.II/doc/news/c-4-0.html b/deal.II/doc/news/c-4-0.html index b63fcd672b..c3aa108f66 100644 --- a/deal.II/doc/news/c-4-0.html +++ b/deal.II/doc/news/c-4-0.html @@ -98,6 +98,15 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

base

    +
  1. Fixed: The constructor of the QGauss class + computed positions and weights of quadrature points in long double accuracy. + However, on machines where long double is the same as double, it therefore + never reached the requested accuracy, in effect leading to an infinite loop. + This is now fixed. +
    + (WB 2002/09/19) +

    +
  2. New: The Function class now exports the value of its template argument through the static member variable dimension. -- 2.39.5