/**
- * Second order Gauss quadrature formula.
+ * 2-Point-Gauss quadrature formula.
*
* Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
*/
template <int dim>
-class QGauss2 : public Quadrature<dim> {
+class QGauss2 : public Quadrature<dim>
+{
public:
QGauss2 ();
};
/**
- * Third order Gauss quadrature formula.
+ * 3-Point-Gauss quadrature formula.
*
* Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
*/
template <int dim>
-class QGauss3 : public Quadrature<dim> {
+class QGauss3 : public Quadrature<dim>
+{
public:
QGauss3 ();
};
/**
- * Fourth order Gauss quadrature formula.
+ * 4-Point-Gauss quadrature formula.
*
* Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
*/
template <int dim>
-class QGauss4 : public Quadrature<dim> {
+class QGauss4 : public Quadrature<dim>
+{
public:
QGauss4 ();
};
/**
- * Fifth order Gauss quadrature formula.
+ * 5-Point-Gauss quadrature formula.
*
* Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
*/
template <int dim>
-class QGauss5 : public Quadrature<dim> {
+class QGauss5 : public Quadrature<dim>
+{
public:
QGauss5 ();
};
/**
- * 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.
* Element Methods
*/
template <int dim>
-class QGauss6 : public Quadrature<dim> {
+class QGauss6 : public Quadrature<dim>
+{
public:
QGauss6 ();
};
/**
- * 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.
* Element Methods
*/
template <int dim>
-class QGauss7 : public Quadrature<dim> {
+class QGauss7 : public Quadrature<dim>
+{
public:
QGauss7 ();
};
/**
- * 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.
* Element Methods
*/
template <int dim>
-class QGauss8 : public Quadrature<dim> {
+class QGauss8 : public Quadrature<dim>
+{
public:
QGauss8 ();
};
/**
* First order midpoint quadrature rule.
+ * For compatibility, this rule may be accessed as #QGauss1#, too.
*/
template <int dim>
-class QMidpoint : public Quadrature<dim> {
+class QMidpoint : public Quadrature<dim>
+{
public:
QMidpoint ();
};
-
+#define QGauss1 QMidpoint
/**
* Simpson quadrature rule.
*/
template <int dim>
-class QSimpson : public Quadrature<dim> {
+class QSimpson : public Quadrature<dim>
+{
public:
QSimpson ();
};
* Trapezoidal quadrature rule.
*/
template <int dim>
-class QTrapez : public Quadrature<dim> {
+class QTrapez : public Quadrature<dim>
+{
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 <int dim>
class QIteratedTrapez :
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 <int dim>
+class QIteratedSimpson :
+ public Quadrature<dim>
+{
+public:
+ QIteratedSimpson(const unsigned intervals);
+};
+
/*---------------------------- quadrature_lib.h ---------------------------*/
/* end of #ifndef __quadrature_lib_H */
#endif
weights[n] = .5/n;
}
+template<>
+QIteratedSimpson<1>::QIteratedSimpson(const unsigned n) :
+ Quadrature<1> (2*n+1)
+{
+ weights.clear();
+
+ for (unsigned int i=0 ; i<n; ++i)
+ {
+ quadrature_points[2*i] = Point<1>(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 <int dim>
QIteratedTrapez<dim>::QIteratedTrapez (const unsigned n) :
Quadrature<dim> (QIteratedTrapez<dim-1>(n), QIteratedTrapez<1>(n)) {};
+template <int dim>
+QIteratedSimpson<dim>::QIteratedSimpson (const unsigned n) :
+ Quadrature<dim> (QIteratedTrapez<dim-1>(n), QIteratedTrapez<1>(n)) {};
+
// explicite specialization
template class QSimpson<2>;
template class QTrapez<2>;
template class QIteratedTrapez<2>;
+template class QIteratedSimpson<2>;
template class QGauss2<3>;
template class QGauss3<3>;
template class QSimpson<3>;
template class QTrapez<3>;
template class QIteratedTrapez<3>;
-
+template class QIteratedSimpson<3>;