From 847dd0e99ffbd631aa7737ef97b342e93dfb7de5 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 30 Aug 2023 10:38:57 +0200 Subject: [PATCH] Delete obsolete code --- .../matrix_free/tensor_product_kernels.h | 417 +++++++----------- 1 file changed, 157 insertions(+), 260 deletions(-) diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 11aa7fefc7..fc9c842945 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -1130,13 +1130,8 @@ namespace internal } 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, + 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, @@ -1268,43 +1263,6 @@ namespace internal 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; @@ -1380,174 +1338,6 @@ namespace internal - template - inline std::enable_if_t - interpolate_to_face(const Number2 *shape_values, - const std::array &n_blocks, - const std::array &steps, - const Number *input, - Number *DEAL_II_RESTRICT output, - const int n_rows_runtime = 0, - const int stride_runtime = 1) - { - const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime; - const int stride = n_rows_template > 0 ? stride_template : stride_runtime; - - Number *output1 = output + n_blocks[0] * n_blocks[1]; - Number *output2 = output1 + n_blocks[0] * n_blocks[1]; - for (int i2 = 0; i2 < n_blocks[1]; ++i2) - { - for (int i1 = 0; i1 < n_blocks[0]; ++i1) - { - Number res0 = shape_values[0] * input[0]; - Number res1, res2; - if (max_derivative > 0) - res1 = shape_values[n_rows] * input[0]; - if (max_derivative > 1) - res2 = shape_values[2 * n_rows] * input[0]; - for (int ind = 1; ind < n_rows; ++ind) - { - res0 += shape_values[ind] * input[stride * ind]; - if (max_derivative > 0) - res1 += shape_values[ind + n_rows] * input[stride * ind]; - if (max_derivative > 1) - res2 += shape_values[ind + 2 * n_rows] * input[stride * ind]; - } - if (add) - { - output[i1] += res0; - if (max_derivative > 0) - output1[i1] += res1; - if (max_derivative > 1) - output2[i2] += res2; - } - else - { - output[i1] = res0; - if (max_derivative > 0) - output1[i1] = res1; - if (max_derivative > 1) - output2[i1] = res2; - } - input += steps[0]; - } - output += n_blocks[0]; - if (max_derivative > 0) - output1 += n_blocks[0]; - if (max_derivative > 1) - output2 += n_blocks[0]; - input += steps[1]; - } - } - - - - template - inline std::enable_if_t - interpolate_to_face(const Number2 *shape_values, - const std::array &n_blocks, - const std::array &steps, - const Number *input, - Number *DEAL_II_RESTRICT output, - const int n_rows_runtime = 0, - const int stride_runtime = 1) - { - const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime; - const int stride = n_rows_template > 0 ? stride_template : stride_runtime; - - const Number *input1 = input + n_blocks[0] * n_blocks[1]; - const Number *input2 = input1 + n_blocks[0] * n_blocks[1]; - for (int i2 = 0; i2 < n_blocks[1]; ++i2) - { - for (int i1 = 0; i1 < n_blocks[0]; ++i1) - { - const Number in = input[i1]; - Number in1, in2; - if (max_derivative > 0) - in1 = input1[i1]; - if (max_derivative > 1) - in2 = input2[i1]; - for (int col = 0; col < n_rows; ++col) - { - Number result = - add ? (output[col * stride] + shape_values[col] * in) : - (shape_values[col] * in); - if (max_derivative > 0) - result += shape_values[col + n_rows] * in1; - if (max_derivative > 1) - result += shape_values[col + 2 * n_rows] * in2; - - output[col * stride] = result; - } - output += steps[0]; - } - input += n_blocks[0]; - if (max_derivative > 0) - input1 += n_blocks[0]; - if (max_derivative > 1) - input2 += n_blocks[0]; - output += steps[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 stride = Utilities::pow(n_rows, face_direction); - std::array steps; - if constexpr (face_direction == 0) - steps = {{n_rows, 0}}; - else if constexpr (face_direction == 1 && dim == 2) - steps = {{1, 0}}; - else if constexpr (face_direction == 1) - // in 3d, the coordinate system is zx, not xz -> switch indices - steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}}; - else if constexpr (face_direction == 2) - steps = {{1, 0}}; - - interpolate_to_face( - this->shape_values, - {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}, - steps, - in, - out); - } - - - /** * Internal evaluator for shape function using the tensor product form * of the basis functions. The same as the other templated class but @@ -1612,13 +1402,8 @@ namespace internal } 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, + 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, @@ -1721,14 +1506,6 @@ namespace internal 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; @@ -1812,45 +1589,165 @@ namespace internal - template - template + interpolate_to_face(const Number2 *shape_values, + const std::array &n_blocks, + const std::array &steps, + const Number *input, + Number *DEAL_II_RESTRICT output, + const int n_rows_runtime = 0, + const int stride_runtime = 1) + { + const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime; + const int stride = n_rows_template > 0 ? stride_template : stride_runtime; + + Number *output1 = output + n_blocks[0] * n_blocks[1]; + Number *output2 = output1 + n_blocks[0] * n_blocks[1]; + for (int i2 = 0; i2 < n_blocks[1]; ++i2) + { + for (int i1 = 0; i1 < n_blocks[0]; ++i1) + { + Number res0 = shape_values[0] * input[0]; + Number res1, res2; + if (max_derivative > 0) + res1 = shape_values[n_rows] * input[0]; + if (max_derivative > 1) + res2 = shape_values[2 * n_rows] * input[0]; + for (int ind = 1; ind < n_rows; ++ind) + { + res0 += shape_values[ind] * input[stride * ind]; + if (max_derivative > 0) + res1 += shape_values[ind + n_rows] * input[stride * ind]; + if (max_derivative > 1) + res2 += shape_values[ind + 2 * n_rows] * input[stride * ind]; + } + if (add) + { + output[i1] += res0; + if (max_derivative > 0) + output1[i1] += res1; + if (max_derivative > 1) + output2[i2] += res2; + } + else + { + output[i1] = res0; + if (max_derivative > 0) + output1[i1] = res1; + if (max_derivative > 1) + output2[i1] = res2; + } + input += steps[0]; + } + output += n_blocks[0]; + if (max_derivative > 0) + output1 += n_blocks[0]; + if (max_derivative > 1) + output2 += n_blocks[0]; + input += steps[1]; + } + } + + + + template - inline void - EvaluatorTensorProduct::apply_face( - const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT out) const + int max_derivative, + typename Number, + typename Number2> + inline std::enable_if_t + interpolate_to_face(const Number2 *shape_values, + const std::array &n_blocks, + const std::array &steps, + const Number *input, + Number *DEAL_II_RESTRICT output, + const int n_rows_runtime = 0, + const int stride_runtime = 1) { - 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 stride = Utilities::pow(n_rows, face_direction); - const int n_rows = this->n_rows; - std::array steps; - if constexpr (face_direction == 0) - steps = {{n_rows, 0}}; - else if constexpr (face_direction == 1 && dim == 2) - steps = {{1, 0}}; - else if constexpr (face_direction == 1) - // in 3d, the coordinate system is zx, not xz -> switch indices - steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}}; - else if constexpr (face_direction == 2) - steps = {{1, 0}}; - - interpolate_to_face<0, 0, contract_to_face, add, max_derivative>( - this->shape_values, - {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}, - steps, - in, - out, - n_rows, - stride); + const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime; + const int stride = n_rows_template > 0 ? stride_template : stride_runtime; + + const Number *input1 = input + n_blocks[0] * n_blocks[1]; + const Number *input2 = input1 + n_blocks[0] * n_blocks[1]; + for (int i2 = 0; i2 < n_blocks[1]; ++i2) + { + for (int i1 = 0; i1 < n_blocks[0]; ++i1) + { + const Number in = input[i1]; + Number in1, in2; + if (max_derivative > 0) + in1 = input1[i1]; + if (max_derivative > 1) + in2 = input2[i1]; + for (int col = 0; col < n_rows; ++col) + { + Number result = + add ? (output[col * stride] + shape_values[col] * in) : + (shape_values[col] * in); + if (max_derivative > 0) + result += shape_values[col + n_rows] * in1; + if (max_derivative > 1) + result += shape_values[col + 2 * n_rows] * in2; + + output[col * stride] = result; + } + output += steps[0]; + } + input += n_blocks[0]; + if (max_derivative > 0) + input1 += n_blocks[0]; + if (max_derivative > 1) + input2 += n_blocks[0]; + output += steps[1]; + } } -- 2.39.5