From: Maximilian Bergbauer Date: Fri, 8 Mar 2024 08:22:01 +0000 (+0100) Subject: Split tensor product points and single point kernels X-Git-Tag: v9.6.0-rc1~496^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16736%2Fhead;p=dealii.git Split tensor product points and single point kernels --- diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 77bd61c7e9..5bc48229fa 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -37,7 +37,7 @@ #include #include #include -#include +#include #include #include diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 84303226aa..e7f1426163 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -31,7 +31,7 @@ #include #include #include -#include +#include #include diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 90cb883ba2..4f2319420b 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2014,1572 +2014,6 @@ namespace internal } } - - - /** - * Struct to avoid using Tensor<1, dim, Point> in - * evaluate_tensor_product_value_and_gradient because a Point cannot be used - * within Tensor. Instead, a specialization of this struct upcasts the point - * to a Tensor<1,dim>. - */ - template - struct ProductTypeNoPoint - { - using type = typename ProductType::type; - }; - - template - struct ProductTypeNoPoint, Number2> - { - using type = typename ProductType, Number2>::type; - }; - - - - /** - * Computes the values and derivatives of the 1d polynomials @p poly at the - * specified point @p p and stores it in @p shapes. - */ - template - inline void - compute_values_of_array( - dealii::ndarray *shapes, - const std::vector> &poly, - const Point &p, - const unsigned int derivative = 1) - { - const int n_shapes = poly.size(); - - // 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, shapes[i].data()); - } - - - - /** - * Specialization of above function for dim = 0. Should not be called. - */ - template - inline void - compute_values_of_array(dealii::ndarray *, - const std::vector> &, - const Point<0, Number> &, - const unsigned int) - { - DEAL_II_ASSERT_UNREACHABLE(); - } - - - - /** - * Interpolate inner dimensions of tensor product shape functions. - */ - template - inline -#ifndef DEBUG - DEAL_II_ALWAYS_INLINE -#endif - std::array::type, - 2 + n_values> - do_interpolate_xy(const Number *values, - const std::vector &renumber, - const dealii::ndarray *shapes, - const int n_shapes_runtime, - int &i) - { - static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); - static_assert(1 <= n_values && n_values <= 2, - "Only n_values=1,2 implemented"); - - const int n_shapes = length > 0 ? length : n_shapes_runtime; - - // If n_values > 1, we want to interpolate from a second array, - // placed in the same array immediately after the main data. This - // is used to interpolate normal derivatives onto faces. - const Number *values_2 = - n_values > 1 ? - values + stride * (length > 0 ? - Utilities::pow(length, dim) : - Utilities::fixed_power(n_shapes_runtime)) : - nullptr; - using Number3 = typename ProductTypeNoPoint::type; - std::array result = {}; - for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) - { - // Interpolation + derivative x direction - std::array inner_result = {}; - - // Distinguish the inner loop based on whether we have a - // renumbering or not - if (do_renumber && !renumber.empty()) - for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - { - // gradient - inner_result[0] += - shapes[i0][1][0] * values[renumber[i] * stride]; - // values - inner_result[1] += - shapes[i0][0][0] * values[renumber[i] * stride]; - if (n_values > 1) - inner_result[2] += - shapes[i0][0][0] * values_2[renumber[i] * stride]; - } - else - for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - { - // gradient - inner_result[0] += shapes[i0][1][0] * values[i * stride]; - // values - inner_result[1] += shapes[i0][0][0] * values[i * stride]; - if (n_values > 1) - inner_result[2] += shapes[i0][0][0] * values_2[i * stride]; - } - - if (dim > 1) - { - // Interpolation + derivative in y direction - // gradient - result[0] += inner_result[0] * shapes[i1][0][1]; - result[1] += inner_result[1] * shapes[i1][1][1]; - // values - result[2] += inner_result[1] * shapes[i1][0][1]; - if (n_values > 1) - result[3] += inner_result[2] * shapes[i1][0][1]; - } - else - { - // gradient - result[0] = inner_result[0]; - // values - result[1] = inner_result[1]; - if (n_values > 1) - result[2] = inner_result[2]; - } - } - return result; - } - - - - /** - * Interpolates the values and gradients into the points specified in - * @p compute_values_of_array() with help of the precomputed @p shapes. - */ - template - inline std::array::type, - dim + n_values> - evaluate_tensor_product_value_and_gradient_shapes( - const dealii::ndarray *shapes, - const int n_shapes, - const Number *values, - const std::vector &renumber = {}) - { - static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); - static_assert(1 <= n_values && n_values <= 2, - "Only n_values=1,2 implemented"); - - using Number3 = typename ProductTypeNoPoint::type; - - std::array result = {}; - if (dim == 0) - { - // We only need the interpolation of the value and normal derivatives on - // faces of a 1d element. As the interpolation is the value at the - // point, simply set the result vector accordingly. - result[0] = values[0]; - if (n_values > 1) - result[1] = values[1 * stride]; - return result; - } - - // Go through the tensor product of shape functions and interpolate - // with optimal algorithm - for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) - { - std::array inner_result; - // Generate separate code with known loop bounds for the most common - // cases - if (n_shapes == 2) - inner_result = - do_interpolate_xy(values, renumber, shapes, n_shapes, i); - else if (n_shapes == 3) - inner_result = - do_interpolate_xy(values, renumber, shapes, n_shapes, i); - else if (n_shapes == 4) - inner_result = - do_interpolate_xy(values, renumber, shapes, n_shapes, i); - else if (n_shapes == 5) - inner_result = - do_interpolate_xy(values, renumber, shapes, n_shapes, i); - else if (n_shapes == 6) - inner_result = - do_interpolate_xy(values, renumber, shapes, n_shapes, i); - else - inner_result = - do_interpolate_xy(values, renumber, shapes, n_shapes, i); - if (dim == 3) - { - // derivative + interpolation in z direction - // gradient - result[0] += inner_result[0] * shapes[i2][0][2]; - result[1] += inner_result[1] * shapes[i2][0][2]; - result[2] += inner_result[2] * shapes[i2][1][2]; - // values - result[3] += inner_result[2] * shapes[i2][0][2]; - if (n_values > 1) - result[4] += inner_result[3] * shapes[i2][0][2]; - } - else if (dim == 2) - { - // gradient - result[0] = inner_result[0]; - result[1] = inner_result[1]; - // values - result[2] = inner_result[2]; - if (n_values > 1) - result[3] = inner_result[3]; - } - else - { - // gradient - result[0] = inner_result[0]; - // values - result[1] = inner_result[1]; - if (n_values > 1) - result[2] = inner_result[2]; - } - } - - return result; - } - - - - /** - * Specializes @p evaluate_tensor_product_value_and_gradient() for linear - * polynomials which massively reduces the necessary instructions. - */ - template - inline std::array::type, - dim + n_values> - evaluate_tensor_product_value_and_gradient_linear( - const Number *values, - const Point &p) - { - static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); - static_assert(1 <= n_values && n_values <= 2, - "Only n_values=1,2 implemented"); - - using Number3 = typename ProductTypeNoPoint::type; - - // If n_values > 1, we want to interpolate from a second array, - // placed in the same array immediately after the main data. This - // is used to interpolate normal derivatives onto faces. - - std::array result; - if (dim == 0) - { - // we only need the value on faces of a 1d element - result[0] = values[0]; - if (n_values > 1) - result[1] = values[1 * stride]; - } - else if (dim == 1) - { - // gradient - result[0] = Number3(values[stride] - values[0]); - // values - result[1] = Number3(values[0]) + p[0] * result[0]; - if (n_values > 1) - result[2] = Number3(values[2 * stride]) + - p[0] * (values[3 * stride] - values[2 * stride]); - } - else if (dim == 2) - { - const Number3 val10 = Number3(values[stride] - values[0]); - const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); - const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; - - // gradient - result[0] = val10 + p[1] * (val32 - val10); - result[1] = tmp1 - tmp0; - - // values - result[2] = tmp0 + p[1] * result[1]; - - if (n_values > 1) - { - const Number3 tmp0_2 = - Number3(values[4 * stride]) + - p[0] * (values[5 * stride] - values[4 * stride]); - const Number3 tmp1_2 = - Number3(values[6 * stride]) + - p[0] * (values[7 * stride] - values[6 * stride]); - result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2); - } - } - else if (dim == 3) - { - const Number3 val10 = Number3(values[stride] - values[0]); - const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); - const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; - const Number3 tmp10 = tmp1 - tmp0; - const Number3 tmpy0 = tmp0 + p[1] * tmp10; - - const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]); - const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]); - const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54; - const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76; - const Number3 tmp32 = tmp3 - tmp2; - const Number3 tmpy1 = tmp2 + p[1] * tmp32; - - // gradient - result[2] = tmpy1 - tmpy0; - result[1] = tmp10 + p[2] * (tmp32 - tmp10); - const Number3 tmpz0 = val10 + p[1] * (val32 - val10); - result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0); - - // value - result[3] = tmpy0 + p[2] * result[2]; - Assert(n_values == 1, ExcNotImplemented()); - } - - return result; - } - - - - /** - * Compute the polynomial interpolation of a tensor product shape function - * $\varphi_i$ given a vector of coefficients $u_i$ in the form - * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape - * functions $\varphi_i(\mathbf{x}) = - * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor - * product. The function returns a pair with the value of the interpolation - * as the first component and the gradient in reference coordinates as the - * second component. Note that for compound types (e.g. the `values` field - * begin a Point argument), the components of the gradient are - * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives - * as the first index; this is a consequence of the generic arguments in the - * function. - * - * @param poly The underlying one-dimensional polynomial basis - * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials. - * - * @param values The expansion coefficients $u_i$ of type `Number` in - * the polynomial interpolation. The coefficients can be simply `double` - * variables but e.g. also Point in case they define arithmetic - * operations with the type `Number2`. - * - * @param p The position in reference coordinates where the interpolation - * should be evaluated. - * - * @param d_linear Flag to specify whether a d-linear (linear in 1d, - * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which - * allows to unroll loops and considerably speed up evaluation. - * - * @param renumber Optional parameter to specify a renumbering in the - * coefficient vector, assuming that `values[renumber[i]]` returns - * the lexicographic (tensor product) entry of the coefficients. If the - * vector is entry, the values are assumed to be sorted lexicographically. - */ - template - inline std::pair< - typename ProductTypeNoPoint::type, - Tensor<1, dim, typename ProductTypeNoPoint::type>> - evaluate_tensor_product_value_and_gradient( - const std::vector> &poly, - const std::vector &values, - const Point &p, - const bool d_linear = false, - const std::vector &renumber = {}) - { - using Number3 = typename ProductTypeNoPoint::type; - - std::array result; - if (d_linear) - { - result = - evaluate_tensor_product_value_and_gradient_linear(values.data(), p); - } - else - { - AssertIndexRange(poly.size(), 200); - std::array, 200> shapes; - compute_values_of_array(shapes.data(), poly, p); - result = evaluate_tensor_product_value_and_gradient_shapes( - shapes.data(), poly.size(), values.data(), renumber); - } - return std::make_pair(result[dim], - Tensor<1, dim, Number3>( - ArrayView(result.data(), dim))); - } - - - - template - inline -#ifndef DEBUG - DEAL_II_ALWAYS_INLINE -#endif - typename ProductTypeNoPoint::type - do_interpolate_xy_value(const Number *values, - const std::vector &renumber, - const dealii::ndarray *shapes, - const int n_shapes_runtime, - int &i) - { - const int n_shapes = length > 0 ? length : n_shapes_runtime; - using Number3 = typename ProductTypeNoPoint::type; - Number3 result = {}; - for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) - { - // Interpolation x direction - Number3 value = {}; - - // Distinguish the inner loop based on whether we have a - // renumbering or not - if (do_renumber && !renumber.empty()) - for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - value += shapes[i0][0][0] * values[renumber[i] * stride]; - else - for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - value += shapes[i0][0][0] * values[i * stride]; - - if (dim > 1) - result += value * shapes[i1][0][1]; - else - result = value; - } - return result; - } - - - - template - inline typename ProductTypeNoPoint::type - evaluate_tensor_product_value_shapes( - const dealii::ndarray *shapes, - const int n_shapes, - const Number *values, - const std::vector &renumber = {}) - { - static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); - - // we only need the value on faces of a 1d element - if (dim == 0) - { - return values[0]; - } - - using Number3 = typename ProductTypeNoPoint::type; - - // Go through the tensor product of shape functions and interpolate - // with optimal algorithm - Number3 result = {}; - for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) - { - Number3 inner_result; - // Generate separate code with known loop bounds for the most common - // cases - if (n_shapes == 2) - inner_result = do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); - else if (n_shapes == 3) - inner_result = do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); - else if (n_shapes == 4) - inner_result = do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); - else if (n_shapes == 5) - inner_result = do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); - else if (n_shapes == 6) - inner_result = do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); - else - inner_result = do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); - if (dim == 3) - { - // Interpolation + derivative in z direction - result += inner_result * shapes[i2][0][2]; - } - else - { - result = inner_result; - } - } - - return result; - } - - - - template - inline typename ProductTypeNoPoint::type - evaluate_tensor_product_value_linear(const Number *values, - const Point &p) - { - static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); - - using Number3 = typename ProductTypeNoPoint::type; - - if (dim == 0) - { - // we only need the value on faces of a 1d element - return values[0]; - } - else if (dim == 1) - { - return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]); - } - else if (dim == 2) - { - const Number3 val10 = Number3(values[stride] - values[0]); - const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); - const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; - return tmp0 + p[1] * (tmp1 - tmp0); - } - else if (dim == 3) - { - const Number3 val10 = Number3(values[stride] - values[0]); - const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); - const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; - const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0); - - const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]); - const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]); - const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54; - const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76; - const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2); - - return tmpy0 + p[2] * (tmpy1 - tmpy0); - } - - // work around a compile error: missing return statement - return Number3(); - } - - - - template - inline typename ProductTypeNoPoint::type - evaluate_tensor_product_value( - const std::vector> &poly, - const std::vector &values, - const Point &p, - const bool d_linear = false, - const std::vector &renumber = {}) - { - typename ProductTypeNoPoint::type result; - if (d_linear) - { - result = evaluate_tensor_product_value_linear(values.data(), p); - } - else - { - AssertIndexRange(poly.size(), 200); - std::array, 200> shapes; - const int n_shapes = poly.size(); - 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, 0, shapes[i].data()); - result = evaluate_tensor_product_value_shapes( - shapes.data(), n_shapes, values.data(), renumber); - } - return result; - } - - - - /** - * 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<1, Number2> &p, - const std::vector &renumber = {}) - { - using Number3 = typename ProductTypeNoPoint::type; - - const int n_shapes = poly.size(); - AssertDimension(n_shapes, values.size()); - Assert(renumber.empty() || renumber.size() == values.size(), - ExcDimensionMismatch(renumber.size(), values.size())); - - 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, derivative_order, &shapes[i][0]); - - Tensor<1, derivative_order + 1, Number3> result; - for (int i1 = 0, i = 0; i1 < n_shapes; ++i1) - { - 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; - - 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) - 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]]; - - // 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) - { - for (unsigned int d = 0; d < 2; ++d) - result[d][2] = hessian[3 + d]; - result[2][2] = hessian[5]; - } - } - - return result; - } - - - - /** - * Test inner dimensions of tensor product shape functions and accumulate. - */ - template - inline -#ifndef DEBUG - DEAL_II_ALWAYS_INLINE -#endif - void - do_apply_test_functions_xy( - Number2 *values, - const dealii::ndarray *shapes, - const std::array &test_grads_value, - const int n_shapes_runtime, - int &i) - { - static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); - static_assert(1 <= n_values && n_values <= 2, - "Only n_values=1,2 implemented"); - - // Note that 'add' is a template argument, so the compiler will remove - // these checks - if (length > 0) - { - constexpr unsigned int array_size = length > 0 ? length : 1; - std::array shape_values_x; - std::array shape_derivs_x; - for (unsigned int j = 0; j < array_size; ++j) - { - shape_values_x[j] = shapes[j][0][0]; - shape_derivs_x[j] = shapes[j][1][0]; - } - for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1) - { - const Number2 test_value_y = - dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + - test_grads_value[1] * shapes[i1][1][1]) : - test_grads_value[2]; - const Number2 test_grad_xy = - dim > 1 ? test_grads_value[0] * shapes[i1][0][1] : - test_grads_value[0]; - Number2 test_value_y_2; - if (n_values > 1) - test_value_y_2 = dim > 1 ? - test_grads_value[3] * shapes[i1][0][1] : - test_grads_value[3]; - - Number2 *values_ptr = values + i + i1 * length; - Number2 *values_ptr_2 = - n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr; - for (int i0 = 0; i0 < length; ++i0) - { - if (add) - values_ptr[i0] += shape_values_x[i0] * test_value_y; - else - values_ptr[i0] = shape_values_x[i0] * test_value_y; - values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy; - if (n_values > 1) - { - if (add) - values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2; - else - values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2; - } - } - } - i += (dim > 1 ? length * length : length); - } - else - { - for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1) - { - const Number2 test_value_y = - dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + - test_grads_value[1] * shapes[i1][1][1]) : - test_grads_value[2]; - const Number2 test_grad_xy = - dim > 1 ? test_grads_value[0] * shapes[i1][0][1] : - test_grads_value[0]; - Number2 test_value_y_2; - if (n_values > 1) - test_value_y_2 = dim > 1 ? - test_grads_value[3] * shapes[i1][0][1] : - test_grads_value[3]; - - Number2 *values_ptr = values + i + i1 * n_shapes_runtime; - Number2 *values_ptr_2 = - n_values > 1 ? - values_ptr + Utilities::fixed_power(n_shapes_runtime) : - nullptr; - for (int i0 = 0; i0 < n_shapes_runtime; ++i0) - { - if (add) - values_ptr[i0] += shapes[i0][0][0] * test_value_y; - else - values_ptr[i0] = shapes[i0][0][0] * test_value_y; - values_ptr[i0] += shapes[i0][1][0] * test_grad_xy; - if (n_values > 1) - { - if (add) - values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2; - else - values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2; - } - } - } - i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime); - } - } - - - - /** - * Same as evaluate_tensor_product_value_and_gradient_shapes() but for - * integration. - */ - template - inline void - integrate_add_tensor_product_value_and_gradient_shapes( - const dealii::ndarray *shapes, - const int n_shapes, - const Number2 *value, - const Tensor<1, dim, Number2> &gradient, - Number2 *values) - { - static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); - static_assert(1 <= n_values && n_values <= 2, - "Only n_values=1,2 implemented"); - - // Note that 'add' is a template argument, so the compiler will remove - // these checks - if (dim == 0) - { - if (add) - values[0] += value[0]; - else - values[0] = value[0]; - if (n_values > 1) - { - if (add) - values[1] += value[1]; - else - values[1] = value[1]; - } - return; - } - - // Implement the transpose of the function above - // as in evaluate, use `int` type to produce better code in this context - std::array test_grads_value; - for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) - { - // test grad x - test_grads_value[0] = - dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0]; - // test grad y - test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] : - (dim > 1 ? gradient[1] : Number2()); - // test value z - test_grads_value[2] = - dim > 2 ? - (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) : - value[0]; - - if (n_values > 1) - test_grads_value[3] = - dim > 2 ? value[1] * shapes[i2][0][2] : value[1]; - // Generate separate code with known loop bounds for the most common - // cases - if (n_shapes == 2) - do_apply_test_functions_xy( - values, shapes, test_grads_value, n_shapes, i); - else if (n_shapes == 3) - do_apply_test_functions_xy( - values, shapes, test_grads_value, n_shapes, i); - else if (n_shapes == 4) - do_apply_test_functions_xy( - values, shapes, test_grads_value, n_shapes, i); - else if (n_shapes == 5) - do_apply_test_functions_xy( - values, shapes, test_grads_value, n_shapes, i); - else if (n_shapes == 6) - do_apply_test_functions_xy( - values, shapes, test_grads_value, n_shapes, i); - else - do_apply_test_functions_xy( - values, shapes, test_grads_value, n_shapes, i); - } - } - - - - /** - * Specializes @p integrate_add_tensor_product_value_and_gradient_shapes() for linear - * polynomials which massively reduces the necessary instructions. - */ - template - inline void - integrate_add_tensor_product_value_and_gradient_linear( - const Number2 *value, - const Tensor<1, dim, Number2> &gradient, - Number2 *values, - const Point &p) - { - static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); - static_assert(1 <= n_values && n_values <= 2, - "Only n_values=1,2 implemented"); - - // Note that 'add' is a template argument, so the compiler will remove - // these checks - if (dim == 0) - { - if (add) - values[0] += value[0]; - else - values[0] = value[0]; - if (n_values > 1) - { - if (add) - values[1] += value[1]; - else - values[1] = value[1]; - } - } - else if (dim == 1) - { - const Number2 difference = value[0] * p[0] + gradient[0]; - if (add) - { - values[0] += value[0] - difference; - values[1] += difference; - } - else - { - values[0] = value[0] - difference; - values[1] = difference; - } - if (n_values > 1) - { - const Number2 product = value[1] * p[0]; - if (add) - { - values[2] += value[1] - product; - values[3] += product; - } - else - { - values[2] = value[1] - product; - values[3] = product; - } - } - } - else if (dim == 2) - { - const Number2 test_value_y1 = value[0] * p[1] + gradient[1]; - const Number2 test_value_y0 = value[0] - test_value_y1; - const Number2 test_grad_xy1 = gradient[0] * p[1]; - const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1; - const Number2 value0 = p[0] * test_value_y0 + test_grad_xy0; - const Number2 value1 = p[0] * test_value_y1 + test_grad_xy1; - - if (add) - { - values[0] += test_value_y0 - value0; - values[1] += value0; - values[2] += test_value_y1 - value1; - values[3] += value1; - } - else - { - values[0] = test_value_y0 - value0; - values[1] = value0; - values[2] = test_value_y1 - value1; - values[3] = value1; - } - - if (n_values > 1) - { - const Number2 test_value_y1_2 = value[1] * p[1]; - const Number2 test_value_y0_2 = value[1] - test_value_y1_2; - const Number2 value0_2 = p[0] * test_value_y0_2; - const Number2 value1_2 = p[0] * test_value_y1_2; - - if (add) - { - values[4] += test_value_y0_2 - value0_2; - values[5] += value0_2; - values[6] += test_value_y1_2 - value1_2; - values[7] += value1_2; - } - else - { - values[4] = test_value_y0_2 - value0_2; - values[5] = value0_2; - values[6] = test_value_y1_2 - value1_2; - values[7] = value1_2; - } - } - } - else if (dim == 3) - { - Assert(n_values == 1, ExcNotImplemented()); - - const Number2 test_value_z1 = value[0] * p[2] + gradient[2]; - const Number2 test_value_z0 = value[0] - test_value_z1; - const Number2 test_grad_x1 = gradient[0] * p[2]; - const Number2 test_grad_x0 = gradient[0] - test_grad_x1; - const Number2 test_grad_y1 = gradient[1] * p[2]; - const Number2 test_grad_y0 = gradient[1] - test_grad_y1; - - const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0; - const Number2 test_value_y00 = test_value_z0 - test_value_y01; - const Number2 test_grad_xy01 = test_grad_x0 * p[1]; - const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01; - const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1; - const Number2 test_value_y10 = test_value_z1 - test_value_y11; - const Number2 test_grad_xy11 = test_grad_x1 * p[1]; - const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11; - - const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00; - const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01; - const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10; - const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11; - - if (add) - { - values[0] += test_value_y00 - value00; - values[1] += value00; - values[2] += test_value_y01 - value01; - values[3] += value01; - values[4] += test_value_y10 - value10; - values[5] += value10; - values[6] += test_value_y11 - value11; - values[7] += value11; - } - else - { - values[0] = test_value_y00 - value00; - values[1] = value00; - values[2] = test_value_y01 - value01; - values[3] = value01; - values[4] = test_value_y10 - value10; - values[5] = value10; - values[6] = test_value_y11 - value11; - values[7] = value11; - } - } - } - - - - /** - * Calls the correct @p integrate_add_tensor_product_value_and_gradient_...() - * function depending on if values should be added to or set and if - * polynomials are linear. - */ - template - inline void - integrate_tensor_product_value_and_gradient( - const dealii::ndarray *shapes, - const unsigned int n_shapes, - const Number2 *value, - const Tensor<1, dim, Number2> &gradient, - Number2 *values, - const Point &p, - const bool do_add) - { - if (do_add) - { - if (is_linear) - internal::integrate_add_tensor_product_value_and_gradient_linear< - dim, - Number, - Number2, - true, - n_values>(value, gradient, values, p); - else - internal::integrate_add_tensor_product_value_and_gradient_shapes< - dim, - Number, - Number2, - true, - n_values>(shapes, n_shapes, value, gradient, values); - } - else - { - if (is_linear) - internal::integrate_add_tensor_product_value_and_gradient_linear< - dim, - Number, - Number2, - false, - n_values>(value, gradient, values, p); - else - internal::integrate_add_tensor_product_value_and_gradient_shapes< - dim, - Number, - Number2, - false, - n_values>(shapes, n_shapes, value, gradient, values); - } - } - - - - /** - * Test inner dimensions of tensor product shape functions and accumulate. - */ - template - inline -#ifndef DEBUG - DEAL_II_ALWAYS_INLINE -#endif - void - do_apply_test_functions_xy_value( - Number2 *values, - const dealii::ndarray *shapes, - const Number2 &test_value, - const int n_shapes_runtime, - int &i) - { - if (length > 0) - { - constexpr unsigned int array_size = length > 0 ? length : 1; - std::array shape_values_x; - for (unsigned int i1 = 0; i1 < array_size; ++i1) - shape_values_x[i1] = shapes[i1][0][0]; - for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1) - { - const Number2 test_value_y = - dim > 1 ? test_value * shapes[i1][0][1] : test_value; - - Number2 *values_ptr = values + i + i1 * length; - for (unsigned int i0 = 0; i0 < length; ++i0) - { - if (add) - values_ptr[i0] += shape_values_x[i0] * test_value_y; - else - values_ptr[i0] = shape_values_x[i0] * test_value_y; - } - } - i += (dim > 1 ? length * length : length); - } - else - { - for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1) - { - const Number2 test_value_y = - dim > 1 ? test_value * shapes[i1][0][1] : test_value; - - Number2 *values_ptr = values + i + i1 * n_shapes_runtime; - for (int i0 = 0; i0 < n_shapes_runtime; ++i0) - { - if (add) - values_ptr[i0] += shapes[i0][0][0] * test_value_y; - else - values_ptr[i0] = shapes[i0][0][0] * test_value_y; - } - } - i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime); - } - } - - - - /** - * Same as evaluate_tensor_product_value_shapes() but for integration. - */ - template - inline void - integrate_add_tensor_product_value_shapes( - const dealii::ndarray *shapes, - const int n_shapes, - const Number2 &value, - Number2 *values) - { - static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); - - // as in evaluate, use `int` type to produce better code in this context - - if (dim == 0) - { - if (add) - values[0] += value; - else - values[0] = value; - return; - } - - // Implement the transpose of the function above - Number2 test_value; - for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) - { - // test value z - test_value = dim > 2 ? value * shapes[i2][0][2] : value; - - // Generate separate code with known loop bounds for the most common - // cases - if (n_shapes == 2) - do_apply_test_functions_xy_value( - values, shapes, test_value, n_shapes, i); - else if (n_shapes == 3) - do_apply_test_functions_xy_value( - values, shapes, test_value, n_shapes, i); - else if (n_shapes == 4) - do_apply_test_functions_xy_value( - values, shapes, test_value, n_shapes, i); - else if (n_shapes == 5) - do_apply_test_functions_xy_value( - values, shapes, test_value, n_shapes, i); - else if (n_shapes == 6) - do_apply_test_functions_xy_value( - values, shapes, test_value, n_shapes, i); - else - do_apply_test_functions_xy_value( - values, shapes, test_value, n_shapes, i); - } - } - - - - /** - * Specializes @p integrate_tensor_product_value_shapes() for linear - * polynomials which massively reduces the necessary instructions. - */ - template - inline void - integrate_add_tensor_product_value_linear(const Number2 &value, - Number2 *values, - const Point &p) - { - static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); - - if (dim == 0) - { - if (add) - values[0] += value; - else - values[0] = value; - } - else if (dim == 1) - { - const auto x0 = 1. - p[0], x1 = p[0]; - - if (add) - { - values[0] += value * x0; - values[1] += value * x1; - } - else - { - values[0] = value * x0; - values[1] = value * x1; - } - } - else if (dim == 2) - { - const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1]; - - const auto test_value_y0 = value * y0; - const auto test_value_y1 = value * y1; - - if (add) - { - values[0] += x0 * test_value_y0; - values[1] += x1 * test_value_y0; - values[2] += x0 * test_value_y1; - values[3] += x1 * test_value_y1; - } - else - { - values[0] = x0 * test_value_y0; - values[1] = x1 * test_value_y0; - values[2] = x0 * test_value_y1; - values[3] = x1 * test_value_y1; - } - } - else if (dim == 3) - { - const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1], - z0 = 1. - p[2], z1 = p[2]; - - const auto test_value_z0 = value * z0; - const auto test_value_z1 = value * z1; - - const auto test_value_y00 = test_value_z0 * y0; - const auto test_value_y01 = test_value_z0 * y1; - const auto test_value_y10 = test_value_z1 * y0; - const auto test_value_y11 = test_value_z1 * y1; - - if (add) - { - values[0] += x0 * test_value_y00; - values[1] += x1 * test_value_y00; - values[2] += x0 * test_value_y01; - values[3] += x1 * test_value_y01; - values[4] += x0 * test_value_y10; - values[5] += x1 * test_value_y10; - values[6] += x0 * test_value_y11; - values[7] += x1 * test_value_y11; - } - else - { - values[0] = x0 * test_value_y00; - values[1] = x1 * test_value_y00; - values[2] = x0 * test_value_y01; - values[3] = x1 * test_value_y01; - values[4] = x0 * test_value_y10; - values[5] = x1 * test_value_y10; - values[6] = x0 * test_value_y11; - values[7] = x1 * test_value_y11; - } - } - } - - - - /** - * Calls the correct @p integrate_add_tensor_product_value_...() - * function depending on if values should be added to or set and if - * polynomials are linear. - */ - template - inline void - integrate_tensor_product_value(const dealii::ndarray *shapes, - const unsigned int n_shapes, - const Number2 &value, - Number2 *values, - const Point &p, - const bool do_add) - { - if (do_add) - { - if (is_linear) - internal::integrate_add_tensor_product_value_linear(value, - values, - p); - else - internal::integrate_add_tensor_product_value_shapes(shapes, - n_shapes, - value, - values); - } - else - { - if (is_linear) - internal::integrate_add_tensor_product_value_linear(value, - values, - p); - else - internal::integrate_add_tensor_product_value_shapes(shapes, - n_shapes, - value, - values); - } - } - - - template inline void weight_fe_q_dofs_by_entity(const Number *weights, diff --git a/include/deal.II/matrix_free/tensor_product_point_kernels.h b/include/deal.II/matrix_free/tensor_product_point_kernels.h new file mode 100644 index 0000000000..4e4e76c0f2 --- /dev/null +++ b/include/deal.II/matrix_free/tensor_product_point_kernels.h @@ -0,0 +1,1601 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2020 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +#ifndef dealii_matrix_free_tensor_product_point_kernels_h +#define dealii_matrix_free_tensor_product_point_kernels_h + +#include + +#include +#include +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + + +namespace internal +{ + /** + * Struct to avoid using Tensor<1, dim, Point> in + * evaluate_tensor_product_value_and_gradient because a Point cannot be used + * within Tensor. Instead, a specialization of this struct upcasts the point + * to a Tensor<1,dim>. + */ + template + struct ProductTypeNoPoint + { + using type = typename ProductType::type; + }; + + template + struct ProductTypeNoPoint, Number2> + { + using type = typename ProductType, Number2>::type; + }; + + + + /** + * Computes the values and derivatives of the 1d polynomials @p poly at the + * specified point @p p and stores it in @p shapes. + */ + template + inline void + compute_values_of_array( + dealii::ndarray *shapes, + const std::vector> &poly, + const Point &p, + const unsigned int derivative = 1) + { + const int n_shapes = poly.size(); + + // 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, shapes[i].data()); + } + + + + /** + * Specialization of above function for dim = 0. Should not be called. + */ + template + inline void + compute_values_of_array(dealii::ndarray *, + const std::vector> &, + const Point<0, Number> &, + const unsigned int) + { + DEAL_II_ASSERT_UNREACHABLE(); + } + + + + /** + * Interpolate inner dimensions of tensor product shape functions. + */ + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + std::array::type, + 2 + n_values> + do_interpolate_xy(const Number *values, + const std::vector &renumber, + const dealii::ndarray *shapes, + const int n_shapes_runtime, + int &i) + { + static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); + static_assert(1 <= n_values && n_values <= 2, + "Only n_values=1,2 implemented"); + + const int n_shapes = length > 0 ? length : n_shapes_runtime; + + // If n_values > 1, we want to interpolate from a second array, + // placed in the same array immediately after the main data. This + // is used to interpolate normal derivatives onto faces. + const Number *values_2 = + n_values > 1 ? + values + stride * (length > 0 ? + Utilities::pow(length, dim) : + Utilities::fixed_power(n_shapes_runtime)) : + nullptr; + using Number3 = typename ProductTypeNoPoint::type; + std::array result = {}; + for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) + { + // Interpolation + derivative x direction + std::array inner_result = {}; + + // Distinguish the inner loop based on whether we have a + // renumbering or not + if (do_renumber && !renumber.empty()) + for (int i0 = 0; i0 < n_shapes; ++i0, ++i) + { + // gradient + inner_result[0] += + shapes[i0][1][0] * values[renumber[i] * stride]; + // values + inner_result[1] += + shapes[i0][0][0] * values[renumber[i] * stride]; + if (n_values > 1) + inner_result[2] += + shapes[i0][0][0] * values_2[renumber[i] * stride]; + } + else + for (int i0 = 0; i0 < n_shapes; ++i0, ++i) + { + // gradient + inner_result[0] += shapes[i0][1][0] * values[i * stride]; + // values + inner_result[1] += shapes[i0][0][0] * values[i * stride]; + if (n_values > 1) + inner_result[2] += shapes[i0][0][0] * values_2[i * stride]; + } + + if (dim > 1) + { + // Interpolation + derivative in y direction + // gradient + result[0] += inner_result[0] * shapes[i1][0][1]; + result[1] += inner_result[1] * shapes[i1][1][1]; + // values + result[2] += inner_result[1] * shapes[i1][0][1]; + if (n_values > 1) + result[3] += inner_result[2] * shapes[i1][0][1]; + } + else + { + // gradient + result[0] = inner_result[0]; + // values + result[1] = inner_result[1]; + if (n_values > 1) + result[2] = inner_result[2]; + } + } + return result; + } + + + + /** + * Interpolates the values and gradients into the points specified in + * @p compute_values_of_array() with help of the precomputed @p shapes. + */ + template + inline std::array::type, + dim + n_values> + evaluate_tensor_product_value_and_gradient_shapes( + const dealii::ndarray *shapes, + const int n_shapes, + const Number *values, + const std::vector &renumber = {}) + { + static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); + static_assert(1 <= n_values && n_values <= 2, + "Only n_values=1,2 implemented"); + + using Number3 = typename ProductTypeNoPoint::type; + + std::array result = {}; + if (dim == 0) + { + // We only need the interpolation of the value and normal derivatives on + // faces of a 1d element. As the interpolation is the value at the + // point, simply set the result vector accordingly. + result[0] = values[0]; + if (n_values > 1) + result[1] = values[1 * stride]; + return result; + } + + // Go through the tensor product of shape functions and interpolate + // with optimal algorithm + for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) + { + std::array inner_result; + // Generate separate code with known loop bounds for the most common + // cases + if (n_shapes == 2) + inner_result = + do_interpolate_xy(values, renumber, shapes, n_shapes, i); + else if (n_shapes == 3) + inner_result = + do_interpolate_xy(values, renumber, shapes, n_shapes, i); + else if (n_shapes == 4) + inner_result = + do_interpolate_xy(values, renumber, shapes, n_shapes, i); + else if (n_shapes == 5) + inner_result = + do_interpolate_xy(values, renumber, shapes, n_shapes, i); + else if (n_shapes == 6) + inner_result = + do_interpolate_xy(values, renumber, shapes, n_shapes, i); + else + inner_result = + do_interpolate_xy(values, renumber, shapes, n_shapes, i); + if (dim == 3) + { + // derivative + interpolation in z direction + // gradient + result[0] += inner_result[0] * shapes[i2][0][2]; + result[1] += inner_result[1] * shapes[i2][0][2]; + result[2] += inner_result[2] * shapes[i2][1][2]; + // values + result[3] += inner_result[2] * shapes[i2][0][2]; + if (n_values > 1) + result[4] += inner_result[3] * shapes[i2][0][2]; + } + else if (dim == 2) + { + // gradient + result[0] = inner_result[0]; + result[1] = inner_result[1]; + // values + result[2] = inner_result[2]; + if (n_values > 1) + result[3] = inner_result[3]; + } + else + { + // gradient + result[0] = inner_result[0]; + // values + result[1] = inner_result[1]; + if (n_values > 1) + result[2] = inner_result[2]; + } + } + + return result; + } + + + + /** + * Specializes @p evaluate_tensor_product_value_and_gradient() for linear + * polynomials which massively reduces the necessary instructions. + */ + template + inline std::array::type, + dim + n_values> + evaluate_tensor_product_value_and_gradient_linear( + const Number *values, + const Point &p) + { + static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); + static_assert(1 <= n_values && n_values <= 2, + "Only n_values=1,2 implemented"); + + using Number3 = typename ProductTypeNoPoint::type; + + // If n_values > 1, we want to interpolate from a second array, + // placed in the same array immediately after the main data. This + // is used to interpolate normal derivatives onto faces. + + std::array result; + if (dim == 0) + { + // we only need the value on faces of a 1d element + result[0] = values[0]; + if (n_values > 1) + result[1] = values[1 * stride]; + } + else if (dim == 1) + { + // gradient + result[0] = Number3(values[stride] - values[0]); + // values + result[1] = Number3(values[0]) + p[0] * result[0]; + if (n_values > 1) + result[2] = Number3(values[2 * stride]) + + p[0] * (values[3 * stride] - values[2 * stride]); + } + else if (dim == 2) + { + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); + const Number3 tmp0 = Number3(values[0]) + p[0] * val10; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; + + // gradient + result[0] = val10 + p[1] * (val32 - val10); + result[1] = tmp1 - tmp0; + + // values + result[2] = tmp0 + p[1] * result[1]; + + if (n_values > 1) + { + const Number3 tmp0_2 = + Number3(values[4 * stride]) + + p[0] * (values[5 * stride] - values[4 * stride]); + const Number3 tmp1_2 = + Number3(values[6 * stride]) + + p[0] * (values[7 * stride] - values[6 * stride]); + result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2); + } + } + else if (dim == 3) + { + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); + const Number3 tmp0 = Number3(values[0]) + p[0] * val10; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; + const Number3 tmp10 = tmp1 - tmp0; + const Number3 tmpy0 = tmp0 + p[1] * tmp10; + + const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]); + const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]); + const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54; + const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76; + const Number3 tmp32 = tmp3 - tmp2; + const Number3 tmpy1 = tmp2 + p[1] * tmp32; + + // gradient + result[2] = tmpy1 - tmpy0; + result[1] = tmp10 + p[2] * (tmp32 - tmp10); + const Number3 tmpz0 = val10 + p[1] * (val32 - val10); + result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0); + + // value + result[3] = tmpy0 + p[2] * result[2]; + Assert(n_values == 1, ExcNotImplemented()); + } + + return result; + } + + + + /** + * Compute the polynomial interpolation of a tensor product shape function + * $\varphi_i$ given a vector of coefficients $u_i$ in the form + * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape + * functions $\varphi_i(\mathbf{x}) = + * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor + * product. The function returns a pair with the value of the interpolation + * as the first component and the gradient in reference coordinates as the + * second component. Note that for compound types (e.g. the `values` field + * begin a Point argument), the components of the gradient are + * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives + * as the first index; this is a consequence of the generic arguments in the + * function. + * + * @param poly The underlying one-dimensional polynomial basis + * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials. + * + * @param values The expansion coefficients $u_i$ of type `Number` in + * the polynomial interpolation. The coefficients can be simply `double` + * variables but e.g. also Point in case they define arithmetic + * operations with the type `Number2`. + * + * @param p The position in reference coordinates where the interpolation + * should be evaluated. + * + * @param d_linear Flag to specify whether a d-linear (linear in 1d, + * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which + * allows to unroll loops and considerably speed up evaluation. + * + * @param renumber Optional parameter to specify a renumbering in the + * coefficient vector, assuming that `values[renumber[i]]` returns + * the lexicographic (tensor product) entry of the coefficients. If the + * vector is entry, the values are assumed to be sorted lexicographically. + */ + template + inline std::pair< + typename ProductTypeNoPoint::type, + Tensor<1, dim, typename ProductTypeNoPoint::type>> + evaluate_tensor_product_value_and_gradient( + const std::vector> &poly, + const std::vector &values, + const Point &p, + const bool d_linear = false, + const std::vector &renumber = {}) + { + using Number3 = typename ProductTypeNoPoint::type; + + std::array result; + if (d_linear) + { + result = + evaluate_tensor_product_value_and_gradient_linear(values.data(), p); + } + else + { + AssertIndexRange(poly.size(), 200); + std::array, 200> shapes; + compute_values_of_array(shapes.data(), poly, p); + result = evaluate_tensor_product_value_and_gradient_shapes( + shapes.data(), poly.size(), values.data(), renumber); + } + return std::make_pair(result[dim], + Tensor<1, dim, Number3>( + ArrayView(result.data(), dim))); + } + + + + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + typename ProductTypeNoPoint::type + do_interpolate_xy_value(const Number *values, + const std::vector &renumber, + const dealii::ndarray *shapes, + const int n_shapes_runtime, + int &i) + { + const int n_shapes = length > 0 ? length : n_shapes_runtime; + using Number3 = typename ProductTypeNoPoint::type; + Number3 result = {}; + for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) + { + // Interpolation x direction + Number3 value = {}; + + // Distinguish the inner loop based on whether we have a + // renumbering or not + if (do_renumber && !renumber.empty()) + for (int i0 = 0; i0 < n_shapes; ++i0, ++i) + value += shapes[i0][0][0] * values[renumber[i] * stride]; + else + for (int i0 = 0; i0 < n_shapes; ++i0, ++i) + value += shapes[i0][0][0] * values[i * stride]; + + if (dim > 1) + result += value * shapes[i1][0][1]; + else + result = value; + } + return result; + } + + + + template + inline typename ProductTypeNoPoint::type + evaluate_tensor_product_value_shapes( + const dealii::ndarray *shapes, + const int n_shapes, + const Number *values, + const std::vector &renumber = {}) + { + static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); + + // we only need the value on faces of a 1d element + if (dim == 0) + { + return values[0]; + } + + using Number3 = typename ProductTypeNoPoint::type; + + // Go through the tensor product of shape functions and interpolate + // with optimal algorithm + Number3 result = {}; + for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) + { + Number3 inner_result; + // Generate separate code with known loop bounds for the most common + // cases + if (n_shapes == 2) + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); + else if (n_shapes == 3) + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); + else if (n_shapes == 4) + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); + else if (n_shapes == 5) + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); + else if (n_shapes == 6) + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); + else + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); + if (dim == 3) + { + // Interpolation + derivative in z direction + result += inner_result * shapes[i2][0][2]; + } + else + { + result = inner_result; + } + } + + return result; + } + + + + template + inline typename ProductTypeNoPoint::type + evaluate_tensor_product_value_linear(const Number *values, + const Point &p) + { + static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); + + using Number3 = typename ProductTypeNoPoint::type; + + if (dim == 0) + { + // we only need the value on faces of a 1d element + return values[0]; + } + else if (dim == 1) + { + return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]); + } + else if (dim == 2) + { + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); + const Number3 tmp0 = Number3(values[0]) + p[0] * val10; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; + return tmp0 + p[1] * (tmp1 - tmp0); + } + else if (dim == 3) + { + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); + const Number3 tmp0 = Number3(values[0]) + p[0] * val10; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; + const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0); + + const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]); + const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]); + const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54; + const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76; + const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2); + + return tmpy0 + p[2] * (tmpy1 - tmpy0); + } + + // work around a compile error: missing return statement + return Number3(); + } + + + + template + inline typename ProductTypeNoPoint::type + evaluate_tensor_product_value( + const std::vector> &poly, + const std::vector &values, + const Point &p, + const bool d_linear = false, + const std::vector &renumber = {}) + { + typename ProductTypeNoPoint::type result; + if (d_linear) + { + result = evaluate_tensor_product_value_linear(values.data(), p); + } + else + { + AssertIndexRange(poly.size(), 200); + std::array, 200> shapes; + const int n_shapes = poly.size(); + 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, 0, shapes[i].data()); + result = evaluate_tensor_product_value_shapes( + shapes.data(), n_shapes, values.data(), renumber); + } + return result; + } + + + + /** + * 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<1, Number2> &p, + const std::vector &renumber = {}) + { + using Number3 = typename ProductTypeNoPoint::type; + + const int n_shapes = poly.size(); + AssertDimension(n_shapes, values.size()); + Assert(renumber.empty() || renumber.size() == values.size(), + ExcDimensionMismatch(renumber.size(), values.size())); + + 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, derivative_order, &shapes[i][0]); + + Tensor<1, derivative_order + 1, Number3> result; + for (int i1 = 0, i = 0; i1 < n_shapes; ++i1) + { + 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; + + 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) + 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]]; + + // 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) + { + for (unsigned int d = 0; d < 2; ++d) + result[d][2] = hessian[3 + d]; + result[2][2] = hessian[5]; + } + } + + return result; + } + + + + /** + * Test inner dimensions of tensor product shape functions and accumulate. + */ + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + void + do_apply_test_functions_xy( + Number2 *values, + const dealii::ndarray *shapes, + const std::array &test_grads_value, + const int n_shapes_runtime, + int &i) + { + static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); + static_assert(1 <= n_values && n_values <= 2, + "Only n_values=1,2 implemented"); + + // Note that 'add' is a template argument, so the compiler will remove + // these checks + if (length > 0) + { + constexpr unsigned int array_size = length > 0 ? length : 1; + std::array shape_values_x; + std::array shape_derivs_x; + for (unsigned int j = 0; j < array_size; ++j) + { + shape_values_x[j] = shapes[j][0][0]; + shape_derivs_x[j] = shapes[j][1][0]; + } + for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1) + { + const Number2 test_value_y = + dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + + test_grads_value[1] * shapes[i1][1][1]) : + test_grads_value[2]; + const Number2 test_grad_xy = + dim > 1 ? test_grads_value[0] * shapes[i1][0][1] : + test_grads_value[0]; + Number2 test_value_y_2; + if (n_values > 1) + test_value_y_2 = dim > 1 ? + test_grads_value[3] * shapes[i1][0][1] : + test_grads_value[3]; + + Number2 *values_ptr = values + i + i1 * length; + Number2 *values_ptr_2 = + n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr; + for (int i0 = 0; i0 < length; ++i0) + { + if (add) + values_ptr[i0] += shape_values_x[i0] * test_value_y; + else + values_ptr[i0] = shape_values_x[i0] * test_value_y; + values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy; + if (n_values > 1) + { + if (add) + values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2; + else + values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2; + } + } + } + i += (dim > 1 ? length * length : length); + } + else + { + for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1) + { + const Number2 test_value_y = + dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + + test_grads_value[1] * shapes[i1][1][1]) : + test_grads_value[2]; + const Number2 test_grad_xy = + dim > 1 ? test_grads_value[0] * shapes[i1][0][1] : + test_grads_value[0]; + Number2 test_value_y_2; + if (n_values > 1) + test_value_y_2 = dim > 1 ? + test_grads_value[3] * shapes[i1][0][1] : + test_grads_value[3]; + + Number2 *values_ptr = values + i + i1 * n_shapes_runtime; + Number2 *values_ptr_2 = + n_values > 1 ? + values_ptr + Utilities::fixed_power(n_shapes_runtime) : + nullptr; + for (int i0 = 0; i0 < n_shapes_runtime; ++i0) + { + if (add) + values_ptr[i0] += shapes[i0][0][0] * test_value_y; + else + values_ptr[i0] = shapes[i0][0][0] * test_value_y; + values_ptr[i0] += shapes[i0][1][0] * test_grad_xy; + if (n_values > 1) + { + if (add) + values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2; + else + values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2; + } + } + } + i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime); + } + } + + + + /** + * Same as evaluate_tensor_product_value_and_gradient_shapes() but for + * integration. + */ + template + inline void + integrate_add_tensor_product_value_and_gradient_shapes( + const dealii::ndarray *shapes, + const int n_shapes, + const Number2 *value, + const Tensor<1, dim, Number2> &gradient, + Number2 *values) + { + static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); + static_assert(1 <= n_values && n_values <= 2, + "Only n_values=1,2 implemented"); + + // Note that 'add' is a template argument, so the compiler will remove + // these checks + if (dim == 0) + { + if (add) + values[0] += value[0]; + else + values[0] = value[0]; + if (n_values > 1) + { + if (add) + values[1] += value[1]; + else + values[1] = value[1]; + } + return; + } + + // Implement the transpose of the function above + // as in evaluate, use `int` type to produce better code in this context + std::array test_grads_value; + for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) + { + // test grad x + test_grads_value[0] = + dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0]; + // test grad y + test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] : + (dim > 1 ? gradient[1] : Number2()); + // test value z + test_grads_value[2] = + dim > 2 ? + (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) : + value[0]; + + if (n_values > 1) + test_grads_value[3] = + dim > 2 ? value[1] * shapes[i2][0][2] : value[1]; + // Generate separate code with known loop bounds for the most common + // cases + if (n_shapes == 2) + do_apply_test_functions_xy( + values, shapes, test_grads_value, n_shapes, i); + else if (n_shapes == 3) + do_apply_test_functions_xy( + values, shapes, test_grads_value, n_shapes, i); + else if (n_shapes == 4) + do_apply_test_functions_xy( + values, shapes, test_grads_value, n_shapes, i); + else if (n_shapes == 5) + do_apply_test_functions_xy( + values, shapes, test_grads_value, n_shapes, i); + else if (n_shapes == 6) + do_apply_test_functions_xy( + values, shapes, test_grads_value, n_shapes, i); + else + do_apply_test_functions_xy( + values, shapes, test_grads_value, n_shapes, i); + } + } + + + + /** + * Specializes @p integrate_add_tensor_product_value_and_gradient_shapes() for linear + * polynomials which massively reduces the necessary instructions. + */ + template + inline void + integrate_add_tensor_product_value_and_gradient_linear( + const Number2 *value, + const Tensor<1, dim, Number2> &gradient, + Number2 *values, + const Point &p) + { + static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); + static_assert(1 <= n_values && n_values <= 2, + "Only n_values=1,2 implemented"); + + // Note that 'add' is a template argument, so the compiler will remove + // these checks + if (dim == 0) + { + if (add) + values[0] += value[0]; + else + values[0] = value[0]; + if (n_values > 1) + { + if (add) + values[1] += value[1]; + else + values[1] = value[1]; + } + } + else if (dim == 1) + { + const Number2 difference = value[0] * p[0] + gradient[0]; + if (add) + { + values[0] += value[0] - difference; + values[1] += difference; + } + else + { + values[0] = value[0] - difference; + values[1] = difference; + } + if (n_values > 1) + { + const Number2 product = value[1] * p[0]; + if (add) + { + values[2] += value[1] - product; + values[3] += product; + } + else + { + values[2] = value[1] - product; + values[3] = product; + } + } + } + else if (dim == 2) + { + const Number2 test_value_y1 = value[0] * p[1] + gradient[1]; + const Number2 test_value_y0 = value[0] - test_value_y1; + const Number2 test_grad_xy1 = gradient[0] * p[1]; + const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1; + const Number2 value0 = p[0] * test_value_y0 + test_grad_xy0; + const Number2 value1 = p[0] * test_value_y1 + test_grad_xy1; + + if (add) + { + values[0] += test_value_y0 - value0; + values[1] += value0; + values[2] += test_value_y1 - value1; + values[3] += value1; + } + else + { + values[0] = test_value_y0 - value0; + values[1] = value0; + values[2] = test_value_y1 - value1; + values[3] = value1; + } + + if (n_values > 1) + { + const Number2 test_value_y1_2 = value[1] * p[1]; + const Number2 test_value_y0_2 = value[1] - test_value_y1_2; + const Number2 value0_2 = p[0] * test_value_y0_2; + const Number2 value1_2 = p[0] * test_value_y1_2; + + if (add) + { + values[4] += test_value_y0_2 - value0_2; + values[5] += value0_2; + values[6] += test_value_y1_2 - value1_2; + values[7] += value1_2; + } + else + { + values[4] = test_value_y0_2 - value0_2; + values[5] = value0_2; + values[6] = test_value_y1_2 - value1_2; + values[7] = value1_2; + } + } + } + else if (dim == 3) + { + Assert(n_values == 1, ExcNotImplemented()); + + const Number2 test_value_z1 = value[0] * p[2] + gradient[2]; + const Number2 test_value_z0 = value[0] - test_value_z1; + const Number2 test_grad_x1 = gradient[0] * p[2]; + const Number2 test_grad_x0 = gradient[0] - test_grad_x1; + const Number2 test_grad_y1 = gradient[1] * p[2]; + const Number2 test_grad_y0 = gradient[1] - test_grad_y1; + + const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0; + const Number2 test_value_y00 = test_value_z0 - test_value_y01; + const Number2 test_grad_xy01 = test_grad_x0 * p[1]; + const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01; + const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1; + const Number2 test_value_y10 = test_value_z1 - test_value_y11; + const Number2 test_grad_xy11 = test_grad_x1 * p[1]; + const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11; + + const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00; + const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01; + const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10; + const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11; + + if (add) + { + values[0] += test_value_y00 - value00; + values[1] += value00; + values[2] += test_value_y01 - value01; + values[3] += value01; + values[4] += test_value_y10 - value10; + values[5] += value10; + values[6] += test_value_y11 - value11; + values[7] += value11; + } + else + { + values[0] = test_value_y00 - value00; + values[1] = value00; + values[2] = test_value_y01 - value01; + values[3] = value01; + values[4] = test_value_y10 - value10; + values[5] = value10; + values[6] = test_value_y11 - value11; + values[7] = value11; + } + } + } + + + + /** + * Calls the correct @p integrate_add_tensor_product_value_and_gradient_...() + * function depending on if values should be added to or set and if + * polynomials are linear. + */ + template + inline void + integrate_tensor_product_value_and_gradient( + const dealii::ndarray *shapes, + const unsigned int n_shapes, + const Number2 *value, + const Tensor<1, dim, Number2> &gradient, + Number2 *values, + const Point &p, + const bool do_add) + { + if (do_add) + { + if (is_linear) + internal::integrate_add_tensor_product_value_and_gradient_linear< + dim, + Number, + Number2, + true, + n_values>(value, gradient, values, p); + else + internal::integrate_add_tensor_product_value_and_gradient_shapes< + dim, + Number, + Number2, + true, + n_values>(shapes, n_shapes, value, gradient, values); + } + else + { + if (is_linear) + internal::integrate_add_tensor_product_value_and_gradient_linear< + dim, + Number, + Number2, + false, + n_values>(value, gradient, values, p); + else + internal::integrate_add_tensor_product_value_and_gradient_shapes< + dim, + Number, + Number2, + false, + n_values>(shapes, n_shapes, value, gradient, values); + } + } + + + + /** + * Test inner dimensions of tensor product shape functions and accumulate. + */ + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + void + do_apply_test_functions_xy_value( + Number2 *values, + const dealii::ndarray *shapes, + const Number2 &test_value, + const int n_shapes_runtime, + int &i) + { + if (length > 0) + { + constexpr unsigned int array_size = length > 0 ? length : 1; + std::array shape_values_x; + for (unsigned int i1 = 0; i1 < array_size; ++i1) + shape_values_x[i1] = shapes[i1][0][0]; + for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1) + { + const Number2 test_value_y = + dim > 1 ? test_value * shapes[i1][0][1] : test_value; + + Number2 *values_ptr = values + i + i1 * length; + for (unsigned int i0 = 0; i0 < length; ++i0) + { + if (add) + values_ptr[i0] += shape_values_x[i0] * test_value_y; + else + values_ptr[i0] = shape_values_x[i0] * test_value_y; + } + } + i += (dim > 1 ? length * length : length); + } + else + { + for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1) + { + const Number2 test_value_y = + dim > 1 ? test_value * shapes[i1][0][1] : test_value; + + Number2 *values_ptr = values + i + i1 * n_shapes_runtime; + for (int i0 = 0; i0 < n_shapes_runtime; ++i0) + { + if (add) + values_ptr[i0] += shapes[i0][0][0] * test_value_y; + else + values_ptr[i0] = shapes[i0][0][0] * test_value_y; + } + } + i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime); + } + } + + + + /** + * Same as evaluate_tensor_product_value_shapes() but for integration. + */ + template + inline void + integrate_add_tensor_product_value_shapes( + const dealii::ndarray *shapes, + const int n_shapes, + const Number2 &value, + Number2 *values) + { + static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); + + // as in evaluate, use `int` type to produce better code in this context + + if (dim == 0) + { + if (add) + values[0] += value; + else + values[0] = value; + return; + } + + // Implement the transpose of the function above + Number2 test_value; + for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) + { + // test value z + test_value = dim > 2 ? value * shapes[i2][0][2] : value; + + // Generate separate code with known loop bounds for the most common + // cases + if (n_shapes == 2) + do_apply_test_functions_xy_value( + values, shapes, test_value, n_shapes, i); + else if (n_shapes == 3) + do_apply_test_functions_xy_value( + values, shapes, test_value, n_shapes, i); + else if (n_shapes == 4) + do_apply_test_functions_xy_value( + values, shapes, test_value, n_shapes, i); + else if (n_shapes == 5) + do_apply_test_functions_xy_value( + values, shapes, test_value, n_shapes, i); + else if (n_shapes == 6) + do_apply_test_functions_xy_value( + values, shapes, test_value, n_shapes, i); + else + do_apply_test_functions_xy_value( + values, shapes, test_value, n_shapes, i); + } + } + + + + /** + * Specializes @p integrate_tensor_product_value_shapes() for linear + * polynomials which massively reduces the necessary instructions. + */ + template + inline void + integrate_add_tensor_product_value_linear(const Number2 &value, + Number2 *values, + const Point &p) + { + static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); + + if (dim == 0) + { + if (add) + values[0] += value; + else + values[0] = value; + } + else if (dim == 1) + { + const auto x0 = 1. - p[0], x1 = p[0]; + + if (add) + { + values[0] += value * x0; + values[1] += value * x1; + } + else + { + values[0] = value * x0; + values[1] = value * x1; + } + } + else if (dim == 2) + { + const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1]; + + const auto test_value_y0 = value * y0; + const auto test_value_y1 = value * y1; + + if (add) + { + values[0] += x0 * test_value_y0; + values[1] += x1 * test_value_y0; + values[2] += x0 * test_value_y1; + values[3] += x1 * test_value_y1; + } + else + { + values[0] = x0 * test_value_y0; + values[1] = x1 * test_value_y0; + values[2] = x0 * test_value_y1; + values[3] = x1 * test_value_y1; + } + } + else if (dim == 3) + { + const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1], + z0 = 1. - p[2], z1 = p[2]; + + const auto test_value_z0 = value * z0; + const auto test_value_z1 = value * z1; + + const auto test_value_y00 = test_value_z0 * y0; + const auto test_value_y01 = test_value_z0 * y1; + const auto test_value_y10 = test_value_z1 * y0; + const auto test_value_y11 = test_value_z1 * y1; + + if (add) + { + values[0] += x0 * test_value_y00; + values[1] += x1 * test_value_y00; + values[2] += x0 * test_value_y01; + values[3] += x1 * test_value_y01; + values[4] += x0 * test_value_y10; + values[5] += x1 * test_value_y10; + values[6] += x0 * test_value_y11; + values[7] += x1 * test_value_y11; + } + else + { + values[0] = x0 * test_value_y00; + values[1] = x1 * test_value_y00; + values[2] = x0 * test_value_y01; + values[3] = x1 * test_value_y01; + values[4] = x0 * test_value_y10; + values[5] = x1 * test_value_y10; + values[6] = x0 * test_value_y11; + values[7] = x1 * test_value_y11; + } + } + } + + + + /** + * Calls the correct @p integrate_add_tensor_product_value_...() + * function depending on if values should be added to or set and if + * polynomials are linear. + */ + template + inline void + integrate_tensor_product_value(const dealii::ndarray *shapes, + const unsigned int n_shapes, + const Number2 &value, + Number2 *values, + const Point &p, + const bool do_add) + { + if (do_add) + { + if (is_linear) + internal::integrate_add_tensor_product_value_linear(value, + values, + p); + else + internal::integrate_add_tensor_product_value_shapes(shapes, + n_shapes, + value, + values); + } + else + { + if (is_linear) + internal::integrate_add_tensor_product_value_linear(value, + values, + p); + else + internal::integrate_add_tensor_product_value_shapes(shapes, + n_shapes, + value, + values); + } + } +} // end of namespace internal + + +DEAL_II_NAMESPACE_CLOSE + +#endif