From: guido Date: Mon, 21 Jun 1999 15:39:42 +0000 (+0000) Subject: iterated simpson rule X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=02bc1626b9d7ab36b5bf3de42c30e885626c926b;p=dealii-svn.git iterated simpson rule git-svn-id: https://svn.dealii.org/trunk@1425 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 1fbe65054b..6982cd08aa 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -10,12 +10,13 @@ /** - * Second order Gauss quadrature formula. + * 2-Point-Gauss quadrature formula. * * Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. */ template -class QGauss2 : public Quadrature { +class QGauss2 : public Quadrature +{ public: QGauss2 (); }; @@ -23,12 +24,13 @@ class QGauss2 : public Quadrature { /** - * Third order Gauss quadrature formula. + * 3-Point-Gauss quadrature formula. * * Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. */ template -class QGauss3 : public Quadrature { +class QGauss3 : public Quadrature +{ public: QGauss3 (); }; @@ -36,12 +38,13 @@ class QGauss3 : public Quadrature { /** - * Fourth order Gauss quadrature formula. + * 4-Point-Gauss quadrature formula. * * Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. */ template -class QGauss4 : public Quadrature { +class QGauss4 : public Quadrature +{ public: QGauss4 (); }; @@ -50,12 +53,13 @@ class QGauss4 : public Quadrature { /** - * Fifth order Gauss quadrature formula. + * 5-Point-Gauss quadrature formula. * * Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. */ template -class QGauss5 : public Quadrature { +class QGauss5 : public Quadrature +{ public: QGauss5 (); }; @@ -63,7 +67,7 @@ class QGauss5 : public Quadrature { /** - * Sixth order Gauss quadrature formula. I have not found explicite + * 6-Point-Gauss quadrature formula. I have not found explicite * representations of the zeros of the Legendre functions of sixth * and higher degree. If anyone finds them, please replace the existing * numbers by these expressions. @@ -72,7 +76,8 @@ class QGauss5 : public Quadrature { * Element Methods */ template -class QGauss6 : public Quadrature { +class QGauss6 : public Quadrature +{ public: QGauss6 (); }; @@ -80,7 +85,7 @@ class QGauss6 : public Quadrature { /** - * Seventh order Gauss quadrature formula. I have not found explicite + * 7-Point-Gauss quadrature formula. I have not found explicite * representations of the zeros of the Legendre functions of sixth * and higher degree. If anyone finds them, please replace the existing * numbers by these expressions. @@ -89,7 +94,8 @@ class QGauss6 : public Quadrature { * Element Methods */ template -class QGauss7 : public Quadrature { +class QGauss7 : public Quadrature +{ public: QGauss7 (); }; @@ -97,7 +103,7 @@ class QGauss7 : public Quadrature { /** - * Eighth order Gauss quadrature formula. I have not found explicite + * 8-Point-Gauss quadrature formula. I have not found explicite * representations of the zeros of the Legendre functions of sixth * and higher degree. If anyone finds them, please replace the existing * numbers by these expressions. @@ -106,7 +112,8 @@ class QGauss7 : public Quadrature { * Element Methods */ template -class QGauss8 : public Quadrature { +class QGauss8 : public Quadrature +{ public: QGauss8 (); }; @@ -117,20 +124,23 @@ class QGauss8 : public Quadrature { /** * First order midpoint quadrature rule. + * For compatibility, this rule may be accessed as #QGauss1#, too. */ template -class QMidpoint : public Quadrature { +class QMidpoint : public Quadrature +{ public: QMidpoint (); }; - +#define QGauss1 QMidpoint /** * Simpson quadrature rule. */ template -class QSimpson : public Quadrature { +class QSimpson : public Quadrature +{ public: QSimpson (); }; @@ -141,22 +151,22 @@ class QSimpson : public Quadrature { * Trapezoidal quadrature rule. */ template -class QTrapez : public Quadrature { +class QTrapez : public Quadrature +{ public: QTrapez (); }; /** - * Iterated trapezoidal rule. The original reason for implementing - * this was gnuplot output for patches. Before this comment is - * removed, please do not rely on the correctness of the quadrature - * weights. + * Iterated trapezoidal rule. The aim of this class is to provide a + * low order formula, where the error constant can be tuned by + * increasing the number of quadrature points. This is useful in + * integrating non-differentiable functions on cells. * * For internal use, it may be worth to know that the points are * ordered in a fashion such that the last coordinate is the one which - * runs fastest, and when overflowing moving ahead the second last - * coordinate, and so on. + * runs fastest and then lexicographically from back to front. */ template class QIteratedTrapez : @@ -166,6 +176,19 @@ public: QIteratedTrapez(const unsigned intervals); }; +/** + * Iterated Simpson rule. + * Like #QIteratedTrapez#, this class provides a lower order formula, + * while the error constant can be tuned by choosing the number of sub-cells. + */ +template +class QIteratedSimpson : + public Quadrature +{ +public: + QIteratedSimpson(const unsigned intervals); +}; + /*---------------------------- quadrature_lib.h ---------------------------*/ /* end of #ifndef __quadrature_lib_H */ #endif diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 27df229525..c8ef626694 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -329,6 +329,24 @@ QIteratedTrapez<1>::QIteratedTrapez (const unsigned n) : weights[n] = .5/n; } +template<> +QIteratedSimpson<1>::QIteratedSimpson(const unsigned n) : + Quadrature<1> (2*n+1) +{ + weights.clear(); + + for (unsigned int i=0 ; i(1.*i/n); + quadrature_points[2*i+1] = Point<1>(1.*i/n+.5/n); + + weights[2*i] += 1./(6*n); + weights[2*i+1] += 4./(6*n); + weights[2*i+2] += 1./(6*n); + } + quadrature_points[2*n] = Point<1>(1.); +} + // construct the quadrature formulae in higher dimensions by // tensor product of lower dimensions template @@ -375,6 +393,10 @@ template QIteratedTrapez::QIteratedTrapez (const unsigned n) : Quadrature (QIteratedTrapez(n), QIteratedTrapez<1>(n)) {}; +template +QIteratedSimpson::QIteratedSimpson (const unsigned n) : + Quadrature (QIteratedTrapez(n), QIteratedTrapez<1>(n)) {}; + // explicite specialization @@ -389,6 +411,7 @@ template class QMidpoint<2>; template class QSimpson<2>; template class QTrapez<2>; template class QIteratedTrapez<2>; +template class QIteratedSimpson<2>; template class QGauss2<3>; template class QGauss3<3>; @@ -401,4 +424,4 @@ template class QMidpoint<3>; template class QSimpson<3>; template class QTrapez<3>; template class QIteratedTrapez<3>; - +template class QIteratedSimpson<3>;