From: Martin Kronbichler Date: Fri, 21 Apr 2023 16:50:38 +0000 (+0200) Subject: Implement arbitrary order derivatives for evaluate_tensor_product X-Git-Tag: v9.5.0-rc1~297^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f94d31cb47fa60af8af0123701c8bcb770281f41;p=dealii.git Implement arbitrary order derivatives for evaluate_tensor_product --- diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 08d09b328e..1c13afb9bf 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -3268,99 +3268,195 @@ namespace internal - template - SymmetricTensor<2, dim, typename ProductTypeNoPoint::type> - evaluate_tensor_product_hessian( + /** + * This function computes derivatives of arbitrary orders in 1d, returning a + * Tensor with the respective derivative + */ + template + inline Tensor<1, 1, typename ProductTypeNoPoint::type> + evaluate_tensor_product_higher_derivatives( const std::vector> &poly, const std::vector & values, - const Point & p, + const Point<1, Number2> & p, const std::vector & renumber = {}) { - static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); - using Number3 = typename ProductTypeNoPoint::type; - // use `int` type for this variable and the loops below to inform the - // compiler that the loops below will never overflow, which allows it to - // generate more optimized code for the variable loop bounds in the - // present context const int n_shapes = poly.size(); - AssertDimension(Utilities::pow(n_shapes, dim), values.size()); + AssertDimension(n_shapes, values.size()); Assert(renumber.empty() || renumber.size() == values.size(), ExcDimensionMismatch(renumber.size(), values.size())); - AssertIndexRange(n_shapes, 200); - dealii::ndarray shapes; + std::array shapes; + Tensor<1, 1, Number3> result; + if (renumber.empty()) + for (int i = 0; i < n_shapes; ++i) + { + poly[i].value(p[0], derivative_order, shapes.data()); + result[0] += shapes[derivative_order] * values[i]; + } + else + for (int i = 0; i < n_shapes; ++i) + { + poly[i].value(p[0], derivative_order, shapes.data()); + result[0] += shapes[derivative_order] * values[renumber[i]]; + } + return result; + } + + + /** + * This function computes derivatives of arbitrary orders in 2d, returning a + * Tensor with the respective derivatives + */ + template + inline Tensor<1, + derivative_order + 1, + typename ProductTypeNoPoint::type> + evaluate_tensor_product_higher_derivatives( + const std::vector> &poly, + const std::vector & values, + const Point<2, Number2> & p, + const std::vector & renumber = {}) + { + using Number3 = typename ProductTypeNoPoint::type; + constexpr int dim = 2; + + const int n_shapes = poly.size(); + AssertDimension(Utilities::pow(n_shapes, 2), values.size()); + Assert(renumber.empty() || renumber.size() == values.size(), + ExcDimensionMismatch(renumber.size(), values.size())); + + AssertIndexRange(n_shapes, 100); + dealii::ndarray shapes; // Evaluate 1d polynomials and their derivatives std::array point; for (unsigned int d = 0; d < dim; ++d) point[d] = p[d]; for (int i = 0; i < n_shapes; ++i) - poly[i].values_of_array(point, 2, &shapes[i][0]); + poly[i].values_of_array(point, derivative_order, &shapes[i][0]); - // Go through the tensor product of shape functions and interpolate - // with optimal algorithm - SymmetricTensor<2, dim, Number3> result; - for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) + Tensor<1, derivative_order + 1, Number3> result; + for (int i1 = 0, i = 0; i1 < n_shapes; ++i1) { - Number3 value_y = {}, deriv_x = {}, deriv_y = {}, deriv_xx = {}, - deriv_xy = {}, deriv_yy = {}; - for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) - { - // Interpolation + derivative x direction - Number3 value = {}, deriv_1 = {}, deriv_2 = {}; + Tensor<1, derivative_order + 1, Number3> result_x; + if (renumber.empty()) + for (int i0 = 0; i0 < n_shapes; ++i0, ++i) + for (unsigned int d = 0; d <= derivative_order; ++d) + result_x[d] += shapes[i0][d][0] * values[i]; + else + for (int i0 = 0; i0 < n_shapes; ++i0, ++i) + for (unsigned int d = 0; d <= derivative_order; ++d) + result_x[d] += shapes[i0][d][0] * values[renumber[i]]; + + for (unsigned int d = 0; d <= derivative_order; ++d) + result[d] += shapes[i1][d][1] * result_x[derivative_order - d]; + } + return result; + } + + + + /** + * This function computes derivatives of arbitrary orders in 3d, returning a + * Tensor with the respective derivatives + */ + template + inline Tensor<1, + ((derivative_order + 1) * (derivative_order + 2)) / 2, + typename ProductTypeNoPoint::type> + evaluate_tensor_product_higher_derivatives( + const std::vector> &poly, + const std::vector & values, + const Point<3, Number2> & p, + const std::vector & renumber = {}) + { + using Number3 = typename ProductTypeNoPoint::type; + constexpr int dim = 3; + constexpr int n_derivatives = + ((derivative_order + 1) * (derivative_order + 2)) / 2; - // Distinguish the inner loop based on whether we have a - // renumbering or not + const int n_shapes = poly.size(); + AssertDimension(Utilities::pow(n_shapes, 3), values.size()); + Assert(renumber.empty() || renumber.size() == values.size(), + ExcDimensionMismatch(renumber.size(), values.size())); + + AssertIndexRange(n_shapes, 100); + dealii::ndarray shapes; + // Evaluate 1d polynomials and their derivatives + std::array point; + for (unsigned int d = 0; d < dim; ++d) + point[d] = p[d]; + for (int i = 0; i < n_shapes; ++i) + poly[i].values_of_array(point, derivative_order, &shapes[i][0]); + + Tensor<1, n_derivatives, Number3> result; + for (int i2 = 0, i = 0; i2 < n_shapes; ++i2) + { + Tensor<1, n_derivatives, Number3> result_xy; + for (int i1 = 0; i1 < n_shapes; ++i1) + { + // apply x derivatives + Tensor<1, derivative_order + 1, Number3> result_x; if (renumber.empty()) for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - { - value += shapes[i0][0][0] * values[i]; - deriv_1 += shapes[i0][1][0] * values[i]; - deriv_2 += shapes[i0][2][0] * values[i]; - } + for (unsigned int d = 0; d <= derivative_order; ++d) + result_x[d] += shapes[i0][d][0] * values[i]; else for (int i0 = 0; i0 < n_shapes; ++i0, ++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 > 1) - { - if (dim > 2) - { - 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[i1][0][1]; - deriv_xy += deriv_1 * shapes[i1][1][1]; - deriv_yy += value * shapes[i1][2][1]; - } - else - { - result[0][0] = deriv_2; - } + for (unsigned int d = 0; d <= derivative_order; ++d) + result_x[d] += shapes[i0][d][0] * values[renumber[i]]; + + // multiply by y derivatives, sorting them in upper triangular + // matrix, starting with highest global derivative order, + // decreasing the combined order of xy derivatives by one in each + // row, to be combined with z derivatives in the next step + for (unsigned int d = 0, c = 0; d <= derivative_order; ++d) + for (unsigned int e = d; e <= derivative_order; ++e, ++c) + result_xy[c] += + shapes[i1][e - d][1] * result_x[derivative_order - e]; } + + // multiply by z derivatives, starting with highest x derivative + for (unsigned int d = 0, c = 0; d <= derivative_order; ++d) + for (unsigned int e = d; e <= derivative_order; ++e, ++c) + result[c] += shapes[i2][d][2] * result_xy[c]; + } + return result; + } + + + + template + SymmetricTensor<2, dim, typename ProductTypeNoPoint::type> + evaluate_tensor_product_hessian( + const std::vector> &poly, + const std::vector & values, + const Point & p, + const std::vector & renumber = {}) + { + static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); + + const auto hessian = + evaluate_tensor_product_higher_derivatives<2>(poly, values, p, renumber); + + using Number3 = typename ProductTypeNoPoint::type; + SymmetricTensor<2, dim, Number3> result; + if (dim == 1) + result[0][0] = hessian[0]; + else if (dim >= 2) + { + // derivatives in Hessian are xx, xy, yy, xz, yz, zz, so must re-order + // them for 3D + for (unsigned int d = 0, c = 0; d < 2; ++d) + for (unsigned int e = d; e < 2; ++e, ++c) + result[d][e] = hessian[c]; if (dim == 3) { - // Interpolation + derivative in z direction - 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) - { - result[0][0] = deriv_xx; - result[1][0] = deriv_xy; - result[1][1] = deriv_yy; + for (unsigned int d = 0; d < 2; ++d) + result[d][2] = hessian[3 + d]; + result[2][2] = hessian[5]; } }