From 264dfb0277687e1136c9007bb5f03091bcfa6ead Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 30 Aug 2023 13:59:15 +0200 Subject: [PATCH] Simplify normal evaluation for Raviart-Thomas --- .../deal.II/matrix_free/evaluation_kernels.h | 262 +++++++++--------- .../matrix_free/tensor_product_kernels.h | 230 +-------------- 2 files changed, 130 insertions(+), 362 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 1401a04180..9e7dbc1a75 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3586,11 +3586,134 @@ namespace internal if (face_direction == face_no / 2 && !increase_max_der) { - interpolate_generic_raviart_thomas_apply_face( - shape_info, face_no, input, output); + constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction); + constexpr int stride0 = Utilities::pow(fe_degree, face_direction); + constexpr int stride2 = fe_degree * (fe_degree + 1); + + const int degree = + fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree; + const int n_rows_n = degree + 1; + const int n_rows_t = degree; + + std::array strides{{1, 1, 1}}; + if (face_direction > 0) + { + strides[0] = + n_rows_n * Utilities::pow(n_rows_t, face_direction - 1); + strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1); + strides[2] = Utilities::pow(n_rows_t, face_direction); + } + const dealii::ndarray dofs_per_direction{ + {{{n_rows_n, n_rows_t, n_rows_t}}, + {{n_rows_t, n_rows_n, n_rows_t}}, + {{n_rows_t, n_rows_t, n_rows_n}}}}; + + std::array steps, n_blocks; + + if constexpr (face_direction == 0) + steps = {{degree + (face_direction == 0), 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 * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}}; + else if constexpr (face_direction == 2) + steps = {{1, 0}}; + + n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim]; + n_blocks[1] = + dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1; + + interpolate_to_face< + (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0), + ((face_direction < 2) ? stride1 : stride2), + do_evaluate, + add_into_output, + max_derivative>(shape_info.data[face_direction != 0] + .shape_data_on_face[face_no % 2] + .begin(), + n_blocks, + steps, + input, + output, + degree + (face_direction == 0), + strides[0]); + + if (do_evaluate) + { + input += n_rows_n * Utilities::pow(n_rows_t, dim - 1); + output += 3 * n_blocks[0] * n_blocks[1]; + } + else + { + output += n_rows_n * Utilities::pow(n_rows_t, dim - 1); + input += 3 * n_blocks[0] * n_blocks[1]; + } + + // must only change steps only for face direction 0 + if constexpr (face_direction == 0) + steps = {{degree, 0}}; + + n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim]; + n_blocks[1] = + dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1; + + interpolate_to_face< + (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0), + ((face_direction < 2) ? stride0 : stride2), + do_evaluate, + add_into_output, + max_derivative>(shape_info.data[face_direction != 1] + .shape_data_on_face[face_no % 2] + .begin(), + n_blocks, + steps, + input, + output, + degree + (face_direction == 1), + strides[1]); + + if constexpr (dim > 2) + { + if (do_evaluate) + { + input += n_rows_n * Utilities::pow(n_rows_t, dim - 1); + output += 3 * n_blocks[0] * n_blocks[1]; + } + else + { + output += n_rows_n * Utilities::pow(n_rows_t, dim - 1); + input += 3 * n_blocks[0] * n_blocks[1]; + } + + if constexpr (face_direction == 0) + steps = {{degree, 0}}; + else if constexpr (face_direction == 1) + // in 3d, the coordinate system is zx, not xz -> switch indices + steps = { + {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}}; + else if constexpr (face_direction == 2) + steps = {{1, 0}}; + + n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim]; + n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim]; + + interpolate_to_face< + (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0), + stride0, + do_evaluate, + add_into_output, + max_derivative>(shape_info.data[face_direction != 2] + .shape_data_on_face[face_no % 2] + .begin(), + n_blocks, + steps, + input, + output, + degree + (face_direction == 2), + strides[2]); + } } else if (face_direction == face_no / 2) { @@ -3623,135 +3746,6 @@ namespace internal } } } - - /* Help function for interpolate_generic_raviart_thomas */ - template - static inline void - interpolate_generic_raviart_thomas_apply_face( - const MatrixFreeFunctions::ShapeInfo &shape_info, - const unsigned int face_no, - const Number *input, - Number *output) - { - // These types are evaluators in either normal or tangential direction - // depending on the face direction, with different normal directions for - // the different components. - using Evalf0 = typename std::conditional< - face_direction == 0, - EvaluatorTensorProductAnisotropic, - EvaluatorTensorProductAnisotropic>::type; - using Evalf1 = typename std::conditional< - face_direction == 1, - EvaluatorTensorProductAnisotropic, - EvaluatorTensorProductAnisotropic>::type; - using Evalf2 = typename std::conditional< - face_direction == 2, - EvaluatorTensorProductAnisotropic, - EvaluatorTensorProductAnisotropic>::type; - - Evalf0 evalf0 = - create_evaluator_tensor_product((face_direction == 0) ? - shape_info.data[0] : - shape_info.data[1], - face_no); - Evalf1 evalf1 = - create_evaluator_tensor_product((face_direction == 1) ? - shape_info.data[0] : - shape_info.data[1], - face_no); - Evalf2 evalf2 = - create_evaluator_tensor_product((face_direction == 2) ? - shape_info.data[0] : - shape_info.data[1], - face_no); - - const unsigned int dofs_per_component_on_cell = - shape_info.dofs_per_component_on_cell; - const unsigned int dofs_per_component_on_face = - 3 * shape_info.dofs_per_component_on_face; - - // NOTE! dofs_per_component_on_face is in the tangent direction, - // i.e (fe.degree+1)*fe.degree. Normal faces are only - // fe.degree*fe.degree - const unsigned int in_stride = - do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face; - const unsigned int out_stride = - do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell; - - const unsigned int in_stride_after_normal = - do_evaluate ? - dofs_per_component_on_cell : - dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2); - const unsigned int out_stride_after_normal = - do_evaluate ? - dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) : - dofs_per_component_on_cell; - - evalf0.template apply_face(input, output); - // stride to next component - input += (face_direction == 0) ? in_stride_after_normal : in_stride; - output += (face_direction == 0) ? out_stride_after_normal : out_stride; - - evalf1.template apply_face(input, output); - - if (dim == 3) - { - // stride to next component - input += (face_direction == 1) ? in_stride_after_normal : in_stride; - output += - (face_direction == 1) ? out_stride_after_normal : out_stride; - - evalf2.template apply_face(input, output); - } - } }; diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index fc9c842945..c9717aab71 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -1928,20 +1928,14 @@ namespace internal const Number *in, Number *out); - 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 - EvaluatorTensorProductAnisotropic< - evaluate_raviart_thomas, - dim, - n_rows, - n_columns, - normal_dir, - Number, - Number2>::apply_face(const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT out) const - { - Assert(dim > 1 && dim < 4, ExcMessage("Only dim=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.")); - - // Determine the number of blocks depending on the face and normaldirection, - // as well as dimension. - constexpr int n_blocks1 = (face_direction == normal_dir) ? (n_rows - 1) : - ((face_direction == 0 && normal_dir == 2) || - (face_direction == 1 && normal_dir == 2) || - (face_direction == 2 && normal_dir == 1)) ? - n_rows : - (n_rows + 1); - constexpr int n_blocks2 = (dim == 2) ? - 1 : - ((face_direction == normal_dir) ? - (n_rows - 1) : - (((face_direction == 0 && normal_dir == 1) || - (face_direction == 1 && normal_dir == 0) || - (face_direction == 2 && normal_dir == 0)) ? - n_rows : - (n_rows + 1))); - - AssertIndexRange(face_direction, dim); - - constexpr int in_stride = - (face_direction == normal_dir) ? - Utilities::pow(n_rows - 1, face_direction) : - ((face_direction == 0) ? - 1 : - ((face_direction == 2) ? - n_rows * (n_rows + 1) : - ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) : - n_rows))); - constexpr int out_stride = n_blocks1 * n_blocks2; - - 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 (normal_dir == 0) - { - if (contract_onto_face) - out += n_rows - 1; - else - in += n_rows - 1; - } - if (normal_dir == 1) - { - if (contract_onto_face) - out += n_rows - 2; - else - in += n_rows - 2; - } - if (normal_dir == 2) - { - if (contract_onto_face) - out += n_rows; - else - in += n_rows; - } - } - break; - - case 2: - ++in; - ++out; - break; - - default: - Assert(false, ExcNotImplemented()); - } - } - if (face_direction == 1 && dim == 3) - { - // adjust for local coordinate system zx - if (contract_onto_face) - { - if (normal_dir == 0) - { - in += (n_rows + 1) * (n_rows - 1); - out -= n_rows * (n_rows + 1) - 1; - } - if (normal_dir == 1) - { - in += (n_rows - 1) * (n_rows - 1); - out -= (n_rows - 1) * (n_rows - 1) - 1; - } - if (normal_dir == 2) - { - in += (n_rows - 1) * (n_rows); - out -= (n_rows) * (n_rows + 1) - 1; - } - } - else - { - if (normal_dir == 0) - { - out += (n_rows + 1) * (n_rows - 1); - in -= n_rows * (n_rows + 1) - 1; - } - if (normal_dir == 1) - { - out += (n_rows - 1) * (n_rows - 1); - in -= (n_rows - 1) * (n_rows - 1) - 1; - } - if (normal_dir == 2) - { - out += (n_rows - 1) * (n_rows); - in -= (n_rows) * (n_rows + 1) - 1; - } - } - } - } - } - /** -- 2.39.5