From 7424f9f91b0dd566c9ade9cbebcf2e9b9fffc887 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 23 May 2022 20:45:20 +0200 Subject: [PATCH] Provide function to evaluate polynomial on array of positions --- include/deal.II/base/polynomial.h | 133 +++++++++++++----- .../matrix_free/tensor_product_kernels.h | 116 +++++++-------- 2 files changed, 158 insertions(+), 91 deletions(-) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 24e9a63418..984e120c15 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -24,6 +24,7 @@ #include #include +#include #include #include @@ -142,6 +143,24 @@ namespace Polynomials 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 + void + values_of_array(const std::array &points, + const unsigned int n_derivatives, + std::array * values) const; + /** * Degree of the polynomial. This is the degree reflected by the number of * coefficients provided by the constructor. Leading non-zero coefficients @@ -827,6 +846,21 @@ namespace Polynomials Polynomial::value(const Number2 x, const unsigned int n_derivatives, Number2 * values) const + { + values_of_array(std::array{{x}}, + n_derivatives, + reinterpret_cast *>(values)); + } + + + + template + template + inline void + Polynomial::values_of_array( + const std::array &x, + const unsigned int n_derivatives, + std::array * values) const { // evaluate Lagrange polynomial and derivatives if (in_lagrange_product_form == true) @@ -839,12 +873,16 @@ namespace Polynomials 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 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 @@ -853,8 +891,10 @@ namespace Polynomials // 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 @@ -868,7 +908,8 @@ namespace Polynomials 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(k + 1); } } @@ -881,46 +922,60 @@ namespace Polynomials // compiler with the pointer aliasing analysis. case 0: { - Number2 value = 1.; + std::array 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 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 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(2) * weight * second[e]; } - values[0] = weight * value; - values[1] = weight * derivative; - values[2] = static_cast(2) * weight * second; break; } } @@ -931,9 +986,12 @@ namespace Polynomials // if derivatives are needed, then do it properly by the full // Horner scheme - const unsigned int m = coefficients.size(); - std::vector a(coefficients.size()); - std::copy(coefficients.begin(), coefficients.end(), a.begin()); + const unsigned int m = coefficients.size(); + std::vector> 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 @@ -943,15 +1001,18 @@ namespace Polynomials for (unsigned int j = 0; j < min_valuessize_m; ++j) { for (int k = m - 2; k >= static_cast(j); --k) - a[k] += x * a[k + 1]; - values[j] = static_cast(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(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.; } diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 3f5233b42a..50f3166800 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -20,6 +20,7 @@ #include #include +#include #include #include @@ -2961,12 +2962,14 @@ namespace internal } AssertIndexRange(n_shapes, 200); - std::array shapes; + dealii::ndarray shapes; // Evaluate 1D polynomials and their derivatives + std::array 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 @@ -2984,22 +2987,22 @@ namespace internal 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 { @@ -3010,10 +3013,10 @@ namespace internal 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) { @@ -3050,12 +3053,14 @@ namespace internal ExcDimensionMismatch(renumber.size(), values.size())); AssertIndexRange(n_shapes, 200); - std::array shapes; + dealii::ndarray shapes; // Evaluate 1D polynomials and their derivatives + std::array 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 @@ -3074,16 +3079,16 @@ namespace internal 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 @@ -3091,13 +3096,13 @@ namespace internal { 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 { @@ -3107,12 +3112,12 @@ namespace internal 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) { @@ -3149,42 +3154,43 @@ namespace internal ExcDimensionMismatch(renumber.size(), values.size())); AssertIndexRange(n_shapes, 200); - std::array shapes; + dealii::ndarray shapes; // Evaluate 1D polynomials and their derivatives + std::array 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; } } } -- 2.39.5