From d566de88673fb5151c5f78d8bf5f0b2073d70787 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 12 Oct 2016 14:09:20 -0600 Subject: [PATCH] Provide paths for static and dynamic memory. --- include/deal.II/base/polynomial.h | 8 ++--- include/deal.II/base/polynomials_piecewise.h | 8 ++--- source/base/polynomial.cc | 30 +++++++++--------- source/base/polynomials_piecewise.cc | 16 +++++----- source/base/tensor_product_polynomials.cc | 32 ++++++++++++++------ 5 files changed, 53 insertions(+), 41 deletions(-) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 22c60b963b..8e4f22a162 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -105,16 +105,16 @@ namespace Polynomials /** * Return the values and the derivatives of the Polynomial at point - * x. values[i], i=0,...,values_size-1 includes the + * x. values[i], i=0,...,n_derivatives 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. + * determined by @p n_derivatives and @p values has to provide sufficient + * space for @p n_derivatives + 1 values. * * This function uses the Horner scheme for numerical stability of the * evaluation. */ void value (const number x, - const unsigned int values_size, + const unsigned int n_derivatives, number *values) const; /** diff --git a/include/deal.II/base/polynomials_piecewise.h b/include/deal.II/base/polynomials_piecewise.h index e070407690..a79933c0ed 100644 --- a/include/deal.II/base/polynomials_piecewise.h +++ b/include/deal.II/base/polynomials_piecewise.h @@ -98,10 +98,10 @@ namespace Polynomials /** * Return the values and the derivatives of the Polynomial at point - * x. values[i], i=0,...,values_size-1 includes the + * x. values[i], i=0,...,n_derivatives 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. + * determined by @p n_derivatives and @p values has to provide sufficient + * space for @p n_derivatives + 1 values. * * Note that all the derivatives evaluate to zero at the border between * intervals (assuming exact arithmetics) in the interior of the unit @@ -112,7 +112,7 @@ namespace Polynomials * make sense. */ void value (const number x, - const unsigned int values_size, + const unsigned int n_derivatives, number *values) const; /** diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index 32b03950eb..cedc1f4936 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -103,7 +103,7 @@ namespace Polynomials { Assert (values.size() > 0, ExcZero()); - value(x,values.size(),&values[0]); + value(x,values.size()-1,&values[0]); } @@ -111,10 +111,10 @@ namespace Polynomials template void Polynomial::value (const number x, - const unsigned int values_size, + const unsigned int n_derivatives, number *values) const { - Assert(values_size > 0, ExcZero()); + Assert(n_derivatives >= 0, ExcZero()); // evaluate Lagrange polynomial and derivatives if (in_lagrange_product_form == true) @@ -123,11 +123,11 @@ namespace Polynomials // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like // automatic differentiation does. const unsigned int n_supp = lagrange_support_points.size(); - switch (values_size) + switch (n_derivatives) { default: values[0] = 1; - for (unsigned int d=1; d0; --k) + for (unsigned int k=n_derivatives; k>0; --k) values[k] = (values[k] * v + values[k-1]); values[0] *= v; } @@ -153,7 +153,7 @@ namespace Polynomials // p^(n)(x)/k! into the actual form of the derivative { number k_faculty = 1; - for (unsigned int k=0; k(k+1); @@ -161,10 +161,10 @@ namespace Polynomials } break; - // manually implement size 1 (values only), size 2 (value + first - // derivative), and size 3 (up to second derivative) since they + // manually implement case 0 (values only), case 1 (value + first + // derivative), and case 2 (up to second derivative) since they // might be called often. then, we can unroll the loop. - case 1: + case 0: values[0] = 1; for (unsigned int i=0; im} are // necessarily zero, as they differentiate the polynomial more often than // the highest power is - const unsigned int min_valuessize_m=std::min(values_size, m); + const unsigned int min_valuessize_m=std::min(n_derivatives+1, m); for (unsigned int j=0; j=static_cast(j); --k) @@ -237,7 +237,7 @@ namespace Polynomials } // fill higher derivatives by zero - for (unsigned int j=min_valuessize_m; j 0, ExcZero()); - value(x,values.size(),&values[0]); + value(x,values.size()-1,&values[0]); } @@ -55,10 +55,10 @@ namespace Polynomials template void PiecewisePolynomial::value (const number x, - const unsigned int values_size, + const unsigned int n_derivatives, number *values) const { - Assert (values_size > 0, ExcZero()); + Assert (n_derivatives >= 0, ExcZero()); // shift polynomial if necessary number y = x; @@ -72,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 - #include DEAL_II_NAMESPACE_OPEN @@ -313,6 +312,9 @@ compute (const Point &p, case 2: index = polynomial*factor + f; break; + default: + ExcNotImplemented(); + break; } const unsigned int i = index_map_inverse[index]; @@ -326,17 +328,27 @@ compute (const Point &p, // Compute the values (and derivatives, if // necessary) of all polynomials at this - // 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 + // evaluation point. To avoid expensive memory allocation, + // we use a small amount of memory on the stack, and store the + // result in an array (that has enough // fields for any evaluation of values and - // derivatives, up to the 4th derivative). - std_cxx11::array, dim> *v = - (std_cxx11::array, dim> *) - alloca(sizeof(std_cxx11::array, dim>) * polynomials.size()); + // derivatives, up to the 4th derivative, for up to 20 polynomials). + // If someone uses a larger number of + // polynomials, we need to allocate more memory on the heap. + std_cxx11::array, dim> *v; + std_cxx11::array, dim>, 20> small_array; + std::vector, dim> > large_array; + + const unsigned int n_polynomials = polynomials.size(); + if (n_polynomials > 20) + { + large_array.resize(n_polynomials); + v = &large_array[0]; + } + else + v = &small_array[0]; - for (unsigned int i=0; i