--- /dev/null
+//-----------------------------------------------------------------------------
+// $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();
+}
--- /dev/null
+
+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::
--- /dev/null
+//-----------------------------------------------------------------------------
+// $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();
+}
--- /dev/null
+
+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::
--- /dev/null
+//-----------------------------------------------------------------------------
+// $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();
+}
--- /dev/null
+
+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::