#include <deal.II/base/exceptions.h>
#include <deal.II/base/table.h>
+#include <boost/container/small_vector.hpp>
#include <array>
DEAL_II_NAMESPACE_OPEN
std::vector<Tensor<3,dim> > &third_derivatives,
std::vector<Tensor<4,dim> > &fourth_derivatives) const
{
+ Assert(dim<=3, ExcNotImplemented());
Assert (values.size()==n_tensor_pols || values.size()==0,
ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
Assert (grads.size()==n_tensor_pols || grads.size()==0,
update_3rd_derivatives = (third_derivatives.size()==n_tensor_pols),
update_4th_derivatives = (fourth_derivatives.size()==n_tensor_pols);
- // check how many
- // values/derivatives we have to
- // compute
+ // check how many values/derivatives we have to compute
unsigned int n_values_and_derivatives = 0;
if (update_values)
n_values_and_derivatives = 1;
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;
- default:
- AssertThrow(false, ExcNotImplemented());
- }
-
- const unsigned int i = index_map_inverse[index];
- values[i] *= value;
- }
- }
- }
- return;
- }
-
-
- // Compute the values (and derivatives, if
- // necessary) of all polynomials at this
- // 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, for up to 20 polynomials).
- // If someone uses a larger number of
- // polynomials, we need to allocate more memory on the heap.
- std::array<std::array<double,5>, dim> *v;
- std::array<std::array<std::array<double,5>, dim>, 20> small_array;
- std::vector<std::array<std::array<double,5>, dim> > large_array;
-
+ // Compute the values (and derivatives, if necessary) of all 1D polynomials
+ // at this evaluation point. We need to compute dim*n_polynomials
+ // evaluations, involving an evaluation of each polynomial for each
+ // coordinate direction. Once we have those values, we perform the
+ // multiplications for the tensor product in the arbitrary dimension.
const unsigned int n_polynomials = polynomials.size();
- if (n_polynomials > 20)
- {
- large_array.resize(n_polynomials);
- v = &large_array[0];
- }
+ boost::container::small_vector<std::array<std::array<double,5>,dim>, 20> values_1d(n_polynomials);
+ if (n_values_and_derivatives == 1)
+ for (unsigned int i=0; i<n_polynomials; ++i)
+ for (unsigned int d=0; d<dim; ++d)
+ values_1d[i][d][0] = polynomials[i].value(p(d));
else
- v = &small_array[0];
-
- 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]);
-
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- {
- // first get the
- // one-dimensional indices of
- // this particular tensor
- // product polynomial
- unsigned int indices[dim];
- compute_index (i, indices);
-
- if (update_values)
- {
- values[i] = 1;
- for (unsigned int x=0; x<dim; ++x)
- values[i] *= v[indices[x]][x][0];
- }
-
- if (update_grads)
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<n_polynomials; ++i)
+ for (unsigned d=0; d<dim; ++d)
+ polynomials[i].value(p(d), n_values_and_derivatives, &values_1d[i][d][0]);
+
+ unsigned int indices[3];
+ unsigned int ind=0;
+ for (indices[2]=0; indices[2]<(dim>2?n_polynomials:1); ++indices[2])
+ for (indices[1]=0; indices[1]<(dim>1?n_polynomials:1); ++indices[1])
+ if (n_values_and_derivatives == 1)
+ for (indices[0]=0; indices[0]<n_polynomials; ++indices[0], ++ind)
{
- grads[i][d] = 1.;
- for (unsigned int x=0; x<dim; ++x)
- grads[i][d] *= v[indices[x]][x][d==x];
+ double value = values_1d[indices[0]][0][0];
+ for (unsigned int d=1; d<dim; ++d)
+ value *= values_1d[indices[d]][d][0];
+ values[index_map_inverse[ind]] = value;
}
+ else
+ for (indices[0]=0; indices[0]<n_polynomials; ++indices[0], ++ind)
+ {
+ unsigned int i = index_map_inverse[ind];
- if (update_grad_grads)
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int d2=0; d2<dim; ++d2)
- {
- grad_grads[i][d1][d2] = 1.;
- for (unsigned int x=0; x<dim; ++x)
- {
- unsigned int derivative=0;
- if (d1==x) ++derivative;
- if (d2==x) ++derivative;
+ if (update_values)
+ {
+ double value = values_1d[indices[0]][0][0];
+ for (unsigned int x=1; x<dim; ++x)
+ value *= values_1d[indices[x]][x][0];
+ values[i] = value;
+ }
- grad_grads[i][d1][d2]
- *= v[indices[x]][x][derivative];
+ if (update_grads)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ double grad = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ grad *= values_1d[indices[x]][x][(d==x)?1:0];
+ grads[i][d] = grad;
}
- }
- if (update_3rd_derivatives)
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int d2=0; d2<dim; ++d2)
- for (unsigned int d3=0; d3<dim; ++d3)
- {
- third_derivatives[i][d1][d2][d3] = 1.;
- for (unsigned int x=0; x<dim; ++x)
+ if (update_grad_grads)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
{
- unsigned int derivative=0;
- if (d1==x) ++derivative;
- if (d2==x) ++derivative;
- if (d3==x) ++derivative;
-
- third_derivatives[i][d1][d2][d3]
- *= v[indices[x]][x][derivative];
+ double der2 = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x) ++derivative;
+ if (d2==x) ++derivative;
+
+ der2 *= values_1d[indices[x]][x][derivative];
+ }
+ grad_grads[i][d1][d2] = der2;
}
- }
- if (update_4th_derivatives)
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int d2=0; d2<dim; ++d2)
- for (unsigned int d3=0; d3<dim; ++d3)
- for (unsigned int d4=0; d4<dim; ++d4)
- {
- fourth_derivatives[i][d1][d2][d3][d4] = 1.;
- for (unsigned int x=0; x<dim; ++x)
+ if (update_3rd_derivatives)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ for (unsigned int d3=0; d3<dim; ++d3)
{
- unsigned int derivative=0;
- if (d1==x) ++derivative;
- if (d2==x) ++derivative;
- if (d3==x) ++derivative;
- if (d4==x) ++derivative;
-
- fourth_derivatives[i][d1][d2][d3][d4]
- *= v[indices[x]][x][derivative];
+ double der3 = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x) ++derivative;
+ if (d2==x) ++derivative;
+ if (d3==x) ++derivative;
+
+ der3 *= values_1d[indices[x]][x][derivative];
+ }
+ third_derivatives[i][d1][d2][d3] = der3;
}
- }
- }
+
+ if (update_4th_derivatives)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ for (unsigned int d3=0; d3<dim; ++d3)
+ for (unsigned int d4=0; d4<dim; ++d4)
+ {
+ double der4 = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x) ++derivative;
+ if (d2==x) ++derivative;
+ if (d3==x) ++derivative;
+ if (d4==x) ++derivative;
+
+ der4 *= values_1d[indices[x]][x][derivative];
+ }
+ fourth_derivatives[i][d1][d2][d3][d4] = der4;
+ }
+ }
}