From da6205dffa384f92996195d96ebfa0d790b5bc89 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 26 Jul 2011 15:12:31 +0000 Subject: [PATCH] Now finally fix LagrangeEquidistant constructor. Add TODO for serialization. git-svn-id: https://svn.dealii.org/trunk@23964 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/base/polynomial.h | 15 ++-- deal.II/source/base/polynomial.cc | 93 ++++++++++++----------- 2 files changed, 57 insertions(+), 51 deletions(-) diff --git a/deal.II/include/deal.II/base/polynomial.h b/deal.II/include/deal.II/base/polynomial.h index 401f040412..d0d56dc493 100644 --- a/deal.II/include/deal.II/base/polynomial.h +++ b/deal.II/include/deal.II/base/polynomial.h @@ -272,7 +272,7 @@ namespace Polynomials * product (x-x_0)*(x-x_1)*...*(x-x_n)/weight, * or not. */ - bool is_lagrange_basis; + bool in_lagrange_product_form; /** * If the polynomial is in Lagrange product @@ -677,8 +677,8 @@ namespace Polynomials inline Polynomial::Polynomial () : - is_lagrange_basis (false), - lagrange_weight (1.) + in_lagrange_product_form (false), + lagrange_weight (1.) {} template @@ -699,7 +699,7 @@ namespace Polynomials { Assert (coefficients.size() > 0, ExcEmptyObject()); - if (is_lagrange_basis == false) + if (in_lagrange_product_form == false) { // Horner scheme const unsigned int m=coefficients.size(); @@ -730,8 +730,13 @@ namespace Polynomials { // forward to serialization // function in the base class. - ar & static_cast(*this); + ar & static_cast(*this); ar & coefficients; + // TODO: adjust tests for including these + // fields + //ar & in_lagrange_product_form; + //ar & lagrange_support_points; + //ar & lagrange_weight; } } diff --git a/deal.II/source/base/polynomial.cc b/deal.II/source/base/polynomial.cc index 4fafca390f..a2557639f3 100644 --- a/deal.II/source/base/polynomial.cc +++ b/deal.II/source/base/polynomial.cc @@ -47,9 +47,9 @@ namespace Polynomials template Polynomial::Polynomial (const std::vector &a) : - coefficients (a), - is_lagrange_basis (false), - lagrange_weight (1.) + coefficients (a), + in_lagrange_product_form (false), + lagrange_weight (1.) {} @@ -57,9 +57,9 @@ namespace Polynomials template Polynomial::Polynomial (const unsigned int n) : - coefficients (n+1, 0.), - is_lagrange_basis (false), - lagrange_weight (1.) + coefficients (n+1, 0.), + in_lagrange_product_form (false), + lagrange_weight (1.) {} @@ -68,7 +68,7 @@ namespace Polynomials Polynomial::Polynomial (const std::vector > &supp, const unsigned int center) : - is_lagrange_basis (true) + in_lagrange_product_form (true) { Assert (supp.size(), ExcEmptyObject()); lagrange_support_points.reserve (supp.size()-1); @@ -124,7 +124,7 @@ namespace Polynomials // evaluate Lagrange polynomial and // derivatives - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { // to compute the value and all derivatives of // a polynomial of the form @@ -270,7 +270,7 @@ namespace Polynomials // to scale (x-x_0)*(x-x_1)*...*(x-x_n), scale // support points by 1./factor and the weight // likewise - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { number inv_fact = number(1.)/factor; number accumulated_fact = 1.; @@ -303,7 +303,7 @@ namespace Polynomials Polynomial& Polynomial::operator *= (const double s) { - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { lagrange_weight *= s; return *this; @@ -323,7 +323,7 @@ namespace Polynomials { // if we are in Lagrange form, just append the // new points - if (is_lagrange_basis == true && p.is_lagrange_basis == true) + if (in_lagrange_product_form == true && p.in_lagrange_product_form == true) { lagrange_weight *= p.lagrange_weight; lagrange_support_points.insert (lagrange_support_points.end(), @@ -333,9 +333,9 @@ namespace Polynomials } // cannot retain Lagrange basis, recompute... - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { - is_lagrange_basis = false; + in_lagrange_product_form = false; lagrange_support_points.clear(); lagrange_weight = 1.; } @@ -361,9 +361,9 @@ namespace Polynomials { // Lagrange product form cannot reasonably be // retained after polynomial addition - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { - is_lagrange_basis = false; + in_lagrange_product_form = false; lagrange_support_points.clear(); lagrange_weight = 1.; } @@ -387,9 +387,9 @@ namespace Polynomials { // Lagrange product form cannot reasonably be // retained after polynomial subtraction - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { - is_lagrange_basis = false; + in_lagrange_product_form = false; lagrange_support_points.clear(); lagrange_weight = 1.; } @@ -496,7 +496,7 @@ namespace Polynomials // shift is simple for a polynomial in product // form, (x-x_0)*(x-x_1)*...*(x-x_n). just add // offset to all shifts - if (is_lagrange_basis == true) + if (in_lagrange_product_form == true) { for (unsigned int i=0; i > + generate_unit_points (const unsigned int n) { - this->coefficients.resize(n+1); - compute_coefficients(n, support_point, this->coefficients); + std::vector > points (n+1); + const double one_over_n = 1./n; + for (unsigned int k=0;k<=n;++k) + points[k](0) = static_cast(k)*one_over_n; + return points; } - else - { - // We have precomputed tables - // up to degree 3. For - // higher order, we have to - // compute by hand. + } + } - // Start with the constant one - this->coefficients.resize(1); - this->coefficients[0] = 1.; - // Then compute the Lagrange - // polynomial as the product - // of linear factors - std::vector two (2, 1.); - for (unsigned int k=0;k<=n;++k) - { - if (k != support_point) - { - two[0] = -1.*k/n; - Polynomial factor(two); - factor.scale(1.*n/(support_point - k)); - (*this) *= factor; - } - } + LagrangeEquidistant::LagrangeEquidistant (const unsigned int n, + const unsigned int support_point) + : + Polynomial (internal::LagrangeEquidistant:: + generate_unit_points (n), + support_point) + { + // For polynomial order up to 3, we have + // precomputed weights. Use these weights + // instead of the product form + if (n <= 3) + { + this->in_lagrange_product_form = false; + this->lagrange_weight = 1.; + this->lagrange_support_points.clear(); + Assert (this->coefficients.size() == n+1, ExcInternalError()); + compute_coefficients(n, support_point, this->coefficients); } } -- 2.39.5