From: Martin Kronbichler Date: Thu, 17 Aug 2023 21:11:16 +0000 (+0200) Subject: Restructure tensor-product kernels via 1d interpolants X-Git-Tag: relicensing~578^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15891%2Fhead;p=dealii.git Restructure tensor-product kernels via 1d interpolants --- diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 98680df2e1..5abcf6507e 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -95,1699 +95,961 @@ namespace internal /** - * Generic evaluator framework that valuates the given shape data in general - * dimensions using the tensor product form. Depending on the particular - * layout in the matrix entries, this corresponds to a usual matrix-matrix - * product or a matrix-matrix product including some symmetries. - * - * @tparam variant Variant of evaluation used for creating template - * specializations - * @tparam dim Dimension of the function - * @tparam n_rows Number of rows in the transformation matrix, which corresponds - * to the number of 1d shape functions in the usual tensor - * contraction setting - * @tparam n_columns Number of columns in the transformation matrix, which - * corresponds to the number of 1d shape functions in the - * usual tensor contraction setting - * @tparam Number Abstract number type for input and output arrays - * @tparam Number2 Abstract number type for coefficient arrays (defaults to - * same type as the input/output arrays); must implement - * operator* with Number to be valid - */ - template - struct EvaluatorTensorProduct - {}; - - /** - * Evaluator framework for anisotropic polynomial spaces that valuates the - * given shape data in general dimensions using the tensor product form. - * - * @tparam variant Variant of evaluation used for creating template - * specializations - * @tparam dim Dimension of the function - * @tparam n_rows Number of rows in the transformation matrix, which corresponds - * to the number of 1d shape functions in the usual tensor - * contraction setting - * @tparam n_columns Number of columns in the transformation matrix, which - * corresponds to the number of 1d shape functions in the - * usual tensor contraction setting - * @tparam Number Abstract number type for input and output arrays - * @tparam Number2 Abstract number type for coefficient arrays (defaults to - * same type as the input/output arrays); must implement - * operator* with Number to be valid - * @tparam normal_dir Indicates the direction of the continuous component for the - * Raviart-Thomas space in terms of the normal onto the - * face, e.g 0 if the is in x-direction, 1 if in y-direction, and 2 if in - * z-direction. - */ - template - struct EvaluatorTensorProductAnisotropic - {}; - - - - /** - * Internal evaluator for shape function in arbitrary dimension using the - * tensor product form of the basis functions. - * - * @tparam dim Space dimension in which this class is applied - * @tparam n_rows Number of rows in the transformation matrix, which corresponds - * to the number of 1d shape functions in the usual tensor - * contraction setting - * @tparam n_columns Number of columns in the transformation matrix, which - * corresponds to the number of 1d shape functions in the - * usual tensor contraction setting - * @tparam Number Abstract number type for input and output arrays - * @tparam Number2 Abstract number type for coefficient arrays (defaults to - * same type as the input/output arrays); must implement - * operator* with Number and produce Number as an output to - * be a valid type - */ - template - struct EvaluatorTensorProduct - { - static constexpr unsigned int n_rows_of_product = - Utilities::pow(n_rows, dim); - static constexpr unsigned int n_columns_of_product = - Utilities::pow(n_columns, dim); - - /** - * Empty constructor. Does nothing. Be careful when using 'values' and - * related methods because they need to be filled with the other pointer - */ - EvaluatorTensorProduct() - : shape_values(nullptr) - , shape_gradients(nullptr) - , shape_hessians(nullptr) - {} - - /** - * Constructor, taking the data from ShapeInfo - */ - EvaluatorTensorProduct(const AlignedVector &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) - : shape_values(shape_values.begin()) - , shape_gradients(shape_gradients.begin()) - , shape_hessians(shape_hessians.begin()) - { - // We can enter this function either for the apply() path that has - // n_rows * n_columns entries or for the apply_face() path that only has - // n_rows * 3 entries in the array. Since we cannot decide about the use - // we must allow for both here. - Assert(shape_values.empty() || - shape_values.size() == n_rows * n_columns || - shape_values.size() == 3 * n_rows, - ExcDimensionMismatch(shape_values.size(), n_rows * n_columns)); - Assert(shape_gradients.empty() || - shape_gradients.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns)); - Assert(shape_hessians.empty() || - shape_hessians.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns)); - (void)dummy1; - (void)dummy2; - } - - /** - * Constructor, taking the data from ShapeInfo via raw pointers - */ - EvaluatorTensorProduct(const Number2 * shape_values, - const Number2 * shape_gradients, - const Number2 * shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) - : shape_values(shape_values) - , shape_gradients(shape_gradients) - , shape_hessians(shape_hessians) - { - (void)dummy1; - (void)dummy2; - } - - template - void - values(const Number in[], Number out[]) const - { - apply(shape_values, in, out); - } - - template - void - gradients(const Number in[], Number out[]) const - { - apply(shape_gradients, in, out); - } - - template - void - hessians(const Number in[], Number out[]) const - { - apply(shape_hessians, in, out); - } - - template - void - values_one_line(const Number in[], Number out[]) const - { - Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); - } - - template - void - gradients_one_line(const Number in[], Number out[]) const - { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, in, out); - } - - template - void - hessians_one_line(const Number in[], Number out[]) const - { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, in, out); - } - - /** - * This function applies the tensor product kernel, corresponding to a - * multiplication of 1d stripes, along the given @p direction of the tensor - * data in the input array. This function allows the @p in and @p out - * arrays to alias for the case n_rows == n_columns, i.e., it is safe to - * perform the contraction in place where @p in and @p out point to the - * same address. For the case n_rows != n_columns, the output is in general - * not correct. - * - * @tparam direction Direction that is evaluated - * @tparam contract_over_rows If true, the tensor contraction sums - * over the rows in the given @p shape_data - * array, otherwise it sums over the columns - * @tparam add If true, the result is added to the output vector, else - * the computed values overwrite the content in the output - * @tparam one_line If true, the kernel is only applied along a single 1d - * stripe within a dim-dimensional tensor, not the full - * n_rows^dim points as in the @p false case. - * - * @param shape_data Transformation matrix with @p n_rows rows and - * @p n_columns columns, stored in row-major format - * @param in Pointer to the start of the input data vector - * @param out Pointer to the start of the output data vector - */ - template - static void - apply(const Number2 *DEAL_II_RESTRICT shape_data, - const Number * in, - Number * out); - - /** - * This function applies the tensor product operation to produce face values - * from cell values. As opposed to the apply method, this method assumes - * that the directions orthogonal to the face have n_rows degrees of - * freedom per direction and not n_columns for those directions lower than - * the one currently applied. In other words, apply_face() must be called - * before calling any interpolation within the face. - * - * @tparam face_direction Direction of the normal vector (0=x, 1=y, etc) - * @tparam contract_onto_face If true, the input vector is of size n_rows^dim - * and interpolation into n_rows^(dim-1) points - * is performed. This is a typical scenario in - * FEFaceEvaluation::evaluate() calls. If false, - * data from n_rows^(dim-1) points is expanded - * into the n_rows^dim points of the higher- - * dimensional data array. Derivatives in the - * case contract_onto_face==false are summed - * together - * @tparam add If true, the result is added to the output vector, else - * the computed values overwrite the content in the output - * @tparam max_derivative Sets the number of derivatives that should be - * computed. 0 means only values, 1 means values and first - * derivatives, 2 second derivates. Note that all the - * derivatives access the data in @p shape_values passed to - * the constructor of the class - * - * @param in address of the input data vector - * @param out address of the output data vector - */ - template - void - apply_face(const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT out) const; - - private: - const Number2 *shape_values; - const Number2 *shape_gradients; - const Number2 *shape_hessians; - }; - - - - template - template - inline void - EvaluatorTensorProduct::apply(const Number2 *DEAL_II_RESTRICT - shape_data, - const Number *in, - Number * out) - { - static_assert(one_line == false || direction == dim - 1, - "Single-line evaluation only works for direction=dim-1."); - Assert(shape_data != nullptr, - ExcMessage( - "The given array shape_data must not be the null pointer!")); - Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || - in != out, - ExcMessage("In-place operation only supported for " - "n_rows==n_columns or single-line interpolation")); - AssertIndexRange(direction, dim); - constexpr int mm = contract_over_rows ? n_rows : n_columns, - nn = contract_over_rows ? n_columns : n_rows; - - constexpr int stride = Utilities::pow(n_columns, direction); - constexpr int n_blocks1 = one_line ? 1 : stride; - constexpr int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); - - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - Number x[mm]; - for (int i = 0; i < mm; ++i) - x[i] = in[stride * i]; - for (int col = 0; col < nn; ++col) - { - Number2 val0; - if (contract_over_rows == true) - val0 = shape_data[col]; - else - val0 = shape_data[col * n_columns]; - Number res0 = val0 * x[0]; - for (int i = 1; i < mm; ++i) - { - if (contract_over_rows == true) - val0 = shape_data[i * n_columns + col]; - else - val0 = shape_data[col * n_columns + i]; - res0 += val0 * x[i]; - } - if (add) - out[stride * col] += res0; - else - out[stride * col] = res0; - } - - if (one_line == false) - { - ++in; - ++out; - } - } - if (one_line == false) - { - in += stride * (mm - 1); - out += stride * (nn - 1); - } - } - } - - - - template - template - inline void - EvaluatorTensorProduct::apply_face(const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT - out) const - { - Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported")); - static_assert(max_derivative >= 0 && max_derivative < 3, - "Only derivative orders 0-2 implemented"); - Assert(shape_values != nullptr, - ExcMessage( - "The given array shape_values must not be the null pointer.")); - - constexpr int n_blocks1 = (dim > 1 ? n_rows : 1); - constexpr int n_blocks2 = (dim > 2 ? n_rows : 1); - - AssertIndexRange(face_direction, dim); - constexpr int in_stride = Utilities::pow(n_rows, face_direction); - constexpr int out_stride = Utilities::pow(n_rows, dim - 1); - const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values; - - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - if (contract_onto_face == true) - { - Number res0 = shape_values[0] * in[0]; - Number res1, res2; - if (max_derivative > 0) - res1 = shape_values[n_rows] * in[0]; - if (max_derivative > 1) - res2 = shape_values[2 * n_rows] * in[0]; - for (int ind = 1; ind < n_rows; ++ind) - { - res0 += shape_values[ind] * in[in_stride * ind]; - if (max_derivative > 0) - res1 += shape_values[ind + n_rows] * in[in_stride * ind]; - if (max_derivative > 1) - res2 += - shape_values[ind + 2 * n_rows] * in[in_stride * ind]; - } - if (add) - { - out[0] += res0; - if (max_derivative > 0) - out[out_stride] += res1; - if (max_derivative > 1) - out[2 * out_stride] += res2; - } - else - { - out[0] = res0; - if (max_derivative > 0) - out[out_stride] = res1; - if (max_derivative > 1) - out[2 * out_stride] = res2; - } - } - else - { - for (int col = 0; col < n_rows; ++col) - { - if (add) - out[col * in_stride] += shape_values[col] * in[0]; - else - out[col * in_stride] = shape_values[col] * in[0]; - if (max_derivative > 0) - out[col * in_stride] += - shape_values[col + n_rows] * in[out_stride]; - if (max_derivative > 1) - out[col * in_stride] += - shape_values[col + 2 * n_rows] * in[2 * out_stride]; - } - } - - // increment: in regular case, just go to the next point in - // x-direction. If we are at the end of one chunk in x-dir, need - // to jump over to the next layer in z-direction - switch (face_direction) - { - case 0: - in += contract_onto_face ? n_rows : 1; - out += contract_onto_face ? 1 : n_rows; - break; - case 1: - ++in; - ++out; - // faces 2 and 3 in 3d use local coordinate system zx, which - // is the other way around compared to the tensor - // product. Need to take that into account. - if (dim == 3) - { - if (contract_onto_face) - out += n_rows - 1; - else - in += n_rows - 1; - } - break; - case 2: - ++in; - ++out; - break; - default: - Assert(false, ExcNotImplemented()); - } - } - - // adjust for local coordinate system zx - if (face_direction == 1 && dim == 3) - { - if (contract_onto_face) - { - in += n_rows * (n_rows - 1); - out -= n_rows * n_rows - 1; - } - else - { - out += n_rows * (n_rows - 1); - in -= n_rows * n_rows - 1; - } - } - } - } - - - - /** - * Internal evaluator for shape function using the tensor product form - * of the basis functions. The same as the other templated class but - * without making use of template arguments and variable loop bounds - * instead. - * - * @tparam dim Space dimension in which this class is applied - * @tparam Number Abstract number type for input and output arrays - * @tparam Number2 Abstract number type for coefficient arrays (defaults to - * same type as the input/output arrays); must implement - * operator* with Number and produce Number as an output to - * be a valid type - */ - template - struct EvaluatorTensorProduct - { - static constexpr unsigned int n_rows_of_product = - numbers::invalid_unsigned_int; - static constexpr unsigned int n_columns_of_product = - numbers::invalid_unsigned_int; - - /** - * Empty constructor. Does nothing. Be careful when using 'values' and - * related methods because they need to be filled with the other constructor - */ - EvaluatorTensorProduct() - : shape_values(nullptr) - , shape_gradients(nullptr) - , shape_hessians(nullptr) - , n_rows(numbers::invalid_unsigned_int) - , n_columns(numbers::invalid_unsigned_int) - {} - - /** - * Constructor, taking the data from ShapeInfo - */ - EvaluatorTensorProduct(const AlignedVector &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int n_rows, - const unsigned int n_columns) - : shape_values(shape_values.begin()) - , shape_gradients(shape_gradients.begin()) - , shape_hessians(shape_hessians.begin()) - , n_rows(n_rows) - , n_columns(n_columns) - { - // We can enter this function either for the apply() path that has - // n_rows * n_columns entries or for the apply_face() path that only has - // n_rows * 3 entries in the array. Since we cannot decide about the use - // we must allow for both here. - Assert(shape_values.empty() || - shape_values.size() == n_rows * n_columns || - shape_values.size() == n_rows * 3, - ExcDimensionMismatch(shape_values.size(), n_rows * n_columns)); - Assert(shape_gradients.empty() || - shape_gradients.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns)); - Assert(shape_hessians.empty() || - shape_hessians.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns)); - } - - /** - * Constructor, taking the data from ShapeInfo - */ - EvaluatorTensorProduct(const Number2 * shape_values, - const Number2 * shape_gradients, - const Number2 * shape_hessians, - const unsigned int n_rows, - const unsigned int n_columns) - : shape_values(shape_values) - , shape_gradients(shape_gradients) - , shape_hessians(shape_hessians) - , n_rows(n_rows) - , n_columns(n_columns) - {} - - template - void - values(const Number *in, Number *out) const - { - apply(shape_values, in, out); - } - - template - void - gradients(const Number *in, Number *out) const - { - apply(shape_gradients, in, out); - } - - template - void - hessians(const Number *in, Number *out) const - { - apply(shape_hessians, in, out); - } - - template - void - values_one_line(const Number in[], Number out[]) const - { - Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); - } - - template - void - gradients_one_line(const Number in[], Number out[]) const - { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, in, out); - } - - template - void - hessians_one_line(const Number in[], Number out[]) const - { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, in, out); - } - - template - void - apply(const Number2 *DEAL_II_RESTRICT shape_data, - const Number * in, - Number * out) const; - - template - void - apply_face(const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT out) const; - - const Number2 * shape_values; - const Number2 * shape_gradients; - const Number2 * shape_hessians; - const unsigned int n_rows; - const unsigned int n_columns; - }; - - - - template - template - inline void - EvaluatorTensorProduct::apply( - const Number2 *DEAL_II_RESTRICT shape_data, - const Number * in, - Number * out) const - { - static_assert(one_line == false || direction == dim - 1, - "Single-line evaluation only works for direction=dim-1."); - Assert(shape_data != nullptr, - ExcMessage( - "The given array shape_data must not be the null pointer!")); - Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || - in != out, - ExcMessage("In-place operation only supported for " - "n_rows==n_columns or single-line interpolation")); - AssertIndexRange(direction, dim); - const int mm = contract_over_rows ? n_rows : n_columns, - nn = contract_over_rows ? n_columns : n_rows; - - const int stride = - direction == 0 ? 1 : Utilities::fixed_power(n_columns); - const int n_blocks1 = one_line ? 1 : stride; - const int n_blocks2 = direction >= dim - 1 ? - 1 : - Utilities::fixed_power(n_rows); - Assert(n_rows <= 128, ExcNotImplemented()); - - // specialization for n_rows = 2 that manually unrolls the innermost loop - // to make the operation perform better (not completely as good as the - // templated one, but much better than the generic version down below, - // because the loop over col can be more effectively unrolled by the - // compiler) - if (contract_over_rows && n_rows == 2) - { - const Number2 *shape_data_1 = shape_data + n_columns; - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - const Number x0 = in[0], x1 = in[stride]; - for (int col = 0; col < nn; ++col) - { - const Number result = - shape_data[col] * x0 + shape_data_1[col] * x1; - if (add) - out[stride * col] += result; - else - out[stride * col] = result; - } - - if (one_line == false) - { - ++in; - ++out; - } - } - if (one_line == false) - { - in += stride * (mm - 1); - out += stride * (nn - 1); - } - } - } - // specialization for n = 3 - else if (contract_over_rows && n_rows == 3) - { - const Number2 *shape_data_1 = shape_data + n_columns; - const Number2 *shape_data_2 = shape_data + 2 * n_columns; - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride]; - for (int col = 0; col < nn; ++col) - { - const Number result = shape_data[col] * x0 + - shape_data_1[col] * x1 + - shape_data_2[col] * x2; - if (add) - out[stride * col] += result; - else - out[stride * col] = result; - } - - if (one_line == false) - { - ++in; - ++out; - } - } - if (one_line == false) - { - in += stride * (mm - 1); - out += stride * (nn - 1); - } - } - } - // general loop for all other cases - else - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - Number x[129]; - for (int i = 0; i < mm; ++i) - x[i] = in[stride * i]; - for (int col = 0; col < nn; ++col) - { - Number2 val0; - if (contract_over_rows == true) - val0 = shape_data[col]; - else - val0 = shape_data[col * n_columns]; - Number res0 = val0 * x[0]; - for (int i = 1; i < mm; ++i) - { - if (contract_over_rows == true) - val0 = shape_data[i * n_columns + col]; - else - val0 = shape_data[col * n_columns + i]; - res0 += val0 * x[i]; - } - if (add) - out[stride * col] += res0; - else - out[stride * col] = res0; - } - - if (one_line == false) - { - ++in; - ++out; - } - } - if (one_line == false) - { - in += stride * (mm - 1); - out += stride * (nn - 1); - } - } - } - - - - template - template - inline void - EvaluatorTensorProduct:: - apply_face(const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT out) const - { - Assert(shape_values != nullptr, - ExcMessage( - "The given array shape_data must not be the null pointer!")); - static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported"); - const int n_blocks1 = dim > 1 ? n_rows : 1; - const int n_blocks2 = dim > 2 ? n_rows : 1; - - AssertIndexRange(face_direction, dim); - const int in_stride = - face_direction > 0 ? Utilities::fixed_power(n_rows) : 1; - const int out_stride = - dim > 1 ? Utilities::fixed_power(n_rows) : 1; - - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - if (contract_onto_face == true) - { - Number res0 = shape_values[0] * in[0]; - Number res1, res2; - if (max_derivative > 0) - res1 = shape_values[n_rows] * in[0]; - if (max_derivative > 1) - res2 = shape_values[2 * n_rows] * in[0]; - for (unsigned int ind = 1; ind < n_rows; ++ind) - { - res0 += shape_values[ind] * in[in_stride * ind]; - if (max_derivative > 0) - res1 += shape_values[ind + n_rows] * in[in_stride * ind]; - if (max_derivative > 1) - res2 += - shape_values[ind + 2 * n_rows] * in[in_stride * ind]; - } - if (add) - { - out[0] += res0; - if (max_derivative > 0) - out[out_stride] += res1; - if (max_derivative > 1) - out[2 * out_stride] += res2; - } - else - { - out[0] = res0; - if (max_derivative > 0) - out[out_stride] = res1; - if (max_derivative > 1) - out[2 * out_stride] = res2; - } - } - else - { - for (unsigned int col = 0; col < n_rows; ++col) - { - if (add) - out[col * in_stride] += shape_values[col] * in[0]; - else - out[col * in_stride] = shape_values[col] * in[0]; - if (max_derivative > 0) - out[col * in_stride] += - shape_values[col + n_rows] * in[out_stride]; - if (max_derivative > 1) - out[col * in_stride] += - shape_values[col + 2 * n_rows] * in[2 * out_stride]; - } - } - - // increment: in regular case, just go to the next point in - // x-direction. If we are at the end of one chunk in x-dir, need - // to jump over to the next layer in z-direction - switch (face_direction) - { - case 0: - in += contract_onto_face ? n_rows : 1; - out += contract_onto_face ? 1 : n_rows; - break; - case 1: - ++in; - ++out; - // faces 2 and 3 in 3d use local coordinate system zx, which - // is the other way around compared to the tensor - // product. Need to take that into account. - if (dim == 3) - { - if (contract_onto_face) - out += n_rows - 1; - else - in += n_rows - 1; - } - break; - case 2: - ++in; - ++out; - break; - default: - Assert(false, ExcNotImplemented()); - } + * One-dimensional kernel for use by the generic tensor product + * interpolation as provided by the class EvaluatorTensorProduct, + * implementing a matrix-vector product along this dimension, controlled by + * the number of rows and columns and the stride in the input and output + * arrays, which are embedded into some lexicographic ordering of unknowns + * in a tensor-product arrangement. + * + * Besides this generic function for templated loop lengths, there are + * several specializations of this class to account for run-time matrix + * sizes as well as some symmetries that reduce the data access or + * arithmetic operations. The specializations are technically realized by + * conditional function overloading with std::enable_if_t based on the first + * template parameter. + */ + template + std::enable_if_t<(variant == evaluate_general), void> + apply_matrix_vector_product(const Number2 *matrix, + const Number * in, + Number * out) + { + // We can only statically assert that one argument is non-zero because + // face evaluation might instantiate some functions, so we need to use the + // run-time assert to verify that we do not end up involuntarily. + static_assert(n_rows > 0 || n_columns > 0, + "Specialization only for n_rows, n_columns > 0"); + Assert(n_rows > 0 && n_columns > 0, + ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " + + std::to_string(n_rows) + ", " + + std::to_string(n_columns) + " was passed!")); + static_assert(quantity == EvaluatorQuantity::value, + "This function should only use EvaluatorQuantity::value"); + + constexpr int mm = transpose_matrix ? n_rows : n_columns, + nn = transpose_matrix ? n_columns : n_rows; + + std::array x; + for (int i = 0; i < mm; ++i) + x[i] = in[stride_in * i]; + for (int col = 0; col < nn; ++col) + { + Number res0; + if (transpose_matrix == true) + { + res0 = matrix[col] * x[0]; + for (int i = 1; i < mm; ++i) + res0 += matrix[i * n_columns + col] * x[i]; } - if (face_direction == 1 && dim == 3) + else { - // adjust for local coordinate system zx - if (contract_onto_face) - { - in += n_rows * (n_rows - 1); - out -= n_rows * n_rows - 1; - } - else - { - out += n_rows * (n_rows - 1); - in -= n_rows * n_rows - 1; - } + res0 = matrix[col * n_columns] * x[0]; + for (int i = 1; i < mm; ++i) + res0 += matrix[col * n_columns + i] * x[i]; } + if (add) + out[stride_out * col] += res0; + else + out[stride_out * col] = res0; } } /** - * Internal evaluator for 1d-3d shape function using the tensor product form - * of the basis functions. This class specializes the general application of - * tensor-product based elements for "symmetric" finite elements, i.e., when - * the shape functions are symmetric about 0.5 and the quadrature points - * are, too. - * - * @tparam dim Space dimension in which this class is applied - * @tparam n_rows Number of rows in the transformation matrix, which corresponds - * to the number of 1d shape functions in the usual tensor - * contraction setting - * @tparam n_columns Number of columns in the transformation matrix, which - * corresponds to the number of 1d shape functions in the - * usual tensor contraction setting - * @tparam Number Abstract number type for input and output arrays - * @tparam Number2 Abstract number type for coefficient arrays (defaults to - * same type as the input/output arrays); must implement - * operator* with Number and produce Number as an output to - * be a valid type + * Specialization of the matrix-vector kernel for run-time loop bounds in + * the generic evaluator. */ - template - struct EvaluatorTensorProduct + std::enable_if_t<(variant == evaluate_general), void> + apply_matrix_vector_product(const Number2 *matrix, + const Number * in, + Number * out, + const int n_rows, + const int n_columns, + const int stride_in, + const int stride_out) { - static constexpr unsigned int n_rows_of_product = - Utilities::pow(n_rows, dim); - static constexpr unsigned int n_columns_of_product = - Utilities::pow(n_columns, dim); - - /** - * Constructor, taking the data from ShapeInfo - */ - EvaluatorTensorProduct(const AlignedVector &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) - : shape_values(shape_values.begin()) - , shape_gradients(shape_gradients.begin()) - , shape_hessians(shape_hessians.begin()) - { - Assert(shape_values.empty() || shape_values.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_values.size(), n_rows * n_columns)); - Assert(shape_gradients.empty() || - shape_gradients.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns)); - Assert(shape_hessians.empty() || - shape_hessians.size() == n_rows * n_columns, - ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns)); - (void)dummy1; - (void)dummy2; - } - - template - void - values(const Number in[], Number out[]) const; + const int mm = transpose_matrix ? n_rows : n_columns, + nn = transpose_matrix ? n_columns : n_rows; + Assert(n_rows <= 128, ExcNotImplemented()); + Assert(n_rows > 0 && n_columns > 0, + ExcInternalError("Empty evaluation task!")); + Assert(n_rows > 0 && n_columns > 0, + ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " + + std::to_string(n_rows) + ", " + + std::to_string(n_columns) + " was passed!")); - template - void - gradients(const Number in[], Number out[]) const; + static_assert(quantity == EvaluatorQuantity::value, + "This function should only use EvaluatorQuantity::value"); - template - void - hessians(const Number in[], Number out[]) const; + // specialization for n_rows = 2 that manually unrolls the innermost loop + // to make the operation perform better (not completely as good as the + // templated one, but much better than the generic version down below, + // because the loop over col can be more effectively unrolled by the + // compiler) + if (transpose_matrix && n_rows == 2) + { + const Number2 *matrix_1 = matrix + n_columns; + const Number x0 = in[0], x1 = in[stride_in]; + for (int col = 0; col < nn; ++col) + { + const Number result = matrix[col] * x0 + matrix_1[col] * x1; + if (add) + out[stride_out * col] += result; + else + out[stride_out * col] = result; + } + } + else if (transpose_matrix && n_rows == 3) + { + const Number2 *matrix_1 = matrix + n_columns; + const Number2 *matrix_2 = matrix_1 + n_columns; + const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in]; + for (int col = 0; col < nn; ++col) + { + const Number result = + matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2; + if (add) + out[stride_out * col] += result; + else + out[stride_out * col] = result; + } + } + else + { + std::array x; + for (int i = 0; i < mm; ++i) + x[i] = in[stride_in * i]; - private: - const Number2 *shape_values; - const Number2 *shape_gradients; - const Number2 *shape_hessians; - }; + Number res0; + for (int col = 0; col < nn; ++col) + { + if (transpose_matrix == true) + { + res0 = matrix[col] * x[0]; + for (int i = 1; i < mm; ++i) + res0 += matrix[i * n_columns + col] * x[i]; + } + else + { + res0 = matrix[col * n_columns] * x[0]; + for (int i = 1; i < mm; ++i) + res0 += matrix[col * n_columns + i] * x[i]; + } + if (add) + out[stride_out * col] += res0; + else + out[stride_out * col] = res0; + } + } + } - // In this case, the 1d shape values read (sorted lexicographically, rows - // run over 1d dofs, columns over quadrature points): - // Q2 --> [ 0.687 0 -0.087 ] - // [ 0.4 1 0.4 ] - // [-0.087 0 0.687 ] - // Q3 --> [ 0.66 0.003 0.002 0.049 ] - // [ 0.521 1.005 -0.01 -0.230 ] - // [-0.230 -0.01 1.005 0.521 ] - // [ 0.049 0.002 0.003 0.66 ] - // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ] - // [ 0.608 1.059 0 0.039 0.176 ] - // [-0.409 -0.113 1 -0.113 -0.409 ] - // [ 0.176 0.039 0 1.059 0.608 ] - // [-0.032 -0.007 0 0.022 0.658 ] - // - // In these matrices, we want to use avoid computations involving zeros and - // ones and in addition use the symmetry in entries to reduce the number of - // read operations. - template - template - inline void - EvaluatorTensorProduct::values(const Number in[], Number out[]) const + std::enable_if_t<(variant == evaluate_symmetric), void> + apply_matrix_vector_product(const Number2 *matrix, + const Number * in, + Number * out) { - Assert(shape_values != nullptr, ExcNotInitialized()); - AssertIndexRange(direction, dim); - constexpr int mm = contract_over_rows ? n_rows : n_columns, - nn = contract_over_rows ? n_columns : n_rows; + // We can only statically assert that one argument is non-zero because + // face evaluation might instantiate some functions, so we need to use the + // run-time assert to verify that we do not end up involuntarily. + static_assert(n_rows > 0 || n_columns > 0, + "Specialization only for n_rows, n_columns > 0"); + Assert(n_rows > 0 && n_columns > 0, + ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " + + std::to_string(n_rows) + ", " + + std::to_string(n_columns) + " was passed!")); + + constexpr int mm = transpose_matrix ? n_rows : n_columns, + nn = transpose_matrix ? n_columns : n_rows; constexpr int n_cols = nn / 2; constexpr int mid = mm / 2; - constexpr int stride = Utilities::pow(n_columns, direction); - constexpr int n_blocks1 = stride; - constexpr int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + std::array x; + for (int i = 0; i < mm; ++i) + x[i] = in[stride_in * i]; - for (int i2 = 0; i2 < n_blocks2; ++i2) + if (quantity == EvaluatorQuantity::value) { - for (int i1 = 0; i1 < n_blocks1; ++i1) + // In this case, the 1d shape values read (sorted lexicographically, + // rows run over 1d dofs, columns over quadrature points): + // Q2 --> [ 0.687 0 -0.087 ] + // [ 0.4 1 0.4 ] + // [-0.087 0 0.687 ] + // Q3 --> [ 0.66 0.003 0.002 0.049 ] + // [ 0.521 1.005 -0.01 -0.230 ] + // [-0.230 -0.01 1.005 0.521 ] + // [ 0.049 0.002 0.003 0.66 ] + // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ] + // [ 0.608 1.059 0 0.039 0.176 ] + // [-0.409 -0.113 1 -0.113 -0.409 ] + // [ 0.176 0.039 0 1.059 0.608 ] + // [-0.032 -0.007 0 0.022 0.658 ] + // + // In these matrices, we want to use avoid computations involving + // zeros and ones and use the symmetry in entries starting from (1,1) + // forward and (N,N) backward, respectively to reduce the number of + // read operations. + for (int col = 0; col < n_cols; ++col) { - for (int col = 0; col < n_cols; ++col) + Number2 val0, val1; + Number res0, res1; + if (transpose_matrix == true) { - Number2 val0, val1; - Number in0, in1, res0, res1; - if (contract_over_rows == true) - { - val0 = shape_values[col]; - val1 = shape_values[nn - 1 - col]; - } - else - { - val0 = shape_values[col * n_columns]; - val1 = shape_values[(col + 1) * n_columns - 1]; - } - if (mid > 0) - { - in0 = in[0]; - in1 = in[stride * (mm - 1)]; - res0 = val0 * in0; - res1 = val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - for (int ind = 1; ind < mid; ++ind) - { - if (contract_over_rows == true) - { - val0 = shape_values[ind * n_columns + col]; - val1 = shape_values[ind * n_columns + nn - 1 - col]; - } - else - { - val0 = shape_values[col * n_columns + ind]; - val1 = - shape_values[(col + 1) * n_columns - 1 - ind]; - } - in0 = in[stride * ind]; - in1 = in[stride * (mm - 1 - ind)]; - res0 += val0 * in0; - res1 += val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - } - } - else - res0 = res1 = Number(); - if (contract_over_rows == true) + val0 = matrix[col]; + val1 = matrix[nn - 1 - col]; + } + else + { + val0 = matrix[col * n_columns]; + val1 = matrix[(col + 1) * n_columns - 1]; + } + if (mid > 0) + { + res0 = val0 * x[0]; + res1 = val1 * x[0]; + res0 += val1 * x[mm - 1]; + res1 += val0 * x[mm - 1]; + for (int ind = 1; ind < mid; ++ind) { - if (mm % 2 == 1) + if (transpose_matrix == true) { - val0 = shape_values[mid * n_columns + col]; - in1 = val0 * in[stride * mid]; - res0 += in1; - res1 += in1; + val0 = matrix[ind * n_columns + col]; + val1 = matrix[ind * n_columns + nn - 1 - col]; } - } - else - { - if (mm % 2 == 1 && nn % 2 == 0) + else { - val0 = shape_values[col * n_columns + mid]; - in1 = val0 * in[stride * mid]; - res0 += in1; - res1 += in1; + val0 = matrix[col * n_columns + ind]; + val1 = matrix[(col + 1) * n_columns - 1 - ind]; } + res0 += val0 * x[ind]; + res1 += val1 * x[ind]; + res0 += val1 * x[mm - 1 - ind]; + res1 += val0 * x[mm - 1 - ind]; } - if (add) + } + else + res0 = res1 = Number(); + if (transpose_matrix == true) + { + if (mm % 2 == 1) { - out[stride * col] += res0; - out[stride * (nn - 1 - col)] += res1; + const Number tmp = matrix[mid * n_columns + col] * x[mid]; + res0 += tmp; + res1 += tmp; } - else + } + else + { + if (mm % 2 == 1 && nn % 2 == 0) { - out[stride * col] = res0; - out[stride * (nn - 1 - col)] = res1; + const Number tmp = matrix[col * n_columns + mid] * x[mid]; + res0 += tmp; + res1 += tmp; } } - if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1) + if (add) { - if (add) - out[stride * n_cols] += in[stride * mid]; - else - out[stride * n_cols] = in[stride * mid]; + out[stride_out * col] += res0; + out[stride_out * (nn - 1 - col)] += res1; + } + else + { + out[stride_out * col] = res0; + out[stride_out * (nn - 1 - col)] = res1; } - else if (contract_over_rows == true && nn % 2 == 1) + } + if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1) + { + if (add) + out[stride_out * n_cols] += x[mid]; + else + out[stride_out * n_cols] = x[mid]; + } + else if (transpose_matrix == true && nn % 2 == 1) + { + Number res0; + if (mid > 0) { - Number res0; - Number2 val0 = shape_values[n_cols]; - if (mid > 0) + res0 = matrix[n_cols] * (x[0] + x[mm - 1]); + for (int ind = 1; ind < mid; ++ind) { - res0 = val0 * (in[0] + in[stride * (mm - 1)]); - for (int ind = 1; ind < mid; ++ind) - { - val0 = shape_values[ind * n_columns + n_cols]; - res0 += val0 * (in[stride * ind] + - in[stride * (mm - 1 - ind)]); - } + const Number2 val0 = matrix[ind * n_columns + n_cols]; + res0 += val0 * (x[ind] + in[mm - 1 - ind]); } - else - res0 = Number(); - if (add) - out[stride * n_cols] += res0; - else - out[stride * n_cols] = res0; } - else if (contract_over_rows == false && nn % 2 == 1) + else + res0 = Number(); + if (add) + out[stride_out * n_cols] += res0; + else + out[stride_out * n_cols] = res0; + } + else if (transpose_matrix == false && nn % 2 == 1) + { + Number res0; + if (mid > 0) + { + res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]); + for (int ind = 1; ind < mid; ++ind) + { + const Number2 val0 = matrix[n_cols * n_columns + ind]; + res0 += val0 * (x[ind] + x[mm - 1 - ind]); + ; + } + if (mm % 2) + res0 += x[mid]; + } + else + res0 = in[0]; + if (add) + out[stride_out * n_cols] += res0; + else + out[stride_out * n_cols] = res0; + } + } + else if (quantity == EvaluatorQuantity::gradient) + { + // For the specialized loop used for gradient computations we again + // exploit symmetries according to the following entries (sorted + // lexicographically, rows run over 1d dofs, columns over quadrature + // points): + // Q2 --> [-2.549 -1 0.549 ] + // [ 3.098 0 -3.098 ] + // [-0.549 1 2.549 ] + // Q3 --> [-4.315 -1.03 0.5 -0.44 ] + // [ 6.07 -1.44 -2.97 2.196 ] + // [-2.196 2.97 1.44 -6.07 ] + // [ 0.44 -0.5 1.03 4.315 ] + // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ] + // [10.111 -2.76 -2.667 2.066 -2.306 ] + // [-5.688 5.773 0 -5.773 5.688 ] + // [ 2.306 -2.066 2.667 2.76 -10.111 ] + // [-0.413 0.353 -0.333 -0.353 0.413 ] + for (int col = 0; col < n_cols; ++col) + { + Number2 val0, val1; + Number res0, res1; + if (transpose_matrix == true) + { + val0 = matrix[col]; + val1 = matrix[nn - 1 - col]; + } + else + { + val0 = matrix[col * n_columns]; + val1 = matrix[(nn - col - 1) * n_columns]; + } + if (mid > 0) { - Number res0; - if (mid > 0) + res0 = val0 * x[0]; + res1 = val1 * x[0]; + res0 -= val1 * x[mm - 1]; + res1 -= val0 * x[mm - 1]; + for (int ind = 1; ind < mid; ++ind) { - Number2 val0 = shape_values[n_cols * n_columns]; - res0 = val0 * (in[0] + in[stride * (mm - 1)]); - for (int ind = 1; ind < mid; ++ind) + if (transpose_matrix == true) + { + val0 = matrix[ind * n_columns + col]; + val1 = matrix[ind * n_columns + nn - 1 - col]; + } + else { - val0 = shape_values[n_cols * n_columns + ind]; - Number in1 = val0 * (in[stride * ind] + - in[stride * (mm - 1 - ind)]); - res0 += in1; + val0 = matrix[col * n_columns + ind]; + val1 = matrix[(nn - col - 1) * n_columns + ind]; } - if (mm % 2) - res0 += in[stride * mid]; + res0 += val0 * x[ind]; + res1 += val1 * x[ind]; + res0 -= val1 * x[mm - 1 - ind]; + res1 -= val0 * x[mm - 1 - ind]; } + } + else + res0 = res1 = Number(); + if (mm % 2 == 1) + { + if (transpose_matrix == true) + val0 = matrix[mid * n_columns + col]; else - res0 = in[0]; - if (add) - out[stride * n_cols] += res0; + val0 = matrix[col * n_columns + mid]; + const Number tmp = val0 * x[mid]; + res0 += tmp; + res1 -= tmp; + } + if (add) + { + out[stride_out * col] += res0; + out[stride_out * (nn - 1 - col)] += res1; + } + else + { + out[stride_out * col] = res0; + out[stride_out * (nn - 1 - col)] = res1; + } + } + if (nn % 2 == 1) + { + Number2 val0; + Number res0; + if (transpose_matrix == true) + val0 = matrix[n_cols]; + else + val0 = matrix[n_cols * n_columns]; + res0 = val0 * (x[0] - x[mm - 1]); + for (int ind = 1; ind < mid; ++ind) + { + if (transpose_matrix == true) + val0 = matrix[ind * n_columns + n_cols]; else - out[stride * n_cols] = res0; + val0 = matrix[n_cols * n_columns + ind]; + Number in1 = val0 * (x[ind] - x[mm - 1 - ind]); + res0 += in1; } - - ++in; - ++out; + if (add) + out[stride_out * n_cols] += res0; + else + out[stride_out * n_cols] = res0; } - in += stride * (mm - 1); - out += stride * (nn - 1); } - } - - - - // For the specialized loop used for the gradient computation in - // here, the 1d shape values read (sorted lexicographically, rows - // run over 1d dofs, columns over quadrature points): - // Q2 --> [-2.549 -1 0.549 ] - // [ 3.098 0 -3.098 ] - // [-0.549 1 2.549 ] - // Q3 --> [-4.315 -1.03 0.5 -0.44 ] - // [ 6.07 -1.44 -2.97 2.196 ] - // [-2.196 2.97 1.44 -6.07 ] - // [ 0.44 -0.5 1.03 4.315 ] - // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ] - // [10.111 -2.76 -2.667 2.066 -2.306 ] - // [-5.688 5.773 0 -5.773 5.688 ] - // [ 2.306 -2.066 2.667 2.76 -10.111 ] - // [-0.413 0.353 -0.333 -0.353 0.413 ] - // - // In these matrices, we want to use avoid computations involving - // zeros and ones and in addition use the symmetry in entries to - // reduce the number of read operations. - template - template - inline void - EvaluatorTensorProduct::gradients(const Number in[], - Number out[]) const - { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - AssertIndexRange(direction, dim); - constexpr int mm = contract_over_rows ? n_rows : n_columns, - nn = contract_over_rows ? n_columns : n_rows; - constexpr int n_cols = nn / 2; - constexpr int mid = mm / 2; - - constexpr int stride = Utilities::pow(n_columns, direction); - constexpr int n_blocks1 = stride; - constexpr int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); - - for (int i2 = 0; i2 < n_blocks2; ++i2) + else { - for (int i1 = 0; i1 < n_blocks1; ++i1) + // Hessians are almost the same as values, apart from some missing '1' + // entries + for (int col = 0; col < n_cols; ++col) { - for (int col = 0; col < n_cols; ++col) + Number2 val0, val1; + Number res0, res1; + if (transpose_matrix == true) { - Number2 val0, val1; - Number in0, in1, res0, res1; - if (contract_over_rows == true) - { - val0 = shape_gradients[col]; - val1 = shape_gradients[nn - 1 - col]; - } - else - { - val0 = shape_gradients[col * n_columns]; - val1 = shape_gradients[(nn - col - 1) * n_columns]; - } - if (mid > 0) + val0 = matrix[col]; + val1 = matrix[nn - 1 - col]; + } + else + { + val0 = matrix[col * n_columns]; + val1 = matrix[(col + 1) * n_columns - 1]; + } + if (mid > 0) + { + res0 = val0 * x[0]; + res1 = val1 * x[0]; + res0 += val1 * x[mm - 1]; + res1 += val0 * x[mm - 1]; + for (int ind = 1; ind < mid; ++ind) { - in0 = in[0]; - in1 = in[stride * (mm - 1)]; - res0 = val0 * in0; - res1 = val1 * in0; - res0 -= val1 * in1; - res1 -= val0 * in1; - for (int ind = 1; ind < mid; ++ind) + if (transpose_matrix == true) { - if (contract_over_rows == true) - { - val0 = shape_gradients[ind * n_columns + col]; - val1 = - shape_gradients[ind * n_columns + nn - 1 - col]; - } - else - { - val0 = shape_gradients[col * n_columns + ind]; - val1 = - shape_gradients[(nn - col - 1) * n_columns + ind]; - } - in0 = in[stride * ind]; - in1 = in[stride * (mm - 1 - ind)]; - res0 += val0 * in0; - res1 += val1 * in0; - res0 -= val1 * in1; - res1 -= val0 * in1; + val0 = matrix[ind * n_columns + col]; + val1 = matrix[ind * n_columns + nn - 1 - col]; } - } - else - res0 = res1 = Number(); - if (mm % 2 == 1) - { - if (contract_over_rows == true) - val0 = shape_gradients[mid * n_columns + col]; else - val0 = shape_gradients[col * n_columns + mid]; - in1 = val0 * in[stride * mid]; - res0 += in1; - res1 -= in1; - } - if (add) - { - out[stride * col] += res0; - out[stride * (nn - 1 - col)] += res1; - } - else - { - out[stride * col] = res0; - out[stride * (nn - 1 - col)] = res1; + { + val0 = matrix[col * n_columns + ind]; + val1 = matrix[(col + 1) * n_columns - 1 - ind]; + } + res0 += val0 * x[ind]; + res1 += val1 * x[ind]; + res0 += val1 * x[mm - 1 - ind]; + res1 += val0 * x[mm - 1 - ind]; } } - if (nn % 2 == 1) + else + res0 = res1 = Number(); + if (mm % 2 == 1) { - Number2 val0; - Number res0; - if (contract_over_rows == true) - val0 = shape_gradients[n_cols]; + if (transpose_matrix == true) + val0 = matrix[mid * n_columns + col]; else - val0 = shape_gradients[n_cols * n_columns]; - res0 = val0 * (in[0] - in[stride * (mm - 1)]); + val0 = matrix[col * n_columns + mid]; + const Number tmp = val0 * x[mid]; + res0 += tmp; + res1 += tmp; + } + if (add) + { + out[stride_out * col] += res0; + out[stride_out * (nn - 1 - col)] += res1; + } + else + { + out[stride_out * col] = res0; + out[stride_out * (nn - 1 - col)] = res1; + } + } + if (nn % 2 == 1) + { + Number2 val0; + Number res0; + if (transpose_matrix == true) + val0 = matrix[n_cols]; + else + val0 = matrix[n_cols * n_columns]; + if (mid > 0) + { + res0 = val0 * (x[0] + x[mm - 1]); for (int ind = 1; ind < mid; ++ind) { - if (contract_over_rows == true) - val0 = shape_gradients[ind * n_columns + n_cols]; + if (transpose_matrix == true) + val0 = matrix[ind * n_columns + n_cols]; else - val0 = shape_gradients[n_cols * n_columns + ind]; - Number in1 = - val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]); + val0 = matrix[n_cols * n_columns + ind]; + Number in1 = val0 * (x[ind] + x[mm - 1 - ind]); res0 += in1; } - if (add) - out[stride * n_cols] += res0; + } + else + res0 = Number(); + if (mm % 2 == 1) + { + if (transpose_matrix == true) + val0 = matrix[mid * n_columns + n_cols]; else - out[stride * n_cols] = res0; + val0 = matrix[n_cols * n_columns + mid]; + res0 += val0 * x[mid]; } - - ++in; - ++out; + if (add) + out[stride_out * n_cols] += res0; + else + out[stride_out * n_cols] = res0; } - in += stride * (mm - 1); - out += stride * (nn - 1); } } - // evaluates the given shape data in 1d-3d using the tensor product - // form assuming the symmetries of unit cell shape hessians for - // finite elements in FEEvaluation - template - template - inline void - EvaluatorTensorProduct::hessians(const Number in[], - Number out[]) const + std::enable_if_t<(variant == evaluate_evenodd), void> + apply_matrix_vector_product(const Number2 *matrix, + const Number * in, + Number * out, + int n_rows_runtime = 0, + int n_columns_runtime = 0, + int stride_in_runtime = 0, + int stride_out_runtime = 0) { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - AssertIndexRange(direction, dim); - constexpr int mm = contract_over_rows ? n_rows : n_columns; - constexpr int nn = contract_over_rows ? n_columns : n_rows; - constexpr int n_cols = nn / 2; - constexpr int mid = mm / 2; + const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static; + const int n_columns = + n_rows_static == 0 ? n_columns_runtime : n_columns_static; + const int stride_in = + n_rows_static == 0 ? stride_in_runtime : stride_in_static; + const int stride_out = + n_rows_static == 0 ? stride_out_runtime : stride_out_static; + + Assert(n_rows > 0 && n_columns > 0, + ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " + + std::to_string(n_rows) + ", " + + std::to_string(n_columns) + " was passed!")); + + const int mm = transpose_matrix ? n_rows : n_columns, + nn = transpose_matrix ? n_columns : n_rows; + const int n_cols = nn / 2; + const int mid = mm / 2; + + constexpr int max_mid = 16; // for non-templated execution + constexpr int static_mid = + n_rows_static == 0 ? 1 : (transpose_matrix ? n_rows : n_columns) / 2; + const int offset = (n_columns + 1) / 2; - constexpr int stride = Utilities::pow(n_columns, direction); - constexpr int n_blocks1 = stride; - constexpr int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + Assert((n_rows_static != 0 && n_columns_static != 0) || mid <= max_mid, + ExcNotImplemented()); - for (int i2 = 0; i2 < n_blocks2; ++i2) + std::array xp, xm; + for (int i = 0; i < mid; ++i) { - for (int i1 = 0; i1 < n_blocks1; ++i1) + if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient) + { + xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)]; + xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)]; + } + else + { + xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)]; + xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)]; + } + } + Number xmid = in[stride_in * mid]; + for (int col = 0; col < n_cols; ++col) + { + Number r0, r1; + if (mid > 0) { - for (int col = 0; col < n_cols; ++col) + if (transpose_matrix == true) { - Number2 val0, val1; - Number in0, in1, res0, res1; - if (contract_over_rows == true) - { - val0 = shape_hessians[col]; - val1 = shape_hessians[nn - 1 - col]; - } - else - { - val0 = shape_hessians[col * n_columns]; - val1 = shape_hessians[(col + 1) * n_columns - 1]; - } - if (mid > 0) - { - in0 = in[0]; - in1 = in[stride * (mm - 1)]; - res0 = val0 * in0; - res1 = val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - for (int ind = 1; ind < mid; ++ind) - { - if (contract_over_rows == true) - { - val0 = shape_hessians[ind * n_columns + col]; - val1 = - shape_hessians[ind * n_columns + nn - 1 - col]; - } - else - { - val0 = shape_hessians[col * n_columns + ind]; - val1 = - shape_hessians[(col + 1) * n_columns - 1 - ind]; - } - in0 = in[stride * ind]; - in1 = in[stride * (mm - 1 - ind)]; - res0 += val0 * in0; - res1 += val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - } - } - else - res0 = res1 = Number(); - if (mm % 2 == 1) - { - if (contract_over_rows == true) - val0 = shape_hessians[mid * n_columns + col]; - else - val0 = shape_hessians[col * n_columns + mid]; - in1 = val0 * in[stride * mid]; - res0 += in1; - res1 += in1; - } - if (add) - { - out[stride * col] += res0; - out[stride * (nn - 1 - col)] += res1; - } - else - { - out[stride * col] = res0; - out[stride * (nn - 1 - col)] = res1; - } + r0 = matrix[col] * xp[0]; + r1 = matrix[(n_rows - 1) * offset + col] * xm[0]; } - if (nn % 2 == 1) + else { - Number2 val0; - Number res0; - if (contract_over_rows == true) - val0 = shape_hessians[n_cols]; - else - val0 = shape_hessians[n_cols * n_columns]; - if (mid > 0) + r0 = matrix[col * offset] * xp[0]; + r1 = matrix[(n_rows - 1 - col) * offset] * xm[0]; + } + for (int ind = 1; ind < mid; ++ind) + { + if (transpose_matrix == true) { - res0 = val0 * (in[0] + in[stride * (mm - 1)]); - for (int ind = 1; ind < mid; ++ind) - { - if (contract_over_rows == true) - val0 = shape_hessians[ind * n_columns + n_cols]; - else - val0 = shape_hessians[n_cols * n_columns + ind]; - Number in1 = val0 * (in[stride * ind] + - in[stride * (mm - 1 - ind)]); - res0 += in1; - } + r0 += matrix[ind * offset + col] * xp[ind]; + r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind]; } else - res0 = Number(); - if (mm % 2 == 1) { - if (contract_over_rows == true) - val0 = shape_hessians[mid * n_columns + n_cols]; - else - val0 = shape_hessians[n_cols * n_columns + mid]; - res0 += val0 * in[stride * mid]; + r0 += matrix[col * offset + ind] * xp[ind]; + r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind]; } - if (add) - out[stride * n_cols] += res0; - else - out[stride * n_cols] = res0; } + } + else + r0 = r1 = Number(); + if (mm % 2 == 1 && transpose_matrix == true) + { + if (quantity == EvaluatorQuantity::gradient) + r1 += matrix[mid * offset + col] * xmid; + else + r0 += matrix[mid * offset + col] * xmid; + } + else if (mm % 2 == 1 && + (nn % 2 == 0 || quantity != EvaluatorQuantity::value || + mm == 3)) + r0 += matrix[col * offset + mid] * xmid; + + if (add) + { + out[stride_out * col] += r0 + r1; + if (quantity == EvaluatorQuantity::gradient && + transpose_matrix == false) + out[stride_out * (nn - 1 - col)] += r1 - r0; + else + out[stride_out * (nn - 1 - col)] += r0 - r1; + } + else + { + out[stride_out * col] = r0 + r1; + if (quantity == EvaluatorQuantity::gradient && + transpose_matrix == false) + out[stride_out * (nn - 1 - col)] = r1 - r0; + else + out[stride_out * (nn - 1 - col)] = r0 - r1; + } + } + if (quantity == EvaluatorQuantity::value && transpose_matrix == true && + nn % 2 == 1 && mm % 2 == 1 && mm > 3) + { + if (add) + out[stride_out * n_cols] += matrix[mid * offset + n_cols] * xmid; + else + out[stride_out * n_cols] = matrix[mid * offset + n_cols] * xmid; + } + else if (transpose_matrix == true && nn % 2 == 1) + { + Number r0; + if (mid > 0) + { + r0 = matrix[n_cols] * xp[0]; + for (int ind = 1; ind < mid; ++ind) + r0 += matrix[ind * offset + n_cols] * xp[ind]; + } + else + r0 = Number(); + if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1) + r0 += matrix[mid * offset + n_cols] * xmid; - ++in; - ++out; + if (add) + out[stride_out * n_cols] += r0; + else + out[stride_out * n_cols] = r0; + } + else if (transpose_matrix == false && nn % 2 == 1) + { + Number r0; + if (mid > 0) + { + if (quantity == EvaluatorQuantity::gradient) + { + r0 = matrix[n_cols * offset] * xm[0]; + for (int ind = 1; ind < mid; ++ind) + r0 += matrix[n_cols * offset + ind] * xm[ind]; + } + else + { + r0 = matrix[n_cols * offset] * xp[0]; + for (int ind = 1; ind < mid; ++ind) + r0 += matrix[n_cols * offset + ind] * xp[ind]; + } } - in += stride * (mm - 1); - out += stride * (nn - 1); + else + r0 = Number(); + + if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1) + r0 += matrix[n_cols * offset + mid] * xmid; + + if (add) + out[stride_out * n_cols] += r0; + else + out[stride_out * n_cols] = r0; } } - template - inline void - even_odd_apply(const int n_rows_in, - const int n_columns_in, - const Number2 *DEAL_II_RESTRICT shapes, - const Number * in, - Number * out) + typename Number2> + std::enable_if_t<(variant == evaluate_evenodd), void> + apply_matrix_vector_product(const Number2 *matrix, + const Number * in, + Number * out, + int n_rows, + int n_columns, + int stride_in, + int stride_out) { - static_assert(type < 3, "Only three variants type=0,1,2 implemented"); - static_assert(one_line == false || direction == dim - 1, - "Single-line evaluation only works for direction=dim-1."); - - const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static; - const int n_columns = - n_columns_static == -1 ? n_columns_in : n_columns_static; - - Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || - in != out, - ExcMessage("In-place operation only supported for " - "n_rows==n_columns or single-line interpolation")); + apply_matrix_vector_product( + matrix, in, out, n_rows, n_columns, stride_in, stride_out); + } - // We cannot statically assert that direction is less than dim, so must do - // an additional dynamic check - AssertIndexRange(direction, dim); - const int nn = contract_over_rows ? n_columns : n_rows; - const int mm = contract_over_rows ? n_rows : n_columns; - constexpr int mm_static = - contract_over_rows ? n_rows_static : n_columns_static; - const int n_cols = nn / 2; - const int mid = mm / 2; - constexpr int mid_static = mm_static / 2; - constexpr int max_mid = 15; // for non-templated execution - Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid, - ExcNotImplemented()); + /** + * Internal evaluator specialized for "symmetric" finite elements in the + * symmetric_hierarchical matrix format. + * + * This class implements an approach similar to the even-odd decomposition + * but with a different type of symmetry. In this case, we assume that a + * single shape function already shows the symmetry over the quadrature + * points, rather than the complete basis that is considered in the even-odd + * case. In particular, we assume that the shape functions are ordered as in + * the Legendre basis, with symmetric shape functions in the even slots + * (rows of the values array) and point-symmetric in the odd slots. Like the + * even-odd decomposition, the number of operations are N^2/2 rather than + * N^2 FMAs (fused multiply-add), where N is the number of 1d dofs. The + * difference is in the way the input and output quantities are symmetrized. + */ + template + std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void> + apply_matrix_vector_product(const Number2 *matrix, + const Number * in, + Number * out) + { + static_assert(n_rows > 0 && n_columns > 0, + "Specialization requires n_rows, n_columns > 0"); - const int stride = Utilities::pow(n_columns, direction); - const int n_blocks1 = one_line ? 1 : stride; - const int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + constexpr bool evaluate_antisymmetric = + (quantity == EvaluatorQuantity::gradient); - const int offset = (n_columns + 1) / 2; + constexpr int mm = transpose_matrix ? n_rows : n_columns, + nn = transpose_matrix ? n_columns : n_rows; + constexpr int n_cols = nn / 2; + constexpr int mid = mm / 2; - // this code may look very inefficient at first sight due to the many - // different cases with if's at the innermost loop part, but all of the - // conditionals can be evaluated at compile time because they are - // templates, so the compiler should optimize everything away - for (int i2 = 0; i2 < n_blocks2; ++i2) + if (transpose_matrix) { - for (int i1 = 0; i1 < n_blocks1; ++i1) + std::array x; + for (unsigned int i = 0; i < mm; ++i) + x[i] = in[stride_in * i]; + for (unsigned int col = 0; col < n_cols; ++col) { - constexpr unsigned int mid_size = - (n_rows_static == -1 || n_columns_static == -1) ? - max_mid : - (mid_static > 0 ? mid_static : 1); - Number xp[mid_size], xm[mid_size]; - for (int i = 0; i < mid; ++i) + Number r0, r1; + if (mid > 0) { - if (contract_over_rows == true && type == 1) - { - xp[i] = in[stride * i] - in[stride * (mm - 1 - i)]; - xm[i] = in[stride * i] + in[stride * (mm - 1 - i)]; - } - else + r0 = matrix[col] * x[0]; + r1 = matrix[col + n_columns] * x[1]; + for (unsigned int ind = 1; ind < mid; ++ind) { - xp[i] = in[stride * i] + in[stride * (mm - 1 - i)]; - xm[i] = in[stride * i] - in[stride * (mm - 1 - i)]; + r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind]; + r1 += + matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1]; } } - Number xmid = in[stride * mid]; - for (int col = 0; col < n_cols; ++col) + else + r0 = r1 = Number(); + if (mm % 2 == 1) + r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1]; + if (add) { - Number r0, r1; - if (mid > 0) - { - if (contract_over_rows == true) - { - r0 = shapes[col] * xp[0]; - r1 = shapes[(n_rows - 1) * offset + col] * xm[0]; - } - else - { - r0 = shapes[col * offset] * xp[0]; - r1 = shapes[(n_rows - 1 - col) * offset] * xm[0]; - } - for (int ind = 1; ind < mid; ++ind) - { - if (contract_over_rows == true) - { - r0 += shapes[ind * offset + col] * xp[ind]; - r1 += shapes[(n_rows - 1 - ind) * offset + col] * - xm[ind]; - } - else - { - r0 += shapes[col * offset + ind] * xp[ind]; - r1 += shapes[(n_rows - 1 - col) * offset + ind] * - xm[ind]; - } - } - } + out[stride_out * col] += r0 + r1; + if (evaluate_antisymmetric) + out[stride_out * (nn - 1 - col)] += r1 - r0; else - r0 = r1 = Number(); - if (mm % 2 == 1 && contract_over_rows == true) - { - if (type == 1) - r1 += shapes[mid * offset + col] * xmid; - else - r0 += shapes[mid * offset + col] * xmid; - } - else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3)) - r0 += shapes[col * offset + mid] * xmid; - - if (add) - { - out[stride * col] += r0 + r1; - if (type == 1 && contract_over_rows == false) - out[stride * (nn - 1 - col)] += r1 - r0; - else - out[stride * (nn - 1 - col)] += r0 - r1; - } + out[stride_out * (nn - 1 - col)] += r0 - r1; + } + else + { + out[stride_out * col] = r0 + r1; + if (evaluate_antisymmetric) + out[stride_out * (nn - 1 - col)] = r1 - r0; else + out[stride_out * (nn - 1 - col)] = r0 - r1; + } + } + if (nn % 2 == 1) + { + Number r0; + const unsigned int shift = evaluate_antisymmetric ? 1 : 0; + if (mid > 0) + { + r0 = matrix[n_cols + shift * n_columns] * x[shift]; + for (unsigned int ind = 1; ind < mid; ++ind) + r0 += matrix[n_cols + (2 * ind + shift) * n_columns] * + x[2 * ind + shift]; + } + else + r0 = 0; + if (!evaluate_antisymmetric && mm % 2 == 1) + r0 += matrix[n_cols + (mm - 1) * n_columns] * x[mm - 1]; + if (add) + out[stride_out * n_cols] += r0; + else + out[stride_out * n_cols] = r0; + } + } + else + { + std::array xp, xm; + for (int i = 0; i < mid; ++i) + if (!evaluate_antisymmetric) + { + xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)]; + xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)]; + } + else + { + xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)]; + xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)]; + } + if (mm % 2 == 1) + xp[mid] = in[stride_in * mid]; + for (unsigned int col = 0; col < n_cols; ++col) + { + Number r0, r1; + if (mid > 0) + { + r0 = matrix[2 * col * n_columns] * xp[0]; + r1 = matrix[(2 * col + 1) * n_columns] * xm[0]; + for (unsigned int ind = 1; ind < mid; ++ind) { - out[stride * col] = r0 + r1; - if (type == 1 && contract_over_rows == false) - out[stride * (nn - 1 - col)] = r1 - r0; - else - out[stride * (nn - 1 - col)] = r0 - r1; + r0 += matrix[2 * col * n_columns + ind] * xp[ind]; + r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind]; } } - if (type == 0 && contract_over_rows == true && nn % 2 == 1 && - mm % 2 == 1 && mm > 3) - { - if (add) - out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid; - else - out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid; - } - else if (contract_over_rows == true && nn % 2 == 1) + else + r0 = r1 = Number(); + if (mm % 2 == 1) { - Number r0; - if (mid > 0) - { - r0 = shapes[n_cols] * xp[0]; - for (int ind = 1; ind < mid; ++ind) - r0 += shapes[ind * offset + n_cols] * xp[ind]; - } - else - r0 = Number(); - if (type != 1 && mm % 2 == 1) - r0 += shapes[mid * offset + n_cols] * xmid; - - if (add) - out[stride * n_cols] += r0; + if (evaluate_antisymmetric) + r1 += matrix[(2 * col + 1) * n_columns + mid] * xp[mid]; else - out[stride * n_cols] = r0; + r0 += matrix[2 * col * n_columns + mid] * xp[mid]; } - else if (contract_over_rows == false && nn % 2 == 1) + if (add) { - Number r0; - if (mid > 0) - { - if (type == 1) - { - r0 = shapes[n_cols * offset] * xm[0]; - for (int ind = 1; ind < mid; ++ind) - r0 += shapes[n_cols * offset + ind] * xm[ind]; - } - else - { - r0 = shapes[n_cols * offset] * xp[0]; - for (int ind = 1; ind < mid; ++ind) - r0 += shapes[n_cols * offset + ind] * xp[ind]; - } - } - else - r0 = Number(); - - if ((type == 0 || type == 2) && mm % 2 == 1) - r0 += shapes[n_cols * offset + mid] * xmid; - - if (add) - out[stride * n_cols] += r0; - else - out[stride * n_cols] = r0; + out[stride_out * (2 * col)] += r0; + out[stride_out * (2 * col + 1)] += r1; } - if (one_line == false) + else { - in += 1; - out += 1; + out[stride_out * (2 * col)] = r0; + out[stride_out * (2 * col + 1)] = r1; } } - if (one_line == false) + if (nn % 2 == 1) { - in += stride * (mm - 1); - out += stride * (nn - 1); + Number r0; + if (mid > 0) + { + r0 = matrix[(nn - 1) * n_columns] * xp[0]; + for (unsigned int ind = 1; ind < mid; ++ind) + r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind]; + } + else + r0 = Number(); + if (mm % 2 == 1 && !evaluate_antisymmetric) + r0 += matrix[(nn - 1) * n_columns + mid] * xp[mid]; + if (add) + out[stride_out * (nn - 1)] += r0; + else + out[stride_out * (nn - 1)] = r0; } } } @@ -1795,22 +1057,13 @@ namespace internal /** - * Internal evaluator for 1d-3d shape function using the tensor product form - * of the basis functions. - * - * This class implements a different approach to the symmetric case for - * values, gradients, and Hessians also treated with the above functions: It - * is possible to reduce the cost per dimension from N^2 to N^2/2, where N - * is the number of 1d dofs (there are only N^2/2 different entries in the - * shape matrix, so this is plausible). The approach is based on the idea of - * applying the operator on the even and odd part of the input vectors - * separately, given that the shape functions evaluated on quadrature points - * are symmetric. This method is presented e.g. in the book "Implementing - * Spectral Methods for Partial Differential Equations" by David A. Kopriva, - * Springer, 2009, section 3.5.3 (Even-Odd-Decomposition). Even though the - * experiments in the book say that the method is not efficient for N<20, it - * is more efficient in the context where the loop bounds are compile-time - * constants (templates). + * Generic evaluator framework that valuates the given shape data in general + * dimensions using the tensor product form. Depending on the particular + * layout in the matrix entries, this corresponds to a usual matrix-matrix + * product or a matrix-matrix product including some symmetries. The actual + * work is implemented by functions of type apply_matrix_vector_product + * working on a single dimension, controlled by suitable strides, using the + * kernel specified via variant. * * @tparam dim Space dimension in which this class is applied * @tparam n_rows Number of rows in the transformation matrix, which corresponds @@ -1825,17 +1078,13 @@ namespace internal * operator* with Number and produce Number as an output to * be a valid type */ - template - struct EvaluatorTensorProduct + typename Number2 = Number> + struct EvaluatorTensorProduct { static constexpr unsigned int n_rows_of_product = Utilities::pow(n_rows, dim); @@ -1844,8 +1093,7 @@ namespace internal /** * Empty constructor. Does nothing. Be careful when using 'values' and - * related methods because they need to be filled with the other - * constructor passing in at least an array for the values. + * related methods because they need to be filled with the other pointer */ EvaluatorTensorProduct() : shape_values(nullptr) @@ -1854,38 +1102,62 @@ namespace internal {} /** - * Constructor, taking the data from ShapeInfo (using the even-odd - * variants stored there) - */ - EvaluatorTensorProduct(const AlignedVector &shape_values) - : shape_values(shape_values.begin()) - , shape_gradients(nullptr) - , shape_hessians(nullptr) - { - AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2)); - } - - /** - * Constructor, taking the data from ShapeInfo (using the even-odd - * variants stored there) + * Constructor, taking the data from ShapeInfo */ EvaluatorTensorProduct(const AlignedVector &shape_values, const AlignedVector &shape_gradients, const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) + const unsigned int = 0, + const unsigned int = 0) : shape_values(shape_values.begin()) , shape_gradients(shape_gradients.begin()) , shape_hessians(shape_hessians.begin()) { - // In this function, we allow for dummy pointers if some of values, - // gradients or hessians should not be computed - if (!shape_values.empty()) - AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2)); - if (!shape_gradients.empty()) - AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2)); - if (!shape_hessians.empty()) - AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2)); + if (variant == evaluate_evenodd) + { + if (!shape_values.empty()) + AssertDimension(shape_values.size(), + n_rows * ((n_columns + 1) / 2)); + if (!shape_gradients.empty()) + AssertDimension(shape_gradients.size(), + n_rows * ((n_columns + 1) / 2)); + if (!shape_hessians.empty()) + AssertDimension(shape_hessians.size(), + n_rows * ((n_columns + 1) / 2)); + } + else + { + // We can enter this function either for the apply() path that has + // n_rows * n_columns entries or for the apply_face() path that only + // has n_rows * 3 entries in the array. Since we cannot decide about + // the use we must allow for both here. + Assert(shape_values.empty() || + shape_values.size() == n_rows * n_columns || + shape_values.size() == 3 * n_rows, + ExcDimensionMismatch(shape_values.size(), n_rows * n_columns)); + Assert(shape_gradients.empty() || + shape_gradients.size() == n_rows * n_columns, + ExcDimensionMismatch(shape_gradients.size(), + n_rows * n_columns)); + Assert(shape_hessians.empty() || + shape_hessians.size() == n_rows * n_columns, + ExcDimensionMismatch(shape_hessians.size(), + n_rows * n_columns)); + } + } + + /** + * Constructor, taking the data from ShapeInfo via raw pointers + */ + EvaluatorTensorProduct(const Number2 * shape_values, + const Number2 * shape_gradients, + const Number2 * shape_hessians, + const unsigned int dummy1 = 0, + const unsigned int dummy2 = 0) + : shape_values(shape_values) + , shape_gradients(shape_gradients) + , shape_hessians(shape_hessians) + { (void)dummy1; (void)dummy2; } @@ -1894,24 +1166,35 @@ namespace internal void values(const Number in[], Number out[]) const { - Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); + apply(shape_values, in, out); } template void gradients(const Number in[], Number out[]) const { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, in, out); + constexpr EvaluatorQuantity gradient_type = + (variant == evaluate_general ? EvaluatorQuantity::value : + EvaluatorQuantity::gradient); + apply( + shape_gradients, in, out); } template void hessians(const Number in[], Number out[]) const { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, in, out); + constexpr EvaluatorQuantity hessian_type = + (((variant == evaluate_general) | + (variant == evaluate_symmetric_hierarchical)) ? + EvaluatorQuantity::value : + EvaluatorQuantity::hessian); + apply( + shape_hessians, in, out); } template @@ -1919,7 +1202,8 @@ namespace internal values_one_line(const Number in[], Number out[]) const { Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); + apply( + shape_values, in, out); } template @@ -1927,9 +1211,11 @@ namespace internal gradients_one_line(const Number in[], Number out[]) const { Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, - in, - out); + constexpr EvaluatorQuantity gradient_type = + (variant == evaluate_general ? EvaluatorQuantity::value : + EvaluatorQuantity::gradient); + apply( + shape_gradients, in, out); } template @@ -1937,9 +1223,13 @@ namespace internal hessians_one_line(const Number in[], Number out[]) const { Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, - in, - out); + constexpr EvaluatorQuantity hessian_type = + (((variant == evaluate_general) | + (variant == evaluate_symmetric_hierarchical)) ? + EvaluatorQuantity::value : + EvaluatorQuantity::hessian); + apply( + shape_hessians, in, out); } /** @@ -1948,8 +1238,8 @@ namespace internal * data in the input array. This function allows the @p in and @p out * arrays to alias for the case n_rows == n_columns, i.e., it is safe to * perform the contraction in place where @p in and @p out point to the - * same address. For the case n_rows != n_columns, the output is only - * correct if @p one_line is set to true. + * same address. For the case n_rows != n_columns, the output is in general + * not correct. * * @tparam direction Direction that is evaluated * @tparam contract_over_rows If true, the tensor contraction sums @@ -1957,10 +1247,6 @@ namespace internal * array, otherwise it sums over the columns * @tparam add If true, the result is added to the output vector, else * the computed values overwrite the content in the output - * @tparam type Determines whether to use the symmetries appearing in - * shape values (type=0), shape gradients (type=1) or - * second derivatives (type=2, similar to type 0 but - * without two additional zero entries) * @tparam one_line If true, the kernel is only applied along a single 1d * stripe within a dim-dimensional tensor, not the full * n_rows^dim points as in the @p false case. @@ -1973,24 +1259,49 @@ namespace internal template + bool one_line = false, + EvaluatorQuantity = EvaluatorQuantity::value> static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number * in, - Number * out) - { - even_odd_apply(n_rows, n_columns, shape_data, in, out); - } + Number * out); + + /** + * This function applies the tensor product operation to produce face values + * from cell values. As opposed to the apply method, this method assumes + * that the directions orthogonal to the face have n_rows degrees of + * freedom per direction and not n_columns for those directions lower than + * the one currently applied. In other words, apply_face() must be called + * before calling any interpolation within the face. + * + * @tparam face_direction Direction of the normal vector (0=x, 1=y, etc) + * @tparam contract_onto_face If true, the input vector is of size n_rows^dim + * and interpolation into n_rows^(dim-1) points + * is performed. This is a typical scenario in + * FEFaceEvaluation::evaluate() calls. If false, + * data from n_rows^(dim-1) points is expanded + * into the n_rows^dim points of the higher- + * dimensional data array. Derivatives in the + * case contract_onto_face==false are summed + * together + * @tparam add If true, the result is added to the output vector, else + * the computed values overwrite the content in the output + * @tparam max_derivative Sets the number of derivatives that should be + * computed. 0 means only values, 1 means values and first + * derivatives, 2 second derivates. Note that all the + * derivatives access the data in @p shape_values passed to + * the constructor of the class + * + * @param in address of the input data vector + * @param out address of the output data vector + */ + template + void + apply_face(const Number *DEAL_II_RESTRICT in, + Number *DEAL_II_RESTRICT out) const; private: const Number2 *shape_values; @@ -1999,243 +1310,336 @@ namespace internal }; - /** - * Internal evaluator for shape function using the tensor product form - * of the basis functions. The same as the other templated class but - * without making use of template arguments and variable loop bounds - * instead. - */ - template - struct EvaluatorTensorProduct + + template + template + inline void + EvaluatorTensorProduct:: + apply(const Number2 *DEAL_II_RESTRICT shape_data, + const Number * in, + Number * out) { - EvaluatorTensorProduct() - : shape_values(nullptr) - , shape_gradients(nullptr) - , shape_hessians(nullptr) - , n_rows(numbers::invalid_unsigned_int) - , n_columns(numbers::invalid_unsigned_int) - {} + static_assert(one_line == false || direction == dim - 1, + "Single-line evaluation only works for direction=dim-1."); + Assert(shape_data != nullptr, + ExcMessage( + "The given array shape_data must not be the null pointer!")); + Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || + in != out, + ExcMessage("In-place operation only supported for " + "n_rows==n_columns or single-line interpolation")); + AssertIndexRange(direction, dim); + constexpr int mm = contract_over_rows ? n_rows : n_columns, + nn = contract_over_rows ? n_columns : n_rows; - EvaluatorTensorProduct(const AlignedVector &shape_values, - const unsigned int n_rows = 0, - const unsigned int n_columns = 0) - : shape_values(shape_values.begin()) - , shape_gradients(nullptr) - , shape_hessians(nullptr) - , n_rows(n_rows) - , n_columns(n_columns) - { - AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2)); - } + constexpr int stride = Utilities::pow(n_columns, direction); + constexpr int n_blocks1 = one_line ? 1 : stride; + constexpr int n_blocks2 = + Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); - EvaluatorTensorProduct(const AlignedVector &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int n_rows = 0, - const unsigned int n_columns = 0) - : shape_values(shape_values.begin()) - , shape_gradients(shape_gradients.begin()) - , shape_hessians(shape_hessians.begin()) - , n_rows(n_rows) - , n_columns(n_columns) - { - if (!shape_values.empty()) - AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2)); - if (!shape_gradients.empty()) - AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2)); - if (!shape_hessians.empty()) - AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2)); - } + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + apply_matrix_vector_product(shape_data, in, out); - template - void - values(const Number in[], Number out[]) const - { - Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); - } + if (one_line == false) + { + ++in; + ++out; + } + } + if (one_line == false) + { + in += stride * (mm - 1); + out += stride * (nn - 1); + } + } + } - template - void - gradients(const Number in[], Number out[]) const - { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, in, out); - } - template - void - hessians(const Number in[], Number out[]) const - { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, in, out); - } - template - void - values_one_line(const Number in[], Number out[]) const - { - Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); - } + template + template + inline void + EvaluatorTensorProduct:: + apply_face(const Number *DEAL_II_RESTRICT in, + Number *DEAL_II_RESTRICT out) const + { + Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported")); + static_assert(max_derivative >= 0 && max_derivative < 3, + "Only derivative orders 0-2 implemented"); + Assert(shape_values != nullptr, + ExcMessage( + "The given array shape_values must not be the null pointer.")); - template - void - gradients_one_line(const Number in[], Number out[]) const - { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, - in, - out); - } + constexpr int n_blocks1 = (dim > 1 ? n_rows : 1); + constexpr int n_blocks2 = (dim > 2 ? n_rows : 1); - template - void - hessians_one_line(const Number in[], Number out[]) const - { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, - in, - out); - } + AssertIndexRange(face_direction, dim); + constexpr int in_stride = Utilities::pow(n_rows, face_direction); + constexpr int out_stride = Utilities::pow(n_rows, dim - 1); + const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values; - template - void - apply(const Number2 *DEAL_II_RESTRICT shape_data, - const Number * in, - Number * out) const - { - even_odd_apply(n_rows, n_columns, shape_data, in, out); - } + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + if (contract_onto_face == true) + { + Number res0 = shape_values[0] * in[0]; + Number res1, res2; + if (max_derivative > 0) + res1 = shape_values[n_rows] * in[0]; + if (max_derivative > 1) + res2 = shape_values[2 * n_rows] * in[0]; + for (int ind = 1; ind < n_rows; ++ind) + { + res0 += shape_values[ind] * in[in_stride * ind]; + if (max_derivative > 0) + res1 += shape_values[ind + n_rows] * in[in_stride * ind]; + if (max_derivative > 1) + res2 += + shape_values[ind + 2 * n_rows] * in[in_stride * ind]; + } + if (add) + { + out[0] += res0; + if (max_derivative > 0) + out[out_stride] += res1; + if (max_derivative > 1) + out[2 * out_stride] += res2; + } + else + { + out[0] = res0; + if (max_derivative > 0) + out[out_stride] = res1; + if (max_derivative > 1) + out[2 * out_stride] = res2; + } + } + else + { + for (int col = 0; col < n_rows; ++col) + { + if (add) + out[col * in_stride] += shape_values[col] * in[0]; + else + out[col * in_stride] = shape_values[col] * in[0]; + if (max_derivative > 0) + out[col * in_stride] += + shape_values[col + n_rows] * in[out_stride]; + if (max_derivative > 1) + out[col * in_stride] += + shape_values[col + 2 * n_rows] * in[2 * out_stride]; + } + } - private: - const Number2 * shape_values; - const Number2 * shape_gradients; - const Number2 * shape_hessians; - const unsigned int n_rows; - const unsigned int n_columns; - }; + // increment: in regular case, just go to the next point in + // x-direction. If we are at the end of one chunk in x-dir, need + // to jump over to the next layer in z-direction + switch (face_direction) + { + case 0: + in += contract_onto_face ? n_rows : 1; + out += contract_onto_face ? 1 : n_rows; + break; + case 1: + ++in; + ++out; + // faces 2 and 3 in 3d use local coordinate system zx, which + // is the other way around compared to the tensor + // product. Need to take that into account. + if (dim == 3) + { + if (contract_onto_face) + out += n_rows - 1; + else + in += n_rows - 1; + } + break; + case 2: + ++in; + ++out; + break; + default: + Assert(false, ExcNotImplemented()); + } + } + + // adjust for local coordinate system zx + if (face_direction == 1 && dim == 3) + { + if (contract_onto_face) + { + in += n_rows * (n_rows - 1); + out -= n_rows * n_rows - 1; + } + else + { + out += n_rows * (n_rows - 1); + in -= n_rows * n_rows - 1; + } + } + } + } /** - * Internal evaluator for 1d-3d shape function using the tensor product form - * of the basis functions. - * - * This class implements an approach similar to the even-odd decomposition - * but with a different type of symmetry. In this case, we assume that a - * single shape function already shows the symmetry over the quadrature - * points, rather than the complete basis that is considered in the even-odd - * case. In particular, we assume that the shape functions are ordered as in - * the Legendre basis, with symmetric shape functions in the even slots - * (rows of the values array) and point-symmetric in the odd slots. Like the - * even-odd decomposition, the number of operations are N^2/2 rather than - * N^2 FMAs (fused multiply-add), where N is the number of 1d dofs. The - * difference is in the way the input and output quantities are symmetrized. + * Internal evaluator for shape function using the tensor product form + * of the basis functions. The same as the other templated class but + * without making use of template arguments and variable loop bounds + * instead. * * @tparam dim Space dimension in which this class is applied - * @tparam n_rows Number of rows in the transformation matrix, which corresponds - * to the number of 1d shape functions in the usual tensor - * contraction setting - * @tparam n_columns Number of columns in the transformation matrix, which - * corresponds to the number of 1d shape functions in the - * usual tensor contraction setting * @tparam Number Abstract number type for input and output arrays * @tparam Number2 Abstract number type for coefficient arrays (defaults to * same type as the input/output arrays); must implement * operator* with Number and produce Number as an output to * be a valid type */ - template - struct EvaluatorTensorProduct + struct EvaluatorTensorProduct { static constexpr unsigned int n_rows_of_product = - Utilities::pow(n_rows, dim); + numbers::invalid_unsigned_int; static constexpr unsigned int n_columns_of_product = - Utilities::pow(n_columns, dim); + numbers::invalid_unsigned_int; /** * Empty constructor. Does nothing. Be careful when using 'values' and - * related methods because they need to be filled with the other - * constructor passing in at least an array for the values. + * related methods because they need to be filled with the other constructor */ EvaluatorTensorProduct() : shape_values(nullptr) , shape_gradients(nullptr) , shape_hessians(nullptr) + , n_rows(numbers::invalid_unsigned_int) + , n_columns(numbers::invalid_unsigned_int) {} /** - * Constructor, taking the data from ShapeInfo (using the even-odd - * variants stored there) - */ - EvaluatorTensorProduct(const AlignedVector &shape_values) - : shape_values(shape_values.begin()) - , shape_gradients(nullptr) - , shape_hessians(nullptr) - {} - - /** - * Constructor, taking the data from ShapeInfo (using the even-odd - * variants stored there) + * Constructor, taking the data from ShapeInfo */ EvaluatorTensorProduct(const AlignedVector &shape_values, const AlignedVector &shape_gradients, const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) + const unsigned int n_rows = 0, + const unsigned int n_columns = 0) : shape_values(shape_values.begin()) , shape_gradients(shape_gradients.begin()) , shape_hessians(shape_hessians.begin()) + , n_rows(n_rows) + , n_columns(n_columns) { - (void)dummy1; - (void)dummy2; + if (variant == evaluate_evenodd) + { + if (!shape_values.empty()) + AssertDimension(shape_values.size(), + n_rows * ((n_columns + 1) / 2)); + if (!shape_gradients.empty()) + AssertDimension(shape_gradients.size(), + n_rows * ((n_columns + 1) / 2)); + if (!shape_hessians.empty()) + AssertDimension(shape_hessians.size(), + n_rows * ((n_columns + 1) / 2)); + } + else + { + // We can enter this function either for the apply() path that has + // n_rows * n_columns entries or for the apply_face() path that only + // has n_rows * 3 entries in the array. Since we cannot decide about + // the use we must allow for both here. + Assert(shape_values.empty() || + shape_values.size() == n_rows * n_columns || + shape_values.size() == n_rows * 3, + ExcDimensionMismatch(shape_values.size(), n_rows * n_columns)); + Assert(shape_gradients.empty() || + shape_gradients.size() == n_rows * n_columns, + ExcDimensionMismatch(shape_gradients.size(), + n_rows * n_columns)); + Assert(shape_hessians.empty() || + shape_hessians.size() == n_rows * n_columns, + ExcDimensionMismatch(shape_hessians.size(), + n_rows * n_columns)); + } } + /** + * Constructor, taking the data from ShapeInfo + */ + EvaluatorTensorProduct(const Number2 * shape_values, + const Number2 * shape_gradients, + const Number2 * shape_hessians, + const unsigned int n_rows = 0, + const unsigned int n_columns = 0) + : shape_values(shape_values) + , shape_gradients(shape_gradients) + , shape_hessians(shape_hessians) + , n_rows(n_rows) + , n_columns(n_columns) + {} + template void - values(const Number in[], Number out[]) const + values(const Number *in, Number *out) const { - Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); + apply(shape_values, in, out); } template void - gradients(const Number in[], Number out[]) const + gradients(const Number *in, Number *out) const { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, in, out); + constexpr EvaluatorQuantity gradient_type = + (variant != evaluate_evenodd ? EvaluatorQuantity::value : + EvaluatorQuantity::gradient); + apply( + shape_gradients, in, out); } template void - hessians(const Number in[], Number out[]) const + hessians(const Number *in, Number *out) const { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, in, out); + constexpr EvaluatorQuantity hessian_type = + (variant != evaluate_evenodd ? EvaluatorQuantity::value : + EvaluatorQuantity::hessian); + apply( + shape_hessians, in, out); } template @@ -2243,274 +1647,297 @@ namespace internal values_one_line(const Number in[], Number out[]) const { Assert(shape_values != nullptr, ExcNotInitialized()); - apply(shape_values, in, out); - } - - template - void - gradients_one_line(const Number in[], Number out[]) const - { - Assert(shape_gradients != nullptr, ExcNotInitialized()); - apply(shape_gradients, - in, - out); + apply( + shape_values, in, out); } template - void - hessians_one_line(const Number in[], Number out[]) const - { - Assert(shape_hessians != nullptr, ExcNotInitialized()); - apply(shape_hessians, - in, - out); - } - - /** - * This function applies the tensor product kernel, corresponding to a - * multiplication of 1d stripes, along the given @p direction of the tensor - * data in the input array. This function allows the @p in and @p out - * arrays to alias for the case n_rows == n_columns, i.e., it is safe to - * perform the contraction in place where @p in and @p out point to the - * same address. For the case n_rows != n_columns, the output is only - * correct if @p one_line is set to true. - * - * @tparam direction Direction that is evaluated - * @tparam contract_over_rows If true, the tensor contraction sums - * over the rows in the given @p shape_data - * array, otherwise it sums over the columns - * @tparam add If true, the result is added to the output vector, else - * the computed values overwrite the content in the output - * @tparam type Determines whether the evaluation is symmetric in even - * rows (type=0) or odd rows (type=1) of @p shape_data and - * skew-symmetric in odd rows (type=0) or even rows (type=1) - * @tparam one_line If true, the kernel is only applied along a single 1d - * stripe within a dim-dimensional tensor, not the full - * n_rows^dim points as in the @p false case. - * - * @param shape_data Transformation matrix with @p n_rows rows and - * @p n_columns columns, stored in row-major format - * @param in Pointer to the start of the input data vector - * @param out Pointer to the start of the output data vector - */ - template - static void + void + gradients_one_line(const Number in[], Number out[]) const + { + Assert(shape_gradients != nullptr, ExcNotInitialized()); + constexpr EvaluatorQuantity gradient_type = + (variant != evaluate_evenodd ? EvaluatorQuantity::value : + EvaluatorQuantity::gradient); + apply( + shape_gradients, in, out); + } + + template + void + hessians_one_line(const Number in[], Number out[]) const + { + Assert(shape_hessians != nullptr, ExcNotInitialized()); + constexpr EvaluatorQuantity hessian_type = + (variant != evaluate_evenodd ? EvaluatorQuantity::value : + EvaluatorQuantity::hessian); + apply( + shape_hessians, in, out); + } + + template + void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number * in, - Number * out); + Number * out) const; - private: - const Number2 *shape_values; - const Number2 *shape_gradients; - const Number2 *shape_hessians; + template + void + apply_face(const Number *DEAL_II_RESTRICT in, + Number *DEAL_II_RESTRICT out) const; + + const Number2 * shape_values; + const Number2 * shape_gradients; + const Number2 * shape_hessians; + const unsigned int n_rows; + const unsigned int n_columns; }; - template - template + template inline void - EvaluatorTensorProduct::apply(const Number2 *DEAL_II_RESTRICT shapes, - const Number * in, - Number * out) + EvaluatorTensorProduct::apply( + const Number2 *DEAL_II_RESTRICT shape_data, + const Number * in, + Number * out) const { static_assert(one_line == false || direction == dim - 1, "Single-line evaluation only works for direction=dim-1."); - static_assert( - type == 0 || type == 1, - "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical."); + Assert(shape_data != nullptr, + ExcMessage( + "The given array shape_data must not be the null pointer!")); Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || in != out, ExcMessage("In-place operation only supported for " "n_rows==n_columns or single-line interpolation")); - - // We cannot statically assert that direction is less than dim, so must do - // an additional dynamic check AssertIndexRange(direction, dim); + const int mm = contract_over_rows ? n_rows : n_columns, + nn = contract_over_rows ? n_columns : n_rows; - constexpr int nn = contract_over_rows ? n_columns : n_rows; - constexpr int mm = contract_over_rows ? n_rows : n_columns; - constexpr int n_cols = nn / 2; - constexpr int mid = mm / 2; + const int stride = + direction == 0 ? 1 : Utilities::fixed_power(n_columns); + const int n_blocks1 = one_line ? 1 : stride; + const int n_blocks2 = direction >= dim - 1 ? + 1 : + Utilities::fixed_power(n_rows); + Assert(n_rows <= 128, ExcNotImplemented()); - constexpr int stride = Utilities::pow(n_columns, direction); - constexpr int n_blocks1 = one_line ? 1 : stride; - constexpr int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + // the empty template case can only run the general evaluator or + // evenodd + constexpr EvaluatorVariant restricted_variant = + variant == evaluate_evenodd ? evaluate_evenodd : evaluate_general; + apply_matrix_vector_product( + shape_data, in, out, n_rows, n_columns, stride, stride); + + if (one_line == false) + { + ++in; + ++out; + } + } + if (one_line == false) + { + in += stride * (mm - 1); + out += stride * (nn - 1); + } + } + } + + + + template + template + inline void + EvaluatorTensorProduct::apply_face( + const Number *DEAL_II_RESTRICT in, + Number *DEAL_II_RESTRICT out) const + { + Assert(shape_values != nullptr, + ExcMessage( + "The given array shape_data must not be the null pointer!")); + static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported"); + const int n_blocks1 = dim > 1 ? n_rows : 1; + const int n_blocks2 = dim > 2 ? n_rows : 1; + + AssertIndexRange(face_direction, dim); + const int in_stride = + face_direction > 0 ? Utilities::fixed_power(n_rows) : 1; + const int out_stride = + dim > 1 ? Utilities::fixed_power(n_rows) : 1; - // this code may look very inefficient at first sight due to the many - // different cases with if's at the innermost loop part, but all of the - // conditionals can be evaluated at compile time because they are - // templates, so the compiler should optimize everything away for (int i2 = 0; i2 < n_blocks2; ++i2) { for (int i1 = 0; i1 < n_blocks1; ++i1) { - if (contract_over_rows) + if (contract_onto_face == true) { - Number x[mm]; - for (unsigned int i = 0; i < mm; ++i) - x[i] = in[stride * i]; - for (unsigned int col = 0; col < n_cols; ++col) + Number res0 = shape_values[0] * in[0]; + Number res1, res2; + if (max_derivative > 0) + res1 = shape_values[n_rows] * in[0]; + if (max_derivative > 1) + res2 = shape_values[2 * n_rows] * in[0]; + for (unsigned int ind = 1; ind < n_rows; ++ind) { - Number r0, r1; - if (mid > 0) - { - r0 = shapes[col] * x[0]; - r1 = shapes[col + n_columns] * x[1]; - for (unsigned int ind = 1; ind < mid; ++ind) - { - r0 += - shapes[col + 2 * ind * n_columns] * x[2 * ind]; - r1 += shapes[col + (2 * ind + 1) * n_columns] * - x[2 * ind + 1]; - } - } - else - r0 = r1 = Number(); - if (mm % 2 == 1) - r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1]; - if (add) - { - out[stride * col] += r0 + r1; - if (type == 1) - out[stride * (nn - 1 - col)] += r1 - r0; - else - out[stride * (nn - 1 - col)] += r0 - r1; - } - else - { - out[stride * col] = r0 + r1; - if (type == 1) - out[stride * (nn - 1 - col)] = r1 - r0; - else - out[stride * (nn - 1 - col)] = r0 - r1; - } + res0 += shape_values[ind] * in[in_stride * ind]; + if (max_derivative > 0) + res1 += shape_values[ind + n_rows] * in[in_stride * ind]; + if (max_derivative > 1) + res2 += + shape_values[ind + 2 * n_rows] * in[in_stride * ind]; + } + if (add) + { + out[0] += res0; + if (max_derivative > 0) + out[out_stride] += res1; + if (max_derivative > 1) + out[2 * out_stride] += res2; } - if (nn % 2 == 1) + else { - Number r0; - const unsigned int shift = type == 1 ? 1 : 0; - if (mid > 0) - { - r0 = shapes[n_cols + shift * n_columns] * x[shift]; - for (unsigned int ind = 1; ind < mid; ++ind) - r0 += shapes[n_cols + (2 * ind + shift) * n_columns] * - x[2 * ind + shift]; - } - else - r0 = 0; - if (type != 1 && mm % 2 == 1) - r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1]; - if (add) - out[stride * n_cols] += r0; - else - out[stride * n_cols] = r0; + out[0] = res0; + if (max_derivative > 0) + out[out_stride] = res1; + if (max_derivative > 1) + out[2 * out_stride] = res2; } } else { - Number xp[mid + 1], xm[mid > 0 ? mid : 1]; - for (int i = 0; i < mid; ++i) - if (type == 0) - { - xp[i] = in[stride * i] + in[stride * (mm - 1 - i)]; - xm[i] = in[stride * i] - in[stride * (mm - 1 - i)]; - } - else - { - xp[i] = in[stride * i] - in[stride * (mm - 1 - i)]; - xm[i] = in[stride * i] + in[stride * (mm - 1 - i)]; - } - if (mm % 2 == 1) - xp[mid] = in[stride * mid]; - for (unsigned int col = 0; col < n_cols; ++col) - { - Number r0, r1; - if (mid > 0) - { - r0 = shapes[2 * col * n_columns] * xp[0]; - r1 = shapes[(2 * col + 1) * n_columns] * xm[0]; - for (unsigned int ind = 1; ind < mid; ++ind) - { - r0 += shapes[2 * col * n_columns + ind] * xp[ind]; - r1 += - shapes[(2 * col + 1) * n_columns + ind] * xm[ind]; - } - } - else - r0 = r1 = Number(); - if (mm % 2 == 1) - { - if (type == 1) - r1 += - shapes[(2 * col + 1) * n_columns + mid] * xp[mid]; - else - r0 += shapes[2 * col * n_columns + mid] * xp[mid]; - } - if (add) - { - out[stride * (2 * col)] += r0; - out[stride * (2 * col + 1)] += r1; - } - else - { - out[stride * (2 * col)] = r0; - out[stride * (2 * col + 1)] = r1; - } - } - if (nn % 2 == 1) + for (unsigned int col = 0; col < n_rows; ++col) { - Number r0; - if (mid > 0) - { - r0 = shapes[(nn - 1) * n_columns] * xp[0]; - for (unsigned int ind = 1; ind < mid; ++ind) - r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind]; - } - else - r0 = Number(); - if (mm % 2 == 1 && type == 0) - r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid]; if (add) - out[stride * (nn - 1)] += r0; + out[col * in_stride] += shape_values[col] * in[0]; else - out[stride * (nn - 1)] = r0; + out[col * in_stride] = shape_values[col] * in[0]; + if (max_derivative > 0) + out[col * in_stride] += + shape_values[col + n_rows] * in[out_stride]; + if (max_derivative > 1) + out[col * in_stride] += + shape_values[col + 2 * n_rows] * in[2 * out_stride]; } } - if (one_line == false) + + // increment: in regular case, just go to the next point in + // x-direction. If we are at the end of one chunk in x-dir, need + // to jump over to the next layer in z-direction + switch (face_direction) { - in += 1; - out += 1; + case 0: + in += contract_onto_face ? n_rows : 1; + out += contract_onto_face ? 1 : n_rows; + break; + case 1: + ++in; + ++out; + // faces 2 and 3 in 3d use local coordinate system zx, which + // is the other way around compared to the tensor + // product. Need to take that into account. + if (dim == 3) + { + if (contract_onto_face) + out += n_rows - 1; + else + in += n_rows - 1; + } + break; + case 2: + ++in; + ++out; + break; + default: + Assert(false, ExcNotImplemented()); } } - if (one_line == false) + if (face_direction == 1 && dim == 3) { - in += stride * (mm - 1); - out += stride * (nn - 1); + // adjust for local coordinate system zx + if (contract_onto_face) + { + in += n_rows * (n_rows - 1); + out -= n_rows * n_rows - 1; + } + else + { + out += n_rows * (n_rows - 1); + in -= n_rows * n_rows - 1; + } } } } + /** + * Generic evaluator framework that valuates the given shape data in general + * dimensions using the tensor product form. Depending on the particular + * layout in the matrix entries, this corresponds to a usual matrix-matrix + * product or a matrix-matrix product including some symmetries. The actual + * work is implemented by functions of type apply_matrix_vector_product + * working on a single dimension, controlled by suitable strides, using the + * kernel specified via variant. + * + * @tparam variant Variant of evaluation used for creating template + * specializations + * @tparam dim Dimension of the function + * @tparam n_rows Number of rows in the transformation matrix, which corresponds + * to the number of 1d shape functions in the usual tensor + * contraction setting + * @tparam n_columns Number of columns in the transformation matrix, which + * corresponds to the number of 1d shape functions in the + * usual tensor contraction setting + * @tparam Number Abstract number type for input and output arrays + * @tparam Number2 Abstract number type for coefficient arrays (defaults to + * same type as the input/output arrays); must implement + * operator* with Number to be valid + * @tparam normal_dir Indicates the direction of the continuous component for the + * Raviart-Thomas space in terms of the normal onto the + * face, e.g 0 if the is in x-direction, 1 if in y-direction, and 2 if in + * z-direction. + */ + template + struct EvaluatorTensorProductAnisotropic + {}; + + + /** * Internal evaluator for shape function in 2d and 3d using the * tensor product form of the anisotropic basis functions of the