From 42956ac212c15d80de460db322e409065bd92401 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 14 Sep 2016 19:19:18 -0600 Subject: [PATCH] Provide shortcut and allocate on stack --- include/deal.II/base/polynomial.h | 14 ++++ include/deal.II/base/polynomials_piecewise.h | 21 ++++- source/base/polynomial.cc | 14 +++- source/base/polynomials_piecewise.cc | 20 ++++- source/base/tensor_product_polynomials.cc | 88 +++++++++++++++----- 5 files changed, 129 insertions(+), 28 deletions(-) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 01888a2292..22c60b963b 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -103,6 +103,20 @@ namespace Polynomials void value (const number x, std::vector &values) const; + /** + * Return the values and the derivatives of the Polynomial at point + * x. values[i], i=0,...,values_size-1 includes the + * ith derivative. The number of derivatives to be computed is + * determined by @p values_size and @p values has to provide sufficient + * space for @p values_size values. + * + * This function uses the Horner scheme for numerical stability of the + * evaluation. + */ + void value (const number x, + const unsigned int values_size, + number *values) const; + /** * Degree of the polynomial. This is the degree reflected by the number of * coefficients provided by the constructor. Leading non-zero coefficients diff --git a/include/deal.II/base/polynomials_piecewise.h b/include/deal.II/base/polynomials_piecewise.h index 3746c1e48d..e070407690 100644 --- a/include/deal.II/base/polynomials_piecewise.h +++ b/include/deal.II/base/polynomials_piecewise.h @@ -83,7 +83,7 @@ namespace Polynomials * Return the values and the derivatives of the Polynomial at point * x. values[i], i=0,...,values.size()-1 includes the * ith derivative. The number of derivatives to be computed is - * thus determined by the size of the array passed. + * thus determined by the size of the vector passed. * * Note that all the derivatives evaluate to zero at the border between * intervals (assuming exact arithmetics) in the interior of the unit @@ -96,6 +96,25 @@ namespace Polynomials void value (const number x, std::vector &values) const; + /** + * Return the values and the derivatives of the Polynomial at point + * x. values[i], i=0,...,values_size-1 includes the + * ith derivative.The number of derivatives to be computed is + * determined by @p values_size and @p values has to provide sufficient + * space for @p values_size values. + * + * Note that all the derivatives evaluate to zero at the border between + * intervals (assuming exact arithmetics) in the interior of the unit + * interval, as there is no unique gradient value in that case for a + * piecewise polynomial. This is not always desired (e.g., when evaluating + * jumps of gradients on the element boundary), but it is the user's + * responsibility to avoid evaluation at these points when it does not + * make sense. + */ + void value (const number x, + const unsigned int values_size, + number *values) const; + /** * Degree of the polynomial. This is the degree of the underlying base * polynomial. diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index cdea89d303..32b03950eb 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -102,7 +102,19 @@ namespace Polynomials std::vector &values) const { Assert (values.size() > 0, ExcZero()); - const unsigned int values_size=values.size(); + + value(x,values.size(),&values[0]); + } + + + + template + void + Polynomial::value (const number x, + const unsigned int values_size, + number *values) const + { + Assert(values_size > 0, ExcZero()); // evaluate Lagrange polynomial and derivatives if (in_lagrange_product_form == true) diff --git a/source/base/polynomials_piecewise.cc b/source/base/polynomials_piecewise.cc index 25656864f9..614b4b3d85 100644 --- a/source/base/polynomials_piecewise.cc +++ b/source/base/polynomials_piecewise.cc @@ -46,7 +46,19 @@ namespace Polynomials std::vector &values) const { Assert (values.size() > 0, ExcZero()); - const unsigned int values_size=values.size(); + + value(x,values.size(),&values[0]); + } + + + + template + void + PiecewisePolynomial::value (const number x, + const unsigned int values_size, + number *values) const + { + Assert (values_size > 0, ExcZero()); // shift polynomial if necessary number y = x; @@ -60,7 +72,7 @@ namespace Polynomials const double offset = step * interval; if (xoffset+step+step) { - for (unsigned int k=0; koffset+step) { - for (unsigned int k=0; k #include +#include + DEAL_II_NAMESPACE_OPEN @@ -254,7 +256,7 @@ compute (const Point &p, Assert (fourth_derivatives.size()==n_tensor_pols|| fourth_derivatives.size()==0, ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0)); - const bool update_values = (values.size() == n_tensor_pols), + const bool update_values = (values.size()==n_tensor_pols), update_grads = (grads.size()==n_tensor_pols), update_grad_grads = (grad_grads.size()==n_tensor_pols), update_3rd_derivatives = (third_derivatives.size()==n_tensor_pols), @@ -275,26 +277,68 @@ compute (const Point &p, if (update_4th_derivatives) n_values_and_derivatives = 5; + // Provide a shortcut if only values are requested. For this case usually the + // temporary memory allocation below does not pay off. Also we loop over all + // tensor polynomials in a peculiar way to avoid all dynamic memory allocation. + // We need to evaluate the polynomial value n_polynomials*dim times, and + // need to compute n_polynomials^dim tensor polynomials. Therefore we can split + // the loop over the tensor polynomials into one loop over n_polynomials, + // evaluate this polynomial for each dimension and multiply it with the + // associated n_polynomials^(dim-1) tensor polynomials before moving to the + // next polynomial. + if (n_values_and_derivatives == 1) + { + const unsigned int n_polynomials = polynomials.size(); + const unsigned int factor = n_tensor_pols/n_polynomials; + + std::fill(values.begin(),values.end(),1.0); + + for (unsigned int polynomial=0; polynomial (that has enough + // evaluation point. To avoid expensive memory allocations, + // use alloca to allocate a (small) amount of memory + // on the stack and store the + // result in a vector of arrays (that has enough // fields for any evaluation of values and - // derivatives, up to the 4th derivative) - Table<2,Tensor<1,5> > v(dim, polynomials.size()); - { - std::vector tmp (n_values_and_derivatives); + // derivatives, up to the 4th derivative). + std_cxx11::array, dim> *v = + (std_cxx11::array, dim> *) + alloca(sizeof(std_cxx11::array, dim>) * polynomials.size()); + + for (unsigned int i=0; i &p, { values[i] = 1; for (unsigned int x=0; x &p, { grads[i][d] = 1.; for (unsigned int x=0; x &p, if (d2==x) ++derivative; grad_grads[i][d1][d2] - *= v(x,indices[x])[derivative]; + *= v[indices[x]][x][derivative]; } } @@ -350,7 +394,7 @@ compute (const Point &p, if (d3==x) ++derivative; third_derivatives[i][d1][d2][d3] - *= v(x,indices[x])[derivative]; + *= v[indices[x]][x][derivative]; } } @@ -370,7 +414,7 @@ compute (const Point &p, if (d4==x) ++derivative; fourth_derivatives[i][d1][d2][d3][d4] - *= v(x,indices[x])[derivative]; + *= v[indices[x]][x][derivative]; } } } -- 2.39.5