From 5a99157295458636a80ee024fce03b20bf834d8c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 6 Jun 2022 00:19:24 +0200 Subject: [PATCH] Remove lex_faces from matrix-free kernels --- .../deal.II/matrix_free/evaluation_kernels.h | 20 ++- .../matrix_free/tensor_product_kernels.h | 118 ++++++------------ 2 files changed, 48 insertions(+), 90 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index ff05bbc5a0..bf83df4ea8 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3536,7 +3536,7 @@ namespace internal }; - template + template struct FEFaceNormalEvaluationImpl { template @@ -3637,20 +3637,17 @@ namespace internal evalf.template apply_face(input, output); + 2>(input, output); else if (flag & EvaluationFlags::gradients) evalf.template apply_face(input, output); + 1>(input, output); else evalf.template apply_face(input, output); + 0>(input, output); input += in_stride; output += out_stride; } @@ -3841,8 +3838,7 @@ namespace internal evalf0.template apply_face(input, output); + max_derivative>(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; @@ -3850,8 +3846,7 @@ namespace internal evalf1.template apply_face(input, output); + max_derivative>(input, output); if (dim == 3) { @@ -3863,8 +3858,7 @@ namespace internal evalf2.template apply_face(input, output); + max_derivative>(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 50f3166800..7a21fcd185 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -334,11 +334,6 @@ namespace internal * derivatives, 2 second derivates. Note that all the * derivatives access the data in @p shape_values passed to * the constructor of the class - * @tparam lex_faces Sets how the evaluation points on the faces should be - * sorted: lexicographically or right-hand-system number - * (special treatment of orientation 1 in 3D). Per default - * right-hand-system number is enabled, which is only - * working for dimensions up to 3. * * @param in address of the input data vector * @param out address of the output data vector @@ -346,8 +341,7 @@ namespace internal template + int max_derivative> void apply_face(const Number *DEAL_II_RESTRICT in, Number *DEAL_II_RESTRICT out) const; @@ -447,8 +441,7 @@ namespace internal template + int max_derivative> inline void EvaluatorTensorProduct 0 && (lex_faces || dim < 4), - ExcMessage("Only dim=1,2,3 supported")); + 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 = - lex_faces ? dealii::Utilities::pow(n_rows, face_direction) : - (dim > 1 ? n_rows : 1); - constexpr int n_blocks2 = - lex_faces ? dealii::Utilities::pow( - n_rows, std::max(dim - face_direction - 1, 0)) : - (dim > 2 ? n_rows : 1); + constexpr int n_blocks1 = (dim > 1 ? n_rows : 1); + constexpr int n_blocks2 = (dim > 2 ? n_rows : 1); AssertIndexRange(face_direction, dim); constexpr int stride = Utilities::pow(n_rows, face_direction); @@ -534,55 +521,41 @@ namespace internal } } - if (lex_faces) + // 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) { - ++out; - ++in; + 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()); } - else - // 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()); - } - } - if (lex_faces) - { - if (contract_onto_face) - in += (dealii::Utilities::pow(n_rows, face_direction + 1) - - n_blocks1); - else - out += (dealii::Utilities::pow(n_rows, face_direction + 1) - - n_blocks1); } - else if (face_direction == 1 && dim == 3) + + // adjust for local coordinate system zx + if (face_direction == 1 && dim == 3) { - // adjust for local coordinate system zx if (contract_onto_face) { in += n_rows * (n_rows - 1); @@ -734,8 +707,7 @@ namespace internal template + int max_derivative> void apply_face(const Number *DEAL_II_RESTRICT in, Number *DEAL_II_RESTRICT out) const; @@ -899,15 +871,12 @@ namespace internal template + int max_derivative> inline void EvaluatorTensorProduct:: apply_face(const Number *DEAL_II_RESTRICT in, Number *DEAL_II_RESTRICT out) const { - static_assert(lex_faces == false, "Not implemented yet."); - Assert(shape_values != nullptr, ExcMessage( "The given array shape_data must not be the null pointer!")); @@ -2505,8 +2474,7 @@ namespace internal template + int max_derivative> void apply_face(const Number *DEAL_II_RESTRICT in, Number *DEAL_II_RESTRICT out) const; @@ -2618,8 +2586,7 @@ namespace internal template + int max_derivative> inline void EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, @@ -2631,12 +2598,9 @@ namespace internal Number2>::apply_face(const Number *DEAL_II_RESTRICT in, Number *DEAL_II_RESTRICT out) const { - Assert(dim > 1 && (lex_faces || dim < 4), - ExcMessage("Only dim=2,3 supported")); + 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"); - static_assert(!lex_faces, - "lex_faces = True is not implemented for Raviart-Thomas"); Assert(shape_values != nullptr, ExcMessage( "The given array shape_values must not be the null pointer.")); -- 2.39.5