]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new class QGauss works
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 29 Jan 2001 07:59:16 +0000 (07:59 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 29 Jan 2001 07:59:16 +0000 (07:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@3815 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/quadrature_lib.cc
tests/base/quadrature_test.cc

index dc3ef1282a4413aa70598d7fa9adabb8648bb003..6e5172ea6acdf07bb3af1fbdfa103894cccc20ba 100644 (file)
 // for all lower dimensions as well. That is why in this file the check
 // is for deal_II_dimension >= any_number and not for ==
 
+#ifndef M_PI
+# define M_PI          3.14159265358979323846  /* pi */
+#endif
+
+template <typename number>
+number abs (const number a)
+{
+  return ((a>0) ? a : -a);
+}
+
+
+template <>
+QGauss<1>::QGauss (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)
+    {
+      z = cos(M_PI * (i-.25)/(n+.5));
+
+                                      // Newton-iteration
+      do
+       {
+                                          // compute L_n (z)
+         p1 = 1.;
+         p2 = 0.;
+         for (unsigned int j=0;j<n;++j)
+           {
+             p3 = p2;
+             p2 = p1;
+             p1 = ((2.*j+1.)*z*p2-j*p3)/(j+1);
+           }
+         pp = n*(z*p1-p2)/(z*z-1);
+         z = z-p1/pp;
+       }
+      while (abs(p1/pp) > 1.e-19);
+
+      double x = .5*z;
+      quadrature_points[i-1] = Point<1>(.5-x);
+      quadrature_points[n-i] = Point<1>(.5+x);
+      
+      double w = 1./((1.-z*z)*pp*pp);
+      weights[i-1] = w;
+      weights[n-i] = w;
+    }
+}
+
+
+
 template <>
 QGauss2<1>::QGauss2 () :
                Quadrature<1> (2)
@@ -303,6 +356,14 @@ QWeddle<1>::QWeddle () :
 
 // construct the quadrature formulae in higher dimensions by
 // tensor product of lower dimensions
+
+template <int dim>
+QGauss<dim>::QGauss (unsigned int n)
+  :  Quadrature<dim> (QGauss<dim-1>(n), QGauss<1>(n))
+{};
+
+
+
 template <int dim>
 QGauss2<dim>::QGauss2 () :  Quadrature<dim> (QGauss2<dim-1>(), QGauss2<1>())  {};
 
@@ -350,6 +411,7 @@ QWeddle<dim>::QWeddle () :
 
 // explicite specialization
 // note that 1d formulae are specialized by implementation above
+template class QGauss<2>;
 template class QGauss2<2>;
 template class QGauss3<2>;
 template class QGauss4<2>;
@@ -362,6 +424,7 @@ template class QSimpson<2>;
 template class QMilne<2>;
 template class QWeddle<2>;
 
+template class QGauss<3>;
 template class QGauss2<3>;
 template class QGauss3<3>;
 template class QGauss4<3>;
index 526a270a37769b41413af47187e3196ac4f36eaa..3223f473396ff8f5af50218b9612d70352998a5e 100644 (file)
@@ -34,6 +34,10 @@ fill_vector (vector<Quadrature<dim> *>& quadratures)
   quadratures.push_back (new QSimpson<dim>());
   quadratures.push_back (new QMilne<dim>());
   quadratures.push_back (new QWeddle<dim>());
+  for (unsigned int i=1;i<9;++i)
+    {
+      quadratures.push_back (new QGauss<dim>(i));
+    }
 }
 
 template <int dim>
@@ -161,7 +165,7 @@ check_faces (vector<Quadrature<dim-1>*>& quadratures, bool sub)
       
          err = fabs(quadrature_int-exact_int);
        }
-      while (err<1e-14);
+      while (err<1e-16);
                                       // Uncomment here for testing
 //      deallog << " (Int " << quadrature_int << ',' << exact_int << ")";
       deallog << " is exact for polynomials of degree " << i-1 << endl;

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.