ofstream logfile("quadrature_test.output");
deallog.attach(logfile);
deallog.depth_console(0);
-
- vector<Quadrature<1> *> quadratures(14);
- quadratures[2]=new QGauss2<1>();
- quadratures[3]=new QGauss3<1>();
- quadratures[4]=new QGauss4<1>();
- quadratures[5]=new QGauss5<1>();
- quadratures[6]=new QGauss6<1>();
- quadratures[7]=new QGauss7<1>();
- quadratures[8]=new QGauss8<1>();
- quadratures[9] = new QMidpoint<1>();
- quadratures[10] = new QTrapez<1>();
- quadratures[11] = new QSimpson<1>();
- quadratures[12] = new QMilne<1>();
- quadratures[13] = new QWeddle<1>();
+ vector<Quadrature<2> *> quadratures(14);
+ quadratures[2]=new QGauss2<2>();
+ quadratures[3]=new QGauss3<2>();
+ quadratures[4]=new QGauss4<2>();
+ quadratures[5]=new QGauss5<2>();
+ quadratures[6]=new QGauss6<2>();
+ quadratures[7]=new QGauss7<2>();
+ quadratures[8]=new QGauss8<2>();
+ quadratures[9] = new QMidpoint<2>();
+ quadratures[10] = new QTrapez<2>();
+ quadratures[11] = new QSimpson<2>();
+ quadratures[12] = new QMilne<2>();
+ quadratures[13] = new QWeddle<2>();
for (unsigned int n=2; n<14; ++n)
{
- Quadrature<1> *quadrature=quadratures[n];
- const vector<Point<1> > &points=quadrature->get_points();
+ Quadrature<2> *quadrature=quadratures[n];
+ const vector<Point<2> > &points=quadrature->get_points();
const vector<double> &weights=quadrature->get_weights();
deallog << "Quadrature no." << n;
++i;
quadrature_int=0;
+
+ // Check the polynomial x^i*y^i
+
for (unsigned int x=0; x<quadrature->n_quadrature_points; ++x)
- quadrature_int+=pow(points[x](0), i)*weights[x];
+ quadrature_int+=pow(points[x](0), i)*pow(points[x](1), i)*weights[x];
// the exact integral is 1/(i+1)
- exact_int=1./(i+1);
+ exact_int=1./(i+1)/(i+1);
err = fabs(quadrature_int-exact_int);
}
while (err<1e-15);