From: kronbichler Date: Tue, 26 Jul 2011 12:58:28 +0000 (+0000) Subject: Polynomial evaluation at high orders: old implementation based on form a_0 + a_1... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9cf8537c6f4817564b91490f2f75024d8cebe972;p=dealii-svn.git Polynomial evaluation at high orders: old implementation based on form a_0 + a_1 x^1 + ... + a_n x^n produces bad results. git-svn-id: https://svn.dealii.org/trunk@23959 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/base/polynomial_lagrange_ops.cc b/tests/base/polynomial_lagrange_ops.cc new file mode 100644 index 0000000000..0494cddb23 --- /dev/null +++ b/tests/base/polynomial_lagrange_ops.cc @@ -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 +#include +#include + +#include +#include +#include + + +using namespace Polynomials; + + +void check_scale (const std::vector > &p) +{ + deallog << "Scale operation"; + for (unsigned int i=0; i 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 > &p) +{ + deallog << "Shift operation"; + for (unsigned int i=0; i 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 > &p) +{ + deallog << "Multiply by scalar"; + for (unsigned int i=0; i 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 > &p) +{ + deallog << "Multiply by polynomial"; + for (unsigned int i=0; i 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 > &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 zero(0, 0.); + for (unsigned int i=0; i 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 > &p) +{ + if (p.size() > 6) + return; + deallog << "Multiply by polynomial expanding"; + Monomial zero(0, 0.); + for (unsigned int i=0; i 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_"< > 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 > 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 index 0000000000..d3b0a5b517 --- /dev/null +++ b/tests/base/polynomial_lagrange_ops/cmp/generic @@ -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 index 0000000000..733e57e1fb --- /dev/null +++ b/tests/base/polynomial_lagrange_order.cc @@ -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 +#include +#include + +#include +#include +#include + + +using namespace Polynomials; + + +void check_interpolation (const std::vector >& p, + const std::vector >& x) +{ + for (unsigned int i=0;i 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 >& 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 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 > 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 > p = LagrangeEquidistant::generate_complete_basis(n); + std::vector > 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 index 0000000000..9b34ee8ef2 --- /dev/null +++ b/tests/base/polynomial_lagrange_order/cmp/generic @@ -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 index 0000000000..b2948ff789 --- /dev/null +++ b/tests/base/polynomial_lagrange_product.cc @@ -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 +#include +#include + +#include +#include +#include + + +using namespace Polynomials; + + +void check_derivatives (const std::vector >& 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 values(n_deriv), values_ref(n_deriv); + Monomial zero (0,0); + for (unsigned int j=0; j q = p[j]; + q += zero; + q.value (x, values_ref); + for (unsigned int i=0; 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 > 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 > 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 index 0000000000..08dfcd9c77 --- /dev/null +++ b/tests/base/polynomial_lagrange_product/cmp/generic @@ -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::