* 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
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.
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;
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;
}
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;
}
}
}
- polynomial.value(y, values);
+ polynomial.value(y, values_size, values);
// change sign if necessary
for (unsigned int j=1; j<values_size; j+=2)
#include <deal.II/base/exceptions.h>
#include <deal.II/base/table.h>
+#include <deal.II/base/std_cxx11/array.h>
+
DEAL_II_NAMESPACE_OPEN
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),
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)
{
{
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)
{
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)
if (d2==x) ++derivative;
grad_grads[i][d1][d2]
- *= v(x,indices[x])[derivative];
+ *= v[indices[x]][x][derivative];
}
}
if (d3==x) ++derivative;
third_derivatives[i][d1][d2][d3]
- *= v(x,indices[x])[derivative];
+ *= v[indices[x]][x][derivative];
}
}
if (d4==x) ++derivative;
fourth_derivatives[i][d1][d2][d3][d4]
- *= v(x,indices[x])[derivative];
+ *= v[indices[x]][x][derivative];
}
}
}
projects / dealii.git / commitdiff