]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix computation of Gauss points.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Sep 2003 15:55:37 +0000 (15:55 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 19 Sep 2003 15:55:37 +0000 (15:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@7989 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/quadrature_lib.cc
deal.II/doc/news/c-4-0.html

index 86d92a06e2532dfad8cd32e320d1d907b637b151..f4709f74b42eae57a0b286f9fd221db71b470e4b 100644 (file)
@@ -14,6 +14,7 @@
 
 #include <base/quadrature_lib.h>
 #include <cmath>
+#include <limits>
 
 
 // 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<long double>(std::numeric_limits<double>::epsilon() / 100),
+                static_cast<long double>(std::numeric_limits<long double>::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);
index b63fcd672ba74b472f7353358aa0008a639f1dbd..c3aa108f66b9e11ab3328aac3698bc94b130625c 100644 (file)
@@ -98,6 +98,15 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>base</h3>
 
 <ol>
+  <li> <p> Fixed: The constructor of the <code class="class">QGauss</code> 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.
+       <br>
+       (WB 2002/09/19)
+       </p>
+
   <li> <p> New: The <code class="class">Function</code> class now
        exports the value of its template argument through the static
        member variable <code class="member">dimension</code>.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.