From: Guido Kanschat Date: Mon, 29 Jan 2001 07:59:16 +0000 (+0000) Subject: new class QGauss works X-Git-Tag: v8.0.0~19770 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c2d19d121770574e717d5c4190dd82b75402113f;p=dealii.git new class QGauss works git-svn-id: https://svn.dealii.org/trunk@3815 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index dc3ef1282a..6e5172ea6a 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -20,6 +20,59 @@ // 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 +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 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 +QGauss::QGauss (unsigned int n) + : Quadrature (QGauss(n), QGauss<1>(n)) +{}; + + + template QGauss2::QGauss2 () : Quadrature (QGauss2(), QGauss2<1>()) {}; @@ -350,6 +411,7 @@ QWeddle::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>; diff --git a/tests/base/quadrature_test.cc b/tests/base/quadrature_test.cc index 526a270a37..3223f47339 100644 --- a/tests/base/quadrature_test.cc +++ b/tests/base/quadrature_test.cc @@ -34,6 +34,10 @@ fill_vector (vector *>& quadratures) quadratures.push_back (new QSimpson()); quadratures.push_back (new QMilne()); quadratures.push_back (new QWeddle()); + for (unsigned int i=1;i<9;++i) + { + quadratures.push_back (new QGauss(i)); + } } template @@ -161,7 +165,7 @@ check_faces (vector*>& 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;