From 2cfb65e814dd5a021cde25c6a88c34826b1f118a Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 14 Jun 2010 20:41:57 +0000 Subject: [PATCH] New test. git-svn-id: https://svn.dealii.org/trunk@21207 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/base/Makefile | 4 +- tests/base/gauss_lobatto.cc | 223 +++++++++++++++++++++++++++ tests/base/gauss_lobatto/cmp/generic | 208 +++++++++++++++++++++++++ 3 files changed, 434 insertions(+), 1 deletion(-) create mode 100644 tests/base/gauss_lobatto.cc create mode 100644 tests/base/gauss_lobatto/cmp/generic diff --git a/tests/base/Makefile b/tests/base/Makefile index 7955fa707b..0eed0491b5 100644 --- a/tests/base/Makefile +++ b/tests/base/Makefile @@ -65,7 +65,9 @@ tests_x = geometry_info_* \ logstream_end \ log_* \ reference \ - quadrature_* qprojector \ + quadrature_* \ + qprojector \ + gauss_lobatto \ path_search \ named_data \ slice_vector \ diff --git a/tests/base/gauss_lobatto.cc b/tests/base/gauss_lobatto.cc new file mode 100644 index 0000000000..91c3d24d0b --- /dev/null +++ b/tests/base/gauss_lobatto.cc @@ -0,0 +1,223 @@ +//---------------------------- gauss_lobatto.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- gauss_lobatto.cc --------------------------- + + +// check accuracy of various quadrature formulas by using them to +// integrate polynomials of increasing degree, and finding the degree +// until which they integrate exactly + + +#include "../tests.h" +#include +#include + +#include +#include +#include +#include + +template +void +fill_vector (std::vector *>& quadratures) +{ + quadratures.push_back (new QMidpoint()); + quadratures.push_back (new QTrapez()); + quadratures.push_back (new QSimpson()); + quadratures.push_back (new QMilne()); + quadratures.push_back (new QWeddle()); + for (unsigned int i=0;i<9;++i) + { + quadratures.push_back (new QGauss(i)); + } + QMilne<1> q1d; + quadratures.push_back (new Quadrature(q1d)); + for (unsigned int i=2;i<8;++i) + { + quadratures.push_back (new QGaussLobatto(i)); + } +} + +template +void +check_cells (std::vector*>& quadratures) +{ + Quadrature quadrature; + for (unsigned int n=0; n > &points=quadrature.get_points(); + const std::vector &weights=quadrature.get_weights(); + + deallog << "Quadrature no." << n; + + unsigned int i=0; + double quadrature_int=0; + double exact_int=0; + double err = 0; + + do + { + ++i; + + quadrature_int=0; + // Check the polynomial x^i*y^i + + for (unsigned int x=0; x(points[x](2)), i*1.0); + case 2: + f *= std::pow(static_cast(points[x](1)), i*1.0); + case 1: + f *= std::pow(static_cast(points[x](0)), i*1.0); + } + quadrature_int+=f*weights[x]; + } + + // the exact integral is 1/(i+1) + exact_int=1./std::pow(static_cast(i+1),dim); + err = std::fabs(quadrature_int-exact_int); + } + while (err<1e-14); + // Uncomment here for testing +// deallog << " (Int " << quadrature_int << ',' << exact_int << ")"; + deallog << " is exact for polynomials of degree " << i-1 << std::endl; + + if (dim==1) + { + // check the ordering of + // the quadrature points + bool in_order=true; + for (unsigned int x=1; x +void +check_faces (const std::vector*>& quadratures, const bool sub) +{ + if (sub) + deallog.push("subfaces"); + else + deallog.push("faces"); + + for (unsigned int n=0; n quadrature (sub == false? + QProjector::project_to_all_faces(*quadratures[n]) : + QProjector::project_to_all_subfaces(*quadratures[n])); + const std::vector > &points=quadrature.get_points(); + const std::vector &weights=quadrature.get_weights(); + + deallog << "Quadrature no." << n; + + unsigned int i=0; + long double quadrature_int=0; + double exact_int=0; + double err = 0; + + do + { + ++i; + + quadrature_int=0; + // Check the polynomial + // x^i*y^i*z^i + + for (unsigned int x=0; x