From 6f8c74d15dbfa36a9857ecfd05fa4fc7201c97bf Mon Sep 17 00:00:00 2001 From: Niklas Wik Date: Mon, 9 May 2022 13:01:06 +0200 Subject: [PATCH] Move template normal_dir from function to struct Creating new EvaluationTensorProductAnisotropic --- .../deal.II/matrix_free/evaluation_kernels.h | 947 +++++++++--------- .../matrix_free/tensor_product_kernels.h | 198 ++-- 2 files changed, 608 insertions(+), 537 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index fe15e2a34f..d5d8e57a4d 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -198,19 +198,6 @@ namespace internal n_q_points_1d, Number> { - using EvalNormal = - EvaluatorTensorProduct; - - using EvalTangent = - EvaluatorTensorProduct; template static void evaluate_or_integrate( @@ -223,22 +210,11 @@ namespace internal template static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::ShapeInfo &shape_info) + const MatrixFreeFunctions::UnivariateShapeData &shape_data) { - // Shape_data is of size 2 with normal direction in [0] and tangential - // direction in [1]. See shape_info.template.h - const MatrixFreeFunctions::UnivariateShapeData *shape_data; - if (std::is_same::value) - { - shape_data = &shape_info.data[0]; - } - else - { - shape_data = &shape_info.data[1]; - } - return EvalType(shape_data->shape_values, - shape_data->shape_gradients, - shape_data->shape_hessians); + return EvalType(shape_data.shape_values, + shape_data.shape_gradients, + shape_data.shape_hessians); } template @@ -1061,6 +1037,22 @@ namespace internal std::integral_constant) { (void)add_into_values_array; + + using EvalNormal = + EvaluatorTensorProductAnisotropic; + + using EvalTangent = + EvaluatorTensorProductAnisotropic; using Eval0 = typename std::conditional::type; using Eval1 = @@ -1068,24 +1060,25 @@ namespace internal using Eval2 = typename std::conditional::type; - Eval0 eval0 = - create_evaluator_tensor_product(fe_eval.get_shape_info()); - Eval1 eval1 = - create_evaluator_tensor_product(fe_eval.get_shape_info()); - Eval2 eval2 = - create_evaluator_tensor_product(fe_eval.get_shape_info()); + const MatrixFreeFunctions::ShapeInfo &shape_info = + fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( + ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1])); + Eval1 eval1 = create_evaluator_tensor_product( + ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1])); + Eval2 eval2 = create_evaluator_tensor_product( + ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1])); Number *temp1 = fe_eval.get_scratch_data().begin(); Number *temp2; - temp2 = temp1 + std::max(Utilities::fixed_power( - fe_eval.get_shape_info().data[0].fe_degree + 1), - Utilities::fixed_power( - fe_eval.get_shape_info().data[0].n_q_points_1d)); + temp2 = + temp1 + + std::max(Utilities::fixed_power(shape_info.data[0].fe_degree + 1), + Utilities::fixed_power(shape_info.data[0].n_q_points_1d)); - const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points; - const std::size_t dofs_per_comp = - fe_eval.get_shape_info().dofs_per_component_on_cell; + const std::size_t n_q_points = shape_info.n_q_points; + const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell; // Initial shift depending on component (normal_dir) Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir; @@ -1101,10 +1094,8 @@ namespace internal case 2: if (evaluation_flag & EvaluationFlags::gradients) { - eval0.template gradients<0, true, false, normal_dir>(values_dofs, - temp1); - eval1.template values<1, true, false, normal_dir>(temp1, - gradients_quad); + eval0.template gradients<0, true, false>(values_dofs, temp1); + eval1.template values<1, true, false>(temp1, gradients_quad); } if (evaluation_flag & EvaluationFlags::hessians) { @@ -1113,44 +1104,39 @@ namespace internal AssertThrow(false, ExcNotImplemented()); // grad xy if (!(evaluation_flag & EvaluationFlags::gradients)) - eval0.template gradients<0, true, false, normal_dir>( - values_dofs, temp1); - eval1.template gradients<1, true, false, normal_dir>( - temp1, hessians_quad + 2 * n_q_points); + eval0.template gradients<0, true, false>(values_dofs, temp1); + eval1.template gradients<1, true, false>(temp1, + hessians_quad + + 2 * n_q_points); // grad xx - eval0.template hessians<0, true, false, normal_dir>(values_dofs, - temp1); - eval1.template values<1, true, false, normal_dir>(temp1, - hessians_quad); + eval0.template hessians<0, true, false>(values_dofs, temp1); + eval1.template values<1, true, false>(temp1, hessians_quad); } // grad y - eval0.template values<0, true, false, normal_dir>(values_dofs, temp1); + eval0.template values<0, true, false>(values_dofs, temp1); if (evaluation_flag & EvaluationFlags::gradients) - eval1.template gradients<1, true, false, normal_dir>( - temp1, gradients_quad + n_q_points); + eval1.template gradients<1, true, false>(temp1, + gradients_quad + + n_q_points); // grad yy if (evaluation_flag & EvaluationFlags::hessians) - eval1.template hessians<1, true, false, normal_dir>(temp1, - hessians_quad + - n_q_points); + eval1.template hessians<1, true, false>(temp1, + hessians_quad + n_q_points); // val: can use values applied in x if (evaluation_flag & EvaluationFlags::values) - eval1.template values<1, true, false, normal_dir>(temp1, - values_quad); + eval1.template values<1, true, false>(temp1, values_quad); break; case 3: if (evaluation_flag & EvaluationFlags::gradients) { // grad x - eval0.template gradients<0, true, false, normal_dir>(values_dofs, - temp1); - eval1.template values<1, true, false, normal_dir>(temp1, temp2); - eval2.template values<2, true, false, normal_dir>(temp2, - gradients_quad); + eval0.template gradients<0, true, false>(values_dofs, temp1); + eval1.template values<1, true, false>(temp1, temp2); + eval2.template values<2, true, false>(temp2, gradients_quad); } if (evaluation_flag & EvaluationFlags::hessians) @@ -1161,73 +1147,69 @@ namespace internal // grad xz if (!(evaluation_flag & EvaluationFlags::gradients)) { - eval0.template gradients<0, true, false, normal_dir>( - values_dofs, temp1); - eval1.template values<1, true, false, normal_dir>(temp1, - temp2); + eval0.template gradients<0, true, false>(values_dofs, temp1); + eval1.template values<1, true, false>(temp1, temp2); } - eval2.template gradients<2, true, false, normal_dir>( - temp2, hessians_quad + 4 * n_q_points); + eval2.template gradients<2, true, false>(temp2, + hessians_quad + + 4 * n_q_points); // grad xy - eval1.template gradients<1, true, false, normal_dir>(temp1, - temp2); - eval2.template values<2, true, false, normal_dir>( - temp2, hessians_quad + 3 * n_q_points); + eval1.template gradients<1, true, false>(temp1, temp2); + eval2.template values<2, true, false>(temp2, + hessians_quad + + 3 * n_q_points); // grad xx - eval0.template hessians<0, true, false, normal_dir>(values_dofs, - temp1); - eval1.template values<1, true, false, normal_dir>(temp1, temp2); - eval2.template values<2, true, false, normal_dir>(temp2, - hessians_quad); + eval0.template hessians<0, true, false>(values_dofs, temp1); + eval1.template values<1, true, false>(temp1, temp2); + eval2.template values<2, true, false>(temp2, hessians_quad); } // grad y - eval0.template values<0, true, false, normal_dir>(values_dofs, temp1); + eval0.template values<0, true, false>(values_dofs, temp1); if (evaluation_flag & EvaluationFlags::gradients) { - eval1.template gradients<1, true, false, normal_dir>(temp1, - temp2); - eval2.template values<2, true, false, normal_dir>(temp2, - gradients_quad + - n_q_points); + eval1.template gradients<1, true, false>(temp1, temp2); + eval2.template values<2, true, false>(temp2, + gradients_quad + + n_q_points); } if (evaluation_flag & EvaluationFlags::hessians) { // grad yz if (!(evaluation_flag & EvaluationFlags::gradients)) - eval1.template gradients<1, true, false, normal_dir>(temp1, - temp2); - eval2.template gradients<2, true, false, normal_dir>( - temp2, hessians_quad + 5 * n_q_points); + eval1.template gradients<1, true, false>(temp1, temp2); + eval2.template gradients<2, true, false>(temp2, + hessians_quad + + 5 * n_q_points); // grad yy - eval1.template hessians<1, true, false, normal_dir>(temp1, temp2); - eval2.template values<2, true, false, normal_dir>(temp2, - hessians_quad + - n_q_points); + eval1.template hessians<1, true, false>(temp1, temp2); + eval2.template values<2, true, false>(temp2, + hessians_quad + n_q_points); } // grad z: can use the values applied in x direction stored in // temp1 - eval1.template values<1, true, false, normal_dir>(temp1, temp2); + eval1.template values<1, true, false>(temp1, temp2); if (evaluation_flag & EvaluationFlags::gradients) - eval2.template gradients<2, true, false, normal_dir>( - temp2, gradients_quad + 2 * n_q_points); + eval2.template gradients<2, true, false>(temp2, + gradients_quad + + 2 * n_q_points); // grad zz: can use the values applied in x and y direction stored // in temp2 if (evaluation_flag & EvaluationFlags::hessians) - eval2.template hessians<2, true, false, normal_dir>( - temp2, hessians_quad + 2 * n_q_points); + eval2.template hessians<2, true, false>(temp2, + hessians_quad + + 2 * n_q_points); // val: can use the values applied in x & y direction stored in // temp2 if (evaluation_flag & EvaluationFlags::values) - eval2.template values<2, true, false, normal_dir>(temp2, - values_quad); + eval2.template values<2, true, false>(temp2, values_quad); break; default: AssertThrow(false, ExcNotImplemented()); @@ -1249,6 +1231,21 @@ namespace internal const bool add_into_values_array, std::integral_constant) { + using EvalNormal = + EvaluatorTensorProductAnisotropic; + + using EvalTangent = + EvaluatorTensorProductAnisotropic; using Eval0 = typename std::conditional::type; using Eval1 = @@ -1256,24 +1253,25 @@ namespace internal using Eval2 = typename std::conditional::type; - Eval0 eval0 = - create_evaluator_tensor_product(fe_eval.get_shape_info()); - Eval1 eval1 = - create_evaluator_tensor_product(fe_eval.get_shape_info()); - Eval2 eval2 = - create_evaluator_tensor_product(fe_eval.get_shape_info()); + const MatrixFreeFunctions::ShapeInfo &shape_info = + fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( + ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1])); + Eval1 eval1 = create_evaluator_tensor_product( + ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1])); + Eval2 eval2 = create_evaluator_tensor_product( + ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1])); Number *temp1 = fe_eval.get_scratch_data().begin(); Number *temp2; - temp2 = temp1 + std::max(Utilities::fixed_power( - fe_eval.get_shape_info().data[0].fe_degree + 1), - Utilities::fixed_power( - fe_eval.get_shape_info().data[0].n_q_points_1d)); + temp2 = + temp1 + + std::max(Utilities::fixed_power(shape_info.data[0].fe_degree + 1), + Utilities::fixed_power(shape_info.data[0].n_q_points_1d)); - const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points; - const std::size_t dofs_per_comp = - fe_eval.get_shape_info().dofs_per_component_on_cell; + const std::size_t n_q_points = shape_info.n_q_points; + const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell; // Initial shift depending on component (normal_dir) Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir; @@ -1291,59 +1289,49 @@ namespace internal if ((evaluation_flag & EvaluationFlags::values) && !(evaluation_flag & EvaluationFlags::gradients)) { - eval1.template values<1, false, false, normal_dir>(values_quad, - temp1); + eval1.template values<1, false, false>(values_quad, temp1); if (add_into_values_array == false) - eval0.template values<0, false, false, normal_dir>(temp1, - values_dofs); + eval0.template values<0, false, false>(temp1, values_dofs); else - eval0.template values<0, false, true, normal_dir>(temp1, - values_dofs); + eval0.template values<0, false, true>(temp1, values_dofs); } if (evaluation_flag & EvaluationFlags::gradients) { - eval1.template gradients<1, false, false, normal_dir>( - gradients_quad + n_q_points, temp1); + eval1.template gradients<1, false, false>(gradients_quad + + n_q_points, + temp1); if ((evaluation_flag & EvaluationFlags::values)) - eval1.template values<1, false, true, normal_dir>(values_quad, - temp1); + eval1.template values<1, false, true>(values_quad, temp1); if (add_into_values_array == false) - eval0.template values<0, false, false, normal_dir>(temp1, - values_dofs); + eval0.template values<0, false, false>(temp1, values_dofs); else - eval0.template values<0, false, true, normal_dir>(temp1, - values_dofs); - eval1.template values<1, false, false, normal_dir>(gradients_quad, - temp1); - eval0.template gradients<0, false, true, normal_dir>(temp1, - values_dofs); + eval0.template values<0, false, true>(temp1, values_dofs); + eval1.template values<1, false, false>(gradients_quad, temp1); + eval0.template gradients<0, false, true>(temp1, values_dofs); } if (evaluation_flag & EvaluationFlags::hessians) { // grad xx - eval1.template values<1, false, false, normal_dir>(hessians_quad, - temp1); + eval1.template values<1, false, false>(hessians_quad, temp1); if ((evaluation_flag & EvaluationFlags::values) || (evaluation_flag & EvaluationFlags::gradients) || add_into_values_array == true) - eval0.template hessians<0, false, true, normal_dir>( - temp1, values_dofs); + eval0.template hessians<0, false, true>(temp1, values_dofs); else - eval0.template hessians<0, false, false, normal_dir>( - temp1, values_dofs); + eval0.template hessians<0, false, false>(temp1, values_dofs); // grad yy - eval1.template hessians<1, false, false, normal_dir>( - hessians_quad + n_q_points, temp1); - eval0.template values<0, false, true, normal_dir>(temp1, - values_dofs); + eval1.template hessians<1, false, false>(hessians_quad + + n_q_points, + temp1); + eval0.template values<0, false, true>(temp1, values_dofs); // grad xy - eval1.template gradients<1, false, false, normal_dir>( - hessians_quad + 2 * n_q_points, temp1); - eval0.template gradients<0, false, true, normal_dir>(temp1, - values_dofs); + eval1.template gradients<1, false, false>(hessians_quad + + 2 * n_q_points, + temp1); + eval0.template gradients<0, false, true>(temp1, values_dofs); } break; @@ -1351,94 +1339,79 @@ namespace internal if ((evaluation_flag & EvaluationFlags::values) && !(evaluation_flag & EvaluationFlags::gradients)) { - eval2.template values<2, false, false, normal_dir>(values_quad, - temp1); - eval1.template values<1, false, false, normal_dir>(temp1, temp2); + eval2.template values<2, false, false>(values_quad, temp1); + eval1.template values<1, false, false>(temp1, temp2); if (add_into_values_array == false) - eval0.template values<0, false, false, normal_dir>(temp2, - values_dofs); + eval0.template values<0, false, false>(temp2, values_dofs); else - eval0.template values<0, false, true, normal_dir>(temp2, - values_dofs); + eval0.template values<0, false, true>(temp2, values_dofs); } if (evaluation_flag & EvaluationFlags::gradients) { - eval2.template gradients<2, false, false, normal_dir>( - gradients_quad + 2 * n_q_points, temp1); + eval2.template gradients<2, false, false>(gradients_quad + + 2 * n_q_points, + temp1); if ((evaluation_flag & EvaluationFlags::values)) - eval2.template values<2, false, true, normal_dir>(values_quad, - temp1); - eval1.template values<1, false, false, normal_dir>(temp1, temp2); - eval2.template values<2, false, false, normal_dir>( - gradients_quad + n_q_points, temp1); - eval1.template gradients<1, false, true, normal_dir>(temp1, - temp2); + eval2.template values<2, false, true>(values_quad, temp1); + eval1.template values<1, false, false>(temp1, temp2); + eval2.template values<2, false, false>(gradients_quad + + n_q_points, + temp1); + eval1.template gradients<1, false, true>(temp1, temp2); if (add_into_values_array == false) - eval0.template values<0, false, false, normal_dir>(temp2, - values_dofs); + eval0.template values<0, false, false>(temp2, values_dofs); else - eval0.template values<0, false, true, normal_dir>(temp2, - values_dofs); - eval2.template values<2, false, false, normal_dir>(gradients_quad, - temp1); - eval1.template values<1, false, false, normal_dir>(temp1, temp2); - eval0.template gradients<0, false, true, normal_dir>(temp2, - values_dofs); + eval0.template values<0, false, true>(temp2, values_dofs); + eval2.template values<2, false, false>(gradients_quad, temp1); + eval1.template values<1, false, false>(temp1, temp2); + eval0.template gradients<0, false, true>(temp2, values_dofs); } if (evaluation_flag & EvaluationFlags::hessians) { // grad xx - eval2.template values<2, false, false, normal_dir>(hessians_quad, - temp1); - eval1.template values<1, false, false, normal_dir>(temp1, temp2); + eval2.template values<2, false, false>(hessians_quad, temp1); + eval1.template values<1, false, false>(temp1, temp2); if ((evaluation_flag & EvaluationFlags::values) || (evaluation_flag & EvaluationFlags::gradients) || add_into_values_array == true) - eval0.template hessians<0, false, true, normal_dir>( - temp2, values_dofs); + eval0.template hessians<0, false, true>(temp2, values_dofs); else - eval0.template hessians<0, false, false, normal_dir>( - temp2, values_dofs); + eval0.template hessians<0, false, false>(temp2, values_dofs); // grad yy - eval2.template values<2, false, false, normal_dir>(hessians_quad + - n_q_points, - temp1); - eval1.template hessians<1, false, false, normal_dir>(temp1, - temp2); - eval0.template values<0, false, true, normal_dir>(temp2, - values_dofs); + eval2.template values<2, false, false>(hessians_quad + n_q_points, + temp1); + eval1.template hessians<1, false, false>(temp1, temp2); + eval0.template values<0, false, true>(temp2, values_dofs); // grad zz - eval2.template hessians<2, false, false, normal_dir>( - hessians_quad + 2 * n_q_points, temp1); - eval1.template values<1, false, false, normal_dir>(temp1, temp2); - eval0.template values<0, false, true, normal_dir>(temp2, - values_dofs); + eval2.template hessians<2, false, false>(hessians_quad + + 2 * n_q_points, + temp1); + eval1.template values<1, false, false>(temp1, temp2); + eval0.template values<0, false, true>(temp2, values_dofs); // grad xy - eval2.template values<2, false, false, normal_dir>( - hessians_quad + 3 * n_q_points, temp1); - eval1.template gradients<1, false, false, normal_dir>(temp1, - temp2); - eval0.template gradients<0, false, true, normal_dir>(temp2, - values_dofs); + eval2.template values<2, false, false>(hessians_quad + + 3 * n_q_points, + temp1); + eval1.template gradients<1, false, false>(temp1, temp2); + eval0.template gradients<0, false, true>(temp2, values_dofs); // grad xz - eval2.template gradients<2, false, false, normal_dir>( - hessians_quad + 4 * n_q_points, temp1); - eval1.template values<1, false, false, normal_dir>(temp1, temp2); - eval0.template gradients<0, false, true, normal_dir>(temp2, - values_dofs); + eval2.template gradients<2, false, false>(hessians_quad + + 4 * n_q_points, + temp1); + eval1.template values<1, false, false>(temp1, temp2); + eval0.template gradients<0, false, true>(temp2, values_dofs); // grad yz - eval2.template gradients<2, false, false, normal_dir>( - hessians_quad + 5 * n_q_points, temp1); - eval1.template gradients<1, false, false, normal_dir>(temp1, - temp2); - eval0.template values<0, false, true, normal_dir>(temp2, - values_dofs); + eval2.template gradients<2, false, false>(hessians_quad + + 5 * n_q_points, + temp1); + eval1.template gradients<1, false, false>(temp1, temp2); + eval0.template values<0, false, true>(temp2, values_dofs); } break; @@ -3035,19 +3008,6 @@ namespace internal template struct FEFaceEvaluationImplRaviartThomas { - using EvalNormal = - EvaluatorTensorProduct; - using EvalTangent = - EvaluatorTensorProduct; - using EvalGeneral = EvaluatorTensorProduct static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::ShapeInfo &shape_info, - const unsigned int subface_index, - const unsigned int direction) + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int subface_index, + const unsigned int direction) { - const MatrixFreeFunctions::UnivariateShapeData &data = - (std::is_same::value) ? shape_info.data.front() : - shape_info.data.back(); if (subface_index >= GeometryInfo::max_children_per_cell) return EvalType(data.shape_values, data.shape_gradients, data.shape_hessians); else { - Assert(false, ExcNotImplemented()); - const unsigned int index = direction == 0 ? subface_index % 2 : subface_index / 2; return EvalType(data.values_within_subface[index], @@ -3091,33 +3046,48 @@ namespace internal { const unsigned int face_direction = face_no / 2; - evaluate_in_face_apply<0, EvalNormal, EvalTangent>( - values_dofs, - fe_eval, - scratch_data, - evaluation_flag, - face_direction, - subface_index, - std::integral_constant()); + // We first evaluate the anisotropic faces, i.e the faces where + // face_direction != component. Note that the call order here is not + // important, since the pointers are shifted accordingly within the + // function. However, this is the order in which the components will be in + // the quadrature points. Furthermore, the isotropic faces have no "normal + // direction" but we still pass in normal_dir = 2 since this is used for + // the pointers. + // ----------------------------------------------------------------------------------- + // | | Anisotropic faces | Isotropic faces| + // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords | + // | --------------------------------------------------------------------------------| + // | 0 | 1, y, 0 | - | 0, y | + // | 1 | 0, x, 0 | - | 1, x | + // | --------------------------------------------------------------------------------| + // | 0 | 1, yz, 0 | 2, yz, 1 | 0, yz | + // | 1 | 2, zx, 0 | 0, zx, 1 | 1, zx | + // | 2 | 0, xy, 0 | 1, xy, 1 | 2, xy | + // ----------------------------------------------------------------------------------- + evaluate_in_face_apply<0>(values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + face_direction, + subface_index, + std::integral_constant()); if (dim == 3) - evaluate_in_face_apply<1, EvalTangent, EvalNormal>( - values_dofs, - fe_eval, - scratch_data, - evaluation_flag, - face_direction, - subface_index, - std::integral_constant()); - - evaluate_in_face_apply<2, EvalGeneral, EvalGeneral>( - values_dofs, - fe_eval, - scratch_data, - evaluation_flag, - face_direction, - subface_index, - std::integral_constant()); + evaluate_in_face_apply<1>(values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + face_direction, + subface_index, + std::integral_constant()); + + evaluate_in_face_apply<2>(values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + face_direction, + subface_index, + std::integral_constant()); } /* @@ -3128,7 +3098,7 @@ namespace internal * integration path. These two functions can be fused together since all * offsets and pointers are the exact same. */ - template + template static inline void evaluate_in_face_apply( Number * values_dofs, @@ -3139,20 +3109,43 @@ namespace internal const unsigned int subface_index, std::integral_constant) { - // TODO. Need to make sure hanging nodes work, i.e call - // create_eval_tensor_product with correct direction. - Eval0 eval0 = - create_evaluator_tensor_product(fe_eval.get_shape_info(), - subface_index, - 0); - Eval1 eval1 = - create_evaluator_tensor_product(fe_eval.get_shape_info(), - subface_index, - 0); + using EvalNormal = + EvaluatorTensorProductAnisotropic; + using EvalTangent = + EvaluatorTensorProductAnisotropic; + + using TempEval0 = typename std:: + conditional::type; + using TempEval1 = typename std:: + conditional::type; + using Eval0 = typename std:: + conditional::type; + using Eval1 = typename std:: + conditional::type; + + const MatrixFreeFunctions::ShapeInfo &shape_info = + fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( + ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]), + subface_index, + 0); + Eval1 eval1 = create_evaluator_tensor_product( + ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]), + subface_index, + 1); constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); - const std::size_t n_dofs_tangent = - fe_eval.get_shape_info().dofs_per_component_on_face; + const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face; const std::size_t n_dofs_normal = n_dofs_tangent - Utilities::pow(fe_degree, dim - 2); const std::size_t dofs_stride = @@ -3188,14 +3181,11 @@ namespace internal switch (dim) { case 3: - eval0.template values<0, true, false, normal_dir>(values_dofs, - values_quad); - eval1.template values<1, true, false, normal_dir>(values_quad, - values_quad); + eval0.template values<0, true, false>(values_dofs, values_quad); + eval1.template values<1, true, false>(values_quad, values_quad); break; case 2: - eval0.template values<0, true, false, normal_dir>(values_dofs, - values_quad); + eval0.template values<0, true, false>(values_dofs, values_quad); break; default: Assert(false, ExcNotImplemented()); @@ -3207,37 +3197,39 @@ namespace internal { case 3: // grad x - eval0.template gradients<0, true, false, normal_dir>( - values_dofs, scratch_data); - eval1.template values<1, true, false, normal_dir>( - scratch_data, gradients_quad); + eval0.template gradients<0, true, false>(values_dofs, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + gradients_quad); // grad y - eval0.template values<0, true, false, normal_dir>(values_dofs, - scratch_data); - eval1.template gradients<1, true, false, normal_dir>( - scratch_data, gradients_quad + n_q_points); + eval0.template values<0, true, false>(values_dofs, + scratch_data); + eval1.template gradients<1, true, false>(scratch_data, + gradients_quad + + n_q_points); if (evaluation_flag & EvaluationFlags::values) - eval1.template values<1, true, false, normal_dir>( - scratch_data, values_quad); + eval1.template values<1, true, false>(scratch_data, + values_quad); // grad z - eval0.template values<0, true, false, normal_dir>(values_dofs + - dofs_stride, - scratch_data); - eval1.template values<1, true, false, normal_dir>( - scratch_data, gradients_quad + 2 * n_q_points); + eval0.template values<0, true, false>(values_dofs + dofs_stride, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + gradients_quad + + 2 * n_q_points); break; case 2: - eval0.template values<0, true, false, normal_dir>( - values_dofs + dofs_stride, gradients_quad + n_q_points); - eval0.template gradients<0, true, false, normal_dir>( - values_dofs, gradients_quad); + eval0.template values<0, true, false>(values_dofs + dofs_stride, + gradients_quad + + n_q_points); + eval0.template gradients<0, true, false>(values_dofs, + gradients_quad); if ((evaluation_flag & EvaluationFlags::values)) - eval0.template values<0, true, false, normal_dir>( - values_dofs, values_quad); + eval0.template values<0, true, false>(values_dofs, + values_quad); break; default: AssertThrow(false, ExcNotImplemented()); @@ -3250,52 +3242,58 @@ namespace internal { case 3: // grad xx - eval0.template hessians<0, true, false, normal_dir>( - values_dofs, scratch_data); - eval1.template values<1, true, false, normal_dir>( - scratch_data, hessians_quad); + eval0.template hessians<0, true, false>(values_dofs, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + hessians_quad); // grad yy - eval0.template values<0, true, false, normal_dir>(values_dofs, - scratch_data); - eval1.template hessians<1, true, false, normal_dir>( - scratch_data, hessians_quad + n_q_points); + eval0.template values<0, true, false>(values_dofs, + scratch_data); + eval1.template hessians<1, true, false>(scratch_data, + hessians_quad + + n_q_points); // grad zz - eval0.template values<0, true, false, normal_dir>( - values_dofs + 2 * dofs_stride, scratch_data); - eval1.template values<1, true, false, normal_dir>( - scratch_data, hessians_quad + 2 * n_q_points); + eval0.template values<0, true, false>(values_dofs + + 2 * dofs_stride, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + hessians_quad + + 2 * n_q_points); // grad xy - eval0.template gradients<0, true, false, normal_dir>( - values_dofs, scratch_data); - eval1.template gradients<1, true, false, normal_dir>( - scratch_data, hessians_quad + 3 * n_q_points); + eval0.template gradients<0, true, false>(values_dofs, + scratch_data); + eval1.template gradients<1, true, false>(scratch_data, + hessians_quad + + 3 * n_q_points); // grad xz - eval0.template gradients<0, true, false, normal_dir>( - values_dofs + dofs_stride, scratch_data); - eval1.template values<1, true, false, normal_dir>( - scratch_data, hessians_quad + 4 * n_q_points); + eval0.template gradients<0, true, false>(values_dofs + + dofs_stride, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + hessians_quad + + 4 * n_q_points); // grad yz - eval0.template values<0, true, false, normal_dir>(values_dofs + - dofs_stride, - scratch_data); - eval1.template gradients<1, true, false, normal_dir>( - scratch_data, hessians_quad + 5 * n_q_points); + eval0.template values<0, true, false>(values_dofs + dofs_stride, + scratch_data); + eval1.template gradients<1, true, false>(scratch_data, + hessians_quad + + 5 * n_q_points); break; case 2: // grad xx - eval0.template hessians<0, true, false, normal_dir>( - values_dofs, hessians_quad); + eval0.template hessians<0, true, false>(values_dofs, + hessians_quad); // grad yy - eval0.template values<0, true, false, normal_dir>( + eval0.template values<0, true, false>( values_dofs + 2 * dofs_stride, hessians_quad + n_q_points); // grad xy - eval0.template gradients<0, true, false, normal_dir>( + eval0.template gradients<0, true, false>( values_dofs + dofs_stride, hessians_quad + 2 * n_q_points); break; default: @@ -3304,7 +3302,7 @@ namespace internal } } - template + template static inline void evaluate_in_face_apply( Number * values_dofs, @@ -3315,20 +3313,43 @@ namespace internal const unsigned int subface_index, std::integral_constant) { - // TODO. Need to make sure hanging nodes work, i.e call - // create_eval_tensor_product with correct direction. - Eval0 eval0 = - create_evaluator_tensor_product(fe_eval.get_shape_info(), - subface_index, - 0); - Eval1 eval1 = - create_evaluator_tensor_product(fe_eval.get_shape_info(), - subface_index, - 0); + using EvalNormal = + EvaluatorTensorProductAnisotropic; + using EvalTangent = + EvaluatorTensorProductAnisotropic; + + using TempEval0 = typename std:: + conditional::type; + using TempEval1 = typename std:: + conditional::type; + using Eval0 = typename std:: + conditional::type; + using Eval1 = typename std:: + conditional::type; + + const MatrixFreeFunctions::ShapeInfo &shape_info = + fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( + ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]), + subface_index, + 0); + Eval1 eval1 = create_evaluator_tensor_product( + ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]), + subface_index, + 1); constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); - const std::size_t n_dofs_tangent = - fe_eval.get_shape_info().dofs_per_component_on_face; + const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face; const std::size_t n_dofs_normal = n_dofs_tangent - Utilities::pow(fe_degree, dim - 2); const std::size_t dofs_stride = @@ -3364,14 +3385,14 @@ namespace internal switch (dim) { case 3: - eval1.template values<1, false, false, normal_dir>(values_quad, - values_quad); - eval0.template values<0, false, false, normal_dir>(values_quad, - values_dofs); + eval1.template values<1, false, false>(values_quad, + values_quad); + eval0.template values<0, false, false>(values_quad, + values_dofs); break; case 2: - eval0.template values<0, false, false, normal_dir>(values_quad, - values_dofs); + eval0.template values<0, false, false>(values_quad, + values_dofs); break; default: Assert(false, ExcNotImplemented()); @@ -3383,42 +3404,45 @@ namespace internal { case 3: // grad z - eval1.template values<1, false, false, normal_dir>( - gradients_quad + 2 * n_q_points, - gradients_quad + 2 * n_q_points); - eval0.template values<0, false, false, normal_dir>( + eval1.template values<1, false, false>(gradients_quad + + 2 * n_q_points, + gradients_quad + + 2 * n_q_points); + eval0.template values<0, false, false>( gradients_quad + 2 * n_q_points, values_dofs + dofs_stride); if (evaluation_flag & EvaluationFlags::values) { - eval1.template values<1, false, false, normal_dir>( - values_quad, scratch_data); - eval1.template gradients<1, false, true, normal_dir>( - gradients_quad + n_q_points, scratch_data); + eval1.template values<1, false, false>(values_quad, + scratch_data); + eval1.template gradients<1, false, true>(gradients_quad + + n_q_points, + scratch_data); } else - eval1.template gradients<1, false, false, normal_dir>( - gradients_quad + n_q_points, scratch_data); + eval1.template gradients<1, false, false>(gradients_quad + + n_q_points, + scratch_data); // grad y - eval0.template values<0, false, false, normal_dir>(scratch_data, - values_dofs); + eval0.template values<0, false, false>(scratch_data, + values_dofs); // grad x - eval1.template values<1, false, false, normal_dir>( - gradients_quad, scratch_data); - eval0.template gradients<0, false, true, normal_dir>( - scratch_data, values_dofs); + eval1.template values<1, false, false>(gradients_quad, + scratch_data); + eval0.template gradients<0, false, true>(scratch_data, + values_dofs); break; case 2: - eval0.template values<0, false, false, normal_dir>( + eval0.template values<0, false, false>( gradients_quad + n_q_points, values_dofs + dofs_stride); - eval0.template gradients<0, false, false, normal_dir>( - gradients_quad, values_dofs); + eval0.template gradients<0, false, false>(gradients_quad, + values_dofs); if (evaluation_flag & EvaluationFlags::values) - eval0.template values<0, false, true, normal_dir>( - values_quad, values_dofs); + eval0.template values<0, false, true>(values_quad, + values_dofs); break; default: AssertThrow(false, ExcNotImplemented()); @@ -3431,70 +3455,79 @@ namespace internal { case 3: // grad xx - eval1.template values<1, false, false, normal_dir>( - hessians_quad, scratch_data); + eval1.template values<1, false, false>(hessians_quad, + scratch_data); if ((evaluation_flag & (EvaluationFlags::values | EvaluationFlags::gradients))) - eval0.template hessians<0, false, true, normal_dir>( - scratch_data, values_dofs); + eval0.template hessians<0, false, true>(scratch_data, + values_dofs); else - eval0.template hessians<0, false, false, normal_dir>( - scratch_data, values_dofs); + eval0.template hessians<0, false, false>(scratch_data, + values_dofs); // grad yy - eval1.template hessians<1, false, false, normal_dir>( - hessians_quad + n_q_points, scratch_data); - eval0.template values<0, false, true, normal_dir>(scratch_data, - values_dofs); + eval1.template hessians<1, false, false>(hessians_quad + + n_q_points, + scratch_data); + eval0.template values<0, false, true>(scratch_data, + values_dofs); // grad zz - eval1.template values<1, false, false, normal_dir>( - hessians_quad + 2 * n_q_points, scratch_data); - eval0.template values<0, false, false, normal_dir>( - scratch_data, values_dofs + 2 * dofs_stride); + eval1.template values<1, false, false>(hessians_quad + + 2 * n_q_points, + scratch_data); + eval0.template values<0, false, false>(scratch_data, + values_dofs + + 2 * dofs_stride); // grad xy - eval1.template gradients<1, false, false, normal_dir>( - hessians_quad + 3 * n_q_points, scratch_data); - eval0.template gradients<0, false, true, normal_dir>( - scratch_data, values_dofs); + eval1.template gradients<1, false, false>(hessians_quad + + 3 * n_q_points, + scratch_data); + eval0.template gradients<0, false, true>(scratch_data, + values_dofs); // grad xz - eval1.template values<1, false, false, normal_dir>( - hessians_quad + 4 * n_q_points, scratch_data); + eval1.template values<1, false, false>(hessians_quad + + 4 * n_q_points, + scratch_data); if ((evaluation_flag & EvaluationFlags::gradients)) - eval0.template gradients<0, false, true, normal_dir>( - scratch_data, values_dofs + dofs_stride); + eval0.template gradients<0, false, true>(scratch_data, + values_dofs + + dofs_stride); else - eval0.template gradients<0, false, false, normal_dir>( - scratch_data, values_dofs + dofs_stride); + eval0.template gradients<0, false, false>(scratch_data, + values_dofs + + dofs_stride); // grad yz - eval1.template gradients<1, false, false, normal_dir>( - hessians_quad + 5 * n_q_points, scratch_data); - eval0.template values<0, false, true, normal_dir>( - scratch_data, values_dofs + dofs_stride); + eval1.template gradients<1, false, false>(hessians_quad + + 5 * n_q_points, + scratch_data); + eval0.template values<0, false, true>(scratch_data, + values_dofs + + dofs_stride); break; case 2: // grad xx if (evaluation_flag & (EvaluationFlags::values | EvaluationFlags::gradients)) - eval0.template hessians<0, false, true, normal_dir>( - hessians_quad, values_dofs); + eval0.template hessians<0, false, true>(hessians_quad, + values_dofs); else - eval0.template hessians<0, false, false, normal_dir>( - hessians_quad, values_dofs); + eval0.template hessians<0, false, false>(hessians_quad, + values_dofs); // grad yy - eval0.template values<0, false, false, normal_dir>( + eval0.template values<0, false, false>( hessians_quad + n_q_points, values_dofs + 2 * dofs_stride); // grad xy if ((evaluation_flag & EvaluationFlags::gradients)) - eval0.template gradients<0, false, true, normal_dir>( + eval0.template gradients<0, false, true>( hessians_quad + 2 * n_q_points, values_dofs + dofs_stride); else - eval0.template gradients<0, false, false, normal_dir>( + eval0.template gradients<0, false, false>( hessians_quad + 2 * n_q_points, values_dofs + dofs_stride); break; default: @@ -3508,19 +3541,6 @@ namespace internal template struct FEFaceNormalEvaluationImpl { - using EvalNormal = - EvaluatorTensorProduct; - using EvalTangent = - EvaluatorTensorProduct; - template static void interpolate(const unsigned int n_components, @@ -3658,21 +3678,10 @@ namespace internal template static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::ShapeInfo &shape_info, - const unsigned int face_no) + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int face_no) { - // Shape_data is of size 2 with normal direction in [0] and tangential - // direction in [1]. See shape_info.template.h - const MatrixFreeFunctions::UnivariateShapeData *shape_data; - if (std::is_same::value) - { - shape_data = &shape_info.data[0]; - } - else - { - shape_data = &shape_info.data[1]; - } - return EvalType(shape_data->shape_data_on_face[face_no % 2], + return EvalType(data.shape_data_on_face[face_no % 2], AlignedVector(), AlignedVector()); } @@ -3747,19 +3756,67 @@ namespace internal const Number * input, Number * output) { - using Evalf0 = typename std:: - conditional::type; - using Evalf1 = typename std:: - conditional::type; - using Evalf2 = typename std:: - conditional::type; + // 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(shape_info, face_no); + create_evaluator_tensor_product((face_direction == 0) ? + shape_info.data[0] : + shape_info.data[1], + face_no); Evalf1 evalf1 = - create_evaluator_tensor_product(shape_info, face_no); + create_evaluator_tensor_product((face_direction == 1) ? + shape_info.data[0] : + shape_info.data[1], + face_no); Evalf2 evalf2 = - create_evaluator_tensor_product(shape_info, face_no); + 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; @@ -3787,8 +3844,7 @@ namespace internal do_evaluate, add_into_output, max_derivative, - lex_faces, - 0>(input, output); + lex_faces>(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; @@ -3797,8 +3853,7 @@ namespace internal do_evaluate, add_into_output, max_derivative, - lex_faces, - 1>(input, output); + lex_faces>(input, output); if (dim == 3) { @@ -3811,8 +3866,7 @@ namespace internal do_evaluate, add_into_output, max_derivative, - lex_faces, - 2>(input, output); + lex_faces>(input, output); } } }; @@ -4209,12 +4263,13 @@ namespace internal (n_q_points_1d < 1) ? 1 : n_q_points_1d, Number>:: - template evaluate_or_integrate_in_face(evaluation_flag, - temp, - fe_eval, - scratch_data, - subface_index, - fe_eval.get_face_no()); + template evaluate_or_integrate_in_face( + evaluation_flag, + temp, + fe_eval, + scratch_data, + subface_index, + fe_eval.get_face_no()); } else if (fe_degree > -1 && subface_index >= GeometryInfo::max_children_per_cell && diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 2e8885f01f..3ff7d54a81 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -122,6 +122,38 @@ namespace internal 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 + {}; + /** @@ -198,30 +230,21 @@ namespace internal (void)dummy2; } - template + template void values(const Number in[], Number out[]) const { apply(shape_values, in, out); } - template + template void gradients(const Number in[], Number out[]) const { apply(shape_gradients, in, out); } - template + template void hessians(const Number in[], Number out[]) const { @@ -2364,13 +2387,15 @@ namespace internal int n_rows, int n_columns, typename Number, + int normal_dir, typename Number2> - struct EvaluatorTensorProduct + struct EvaluatorTensorProductAnisotropic { static constexpr unsigned int n_rows_of_product = numbers::invalid_unsigned_int; @@ -2381,7 +2406,7 @@ namespace internal * Empty constructor. Does nothing. Be careful when using 'values' and * related methods because they need to be filled with the other pointer */ - EvaluatorTensorProduct() + EvaluatorTensorProductAnisotropic() : shape_values(nullptr) , shape_gradients(nullptr) , shape_hessians(nullptr) @@ -2390,11 +2415,12 @@ namespace internal /** * 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) + EvaluatorTensorProductAnisotropic( + 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()) @@ -2417,31 +2443,25 @@ namespace internal (void)dummy2; } - template + template void values(const Number in[], Number out[]) const { - apply(shape_values, - in, - out); + apply(shape_values, in, out); } - template + template void gradients(const Number in[], Number out[]) const { - apply(shape_gradients, - in, - out); + apply(shape_gradients, in, out); } - template + template void hessians(const Number in[], Number out[]) const { - apply(shape_hessians, - in, - out); + apply(shape_hessians, in, out); } /** @@ -2475,7 +2495,6 @@ namespace internal template static void apply(const Number2 *DEAL_II_RESTRICT shape_data, @@ -2486,8 +2505,7 @@ namespace internal bool contract_onto_face, bool add, int max_derivative, - bool lex_faces = false, - int normal_direction> + bool lex_faces = false> void apply_face(const Number *DEAL_II_RESTRICT in, Number *DEAL_II_RESTRICT out) const; @@ -2501,22 +2519,20 @@ namespace internal int n_rows, int n_columns, typename Number, + int normal_dir, typename Number2> - template + template inline void - EvaluatorTensorProduct::apply(const Number2 *DEAL_II_RESTRICT - shape_data, - const Number *in, - Number * out) + EvaluatorTensorProductAnisotropic< + evaluate_raviart_thomas, + dim, + n_rows, + n_columns, + Number, + normal_dir, + Number2>::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."); @@ -2596,22 +2612,23 @@ namespace internal int n_rows, int n_columns, typename Number, + int normal_dir, typename Number2> template + bool lex_faces> inline void - EvaluatorTensorProduct::apply_face(const Number *DEAL_II_RESTRICT in, - Number *DEAL_II_RESTRICT - out) const + EvaluatorTensorProductAnisotropic< + evaluate_raviart_thomas, + dim, + n_rows, + n_columns, + Number, + normal_dir, + 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")); @@ -2625,34 +2642,33 @@ namespace internal // Determine the number of blocks depending on the face and normaldirection, // as well as dimension. - constexpr int n_blocks1 = (face_direction == normal_direction) ? - (n_rows - 1) : - ((face_direction == 0 && normal_direction == 2) || - (face_direction == 1 && normal_direction == 2) || - (face_direction == 2 && normal_direction == 1)) ? - n_rows : - (n_rows + 1); - constexpr int n_blocks2 = - (dim == 2) ? 1 : - ((face_direction == normal_direction) ? - (n_rows - 1) : - (((face_direction == 0 && normal_direction == 1) || - (face_direction == 1 && normal_direction == 0) || - (face_direction == 2 && normal_direction == 0)) ? - n_rows : - (n_rows + 1))); + 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 stride = - (face_direction == normal_direction) ? + (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_direction == 0) ? (n_rows + 1) : - n_rows))); + ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) : + n_rows))); constexpr int out_stride = n_blocks1 * n_blocks2; const Number *DEAL_II_RESTRICT shape_values = this->shape_values; @@ -2739,21 +2755,21 @@ namespace internal // product. Need to take that into account. if (dim == 3) { - if (normal_direction == 0) + if (normal_dir == 0) { if (contract_onto_face) out += n_rows - 1; else in += n_rows - 1; } - if (normal_direction == 1) + if (normal_dir == 1) { if (contract_onto_face) out += n_rows - 2; else in += n_rows - 2; } - if (normal_direction == 2) + if (normal_dir == 2) { if (contract_onto_face) out += n_rows; @@ -2777,17 +2793,17 @@ namespace internal // adjust for local coordinate system zx if (contract_onto_face) { - if (normal_direction == 0) + if (normal_dir == 0) { in += (n_rows + 1) * (n_rows - 1); out -= n_rows * (n_rows + 1) - 1; } - if (normal_direction == 1) + if (normal_dir == 1) { in += (n_rows - 1) * (n_rows - 1); out -= (n_rows - 1) * (n_rows - 1) - 1; } - if (normal_direction == 2) + if (normal_dir == 2) { in += (n_rows - 1) * (n_rows); out -= (n_rows) * (n_rows + 1) - 1; @@ -2795,17 +2811,17 @@ namespace internal } else { - if (normal_direction == 0) + if (normal_dir == 0) { out += (n_rows + 1) * (n_rows - 1); in -= n_rows * (n_rows + 1) - 1; } - if (normal_direction == 1) + if (normal_dir == 1) { out += (n_rows - 1) * (n_rows - 1); in -= (n_rows - 1) * (n_rows - 1) - 1; } - if (normal_direction == 2) + if (normal_dir == 2) { out += (n_rows - 1) * (n_rows); in -= (n_rows) * (n_rows + 1) - 1; -- 2.39.5