From d335a2ab83996acd377265ab70742afe75e04ac7 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 12 Sep 2023 14:26:06 +0200 Subject: [PATCH] Split up evaluation_kernels.h into cell and face parts --- .../deal.II/matrix_free/evaluation_kernels.h | 3410 +---------------- .../matrix_free/evaluation_kernels_face.h | 2920 ++++++++++++++ ...aluation_template_face_factory.templates.h | 3 +- .../evaluation_template_factory_internal.h | 15 + include/deal.II/matrix_free/fe_evaluation.h | 1 + .../deal.II/matrix_free/fe_point_evaluation.h | 2 +- .../matrix_free/tensor_product_kernels.h | 186 + 7 files changed, 3297 insertions(+), 3240 deletions(-) create mode 100644 include/deal.II/matrix_free/evaluation_kernels_face.h diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index e63dc157b9..afdf62854a 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -19,12 +19,9 @@ #include -#include -#include #include #include -#include #include #include #include @@ -1825,164 +1822,6 @@ namespace internal Number *values_dofs_actual, FEEvaluationData &fe_eval, const bool add_into_values_array = false); - - template - static void - work_normal(const MatrixFreeFunctions::UnivariateShapeData &data, - const Number *in, - Number *out, - const bool add_into_result = false, - const int subface_index_1d = 0) - { - AssertIndexRange(direction, dim); - AssertDimension(fe_degree, data.fe_degree); - AssertDimension(n_q_points_1d, data.n_q_points_1d); - constexpr int n_rows = fe_degree + 1; - constexpr int n_columns = n_q_points_1d; - constexpr int mm = contract_over_rows ? n_rows : n_columns; - constexpr int nn = contract_over_rows ? n_columns : n_rows; - const Number2 *shape_data = - symmetric_evaluate ? - data.shape_values_eo.data() : - data.values_within_subface[subface_index_1d].data(); - Assert(shape_data != nullptr, ExcNotInitialized()); - Assert(contract_over_rows == false || !add_into_result, - ExcMessage("Cannot add into result if contract_over_rows = true")); - - constexpr int n_blocks1 = Utilities::pow(fe_degree, direction); - constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1); - constexpr int stride_in = contract_over_rows ? 1 : extra_stride; - constexpr int stride_out = contract_over_rows ? extra_stride : 1; - constexpr EvaluatorVariant variant = - symmetric_evaluate ? evaluate_evenodd : evaluate_general; - - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - if (contract_over_rows == false && add_into_result) - apply_matrix_vector_product(shape_data, in, out); - else - apply_matrix_vector_product(shape_data, in, out); - - in += stride_in; - out += stride_out; - } - in += n_blocks1 * (mm - 1) * stride_in; - out += n_blocks1 * (nn - 1) * stride_out; - } - } - - template - static void - work_tangential( - const MatrixFreeFunctions::UnivariateShapeData &data, - const Number *in, - Number *out, - const int subface_index_1d = 0) - { - AssertIndexRange(direction, dim); - AssertDimension(fe_degree - 1, data.fe_degree); - AssertDimension(n_q_points_1d, data.n_q_points_1d); - static_assert(direction != normal_direction, - "Cannot interpolate tangentially in normal direction"); - - constexpr int n_rows = fe_degree; - constexpr int n_columns = n_q_points_1d; - const Number2 *shape_data = - symmetric_evaluate ? - data.shape_values_eo.data() : - data.values_within_subface[subface_index_1d].data(); - Assert(shape_data != nullptr, ExcNotInitialized()); - - constexpr int n_blocks1 = - (direction > normal_direction) ? - Utilities::pow(n_q_points_1d, direction) : - (direction > 0 ? - (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) : - 1); - constexpr int n_blocks2 = - (direction > normal_direction) ? - Utilities::pow(fe_degree, dim - 1 - direction) : - ((direction + 1 < dim) ? - (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) : - 1); - - constexpr EvaluatorVariant variant = - symmetric_evaluate ? evaluate_evenodd : evaluate_general; - - // Since we may perform an in-place interpolation, we must run the step - // expanding the size of the basis backward ('contract_over_rows' aka - // 'evaluate' case), so shift the pointers and decrement during the loop - if (contract_over_rows) - { - in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1; - out += extra_stride * - ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1); - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - apply_matrix_vector_product(shape_data, in, out); - - --in; - out -= extra_stride; - } - in -= n_blocks1 * (n_rows - 1); - out -= n_blocks1 * (n_columns - 1) * extra_stride; - } - } - else - { - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - apply_matrix_vector_product(shape_data, in, out); - - in += extra_stride; - ++out; - } - in += n_blocks1 * (n_columns - 1) * extra_stride; - out += n_blocks1 * (n_rows - 1); - } - } - } }; @@ -1999,7 +1838,7 @@ namespace internal const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData &fe_eval, - const bool add_into_values_array) + const bool add) { Assert(dim == 2 || dim == 3, ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas " @@ -2025,6 +1864,9 @@ namespace internal if (integrate) { + EvaluatorTensorProductAnisotropic + eval; + const bool do_values = evaluation_flag & EvaluationFlags::values; if ((evaluation_flag & EvaluationFlags::gradients) != 0u) integrate_gradients_collocation(shape_data[0], @@ -2032,12 +1874,9 @@ namespace internal gradients, do_values); if constexpr (dim > 2) - work_tangential<2, 0, false>(shape_data[1], values, values); - work_tangential<1, 0, false>(shape_data[1], values, values); - work_normal<0, false>(shape_data[0], - values, - values_dofs, - add_into_values_array); + eval.template tangential<2, 0>(shape_data[1], values, values); + eval.template tangential<1, 0>(shape_data[1], values, values); + eval.template normal<0>(shape_data[0], values, values_dofs, add); values += n_points; gradients += n_points * dim; @@ -2049,12 +1888,9 @@ namespace internal gradients, do_values); if constexpr (dim > 2) - work_tangential<2, 1, false>(shape_data[1], values, values); - work_tangential<0, 1, false>(shape_data[1], values, values); - work_normal<1, false>(shape_data[0], - values, - values_dofs, - add_into_values_array); + eval.template tangential<2, 1>(shape_data[1], values, values); + eval.template tangential<0, 1>(shape_data[1], values, values); + eval.template normal<1>(shape_data[0], values, values_dofs, add); if constexpr (dim > 2) { @@ -2067,20 +1903,19 @@ namespace internal values, gradients, do_values); - work_tangential<1, 2, false>(shape_data[1], values, values); - work_tangential<0, 2, false>(shape_data[1], values, values); - work_normal<2, false>(shape_data[0], - values, - values_dofs, - add_into_values_array); + eval.template tangential<1, 2>(shape_data[1], values, values); + eval.template tangential<0, 2>(shape_data[1], values, values); + eval.template normal<0>(shape_data[0], values, values_dofs, add); } } else { - work_normal<0, true>(shape_data[0], values_dofs, values); - work_tangential<1, 0, true>(shape_data[1], values, values); + EvaluatorTensorProductAnisotropic + eval; + eval.template normal<0>(shape_data[0], values_dofs, values); + eval.template tangential<1, 0>(shape_data[1], values, values); if constexpr (dim > 2) - work_tangential<2, 0, true>(shape_data[1], values, values); + eval.template tangential<2, 0>(shape_data[1], values, values); if ((evaluation_flag & EvaluationFlags::gradients) != 0u) evaluate_gradients_collocation(shape_data[0], values, @@ -2090,10 +1925,10 @@ namespace internal gradients += n_points * dim; values_dofs += dofs_per_component; - work_normal<1, true>(shape_data[0], values_dofs, values); - work_tangential<0, 1, true>(shape_data[1], values, values); + eval.template normal<1>(shape_data[0], values_dofs, values); + eval.template tangential<0, 1>(shape_data[1], values, values); if constexpr (dim > 2) - work_tangential<2, 1, true>(shape_data[1], values, values); + eval.template tangential<2, 1>(shape_data[1], values, values); if ((evaluation_flag & EvaluationFlags::gradients) != 0u) evaluate_gradients_collocation(shape_data[0], values, @@ -2105,9 +1940,9 @@ namespace internal gradients += n_points * dim; values_dofs += dofs_per_component; - work_normal<2, true>(shape_data[0], values_dofs, values); - work_tangential<0, 2, true>(shape_data[1], values, values); - work_tangential<1, 2, true>(shape_data[1], values, values); + eval.template normal<2>(shape_data[0], values_dofs, values); + eval.template tangential<0, 2>(shape_data[1], values, values); + eval.template tangential<1, 2>(shape_data[1], values, values); if ((evaluation_flag & EvaluationFlags::gradients) != 0u) evaluate_gradients_collocation(shape_data[0], values, @@ -2118,22 +1953,6 @@ namespace internal - /** - * Helper function to specify whether transformation to collocation should - * be used: It should give correct results (first condition), we need to be - * able to initialize the fields in shape_info.templates.h from the - * polynomials (second condition), and it should be the most efficient - * choice in terms of operation counts (third condition). - */ - constexpr bool - use_collocation_evaluation(const unsigned int fe_degree, - const unsigned int n_q_points_1d) - { - return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) && - (n_q_points_1d <= 3 * fe_degree / 2 + 1); - } - - /** * This class chooses an appropriate evaluation/integration strategy based on * the template parameters and the shape_info variable which contains runtime @@ -2384,3055 +2203,187 @@ namespace internal - template - struct FEFaceEvaluationImpl + /** + * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation, + * using an FEEvaluationData argument. + */ + template + struct CellwiseInverseMassMatrixImplBasic { - // We enable a transformation to collocation for derivatives if it gives - // correct results (first two conditions), if it is the most efficient - // choice in terms of operation counts (third condition) and if we were - // able to initialize the fields in shape_info.templates.h from the - // polynomials (fourth condition). using Number2 = - typename FEEvaluationData::shape_info_number_type; - - using Eval = EvaluatorTensorProduct; - - static Eval - create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &data, - const unsigned int subface_index, - const unsigned int direction) - { - if (symmetric_evaluate) - return Eval(data.shape_values_eo, - data.shape_gradients_eo, - data.shape_hessians_eo, - data.fe_degree + 1, - data.n_q_points_1d); - else if (subface_index >= GeometryInfo::max_children_per_cell) - return Eval(data.shape_values, - data.shape_gradients, - data.shape_hessians, - data.fe_degree + 1, - data.n_q_points_1d); - else - { - const unsigned int index = - direction == 0 ? subface_index % 2 : subface_index / 2; - return Eval(data.values_within_subface[index], - data.gradients_within_subface[index], - data.hessians_within_subface[index], - data.fe_degree + 1, - data.n_q_points_1d); - } - } + typename FEEvaluationData::shape_info_number_type; - static void - evaluate_in_face( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::UnivariateShapeData &data, - Number *values_dofs, - Number *values_quad, - Number *gradients_quad, - Number *hessians_quad, - Number *scratch_data, - const unsigned int subface_index) + template + static bool + run(const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number *in_array, + Number *out_array) { - Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); - Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); + const unsigned int given_degree = + (fe_degree > -1) ? fe_degree : + fe_eval.get_shape_info().data.front().fe_degree; - const std::size_t n_dofs = fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(data.fe_degree + 1, dim - 1); - const std::size_t n_q_points = - fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : - Utilities::pow(data.n_q_points_1d, dim - 1); + const unsigned int dofs_per_component = + Utilities::pow(given_degree + 1, dim); - // keep a copy of the original pointer for the case of the Hessians - Number *values_dofs_ptr = values_dofs; + Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); + Assert(fe_eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric_no_collocation, + ExcNotImplemented()); - if ((evaluation_flag & EvaluationFlags::values) != 0u && - ((evaluation_flag & EvaluationFlags::gradients) == 0u)) - for (unsigned int c = 0; c < n_components; ++c) - { - switch (dim) - { - case 3: - 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>(values_dofs, - values_quad); - break; - case 1: - values_quad[0] = values_dofs[0]; - break; - default: - Assert(false, ExcNotImplemented()); - } - // Note: we always keep storage of values, 1st and 2nd derivatives - // in an array - values_dofs += 3 * n_dofs; - values_quad += n_q_points; - } - else if ((evaluation_flag & EvaluationFlags::gradients) != 0u) - for (unsigned int c = 0; c < n_components; ++c) - { - switch (dim) - { - case 3: - if (symmetric_evaluate && - use_collocation_evaluation(fe_degree, n_q_points_1d)) - { - eval0.template values<0, true, false>(values_dofs, - values_quad); - eval0.template values<1, true, false>(values_quad, - values_quad); - EvaluatorTensorProduct - eval_grad({}, data.shape_gradients_collocation_eo, {}); - eval_grad.template gradients<0, true, false, 3>( - values_quad, gradients_quad); - eval_grad.template gradients<1, true, false, 3>( - values_quad, gradients_quad + 1); - } - else - { - // grad x - eval0.template gradients<0, true, false>(values_dofs, - scratch_data); - eval1.template values<1, true, false, 3>(scratch_data, - gradients_quad); - - // grad y - eval0.template values<0, true, false>(values_dofs, - scratch_data); - eval1.template gradients<1, true, false, 3>( - scratch_data, gradients_quad + 1); - - if ((evaluation_flag & EvaluationFlags::values) != 0u) - eval1.template values<1, true, false>(scratch_data, - values_quad); - } - // grad z - eval0.template values<0, true, false>(values_dofs + n_dofs, - scratch_data); - eval1.template values<1, true, false, 3>(scratch_data, - gradients_quad + 2); - - break; - case 2: - eval0.template values<0, true, false, 2>(values_dofs + n_dofs, - gradients_quad + 1); - eval0.template gradients<0, true, false, 2>(values_dofs, - gradients_quad); - if ((evaluation_flag & EvaluationFlags::values) != 0u) - eval0.template values<0, true, false>(values_dofs, - values_quad); - break; - case 1: - values_quad[0] = values_dofs[0]; - gradients_quad[0] = values_dofs[1]; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - values_dofs += 3 * n_dofs; - values_quad += n_q_points; - gradients_quad += dim * n_q_points; - } + EvaluatorTensorProduct + evaluator({}, + {}, + fe_eval.get_shape_info().data.front().inverse_shape_values_eo, + given_degree + 1, + given_degree + 1); - if ((evaluation_flag & EvaluationFlags::hessians) != 0u) + for (unsigned int d = 0; d < n_components; ++d) { - values_dofs = values_dofs_ptr; - for (unsigned int c = 0; c < n_components; ++c) - { - switch (dim) - { - case 3: - // grad xx - 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>(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>(values_dofs + - 2 * n_dofs, - scratch_data); - eval1.template values<1, true, false>(scratch_data, - hessians_quad + - 2 * n_q_points); - - // grad xy - 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>(values_dofs + - n_dofs, - scratch_data); - eval1.template values<1, true, false>(scratch_data, - hessians_quad + - 4 * n_q_points); - - // grad yz - eval0.template values<0, true, false>(values_dofs + n_dofs, - 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>(values_dofs, - hessians_quad); - // grad yy - eval0.template values<0, true, false>( - values_dofs + 2 * n_dofs, hessians_quad + n_q_points); - // grad xy - eval0.template gradients<0, true, false>( - values_dofs + n_dofs, hessians_quad + 2 * n_q_points); - break; - case 1: - hessians_quad[0] = values_dofs[2]; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - values_dofs += 3 * n_dofs; - hessians_quad += dim * (dim + 1) / 2 * n_q_points; - } + const Number *in = in_array + d * dofs_per_component; + Number *out = out_array + d * dofs_per_component; + // Need to select 'apply' method with hessian slot because values + // assume symmetries that do not exist in the inverse shapes + evaluator.template hessians<0, true, false>(in, out); + if (dim > 1) + evaluator.template hessians<1, true, false>(out, out); + if (dim > 2) + evaluator.template hessians<2, true, false>(out, out); } - } - - static void - integrate_in_face( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::UnivariateShapeData &data, - Number *values_dofs, - Number *values_quad, - Number *gradients_quad, - Number *hessians_quad, - Number *scratch_data, - const unsigned int subface_index) - { - Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); - Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); - - const std::size_t n_dofs = - fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - (dim > 1 ? Utilities::fixed_power(data.fe_degree + 1) : 1); - const std::size_t n_q_points = - fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : - Utilities::pow(data.n_q_points_1d, dim - 1); - - // keep a copy of the original pointer for the case of the Hessians - Number *values_dofs_ptr = values_dofs; - - if ((integration_flag & EvaluationFlags::values) != 0u && - (integration_flag & EvaluationFlags::gradients) == 0u) - for (unsigned int c = 0; c < n_components; ++c) - { - switch (dim) - { - case 3: - 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>(values_quad, - values_dofs); - break; - case 1: - values_dofs[0] = values_quad[0]; - break; - default: - Assert(false, ExcNotImplemented()); - } - values_dofs += 3 * n_dofs; - values_quad += n_q_points; - } - else if ((integration_flag & EvaluationFlags::gradients) != 0u) - for (unsigned int c = 0; c < n_components; ++c) - { - switch (dim) - { - case 3: - // grad z - eval1.template values<1, false, false, 3>(gradients_quad + 2, - scratch_data); - eval0.template values<0, false, false>(scratch_data, - values_dofs + n_dofs); - if (symmetric_evaluate && - use_collocation_evaluation(fe_degree, n_q_points_1d)) - { - EvaluatorTensorProduct - eval_grad({}, data.shape_gradients_collocation_eo, {}); - if ((integration_flag & EvaluationFlags::values) != 0u) - eval_grad.template gradients<1, false, true, 3>( - gradients_quad + 1, values_quad); - else - eval_grad.template gradients<1, false, false, 3>( - gradients_quad + 1, values_quad); - eval_grad.template gradients<0, false, true, 3>( - gradients_quad, values_quad); - eval0.template values<1, false, false>(values_quad, - values_quad); - eval0.template values<0, false, false>(values_quad, - values_dofs); - } - else - { - if ((integration_flag & EvaluationFlags::values) != 0u) - { - eval1.template values<1, false, false>(values_quad, - scratch_data); - eval1.template gradients<1, false, true, 3>( - gradients_quad + 1, scratch_data); - } - else - eval1.template gradients<1, false, false, 3>( - gradients_quad + 1, scratch_data); - - // grad y - eval0.template values<0, false, false>(scratch_data, - values_dofs); - - // grad x - eval1.template values<1, false, false, 3>(gradients_quad, - scratch_data); - eval0.template gradients<0, false, true>(scratch_data, - values_dofs); - } - break; - case 2: - eval0.template values<0, false, false, 2>(gradients_quad + 1, - values_dofs + - n_dofs); - eval0.template gradients<0, false, false, 2>(gradients_quad, - values_dofs); - if ((integration_flag & EvaluationFlags::values) != 0u) - eval0.template values<0, false, true>(values_quad, - values_dofs); - break; - case 1: - values_dofs[0] = values_quad[0]; - values_dofs[1] = gradients_quad[0]; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - values_dofs += 3 * n_dofs; - values_quad += n_q_points; - gradients_quad += dim * n_q_points; - } - - if ((integration_flag & EvaluationFlags::hessians) != 0u) + for (unsigned int q = 0; q < dofs_per_component; ++q) { - values_dofs = values_dofs_ptr; - for (unsigned int c = 0; c < n_components; ++c) - { - switch (dim) - { - case 3: - // grad xx - eval1.template values<1, false, false>(hessians_quad, - scratch_data); - if ((integration_flag & (EvaluationFlags::values | - EvaluationFlags::gradients)) != 0u) - eval0.template hessians<0, false, true>(scratch_data, - values_dofs); - else - eval0.template hessians<0, false, false>(scratch_data, - values_dofs); - - // grad yy - 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>(hessians_quad + - 2 * n_q_points, - scratch_data); - eval0.template values<0, false, false>(scratch_data, - values_dofs + - 2 * n_dofs); - - // grad xy - 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>(hessians_quad + - 4 * n_q_points, - scratch_data); - if ((integration_flag & EvaluationFlags::gradients) != 0u) - eval0.template gradients<0, false, true>(scratch_data, - values_dofs + - n_dofs); - else - eval0.template gradients<0, false, false>(scratch_data, - values_dofs + - n_dofs); - - // grad yz - eval1.template gradients<1, false, false>(hessians_quad + - 5 * n_q_points, - scratch_data); - eval0.template values<0, false, true>(scratch_data, - values_dofs + n_dofs); - - break; - case 2: - // grad xx - if ((integration_flag & (EvaluationFlags::values | - EvaluationFlags::gradients)) != 0u) - eval0.template hessians<0, false, true>(hessians_quad, - values_dofs); - else - eval0.template hessians<0, false, false>(hessians_quad, - values_dofs); - - // grad yy - eval0.template values<0, false, false>( - hessians_quad + n_q_points, values_dofs + 2 * n_dofs); - // grad xy - if ((integration_flag & EvaluationFlags::gradients) != 0u) - eval0.template gradients<0, false, true>( - hessians_quad + 2 * n_q_points, values_dofs + n_dofs); - else - eval0.template gradients<0, false, false>( - hessians_quad + 2 * n_q_points, values_dofs + n_dofs); - break; - case 1: - values_dofs[2] = hessians_quad[0]; - if ((integration_flag & EvaluationFlags::values) == 0u) - values_dofs[0] = 0; - if ((integration_flag & EvaluationFlags::gradients) == 0u) - values_dofs[1] = 0; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - values_dofs += 3 * n_dofs; - hessians_quad += dim * (dim + 1) / 2 * n_q_points; - } + const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q); + for (unsigned int d = 0; d < n_components; ++d) + out_array[q + d * dofs_per_component] *= inverse_JxW_q; } + for (unsigned int d = 0; d < n_components; ++d) + { + Number *out = out_array + d * dofs_per_component; + if (dim > 2) + evaluator.template hessians<2, false, false>(out, out); + if (dim > 1) + evaluator.template hessians<1, false, false>(out, out); + evaluator.template hessians<0, false, false>(out, out); + } + return false; } }; - template - struct FEFaceEvaluationImplRaviartThomas + /** + * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation + * with user-provided coefficients at quadrature points (in contrast to + * CellwiseInverseMassMatrixImplBasic, which implicitly uses `1/(|J|xW)' as + * coefficient). + */ + template + struct CellwiseInverseMassMatrixImplFlexible { using Number2 = - typename FEEvaluationData::shape_info_number_type; + typename FEEvaluationData::shape_info_number_type; - /** - * Apply the sum factorization kernels within the face for Raviart-Thomas - * elements for either evaluation or integration - */ - template - static inline void - evaluate_or_integrate_in_face( - const EvaluationFlags::EvaluationFlags evaluation_flag, - const std::vector> - &shape_data, - Number *values_dofs_in, - Number *values_quad, - Number *gradients_quad, - Number *scratch_data, - const unsigned int subface_index, - const unsigned int face_direction) + template + static bool + run(const unsigned int n_desired_components, + const FEEvaluationData &fe_eval, + const ArrayView &inverse_coefficients, + const bool dyadic_coefficients, + const Number *in_array, + Number *out_array) { - AssertDimension(shape_data.size(), 2); - - const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree; - const int n_rows_n = degree + 1; - const int n_rows_t = degree; - 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}}}}; - - (void)scratch_data; - (void)subface_index; - // TODO: This is currently not implemented, but the test - // matrix_vector_rt_face_03 apparently works without it -> check - // if (subface_index < GeometryInfo::max_children_per_cell) - // Assert(false, ExcNotImplemented()); - - using EvalAniso = - FEEvaluationImpl; - using Eval = EvaluatorTensorProduct; - - std::array values_dofs_offsets = {}; - for (unsigned int comp = 0; comp < dim - 1; ++comp) + const unsigned int given_degree = + (fe_degree > -1) ? fe_degree : + fe_eval.get_shape_info().data.front().fe_degree; + + const unsigned int dofs_per_component = + Utilities::pow(given_degree + 1, dim); + + Assert(inverse_coefficients.size() > 0 && + inverse_coefficients.size() % dofs_per_component == 0, + ExcMessage( + "Expected diagonal to be a multiple of scalar dof per cells")); + + if (!dyadic_coefficients) { - if (dim == 2) - values_dofs_offsets[comp + 1] = - values_dofs_offsets[comp] + - 3 * dofs_per_direction[comp][(face_direction + 1) % dim]; - else - values_dofs_offsets[comp + 1] = - values_dofs_offsets[comp] + - 3 * dofs_per_direction[comp][(face_direction + 1) % dim] * - dofs_per_direction[comp][(face_direction + 2) % dim]; + if (inverse_coefficients.size() != dofs_per_component) + AssertDimension(n_desired_components * dofs_per_component, + inverse_coefficients.size()); + } + else + { + AssertDimension(n_desired_components * n_desired_components * + dofs_per_component, + inverse_coefficients.size()); } - // Jacobians on faces are reordered to enable simple access with the - // regular evaluators; to get the RT Piola transform right, we need to - // pass through the values_dofs array in a permuted right order - std::array components; - for (unsigned int comp = 0; comp < dim; ++comp) - components[comp] = (face_direction + comp + 1) % dim; + Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); + Assert(fe_eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric_no_collocation, + ExcNotImplemented()); - for (const unsigned int comp : components) - { - Number *values_dofs = values_dofs_in + values_dofs_offsets[comp]; + EvaluatorTensorProduct + evaluator({}, + {}, + fe_eval.get_shape_info().data.front().inverse_shape_values_eo, + given_degree + 1, + given_degree + 1); + + const Number *in = in_array; + Number *out = out_array; + + const Number *inv_coefficient = inverse_coefficients.data(); + + const unsigned int shift_coefficient = + inverse_coefficients.size() > dofs_per_component ? dofs_per_component : + 0; - std::array n_blocks{ - {dofs_per_direction[comp][(face_direction + 1) % dim], - (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] : - 1)}}; + const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components; + const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1; - if constexpr (dim == 3) + for (unsigned int d = 0; d < n_comp_outer; ++d) + { + for (unsigned int di = 0; di < n_comp_inner; ++di) + { + const Number *in_ = in + di * dofs_per_component; + Number *out_ = out + di * dofs_per_component; + evaluator.template hessians<0, true, false>(in_, out_); + if (dim > 1) + evaluator.template hessians<1, true, false>(out_, out_); + if (dim > 2) + evaluator.template hessians<2, true, false>(out_, out_); + } + if (dyadic_coefficients) { - EvaluatorTensorProduct - eval_g({}, - shape_data[0].shape_gradients_collocation_eo.data(), - {}); - if (!do_integrate) + const auto n_coeff_components = + n_desired_components * n_desired_components; + if (n_desired_components == dim) { - // Evaluate in 3d - if (n_blocks[0] == n_rows_n) - { - EvalAniso::template work_normal<0, true, true, 1>( - shape_data[0], values_dofs, values_quad); - EvalAniso::template work_tangential<1, 0, true, true, 1>( - shape_data[1], values_quad, values_quad); - - if (evaluation_flag & EvaluationFlags::gradients) - { - EvalAniso::template work_normal<0, true, true, 1>( - shape_data[0], - values_dofs + n_blocks[0] * n_blocks[1], - scratch_data); - EvalAniso:: - template work_tangential<1, 0, true, true, dim>( - shape_data[1], scratch_data, gradients_quad + 2); - } - } - else if (n_blocks[1] == n_rows_n) - { - EvalAniso::template work_normal<1, true, true, 1>( - shape_data[0], values_dofs, values_quad); - EvalAniso::template work_tangential<0, 1, true, true, 1>( - shape_data[1], values_quad, values_quad); - - if (evaluation_flag & EvaluationFlags::gradients) - { - EvalAniso::template work_normal<1, true, true, 1>( - shape_data[0], - values_dofs + n_blocks[0] * n_blocks[1], - scratch_data); - EvalAniso:: - template work_tangential<0, 1, true, true, dim>( - shape_data[1], scratch_data, gradients_quad + 2); - } - } - else - { - Eval eval(shape_data[1].shape_values_eo.data(), {}, {}); - eval.template values<0, true, false>(values_dofs, - values_quad); - eval.template values<1, true, false>(values_quad, - values_quad); - if (evaluation_flag & EvaluationFlags::gradients) - { - eval.template values<0, true, false>(values_dofs + - n_blocks[0] * - n_blocks[1], - scratch_data); - eval.template values<1, true, false, dim>( - scratch_data, gradients_quad + 2); - } - } - if (evaluation_flag & EvaluationFlags::gradients) - { - eval_g.template gradients<0, true, false, dim>( - values_quad, gradients_quad); - eval_g.template gradients<1, true, false, dim>( - values_quad, gradients_quad + 1); - } + for (unsigned int q = 0; q < dofs_per_component; ++q) + vmult(&inv_coefficient[q * n_coeff_components], + &in[q], + &out[q], + dofs_per_component); } else { - // Integrate in 3d - if (evaluation_flag & EvaluationFlags::gradients) - { - if (evaluation_flag & EvaluationFlags::values) - eval_g.template gradients<0, false, true, dim>( - gradients_quad, values_quad); - else - eval_g.template gradients<0, false, false, dim>( - gradients_quad, values_quad); - eval_g.template gradients<1, false, true, dim>( - gradients_quad + 1, values_quad); - } - if (n_blocks[0] == n_rows_n) - { - EvalAniso::template work_tangential<1, 0, false, true, 1>( - shape_data[1], values_quad, values_quad); - EvalAniso::template work_normal<0, false, true, 1>( - shape_data[0], values_quad, values_dofs); - - if (evaluation_flag & EvaluationFlags::gradients) - { - EvalAniso:: - template work_tangential<1, 0, false, true, dim>( - shape_data[1], gradients_quad + 2, scratch_data); - EvalAniso::template work_normal<0, false, true, 1>( - shape_data[0], - scratch_data, - values_dofs + n_blocks[0] * n_blocks[1]); - } - } - else if (n_blocks[1] == n_rows_n) - { - EvalAniso::template work_tangential<0, 1, false, true, 1>( - shape_data[1], values_quad, values_quad); - EvalAniso::template work_normal<1, false, true, 1>( - shape_data[0], values_quad, values_dofs); - - if (evaluation_flag & EvaluationFlags::gradients) - { - EvalAniso:: - template work_tangential<0, 1, false, true, dim>( - shape_data[1], gradients_quad + 2, scratch_data); - EvalAniso::template work_normal<1, false, true, 1>( - shape_data[0], - scratch_data, - values_dofs + n_blocks[0] * n_blocks[1]); - } - } - else - { - Eval eval(shape_data[1].shape_values_eo.data(), {}, {}); - eval.template values<1, false, false>(values_quad, - values_quad); - eval.template values<0, false, false>(values_quad, - values_dofs); - if (evaluation_flag & EvaluationFlags::gradients) - { - eval.template values<1, false, false, dim>( - gradients_quad + 2, scratch_data); - eval.template values<0, false, false>( - scratch_data, - values_dofs + n_blocks[0] * n_blocks[1]); - } - } - } - } - else - { - using EvalN = EvaluatorTensorProduct; - if (!do_integrate) - { - // Evaluate in 2d - if (n_blocks[0] == n_rows_n) - { - EvalN eval(shape_data[0].shape_values_eo, - shape_data[0].shape_gradients_eo, - {}); - eval.template values<0, true, false>(values_dofs, - values_quad); - if (evaluation_flag & EvaluationFlags::gradients) - { - eval.template gradients<0, true, false, dim>( - values_dofs, gradients_quad); - eval.template values<0, true, false, dim>( - values_dofs + n_rows_n, gradients_quad + 1); - } - } - else - { - Eval eval(shape_data[1].shape_values_eo, - shape_data[1].shape_gradients_eo, - {}); - eval.template values<0, true, false>(values_dofs, - values_quad); - if (evaluation_flag & EvaluationFlags::gradients) - { - eval.template gradients<0, true, false, dim>( - values_dofs, gradients_quad); - eval.template values<0, true, false, dim>( - values_dofs + n_rows_t, gradients_quad + 1); - } - } - } - else - { - // Integrate in 2d - if (n_blocks[0] == n_rows_n) - { - EvalN eval(shape_data[0].shape_values_eo, - shape_data[0].shape_gradients_eo, - {}); - if (evaluation_flag & EvaluationFlags::values) - eval.template values<0, false, false>(values_quad, - values_dofs); - if (evaluation_flag & EvaluationFlags::gradients) - { - if (evaluation_flag & EvaluationFlags::values) - eval.template gradients<0, false, true, dim>( - gradients_quad, values_dofs); - else - eval.template gradients<0, false, false, dim>( - gradients_quad, values_dofs); - eval.template values<0, false, false, dim>( - gradients_quad + 1, values_dofs + n_rows_n); - } - } - else - { - Eval eval(shape_data[1].shape_values_eo, - shape_data[1].shape_gradients_eo, - {}); - if (evaluation_flag & EvaluationFlags::values) - eval.template values<0, false, false>(values_quad, - values_dofs); - if (evaluation_flag & EvaluationFlags::gradients) - { - if (evaluation_flag & EvaluationFlags::values) - eval.template gradients<0, false, true, dim>( - gradients_quad, values_dofs); - else - eval.template gradients<0, false, false, dim>( - gradients_quad, values_dofs); - eval.template values<0, false, false, dim>( - gradients_quad + 1, values_dofs + n_rows_t); - } - } - } - } - values_quad += Utilities::pow(n_q_points_1d, dim - 1); - gradients_quad += dim * Utilities::pow(n_q_points_1d, dim - 1); - } - } - }; - - - - template - struct FEFaceNormalEvaluationImpl - { - using Number2 = - typename FEEvaluationData::shape_info_number_type; - - template - static void - interpolate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags flags, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number *input, - Number *output, - const unsigned int face_no) - { - Assert(static_cast(fe_degree) == - shape_info.data.front().fe_degree || - fe_degree == -1, - ExcInternalError()); - if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) - interpolate_raviart_thomas( - n_components, input, output, flags, face_no, shape_info); - else - interpolate_generic( - n_components, - input, - output, - flags, - face_no, - shape_info.data.front().fe_degree + 1, - shape_info.data.front().shape_data_on_face, - shape_info.dofs_per_component_on_cell, - 3 * shape_info.dofs_per_component_on_face); - } - - /** - * Interpolate the values on the cell quadrature points onto a face. - */ - template - static void - interpolate_quadrature( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags flags, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number *input, - Number *output, - const unsigned int face_no) - { - Assert(static_cast(fe_degree + 1) == - shape_info.data.front().n_q_points_1d || - fe_degree == -1, - ExcInternalError()); - - interpolate_generic( - n_components, - input, - output, - flags, - face_no, - shape_info.data.front().quadrature.size(), - shape_info.data.front().quadrature_data_on_face, - shape_info.n_q_points, - shape_info.n_q_points_face); - } - - private: - template - static void - interpolate_generic(const unsigned int n_components, - const Number *input, - Number *output, - const EvaluationFlags::EvaluationFlags flag, - const unsigned int face_no, - const unsigned int n_points_1d, - const std::array, 2> &shape_data, - const unsigned int dofs_per_component_on_cell, - const unsigned int dofs_per_component_on_face) - { - if (face_direction == face_no / 2) - { - constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction); - - const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d; - const int stride = Utilities::pow(n_rows, face_direction); - const std::array n_blocks{ - {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}; - 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}}; - - for (unsigned int c = 0; c < n_components; ++c) - { - if (flag & EvaluationFlags::hessians) - interpolate_to_face(shape_data[face_no % 2].begin(), - n_blocks, - steps, - input, - output, - n_rows, - stride); - else if (flag & EvaluationFlags::gradients) - interpolate_to_face(shape_data[face_no % 2].begin(), - n_blocks, - steps, - input, - output, - n_rows, - stride); - else - interpolate_to_face(shape_data[face_no % 2].begin(), - n_blocks, - steps, - input, - output, - n_rows, - stride); - if (do_evaluate) - { - input += dofs_per_component_on_cell; - output += dofs_per_component_on_face; - } - else - { - output += dofs_per_component_on_cell; - input += dofs_per_component_on_face; - } - } - } - else if (face_direction < dim) - { - interpolate_generic( - n_components, - input, - output, - flag, - face_no, - n_points_1d, - shape_data, - dofs_per_component_on_cell, - dofs_per_component_on_face); - } - } - - template - static void - interpolate_raviart_thomas( - const unsigned int n_components, - const Number *input, - Number *output, - const EvaluationFlags::EvaluationFlags flag, - const unsigned int face_no, - const MatrixFreeFunctions::ShapeInfo &shape_info) - { - if (dim == 1) - { - // This should never happen since the FE_RaviartThomasNodal is not - // defined for dim = 1. It prevents compiler warnings of infinite - // recursion. - Assert(false, ExcInternalError()); - return; - } - - bool increase_max_der = false; - if ((flag & EvaluationFlags::hessians && max_derivative < 2) || - (flag & EvaluationFlags::gradients && max_derivative < 1)) - increase_max_der = true; - - if (face_direction == face_no / 2 && !increase_max_der) - { - 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) - { - // Only increase max_derivative - interpolate_raviart_thomas( - n_components, input, output, flag, face_no, shape_info); - } - else if (face_direction < dim) - { - if (increase_max_der) - { - interpolate_raviart_thomas( - n_components, input, output, flag, face_no, shape_info); - } - else - { - interpolate_raviart_thomas( - n_components, input, output, flag, face_no, shape_info); - } - } - } - }; - - - - // internal helper function for reading data; base version of different types - template - void - do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst) - { - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - dst[v] = src_ptr[v]; - } - - - - // internal helper function for reading data; specialized version where we - // can use a dedicated load function - template - void - do_vectorized_read(const Number *src_ptr, VectorizedArray &dst) - { - dst.load(src_ptr); - } - - - - // internal helper function for reading data; base version of different types - template - void - do_vectorized_gather(const Number2 *src_ptr, - const unsigned int *indices, - VectorizedArrayType &dst) - { - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - dst[v] = src_ptr[indices[v]]; - } - - - - // internal helper function for reading data; specialized version where we - // can use a dedicated gather function - template - void - do_vectorized_gather(const Number *src_ptr, - const unsigned int *indices, - VectorizedArray &dst) - { - dst.gather(src_ptr, indices); - } - - - - // internal helper function for reading data; base version of different types - template - void - do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr) - { - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - dst_ptr[v] += src[v]; - } - - - - // internal helper function for reading data; specialized version where we - // can use a dedicated load function - template - void - do_vectorized_add(const VectorizedArray src, Number *dst_ptr) - { - VectorizedArray tmp; - tmp.load(dst_ptr); - (tmp + src).store(dst_ptr); - } - - - - // internal helper function for reading data; base version of different types - template - void - do_vectorized_scatter_add(const VectorizedArrayType src, - const unsigned int *indices, - Number2 *dst_ptr) - { - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - dst_ptr[indices[v]] += src[v]; - } - - - - // internal helper function for reading data; specialized version where we - // can use a dedicated gather function - template - void - do_vectorized_scatter_add(const VectorizedArray src, - const unsigned int *indices, - Number *dst_ptr) - { -#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512 - for (unsigned int v = 0; v < width; ++v) - dst_ptr[indices[v]] += src[v]; -#else - VectorizedArray tmp; - tmp.gather(dst_ptr, indices); - (tmp + src).scatter(indices, dst_ptr); -#endif - } - - - - template - void - adjust_for_face_orientation(const unsigned int dim, - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags flag, - const unsigned int *orientation, - const bool integrate, - const std::size_t n_q_points, - Number *tmp_values, - Number *values_quad, - Number *gradients_quad, - Number *hessians_quad) - { - for (unsigned int c = 0; c < n_components; ++c) - { - if (flag & EvaluationFlags::values) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = values_quad[c * n_q_points + orientation[q]]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = values_quad[c * n_q_points + q]; - for (unsigned int q = 0; q < n_q_points; ++q) - values_quad[c * n_q_points + q] = tmp_values[q]; - } - if (flag & EvaluationFlags::gradients) - for (unsigned int d = 0; d < dim; ++d) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = - gradients_quad[(c * n_q_points + orientation[q]) * dim + d]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = - gradients_quad[(c * n_q_points + q) * dim + d]; - for (unsigned int q = 0; q < n_q_points; ++q) - gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q]; - } - if (flag & EvaluationFlags::hessians) - { - const unsigned int hdim = (dim * (dim + 1)) / 2; - for (unsigned int d = 0; d < hdim; ++d) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points + - orientation[q]]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = - hessians_quad[(c * hdim + d) * n_q_points + q]; - for (unsigned int q = 0; q < n_q_points; ++q) - hessians_quad[(c * hdim + d) * n_q_points + q] = - tmp_values[q]; - } - } - } - } - - - - template - void - adjust_for_face_orientation_per_lane( - const unsigned int dim, - const unsigned int n_components, - const unsigned int v, - const EvaluationFlags::EvaluationFlags flag, - const unsigned int *orientation, - const bool integrate, - const std::size_t n_q_points, - Number *tmp_values, - VectorizedArrayType *values_quad, - VectorizedArrayType *gradients_quad = nullptr, - VectorizedArrayType *hessians_quad = nullptr) - { - for (unsigned int c = 0; c < n_components; ++c) - { - if (flag & EvaluationFlags::values) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v]; - for (unsigned int q = 0; q < n_q_points; ++q) - values_quad[c * n_q_points + q][v] = tmp_values[q]; - } - if (flag & EvaluationFlags::gradients) - for (unsigned int d = 0; d < dim; ++d) - { - Assert(gradients_quad != nullptr, ExcInternalError()); - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = - gradients_quad[(c * n_q_points + orientation[q]) * dim + d] - [v]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = - gradients_quad[(c * n_q_points + q) * dim + d][v]; - for (unsigned int q = 0; q < n_q_points; ++q) - gradients_quad[(c * n_q_points + q) * dim + d][v] = - tmp_values[q]; - } - if (flag & EvaluationFlags::hessians) - { - Assert(hessians_quad != nullptr, ExcInternalError()); - const unsigned int hdim = (dim * (dim + 1)) / 2; - for (unsigned int d = 0; d < hdim; ++d) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points + - orientation[q]][v]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = - hessians_quad[(c * hdim + d) * n_q_points + q][v]; - for (unsigned int q = 0; q < n_q_points; ++q) - hessians_quad[(c * hdim + d) * n_q_points + q][v] = - tmp_values[q]; - } - } - } - } - - - - template - struct FEFaceEvaluationImplEvaluateSelector - { - template - static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - FEEvaluationData &fe_eval) - { - const auto &shape_info = fe_eval.get_shape_info(); - const auto &shape_data = shape_info.data.front(); - using Number2 = - typename FEEvaluationData::shape_info_number_type; - - if (shape_info.element_type == MatrixFreeFunctions::tensor_none) - { - Assert((fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) == false, - ExcNotImplemented()); - - const unsigned int face_no = fe_eval.get_face_no(); - const unsigned int face_orientation = fe_eval.get_face_orientation(); - const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; - const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; - - using Eval = - EvaluatorTensorProduct; - - if (evaluation_flag & EvaluationFlags::values) - { - const auto *const shape_values = - &shape_data.shape_values_face(face_no, face_orientation, 0); - - auto *values_quad_ptr = fe_eval.begin_values(); - auto *values_dofs_actual_ptr = values_dofs; - - Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); - for (unsigned int c = 0; c < n_components; ++c) - { - eval.template values<0, true, false>(values_dofs_actual_ptr, - values_quad_ptr); - - values_quad_ptr += n_q_points; - values_dofs_actual_ptr += n_dofs; - } - } - - if (evaluation_flag & EvaluationFlags::gradients) - { - auto *gradients_quad_ptr = fe_eval.begin_gradients(); - const auto *values_dofs_actual_ptr = values_dofs; - - std::array shape_gradients; - for (unsigned int d = 0; d < dim; ++d) - shape_gradients[d] = &shape_data.shape_gradients_face( - face_no, face_orientation, d, 0); - - for (unsigned int c = 0; c < n_components; ++c) - { - for (unsigned int d = 0; d < dim; ++d) - { - Eval eval(nullptr, - shape_gradients[d], - nullptr, - n_dofs, - n_q_points); - - eval.template gradients<0, true, false, dim>( - values_dofs_actual_ptr, gradients_quad_ptr + d); - } - gradients_quad_ptr += n_q_points * dim; - values_dofs_actual_ptr += n_dofs; - } - } - - Assert(!(evaluation_flag & EvaluationFlags::hessians), - ExcNotImplemented()); - - return true; - } - - const unsigned int dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(shape_data.fe_degree + 1, dim - 1); - - // Note: we always keep storage of values, 1st and 2nd derivatives in an - // array, so reserve space for all three here - Number *temp = fe_eval.get_scratch_data().begin(); - Number *scratch_data = temp + 3 * n_components * dofs_per_face; - - bool use_vectorization = true; - - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) // exterior faces in the ECL loop - for (unsigned int v = 0; v < Number::size(); ++v) - if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int && - fe_eval.get_face_no(v) != fe_eval.get_face_no(0)) - use_vectorization = false; - - if (use_vectorization == false) - { - for (unsigned int v = 0; v < Number::size(); ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - { - for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; - ++i) - temp[i][v] = 0; - continue; - } - - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - evaluation_flag, - shape_info, - values_dofs, - scratch_data, - fe_eval.get_face_no(v)); - - for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; - ++i) - temp[i][v] = scratch_data[i][v]; - } - } - else - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - evaluation_flag, - shape_info, - values_dofs, - temp, - fe_eval.get_face_no()); - - const unsigned int subface_index = fe_eval.get_subface_index(); - constexpr unsigned int n_q_points_1d_actual = - fe_degree > -1 ? n_q_points_1d : 0; - - if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) - { - FEFaceEvaluationImplRaviartThomas:: - template evaluate_or_integrate_in_face( - evaluation_flag, - fe_eval.get_shape_info().data, - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - scratch_data, - subface_index, - fe_eval.get_face_no() / 2); - } - else if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl::evaluate_in_face(n_components, - evaluation_flag, - shape_data, - temp, - fe_eval.begin_values(), - fe_eval - .begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - else - FEFaceEvaluationImpl::evaluate_in_face(n_components, - evaluation_flag, - shape_data, - temp, - fe_eval.begin_values(), - fe_eval - .begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - - if (use_vectorization == false) - { - for (unsigned int v = 0; v < Number::size(); ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - if (fe_eval.get_face_orientation(v) != 0) - adjust_for_face_orientation_per_lane( - dim, - n_components, - v, - evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(v), 0), - false, - shape_info.n_q_points_face, - &temp[0][0], - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - } - } - else if (fe_eval.get_face_orientation() != 0) - adjust_for_face_orientation( - dim, - n_components, - evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(), 0), - false, - shape_info.n_q_points_face, - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - - return false; - } - }; - - - - template - struct FEFaceEvaluationImplIntegrateSelector - { - template - static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, - FEEvaluationData &fe_eval) - { - const auto &shape_info = fe_eval.get_shape_info(); - const auto &shape_data = shape_info.data.front(); - using Number2 = - typename FEEvaluationData::shape_info_number_type; - - if (shape_info.element_type == MatrixFreeFunctions::tensor_none) - { - Assert((fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) == false, - ExcNotImplemented()); - - const unsigned int face_no = fe_eval.get_face_no(); - const unsigned int face_orientation = fe_eval.get_face_orientation(); - const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; - const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; - - using Eval = - EvaluatorTensorProduct; - - if (integration_flag & EvaluationFlags::values) - { - const auto *const shape_values = - &shape_data.shape_values_face(face_no, face_orientation, 0); - - auto *values_quad_ptr = fe_eval.begin_values(); - auto *values_dofs_actual_ptr = values_dofs; - - Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); - for (unsigned int c = 0; c < n_components; ++c) - { - eval.template values<0, false, false>(values_quad_ptr, - values_dofs_actual_ptr); - - values_quad_ptr += n_q_points; - values_dofs_actual_ptr += n_dofs; - } - } - - if (integration_flag & EvaluationFlags::gradients) - { - auto *gradients_quad_ptr = fe_eval.begin_gradients(); - auto *values_dofs_actual_ptr = values_dofs; - - std::array shape_gradients; - for (unsigned int d = 0; d < dim; ++d) - shape_gradients[d] = &shape_data.shape_gradients_face( - face_no, face_orientation, d, 0); - - for (unsigned int c = 0; c < n_components; ++c) - { - for (unsigned int d = 0; d < dim; ++d) - { - Eval eval(nullptr, - shape_gradients[d], - nullptr, - n_dofs, - n_q_points); - - if (!(integration_flag & EvaluationFlags::values) && - d == 0) - eval.template gradients<0, false, false, dim>( - gradients_quad_ptr + d, values_dofs_actual_ptr); - else - eval.template gradients<0, false, true, dim>( - gradients_quad_ptr + d, values_dofs_actual_ptr); - } - gradients_quad_ptr += n_q_points * dim; - values_dofs_actual_ptr += n_dofs; - } - } - - Assert(!(integration_flag & EvaluationFlags::hessians), - ExcNotImplemented()); - - return true; - } - - const unsigned int dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(shape_data.fe_degree + 1, dim - 1); - - Number *temp = fe_eval.get_scratch_data().begin(); - Number *scratch_data = temp + 3 * n_components * dofs_per_face; - - bool use_vectorization = true; - - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) // exterior faces in the ECL loop - use_vectorization = - fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int && - std::all_of(fe_eval.get_cell_ids().begin() + 1, - fe_eval.get_cell_ids().end(), - [&](const auto &v) { - return v == fe_eval.get_cell_ids()[0] || - v == numbers::invalid_unsigned_int; - }); - - if (use_vectorization == false) - { - for (unsigned int v = 0; v < Number::size(); ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - if (fe_eval.get_face_orientation(v) != 0) - adjust_for_face_orientation_per_lane( - dim, - n_components, - v, - integration_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(v), 0), - true, - shape_info.n_q_points_face, - &temp[0][0], - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - } - } - else if (fe_eval.get_face_orientation() != 0) - adjust_for_face_orientation( - dim, - n_components, - integration_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(), 0), - true, - shape_info.n_q_points_face, - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - - const unsigned int n_q_points_1d_actual = - fe_degree > -1 ? n_q_points_1d : 0; - const unsigned int subface_index = fe_eval.get_subface_index(); - - if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) - { - FEFaceEvaluationImplRaviartThomas:: - template evaluate_or_integrate_in_face( - integration_flag, - fe_eval.get_shape_info().data, - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - scratch_data, - subface_index, - fe_eval.get_face_no() / 2); - } - else if (fe_degree > -1 && - fe_eval.get_subface_index() >= - GeometryInfo::max_children_per_cell && - shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d_actual, - Number>::integrate_in_face(n_components, - integration_flag, - shape_data, - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - else - FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d_actual, - Number>::integrate_in_face(n_components, - integration_flag, - shape_data, - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - - if (use_vectorization == false) - { - for (unsigned int v = 0; v < Number::size(); ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - integration_flag, - shape_info, - values_dofs, - scratch_data, - fe_eval.get_face_no(v)); - - for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; - ++i) - temp[i][v] = scratch_data[i][v]; - } - } - else - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - integration_flag, - shape_info, - temp, - values_dofs, - fe_eval.get_face_no()); - return false; - } - }; - - - - template - void - fe_face_evaluation_process_and_io( - Processor &proc, - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - typename Processor::Number2_ *global_vector_ptr, - const std::vector> *sm_ptr, - const EvaluationData &fe_eval, - typename Processor::VectorizedArrayType_ *temp1) - { - constexpr int dim = Processor::dim_; - constexpr int fe_degree = Processor::fe_degree_; - using VectorizedArrayType = typename Processor::VectorizedArrayType_; - constexpr int n_lanes = VectorizedArrayType::size(); - - using Number = typename Processor::Number_; - using Number2_ = typename Processor::Number2_; - - const auto &shape_data = fe_eval.get_shape_info().data.front(); - constexpr bool integrate = Processor::do_integrate; - const unsigned int face_no = fe_eval.get_face_no(); - const auto &dof_info = fe_eval.get_dof_info(); - const unsigned int cell = fe_eval.get_cell_or_face_batch_id(); - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index = - fe_eval.get_dof_access_index(); - AssertIndexRange(cell, - dof_info.index_storage_variants[dof_access_index].size()); - constexpr unsigned int dofs_per_face = - Utilities::pow(fe_degree + 1, dim - 1); - const unsigned int subface_index = fe_eval.get_subface_index(); - - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - bool all_faces_are_same = n_filled_lanes == n_lanes; - if (n_face_orientations == n_lanes) - for (unsigned int v = 1; v < n_lanes; ++v) - if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) || - fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0)) - { - all_faces_are_same = false; - break; - } - - // check for re-orientation ... - std::array orientation = {}; - - if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same && - fe_eval.is_interior_face() == 0) - for (unsigned int v = 0; v < n_lanes; ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - if (shape_data.nodal_at_cell_boundaries && - fe_eval.get_face_orientation(v) != 0) - { - // ... and in case we detect a re-orientation, go to the other - // version of this function that actually allows for this - if (subface_index == GeometryInfo::max_children_per_cell && - check_face_orientations == false) - { - fe_face_evaluation_process_and_io(proc, - n_components, - evaluation_flag, - global_vector_ptr, - sm_ptr, - fe_eval, - temp1); - return; - } - orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs( - fe_eval.get_face_orientation(v), 0); - } - } - else if (dim == 3 && fe_eval.get_face_orientation() != 0) - { - // go to the other version of this function - if (subface_index == GeometryInfo::max_children_per_cell && - check_face_orientations == false) - { - fe_face_evaluation_process_and_io(proc, - n_components, - evaluation_flag, - global_vector_ptr, - sm_ptr, - fe_eval, - temp1); - return; - } - for (unsigned int v = 0; v < n_face_orientations; ++v) - orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs( - fe_eval.get_face_orientation(), 0); - } - - // we know that the gradient weights for the Hermite case on the - // right (side==1) are the negative from the value at the left - // (side==0), so we only read out one of them. - VectorizedArrayType grad_weight = - shape_data - .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) : - (1 + face_no % 2))]; - - // face_to_cell_index_hermite - std::array index_array_hermite = - {}; - if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) - { - if (n_face_orientations == 1) - index_array_hermite[0] = - &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0); - else - { - for (unsigned int v = 0; v < n_lanes; ++v) - { - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - const auto face_no = fe_eval.get_face_no(v); - - grad_weight[v] = - shape_data.shape_data_on_face[0][fe_degree + - (integrate ? - (2 - (face_no % 2)) : - (1 + (face_no % 2)))][0]; - - index_array_hermite[v] = - &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, - 0); - } - } - } - - // face_to_cell_index_nodal - std::array index_array_nodal = - {}; - if (shape_data.nodal_at_cell_boundaries == true) - { - if (n_face_orientations == 1) - index_array_nodal[0] = - &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0); - else - { - for (unsigned int v = 0; v < n_lanes; ++v) - { - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - const auto face_no = fe_eval.get_face_no(v); - - index_array_nodal[v] = - &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, - 0); - } - } - } - - - const auto reorientate = [&](const unsigned int v, const unsigned int i) { - return (!check_face_orientations || orientation[v] == nullptr) ? - i : - orientation[v][i]; - }; - - const unsigned int cell_index = - dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ? - fe_eval.get_cell_ids()[0] : - cell * n_lanes; - const unsigned int *dof_indices = - &dof_info.dof_indices_contiguous[dof_access_index][cell_index]; - - for (unsigned int comp = 0; comp < n_components; ++comp) - { - const std::size_t index_offset = - dof_info.component_dof_indices_offset - [fe_eval.get_active_fe_index()] - [fe_eval.get_first_selected_component()] + - comp * Utilities::pow(fe_degree + 1, dim); - - // case 1: contiguous and interleaved indices - if (n_face_orientations == 1 && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous) - { - AssertDimension( - dof_info.n_vectorization_lanes_filled[dof_access_index][cell], - n_lanes); - Number2_ *vector_ptr = - global_vector_ptr + dof_indices[0] + index_offset * n_lanes; - - if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - Assert(n_face_orientations == 1, ExcNotImplemented()); - - const unsigned int ind1 = index_array_hermite[0][2 * i]; - const unsigned int ind2 = index_array_hermite[0][2 * i + 1]; - const unsigned int i_ = reorientate(0, i); - proc.hermite_grad_vectorized(temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1 * n_lanes, - vector_ptr + ind2 * n_lanes, - grad_weight); - } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - Assert(n_face_orientations == 1, ExcNotImplemented()); - - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = index_array_nodal[0][i]; - proc.value_vectorized(temp1[i_], - vector_ptr + ind * n_lanes); - } - } - } - - // case 2: contiguous and interleaved indices with fixed stride - else if (n_face_orientations == 1 && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_strided) - { - AssertDimension( - dof_info.n_vectorization_lanes_filled[dof_access_index][cell], - n_lanes); - Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes; - if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - Assert(n_face_orientations == 1, ExcNotImplemented()); - - const unsigned int i_ = reorientate(0, i); - const unsigned int ind1 = - index_array_hermite[0][2 * i] * n_lanes; - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1] * n_lanes; - proc.hermite_grad_vectorized_indexed( - temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1, - vector_ptr + ind2, - grad_weight, - dof_indices, - dof_indices); - } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - Assert(n_face_orientations == 1, ExcNotImplemented()); - - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = index_array_nodal[0][i] * n_lanes; - proc.value_vectorized_indexed(temp1[i_], - vector_ptr + ind, - dof_indices); - } - } - } - - // case 3: contiguous and interleaved indices with mixed stride - else if (n_face_orientations == 1 && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_mixed_strides) - { - const unsigned int *strides = - &dof_info.dof_indices_interleave_strides[dof_access_index] - [cell * n_lanes]; - unsigned int indices[n_lanes]; - for (unsigned int v = 0; v < n_lanes; ++v) - indices[v] = dof_indices[v] + index_offset * strides[v]; - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) - { - if (n_filled_lanes == n_lanes) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - Assert(n_face_orientations == 1, ExcNotImplemented()); - - const unsigned int i_ = reorientate(0, i); - unsigned int ind1[n_lanes]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; v < n_lanes; ++v) - ind1[v] = indices[v] + - index_array_hermite[0][2 * i] * strides[v]; - unsigned int ind2[n_lanes]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; v < n_lanes; ++v) - ind2[v] = - indices[v] + - // TODO - index_array_hermite[0][2 * i + 1] * strides[v]; - proc.hermite_grad_vectorized_indexed( - temp1[i_], - temp1[i_ + dofs_per_face], - global_vector_ptr, - global_vector_ptr, - grad_weight, - ind1, - ind2); - } - else - { - if (integrate == false) - for (unsigned int i = 0; i < 2 * dofs_per_face; ++i) - temp1[i] = VectorizedArrayType(); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int i_ = - reorientate(n_face_orientations == 1 ? 0 : v, i); - proc.hermite_grad( - temp1[i_][v], - temp1[i_ + dofs_per_face][v], - global_vector_ptr - [indices[v] + - index_array_hermite - [n_face_orientations == 1 ? 0 : v][2 * i] * - strides[v]], - global_vector_ptr - [indices[v] + - index_array_hermite[n_face_orientations == 1 ? - 0 : - v][2 * i + 1] * - strides[v]], - grad_weight[n_face_orientations == 1 ? 0 : v]); - } - } - } - else - { - if (n_filled_lanes == n_lanes) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - Assert(n_face_orientations == 1, ExcInternalError()); - unsigned int ind[n_lanes]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; v < n_lanes; ++v) - ind[v] = - indices[v] + index_array_nodal[0][i] * strides[v]; - const unsigned int i_ = reorientate(0, i); - proc.value_vectorized_indexed(temp1[i_], - global_vector_ptr, - ind); - } - else - { - if (integrate == false) - for (unsigned int i = 0; i < dofs_per_face; ++i) - temp1[i] = VectorizedArrayType(); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int i = 0; i < dofs_per_face; ++i) - proc.value( - temp1[reorientate(n_face_orientations == 1 ? 0 : v, - i)][v], - global_vector_ptr - [indices[v] + - index_array_nodal[n_face_orientations == 1 ? 0 : v] - [i] * - strides[v]]); - } - } - } - - // case 4: contiguous indices without interleaving - else if (n_face_orientations > 1 || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - contiguous) - { - Number2_ *vector_ptr = global_vector_ptr + index_offset; - - const bool vectorization_possible = - all_faces_are_same && (sm_ptr == nullptr); - - std::array vector_ptrs; - std::array reordered_indices; - - if (vectorization_possible == false) - { - vector_ptrs = {}; - if (n_face_orientations == 1) - { - for (unsigned int v = 0; v < n_filled_lanes; ++v) - if (sm_ptr == nullptr) - { - vector_ptrs[v] = vector_ptr + dof_indices[v]; - } - else - { - const auto &temp = - dof_info - .dof_indices_contiguous_sm[dof_access_index] - [cell * n_lanes + v]; - vector_ptrs[v] = const_cast( - sm_ptr->operator[](temp.first).data() + - temp.second + index_offset); - } - } - else if (n_face_orientations == n_lanes) - { - const auto &cells = fe_eval.get_cell_ids(); - for (unsigned int v = 0; v < n_lanes; ++v) - if (cells[v] != numbers::invalid_unsigned_int) - { - if (sm_ptr == nullptr) - { - vector_ptrs[v] = - vector_ptr + - dof_info - .dof_indices_contiguous[dof_access_index] - [cells[v]]; - } - else - { - const auto &temp = - dof_info - .dof_indices_contiguous_sm[dof_access_index] - [cells[v]]; - vector_ptrs[v] = const_cast( - sm_ptr->operator[](temp.first).data() + - temp.second + index_offset); - } - } - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (n_face_orientations == n_lanes) - { - for (unsigned int v = 0; v < n_lanes; ++v) - reordered_indices[v] = - dof_info.dof_indices_contiguous[dof_access_index] - [fe_eval.get_cell_ids()[v]]; - dof_indices = reordered_indices.data(); - } - - if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) - { - if (vectorization_possible) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = index_array_hermite[0][2 * i]; - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1]; - const unsigned int i_ = reorientate(0, i); - - proc.hermite_grad_vectorized_indexed( - temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1, - vector_ptr + ind2, - grad_weight, - dof_indices, - dof_indices); - } - else if (n_face_orientations == 1) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = index_array_hermite[0][2 * i]; - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1]; - const unsigned int i_ = reorientate(0, i); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - proc.hermite_grad(temp1[i_][v], - temp1[i_ + dofs_per_face][v], - vector_ptrs[v][ind1], - vector_ptrs[v][ind2], - grad_weight[v]); - - if (integrate == false) - for (unsigned int v = n_filled_lanes; v < n_lanes; ++v) - { - temp1[i][v] = 0.0; - temp1[i + dofs_per_face][v] = 0.0; - } - } - else - { - if (integrate == false && n_filled_lanes < n_lanes) - for (unsigned int i = 0; i < dofs_per_face; ++i) - temp1[i] = temp1[i + dofs_per_face] = Number(); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int i = 0; i < dofs_per_face; ++i) - proc.hermite_grad( - temp1[reorientate(v, i)][v], - temp1[reorientate(v, i) + dofs_per_face][v], - vector_ptrs[v][index_array_hermite[v][2 * i]], - vector_ptrs[v][index_array_hermite[v][2 * i + 1]], - grad_weight[v]); - } - } - else - { - if (vectorization_possible) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind = index_array_nodal[0][i]; - const unsigned int i_ = reorientate(0, i); - - proc.value_vectorized_indexed(temp1[i_], - vector_ptr + ind, - dof_indices); - } - else if (n_face_orientations == 1) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind = index_array_nodal[0][i]; - const unsigned int i_ = reorientate(0, i); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - proc.value(temp1[i_][v], vector_ptrs[v][ind]); - - if (integrate == false) - for (unsigned int v = n_filled_lanes; v < n_lanes; ++v) - temp1[i_][v] = 0.0; - } - else - { - if (integrate == false && n_filled_lanes < n_lanes) - for (unsigned int i = 0; i < dofs_per_face; ++i) - temp1[i] = Number(); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int i = 0; i < dofs_per_face; ++i) - proc.value(temp1[reorientate(v, i)][v], - vector_ptrs[v][index_array_nodal[v][i]]); - } - } - } - else - { - // We should not end up here, this should be caught by - // FEFaceEvaluationImplGatherEvaluateSelector::supports() - Assert(false, ExcInternalError()); - } - temp1 += 3 * dofs_per_face; - } - } - - - - template - struct FEFaceEvaluationImplGatherEvaluateSelector - { - using Number = typename VectorizedArrayType::value_type; - - template - static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number2 *src_ptr, - const std::vector> *sm_ptr, - FEEvaluationData &fe_eval) - { - Assert(fe_degree > -1, ExcInternalError()); - Assert(fe_eval.get_shape_info().element_type <= - MatrixFreeFunctions::tensor_symmetric_no_collocation, - ExcInternalError()); - - const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1); - - VectorizedArrayType *temp = fe_eval.get_scratch_data().begin(); - VectorizedArrayType *scratch_data = - temp + 3 * n_components * dofs_per_face; - - Processor p; - - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) - fe_face_evaluation_process_and_io( - p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp); - else - fe_face_evaluation_process_and_io<1>( - p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp); - - const unsigned int subface_index = fe_eval.get_subface_index(); - - if (subface_index >= GeometryInfo::max_children_per_cell) - FEFaceEvaluationImpl:: - evaluate_in_face(n_components, - evaluation_flag, - fe_eval.get_shape_info().data.front(), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - else - FEFaceEvaluationImpl:: - evaluate_in_face(n_components, - evaluation_flag, - fe_eval.get_shape_info().data.front(), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - - // re-orientation for cases not possible with above algorithm - if (subface_index < GeometryInfo::max_children_per_cell) - { - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) - { - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == - numbers::invalid_unsigned_int) - continue; - - if (fe_eval.get_face_orientation(v) != 0) - adjust_for_face_orientation_per_lane( - dim, - n_components, - v, - evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(v), 0), - false, - Utilities::pow(n_q_points_1d, dim - 1), - &temp[0][0], - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - } - } - else if (fe_eval.get_face_orientation() != 0) - adjust_for_face_orientation( - dim, - n_components, - evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(), 0), - false, - Utilities::pow(n_q_points_1d, dim - 1), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - } - - return false; - } - - template - static bool - supports(const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number2 *vector_ptr, - MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage) - { - const unsigned int fe_degree = shape_info.data.front().fe_degree; - if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries || - (evaluation_flag & EvaluationFlags::gradients && - (fe_degree < 2 || - shape_info.data.front().element_type != - MatrixFreeFunctions::tensor_symmetric_hermite)) || - (evaluation_flag & EvaluationFlags::hessians) || - vector_ptr == nullptr || - shape_info.data.front().element_type > - MatrixFreeFunctions::tensor_symmetric_no_collocation || - storage < - MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) - return false; - else - return true; - } - - private: - template - struct Processor - { - static const bool do_integrate = false; - static const int dim_ = dim; - static const int fe_degree_ = fe_degree; - using VectorizedArrayType_ = VectorizedArrayType; - using Number_ = Number; - using Number2_ = const Number2; - - template - void - hermite_grad_vectorized(T0 &temp_1, - T0 &temp_2, - const T1 src_ptr_1, - const T1 src_ptr_2, - const T2 &grad_weight) - { - do_vectorized_read(src_ptr_1, temp_1); - do_vectorized_read(src_ptr_2, temp_2); - temp_2 = grad_weight * (temp_1 - temp_2); - } - - template - void - value_vectorized(T1 &temp, const T2 src_ptr) - { - do_vectorized_read(src_ptr, temp); - } - - template - void - hermite_grad_vectorized_indexed(T0 &temp_1, - T0 &temp_2, - const T1 src_ptr_1, - const T1 src_ptr_2, - const T2 &grad_weight, - const T3 &indices_1, - const T3 &indices_2) - { - do_vectorized_gather(src_ptr_1, indices_1, temp_1); - do_vectorized_gather(src_ptr_2, indices_2, temp_2); - temp_2 = grad_weight * (temp_1 - temp_2); - } - - template - void - value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices) - { - do_vectorized_gather(src_ptr, indices, temp); - } - - template - void - hermite_grad(T0 &temp_1, - T0 &temp_2, - const T1 &src_ptr_1, - const T1 &src_ptr_2, - const T2 &grad_weight) - { - // case 3a) - temp_1 = src_ptr_1; - temp_2 = grad_weight * (temp_1 - src_ptr_2); - } - - template - void - value(T1 &temp, const T2 &src_ptr) - { - // case 3b) - temp = src_ptr; - } - }; - }; - - - - template - struct FEFaceEvaluationImplIntegrateScatterSelector - { - using Number = typename VectorizedArrayType::value_type; - - template - static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number2 *dst_ptr, - const std::vector> *sm_ptr, - FEEvaluationData &fe_eval) - { - Assert(fe_degree > -1, ExcInternalError()); - Assert(fe_eval.get_shape_info().element_type <= - MatrixFreeFunctions::tensor_symmetric_no_collocation, - ExcInternalError()); - - const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1); - - VectorizedArrayType *temp = fe_eval.get_scratch_data().begin(); - VectorizedArrayType *scratch_data = - temp + 3 * n_components * dofs_per_face; - - const unsigned int subface_index = fe_eval.get_subface_index(); - - // re-orientation for cases not possible with the io function below - if (subface_index < GeometryInfo::max_children_per_cell) - { - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - { - // the loop breaks once an invalid_unsigned_int is hit for - // all cases except the exterior faces in the ECL loop (where - // some faces might be at the boundaries but others not) - if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) - continue; - - if (fe_eval.get_face_orientation(v) != 0) - adjust_for_face_orientation_per_lane( - dim, - n_components, - v, - integration_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(v), 0), - true, - Utilities::pow(n_q_points_1d, dim - 1), - &temp[0][0], - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - } - else if (fe_eval.get_face_orientation() != 0) - adjust_for_face_orientation( - dim, - n_components, - integration_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(), 0), - true, - Utilities::pow(n_q_points_1d, dim - 1), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); - } - - if (fe_degree > -1 && fe_eval.get_subface_index() >= - GeometryInfo::max_children_per_cell) - FEFaceEvaluationImpl:: - integrate_in_face(n_components, - integration_flag, - fe_eval.get_shape_info().data.front(), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - else - FEFaceEvaluationImpl:: - integrate_in_face(n_components, - integration_flag, - fe_eval.get_shape_info().data.front(), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians(), - scratch_data, - subface_index); - - Processor p; - - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) - fe_face_evaluation_process_and_io( - p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp); - else - fe_face_evaluation_process_and_io<1>( - p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp); - - return false; - } - - private: - template - struct Processor - { - static const bool do_integrate = true; - static const int dim_ = dim; - static const int fe_degree_ = fe_degree; - using VectorizedArrayType_ = VectorizedArrayType; - using Number_ = Number; - using Number2_ = Number2; - - template - void - hermite_grad_vectorized(const T0 &temp_1, - const T1 &temp_2, - T2 dst_ptr_1, - T3 dst_ptr_2, - const T4 &grad_weight) - { - // case 1a) - const VectorizedArrayType val = temp_1 - grad_weight * temp_2; - const VectorizedArrayType grad = grad_weight * temp_2; - do_vectorized_add(val, dst_ptr_1); - do_vectorized_add(grad, dst_ptr_2); - } - - template - void - value_vectorized(const T0 &temp, T1 dst_ptr) - { - // case 1b) - do_vectorized_add(temp, dst_ptr); - } - - template - void - hermite_grad_vectorized_indexed(const T0 &temp_1, - const T0 &temp_2, - T1 dst_ptr_1, - T1 dst_ptr_2, - const T2 &grad_weight, - const T3 &indices_1, - const T3 &indices_2) - { - // case 2a) - const VectorizedArrayType val = temp_1 - grad_weight * temp_2; - const VectorizedArrayType grad = grad_weight * temp_2; - do_vectorized_scatter_add(val, indices_1, dst_ptr_1); - do_vectorized_scatter_add(grad, indices_2, dst_ptr_2); - } - - template - void - value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices) - { - // case 2b) - do_vectorized_scatter_add(temp, indices, dst_ptr); - } - - template - void - hermite_grad(const T0 &temp_1, - const T0 &temp_2, - T1 &dst_ptr_1, - T1 &dst_ptr_2, - const T2 &grad_weight) - { - // case 3a) - const Number val = temp_1 - grad_weight * temp_2; - const Number grad = grad_weight * temp_2; - dst_ptr_1 += val; - dst_ptr_2 += grad; - } - - template - void - value(const T0 &temp, T1 &dst_ptr) - { - // case 3b) - dst_ptr += temp; - } - }; - }; - - - - /** - * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation, - * using an FEEvaluationData argument. - */ - template - struct CellwiseInverseMassMatrixImplBasic - { - using Number2 = - typename FEEvaluationData::shape_info_number_type; - - template - static bool - run(const unsigned int n_components, - const FEEvaluationData &fe_eval, - const Number *in_array, - Number *out_array) - { - const unsigned int given_degree = - (fe_degree > -1) ? fe_degree : - fe_eval.get_shape_info().data.front().fe_degree; - - const unsigned int dofs_per_component = - Utilities::pow(given_degree + 1, dim); - - Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); - Assert(fe_eval.get_shape_info().element_type <= - MatrixFreeFunctions::tensor_symmetric_no_collocation, - ExcNotImplemented()); - - EvaluatorTensorProduct - evaluator({}, - {}, - fe_eval.get_shape_info().data.front().inverse_shape_values_eo, - given_degree + 1, - given_degree + 1); - - for (unsigned int d = 0; d < n_components; ++d) - { - const Number *in = in_array + d * dofs_per_component; - Number *out = out_array + d * dofs_per_component; - // Need to select 'apply' method with hessian slot because values - // assume symmetries that do not exist in the inverse shapes - evaluator.template hessians<0, true, false>(in, out); - if (dim > 1) - evaluator.template hessians<1, true, false>(out, out); - if (dim > 2) - evaluator.template hessians<2, true, false>(out, out); - } - for (unsigned int q = 0; q < dofs_per_component; ++q) - { - const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q); - for (unsigned int d = 0; d < n_components; ++d) - out_array[q + d * dofs_per_component] *= inverse_JxW_q; - } - for (unsigned int d = 0; d < n_components; ++d) - { - Number *out = out_array + d * dofs_per_component; - if (dim > 2) - evaluator.template hessians<2, false, false>(out, out); - if (dim > 1) - evaluator.template hessians<1, false, false>(out, out); - evaluator.template hessians<0, false, false>(out, out); - } - return false; - } - }; - - - - /** - * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation - * with user-provided coefficients at quadrature points (in contrast to - * CellwiseInverseMassMatrixImplBasic, which implicitly uses `1/(|J|xW)' as - * coefficient). - */ - template - struct CellwiseInverseMassMatrixImplFlexible - { - using Number2 = - typename FEEvaluationData::shape_info_number_type; - - template - static bool - run(const unsigned int n_desired_components, - const FEEvaluationData &fe_eval, - const ArrayView &inverse_coefficients, - const bool dyadic_coefficients, - const Number *in_array, - Number *out_array) - { - const unsigned int given_degree = - (fe_degree > -1) ? fe_degree : - fe_eval.get_shape_info().data.front().fe_degree; - - const unsigned int dofs_per_component = - Utilities::pow(given_degree + 1, dim); - - Assert(inverse_coefficients.size() > 0 && - inverse_coefficients.size() % dofs_per_component == 0, - ExcMessage( - "Expected diagonal to be a multiple of scalar dof per cells")); - - if (!dyadic_coefficients) - { - if (inverse_coefficients.size() != dofs_per_component) - AssertDimension(n_desired_components * dofs_per_component, - inverse_coefficients.size()); - } - else - { - AssertDimension(n_desired_components * n_desired_components * - dofs_per_component, - inverse_coefficients.size()); - } - - Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); - Assert(fe_eval.get_shape_info().element_type <= - MatrixFreeFunctions::tensor_symmetric_no_collocation, - ExcNotImplemented()); - - EvaluatorTensorProduct - evaluator({}, - {}, - fe_eval.get_shape_info().data.front().inverse_shape_values_eo, - given_degree + 1, - given_degree + 1); - - const Number *in = in_array; - Number *out = out_array; - - const Number *inv_coefficient = inverse_coefficients.data(); - - const unsigned int shift_coefficient = - inverse_coefficients.size() > dofs_per_component ? dofs_per_component : - 0; - - const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components; - const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1; - - for (unsigned int d = 0; d < n_comp_outer; ++d) - { - for (unsigned int di = 0; di < n_comp_inner; ++di) - { - const Number *in_ = in + di * dofs_per_component; - Number *out_ = out + di * dofs_per_component; - evaluator.template hessians<0, true, false>(in_, out_); - if (dim > 1) - evaluator.template hessians<1, true, false>(out_, out_); - if (dim > 2) - evaluator.template hessians<2, true, false>(out_, out_); - } - if (dyadic_coefficients) - { - const auto n_coeff_components = - n_desired_components * n_desired_components; - if (n_desired_components == dim) - { - for (unsigned int q = 0; q < dofs_per_component; ++q) - vmult(&inv_coefficient[q * n_coeff_components], - &in[q], - &out[q], - dofs_per_component); - } - else - { - for (unsigned int q = 0; q < dofs_per_component; ++q) - vmult<-1>(&inv_coefficient[q * n_coeff_components], - &in[q], - &out[q], - dofs_per_component, - n_desired_components); + for (unsigned int q = 0; q < dofs_per_component; ++q) + vmult<-1>(&inv_coefficient[q * n_coeff_components], + &in[q], + &out[q], + dofs_per_component, + n_desired_components); } } else @@ -5567,21 +2518,6 @@ namespace internal return false; } }; - - /** - * This struct is used to implement - * FEEvaluation::fast_evaluation_supported() and - * FEFaceEvaluation::fast_evaluation_supported(). - */ - struct FastEvaluationSupported - { - template - static bool - run() - { - return fe_degree != -1; - } - }; } // end of namespace internal diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h new file mode 100644 index 0000000000..09a01d423c --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_kernels_face.h @@ -0,0 +1,2920 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_evaluation_kernels_face_h +#define dealii_matrix_free_evaluation_kernels_face_h + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + template + struct FEFaceEvaluationImpl + { + // We enable a transformation to collocation for derivatives if it gives + // correct results (first two conditions), if it is the most efficient + // choice in terms of operation counts (third condition) and if we were + // able to initialize the fields in shape_info.templates.h from the + // polynomials (fourth condition). + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + using Eval = EvaluatorTensorProduct; + + static Eval + create_evaluator_tensor_product( + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int subface_index, + const unsigned int direction) + { + if (symmetric_evaluate) + return Eval(data.shape_values_eo, + data.shape_gradients_eo, + data.shape_hessians_eo, + data.fe_degree + 1, + data.n_q_points_1d); + else if (subface_index >= GeometryInfo::max_children_per_cell) + return Eval(data.shape_values, + data.shape_gradients, + data.shape_hessians, + data.fe_degree + 1, + data.n_q_points_1d); + else + { + const unsigned int index = + direction == 0 ? subface_index % 2 : subface_index / 2; + return Eval(data.values_within_subface[index], + data.gradients_within_subface[index], + data.hessians_within_subface[index], + data.fe_degree + 1, + data.n_q_points_1d); + } + } + + static void + evaluate_in_face( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::UnivariateShapeData &data, + Number *values_dofs, + Number *values_quad, + Number *gradients_quad, + Number *hessians_quad, + Number *scratch_data, + const unsigned int subface_index) + { + Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); + Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); + + const std::size_t n_dofs = fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(data.fe_degree + 1, dim - 1); + const std::size_t n_q_points = + fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : + Utilities::pow(data.n_q_points_1d, dim - 1); + + // keep a copy of the original pointer for the case of the Hessians + Number *values_dofs_ptr = values_dofs; + + if ((evaluation_flag & EvaluationFlags::values) != 0u && + ((evaluation_flag & EvaluationFlags::gradients) == 0u)) + for (unsigned int c = 0; c < n_components; ++c) + { + switch (dim) + { + case 3: + 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>(values_dofs, + values_quad); + break; + case 1: + values_quad[0] = values_dofs[0]; + break; + default: + Assert(false, ExcNotImplemented()); + } + // Note: we always keep storage of values, 1st and 2nd derivatives + // in an array + values_dofs += 3 * n_dofs; + values_quad += n_q_points; + } + else if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + for (unsigned int c = 0; c < n_components; ++c) + { + switch (dim) + { + case 3: + if (symmetric_evaluate && + use_collocation_evaluation(fe_degree, n_q_points_1d)) + { + eval0.template values<0, true, false>(values_dofs, + values_quad); + eval0.template values<1, true, false>(values_quad, + values_quad); + EvaluatorTensorProduct + eval_grad({}, data.shape_gradients_collocation_eo, {}); + eval_grad.template gradients<0, true, false, 3>( + values_quad, gradients_quad); + eval_grad.template gradients<1, true, false, 3>( + values_quad, gradients_quad + 1); + } + else + { + // grad x + eval0.template gradients<0, true, false>(values_dofs, + scratch_data); + eval1.template values<1, true, false, 3>(scratch_data, + gradients_quad); + + // grad y + eval0.template values<0, true, false>(values_dofs, + scratch_data); + eval1.template gradients<1, true, false, 3>( + scratch_data, gradients_quad + 1); + + if ((evaluation_flag & EvaluationFlags::values) != 0u) + eval1.template values<1, true, false>(scratch_data, + values_quad); + } + // grad z + eval0.template values<0, true, false>(values_dofs + n_dofs, + scratch_data); + eval1.template values<1, true, false, 3>(scratch_data, + gradients_quad + 2); + + break; + case 2: + eval0.template values<0, true, false, 2>(values_dofs + n_dofs, + gradients_quad + 1); + eval0.template gradients<0, true, false, 2>(values_dofs, + gradients_quad); + if ((evaluation_flag & EvaluationFlags::values) != 0u) + eval0.template values<0, true, false>(values_dofs, + values_quad); + break; + case 1: + values_quad[0] = values_dofs[0]; + gradients_quad[0] = values_dofs[1]; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + values_dofs += 3 * n_dofs; + values_quad += n_q_points; + gradients_quad += dim * n_q_points; + } + + if ((evaluation_flag & EvaluationFlags::hessians) != 0u) + { + values_dofs = values_dofs_ptr; + for (unsigned int c = 0; c < n_components; ++c) + { + switch (dim) + { + case 3: + // grad xx + 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>(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>(values_dofs + + 2 * n_dofs, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + hessians_quad + + 2 * n_q_points); + + // grad xy + 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>(values_dofs + + n_dofs, + scratch_data); + eval1.template values<1, true, false>(scratch_data, + hessians_quad + + 4 * n_q_points); + + // grad yz + eval0.template values<0, true, false>(values_dofs + n_dofs, + 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>(values_dofs, + hessians_quad); + // grad yy + eval0.template values<0, true, false>( + values_dofs + 2 * n_dofs, hessians_quad + n_q_points); + // grad xy + eval0.template gradients<0, true, false>( + values_dofs + n_dofs, hessians_quad + 2 * n_q_points); + break; + case 1: + hessians_quad[0] = values_dofs[2]; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + values_dofs += 3 * n_dofs; + hessians_quad += dim * (dim + 1) / 2 * n_q_points; + } + } + } + + static void + integrate_in_face( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const MatrixFreeFunctions::UnivariateShapeData &data, + Number *values_dofs, + Number *values_quad, + Number *gradients_quad, + Number *hessians_quad, + Number *scratch_data, + const unsigned int subface_index) + { + Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); + Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); + + const std::size_t n_dofs = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + (dim > 1 ? Utilities::fixed_power(data.fe_degree + 1) : 1); + const std::size_t n_q_points = + fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : + Utilities::pow(data.n_q_points_1d, dim - 1); + + // keep a copy of the original pointer for the case of the Hessians + Number *values_dofs_ptr = values_dofs; + + if ((integration_flag & EvaluationFlags::values) != 0u && + (integration_flag & EvaluationFlags::gradients) == 0u) + for (unsigned int c = 0; c < n_components; ++c) + { + switch (dim) + { + case 3: + 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>(values_quad, + values_dofs); + break; + case 1: + values_dofs[0] = values_quad[0]; + break; + default: + Assert(false, ExcNotImplemented()); + } + values_dofs += 3 * n_dofs; + values_quad += n_q_points; + } + else if ((integration_flag & EvaluationFlags::gradients) != 0u) + for (unsigned int c = 0; c < n_components; ++c) + { + switch (dim) + { + case 3: + // grad z + eval1.template values<1, false, false, 3>(gradients_quad + 2, + scratch_data); + eval0.template values<0, false, false>(scratch_data, + values_dofs + n_dofs); + if (symmetric_evaluate && + use_collocation_evaluation(fe_degree, n_q_points_1d)) + { + EvaluatorTensorProduct + eval_grad({}, data.shape_gradients_collocation_eo, {}); + if ((integration_flag & EvaluationFlags::values) != 0u) + eval_grad.template gradients<1, false, true, 3>( + gradients_quad + 1, values_quad); + else + eval_grad.template gradients<1, false, false, 3>( + gradients_quad + 1, values_quad); + eval_grad.template gradients<0, false, true, 3>( + gradients_quad, values_quad); + eval0.template values<1, false, false>(values_quad, + values_quad); + eval0.template values<0, false, false>(values_quad, + values_dofs); + } + else + { + if ((integration_flag & EvaluationFlags::values) != 0u) + { + eval1.template values<1, false, false>(values_quad, + scratch_data); + eval1.template gradients<1, false, true, 3>( + gradients_quad + 1, scratch_data); + } + else + eval1.template gradients<1, false, false, 3>( + gradients_quad + 1, scratch_data); + + // grad y + eval0.template values<0, false, false>(scratch_data, + values_dofs); + + // grad x + eval1.template values<1, false, false, 3>(gradients_quad, + scratch_data); + eval0.template gradients<0, false, true>(scratch_data, + values_dofs); + } + break; + case 2: + eval0.template values<0, false, false, 2>(gradients_quad + 1, + values_dofs + + n_dofs); + eval0.template gradients<0, false, false, 2>(gradients_quad, + values_dofs); + if ((integration_flag & EvaluationFlags::values) != 0u) + eval0.template values<0, false, true>(values_quad, + values_dofs); + break; + case 1: + values_dofs[0] = values_quad[0]; + values_dofs[1] = gradients_quad[0]; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + values_dofs += 3 * n_dofs; + values_quad += n_q_points; + gradients_quad += dim * n_q_points; + } + + if ((integration_flag & EvaluationFlags::hessians) != 0u) + { + values_dofs = values_dofs_ptr; + for (unsigned int c = 0; c < n_components; ++c) + { + switch (dim) + { + case 3: + // grad xx + eval1.template values<1, false, false>(hessians_quad, + scratch_data); + if ((integration_flag & (EvaluationFlags::values | + EvaluationFlags::gradients)) != 0u) + eval0.template hessians<0, false, true>(scratch_data, + values_dofs); + else + eval0.template hessians<0, false, false>(scratch_data, + values_dofs); + + // grad yy + 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>(hessians_quad + + 2 * n_q_points, + scratch_data); + eval0.template values<0, false, false>(scratch_data, + values_dofs + + 2 * n_dofs); + + // grad xy + 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>(hessians_quad + + 4 * n_q_points, + scratch_data); + if ((integration_flag & EvaluationFlags::gradients) != 0u) + eval0.template gradients<0, false, true>(scratch_data, + values_dofs + + n_dofs); + else + eval0.template gradients<0, false, false>(scratch_data, + values_dofs + + n_dofs); + + // grad yz + eval1.template gradients<1, false, false>(hessians_quad + + 5 * n_q_points, + scratch_data); + eval0.template values<0, false, true>(scratch_data, + values_dofs + n_dofs); + + break; + case 2: + // grad xx + if ((integration_flag & (EvaluationFlags::values | + EvaluationFlags::gradients)) != 0u) + eval0.template hessians<0, false, true>(hessians_quad, + values_dofs); + else + eval0.template hessians<0, false, false>(hessians_quad, + values_dofs); + + // grad yy + eval0.template values<0, false, false>( + hessians_quad + n_q_points, values_dofs + 2 * n_dofs); + // grad xy + if ((integration_flag & EvaluationFlags::gradients) != 0u) + eval0.template gradients<0, false, true>( + hessians_quad + 2 * n_q_points, values_dofs + n_dofs); + else + eval0.template gradients<0, false, false>( + hessians_quad + 2 * n_q_points, values_dofs + n_dofs); + break; + case 1: + values_dofs[2] = hessians_quad[0]; + if ((integration_flag & EvaluationFlags::values) == 0u) + values_dofs[0] = 0; + if ((integration_flag & EvaluationFlags::gradients) == 0u) + values_dofs[1] = 0; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + values_dofs += 3 * n_dofs; + hessians_quad += dim * (dim + 1) / 2 * n_q_points; + } + } + } + }; + + + + template + struct FEFaceEvaluationImplRaviartThomas + { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + /** + * Apply the sum factorization kernels within the face for Raviart-Thomas + * elements for either evaluation or integration + */ + template + static inline void + evaluate_or_integrate_in_face( + const EvaluationFlags::EvaluationFlags evaluation_flag, + const std::vector> + &shape_data, + Number *values_dofs_in, + Number *values, + Number *gradients, + Number *scratch_data, + const unsigned int subface_index, + const unsigned int face_direction) + { + AssertDimension(shape_data.size(), 2); + + const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree; + const int n_rows_n = degree + 1; + const int n_rows_t = degree; + 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}}}}; + + (void)scratch_data; + (void)subface_index; + // TODO: This is currently not implemented, but the test + // matrix_vector_rt_face_03 apparently works without it -> check + // if (subface_index < GeometryInfo::max_children_per_cell) + // Assert(false, ExcNotImplemented()); + + using Eval = EvaluatorTensorProduct; + + std::array values_dofs_offsets = {}; + for (unsigned int comp = 0; comp < dim - 1; ++comp) + { + if (dim == 2) + values_dofs_offsets[comp + 1] = + values_dofs_offsets[comp] + + 3 * dofs_per_direction[comp][(face_direction + 1) % dim]; + else + values_dofs_offsets[comp + 1] = + values_dofs_offsets[comp] + + 3 * dofs_per_direction[comp][(face_direction + 1) % dim] * + dofs_per_direction[comp][(face_direction + 2) % dim]; + } + + // Jacobians on faces are reordered to enable simple access with the + // regular evaluators; to get the RT Piola transform right, we need to + // pass through the values_dofs array in a permuted right order + std::array components; + for (unsigned int comp = 0; comp < dim; ++comp) + components[comp] = (face_direction + comp + 1) % dim; + + for (const unsigned int comp : components) + { + Number *values_dofs = values_dofs_in + values_dofs_offsets[comp]; + + std::array n_blocks{ + {dofs_per_direction[comp][(face_direction + 1) % dim], + (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] : + 1)}}; + + if constexpr (dim == 3) + { + EvaluatorTensorProduct + eval_g({}, + shape_data[0].shape_gradients_collocation_eo.data(), + {}); + if (!do_integrate) + { + EvaluatorTensorProductAnisotropic + eval; + // Evaluate in 3d + if (n_blocks[0] == n_rows_n) + { + eval.template normal<0>(shape_data[0], + values_dofs, + values); + eval.template tangential<1, 0>(shape_data[1], + values, + values); + + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template normal<0>(shape_data[0], + values_dofs + + n_blocks[0] * n_blocks[1], + scratch_data); + eval.template tangential<1, 0, dim>(shape_data[1], + scratch_data, + gradients + 2); + } + } + else if (n_blocks[1] == n_rows_n) + { + eval.template normal<1>(shape_data[0], + values_dofs, + values); + eval.template tangential<0, 1>(shape_data[1], + values, + values); + + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template normal<1>(shape_data[0], + values_dofs + + n_blocks[0] * n_blocks[1], + scratch_data); + eval.template tangential<0, 1, dim>(shape_data[1], + scratch_data, + gradients + 2); + } + } + else + { + Eval eval(shape_data[1].shape_values_eo.data(), {}, {}); + eval.template values<0, true, false>(values_dofs, values); + eval.template values<1, true, false>(values, values); + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template values<0, true, false>(values_dofs + + n_blocks[0] * + n_blocks[1], + scratch_data); + eval.template values<1, true, false, dim>( + scratch_data, gradients + 2); + } + } + if (evaluation_flag & EvaluationFlags::gradients) + { + eval_g.template gradients<0, true, false, dim>(values, + gradients); + eval_g.template gradients<1, true, false, dim>(values, + gradients + + 1); + } + } + else + { + EvaluatorTensorProductAnisotropic + eval; + // Integrate in 3d + if (evaluation_flag & EvaluationFlags::gradients) + { + if (evaluation_flag & EvaluationFlags::values) + eval_g.template gradients<0, false, true, dim>( + gradients, values); + else + eval_g.template gradients<0, false, false, dim>( + gradients, values); + eval_g.template gradients<1, false, true, dim>(gradients + + 1, + values); + } + if (n_blocks[0] == n_rows_n) + { + eval.template tangential<1, 0>(shape_data[1], + values, + values); + eval.template normal<0>(shape_data[0], + values, + values_dofs); + + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template tangential<1, 0, dim>(shape_data[1], + gradients + 2, + scratch_data); + eval.template normal<0>(shape_data[0], + scratch_data, + values_dofs + + n_blocks[0] * n_blocks[1]); + } + } + else if (n_blocks[1] == n_rows_n) + { + eval.template tangential<0, 1>(shape_data[1], + values, + values); + eval.template normal<1>(shape_data[0], + values, + values_dofs); + + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template tangential<0, 1, dim>(shape_data[1], + gradients + 2, + scratch_data); + eval.template normal<1>(shape_data[0], + scratch_data, + values_dofs + + n_blocks[0] * n_blocks[1]); + } + } + else + { + Eval eval_iso(shape_data[1].shape_values_eo.data(), + {}, + {}); + eval_iso.template values<1, false, false>(values, values); + eval_iso.template values<0, false, false>(values, + values_dofs); + if (evaluation_flag & EvaluationFlags::gradients) + { + eval_iso.template values<1, false, false, dim>( + gradients + 2, scratch_data); + eval_iso.template values<0, false, false>( + scratch_data, + values_dofs + n_blocks[0] * n_blocks[1]); + } + } + } + } + else + { + using EvalN = EvaluatorTensorProduct; + if (!do_integrate) + { + // Evaluate in 2d + if (n_blocks[0] == n_rows_n) + { + EvalN eval(shape_data[0].shape_values_eo, + shape_data[0].shape_gradients_eo, + {}); + eval.template values<0, true, false>(values_dofs, values); + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template gradients<0, true, false, dim>( + values_dofs, gradients); + eval.template values<0, true, false, dim>( + values_dofs + n_rows_n, gradients + 1); + } + } + else + { + Eval eval(shape_data[1].shape_values_eo, + shape_data[1].shape_gradients_eo, + {}); + eval.template values<0, true, false>(values_dofs, values); + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template gradients<0, true, false, dim>( + values_dofs, gradients); + eval.template values<0, true, false, dim>( + values_dofs + n_rows_t, gradients + 1); + } + } + } + else + { + // Integrate in 2d + if (n_blocks[0] == n_rows_n) + { + EvalN eval(shape_data[0].shape_values_eo, + shape_data[0].shape_gradients_eo, + {}); + if (evaluation_flag & EvaluationFlags::values) + eval.template values<0, false, false>(values, + values_dofs); + if (evaluation_flag & EvaluationFlags::gradients) + { + if (evaluation_flag & EvaluationFlags::values) + eval.template gradients<0, false, true, dim>( + gradients, values_dofs); + else + eval.template gradients<0, false, false, dim>( + gradients, values_dofs); + eval.template values<0, false, false, dim>( + gradients + 1, values_dofs + n_rows_n); + } + } + else + { + Eval eval(shape_data[1].shape_values_eo, + shape_data[1].shape_gradients_eo, + {}); + if (evaluation_flag & EvaluationFlags::values) + eval.template values<0, false, false>(values, + values_dofs); + if (evaluation_flag & EvaluationFlags::gradients) + { + if (evaluation_flag & EvaluationFlags::values) + eval.template gradients<0, false, true, dim>( + gradients, values_dofs); + else + eval.template gradients<0, false, false, dim>( + gradients, values_dofs); + eval.template values<0, false, false, dim>( + gradients + 1, values_dofs + n_rows_t); + } + } + } + } + values += Utilities::pow(n_q_points_1d, dim - 1); + gradients += dim * Utilities::pow(n_q_points_1d, dim - 1); + } + } + }; + + + + template + struct FEFaceNormalEvaluationImpl + { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + template + static void + interpolate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flags, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number *input, + Number *output, + const unsigned int face_no) + { + Assert(static_cast(fe_degree) == + shape_info.data.front().fe_degree || + fe_degree == -1, + ExcInternalError()); + if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) + interpolate_raviart_thomas( + n_components, input, output, flags, face_no, shape_info); + else + interpolate_generic( + n_components, + input, + output, + flags, + face_no, + shape_info.data.front().fe_degree + 1, + shape_info.data.front().shape_data_on_face, + shape_info.dofs_per_component_on_cell, + 3 * shape_info.dofs_per_component_on_face); + } + + /** + * Interpolate the values on the cell quadrature points onto a face. + */ + template + static void + interpolate_quadrature( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flags, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number *input, + Number *output, + const unsigned int face_no) + { + Assert(static_cast(fe_degree + 1) == + shape_info.data.front().n_q_points_1d || + fe_degree == -1, + ExcInternalError()); + + interpolate_generic( + n_components, + input, + output, + flags, + face_no, + shape_info.data.front().quadrature.size(), + shape_info.data.front().quadrature_data_on_face, + shape_info.n_q_points, + shape_info.n_q_points_face); + } + + private: + template + static void + interpolate_generic(const unsigned int n_components, + const Number *input, + Number *output, + const EvaluationFlags::EvaluationFlags flag, + const unsigned int face_no, + const unsigned int n_points_1d, + const std::array, 2> &shape_data, + const unsigned int dofs_per_component_on_cell, + const unsigned int dofs_per_component_on_face) + { + if (face_direction == face_no / 2) + { + constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction); + + const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d; + const int stride = Utilities::pow(n_rows, face_direction); + const std::array n_blocks{ + {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}; + 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}}; + + for (unsigned int c = 0; c < n_components; ++c) + { + if (flag & EvaluationFlags::hessians) + interpolate_to_face(shape_data[face_no % 2].begin(), + n_blocks, + steps, + input, + output, + n_rows, + stride); + else if (flag & EvaluationFlags::gradients) + interpolate_to_face(shape_data[face_no % 2].begin(), + n_blocks, + steps, + input, + output, + n_rows, + stride); + else + interpolate_to_face(shape_data[face_no % 2].begin(), + n_blocks, + steps, + input, + output, + n_rows, + stride); + if (do_evaluate) + { + input += dofs_per_component_on_cell; + output += dofs_per_component_on_face; + } + else + { + output += dofs_per_component_on_cell; + input += dofs_per_component_on_face; + } + } + } + else if (face_direction < dim) + { + interpolate_generic( + n_components, + input, + output, + flag, + face_no, + n_points_1d, + shape_data, + dofs_per_component_on_cell, + dofs_per_component_on_face); + } + } + + template + static void + interpolate_raviart_thomas( + const unsigned int n_components, + const Number *input, + Number *output, + const EvaluationFlags::EvaluationFlags flag, + const unsigned int face_no, + const MatrixFreeFunctions::ShapeInfo &shape_info) + { + if (dim == 1) + { + // This should never happen since the FE_RaviartThomasNodal is not + // defined for dim = 1. It prevents compiler warnings of infinite + // recursion. + Assert(false, ExcInternalError()); + return; + } + + bool increase_max_der = false; + if ((flag & EvaluationFlags::hessians && max_derivative < 2) || + (flag & EvaluationFlags::gradients && max_derivative < 1)) + increase_max_der = true; + + if (face_direction == face_no / 2 && !increase_max_der) + { + 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) + { + // Only increase max_derivative + interpolate_raviart_thomas( + n_components, input, output, flag, face_no, shape_info); + } + else if (face_direction < dim) + { + if (increase_max_der) + { + interpolate_raviart_thomas( + n_components, input, output, flag, face_no, shape_info); + } + else + { + interpolate_raviart_thomas( + n_components, input, output, flag, face_no, shape_info); + } + } + } + }; + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + dst[v] = src_ptr[v]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated load function + template + void + do_vectorized_read(const Number *src_ptr, VectorizedArray &dst) + { + dst.load(src_ptr); + } + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_gather(const Number2 *src_ptr, + const unsigned int *indices, + VectorizedArrayType &dst) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + dst[v] = src_ptr[indices[v]]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated gather function + template + void + do_vectorized_gather(const Number *src_ptr, + const unsigned int *indices, + VectorizedArray &dst) + { + dst.gather(src_ptr, indices); + } + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + dst_ptr[v] += src[v]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated load function + template + void + do_vectorized_add(const VectorizedArray src, Number *dst_ptr) + { + VectorizedArray tmp; + tmp.load(dst_ptr); + (tmp + src).store(dst_ptr); + } + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_scatter_add(const VectorizedArrayType src, + const unsigned int *indices, + Number2 *dst_ptr) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + dst_ptr[indices[v]] += src[v]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated gather function + template + void + do_vectorized_scatter_add(const VectorizedArray src, + const unsigned int *indices, + Number *dst_ptr) + { +#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512 + for (unsigned int v = 0; v < width; ++v) + dst_ptr[indices[v]] += src[v]; +#else + VectorizedArray tmp; + tmp.gather(dst_ptr, indices); + (tmp + src).scatter(indices, dst_ptr); +#endif + } + + + + template + void + adjust_for_face_orientation(const unsigned int dim, + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flag, + const unsigned int *orientation, + const bool integrate, + const std::size_t n_q_points, + Number *tmp_values, + Number *values_quad, + Number *gradients_quad, + Number *hessians_quad) + { + for (unsigned int c = 0; c < n_components; ++c) + { + if (flag & EvaluationFlags::values) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = values_quad[c * n_q_points + orientation[q]]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = values_quad[c * n_q_points + q]; + for (unsigned int q = 0; q < n_q_points; ++q) + values_quad[c * n_q_points + q] = tmp_values[q]; + } + if (flag & EvaluationFlags::gradients) + for (unsigned int d = 0; d < dim; ++d) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = + gradients_quad[(c * n_q_points + orientation[q]) * dim + d]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = + gradients_quad[(c * n_q_points + q) * dim + d]; + for (unsigned int q = 0; q < n_q_points; ++q) + gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q]; + } + if (flag & EvaluationFlags::hessians) + { + const unsigned int hdim = (dim * (dim + 1)) / 2; + for (unsigned int d = 0; d < hdim; ++d) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points + + orientation[q]]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = + hessians_quad[(c * hdim + d) * n_q_points + q]; + for (unsigned int q = 0; q < n_q_points; ++q) + hessians_quad[(c * hdim + d) * n_q_points + q] = + tmp_values[q]; + } + } + } + } + + + + template + void + adjust_for_face_orientation_per_lane( + const unsigned int dim, + const unsigned int n_components, + const unsigned int v, + const EvaluationFlags::EvaluationFlags flag, + const unsigned int *orientation, + const bool integrate, + const std::size_t n_q_points, + Number *tmp_values, + VectorizedArrayType *values_quad, + VectorizedArrayType *gradients_quad = nullptr, + VectorizedArrayType *hessians_quad = nullptr) + { + for (unsigned int c = 0; c < n_components; ++c) + { + if (flag & EvaluationFlags::values) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v]; + for (unsigned int q = 0; q < n_q_points; ++q) + values_quad[c * n_q_points + q][v] = tmp_values[q]; + } + if (flag & EvaluationFlags::gradients) + for (unsigned int d = 0; d < dim; ++d) + { + Assert(gradients_quad != nullptr, ExcInternalError()); + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = + gradients_quad[(c * n_q_points + orientation[q]) * dim + d] + [v]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = + gradients_quad[(c * n_q_points + q) * dim + d][v]; + for (unsigned int q = 0; q < n_q_points; ++q) + gradients_quad[(c * n_q_points + q) * dim + d][v] = + tmp_values[q]; + } + if (flag & EvaluationFlags::hessians) + { + Assert(hessians_quad != nullptr, ExcInternalError()); + const unsigned int hdim = (dim * (dim + 1)) / 2; + for (unsigned int d = 0; d < hdim; ++d) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points + + orientation[q]][v]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = + hessians_quad[(c * hdim + d) * n_q_points + q][v]; + for (unsigned int q = 0; q < n_q_points; ++q) + hessians_quad[(c * hdim + d) * n_q_points + q][v] = + tmp_values[q]; + } + } + } + } + + + + template + struct FEFaceEvaluationImplEvaluateSelector + { + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + if (shape_info.element_type == MatrixFreeFunctions::tensor_none) + { + Assert((fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) == false, + ExcNotImplemented()); + + const unsigned int face_no = fe_eval.get_face_no(); + const unsigned int face_orientation = fe_eval.get_face_orientation(); + const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; + const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; + + using Eval = + EvaluatorTensorProduct; + + if (evaluation_flag & EvaluationFlags::values) + { + const auto *const shape_values = + &shape_data.shape_values_face(face_no, face_orientation, 0); + + auto *values_quad_ptr = fe_eval.begin_values(); + auto *values_dofs_actual_ptr = values_dofs; + + Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); + for (unsigned int c = 0; c < n_components; ++c) + { + eval.template values<0, true, false>(values_dofs_actual_ptr, + values_quad_ptr); + + values_quad_ptr += n_q_points; + values_dofs_actual_ptr += n_dofs; + } + } + + if (evaluation_flag & EvaluationFlags::gradients) + { + auto *gradients_quad_ptr = fe_eval.begin_gradients(); + const auto *values_dofs_actual_ptr = values_dofs; + + std::array shape_gradients; + for (unsigned int d = 0; d < dim; ++d) + shape_gradients[d] = &shape_data.shape_gradients_face( + face_no, face_orientation, d, 0); + + for (unsigned int c = 0; c < n_components; ++c) + { + for (unsigned int d = 0; d < dim; ++d) + { + Eval eval(nullptr, + shape_gradients[d], + nullptr, + n_dofs, + n_q_points); + + eval.template gradients<0, true, false, dim>( + values_dofs_actual_ptr, gradients_quad_ptr + d); + } + gradients_quad_ptr += n_q_points * dim; + values_dofs_actual_ptr += n_dofs; + } + } + + Assert(!(evaluation_flag & EvaluationFlags::hessians), + ExcNotImplemented()); + + return true; + } + + const unsigned int dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(shape_data.fe_degree + 1, dim - 1); + + // Note: we always keep storage of values, 1st and 2nd derivatives in an + // array, so reserve space for all three here + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_face; + + bool use_vectorization = true; + + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) // exterior faces in the ECL loop + for (unsigned int v = 0; v < Number::size(); ++v) + if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int && + fe_eval.get_face_no(v) != fe_eval.get_face_no(0)) + use_vectorization = false; + + if (use_vectorization == false) + { + for (unsigned int v = 0; v < Number::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + { + for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; + ++i) + temp[i][v] = 0; + continue; + } + + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + evaluation_flag, + shape_info, + values_dofs, + scratch_data, + fe_eval.get_face_no(v)); + + for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; + ++i) + temp[i][v] = scratch_data[i][v]; + } + } + else + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + evaluation_flag, + shape_info, + values_dofs, + temp, + fe_eval.get_face_no()); + + const unsigned int subface_index = fe_eval.get_subface_index(); + constexpr unsigned int n_q_points_1d_actual = + fe_degree > -1 ? n_q_points_1d : 0; + + if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) + { + FEFaceEvaluationImplRaviartThomas:: + template evaluate_or_integrate_in_face( + evaluation_flag, + fe_eval.get_shape_info().data, + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + scratch_data, + subface_index, + fe_eval.get_face_no() / 2); + } + else if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl::evaluate_in_face(n_components, + evaluation_flag, + shape_data, + temp, + fe_eval.begin_values(), + fe_eval + .begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + else + FEFaceEvaluationImpl::evaluate_in_face(n_components, + evaluation_flag, + shape_data, + temp, + fe_eval.begin_values(), + fe_eval + .begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + + if (use_vectorization == false) + { + for (unsigned int v = 0; v < Number::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + if (fe_eval.get_face_orientation(v) != 0) + adjust_for_face_orientation_per_lane( + dim, + n_components, + v, + evaluation_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(v), 0), + false, + shape_info.n_q_points_face, + &temp[0][0], + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + } + } + else if (fe_eval.get_face_orientation() != 0) + adjust_for_face_orientation( + dim, + n_components, + evaluation_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(), 0), + false, + shape_info.n_q_points_face, + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + + return false; + } + }; + + + + template + struct FEFaceEvaluationImplIntegrateSelector + { + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + if (shape_info.element_type == MatrixFreeFunctions::tensor_none) + { + Assert((fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) == false, + ExcNotImplemented()); + + const unsigned int face_no = fe_eval.get_face_no(); + const unsigned int face_orientation = fe_eval.get_face_orientation(); + const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; + const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; + + using Eval = + EvaluatorTensorProduct; + + if (integration_flag & EvaluationFlags::values) + { + const auto *const shape_values = + &shape_data.shape_values_face(face_no, face_orientation, 0); + + auto *values_quad_ptr = fe_eval.begin_values(); + auto *values_dofs_actual_ptr = values_dofs; + + Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); + for (unsigned int c = 0; c < n_components; ++c) + { + eval.template values<0, false, false>(values_quad_ptr, + values_dofs_actual_ptr); + + values_quad_ptr += n_q_points; + values_dofs_actual_ptr += n_dofs; + } + } + + if (integration_flag & EvaluationFlags::gradients) + { + auto *gradients_quad_ptr = fe_eval.begin_gradients(); + auto *values_dofs_actual_ptr = values_dofs; + + std::array shape_gradients; + for (unsigned int d = 0; d < dim; ++d) + shape_gradients[d] = &shape_data.shape_gradients_face( + face_no, face_orientation, d, 0); + + for (unsigned int c = 0; c < n_components; ++c) + { + for (unsigned int d = 0; d < dim; ++d) + { + Eval eval(nullptr, + shape_gradients[d], + nullptr, + n_dofs, + n_q_points); + + if (!(integration_flag & EvaluationFlags::values) && + d == 0) + eval.template gradients<0, false, false, dim>( + gradients_quad_ptr + d, values_dofs_actual_ptr); + else + eval.template gradients<0, false, true, dim>( + gradients_quad_ptr + d, values_dofs_actual_ptr); + } + gradients_quad_ptr += n_q_points * dim; + values_dofs_actual_ptr += n_dofs; + } + } + + Assert(!(integration_flag & EvaluationFlags::hessians), + ExcNotImplemented()); + + return true; + } + + const unsigned int dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(shape_data.fe_degree + 1, dim - 1); + + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_face; + + bool use_vectorization = true; + + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) // exterior faces in the ECL loop + use_vectorization = + fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int && + std::all_of(fe_eval.get_cell_ids().begin() + 1, + fe_eval.get_cell_ids().end(), + [&](const auto &v) { + return v == fe_eval.get_cell_ids()[0] || + v == numbers::invalid_unsigned_int; + }); + + if (use_vectorization == false) + { + for (unsigned int v = 0; v < Number::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + if (fe_eval.get_face_orientation(v) != 0) + adjust_for_face_orientation_per_lane( + dim, + n_components, + v, + integration_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(v), 0), + true, + shape_info.n_q_points_face, + &temp[0][0], + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + } + } + else if (fe_eval.get_face_orientation() != 0) + adjust_for_face_orientation( + dim, + n_components, + integration_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(), 0), + true, + shape_info.n_q_points_face, + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + + const unsigned int n_q_points_1d_actual = + fe_degree > -1 ? n_q_points_1d : 0; + const unsigned int subface_index = fe_eval.get_subface_index(); + + if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) + { + FEFaceEvaluationImplRaviartThomas:: + template evaluate_or_integrate_in_face( + integration_flag, + fe_eval.get_shape_info().data, + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + scratch_data, + subface_index, + fe_eval.get_face_no() / 2); + } + else if (fe_degree > -1 && + fe_eval.get_subface_index() >= + GeometryInfo::max_children_per_cell && + shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d_actual, + Number>::integrate_in_face(n_components, + integration_flag, + shape_data, + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + else + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d_actual, + Number>::integrate_in_face(n_components, + integration_flag, + shape_data, + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + + if (use_vectorization == false) + { + for (unsigned int v = 0; v < Number::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flag, + shape_info, + values_dofs, + scratch_data, + fe_eval.get_face_no(v)); + + for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; + ++i) + temp[i][v] = scratch_data[i][v]; + } + } + else + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flag, + shape_info, + temp, + values_dofs, + fe_eval.get_face_no()); + return false; + } + }; + + + + template + void + fe_face_evaluation_process_and_io( + Processor &proc, + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + typename Processor::Number2_ *global_vector_ptr, + const std::vector> *sm_ptr, + const EvaluationData &fe_eval, + typename Processor::VectorizedArrayType_ *temp1) + { + constexpr int dim = Processor::dim_; + constexpr int fe_degree = Processor::fe_degree_; + using VectorizedArrayType = typename Processor::VectorizedArrayType_; + constexpr int n_lanes = VectorizedArrayType::size(); + + using Number = typename Processor::Number_; + using Number2_ = typename Processor::Number2_; + + const auto &shape_data = fe_eval.get_shape_info().data.front(); + constexpr bool integrate = Processor::do_integrate; + const unsigned int face_no = fe_eval.get_face_no(); + const auto &dof_info = fe_eval.get_dof_info(); + const unsigned int cell = fe_eval.get_cell_or_face_batch_id(); + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index = + fe_eval.get_dof_access_index(); + AssertIndexRange(cell, + dof_info.index_storage_variants[dof_access_index].size()); + constexpr unsigned int dofs_per_face = + Utilities::pow(fe_degree + 1, dim - 1); + const unsigned int subface_index = fe_eval.get_subface_index(); + + const unsigned int n_filled_lanes = + dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; + + bool all_faces_are_same = n_filled_lanes == n_lanes; + if (n_face_orientations == n_lanes) + for (unsigned int v = 1; v < n_lanes; ++v) + if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) || + fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0)) + { + all_faces_are_same = false; + break; + } + + // check for re-orientation ... + std::array orientation = {}; + + if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same && + fe_eval.is_interior_face() == 0) + for (unsigned int v = 0; v < n_lanes; ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + if (shape_data.nodal_at_cell_boundaries && + fe_eval.get_face_orientation(v) != 0) + { + // ... and in case we detect a re-orientation, go to the other + // version of this function that actually allows for this + if (subface_index == GeometryInfo::max_children_per_cell && + check_face_orientations == false) + { + fe_face_evaluation_process_and_io(proc, + n_components, + evaluation_flag, + global_vector_ptr, + sm_ptr, + fe_eval, + temp1); + return; + } + orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs( + fe_eval.get_face_orientation(v), 0); + } + } + else if (dim == 3 && fe_eval.get_face_orientation() != 0) + { + // go to the other version of this function + if (subface_index == GeometryInfo::max_children_per_cell && + check_face_orientations == false) + { + fe_face_evaluation_process_and_io(proc, + n_components, + evaluation_flag, + global_vector_ptr, + sm_ptr, + fe_eval, + temp1); + return; + } + for (unsigned int v = 0; v < n_face_orientations; ++v) + orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs( + fe_eval.get_face_orientation(), 0); + } + + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + VectorizedArrayType grad_weight = + shape_data + .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) : + (1 + face_no % 2))]; + + // face_to_cell_index_hermite + std::array index_array_hermite = + {}; + if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) + { + if (n_face_orientations == 1) + index_array_hermite[0] = + &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0); + else + { + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + const auto face_no = fe_eval.get_face_no(v); + + grad_weight[v] = + shape_data.shape_data_on_face[0][fe_degree + + (integrate ? + (2 - (face_no % 2)) : + (1 + (face_no % 2)))][0]; + + index_array_hermite[v] = + &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, + 0); + } + } + } + + // face_to_cell_index_nodal + std::array index_array_nodal = + {}; + if (shape_data.nodal_at_cell_boundaries == true) + { + if (n_face_orientations == 1) + index_array_nodal[0] = + &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0); + else + { + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + const auto face_no = fe_eval.get_face_no(v); + + index_array_nodal[v] = + &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, + 0); + } + } + } + + + const auto reorientate = [&](const unsigned int v, const unsigned int i) { + return (!check_face_orientations || orientation[v] == nullptr) ? + i : + orientation[v][i]; + }; + + const unsigned int cell_index = + dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ? + fe_eval.get_cell_ids()[0] : + cell * n_lanes; + const unsigned int *dof_indices = + &dof_info.dof_indices_contiguous[dof_access_index][cell_index]; + + for (unsigned int comp = 0; comp < n_components; ++comp) + { + const std::size_t index_offset = + dof_info.component_dof_indices_offset + [fe_eval.get_active_fe_index()] + [fe_eval.get_first_selected_component()] + + comp * Utilities::pow(fe_degree + 1, dim); + + // case 1: contiguous and interleaved indices + if (n_face_orientations == 1 && + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous) + { + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + n_lanes); + Number2_ *vector_ptr = + global_vector_ptr + dof_indices[0] + index_offset * n_lanes; + + if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + Assert(n_face_orientations == 1, ExcNotImplemented()); + + const unsigned int ind1 = index_array_hermite[0][2 * i]; + const unsigned int ind2 = index_array_hermite[0][2 * i + 1]; + const unsigned int i_ = reorientate(0, i); + proc.hermite_grad_vectorized(temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1 * n_lanes, + vector_ptr + ind2 * n_lanes, + grad_weight); + } + } + else + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + Assert(n_face_orientations == 1, ExcNotImplemented()); + + const unsigned int i_ = reorientate(0, i); + const unsigned int ind = index_array_nodal[0][i]; + proc.value_vectorized(temp1[i_], + vector_ptr + ind * n_lanes); + } + } + } + + // case 2: contiguous and interleaved indices with fixed stride + else if (n_face_orientations == 1 && + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided) + { + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + n_lanes); + Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes; + if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + Assert(n_face_orientations == 1, ExcNotImplemented()); + + const unsigned int i_ = reorientate(0, i); + const unsigned int ind1 = + index_array_hermite[0][2 * i] * n_lanes; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1] * n_lanes; + proc.hermite_grad_vectorized_indexed( + temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1, + vector_ptr + ind2, + grad_weight, + dof_indices, + dof_indices); + } + } + else + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + Assert(n_face_orientations == 1, ExcNotImplemented()); + + const unsigned int i_ = reorientate(0, i); + const unsigned int ind = index_array_nodal[0][i] * n_lanes; + proc.value_vectorized_indexed(temp1[i_], + vector_ptr + ind, + dof_indices); + } + } + } + + // case 3: contiguous and interleaved indices with mixed stride + else if (n_face_orientations == 1 && + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_mixed_strides) + { + const unsigned int *strides = + &dof_info.dof_indices_interleave_strides[dof_access_index] + [cell * n_lanes]; + unsigned int indices[n_lanes]; + for (unsigned int v = 0; v < n_lanes; ++v) + indices[v] = dof_indices[v] + index_offset * strides[v]; + const unsigned int n_filled_lanes = + dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; + + if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) + { + if (n_filled_lanes == n_lanes) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + Assert(n_face_orientations == 1, ExcNotImplemented()); + + const unsigned int i_ = reorientate(0, i); + unsigned int ind1[n_lanes]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; v < n_lanes; ++v) + ind1[v] = indices[v] + + index_array_hermite[0][2 * i] * strides[v]; + unsigned int ind2[n_lanes]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; v < n_lanes; ++v) + ind2[v] = + indices[v] + + // TODO + index_array_hermite[0][2 * i + 1] * strides[v]; + proc.hermite_grad_vectorized_indexed( + temp1[i_], + temp1[i_ + dofs_per_face], + global_vector_ptr, + global_vector_ptr, + grad_weight, + ind1, + ind2); + } + else + { + if (integrate == false) + for (unsigned int i = 0; i < 2 * dofs_per_face; ++i) + temp1[i] = VectorizedArrayType(); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int i_ = + reorientate(n_face_orientations == 1 ? 0 : v, i); + proc.hermite_grad( + temp1[i_][v], + temp1[i_ + dofs_per_face][v], + global_vector_ptr + [indices[v] + + index_array_hermite + [n_face_orientations == 1 ? 0 : v][2 * i] * + strides[v]], + global_vector_ptr + [indices[v] + + index_array_hermite[n_face_orientations == 1 ? + 0 : + v][2 * i + 1] * + strides[v]], + grad_weight[n_face_orientations == 1 ? 0 : v]); + } + } + } + else + { + if (n_filled_lanes == n_lanes) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + Assert(n_face_orientations == 1, ExcInternalError()); + unsigned int ind[n_lanes]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; v < n_lanes; ++v) + ind[v] = + indices[v] + index_array_nodal[0][i] * strides[v]; + const unsigned int i_ = reorientate(0, i); + proc.value_vectorized_indexed(temp1[i_], + global_vector_ptr, + ind); + } + else + { + if (integrate == false) + for (unsigned int i = 0; i < dofs_per_face; ++i) + temp1[i] = VectorizedArrayType(); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int i = 0; i < dofs_per_face; ++i) + proc.value( + temp1[reorientate(n_face_orientations == 1 ? 0 : v, + i)][v], + global_vector_ptr + [indices[v] + + index_array_nodal[n_face_orientations == 1 ? 0 : v] + [i] * + strides[v]]); + } + } + } + + // case 4: contiguous indices without interleaving + else if (n_face_orientations > 1 || + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous) + { + Number2_ *vector_ptr = global_vector_ptr + index_offset; + + const bool vectorization_possible = + all_faces_are_same && (sm_ptr == nullptr); + + std::array vector_ptrs; + std::array reordered_indices; + + if (vectorization_possible == false) + { + vector_ptrs = {}; + if (n_face_orientations == 1) + { + for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (sm_ptr == nullptr) + { + vector_ptrs[v] = vector_ptr + dof_indices[v]; + } + else + { + const auto &temp = + dof_info + .dof_indices_contiguous_sm[dof_access_index] + [cell * n_lanes + v]; + vector_ptrs[v] = const_cast( + sm_ptr->operator[](temp.first).data() + + temp.second + index_offset); + } + } + else if (n_face_orientations == n_lanes) + { + const auto &cells = fe_eval.get_cell_ids(); + for (unsigned int v = 0; v < n_lanes; ++v) + if (cells[v] != numbers::invalid_unsigned_int) + { + if (sm_ptr == nullptr) + { + vector_ptrs[v] = + vector_ptr + + dof_info + .dof_indices_contiguous[dof_access_index] + [cells[v]]; + } + else + { + const auto &temp = + dof_info + .dof_indices_contiguous_sm[dof_access_index] + [cells[v]]; + vector_ptrs[v] = const_cast( + sm_ptr->operator[](temp.first).data() + + temp.second + index_offset); + } + } + } + else + { + Assert(false, ExcNotImplemented()); + } + } + else if (n_face_orientations == n_lanes) + { + for (unsigned int v = 0; v < n_lanes; ++v) + reordered_indices[v] = + dof_info.dof_indices_contiguous[dof_access_index] + [fe_eval.get_cell_ids()[v]]; + dof_indices = reordered_indices.data(); + } + + if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) + { + if (vectorization_possible) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = index_array_hermite[0][2 * i]; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1]; + const unsigned int i_ = reorientate(0, i); + + proc.hermite_grad_vectorized_indexed( + temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1, + vector_ptr + ind2, + grad_weight, + dof_indices, + dof_indices); + } + else if (n_face_orientations == 1) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = index_array_hermite[0][2 * i]; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1]; + const unsigned int i_ = reorientate(0, i); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + proc.hermite_grad(temp1[i_][v], + temp1[i_ + dofs_per_face][v], + vector_ptrs[v][ind1], + vector_ptrs[v][ind2], + grad_weight[v]); + + if (integrate == false) + for (unsigned int v = n_filled_lanes; v < n_lanes; ++v) + { + temp1[i][v] = 0.0; + temp1[i + dofs_per_face][v] = 0.0; + } + } + else + { + if (integrate == false && n_filled_lanes < n_lanes) + for (unsigned int i = 0; i < dofs_per_face; ++i) + temp1[i] = temp1[i + dofs_per_face] = Number(); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int i = 0; i < dofs_per_face; ++i) + proc.hermite_grad( + temp1[reorientate(v, i)][v], + temp1[reorientate(v, i) + dofs_per_face][v], + vector_ptrs[v][index_array_hermite[v][2 * i]], + vector_ptrs[v][index_array_hermite[v][2 * i + 1]], + grad_weight[v]); + } + } + else + { + if (vectorization_possible) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = index_array_nodal[0][i]; + const unsigned int i_ = reorientate(0, i); + + proc.value_vectorized_indexed(temp1[i_], + vector_ptr + ind, + dof_indices); + } + else if (n_face_orientations == 1) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = index_array_nodal[0][i]; + const unsigned int i_ = reorientate(0, i); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + proc.value(temp1[i_][v], vector_ptrs[v][ind]); + + if (integrate == false) + for (unsigned int v = n_filled_lanes; v < n_lanes; ++v) + temp1[i_][v] = 0.0; + } + else + { + if (integrate == false && n_filled_lanes < n_lanes) + for (unsigned int i = 0; i < dofs_per_face; ++i) + temp1[i] = Number(); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int i = 0; i < dofs_per_face; ++i) + proc.value(temp1[reorientate(v, i)][v], + vector_ptrs[v][index_array_nodal[v][i]]); + } + } + } + else + { + // We should not end up here, this should be caught by + // FEFaceEvaluationImplGatherEvaluateSelector::supports() + Assert(false, ExcInternalError()); + } + temp1 += 3 * dofs_per_face; + } + } + + + + template + struct FEFaceEvaluationImplGatherEvaluateSelector + { + using Number = typename VectorizedArrayType::value_type; + + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number2 *src_ptr, + const std::vector> *sm_ptr, + FEEvaluationData &fe_eval) + { + Assert(fe_degree > -1, ExcInternalError()); + Assert(fe_eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric_no_collocation, + ExcInternalError()); + + const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1); + + VectorizedArrayType *temp = fe_eval.get_scratch_data().begin(); + VectorizedArrayType *scratch_data = + temp + 3 * n_components * dofs_per_face; + + Processor p; + + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) + fe_face_evaluation_process_and_io( + p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp); + else + fe_face_evaluation_process_and_io<1>( + p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp); + + const unsigned int subface_index = fe_eval.get_subface_index(); + + if (subface_index >= GeometryInfo::max_children_per_cell) + FEFaceEvaluationImpl:: + evaluate_in_face(n_components, + evaluation_flag, + fe_eval.get_shape_info().data.front(), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + else + FEFaceEvaluationImpl:: + evaluate_in_face(n_components, + evaluation_flag, + fe_eval.get_shape_info().data.front(), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + + // re-orientation for cases not possible with above algorithm + if (subface_index < GeometryInfo::max_children_per_cell) + { + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == + numbers::invalid_unsigned_int) + continue; + + if (fe_eval.get_face_orientation(v) != 0) + adjust_for_face_orientation_per_lane( + dim, + n_components, + v, + evaluation_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(v), 0), + false, + Utilities::pow(n_q_points_1d, dim - 1), + &temp[0][0], + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + } + } + else if (fe_eval.get_face_orientation() != 0) + adjust_for_face_orientation( + dim, + n_components, + evaluation_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(), 0), + false, + Utilities::pow(n_q_points_1d, dim - 1), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + } + + return false; + } + + template + static bool + supports(const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number2 *vector_ptr, + MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage) + { + const unsigned int fe_degree = shape_info.data.front().fe_degree; + if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries || + (evaluation_flag & EvaluationFlags::gradients && + (fe_degree < 2 || + shape_info.data.front().element_type != + MatrixFreeFunctions::tensor_symmetric_hermite)) || + (evaluation_flag & EvaluationFlags::hessians) || + vector_ptr == nullptr || + shape_info.data.front().element_type > + MatrixFreeFunctions::tensor_symmetric_no_collocation || + storage < + MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) + return false; + else + return true; + } + + private: + template + struct Processor + { + static const bool do_integrate = false; + static const int dim_ = dim; + static const int fe_degree_ = fe_degree; + using VectorizedArrayType_ = VectorizedArrayType; + using Number_ = Number; + using Number2_ = const Number2; + + template + void + hermite_grad_vectorized(T0 &temp_1, + T0 &temp_2, + const T1 src_ptr_1, + const T1 src_ptr_2, + const T2 &grad_weight) + { + do_vectorized_read(src_ptr_1, temp_1); + do_vectorized_read(src_ptr_2, temp_2); + temp_2 = grad_weight * (temp_1 - temp_2); + } + + template + void + value_vectorized(T1 &temp, const T2 src_ptr) + { + do_vectorized_read(src_ptr, temp); + } + + template + void + hermite_grad_vectorized_indexed(T0 &temp_1, + T0 &temp_2, + const T1 src_ptr_1, + const T1 src_ptr_2, + const T2 &grad_weight, + const T3 &indices_1, + const T3 &indices_2) + { + do_vectorized_gather(src_ptr_1, indices_1, temp_1); + do_vectorized_gather(src_ptr_2, indices_2, temp_2); + temp_2 = grad_weight * (temp_1 - temp_2); + } + + template + void + value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices) + { + do_vectorized_gather(src_ptr, indices, temp); + } + + template + void + hermite_grad(T0 &temp_1, + T0 &temp_2, + const T1 &src_ptr_1, + const T1 &src_ptr_2, + const T2 &grad_weight) + { + // case 3a) + temp_1 = src_ptr_1; + temp_2 = grad_weight * (temp_1 - src_ptr_2); + } + + template + void + value(T1 &temp, const T2 &src_ptr) + { + // case 3b) + temp = src_ptr; + } + }; + }; + + + + template + struct FEFaceEvaluationImplIntegrateScatterSelector + { + using Number = typename VectorizedArrayType::value_type; + + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number2 *dst_ptr, + const std::vector> *sm_ptr, + FEEvaluationData &fe_eval) + { + Assert(fe_degree > -1, ExcInternalError()); + Assert(fe_eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric_no_collocation, + ExcInternalError()); + + const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1); + + VectorizedArrayType *temp = fe_eval.get_scratch_data().begin(); + VectorizedArrayType *scratch_data = + temp + 3 * n_components * dofs_per_face; + + const unsigned int subface_index = fe_eval.get_subface_index(); + + // re-orientation for cases not possible with the io function below + if (subface_index < GeometryInfo::max_children_per_cell) + { + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) + continue; + + if (fe_eval.get_face_orientation(v) != 0) + adjust_for_face_orientation_per_lane( + dim, + n_components, + v, + integration_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(v), 0), + true, + Utilities::pow(n_q_points_1d, dim - 1), + &temp[0][0], + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + } + else if (fe_eval.get_face_orientation() != 0) + adjust_for_face_orientation( + dim, + n_components, + integration_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(), 0), + true, + Utilities::pow(n_q_points_1d, dim - 1), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); + } + + if (fe_degree > -1 && fe_eval.get_subface_index() >= + GeometryInfo::max_children_per_cell) + FEFaceEvaluationImpl:: + integrate_in_face(n_components, + integration_flag, + fe_eval.get_shape_info().data.front(), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + else + FEFaceEvaluationImpl:: + integrate_in_face(n_components, + integration_flag, + fe_eval.get_shape_info().data.front(), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians(), + scratch_data, + subface_index); + + Processor p; + + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) + fe_face_evaluation_process_and_io( + p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp); + else + fe_face_evaluation_process_and_io<1>( + p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp); + + return false; + } + + private: + template + struct Processor + { + static const bool do_integrate = true; + static const int dim_ = dim; + static const int fe_degree_ = fe_degree; + using VectorizedArrayType_ = VectorizedArrayType; + using Number_ = Number; + using Number2_ = Number2; + + template + void + hermite_grad_vectorized(const T0 &temp_1, + const T1 &temp_2, + T2 dst_ptr_1, + T3 dst_ptr_2, + const T4 &grad_weight) + { + // case 1a) + const VectorizedArrayType val = temp_1 - grad_weight * temp_2; + const VectorizedArrayType grad = grad_weight * temp_2; + do_vectorized_add(val, dst_ptr_1); + do_vectorized_add(grad, dst_ptr_2); + } + + template + void + value_vectorized(const T0 &temp, T1 dst_ptr) + { + // case 1b) + do_vectorized_add(temp, dst_ptr); + } + + template + void + hermite_grad_vectorized_indexed(const T0 &temp_1, + const T0 &temp_2, + T1 dst_ptr_1, + T1 dst_ptr_2, + const T2 &grad_weight, + const T3 &indices_1, + const T3 &indices_2) + { + // case 2a) + const VectorizedArrayType val = temp_1 - grad_weight * temp_2; + const VectorizedArrayType grad = grad_weight * temp_2; + do_vectorized_scatter_add(val, indices_1, dst_ptr_1); + do_vectorized_scatter_add(grad, indices_2, dst_ptr_2); + } + + template + void + value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices) + { + // case 2b) + do_vectorized_scatter_add(temp, indices, dst_ptr); + } + + template + void + hermite_grad(const T0 &temp_1, + const T0 &temp_2, + T1 &dst_ptr_1, + T1 &dst_ptr_2, + const T2 &grad_weight) + { + // case 3a) + const Number val = temp_1 - grad_weight * temp_2; + const Number grad = grad_weight * temp_2; + dst_ptr_1 += val; + dst_ptr_2 += grad; + } + + template + void + value(const T0 &temp, T1 &dst_ptr) + { + // case 3b) + dst_ptr += temp; + } + }; + }; +} // end of namespace internal + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h index a429df6208..74c7c4899a 100644 --- a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h @@ -20,8 +20,7 @@ #include -#include -#include +#include #include #include diff --git a/include/deal.II/matrix_free/evaluation_template_factory_internal.h b/include/deal.II/matrix_free/evaluation_template_factory_internal.h index c2378a9812..ea52c477dd 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory_internal.h +++ b/include/deal.II/matrix_free/evaluation_template_factory_internal.h @@ -28,6 +28,21 @@ DEAL_II_NAMESPACE_OPEN namespace internal { + /** + * This struct is used to implement + * FEEvaluation::fast_evaluation_supported() and + * FEFaceEvaluation::fast_evaluation_supported(). + */ + struct FastEvaluationSupported + { + template + static bool + run() + { + return fe_degree != -1; + } + }; + template bool instantiation_helper_run(const unsigned int given_degree, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index cf7d1c4057..9c3bcfdb33 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -31,6 +31,7 @@ #include #include +#include #include #include #include diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 696a16c8b2..9c7a568542 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 0f8f660cc6..837510b986 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -24,6 +24,8 @@ #include #include +#include + DEAL_II_NAMESPACE_OPEN @@ -1653,6 +1655,173 @@ namespace internal + template + struct EvaluatorTensorProductAnisotropic + { + template + static void + normal(const MatrixFreeFunctions::UnivariateShapeData &data, + const Number *in, + Number *out, + const bool add_into_result = false, + const int subface_index_1d = 0) + { + AssertIndexRange(direction, dim); + AssertDimension(fe_degree, data.fe_degree); + AssertDimension(n_q_points_1d, data.n_q_points_1d); + constexpr int n_rows = fe_degree + 1; + constexpr int n_columns = n_q_points_1d; + constexpr int mm = contract_over_rows ? n_rows : n_columns; + constexpr int nn = contract_over_rows ? n_columns : n_rows; + const Number2 *shape_data = + symmetric_evaluate ? + data.shape_values_eo.data() : + data.values_within_subface[subface_index_1d].data(); + Assert(shape_data != nullptr, ExcNotInitialized()); + Assert(contract_over_rows == false || !add_into_result, + ExcMessage("Cannot add into result if contract_over_rows = true")); + + constexpr int n_blocks1 = Utilities::pow(fe_degree, direction); + constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1); + constexpr int stride_in = contract_over_rows ? 1 : stride; + constexpr int stride_out = contract_over_rows ? stride : 1; + constexpr EvaluatorVariant variant = + symmetric_evaluate ? evaluate_evenodd : evaluate_general; + + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + if (contract_over_rows == false && add_into_result) + apply_matrix_vector_product(shape_data, in, out); + else + apply_matrix_vector_product(shape_data, in, out); + + in += stride_in; + out += stride_out; + } + in += n_blocks1 * (mm - 1) * stride_in; + out += n_blocks1 * (nn - 1) * stride_out; + } + } + + template + static void + tangential(const MatrixFreeFunctions::UnivariateShapeData &data, + const Number *in, + Number *out, + const int subface_index_1d = 0) + { + AssertIndexRange(direction, dim); + AssertDimension(fe_degree - 1, data.fe_degree); + AssertDimension(n_q_points_1d, data.n_q_points_1d); + static_assert(direction != normal_direction, + "Cannot interpolate tangentially in normal direction"); + + constexpr int n_rows = fe_degree; + constexpr int n_columns = n_q_points_1d; + const Number2 *shape_data = + symmetric_evaluate ? + data.shape_values_eo.data() : + data.values_within_subface[subface_index_1d].data(); + Assert(shape_data != nullptr, ExcNotInitialized()); + + constexpr int n_blocks1 = + (direction > normal_direction) ? + Utilities::pow(n_q_points_1d, direction) : + (direction > 0 ? + (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) : + 1); + constexpr int n_blocks2 = + (direction > normal_direction) ? + Utilities::pow(fe_degree, dim - 1 - direction) : + ((direction + 1 < dim) ? + (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) : + 1); + + constexpr EvaluatorVariant variant = + symmetric_evaluate ? evaluate_evenodd : evaluate_general; + + // Since we may perform an in-place interpolation, we must run the step + // expanding the size of the basis backward ('contract_over_rows' aka + // 'evaluate' case), so shift the pointers and decrement during the loop + if (contract_over_rows) + { + in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1; + out += + stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1); + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + apply_matrix_vector_product(shape_data, in, out); + + --in; + out -= stride; + } + in -= n_blocks1 * (n_rows - 1); + out -= n_blocks1 * (n_columns - 1) * stride; + } + } + else + { + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + apply_matrix_vector_product(shape_data, in, out); + + in += stride; + ++out; + } + in += n_blocks1 * (n_columns - 1) * stride; + out += n_blocks1 * (n_rows - 1); + } + } + } + }; + + + /** * This function applies the tensor product operation to produce face values * from cell values. The algorithm involved here can be interpreted as the @@ -1767,6 +1936,23 @@ namespace internal + /** + * Helper function to specify whether a transformation to collocation should + * be used: It should give correct results (first condition), we need to be + * able to initialize the fields in shape_info.templates.h from the + * polynomials (second condition), and it should be the most efficient + * choice in terms of operation counts (third condition). + */ + constexpr bool + use_collocation_evaluation(const unsigned int fe_degree, + const unsigned int n_q_points_1d) + { + return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) && + (n_q_points_1d <= 3 * fe_degree / 2 + 1); + } + + + /** * This function performs the opposite operation to the interpolate_to_face * function, done as the last step in sum factorization to embed face values -- 2.39.5