New test.

git-svn-id: https://svn.dealii.org/trunk@21207 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/base/Makefile b/tests/base/Makefile
index 7955fa707bd311e145861abd112de2e32bcb1cd6..0eed0491b58c6fc4a41b5124b0b533ab7f9f0ddb 100644 (file)
--- 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 (file)
index 0000000..91c3d24
--- /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 <iomanip>
+#include <fstream>
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <base/qprojector.h>
+#include <cmath>
+
+template <int dim>
+void
+fill_vector (std::vector<Quadrature<dim> *>& quadratures)
+{
+  quadratures.push_back (new QMidpoint<dim>());
+  quadratures.push_back (new QTrapez<dim>());
+  quadratures.push_back (new QSimpson<dim>());
+  quadratures.push_back (new QMilne<dim>());
+  quadratures.push_back (new QWeddle<dim>());
+  for (unsigned int i=0;i<9;++i)
+    {
+      quadratures.push_back (new QGauss<dim>(i));
+    }
+  QMilne<1> q1d;
+  quadratures.push_back (new Quadrature<dim>(q1d));
+  for (unsigned int i=2;i<8;++i)
+    {
+      quadratures.push_back (new QGaussLobatto<dim>(i));
+    }
+}
+
+template <int dim>
+void
+check_cells (std::vector<Quadrature<dim>*>& quadratures)
+{
+  Quadrature<dim> quadrature;
+  for (unsigned int n=0; n<quadratures.size(); ++n)
+    {
+      quadrature = *quadratures[n];
+      const std::vector<Point<dim> > &points=quadrature.get_points();
+      const std::vector<double> &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<quadrature.n_quadrature_points; ++x)
+           {
+             double f=1.;
+             switch (dim)
+               {
+               case 3:
+                 f *= std::pow(static_cast<double>(points[x](2)), i*1.0);
+               case 2:
+                 f *= std::pow(static_cast<double>(points[x](1)), i*1.0);
+               case 1:
+                 f *= std::pow(static_cast<double>(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<double>(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<quadrature.n_quadrature_points; ++x)
+           {
+             if (points[x](0)<=points[x-1](0))
+               in_order=false;
+           }
+         if (!in_order)
+           for (unsigned int x=0; x<quadrature.n_quadrature_points; ++x)
+             deallog << points[x] << std::endl;
+       }
+    }
+}
+
+
+template <int dim>
+void
+check_faces (const std::vector<Quadrature<dim-1>*>& quadratures, const bool sub)
+{
+  if (sub)
+    deallog.push("subfaces");
+  else
+    deallog.push("faces");
+
+  for (unsigned int n=0; n<quadratures.size(); ++n)
+    {
+      Quadrature<dim> quadrature (sub == false?
+                                 QProjector<dim>::project_to_all_faces(*quadratures[n]) :
+                                 QProjector<dim>::project_to_all_subfaces(*quadratures[n]));
+      const std::vector<Point<dim> > &points=quadrature.get_points();
+      const std::vector<double> &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<quadrature.n_quadrature_points; ++x)
+           {
+             long double f=1.;
+             switch (dim)
+               {
+               case 3:
+                 f *= std::pow((long double) points[x](2), i*1.0L);
+               case 2:
+                 f *= std::pow((long double) points[x](1), i*1.0L);
+               case 1:
+                 f *= std::pow((long double) points[x](0), i*1.0L);
+               }
+             quadrature_int+=f*weights[x];
+           }
+
+                                          // the exact integral is
+                                          // 1/(i+1)^(dim-1)
+         switch (dim)
+           {
+           case 2:
+             exact_int = 2 * (sub ? 2:1) / (double) (i+1);
+             break;
+           case 3:
+             exact_int = 3 * (sub ? (4+2+2):1)*8 / (double) (i+1)/(i+1);
+             break;
+           }
+
+         err = std::fabs(quadrature_int-exact_int);
+       }
+                                      // for comparison: use factor 8 in case
+                                      // of dim==3, as we integrate 8 times
+                                      // over the whole surface (all
+                                      // combinations of face_orientation,
+                                      // face_flip and face_rotation)
+      while (err < (dim==3 ? 8 : 1) * 2e-14);
+                                      // Uncomment here for testing
+      //      deallog << " (Int " << quadrature_int << '-' << exact_int << '=' << err << ")";
+      deallog << " is exact for polynomials of degree " << i-1 << std::endl;
+    }
+  deallog.pop();
+}
+
+int main()
+{
+  std::ofstream logfile("gauss_lobatto/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  for (unsigned int n=2; n<20; ++n)
+    {
+      deallog << "QGaussLobatto(" << n << ")" << std::endl;
+
+      QGaussLobatto<1> q(n);
+      for (unsigned int i=0; i<q.n_quadrature_points; ++i)
+       deallog << q.point(i) << ' ' << q.weight(i) << std::endl;
+
+                                      // the points must be
+                                      // symmetrically located around
+                                      // 0.5
+      double p = 0;
+      for (unsigned int i=0; i<q.n_quadrature_points; ++i)
+       p += (q.point(i)[0] - 0.5);
+      Assert (std::fabs(p) < 1e-12, ExcInternalError());
+
+                                      // the sum of weights must be
+                                      // one
+      double w = 0;
+      for (unsigned int i=0; i<q.n_quadrature_points; ++i)
+       w += q.weight(i);
+      Assert (std::fabs(w-1) < 1e-12, ExcInternalError());
+    }
+}
+
+

diff --git a/tests/base/gauss_lobatto/cmp/generic b/tests/base/gauss_lobatto/cmp/generic
new file mode 100644 (file)
index 0000000..9cec587
--- /dev/null
+++ b/tests/base/gauss_lobatto/cmp/generic
@@ -0,0 +1,208 @@
+
+DEAL::QGaussLobatto(2)
+DEAL::0.00000 0.500000
+DEAL::1.00000 0.500000
+DEAL::QGaussLobatto(3)
+DEAL::0.00000 0.166667
+DEAL::0.500000 0.666667
+DEAL::1.00000 0.166667
+DEAL::QGaussLobatto(4)
+DEAL::0.00000 0.0833333
+DEAL::0.276393 0.416667
+DEAL::0.723607 0.416667
+DEAL::1.00000 0.0833333
+DEAL::QGaussLobatto(5)
+DEAL::0.00000 0.0500000
+DEAL::0.172673 0.272222
+DEAL::0.500000 0.355556
+DEAL::0.827327 0.272222
+DEAL::1.00000 0.0500000
+DEAL::QGaussLobatto(6)
+DEAL::0.00000 0.0333333
+DEAL::0.117472 0.189237
+DEAL::0.357384 0.277429
+DEAL::0.642616 0.277429
+DEAL::0.882528 0.189237
+DEAL::1.00000 0.0333333
+DEAL::QGaussLobatto(7)
+DEAL::0.00000 0.0238095
+DEAL::0.0848881 0.138413
+DEAL::0.265576 0.215873
+DEAL::0.500000 0.243810
+DEAL::0.734424 0.215873
+DEAL::0.915112 0.138413
+DEAL::1.00000 0.0238095
+DEAL::QGaussLobatto(8)
+DEAL::0.00000 0.0178571
+DEAL::0.0641299 0.105352
+DEAL::0.204150 0.170561
+DEAL::0.395350 0.206229
+DEAL::0.604650 0.206229
+DEAL::0.795850 0.170561
+DEAL::0.935870 0.105352
+DEAL::1.00000 0.0178571
+DEAL::QGaussLobatto(9)
+DEAL::0.00000 0.0138889
+DEAL::0.0501210 0.0827477
+DEAL::0.161407 0.137269
+DEAL::0.318441 0.173214
+DEAL::0.500000 0.185760
+DEAL::0.681559 0.173214
+DEAL::0.838593 0.137269
+DEAL::0.949879 0.0827477
+DEAL::1.00000 0.0138889
+DEAL::QGaussLobatto(10)
+DEAL::0.00000 0.0111111
+DEAL::0.0402330 0.0666530
+DEAL::0.130613 0.112445
+DEAL::0.261038 0.146021
+DEAL::0.417361 0.163770
+DEAL::0.582639 0.163770
+DEAL::0.738962 0.146021
+DEAL::0.869387 0.112445
+DEAL::0.959767 0.0666530
+DEAL::1.00000 0.0111111
+DEAL::QGaussLobatto(11)
+DEAL::0.00000 0.00909091
+DEAL::0.0329993 0.0548061
+DEAL::0.107758 0.0935849
+DEAL::0.217382 0.124024
+DEAL::0.352121 0.143440
+DEAL::0.500000 0.150109
+DEAL::0.647879 0.143440
+DEAL::0.782618 0.124024
+DEAL::0.892242 0.0935849
+DEAL::0.967001 0.0548061
+DEAL::1.00000 0.00909091
+DEAL::QGaussLobatto(12)
+DEAL::0.00000 0.00757576
+DEAL::0.0275504 0.0458423
+DEAL::0.0903603 0.0789874
+DEAL::0.183562 0.106254
+DEAL::0.300235 0.125638
+DEAL::0.431724 0.135703
+DEAL::0.568276 0.135703
+DEAL::0.699765 0.125638
+DEAL::0.816438 0.106254
+DEAL::0.909640 0.0789874
+DEAL::0.972450 0.0458423
+DEAL::1.00000 0.00757576
+DEAL::QGaussLobatto(13)
+DEAL::0.00000 0.00641026
+DEAL::0.0233451 0.0389008
+DEAL::0.0768262 0.0674910
+DEAL::0.156906 0.0918234
+DEAL::0.258545 0.110384
+DEAL::0.375357 0.122008
+DEAL::0.500000 0.125965
+DEAL::0.624643 0.122008
+DEAL::0.741455 0.110384
+DEAL::0.843094 0.0918234
+DEAL::0.923174 0.0674910
+DEAL::0.976655 0.0389008
+DEAL::1.00000 0.00641026
+DEAL::QGaussLobatto(14)
+DEAL::0.00000 0.00549451
+DEAL::0.0200325 0.0334186
+DEAL::0.0660995 0.0582933
+DEAL::0.135566 0.0800109
+DEAL::0.224680 0.0974131
+DEAL::0.328638 0.109563
+DEAL::0.441834 0.115806
+DEAL::0.558166 0.115806
+DEAL::0.671362 0.109563
+DEAL::0.775320 0.0974131
+DEAL::0.864434 0.0800109
+DEAL::0.933901 0.0582933
+DEAL::0.979968 0.0334186
+DEAL::1.00000 0.00549451
+DEAL::QGaussLobatto(15)
+DEAL::0.00000 0.00476190
+DEAL::0.0173770 0.0290149
+DEAL::0.0574590 0.0508300
+DEAL::0.118240 0.0702558
+DEAL::0.196873 0.0863948
+DEAL::0.289681 0.0984936
+DEAL::0.392323 0.105987
+DEAL::0.500000 0.108524
+DEAL::0.607677 0.105987
+DEAL::0.710319 0.0984936
+DEAL::0.803127 0.0863948
+DEAL::0.881760 0.0702558
+DEAL::0.942541 0.0508300
+DEAL::0.982623 0.0290149
+DEAL::1.00000 0.00476190
+DEAL::QGaussLobatto(16)
+DEAL::0.00000 0.00416667
+DEAL::0.0152160 0.0254252
+DEAL::0.0503997 0.0446968
+DEAL::0.103996 0.0621277
+DEAL::0.173806 0.0770135
+DEAL::0.256970 0.0887460
+DEAL::0.350085 0.0968450
+DEAL::0.449337 0.100979
+DEAL::0.550663 0.100979
+DEAL::0.649915 0.0968450
+DEAL::0.743030 0.0887460
+DEAL::0.826194 0.0770135
+DEAL::0.896004 0.0621277
+DEAL::0.949600 0.0446968
+DEAL::0.984784 0.0254252
+DEAL::1.00000 0.00416667
+DEAL::QGaussLobatto(17)
+DEAL::0.00000 0.00367647
+DEAL::0.0134339 0.0224610
+DEAL::0.0445600 0.0395991
+DEAL::0.0921519 0.0552965
+DEAL::0.154486 0.0689939
+DEAL::0.229307 0.0801973
+DEAL::0.313913 0.0885021
+DEAL::0.405244 0.0936082
+DEAL::0.500000 0.0953309
+DEAL::0.594756 0.0936082
+DEAL::0.686087 0.0885021
+DEAL::0.770693 0.0801973
+DEAL::0.845514 0.0689939
+DEAL::0.907848 0.0552965
+DEAL::0.955440 0.0395991
+DEAL::0.986566 0.0224610
+DEAL::1.00000 0.00367647
+DEAL::QGaussLobatto(18)
+DEAL::0.00000 0.00326797
+DEAL::0.0119472 0.0199853
+DEAL::0.0396754 0.0353186
+DEAL::0.0822032 0.0495081
+DEAL::0.138160 0.0621053
+DEAL::0.205748 0.0727060
+DEAL::0.282792 0.0809698
+DEAL::0.366819 0.0866311
+DEAL::0.455125 0.0895079
+DEAL::0.544875 0.0895079
+DEAL::0.633181 0.0866311
+DEAL::0.717208 0.0809698
+DEAL::0.794252 0.0727060
+DEAL::0.861840 0.0621053
+DEAL::0.917797 0.0495081
+DEAL::0.960325 0.0353186
+DEAL::0.988053 0.0199853
+DEAL::1.00000 0.00326797
+DEAL::QGaussLobatto(19)
+DEAL::0.00000 0.00292398
+DEAL::0.0106941 0.0178967
+DEAL::0.0355492 0.0316909
+DEAL::0.0737697 0.0445659
+DEAL::0.124253 0.0561577
+DEAL::0.185546 0.0661336
+DEAL::0.255885 0.0742070
+DEAL::0.333248 0.0801455
+DEAL::0.415407 0.0837783
+DEAL::0.500000 0.0850010
+DEAL::0.584593 0.0837783
+DEAL::0.666752 0.0801455
+DEAL::0.744115 0.0742070
+DEAL::0.814454 0.0661336
+DEAL::0.875747 0.0561577
+DEAL::0.926230 0.0445659
+DEAL::0.964451 0.0316909
+DEAL::0.989306 0.0178967
+DEAL::1.00000 0.00292398