/**
* 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>x</tt>. <tt>values[i], i=0,...,n_derivatives</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.
+ * 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;
/**
/**
* 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>x</tt>. <tt>values[i], i=0,...,n_derivatives</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.
+ * 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
* make sense.
*/
void value (const number x,
- const unsigned int values_size,
+ const unsigned int n_derivatives,
number *values) const;
/**
{
Assert (values.size() > 0, ExcZero());
- value(x,values.size(),&values[0]);
+ value(x,values.size()-1,&values[0]);
}
template <typename number>
void
Polynomial<number>::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)
// 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; d<values_size; ++d)
+ for (unsigned int d=1; d<=n_derivatives; ++d)
values[d] = 0;
for (unsigned int i=0; i<n_supp; ++i)
{
// i.e., expand value v and derivative one). since we reuse a
// value from the next lower derivative from the steps before,
// need to start from the highest derivative
- for (unsigned int k=values_size-1; k>0; --k)
+ for (unsigned int k=n_derivatives; k>0; --k)
values[k] = (values[k] * v + values[k-1]);
values[0] *= v;
}
// p^(n)(x)/k! into the actual form of the derivative
{
number k_faculty = 1;
- for (unsigned int k=0; k<values_size; ++k)
+ for (unsigned int k=0; k<=n_derivatives; ++k)
{
values[k] *= k_faculty * lagrange_weight;
k_faculty *= static_cast<number>(k+1);
}
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; i<n_supp; ++i)
{
values[0] *= lagrange_weight;
break;
- case 2:
+ case 1:
values[0] = 1;
values[1] = 0;
for (unsigned int i=0; i<n_supp; ++i)
values[1] *= lagrange_weight;
break;
- case 3:
+ case 2:
values[0] = 1;
values[1] = 0;
values[2] = 0;
// if we only need the value, then call the other function since that is
// significantly faster (there is no need to allocate and free memory,
// which is really expensive compared to all the other operations!)
- if (values_size == 1)
+ if (n_derivatives == 0)
{
values[0] = value(x);
return;
// loop over all requested derivatives. note that derivatives @p{j>m} 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<min_valuessize_m; ++j)
{
for (int k=m-2; k>=static_cast<int>(j); --k)
}
// fill higher derivatives by zero
- for (unsigned int j=min_valuessize_m; j<values_size; ++j)
+ for (unsigned int j=min_valuessize_m; j<=n_derivatives; ++j)
values[j] = 0;
}
{
Assert (values.size() > 0, ExcZero());
- value(x,values.size(),&values[0]);
+ value(x,values.size()-1,&values[0]);
}
template <typename number>
void
PiecewisePolynomial<number>::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;
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<=n_derivatives; ++k)
values[k] = 0;
return;
}
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<=n_derivatives; ++k)
values[k] = 0;
return;
}
(interval < n_intervals-1 || derivative_change_sign == -1.)))
{
values[0] = value(x);
- for (unsigned int d=1; d<values_size; ++d)
+ for (unsigned int d=1; d<=n_derivatives; ++d)
values[d] = 0;
return;
}
}
- polynomial.value(y, values_size, values);
+ polynomial.value(y, n_derivatives, values);
// change sign if necessary
- for (unsigned int j=1; j<values_size; j+=2)
+ for (unsigned int j=1; j<=n_derivatives; j+=2)
values[j] *= derivative_change_sign;
}
#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/table.h>
-
#include <deal.II/base/std_cxx11/array.h>
DEAL_II_NAMESPACE_OPEN
case 2:
index = polynomial*factor + f;
break;
+ default:
+ ExcNotImplemented();
+ break;
}
const unsigned int i = index_map_inverse[index];
// 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<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());
+ // 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<std_cxx11::array<double,5>, dim> *v;
+ std_cxx11::array<std_cxx11::array<std_cxx11::array<double,5>, dim>, 20> small_array;
+ std::vector<std_cxx11::array<std_cxx11::array<double,5>, 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<polynomials.size(); ++i)
+ for (unsigned int i=0; i<n_polynomials; ++i)
for (unsigned int d=0; d<dim; ++d)
polynomials[i].value(p(d), n_values_and_derivatives, &v[i][d][0]);
projects / dealii.git / commitdiff