#include <deal.II/base/point.h>
#include <deal.II/base/subscriptor.h>
+#include <array>
#include <memory>
#include <vector>
const unsigned int n_derivatives,
Number2 * values) const;
+ /**
+ * Similar to the function above, but evaluate the polynomials on several
+ * positions at once, as described by the array argument @p points. This
+ * function is can be faster than the other function when the same
+ * polynomial should be evaluated on several positions at once, e.g., the
+ * x,y,z coordinates of a point for tensor-product polynomials.
+ *
+ * The template type `Number2` must implement arithmetic
+ * operations such as additions or multiplication with the type
+ * `number` of the polynomial, and must be convertible from
+ * `number` by `operator=`.
+ */
+ template <std::size_t n_entries, typename Number2>
+ void
+ values_of_array(const std::array<Number2, n_entries> &points,
+ const unsigned int n_derivatives,
+ std::array<Number2, n_entries> * values) const;
+
/**
* Degree of the polynomial. This is the degree reflected by the number of
* coefficients provided by the constructor. Leading non-zero coefficients
Polynomial<number>::value(const Number2 x,
const unsigned int n_derivatives,
Number2 * values) const
+ {
+ values_of_array(std::array<Number2, 1ul>{{x}},
+ n_derivatives,
+ reinterpret_cast<std::array<Number2, 1ul> *>(values));
+ }
+
+
+
+ template <typename number>
+ template <std::size_t n_entries, typename Number2>
+ inline void
+ Polynomial<number>::values_of_array(
+ const std::array<Number2, n_entries> &x,
+ const unsigned int n_derivatives,
+ std::array<Number2, n_entries> * values) const
{
// evaluate Lagrange polynomial and derivatives
if (in_lagrange_product_form == true)
switch (n_derivatives)
{
default:
- values[0] = 1.;
- for (unsigned int d = 1; d <= n_derivatives; ++d)
- values[d] = 0.;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[0][e] = 1.;
+ for (unsigned int k = 1; k <= n_derivatives; ++k)
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[k][e] = 0.;
for (unsigned int i = 0; i < n_supp; ++i)
{
- const Number2 v = x - lagrange_support_points[i];
+ std::array<Number2, n_entries> v = x;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ v[e] -= lagrange_support_points[i];
// multiply by (x-x_i) and compute action on all derivatives,
// too (inspired from automatic differentiation: implement the
// value from the next lower derivative from the steps before,
// need to start from the highest derivative
for (unsigned int k = n_derivatives; k > 0; --k)
- values[k] = (values[k] * v + values[k - 1]);
- values[0] *= v;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[0][e] *= v[e];
}
// finally, multiply by the weight in the Lagrange
// denominator. Could be done instead of setting values[0] = 1
number k_factorial = 1;
for (unsigned int k = 0; k <= n_derivatives; ++k)
{
- values[k] *= k_factorial * weight;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[k][e] *= k_factorial * weight;
k_factorial *= static_cast<number>(k + 1);
}
}
// compiler with the pointer aliasing analysis.
case 0:
{
- Number2 value = 1.;
+ std::array<Number2, n_entries> value;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ value[e] = 1.;
for (unsigned int i = 0; i < n_supp; ++i)
- {
- const Number2 v = x - lagrange_support_points[i];
- value *= v;
- }
- values[0] = weight * value;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ value[e] *= (x[e] - lagrange_support_points[i]);
+
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[0][e] = weight * value[e];
break;
}
case 1:
{
- Number2 value = 1.;
- Number2 derivative = 0.;
+ std::array<Number2, n_entries> value, derivative = {};
+ for (unsigned int e = 0; e < n_entries; ++e)
+ value[e] = 1.;
for (unsigned int i = 0; i < n_supp; ++i)
+ for (unsigned int e = 0; e < n_entries; ++e)
+ {
+ const Number2 v = x[e] - lagrange_support_points[i];
+ derivative[e] = derivative[e] * v + value[e];
+ value[e] *= v;
+ }
+
+ for (unsigned int e = 0; e < n_entries; ++e)
{
- const Number2 v = x - lagrange_support_points[i];
- derivative = derivative * v + value;
- value *= v;
+ values[0][e] = weight * value[e];
+ values[1][e] = weight * derivative[e];
}
- values[0] = weight * value;
- values[1] = weight * derivative;
break;
}
case 2:
{
- Number2 value = 1.;
- Number2 derivative = 0.;
- Number2 second = 0.;
+ std::array<Number2, n_entries> value, derivative = {},
+ second = {};
+ for (unsigned int e = 0; e < n_entries; ++e)
+ value[e] = 1.;
for (unsigned int i = 0; i < n_supp; ++i)
+ for (unsigned int e = 0; e < n_entries; ++e)
+ {
+ const Number2 v = x[e] - lagrange_support_points[i];
+ second[e] = second[e] * v + derivative[e];
+ derivative[e] = derivative[e] * v + value[e];
+ value[e] *= v;
+ }
+
+ for (unsigned int e = 0; e < n_entries; ++e)
{
- const Number2 v = x - lagrange_support_points[i];
- second = second * v + derivative;
- derivative = derivative * v + value;
- value *= v;
+ values[0][e] = weight * value[e];
+ values[1][e] = weight * derivative[e];
+ values[2][e] = static_cast<number>(2) * weight * second[e];
}
- values[0] = weight * value;
- values[1] = weight * derivative;
- values[2] = static_cast<number>(2) * weight * second;
break;
}
}
// if derivatives are needed, then do it properly by the full
// Horner scheme
- const unsigned int m = coefficients.size();
- std::vector<Number2> a(coefficients.size());
- std::copy(coefficients.begin(), coefficients.end(), a.begin());
+ const unsigned int m = coefficients.size();
+ std::vector<std::array<Number2, n_entries>> a(coefficients.size());
+ for (unsigned int i = 0; i < coefficients.size(); ++i)
+ for (unsigned int e = 0; e < n_entries; ++e)
+ a[i][e] = coefficients[i];
+
unsigned int j_factorial = 1;
// loop over all requested derivatives. note that derivatives @p{j>m} are
for (unsigned int j = 0; j < min_valuessize_m; ++j)
{
for (int k = m - 2; k >= static_cast<int>(j); --k)
- a[k] += x * a[k + 1];
- values[j] = static_cast<number>(j_factorial) * a[j];
+ for (unsigned int e = 0; e < n_entries; ++e)
+ a[k][e] += x[e] * a[k + 1][e];
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[j][e] = static_cast<number>(j_factorial) * a[j][e];
j_factorial *= j + 1;
}
// fill higher derivatives by zero
for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
- values[j] = 0.;
+ for (unsigned int e = 0; e < n_entries; ++e)
+ values[j][e] = 0.;
}
#include <deal.II/base/config.h>
#include <deal.II/base/aligned_vector.h>
+#include <deal.II/base/ndarray.h>
#include <deal.II/base/polynomial.h>
#include <deal.II/base/utilities.h>
}
AssertIndexRange(n_shapes, 200);
- std::array<Number2, 2 * dim * 200> shapes;
+ dealii::ndarray<Number2, 200, 2, dim> shapes;
// Evaluate 1D polynomials and their derivatives
+ std::array<Number2, dim> point;
for (unsigned int d = 0; d < dim; ++d)
- for (int i = 0; i < n_shapes; ++i)
- poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
+ point[d] = p[d];
+ for (int i = 0; i < n_shapes; ++i)
+ poly[i].values_of_array(point, 1, &shapes[i][0]);
// Go through the tensor product of shape functions and interpolate
// with optimal algorithm
if (renumber.empty())
for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
{
- value += shapes[2 * i0] * values[i];
- deriv += shapes[2 * i0 + 1] * values[i];
+ value += shapes[i0][0][0] * values[i];
+ deriv += shapes[i0][1][0] * values[i];
}
else
for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
{
- value += shapes[2 * i0] * values[renumber[i]];
- deriv += shapes[2 * i0 + 1] * values[renumber[i]];
+ value += shapes[i0][0][0] * values[renumber[i]];
+ deriv += shapes[i0][1][0] * values[renumber[i]];
}
// Interpolation + derivative in y direction
if (dim > 1)
{
- value_y += value * shapes[2 * n_shapes + 2 * i1];
- deriv_x += deriv * shapes[2 * n_shapes + 2 * i1];
- deriv_y += value * shapes[2 * n_shapes + 2 * i1 + 1];
+ value_y += value * shapes[i1][0][1];
+ deriv_x += deriv * shapes[i1][0][1];
+ deriv_y += value * shapes[i1][1][1];
}
else
{
if (dim == 3)
{
// Interpolation + derivative in z direction
- result.first += value_y * shapes[4 * n_shapes + 2 * i2];
- result.second[0] += deriv_x * shapes[4 * n_shapes + 2 * i2];
- result.second[1] += deriv_y * shapes[4 * n_shapes + 2 * i2];
- result.second[2] += value_y * shapes[4 * n_shapes + 2 * i2 + 1];
+ result.first += value_y * shapes[i2][0][2];
+ result.second[0] += deriv_x * shapes[i2][0][2];
+ result.second[1] += deriv_y * shapes[i2][0][2];
+ result.second[2] += value_y * shapes[i2][1][2];
}
else if (dim == 2)
{
ExcDimensionMismatch(renumber.size(), values.size()));
AssertIndexRange(n_shapes, 200);
- std::array<Number2, 3 * dim * 200> shapes;
+ dealii::ndarray<Number2, 200, 3, dim> shapes;
// Evaluate 1D polynomials and their derivatives
+ std::array<Number2, dim> point;
for (unsigned int d = 0; d < dim; ++d)
- for (int i = 0; i < n_shapes; ++i)
- poly[i].value(p[d], 2, shapes.data() + 3 * (d * n_shapes + i));
+ point[d] = p[d];
+ for (int i = 0; i < n_shapes; ++i)
+ poly[i].values_of_array(point, 2, &shapes[i][0]);
// Go through the tensor product of shape functions and interpolate
// with optimal algorithm
if (renumber.empty())
for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
{
- value += shapes[3 * i0] * values[i];
- deriv_1 += shapes[3 * i0 + 1] * values[i];
- deriv_2 += shapes[3 * i0 + 2] * values[i];
+ value += shapes[i0][0][0] * values[i];
+ deriv_1 += shapes[i0][1][0] * values[i];
+ deriv_2 += shapes[i0][2][0] * values[i];
}
else
for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
{
- value += shapes[3 * i0] * values[renumber[i]];
- deriv_1 += shapes[3 * i0 + 1] * values[renumber[i]];
- deriv_2 += shapes[3 * i0 + 2] * values[renumber[i]];
+ value += shapes[i0][0][0] * values[renumber[i]];
+ deriv_1 += shapes[i0][1][0] * values[renumber[i]];
+ deriv_2 += shapes[i0][2][0] * values[renumber[i]];
}
// Interpolation + derivative in y direction
{
if (dim > 2)
{
- value_y += value * shapes[3 * n_shapes + 3 * i1];
- deriv_x += deriv_1 * shapes[3 * n_shapes + 3 * i1];
- deriv_y += value * shapes[3 * n_shapes + 3 * i1 + 1];
+ value_y += value * shapes[i1][0][1];
+ deriv_x += deriv_1 * shapes[i1][0][1];
+ deriv_y += value * shapes[i1][1][1];
}
- deriv_xx += deriv_2 * shapes[3 * n_shapes + 3 * i1];
- deriv_xy += deriv_1 * shapes[3 * n_shapes + 3 * i1 + 1];
- deriv_yy += value * shapes[3 * n_shapes + 3 * i1 + 2];
+ deriv_xx += deriv_2 * shapes[i1][0][1];
+ deriv_xy += deriv_1 * shapes[i1][1][1];
+ deriv_yy += value * shapes[i1][2][1];
}
else
{
if (dim == 3)
{
// Interpolation + derivative in z direction
- result[0][0] += deriv_xx * shapes[6 * n_shapes + 3 * i2];
- result[0][1] += deriv_xy * shapes[6 * n_shapes + 3 * i2];
- result[0][2] += deriv_x * shapes[6 * n_shapes + 3 * i2 + 1];
- result[1][1] += deriv_yy * shapes[6 * n_shapes + 3 * i2];
- result[1][2] += deriv_y * shapes[6 * n_shapes + 3 * i2 + 1];
- result[2][2] += value_y * shapes[6 * n_shapes + 3 * i2 + 2];
+ result[0][0] += deriv_xx * shapes[i2][0][2];
+ result[0][1] += deriv_xy * shapes[i2][0][2];
+ result[0][2] += deriv_x * shapes[i2][1][2];
+ result[1][1] += deriv_yy * shapes[i2][0][2];
+ result[1][2] += deriv_y * shapes[i2][1][2];
+ result[2][2] += value_y * shapes[i2][2][2];
}
else if (dim == 2)
{
ExcDimensionMismatch(renumber.size(), values.size()));
AssertIndexRange(n_shapes, 200);
- std::array<Number, 2 * dim * 200> shapes;
+ dealii::ndarray<Number, 200, 2, dim> shapes;
// Evaluate 1D polynomials and their derivatives
+ std::array<Number, dim> point;
for (unsigned int d = 0; d < dim; ++d)
- for (int i = 0; i < n_shapes; ++i)
- poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
+ point[d] = p[d];
+ for (int i = 0; i < n_shapes; ++i)
+ poly[i].values_of_array(point, 1, &shapes[i][0]);
// Implement the transpose of the function above
for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
{
const Number2 test_value_z =
- dim > 2 ? (value * shapes[4 * n_shapes + 2 * i2] +
- gradient[2] * shapes[4 * n_shapes + 2 * i2 + 1]) :
- value;
+ dim > 2 ?
+ (value * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
+ value;
const Number2 test_grad_x =
- dim > 2 ? gradient[0] * shapes[4 * n_shapes + 2 * i2] : gradient[0];
- const Number2 test_grad_y =
- dim > 2 ? gradient[1] * shapes[4 * n_shapes + 2 * i2] :
- (dim > 1 ? gradient[1] : Number2());
+ dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
+ const Number2 test_grad_y = dim > 2 ?
+ gradient[1] * shapes[i2][0][2] :
+ (dim > 1 ? gradient[1] : Number2());
for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
{
- const Number2 test_value_y =
- dim > 1 ? (test_value_z * shapes[2 * n_shapes + 2 * i1] +
- test_grad_y * shapes[2 * n_shapes + 2 * i1 + 1]) :
- test_value_z;
+ const Number2 test_value_y = dim > 1 ?
+ (test_value_z * shapes[i1][0][1] +
+ test_grad_y * shapes[i1][1][1]) :
+ test_value_z;
const Number2 test_grad_xy =
- dim > 1 ? test_grad_x * shapes[2 * n_shapes + 2 * i1] :
- test_grad_x;
+ dim > 1 ? test_grad_x * shapes[i1][0][1] : test_grad_x;
if (renumber.empty())
for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
- values[i] += shapes[2 * i0] * test_value_y +
- shapes[2 * i0 + 1] * test_grad_xy;
+ values[i] += shapes[i0][0][0] * test_value_y +
+ shapes[i0][1][0] * test_grad_xy;
else
for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
- values[renumber[i]] += shapes[2 * i0] * test_value_y +
- shapes[2 * i0 + 1] * test_grad_xy;
+ values[renumber[i]] += shapes[i0][0][0] * test_value_y +
+ shapes[i0][1][0] * test_grad_xy;
}
}
}