new high order closed Newton-Cotes formulae

git-svn-id: https://svn.dealii.org/trunk@3142 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 4446712a0b528b8bae00ad172d6a1962db4526ac..a2868f65e30c82352c97acb7493823cdea4a9c4e 100644 (file)
--- a/deal.II/base/include/base/quadrature_lib.h
+++ b/deal.II/base/include/base/quadrature_lib.h
@@ -169,5 +169,28 @@ class QTrapez : public Quadrature<dim>
     QTrapez ();
 };
 
+/**
+ * Milne-rule. Closed Newton-Cotes formula, exact for polynomials of degree 5.
+ * See Stoer: Einführung in die Numerische Mathematik I, p. 102
+ */
+template <int dim>
+class QMilne : public Quadrature<dim>
+{
+  public:
+    QMilne ();
+};
+
+
+/**
+ * Weddle-rule. Closed Newton-Cotes formula, exact for polynomials of degree 7.
+ * See Stoer: Einführung in die Numerische Mathematik I, p. 102
+ */
+template <int dim>
+class QWeddle : public Quadrature<dim>
+{
+  public:
+    QWeddle ();
+};
+
 
 #endif

diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc
index b4bc6549f3af583a1745235ba1134e77aef53997..2a0bcdce5541908d9770d790a44be2e0a11073f9 100644 (file)
--- a/deal.II/base/source/quadrature_lib.cc
+++ b/deal.II/base/source/quadrature_lib.cc
@@ -286,6 +286,21 @@ QMidpoint<1>::QMidpoint () :
 };
 
 
+template <>
+QTrapez<1>::QTrapez () :
+               Quadrature<1> (2)
+{
+  static const double xpts[] = { 0.0, 1.0 };
+  static const double wts[]  = { 0.5, 0.5 };
+
+  for (unsigned int i=0; i<n_quadrature_points; ++i) 
+    {
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
+    };
+};
+
+
 template <>
 QSimpson<1>::QSimpson () :
                Quadrature<1> (3)
@@ -302,11 +317,28 @@ QSimpson<1>::QSimpson () :
 
 
 template <>
-QTrapez<1>::QTrapez () :
-               Quadrature<1> (2)
+QMilne<1>::QMilne () :
+               Quadrature<1> (5)
 {
-  static const double xpts[] = { 0.0, 1.0 };
-  static const double wts[]  = { 0.5, 0.5 };
+  static const double xpts[] = { 0.0, .25, .5, .75, 1.0 };
+  static const double wts[]  = { 7./90., 32./90., 12./90., 32./90., 7./90. };
+
+  for (unsigned int i=0; i<n_quadrature_points; ++i) 
+    {
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
+    };
+};
+
+
+template <>
+QWeddle<1>::QWeddle () :
+               Quadrature<1> (7)
+{
+  static const double xpts[] = { 0.0, 1./6., 1./3., .5, 2./3., 5./6., 1.0 };
+  static const double wts[]  = { 41./840., 216./840., 27./840., 272./840.,
+                                27./840., 216./840., 41./840.
+  };
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
@@ -316,6 +348,7 @@ QTrapez<1>::QTrapez () :
 };
 
 
+
 // construct the quadrature formulae in higher dimensions by
 // tensor product of lower dimensions
 template <int dim>
@@ -369,8 +402,10 @@ template class QGauss6<2>;
 template class QGauss7<2>;
 template class QGauss8<2>;
 template class QMidpoint<2>;
-template class QSimpson<2>;
 template class QTrapez<2>;
+template class QSimpson<2>;
+template class QMilne<2>;
+template class QWeddle<2>;
 
 template class QGauss2<3>;
 template class QGauss3<3>;
@@ -380,5 +415,7 @@ template class QGauss6<3>;
 template class QGauss7<3>;
 template class QGauss8<3>;
 template class QMidpoint<3>;
-template class QSimpson<3>;
 template class QTrapez<3>;
+template class QSimpson<3>;
+template class QMilne<3>;
+template class QWeddle<3>;