]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Polynomial evaluation at high orders: old implementation based on form a_0 + a_1...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Jul 2011 12:58:28 +0000 (12:58 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Jul 2011 12:58:28 +0000 (12:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@23959 0785d39b-7218-0410-832d-ea1e28bc413d

tests/base/polynomial_lagrange_ops.cc [new file with mode: 0644]
tests/base/polynomial_lagrange_ops/cmp/generic [new file with mode: 0644]
tests/base/polynomial_lagrange_order.cc [new file with mode: 0644]
tests/base/polynomial_lagrange_order/cmp/generic [new file with mode: 0644]
tests/base/polynomial_lagrange_product.cc [new file with mode: 0644]
tests/base/polynomial_lagrange_product/cmp/generic [new file with mode: 0644]

diff --git a/tests/base/polynomial_lagrange_ops.cc b/tests/base/polynomial_lagrange_ops.cc
new file mode 100644 (file)
index 0000000..0494cdd
--- /dev/null
@@ -0,0 +1,229 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2011 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.
+//
+//-----------------------------------------------------------------------------
+
+// Tests operations similar to polynomial_lagrange_order when the polynomial
+// is modified
+
+#include "../tests.h"
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+
+
+using namespace Polynomials;
+
+
+void check_scale (const std::vector<Polynomial<double> > &p)
+{
+  deallog << "Scale operation";
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      Polynomial<double> q = p[i];
+      double x = (double)rand()/RAND_MAX;
+      double factor = 5.*(double)rand()/RAND_MAX;
+      q.scale (factor);
+      double value1 = p[i].value (factor * x);
+      double value2 = q   .value (x);
+      if (std::fabs(value1-value2) > std::max(1e-13,1e-13*std::fabs(value1)))
+       deallog << "Error scale at x=" << x << ": p(t)=" << value1
+               << ", q(x)=" << value2 << std::endl;
+      deallog << ".";
+    }
+  deallog << std::endl;
+}
+
+
+
+void check_shift (const std::vector<Polynomial<double> > &p)
+{
+  deallog << "Shift operation";
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      Polynomial<double> q = p[i];
+      double x = (double)rand()/RAND_MAX;
+      double a = 10.*(-1.+2.*(double)rand()/RAND_MAX);
+      q.shift (a);
+      double value1 = p[i].value (x+a);
+      double value2 = q   .value (x);
+      if (std::fabs(value1-value2) > std::max(1e-13,1e-13*std::fabs(value1)))
+       deallog << "Error shift at x=" << x << ": p(t)=" << value1
+               << ", q(x)=" << value2 << std::endl;
+      deallog << ".";
+    }
+  deallog << std::endl;
+}
+
+
+
+void check_mult_scalar (const std::vector<Polynomial<double> > &p)
+{
+  deallog << "Multiply by scalar";
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      Polynomial<double> q = p[i];
+      double x = (double)rand()/RAND_MAX;
+      double a = (double)rand()/RAND_MAX;
+      q *= a;
+      double value1 = p[i].value (x) * a;
+      double value2 = q   .value (x);
+      if (std::fabs(value1-value2) > std::max(1e-13,1e-13*std::fabs(value1)))
+       deallog << "Error multiply at x=" << x << ": a*p(x)=" << value1
+               << ", q(x)=" << value2 << std::endl;
+      deallog << ".";
+    }
+  deallog << std::endl;
+}
+
+
+
+void check_mult (const std::vector<Polynomial<double> > &p)
+{
+  deallog << "Multiply by polynomial";
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      for (unsigned int j=0; j<p.size(); ++j)
+       {
+         Polynomial<double> q = p[i];
+         q *= p[j];
+         double x = (double)rand()/RAND_MAX;
+         double value1 = p[i].value (x) * p[j].value(x);
+         double value2 = q   .value (x);
+         if (std::fabs(value1-value2) > std::max(1e-13,1e-13*std::fabs(value1)))
+           deallog << "Error multiply at x=" << x << ": p_1(x)*p_2(x)=" << value1
+                   << ", q(x)=" << value2 << std::endl;
+       }
+      deallog << ".";
+    }
+  deallog << std::endl;
+}
+
+
+
+void check_expand (const std::vector<Polynomial<double> > &p)
+{
+  if (p.size() > 10)
+    return;
+                               // this checks whether the Lagrange product
+                               // form and the usual form with factors for
+                               // different powers does the same
+                               // thing. Realize this by adding the
+                               // polynomial '0' to the current
+                               // polynomial. That destroys the product form
+                               // and triggers the usual form. do not do this
+                               // for higher order because then the standard
+                               // form is unstable
+  deallog << "Expansion operation";
+  Monomial<double> zero(0, 0.);
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      Polynomial<double> q = p[i];
+      double x = (double)rand()/RAND_MAX;
+      q += zero;
+      double value1 = p[i].value (x);
+      double value2 = q   .value (x);
+      if (std::fabs(value1-value2) > std::max(1e-10,1e-10*std::fabs(value1)))
+       deallog << "Error expansion at x=" << x << ": p(x)=" << value1
+               << ", q(x)=" << value2 << std::endl;
+      deallog << ".";
+    }
+  deallog << std::endl;
+}
+
+
+
+void check_mult_expand (const std::vector<Polynomial<double> > &p)
+{
+  if (p.size() > 6)
+    return;
+  deallog << "Multiply by polynomial expanding";
+  Monomial<double> zero(0, 0.);
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      for (unsigned int j=0; j<p.size(); ++j)
+       {
+         Polynomial<double> q = p[i];
+         q += zero;
+         q *= p[j];
+         double x = (double)rand()/RAND_MAX;
+         double value1 = p[i].value (x) * p[j].value(x);
+         double value2 = q   .value (x);
+         if (std::fabs(value1-value2) > std::max(1e-10,1e-10*std::fabs(value1)))
+           deallog << "Error multiply at x=" << x
+                   << ": p_"<<i<< "(x)*p_"<<j<<"(x)=" << value1
+                   << ", q(x)=" << value2 << std::endl;
+       }
+      deallog << ".";
+    }
+  deallog << std::endl;
+}
+
+
+
+void
+check_lge (unsigned int n)
+{
+  deallog << "Points: " << n+1 << std::endl;
+  std::vector<Polynomial<double> > p = LagrangeEquidistant::generate_complete_basis(n);
+  check_scale (p);
+  check_shift (p);
+  check_mult_scalar(p);
+  check_mult  (p);
+  check_expand(p);
+  check_mult_expand(p);
+  deallog << std::endl;
+}
+
+
+
+void
+check_poly (const Quadrature<1>& q)
+{
+  deallog << "Points: " << q.size() << std::endl;
+  std::vector<Polynomial<double> > p = generate_complete_Lagrange_basis(q.get_points());
+  check_scale (p);
+  check_shift (p);
+  check_mult_scalar(p);
+  check_mult  (p);
+  check_expand(p);
+  check_mult_expand(p);
+  deallog << std::endl;
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("polynomial_lagrange_ops/output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog.push("LagrangeEquidistant");
+  for (unsigned i=1; i<10; i+=2)
+    check_lge (i);
+  deallog.pop();
+  deallog << std::endl;
+
+                               // Lagrange elements on GL points have good
+                               // conditioning, so test to some very high
+                               // orders
+  deallog.push("GaussLobatto");
+  for (unsigned i=1; i<40; i+=3)
+    check_poly (QGaussLobatto<1>(i+1));
+  deallog.pop();
+}
diff --git a/tests/base/polynomial_lagrange_ops/cmp/generic b/tests/base/polynomial_lagrange_ops/cmp/generic
new file mode 100644 (file)
index 0000000..d3b0a5b
--- /dev/null
@@ -0,0 +1,123 @@
+
+DEAL:LagrangeEquidistant::Points: 2
+DEAL:LagrangeEquidistant::Scale operation..
+DEAL:LagrangeEquidistant::Shift operation..
+DEAL:LagrangeEquidistant::Multiply by scalar..
+DEAL:LagrangeEquidistant::Multiply by polynomial..
+DEAL:LagrangeEquidistant::Expansion operation..
+DEAL:LagrangeEquidistant::Multiply by polynomial expanding..
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 4
+DEAL:LagrangeEquidistant::Scale operation....
+DEAL:LagrangeEquidistant::Shift operation....
+DEAL:LagrangeEquidistant::Multiply by scalar....
+DEAL:LagrangeEquidistant::Multiply by polynomial....
+DEAL:LagrangeEquidistant::Expansion operation....
+DEAL:LagrangeEquidistant::Multiply by polynomial expanding....
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 6
+DEAL:LagrangeEquidistant::Scale operation......
+DEAL:LagrangeEquidistant::Shift operation......
+DEAL:LagrangeEquidistant::Multiply by scalar......
+DEAL:LagrangeEquidistant::Multiply by polynomial......
+DEAL:LagrangeEquidistant::Expansion operation......
+DEAL:LagrangeEquidistant::Multiply by polynomial expanding......
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 8
+DEAL:LagrangeEquidistant::Scale operation........
+DEAL:LagrangeEquidistant::Shift operation........
+DEAL:LagrangeEquidistant::Multiply by scalar........
+DEAL:LagrangeEquidistant::Multiply by polynomial........
+DEAL:LagrangeEquidistant::Expansion operation........
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 10
+DEAL:LagrangeEquidistant::Scale operation..........
+DEAL:LagrangeEquidistant::Shift operation..........
+DEAL:LagrangeEquidistant::Multiply by scalar..........
+DEAL:LagrangeEquidistant::Multiply by polynomial..........
+DEAL:LagrangeEquidistant::Expansion operation..........
+DEAL:LagrangeEquidistant::
+DEAL::
+DEAL:GaussLobatto::Points: 2
+DEAL:GaussLobatto::Scale operation..
+DEAL:GaussLobatto::Shift operation..
+DEAL:GaussLobatto::Multiply by scalar..
+DEAL:GaussLobatto::Multiply by polynomial..
+DEAL:GaussLobatto::Expansion operation..
+DEAL:GaussLobatto::Multiply by polynomial expanding..
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 5
+DEAL:GaussLobatto::Scale operation.....
+DEAL:GaussLobatto::Shift operation.....
+DEAL:GaussLobatto::Multiply by scalar.....
+DEAL:GaussLobatto::Multiply by polynomial.....
+DEAL:GaussLobatto::Expansion operation.....
+DEAL:GaussLobatto::Multiply by polynomial expanding.....
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 8
+DEAL:GaussLobatto::Scale operation........
+DEAL:GaussLobatto::Shift operation........
+DEAL:GaussLobatto::Multiply by scalar........
+DEAL:GaussLobatto::Multiply by polynomial........
+DEAL:GaussLobatto::Expansion operation........
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 11
+DEAL:GaussLobatto::Scale operation...........
+DEAL:GaussLobatto::Shift operation...........
+DEAL:GaussLobatto::Multiply by scalar...........
+DEAL:GaussLobatto::Multiply by polynomial...........
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 14
+DEAL:GaussLobatto::Scale operation..............
+DEAL:GaussLobatto::Shift operation..............
+DEAL:GaussLobatto::Multiply by scalar..............
+DEAL:GaussLobatto::Multiply by polynomial..............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 17
+DEAL:GaussLobatto::Scale operation.................
+DEAL:GaussLobatto::Shift operation.................
+DEAL:GaussLobatto::Multiply by scalar.................
+DEAL:GaussLobatto::Multiply by polynomial.................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 20
+DEAL:GaussLobatto::Scale operation....................
+DEAL:GaussLobatto::Shift operation....................
+DEAL:GaussLobatto::Multiply by scalar....................
+DEAL:GaussLobatto::Multiply by polynomial....................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 23
+DEAL:GaussLobatto::Scale operation.......................
+DEAL:GaussLobatto::Shift operation.......................
+DEAL:GaussLobatto::Multiply by scalar.......................
+DEAL:GaussLobatto::Multiply by polynomial.......................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 26
+DEAL:GaussLobatto::Scale operation..........................
+DEAL:GaussLobatto::Shift operation..........................
+DEAL:GaussLobatto::Multiply by scalar..........................
+DEAL:GaussLobatto::Multiply by polynomial..........................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 29
+DEAL:GaussLobatto::Scale operation.............................
+DEAL:GaussLobatto::Shift operation.............................
+DEAL:GaussLobatto::Multiply by scalar.............................
+DEAL:GaussLobatto::Multiply by polynomial.............................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 32
+DEAL:GaussLobatto::Scale operation................................
+DEAL:GaussLobatto::Shift operation................................
+DEAL:GaussLobatto::Multiply by scalar................................
+DEAL:GaussLobatto::Multiply by polynomial................................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 35
+DEAL:GaussLobatto::Scale operation...................................
+DEAL:GaussLobatto::Shift operation...................................
+DEAL:GaussLobatto::Multiply by scalar...................................
+DEAL:GaussLobatto::Multiply by polynomial...................................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 38
+DEAL:GaussLobatto::Scale operation......................................
+DEAL:GaussLobatto::Shift operation......................................
+DEAL:GaussLobatto::Multiply by scalar......................................
+DEAL:GaussLobatto::Multiply by polynomial......................................
+DEAL:GaussLobatto::
diff --git a/tests/base/polynomial_lagrange_order.cc b/tests/base/polynomial_lagrange_order.cc
new file mode 100644 (file)
index 0000000..733e57e
--- /dev/null
@@ -0,0 +1,133 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2011 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.
+//
+//-----------------------------------------------------------------------------
+
+// Similar to polyomial_lagrange, but test Lagrange interpolation for high
+// order with tighter tolerances, in particular the effect of stability of the
+// polynomial evaluation at random points
+
+#include "../tests.h"
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+
+
+using namespace Polynomials;
+
+
+void check_interpolation (const std::vector<Polynomial<double> >& p,
+                         const std::vector<Point<1> >& x)
+{
+  for (unsigned int i=0;i<p.size();++i)
+    {
+      deallog << i;
+      for (unsigned int k=0;k<x.size();++k)
+       {
+         deallog << '.';
+         const double y = p[i].value(x[k](0));
+         if (i == k)
+           {
+             if (std::fabs(y-1.) > 1e-13)
+               deallog << "Error1  lg y=" << std::log10(std::fabs(y-1.))
+                       << std::endl;
+           }
+         else
+           {
+             if (std::fabs(y) > 1e-13)
+               deallog << "Error0  lg y=" << std::log10(std::fabs(y))
+                       << std::endl;
+           }
+       }
+      deallog << std::endl;
+    }
+}
+
+
+
+void check_constant (const std::vector<Polynomial<double> >& p)
+{
+                               // check whether the sum of all polynomials in
+                               // the basis gives one for a given point
+  deallog << "Representation of one at random points";
+  for (unsigned int j=0; j<12; ++j)
+    {
+      double x = (double)rand()/RAND_MAX;
+      double value = 0;
+      for (unsigned int i=0;i<p.size();++i)
+       {
+         value += p[i].value(x);
+       }
+      deallog << ".";
+      if (std::fabs (1.-value) > 1e-13)
+       deallog << "Error1  lg y=" << std::log10(std::fabs(1.-value))
+               << std::endl;
+    }
+  deallog << std::endl;
+
+}
+
+
+
+void
+check_poly (const Quadrature<1>& q)
+{
+  deallog << "Points: " << q.size() << std::endl;
+  std::vector<Polynomial<double> > p = generate_complete_Lagrange_basis(q.get_points());
+  check_interpolation(p, q.get_points());
+  check_constant (p);
+  deallog << std::endl;
+}
+
+
+
+void
+check_lge (unsigned int n)
+{
+  deallog << "Points: " << n+1 << std::endl;
+  std::vector<Polynomial<double> > p = LagrangeEquidistant::generate_complete_basis(n);
+  std::vector<Point<1> > x(n+1);
+  const double h = 1./n;
+  for (unsigned int i=0;i<=n;++i)
+    x[i](0) = h*i;
+  check_interpolation(p, x);
+  check_constant(p);
+  deallog << std::endl;
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("polynomial_lagrange_order/output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog.push("LagrangeEquidistant");
+  for (unsigned i=8; i<18; i+=2)
+    check_lge (i);
+  deallog.pop();
+  deallog << std::endl;
+
+                               // Lagrange elements on GL points have good
+                               // conditioning, so test to some very high
+                               // orders
+  deallog.push("GaussLobatto");
+  for (unsigned i=8; i<40; i+=3)
+    check_poly (QGaussLobatto<1>(i+1));
+  deallog.pop();
+}
diff --git a/tests/base/polynomial_lagrange_order/cmp/generic b/tests/base/polynomial_lagrange_order/cmp/generic
new file mode 100644 (file)
index 0000000..9b34ee8
--- /dev/null
@@ -0,0 +1,379 @@
+
+DEAL:LagrangeEquidistant::Points: 9
+DEAL:LagrangeEquidistant::0.........
+DEAL:LagrangeEquidistant::1.........
+DEAL:LagrangeEquidistant::2.........
+DEAL:LagrangeEquidistant::3.........
+DEAL:LagrangeEquidistant::4.........
+DEAL:LagrangeEquidistant::5.........
+DEAL:LagrangeEquidistant::6.........
+DEAL:LagrangeEquidistant::7.........
+DEAL:LagrangeEquidistant::8.........
+DEAL:LagrangeEquidistant::Representation of one at random points............
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 11
+DEAL:LagrangeEquidistant::0...........
+DEAL:LagrangeEquidistant::1...........
+DEAL:LagrangeEquidistant::2...........
+DEAL:LagrangeEquidistant::3...........
+DEAL:LagrangeEquidistant::4...........
+DEAL:LagrangeEquidistant::5...........
+DEAL:LagrangeEquidistant::6...........
+DEAL:LagrangeEquidistant::7...........
+DEAL:LagrangeEquidistant::8...........
+DEAL:LagrangeEquidistant::9...........
+DEAL:LagrangeEquidistant::10...........
+DEAL:LagrangeEquidistant::Representation of one at random points............
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 13
+DEAL:LagrangeEquidistant::0.............
+DEAL:LagrangeEquidistant::1.............
+DEAL:LagrangeEquidistant::2.............
+DEAL:LagrangeEquidistant::3.............
+DEAL:LagrangeEquidistant::4.............
+DEAL:LagrangeEquidistant::5.............
+DEAL:LagrangeEquidistant::6.............
+DEAL:LagrangeEquidistant::7.............
+DEAL:LagrangeEquidistant::8.............
+DEAL:LagrangeEquidistant::9.............
+DEAL:LagrangeEquidistant::10.............
+DEAL:LagrangeEquidistant::11.............
+DEAL:LagrangeEquidistant::12.............
+DEAL:LagrangeEquidistant::Representation of one at random points............
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 15
+DEAL:LagrangeEquidistant::0...............
+DEAL:LagrangeEquidistant::1...............
+DEAL:LagrangeEquidistant::2...............
+DEAL:LagrangeEquidistant::3...............
+DEAL:LagrangeEquidistant::4...............
+DEAL:LagrangeEquidistant::5...............
+DEAL:LagrangeEquidistant::6...............
+DEAL:LagrangeEquidistant::7...............
+DEAL:LagrangeEquidistant::8...............
+DEAL:LagrangeEquidistant::9...............
+DEAL:LagrangeEquidistant::10...............
+DEAL:LagrangeEquidistant::11...............
+DEAL:LagrangeEquidistant::12...............
+DEAL:LagrangeEquidistant::13...............
+DEAL:LagrangeEquidistant::14...............
+DEAL:LagrangeEquidistant::Representation of one at random points............
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 17
+DEAL:LagrangeEquidistant::0.................
+DEAL:LagrangeEquidistant::1.................
+DEAL:LagrangeEquidistant::2.................
+DEAL:LagrangeEquidistant::3.................
+DEAL:LagrangeEquidistant::4.................
+DEAL:LagrangeEquidistant::5.................
+DEAL:LagrangeEquidistant::6.................
+DEAL:LagrangeEquidistant::7.................
+DEAL:LagrangeEquidistant::8.................
+DEAL:LagrangeEquidistant::9.................
+DEAL:LagrangeEquidistant::10.................
+DEAL:LagrangeEquidistant::11.................
+DEAL:LagrangeEquidistant::12.................
+DEAL:LagrangeEquidistant::13.................
+DEAL:LagrangeEquidistant::14.................
+DEAL:LagrangeEquidistant::15.................
+DEAL:LagrangeEquidistant::16.................
+DEAL:LagrangeEquidistant::Representation of one at random points............
+DEAL:LagrangeEquidistant::
+DEAL::
+DEAL:GaussLobatto::Points: 9
+DEAL:GaussLobatto::0.........
+DEAL:GaussLobatto::1.........
+DEAL:GaussLobatto::2.........
+DEAL:GaussLobatto::3.........
+DEAL:GaussLobatto::4.........
+DEAL:GaussLobatto::5.........
+DEAL:GaussLobatto::6.........
+DEAL:GaussLobatto::7.........
+DEAL:GaussLobatto::8.........
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 12
+DEAL:GaussLobatto::0............
+DEAL:GaussLobatto::1............
+DEAL:GaussLobatto::2............
+DEAL:GaussLobatto::3............
+DEAL:GaussLobatto::4............
+DEAL:GaussLobatto::5............
+DEAL:GaussLobatto::6............
+DEAL:GaussLobatto::7............
+DEAL:GaussLobatto::8............
+DEAL:GaussLobatto::9............
+DEAL:GaussLobatto::10............
+DEAL:GaussLobatto::11............
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 15
+DEAL:GaussLobatto::0...............
+DEAL:GaussLobatto::1...............
+DEAL:GaussLobatto::2...............
+DEAL:GaussLobatto::3...............
+DEAL:GaussLobatto::4...............
+DEAL:GaussLobatto::5...............
+DEAL:GaussLobatto::6...............
+DEAL:GaussLobatto::7...............
+DEAL:GaussLobatto::8...............
+DEAL:GaussLobatto::9...............
+DEAL:GaussLobatto::10...............
+DEAL:GaussLobatto::11...............
+DEAL:GaussLobatto::12...............
+DEAL:GaussLobatto::13...............
+DEAL:GaussLobatto::14...............
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 18
+DEAL:GaussLobatto::0..................
+DEAL:GaussLobatto::1..................
+DEAL:GaussLobatto::2..................
+DEAL:GaussLobatto::3..................
+DEAL:GaussLobatto::4..................
+DEAL:GaussLobatto::5..................
+DEAL:GaussLobatto::6..................
+DEAL:GaussLobatto::7..................
+DEAL:GaussLobatto::8..................
+DEAL:GaussLobatto::9..................
+DEAL:GaussLobatto::10..................
+DEAL:GaussLobatto::11..................
+DEAL:GaussLobatto::12..................
+DEAL:GaussLobatto::13..................
+DEAL:GaussLobatto::14..................
+DEAL:GaussLobatto::15..................
+DEAL:GaussLobatto::16..................
+DEAL:GaussLobatto::17..................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 21
+DEAL:GaussLobatto::0.....................
+DEAL:GaussLobatto::1.....................
+DEAL:GaussLobatto::2.....................
+DEAL:GaussLobatto::3.....................
+DEAL:GaussLobatto::4.....................
+DEAL:GaussLobatto::5.....................
+DEAL:GaussLobatto::6.....................
+DEAL:GaussLobatto::7.....................
+DEAL:GaussLobatto::8.....................
+DEAL:GaussLobatto::9.....................
+DEAL:GaussLobatto::10.....................
+DEAL:GaussLobatto::11.....................
+DEAL:GaussLobatto::12.....................
+DEAL:GaussLobatto::13.....................
+DEAL:GaussLobatto::14.....................
+DEAL:GaussLobatto::15.....................
+DEAL:GaussLobatto::16.....................
+DEAL:GaussLobatto::17.....................
+DEAL:GaussLobatto::18.....................
+DEAL:GaussLobatto::19.....................
+DEAL:GaussLobatto::20.....................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 24
+DEAL:GaussLobatto::0........................
+DEAL:GaussLobatto::1........................
+DEAL:GaussLobatto::2........................
+DEAL:GaussLobatto::3........................
+DEAL:GaussLobatto::4........................
+DEAL:GaussLobatto::5........................
+DEAL:GaussLobatto::6........................
+DEAL:GaussLobatto::7........................
+DEAL:GaussLobatto::8........................
+DEAL:GaussLobatto::9........................
+DEAL:GaussLobatto::10........................
+DEAL:GaussLobatto::11........................
+DEAL:GaussLobatto::12........................
+DEAL:GaussLobatto::13........................
+DEAL:GaussLobatto::14........................
+DEAL:GaussLobatto::15........................
+DEAL:GaussLobatto::16........................
+DEAL:GaussLobatto::17........................
+DEAL:GaussLobatto::18........................
+DEAL:GaussLobatto::19........................
+DEAL:GaussLobatto::20........................
+DEAL:GaussLobatto::21........................
+DEAL:GaussLobatto::22........................
+DEAL:GaussLobatto::23........................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 27
+DEAL:GaussLobatto::0...........................
+DEAL:GaussLobatto::1...........................
+DEAL:GaussLobatto::2...........................
+DEAL:GaussLobatto::3...........................
+DEAL:GaussLobatto::4...........................
+DEAL:GaussLobatto::5...........................
+DEAL:GaussLobatto::6...........................
+DEAL:GaussLobatto::7...........................
+DEAL:GaussLobatto::8...........................
+DEAL:GaussLobatto::9...........................
+DEAL:GaussLobatto::10...........................
+DEAL:GaussLobatto::11...........................
+DEAL:GaussLobatto::12...........................
+DEAL:GaussLobatto::13...........................
+DEAL:GaussLobatto::14...........................
+DEAL:GaussLobatto::15...........................
+DEAL:GaussLobatto::16...........................
+DEAL:GaussLobatto::17...........................
+DEAL:GaussLobatto::18...........................
+DEAL:GaussLobatto::19...........................
+DEAL:GaussLobatto::20...........................
+DEAL:GaussLobatto::21...........................
+DEAL:GaussLobatto::22...........................
+DEAL:GaussLobatto::23...........................
+DEAL:GaussLobatto::24...........................
+DEAL:GaussLobatto::25...........................
+DEAL:GaussLobatto::26...........................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 30
+DEAL:GaussLobatto::0..............................
+DEAL:GaussLobatto::1..............................
+DEAL:GaussLobatto::2..............................
+DEAL:GaussLobatto::3..............................
+DEAL:GaussLobatto::4..............................
+DEAL:GaussLobatto::5..............................
+DEAL:GaussLobatto::6..............................
+DEAL:GaussLobatto::7..............................
+DEAL:GaussLobatto::8..............................
+DEAL:GaussLobatto::9..............................
+DEAL:GaussLobatto::10..............................
+DEAL:GaussLobatto::11..............................
+DEAL:GaussLobatto::12..............................
+DEAL:GaussLobatto::13..............................
+DEAL:GaussLobatto::14..............................
+DEAL:GaussLobatto::15..............................
+DEAL:GaussLobatto::16..............................
+DEAL:GaussLobatto::17..............................
+DEAL:GaussLobatto::18..............................
+DEAL:GaussLobatto::19..............................
+DEAL:GaussLobatto::20..............................
+DEAL:GaussLobatto::21..............................
+DEAL:GaussLobatto::22..............................
+DEAL:GaussLobatto::23..............................
+DEAL:GaussLobatto::24..............................
+DEAL:GaussLobatto::25..............................
+DEAL:GaussLobatto::26..............................
+DEAL:GaussLobatto::27..............................
+DEAL:GaussLobatto::28..............................
+DEAL:GaussLobatto::29..............................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 33
+DEAL:GaussLobatto::0.................................
+DEAL:GaussLobatto::1.................................
+DEAL:GaussLobatto::2.................................
+DEAL:GaussLobatto::3.................................
+DEAL:GaussLobatto::4.................................
+DEAL:GaussLobatto::5.................................
+DEAL:GaussLobatto::6.................................
+DEAL:GaussLobatto::7.................................
+DEAL:GaussLobatto::8.................................
+DEAL:GaussLobatto::9.................................
+DEAL:GaussLobatto::10.................................
+DEAL:GaussLobatto::11.................................
+DEAL:GaussLobatto::12.................................
+DEAL:GaussLobatto::13.................................
+DEAL:GaussLobatto::14.................................
+DEAL:GaussLobatto::15.................................
+DEAL:GaussLobatto::16.................................
+DEAL:GaussLobatto::17.................................
+DEAL:GaussLobatto::18.................................
+DEAL:GaussLobatto::19.................................
+DEAL:GaussLobatto::20.................................
+DEAL:GaussLobatto::21.................................
+DEAL:GaussLobatto::22.................................
+DEAL:GaussLobatto::23.................................
+DEAL:GaussLobatto::24.................................
+DEAL:GaussLobatto::25.................................
+DEAL:GaussLobatto::26.................................
+DEAL:GaussLobatto::27.................................
+DEAL:GaussLobatto::28.................................
+DEAL:GaussLobatto::29.................................
+DEAL:GaussLobatto::30.................................
+DEAL:GaussLobatto::31.................................
+DEAL:GaussLobatto::32.................................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 36
+DEAL:GaussLobatto::0....................................
+DEAL:GaussLobatto::1....................................
+DEAL:GaussLobatto::2....................................
+DEAL:GaussLobatto::3....................................
+DEAL:GaussLobatto::4....................................
+DEAL:GaussLobatto::5....................................
+DEAL:GaussLobatto::6....................................
+DEAL:GaussLobatto::7....................................
+DEAL:GaussLobatto::8....................................
+DEAL:GaussLobatto::9....................................
+DEAL:GaussLobatto::10....................................
+DEAL:GaussLobatto::11....................................
+DEAL:GaussLobatto::12....................................
+DEAL:GaussLobatto::13....................................
+DEAL:GaussLobatto::14....................................
+DEAL:GaussLobatto::15....................................
+DEAL:GaussLobatto::16....................................
+DEAL:GaussLobatto::17....................................
+DEAL:GaussLobatto::18....................................
+DEAL:GaussLobatto::19....................................
+DEAL:GaussLobatto::20....................................
+DEAL:GaussLobatto::21....................................
+DEAL:GaussLobatto::22....................................
+DEAL:GaussLobatto::23....................................
+DEAL:GaussLobatto::24....................................
+DEAL:GaussLobatto::25....................................
+DEAL:GaussLobatto::26....................................
+DEAL:GaussLobatto::27....................................
+DEAL:GaussLobatto::28....................................
+DEAL:GaussLobatto::29....................................
+DEAL:GaussLobatto::30....................................
+DEAL:GaussLobatto::31....................................
+DEAL:GaussLobatto::32....................................
+DEAL:GaussLobatto::33....................................
+DEAL:GaussLobatto::34....................................
+DEAL:GaussLobatto::35....................................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 39
+DEAL:GaussLobatto::0.......................................
+DEAL:GaussLobatto::1.......................................
+DEAL:GaussLobatto::2.......................................
+DEAL:GaussLobatto::3.......................................
+DEAL:GaussLobatto::4.......................................
+DEAL:GaussLobatto::5.......................................
+DEAL:GaussLobatto::6.......................................
+DEAL:GaussLobatto::7.......................................
+DEAL:GaussLobatto::8.......................................
+DEAL:GaussLobatto::9.......................................
+DEAL:GaussLobatto::10.......................................
+DEAL:GaussLobatto::11.......................................
+DEAL:GaussLobatto::12.......................................
+DEAL:GaussLobatto::13.......................................
+DEAL:GaussLobatto::14.......................................
+DEAL:GaussLobatto::15.......................................
+DEAL:GaussLobatto::16.......................................
+DEAL:GaussLobatto::17.......................................
+DEAL:GaussLobatto::18.......................................
+DEAL:GaussLobatto::19.......................................
+DEAL:GaussLobatto::20.......................................
+DEAL:GaussLobatto::21.......................................
+DEAL:GaussLobatto::22.......................................
+DEAL:GaussLobatto::23.......................................
+DEAL:GaussLobatto::24.......................................
+DEAL:GaussLobatto::25.......................................
+DEAL:GaussLobatto::26.......................................
+DEAL:GaussLobatto::27.......................................
+DEAL:GaussLobatto::28.......................................
+DEAL:GaussLobatto::29.......................................
+DEAL:GaussLobatto::30.......................................
+DEAL:GaussLobatto::31.......................................
+DEAL:GaussLobatto::32.......................................
+DEAL:GaussLobatto::33.......................................
+DEAL:GaussLobatto::34.......................................
+DEAL:GaussLobatto::35.......................................
+DEAL:GaussLobatto::36.......................................
+DEAL:GaussLobatto::37.......................................
+DEAL:GaussLobatto::38.......................................
+DEAL:GaussLobatto::Representation of one at random points............
+DEAL:GaussLobatto::
diff --git a/tests/base/polynomial_lagrange_product.cc b/tests/base/polynomial_lagrange_product.cc
new file mode 100644 (file)
index 0000000..b2948ff
--- /dev/null
@@ -0,0 +1,112 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2011 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.
+//
+//-----------------------------------------------------------------------------
+
+// Tests correctness of values and derivatives for polynomials derived from
+// Lagrange product form
+
+#include "../tests.h"
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+
+
+using namespace Polynomials;
+
+
+void check_derivatives (const std::vector<Polynomial<double> >& p,
+                       const unsigned int                      n_deriv)
+{
+                               // check whether the values and derivatives
+                               // are evaluated correctly some randomly
+                               // generated points. compare with a polynomial
+                               // that is not in product form (we get the
+                               // expanded form by adding a dummy polynomial;
+                               // addition of polynomials destroys the
+                               // product form in the current implementation)
+  deallog << "Representation of derivatives up to order " << n_deriv-1 << std::endl;
+  std::vector<double> values(n_deriv), values_ref(n_deriv);
+  Monomial<double> zero (0,0);
+  for (unsigned int j=0; j<p.size(); ++j)
+    {
+      double x = (double)rand()/RAND_MAX;
+      p[j].value (x, values);
+      Polynomial<double> q = p[j];
+      q += zero;
+      q.value (x, values_ref);
+      for (unsigned int i=0; i<n_deriv; ++i)
+       {
+         deallog << ".";
+         if (std::fabs (values[i]-values_ref[i]) >
+             std::max(1e-11,1e-11*std::fabs(values[i])))
+           deallog << "Error deriv" << i << "  lg y="
+                   << std::log10(std::fabs(values[i]-values_ref[i]))
+                   << ", is: " << values[i] << ", should be: " << values_ref[i]
+                   << std::endl;
+       }
+    }
+  deallog << std::endl;
+
+}
+
+
+
+void
+check_poly (const Quadrature<1>& q)
+{
+  deallog << "Points: " << q.size() << std::endl;
+  std::vector<Polynomial<double> > p = generate_complete_Lagrange_basis(q.get_points());
+  for (unsigned int i=1; i<6; ++i)
+    check_derivatives(p, i);
+  deallog << std::endl;
+}
+
+
+
+void
+check_lge (unsigned int n)
+{
+  deallog << "Points: " << n+1 << std::endl;
+  std::vector<Polynomial<double> > p = LagrangeEquidistant::generate_complete_basis(n);
+  for (unsigned int i=1; i<6; ++i)
+    check_derivatives(p, i);
+  deallog << std::endl;
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("polynomial_lagrange_product/output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog.push("LagrangeEquidistant");
+  for (unsigned i=1; i<8; i+=2)
+    check_lge (i);
+  deallog.pop();
+  deallog << std::endl;
+
+                               // Lagrange elements on GL points have good
+                               // conditioning, so test to some very high
+                               // orders
+  deallog.push("GaussLobatto");
+  for (unsigned i=1; i<8; i+=2)
+    check_poly (QGaussLobatto<1>(i+1));
+  deallog.pop();
+}
diff --git a/tests/base/polynomial_lagrange_product/cmp/generic b/tests/base/polynomial_lagrange_product/cmp/generic
new file mode 100644 (file)
index 0000000..08dfcd9
--- /dev/null
@@ -0,0 +1,98 @@
+
+DEAL:LagrangeEquidistant::Points: 2
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 0
+DEAL:LagrangeEquidistant::..
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 1
+DEAL:LagrangeEquidistant::....
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 2
+DEAL:LagrangeEquidistant::......
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 3
+DEAL:LagrangeEquidistant::........
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 4
+DEAL:LagrangeEquidistant::..........
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 4
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 0
+DEAL:LagrangeEquidistant::....
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 1
+DEAL:LagrangeEquidistant::........
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 2
+DEAL:LagrangeEquidistant::............
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 3
+DEAL:LagrangeEquidistant::................
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 4
+DEAL:LagrangeEquidistant::....................
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 6
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 0
+DEAL:LagrangeEquidistant::......
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 1
+DEAL:LagrangeEquidistant::............
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 2
+DEAL:LagrangeEquidistant::..................
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 3
+DEAL:LagrangeEquidistant::........................
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 4
+DEAL:LagrangeEquidistant::..............................
+DEAL:LagrangeEquidistant::
+DEAL:LagrangeEquidistant::Points: 8
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 0
+DEAL:LagrangeEquidistant::........
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 1
+DEAL:LagrangeEquidistant::................
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 2
+DEAL:LagrangeEquidistant::........................
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 3
+DEAL:LagrangeEquidistant::................................
+DEAL:LagrangeEquidistant::Representation of derivatives up to order 4
+DEAL:LagrangeEquidistant::........................................
+DEAL:LagrangeEquidistant::
+DEAL::
+DEAL:GaussLobatto::Points: 2
+DEAL:GaussLobatto::Representation of derivatives up to order 0
+DEAL:GaussLobatto::..
+DEAL:GaussLobatto::Representation of derivatives up to order 1
+DEAL:GaussLobatto::....
+DEAL:GaussLobatto::Representation of derivatives up to order 2
+DEAL:GaussLobatto::......
+DEAL:GaussLobatto::Representation of derivatives up to order 3
+DEAL:GaussLobatto::........
+DEAL:GaussLobatto::Representation of derivatives up to order 4
+DEAL:GaussLobatto::..........
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 4
+DEAL:GaussLobatto::Representation of derivatives up to order 0
+DEAL:GaussLobatto::....
+DEAL:GaussLobatto::Representation of derivatives up to order 1
+DEAL:GaussLobatto::........
+DEAL:GaussLobatto::Representation of derivatives up to order 2
+DEAL:GaussLobatto::............
+DEAL:GaussLobatto::Representation of derivatives up to order 3
+DEAL:GaussLobatto::................
+DEAL:GaussLobatto::Representation of derivatives up to order 4
+DEAL:GaussLobatto::....................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 6
+DEAL:GaussLobatto::Representation of derivatives up to order 0
+DEAL:GaussLobatto::......
+DEAL:GaussLobatto::Representation of derivatives up to order 1
+DEAL:GaussLobatto::............
+DEAL:GaussLobatto::Representation of derivatives up to order 2
+DEAL:GaussLobatto::..................
+DEAL:GaussLobatto::Representation of derivatives up to order 3
+DEAL:GaussLobatto::........................
+DEAL:GaussLobatto::Representation of derivatives up to order 4
+DEAL:GaussLobatto::..............................
+DEAL:GaussLobatto::
+DEAL:GaussLobatto::Points: 8
+DEAL:GaussLobatto::Representation of derivatives up to order 0
+DEAL:GaussLobatto::........
+DEAL:GaussLobatto::Representation of derivatives up to order 1
+DEAL:GaussLobatto::................
+DEAL:GaussLobatto::Representation of derivatives up to order 2
+DEAL:GaussLobatto::........................
+DEAL:GaussLobatto::Representation of derivatives up to order 3
+DEAL:GaussLobatto::................................
+DEAL:GaussLobatto::Representation of derivatives up to order 4
+DEAL:GaussLobatto::........................................
+DEAL:GaussLobatto::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.