From 42956ac212c15d80de460db322e409065bd92401 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller <rene.gassmoeller@mailbox.org> 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<number> &values) const; + /** + * Return the values and the derivatives of the Polynomial at point + * <tt>x</tt>. <tt>values[i], i=0,...,values_size-1</tt> includes the + * <tt>i</tt>th 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 * <tt>x</tt>. <tt>values[i], i=0,...,values.size()-1</tt> includes the * <tt>i</tt>th 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<number> &values) const; + /** + * Return the values and the derivatives of the Polynomial at point + * <tt>x</tt>. <tt>values[i], i=0,...,values_size-1</tt> includes the + * <tt>i</tt>th 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<number> &values) const { Assert (values.size() > 0, ExcZero()); - const unsigned int values_size=values.size(); + + value(x,values.size(),&values[0]); + } + + + + template <typename number> + void + Polynomial<number>::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<number> &values) const { Assert (values.size() > 0, ExcZero()); - const unsigned int values_size=values.size(); + + value(x,values.size(),&values[0]); + } + + + + template <typename number> + void + PiecewisePolynomial<number>::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 (x<offset || x>offset+step+step) { - for (unsigned int k=0; k<values.size(); ++k) + for (unsigned int k=0; k<values_size; ++k) values[k] = 0; return; } @@ -77,7 +89,7 @@ namespace Polynomials const double offset = step * interval; if (x<offset || x>offset+step) { - for (unsigned int k=0; k<values.size(); ++k) + for (unsigned int k=0; k<values_size; ++k) values[k] = 0; return; } @@ -100,7 +112,7 @@ namespace Polynomials } } - polynomial.value(y, values); + polynomial.value(y, values_size, values); // change sign if necessary for (unsigned int j=1; j<values_size; j+=2) diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc index 9ee229c6b2..fa382021db 100644 --- a/source/base/tensor_product_polynomials.cc +++ b/source/base/tensor_product_polynomials.cc @@ -18,6 +18,8 @@ #include <deal.II/base/exceptions.h> #include <deal.II/base/table.h> +#include <deal.II/base/std_cxx11/array.h> + DEAL_II_NAMESPACE_OPEN @@ -254,7 +256,7 @@ compute (const Point<dim> &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<dim> &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<n_polynomials; ++polynomial) + { + for (unsigned int d=0; d<dim; ++d) + { + const double value = polynomials[polynomial].value(p(d)); + + for (unsigned int f=0; f<factor; ++f) + { + unsigned int index = 0; + switch (d) + { + case 0: + index = polynomial+f*n_polynomials; + break; + case 1: + index = polynomial*n_polynomials + (f/n_polynomials)*n_polynomials*n_polynomials + f%n_polynomials; + break; + case 2: + index = polynomial*factor + f; + break; + } + + const unsigned int i = index_map_inverse[index]; + values[i] *= value; + } + } + } + return; + } + - // compute the values (and derivatives, if + // Compute the values (and derivatives, if // necessary) of all polynomials at this - // evaluation point. to avoid many - // reallocation, use one std::vector for - // polynomial evaluation and store the - // result as Tensor<1,5> (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<double> tmp (n_values_and_derivatives); + // derivatives, up to the 4th derivative). + std_cxx11::array<std_cxx11::array<double, 5>, dim> *v = + (std_cxx11::array<std_cxx11::array<double, 5>, dim> *) + alloca(sizeof(std_cxx11::array<std_cxx11::array<double, 5>, dim>) * polynomials.size()); + + for (unsigned int i=0; i<polynomials.size(); ++i) for (unsigned int d=0; d<dim; ++d) - for (unsigned int i=0; i<polynomials.size(); ++i) - { - polynomials[i].value(p(d), tmp); - for (unsigned int e=0; e<n_values_and_derivatives; ++e) - v(d,i)[e] = tmp[e]; - }; - } + polynomials[i].value(p(d), n_values_and_derivatives, &v[i][d][0]); for (unsigned int i=0; i<n_tensor_pols; ++i) { @@ -309,7 +353,7 @@ compute (const Point<dim> &p, { values[i] = 1; for (unsigned int x=0; x<dim; ++x) - values[i] *= v(x,indices[x])[0]; + values[i] *= v[indices[x]][x][0]; } if (update_grads) @@ -317,7 +361,7 @@ compute (const Point<dim> &p, { grads[i][d] = 1.; for (unsigned int x=0; x<dim; ++x) - grads[i][d] *= v(x,indices[x])[d==x]; + grads[i][d] *= v[indices[x]][x][d==x]; } if (update_grad_grads) @@ -332,7 +376,7 @@ compute (const Point<dim> &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<dim> &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<dim> &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