From: Martin Kronbichler Date: Thu, 9 Dec 2021 17:13:08 +0000 (+0100) Subject: Simplify calls to FEEvaluation kernels X-Git-Tag: v9.4.0-rc1~740^2~15 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3bbcc17dedb564b459be940c57a3757b633e3afd;p=dealii.git Simplify calls to FEEvaluation kernels --- diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 306b257d5c..ddb9a63352 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -506,30 +506,6 @@ public: */ mutable AlignedVector scratch; - /** - * In case the quadrature rule given represents a tensor product - * the values at the mapped support points are stored in this object. - */ - mutable AlignedVector values_dofs; - - /** - * In case the quadrature rule given represents a tensor product - * the values at the quadrature points are stored in this object. - */ - mutable AlignedVector values_quad; - - /** - * In case the quadrature rule given represents a tensor product - * the gradients at the quadrature points are stored in this object. - */ - mutable AlignedVector gradients_quad; - - /** - * In case the quadrature rule given represents a tensor product - * the hessians at the quadrature points are stored in this object. - */ - mutable AlignedVector hessians_quad; - /** * Indicates whether the given Quadrature object is a tensor product. */ diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 9df99a0404..3d40e5154d 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -35,6 +35,8 @@ #include #include +#include +#include #include #include @@ -1128,16 +1130,18 @@ namespace internal return; } + internal::MatrixFreeFunctions:: + MappingInfoStorage + temp_data; + temp_data.descriptor.resize(1); + temp_data.descriptor[0].n_q_points = n_q_points; + FEEvaluationBaseData eval( + {&data.shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0); + // prepare arrays if (evaluation_flag != EvaluationFlags::nothing) { - data.values_dofs.resize(n_comp * n_shape_values); - data.values_quad.resize(n_comp * n_q_points); - data.gradients_quad.resize(n_comp * n_q_points * dim); - data.scratch.resize(2 * std::max(n_q_points, n_shape_values)); - - if (evaluation_flag & EvaluationFlags::hessians) - data.hessians_quad.resize(n_comp * n_q_points * n_hessians); + eval.set_data_pointers(&data.scratch, n_comp); const std::vector &renumber_to_lexicographic = data.shape_info.lexicographic_numbering; @@ -1146,20 +1150,14 @@ namespace internal { const unsigned int in_comp = d % n_lanes; const unsigned int out_comp = d / n_lanes; - data.values_dofs[out_comp * n_shape_values + i][in_comp] = + eval + .begin_dof_values()[out_comp * n_shape_values + i][in_comp] = data.mapping_support_points[renumber_to_lexicographic[i]][d]; } // do the actual tensorized evaluation internal::FEEvaluationFactory:: - evaluate(n_comp, - evaluation_flag, - data.shape_info, - data.values_dofs.begin(), - data.values_quad.begin(), - data.gradients_quad.begin(), - data.hessians_quad.begin(), - data.scratch.begin()); + evaluate(n_comp, evaluation_flag, eval.begin_dof_values(), eval); } // do the postprocessing @@ -1171,7 +1169,7 @@ namespace internal in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes; ++in_comp) quadrature_points[i][out_comp * n_lanes + in_comp] = - data.values_quad[out_comp * n_q_points + i][in_comp]; + eval.begin_values()[out_comp * n_q_points + i][in_comp]; } if (evaluation_flag & EvaluationFlags::gradients) @@ -1193,9 +1191,9 @@ namespace internal const unsigned int new_point = total_number % n_q_points; data.contravariant[new_point][out_comp * n_lanes + in_comp] [new_comp] = - data - .gradients_quad[(out_comp * n_q_points + point) * dim + - j][in_comp]; + eval.begin_gradients()[(out_comp * n_q_points + point) * + dim + + j][in_comp]; } } if (update_flags & update_covariant_transformation) @@ -1236,9 +1234,9 @@ namespace internal dim == 2 ? desymmetrize_2d[new_hessian_comp][1] : desymmetrize_3d[new_hessian_comp][1]; const double value = - data.hessians_quad[(out_comp * n_q_points + point) * - n_hessians + - j][in_comp]; + eval.begin_hessians()[(out_comp * n_q_points + point) * + n_hessians + + j][in_comp]; jacobian_grads[new_point][out_comp * n_lanes + in_comp] [new_hessian_comp_i][new_hessian_comp_j] = value; diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 81ff490f76..92001b04db 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -121,26 +121,21 @@ namespace internal n_q_points_1d, Number>; + using EvaluationData = + FEEvaluationBaseData; + static void - evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data); + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs_actual, + EvaluationData & eval); static void - integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array); + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs_actual, + EvaluationData & eval, + const bool add_into_values_array); static Eval create_evaluator_tensor_product( @@ -175,26 +170,21 @@ namespace internal n_q_points_1d, Number> { + using EvaluationData = + FEEvaluationBaseData; + static void - evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data); + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs_actual, + EvaluationData & eval); static void - integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array); + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs_actual, + EvaluationData & eval, + const bool add_into_values_array); }; @@ -206,14 +196,10 @@ namespace internal typename Number> inline void FEEvaluationImpl::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs_actual, + EvaluationData & eval) { if (evaluation_flag == EvaluationFlags::nothing) return; @@ -221,11 +207,13 @@ namespace internal std::array *, 3> univariate_shape_data; - univariate_shape_data.fill(&shape_info.data[0]); + const auto &shape_data = eval.get_shape_info().data; - if (shape_info.data.size() == dim) + univariate_shape_data.fill(&shape_data[0]); + + if (shape_data.size() == dim) for (int i = 1; i < dim; ++i) - univariate_shape_data[i] = &shape_info.data[i]; + univariate_shape_data[i] = &shape_data[i]; Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]); Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]); @@ -237,34 +225,36 @@ namespace internal (Eval::n_rows_of_product > Eval::n_columns_of_product ? Eval::n_rows_of_product : Eval::n_columns_of_product); - Number *temp1 = scratch_data; + Number *temp1 = eval.get_scratch_data().begin(); Number *temp2; if (temp_size == 0) { - temp2 = temp1 + std::max(Utilities::fixed_power( - shape_info.data.front().fe_degree + 1), - Utilities::fixed_power( - shape_info.data.front().n_q_points_1d)); + temp2 = + temp1 + + std::max(Utilities::fixed_power(shape_data[0].fe_degree + 1), + Utilities::fixed_power(shape_data[0].n_q_points_1d)); } else { temp2 = temp1 + temp_size; } - const unsigned int n_q_points = - temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product; - const unsigned int dofs_per_comp = + const std::size_t n_q_points = temp_size == 0 ? + eval.get_shape_info().n_q_points : + Eval::n_columns_of_product; + const std::size_t dofs_per_comp = (type == MatrixFreeFunctions::truncated_tensor) ? - Utilities::fixed_power(shape_info.data.front().fe_degree + 1) : - shape_info.dofs_per_component_on_cell; + Utilities::pow(shape_data[0].fe_degree + 1, dim) : + eval.get_shape_info().dofs_per_component_on_cell; const Number *values_dofs = values_dofs_actual; if (type == MatrixFreeFunctions::truncated_tensor) { + const std::size_t stored_dofs_per_comp = + eval.get_shape_info().dofs_per_component_on_cell; Number *values_dofs_tmp = - scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell, - shape_info.n_q_points)); + temp1 + 2 * (std::max(stored_dofs_per_comp, n_q_points)); const int degree = - fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; + fe_degree != -1 ? fe_degree : shape_data[0].fe_degree; for (unsigned int c = 0; c < n_components; ++c) for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1); @@ -275,8 +265,7 @@ namespace internal for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q) values_dofs_tmp[c * dofs_per_comp + count_q] = - values_dofs_actual - [c * shape_info.dofs_per_component_on_cell + count_p]; + values_dofs_actual[c * stored_dofs_per_comp + count_p]; for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q) values_dofs_tmp[c * dofs_per_comp + count_q] = Number(); @@ -288,6 +277,10 @@ namespace internal values_dofs = values_dofs_tmp; } + Number *values_quad = eval.begin_values(); + Number *gradients_quad = eval.begin_gradients(); + Number *hessians_quad = eval.begin_hessians(); + switch (dim) { case 1: @@ -460,10 +453,10 @@ namespace internal { values_quad -= n_components * n_q_points; values_dofs -= n_components * dofs_per_comp; - for (unsigned int c = 0; c < n_components; ++c) - for (unsigned int q = 0; q < shape_info.n_q_points; ++q) - values_quad[c * shape_info.n_q_points + q] += - values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1]; + for (std::size_t c = 0; c < n_components; ++c) + for (std::size_t q = 0; q < n_q_points; ++q) + values_quad[c * n_q_points + q] += + values_dofs[(c + 1) * dofs_per_comp - 1]; } } @@ -476,24 +469,21 @@ namespace internal typename Number> inline void FEEvaluationImpl::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs_actual, + EvaluationData & eval, + const bool add_into_values_array) { std::array *, 3> univariate_shape_data; - univariate_shape_data.fill(&shape_info.data[0]); + const auto &shape_data = eval.get_shape_info().data; + univariate_shape_data.fill(&shape_data[0]); - if (shape_info.data.size() == dim) + if (shape_data.size() == dim) for (int i = 1; i < dim; ++i) - univariate_shape_data[i] = &shape_info.data[i]; + univariate_shape_data[i] = &shape_data[i]; Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]); Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]); @@ -505,33 +495,38 @@ namespace internal (Eval::n_rows_of_product > Eval::n_columns_of_product ? Eval::n_rows_of_product : Eval::n_columns_of_product); - Number *temp1 = scratch_data; + Number *temp1 = eval.get_scratch_data().begin(); Number *temp2; if (temp_size == 0) { - temp2 = temp1 + std::max(Utilities::fixed_power( - shape_info.data.front().fe_degree + 1), - Utilities::fixed_power( - shape_info.data.front().n_q_points_1d)); + temp2 = + temp1 + + std::max(Utilities::fixed_power(shape_data[0].fe_degree + 1), + Utilities::fixed_power(shape_data[0].n_q_points_1d)); } else { temp2 = temp1 + temp_size; } - const unsigned int n_q_points = - temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product; + const unsigned int n_q_points = temp_size == 0 ? + eval.get_shape_info().n_q_points : + Eval::n_columns_of_product; const unsigned int dofs_per_comp = (type == MatrixFreeFunctions::truncated_tensor) ? - Utilities::fixed_power(shape_info.data.front().fe_degree + 1) : - shape_info.dofs_per_component_on_cell; + Utilities::fixed_power(shape_data[0].fe_degree + 1) : + eval.get_shape_info().dofs_per_component_on_cell; // expand dof_values to tensor product for truncated tensor products Number *values_dofs = (type == MatrixFreeFunctions::truncated_tensor) ? - scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell, - shape_info.n_q_points)) : + temp1 + 2 * (std::max(eval.get_shape_info().dofs_per_component_on_cell, + n_q_points)) : values_dofs_actual; + Number *values_quad = eval.begin_values(); + Number *gradients_quad = eval.begin_gradients(); + Number *hessians_quad = eval.begin_hessians(); + switch (dim) { case 1: @@ -734,14 +729,13 @@ namespace internal // case FE_Q_DG0: add values, gradients and second derivatives are zero if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0) { - values_dofs -= n_components * dofs_per_comp - - shape_info.dofs_per_component_on_cell + 1; + values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1; values_quad -= n_components * n_q_points; if (integration_flag & EvaluationFlags::values) for (unsigned int c = 0; c < n_components; ++c) { values_dofs[0] = values_quad[0]; - for (unsigned int q = 1; q < shape_info.n_q_points; ++q) + for (unsigned int q = 1; q < n_q_points; ++q) values_dofs[0] += values_quad[q]; values_dofs += dofs_per_comp; values_quad += n_q_points; @@ -749,16 +743,18 @@ namespace internal else { for (unsigned int c = 0; c < n_components; ++c) - values_dofs[c * shape_info.dofs_per_component_on_cell] = Number(); - values_dofs += n_components * shape_info.dofs_per_component_on_cell; + values_dofs[c * dofs_per_comp] = Number(); + values_dofs += n_components * dofs_per_comp; } } if (type == MatrixFreeFunctions::truncated_tensor) { + const std::size_t stored_dofs_per_comp = + eval.get_shape_info().dofs_per_component_on_cell; values_dofs -= dofs_per_comp * n_components; const int degree = - fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; + fe_degree != -1 ? fe_degree : shape_data[0].fe_degree; for (unsigned int c = 0; c < n_components; ++c) for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1); @@ -768,9 +764,7 @@ namespace internal { for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q) - values_dofs_actual[c * - shape_info.dofs_per_component_on_cell + - count_p] = + values_dofs_actual[c * stored_dofs_per_comp + count_p] = values_dofs[c * dofs_per_comp + count_q]; count_q += j + i; } @@ -783,32 +777,29 @@ namespace internal template inline void - FEEvaluationImpl< - MatrixFreeFunctions::tensor_none, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number *values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data) + FEEvaluationImpl::evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags + evaluation_flag, + const Number * values_dofs_actual, + EvaluationData &eval) { - (void)scratch_data; + const unsigned int n_dofs = + eval.get_shape_info().dofs_per_component_on_cell; + const unsigned int n_q_points = eval.get_shape_info().n_q_points; - const unsigned int n_dofs = shape_info.dofs_per_component_on_cell; - const unsigned int n_q_points = shape_info.n_q_points; + const auto &shape_data = eval.get_shape_info().data; using Eval = EvaluatorTensorProduct; if (evaluation_flag & EvaluationFlags::values) { - const auto shape_values = shape_info.data.front().shape_values.data(); - auto values_quad_ptr = values_quad; + const auto shape_values = shape_data[0].shape_values.data(); + auto values_quad_ptr = eval.begin_values(); auto values_dofs_actual_ptr = values_dofs_actual; Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); @@ -824,10 +815,9 @@ namespace internal if (evaluation_flag & EvaluationFlags::gradients) { - const auto shape_gradients = - shape_info.data.front().shape_gradients.data(); - auto gradients_quad_ptr = gradients_quad; - auto values_dofs_actual_ptr = values_dofs_actual; + const auto shape_gradients = shape_data[0].shape_gradients.data(); + auto gradients_quad_ptr = eval.begin_gradients(); + auto values_dofs_actual_ptr = values_dofs_actual; for (unsigned int c = 0; c < n_components; ++c) { @@ -849,47 +839,41 @@ namespace internal } if (evaluation_flag & EvaluationFlags::hessians) - { - Assert(false, ExcNotImplemented()); - (void)hessians_quad; - } + Assert(false, ExcNotImplemented()); } template inline void - FEEvaluationImpl< - MatrixFreeFunctions::tensor_none, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array) + FEEvaluationImpl::integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags + integration_flag, + Number * values_dofs_actual, + EvaluationData &eval, + const bool add_into_values_array) { // TODO: implement hessians - (void)hessians_quad; AssertThrow(!(integration_flag & EvaluationFlags::hessians), ExcNotImplemented()); - (void)scratch_data; - const unsigned int n_dofs = shape_info.dofs_per_component_on_cell; - const unsigned int n_q_points = shape_info.n_q_points; + const unsigned int n_dofs = + eval.get_shape_info().dofs_per_component_on_cell; + const unsigned int n_q_points = eval.get_shape_info().n_q_points; + + const auto &shape_data = eval.get_shape_info().data; using Eval = EvaluatorTensorProduct; if (integration_flag & EvaluationFlags::values) { - const auto shape_values = shape_info.data.front().shape_values.data(); - auto values_quad_ptr = values_quad; + const auto shape_values = shape_data[0].shape_values.data(); + auto values_quad_ptr = eval.begin_values(); auto values_dofs_actual_ptr = values_dofs_actual; Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); @@ -909,10 +893,9 @@ namespace internal if (integration_flag & EvaluationFlags::gradients) { - const auto shape_gradients = - shape_info.data.front().shape_gradients.data(); - auto gradients_quad_ptr = gradients_quad; - auto values_dofs_actual_ptr = values_dofs_actual; + const auto shape_gradients = shape_data[0].shape_gradients.data(); + auto gradients_quad_ptr = eval.begin_gradients(); + auto values_dofs_actual_ptr = values_dofs_actual; for (unsigned int c = 0; c < n_components; ++c) { @@ -1357,26 +1340,36 @@ namespace internal template struct FEEvaluationImplCollocation { + using EvaluationData = + FEEvaluationBaseData; + + static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + EvaluationData & eval); + + static void + do_evaluate(const MatrixFreeFunctions::UnivariateShapeData &shape, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + Number * gradients_quad, + Number * hessians_quad); + static void - evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data); + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + EvaluationData & eval, + const bool add_into_values_array); static void - integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array); + do_integrate(const MatrixFreeFunctions::UnivariateShapeData &shape, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + Number * gradients_quad, + const Number * hessians_quad, + const bool add_into_values_array); }; @@ -1384,76 +1377,86 @@ namespace internal template inline void FEEvaluationImplCollocation::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number *) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + EvaluationData & eval) + { + constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); + + for (unsigned int c = 0; c < n_components; ++c) + { + if (evaluation_flag & EvaluationFlags::values) + for (unsigned int i = 0; i < n_points; ++i) + eval.begin_values()[n_points * c + i] = + values_dofs[n_points * c + i]; + + do_evaluate(eval.get_shape_info().data[0], + evaluation_flag, + values_dofs + c * n_points, + eval.begin_gradients() + c * dim * n_points, + eval.begin_hessians() + c * dim * (dim + 1) / 2 * n_points); + } + } + + + + template + inline void + FEEvaluationImplCollocation::do_evaluate( + const MatrixFreeFunctions::UnivariateShapeData &shape, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + Number * gradients_quad, + Number * hessians_quad) { - AssertDimension( - shape_info.data.front().shape_gradients_collocation_eo.size(), - (fe_degree + 2) / 2 * (fe_degree + 1)); + AssertDimension(shape.shape_gradients_collocation_eo.size(), + (fe_degree + 2) / 2 * (fe_degree + 1)); + constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); EvaluatorTensorProduct - eval(AlignedVector(), - shape_info.data.front().shape_gradients_collocation_eo, - shape_info.data.front().shape_hessians_collocation_eo); - constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim); - - for (unsigned int c = 0; c < n_components; ++c) + eval(AlignedVector(), + shape.shape_gradients_collocation_eo, + shape.shape_hessians_collocation_eo); + if (evaluation_flag & + (EvaluationFlags::gradients | EvaluationFlags::hessians)) { - if (evaluation_flag & EvaluationFlags::values) - for (unsigned int i = 0; i < n_q_points; ++i) - values_quad[i] = values_dofs[i]; - if (evaluation_flag & - (EvaluationFlags::gradients | EvaluationFlags::hessians)) + eval.template gradients<0, true, false>(values_dofs, gradients_quad); + if (dim > 1) + eval.template gradients<1, true, false>(values_dofs, + gradients_quad + n_points); + if (dim > 2) + eval.template gradients<2, true, false>(values_dofs, + gradients_quad + + 2 * n_points); + } + if (evaluation_flag & EvaluationFlags::hessians) + { + eval.template hessians<0, true, false>(values_dofs, hessians_quad); + if (dim > 1) { - eval.template gradients<0, true, false>(values_dofs, - gradients_quad); - if (dim > 1) - eval.template gradients<1, true, false>(values_dofs, - gradients_quad + - n_q_points); - if (dim > 2) - eval.template gradients<2, true, false>(values_dofs, - gradients_quad + - 2 * n_q_points); + eval.template gradients<1, true, false>(gradients_quad, + hessians_quad + + dim * n_points); + eval.template hessians<1, true, false>(values_dofs, + hessians_quad + n_points); } - if (evaluation_flag & EvaluationFlags::hessians) + if (dim > 2) { - eval.template hessians<0, true, false>(values_dofs, hessians_quad); - if (dim > 1) - { - eval.template gradients<1, true, false>(gradients_quad, - hessians_quad + - dim * n_q_points); - eval.template hessians<1, true, false>(values_dofs, - hessians_quad + - n_q_points); - } - if (dim > 2) - { - eval.template gradients<2, true, false>(gradients_quad, - hessians_quad + - 4 * n_q_points); - eval.template gradients<2, true, false>( - gradients_quad + n_q_points, hessians_quad + 5 * n_q_points); - eval.template hessians<2, true, false>(values_dofs, - hessians_quad + - 2 * n_q_points); - } - hessians_quad += (dim * (dim + 1)) / 2 * n_q_points; + eval.template gradients<2, true, false>(gradients_quad, + hessians_quad + + 4 * n_points); + eval.template gradients<2, true, false>(gradients_quad + n_points, + hessians_quad + + 5 * n_points); + eval.template hessians<2, true, false>(values_dofs, + hessians_quad + + 2 * n_points); } - gradients_quad += dim * n_q_points; - values_quad += n_q_points; - values_dofs += n_q_points; } } @@ -1462,117 +1465,140 @@ namespace internal template inline void FEEvaluationImplCollocation::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + EvaluationData & eval, + const bool add_into_values_array) + { + constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); + + for (unsigned int c = 0; c < n_components; ++c) + { + if (integration_flag & EvaluationFlags::values) + { + if (add_into_values_array) + for (unsigned int i = 0; i < n_points; ++i) + values_dofs[n_points * c + i] += + eval.begin_values()[n_points * c + i]; + else + for (unsigned int i = 0; i < n_points; ++i) + values_dofs[n_points * c + i] = + eval.begin_values()[n_points * c + i]; + } + + do_integrate(eval.get_shape_info().data[0], + integration_flag, + values_dofs + c * n_points, + eval.begin_gradients() + c * dim * n_points, + eval.begin_hessians() + c * dim * (dim + 1) / 2 * n_points, + add_into_values_array || + (integration_flag & EvaluationFlags::values)); + } + } + + + + template + inline void + FEEvaluationImplCollocation::do_integrate( + const MatrixFreeFunctions::UnivariateShapeData &shape, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + Number * gradients_quad, + const Number * hessians_quad, + const bool add_into_values_array) { - AssertDimension( - shape_info.data.front().shape_gradients_collocation_eo.size(), - (fe_degree + 2) / 2 * (fe_degree + 1)); + AssertDimension(shape.shape_gradients_collocation_eo.size(), + (fe_degree + 2) / 2 * (fe_degree + 1)); EvaluatorTensorProduct - eval(AlignedVector(), - shape_info.data.front().shape_gradients_collocation_eo, - shape_info.data.front().shape_hessians_collocation_eo); - constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim); - constexpr unsigned int hdim = (dim * (dim + 1)) / 2; + eval(AlignedVector(), + shape.shape_gradients_collocation_eo, + shape.shape_hessians_collocation_eo); + constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); - for (unsigned int c = 0; c < n_components; ++c) + if (integration_flag & EvaluationFlags::hessians) { - if (integration_flag & EvaluationFlags::values) + // diagonal + // grad xx + if (add_into_values_array == true) + eval.template hessians<0, false, true>(hessians_quad, values_dofs); + else + eval.template hessians<0, false, false>(hessians_quad, values_dofs); + // grad yy + if (dim > 1) + eval.template hessians<1, false, true>(hessians_quad + n_points, + values_dofs); + // grad zz + if (dim > 2) + eval.template hessians<2, false, true>(hessians_quad + 2 * n_points, + values_dofs); + // off-diagonal + if (dim == 2) { - if (add_into_values_array == false) - for (unsigned int i = 0; i < n_q_points; ++i) - values_dofs[i] = values_quad[i]; + // grad xy, queue into gradient + if (integration_flag & EvaluationFlags::gradients) + eval.template gradients<1, false, true>(hessians_quad + + 2 * n_points, + gradients_quad); else - for (unsigned int i = 0; i < n_q_points; ++i) - values_dofs[i] += values_quad[i]; + eval.template gradients<1, false, false>(hessians_quad + + 2 * n_points, + gradients_quad); } - if (integration_flag & EvaluationFlags::gradients) + if (dim == 3) { - if (integration_flag & EvaluationFlags::values || - add_into_values_array == true) - eval.template gradients<0, false, true>(gradients_quad, - values_dofs); + // grad xy, queue into gradient + if (integration_flag & EvaluationFlags::gradients) + eval.template gradients<1, false, true>(hessians_quad + + 3 * n_points, + gradients_quad); else - eval.template gradients<0, false, false>(gradients_quad, - values_dofs); - if (dim > 1) - eval.template gradients<1, false, true>(gradients_quad + - n_q_points, - values_dofs); - if (dim > 2) - eval.template gradients<2, false, true>(gradients_quad + - 2 * n_q_points, - values_dofs); - } - if (integration_flag & EvaluationFlags::hessians) - { - // diagonal - // grad xx - if (integration_flag & EvaluationFlags::values || - integration_flag & EvaluationFlags::gradients || - add_into_values_array == true) - eval.template hessians<0, false, true>(hessians_quad, - values_dofs); + eval.template gradients<1, false, false>(hessians_quad + + 3 * n_points, + gradients_quad); + + // grad xz + eval.template gradients<2, false, true>(hessians_quad + + 4 * n_points, + gradients_quad); + + // grad yz + if (integration_flag & EvaluationFlags::gradients) + eval.template gradients<2, false, true>( + hessians_quad + 5 * n_points, gradients_quad + n_points); else - eval.template hessians<0, false, false>(hessians_quad, - values_dofs); - // grad yy - if (dim > 1) - eval.template hessians<1, false, true>(hessians_quad + n_q_points, - values_dofs); - // grad zz - if (dim > 2) - eval.template hessians<2, false, true>(hessians_quad + - 2 * n_q_points, - values_dofs); - // off-diagonal - if (dim == 2) - { - // grad xy - eval.template gradients<0, false, false>(hessians_quad + - 2 * n_q_points, - scratch_data); - eval.template gradients<1, false, true>(scratch_data, - values_dofs); - } - if (dim == 3) - { - // grad xy - eval.template gradients<0, false, false>(hessians_quad + - 3 * n_q_points, - scratch_data); - eval.template gradients<1, false, true>(scratch_data, - values_dofs); - // grad xz - eval.template gradients<0, false, false>(hessians_quad + - 4 * n_q_points, - scratch_data); - eval.template gradients<2, false, true>(scratch_data, - values_dofs); - // grad yz - eval.template gradients<1, false, false>(hessians_quad + - 5 * n_q_points, - scratch_data); - eval.template gradients<2, false, true>(scratch_data, - values_dofs); - } - hessians_quad += hdim * n_q_points; + eval.template gradients<2, false, false>( + hessians_quad + 5 * n_points, gradients_quad + n_points); } - gradients_quad += dim * n_q_points; - values_quad += n_q_points; - values_dofs += n_q_points; + + // if we did not integrate gradients, set the last slot to zero + // which was not touched before, in order to avoid the if + // statement in the gradients loop below + if (!(integration_flag & EvaluationFlags::gradients)) + for (unsigned int q = 0; q < n_points; ++q) + gradients_quad[(dim - 1) * n_points + q] = Number(); + } + + if (integration_flag & + (EvaluationFlags::gradients | EvaluationFlags::hessians)) + { + if (add_into_values_array || + integration_flag & EvaluationFlags::hessians) + eval.template gradients<0, false, true>(gradients_quad, values_dofs); + else + eval.template gradients<0, false, false>(gradients_quad, values_dofs); + if (dim > 1) + eval.template gradients<1, false, true>(gradients_quad + n_points, + values_dofs); + if (dim > 2) + eval.template gradients<2, false, true>(gradients_quad + 2 * n_points, + values_dofs); } } @@ -1591,26 +1617,21 @@ namespace internal template struct FEEvaluationImplTransformToCollocation { + using EvaluationData = + FEEvaluationBaseData; + static void - evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data); + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + EvaluationData & eval); static void - integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool add_into_values_array); + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs, + EvaluationData & eval, + const bool add_into_values_array); }; @@ -1623,19 +1644,17 @@ namespace internal n_q_points_1d, Number>::evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * values_dofs, - Number * values_quad, - Number *gradients_quad, - Number *hessians_quad, - Number *) + const Number * values_dofs, + EvaluationData & eval) { + const auto &shape_data = eval.get_shape_info().data[0]; + Assert(n_q_points_1d > fe_degree, ExcMessage("You lose information when going to a collocation space " "of lower degree, so the evaluation results would be " "wrong. Thus, this class does not permit the desired " "operation.")); - constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim); + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim); for (unsigned int c = 0; c < n_components; ++c) { @@ -1647,28 +1666,22 @@ namespace internal n_q_points_1d, Number, Number>::do_forward(1, - shape_info.data.front().shape_values_eo, - values_dofs, - values_quad); + shape_data.shape_values_eo, + values_dofs + + c * Utilities::pow(fe_degree + 1, dim), + eval.begin_values() + c * n_q_points); // apply derivatives in the collocation space if (evaluation_flag & (EvaluationFlags::gradients | EvaluationFlags::hessians)) - FEEvaluationImplCollocation::evaluate( - 1, - evaluation_flag & - (EvaluationFlags::gradients | EvaluationFlags::hessians), - shape_info, - values_quad, - nullptr, - gradients_quad, - hessians_quad, - nullptr); - - values_dofs += shape_info.dofs_per_component_on_cell; - values_quad += n_q_points; - gradients_quad += dim * n_q_points; - hessians_quad += (dim * (dim + 1)) / 2 * n_q_points; + FEEvaluationImplCollocation:: + do_evaluate(shape_data, + evaluation_flag & (EvaluationFlags::gradients | + EvaluationFlags::hessians), + eval.begin_values() + c * n_q_points, + eval.begin_gradients() + c * dim * n_q_points, + eval.begin_hessians() + + c * dim * (dim + 1) / 2 * n_q_points); } } @@ -1682,24 +1695,18 @@ namespace internal n_q_points_1d, Number>::integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, + Number * values_dofs, + EvaluationData & eval, const bool add_into_values_array) { + const auto &shape_data = eval.get_shape_info().data[0]; + Assert(n_q_points_1d > fe_degree, ExcMessage("You lose information when going to a collocation space " "of lower degree, so the evaluation results would be " "wrong. Thus, this class does not permit the desired " "operation.")); - AssertDimension( - shape_info.data.front().shape_gradients_collocation_eo.size(), - (n_q_points_1d + 1) / 2 * n_q_points_1d); constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim); - constexpr unsigned int hdim = (dim * (dim + 1)) / 2; for (unsigned int c = 0; c < n_components; ++c) { @@ -1707,17 +1714,15 @@ namespace internal if (integration_flag & (EvaluationFlags::gradients | EvaluationFlags::hessians)) FEEvaluationImplCollocation:: - integrate(1, - integration_flag & (EvaluationFlags::gradients | - EvaluationFlags::hessians), - shape_info, - values_quad, - nullptr, - gradients_quad, - hessians_quad, - scratch_data, - /*add_into_values_array=*/integration_flag & - EvaluationFlags::values); + do_integrate(shape_data, + integration_flag & (EvaluationFlags::gradients | + EvaluationFlags::hessians), + eval.begin_values() + c * n_q_points, + eval.begin_gradients() + c * dim * n_q_points, + eval.begin_hessians() + + c * dim * (dim + 1) / 2 * n_q_points, + /*add_into_values_array=*/ + integration_flag & EvaluationFlags::values); // transform back to the original space FEEvaluationImplBasisChange< @@ -1728,14 +1733,11 @@ namespace internal n_q_points_1d, Number, Number>::do_backward(1, - shape_info.data.front().shape_values_eo, + shape_data.shape_values_eo, add_into_values_array, - values_quad, - values_dofs); - hessians_quad += hdim * n_q_points; - gradients_quad += dim * n_q_points; - values_quad += n_q_points; - values_dofs += shape_info.dofs_per_component_on_cell; + eval.begin_values() + c * n_q_points, + values_dofs + + c * Utilities::pow(fe_degree + 1, dim)); } } @@ -1761,146 +1763,100 @@ namespace internal { template static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const internal::MatrixFreeFunctions::ShapeInfo &shape_info, - Number *values_dofs_actual, - Number *values_quad, - Number *gradients_quad, - Number *hessians_quad, - Number *scratch_data) + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + FEEvaluationBaseData + &eval) { + const auto element_type = eval.get_shape_info().element_type; + using ElementType = internal::MatrixFreeFunctions::ElementType; + // We enable a transformation to collocation for derivatives if it gives // correct results (first condition), if it is the most efficient choice // in terms of operation counts (second condition) and if we were able to // initialize the fields in shape_info.templates.h from the polynomials // (third condition). - static constexpr bool use_collocation = - n_q_points_1d > fe_degree && n_q_points_1d <= 3 * fe_degree / 2 + 1 && - n_q_points_1d < 200; - - Assert(shape_info.data.size() == 1 || - (shape_info.data.size() == dim && - shape_info.element_type == - internal::MatrixFreeFunctions::tensor_general), + constexpr bool use_collocation = (n_q_points_1d > fe_degree) && + (n_q_points_1d < 200) && + (n_q_points_1d <= 3 * fe_degree / 2 + 1); + + Assert(eval.get_shape_info().data.size() == 1 || + (eval.get_shape_info().data.size() == dim && + element_type == ElementType::tensor_general), ExcNotImplemented()); if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && - shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_collocation) + element_type == ElementType::tensor_symmetric_collocation) { internal::FEEvaluationImplCollocation:: - evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + evaluate(n_components, evaluation_flag, values_dofs, eval); } // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see // shape_info.h for more details else if (fe_degree >= 0 && use_collocation && - shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) + element_type <= ElementType::tensor_symmetric) { internal::FEEvaluationImplTransformToCollocation< dim, fe_degree, n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + Number>::evaluate(n_components, evaluation_flag, values_dofs, eval); } - else if (fe_degree >= 0 && - shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) + else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric) { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + internal::FEEvaluationImpl::evaluate(n_components, + evaluation_flag, + values_dofs, + eval); } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + else if (element_type == ElementType::tensor_symmetric_plus_dg0) { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + internal::FEEvaluationImpl::evaluate(n_components, + evaluation_flag, + values_dofs, + eval); } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::truncated_tensor) + else if (element_type == ElementType::truncated_tensor) { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::truncated_tensor, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + internal::FEEvaluationImpl::evaluate(n_components, + evaluation_flag, + values_dofs, + eval); } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_none) + else if (element_type == ElementType::tensor_none) { - internal::FEEvaluationImpl::evaluate(n_components, evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + values_dofs, + eval); } else { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_general, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + internal::FEEvaluationImpl::evaluate(n_components, + evaluation_flag, + values_dofs, + eval); } return false; @@ -1931,49 +1887,42 @@ namespace internal static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - const internal::MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, + Number * values_dofs, + FEEvaluationBaseData + & eval, const bool sum_into_values_array) { + const auto element_type = eval.get_shape_info().element_type; + using ElementType = internal::MatrixFreeFunctions::ElementType; + // We enable a transformation to collocation for derivatives if it gives // correct results (first condition), if it is the most efficient choice // in terms of operation counts (second condition) and if we were able to // initialize the fields in shape_info.templates.h from the polynomials // (third condition). - constexpr bool use_collocation = n_q_points_1d > fe_degree && - n_q_points_1d <= 3 * fe_degree / 2 + 1 && - n_q_points_1d < 200; - - Assert(shape_info.data.size() == 1 || - (shape_info.data.size() == dim && - shape_info.element_type == - internal::MatrixFreeFunctions::tensor_general), + constexpr bool use_collocation = (n_q_points_1d > fe_degree) && + (n_q_points_1d < 200) && + (n_q_points_1d <= 3 * fe_degree / 2 + 1); + + Assert(eval.get_shape_info().data.size() == 1 || + (eval.get_shape_info().data.size() == dim && + element_type == ElementType::tensor_general), ExcNotImplemented()); if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && - shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_collocation) + element_type == ElementType::tensor_symmetric_collocation) { internal::FEEvaluationImplCollocation:: integrate(n_components, integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, + values_dofs, + eval, sum_into_values_array); } // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see // shape_info.h for more details else if (fe_degree >= 0 && use_collocation && - shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) + element_type <= ElementType::tensor_symmetric) { internal::FEEvaluationImplTransformToCollocation< dim, @@ -1981,102 +1930,69 @@ namespace internal n_q_points_1d, Number>::integrate(n_components, integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, + values_dofs, + eval, sum_into_values_array); } - else if (fe_degree >= 0 && - shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) + else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric) { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - sum_into_values_array); + internal::FEEvaluationImpl::integrate(n_components, + integration_flag, + values_dofs, + eval, + sum_into_values_array); } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + else if (element_type == ElementType::tensor_symmetric_plus_dg0) { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - sum_into_values_array); + internal::FEEvaluationImpl::integrate(n_components, + integration_flag, + values_dofs, + eval, + sum_into_values_array); } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::truncated_tensor) + else if (element_type == ElementType::truncated_tensor) { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::truncated_tensor, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - sum_into_values_array); + internal::FEEvaluationImpl::integrate(n_components, + integration_flag, + values_dofs, + eval, + sum_into_values_array); } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_none) + else if (element_type == ElementType::tensor_none) { - internal::FEEvaluationImpl::integrate(n_components, integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, + values_dofs, + eval, sum_into_values_array); } else { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_general, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - sum_into_values_array); + internal::FEEvaluationImpl::integrate(n_components, + integration_flag, + values_dofs, + eval, + sum_into_values_array); } return false; @@ -2098,96 +2014,76 @@ namespace internal // able to initialize the fields in shape_info.templates.h from the // polynomials (fourth condition). static constexpr bool use_collocation = - symmetric_evaluate && - n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 && - n_q_points_1d < 200; + symmetric_evaluate && (n_q_points_1d > fe_degree) && + (n_q_points_1d < 200) && (n_q_points_1d <= 3 * fe_degree / 2 + 1); + + using Eval = + EvaluatorTensorProduct; + + using EvaluationData = + FEEvaluationBaseData; + + static Eval + create_evaluator_tensor_product( + const internal::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 MatrixFreeFunctions::ShapeInfo &data, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_hessian, - const unsigned int subface_index) + 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) { - const AlignedVector &val0 = - symmetric_evaluate ? - data.data.front().shape_values_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_values : - data.data.front().values_within_subface[subface_index % 2]); - const AlignedVector &val1 = - symmetric_evaluate ? - data.data.front().shape_values_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_values : - data.data.front().values_within_subface[subface_index / 2]); - - const AlignedVector &grad0 = - symmetric_evaluate ? - data.data.front().shape_gradients_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_gradients : - data.data.front().gradients_within_subface[subface_index % 2]); - const AlignedVector &grad1 = - symmetric_evaluate ? - data.data.front().shape_gradients_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_gradients : - data.data.front().gradients_within_subface[subface_index / 2]); - - const AlignedVector &hessian0 = - symmetric_evaluate ? - data.data.front().shape_hessians_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_hessians : - data.data.front().hessians_within_subface[subface_index % 2]); - const AlignedVector &hessian1 = - symmetric_evaluate ? - data.data.front().shape_hessians_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_hessians : - data.data.front().hessians_within_subface[subface_index / 2]); - - using Eval = - internal::EvaluatorTensorProduct; - Eval eval0(val0, - grad0, - hessian0, - data.data.front().fe_degree + 1, - data.data.front().n_q_points_1d); - Eval eval1(val1, - grad1, - hessian1, - data.data.front().fe_degree + 1, - data.data.front().n_q_points_1d); - - const unsigned int size_deg = - fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - (dim > 1 ? - Utilities::fixed_power(data.data.front().fe_degree + 1) : - 1); + Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); + Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); - const unsigned int n_q_points = fe_degree > -1 ? - Utilities::pow(n_q_points_1d, dim - 1) : - data.n_q_points_face; + const std::size_t size_deg = + 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 (evaluate_val && !evaluate_grad) + if (evaluation_flag == EvaluationFlags::values) for (unsigned int c = 0; c < n_components; ++c) { switch (dim) @@ -2211,7 +2107,7 @@ namespace internal values_dofs += 3 * size_deg; values_quad += n_q_points; } - else if (evaluate_grad) + else if (evaluation_flag & EvaluationFlags::gradients) for (unsigned int c = 0; c < n_components; ++c) { switch (dim) @@ -2229,10 +2125,9 @@ namespace internal n_q_points_1d, n_q_points_1d, Number> - eval_grad( - AlignedVector(), - data.data.front().shape_gradients_collocation_eo, - AlignedVector()); + eval_grad(AlignedVector(), + data.shape_gradients_collocation_eo, + AlignedVector()); eval_grad.template gradients<0, true, false>( values_quad, gradients_quad); eval_grad.template gradients<1, true, false>( @@ -2253,7 +2148,7 @@ namespace internal gradients_quad + n_q_points); - if (evaluate_val == true) + if (evaluation_flag & EvaluationFlags::values) eval1.template values<1, true, false>(scratch_data, values_quad); } @@ -2267,11 +2162,10 @@ namespace internal case 2: eval0.template values<0, true, false>(values_dofs + size_deg, gradients_quad + - (dim - 1) * - n_q_points); + n_q_points); eval0.template gradients<0, true, false>(values_dofs, gradients_quad); - if (evaluate_val == true) + if (evaluation_flag & EvaluationFlags::values) eval0.template values<0, true, false>(values_dofs, values_quad); break; @@ -2287,10 +2181,9 @@ namespace internal gradients_quad += dim * n_q_points; } - if (evaluate_hessian) + if (evaluation_flag & EvaluationFlags::hessians) { - values_dofs = values_dofs_ptr; - constexpr unsigned int hdim = (dim * (dim + 1)) / 2; + values_dofs = values_dofs_ptr; for (unsigned int c = 0; c < n_components; ++c) { switch (dim) @@ -2359,97 +2252,38 @@ namespace internal AssertThrow(false, ExcNotImplemented()); } values_dofs += 3 * size_deg; - hessians_quad += hdim * n_q_points; + hessians_quad += dim * (dim + 1) / 2 * n_q_points; } } } static void - integrate_in_face(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const bool integrate_val, - const bool integrate_grad, - const bool integrate_hessian, - const unsigned int subface_index) + 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) { - const AlignedVector &val0 = - symmetric_evaluate ? - data.data.front().shape_values_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_values : - data.data.front().values_within_subface[subface_index % 2]); - const AlignedVector &val1 = - symmetric_evaluate ? - data.data.front().shape_values_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_values : - data.data.front().values_within_subface[subface_index / 2]); - - const AlignedVector &grad0 = - symmetric_evaluate ? - data.data.front().shape_gradients_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_gradients : - data.data.front().gradients_within_subface[subface_index % 2]); - const AlignedVector &grad1 = - symmetric_evaluate ? - data.data.front().shape_gradients_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_gradients : - data.data.front().gradients_within_subface[subface_index / 2]); - - const AlignedVector &hessian0 = - symmetric_evaluate ? - data.data.front().shape_hessians_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_hessians : - data.data.front().hessians_within_subface[subface_index % 2]); - const AlignedVector &hessian1 = - symmetric_evaluate ? - data.data.front().shape_hessians_eo : - (subface_index >= GeometryInfo::max_children_per_cell ? - data.data.front().shape_hessians : - data.data.front().hessians_within_subface[subface_index / 2]); - - using Eval = - internal::EvaluatorTensorProduct; - Eval eval0(val0, - grad0, - hessian0, - data.data.front().fe_degree + 1, - data.data.front().n_q_points_1d); - Eval eval1(val1, - grad1, - hessian1, - data.data.front().fe_degree + 1, - data.data.front().n_q_points_1d); - - const unsigned int size_deg = + Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); + Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); + + const std::size_t size_deg = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - (dim > 1 ? - Utilities::fixed_power(data.data.front().fe_degree + 1) : - 1); - - const unsigned int n_q_points = fe_degree > -1 ? - Utilities::pow(n_q_points_1d, dim - 1) : - data.n_q_points_face; + (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 (integrate_val && !integrate_grad) + if (integration_flag == EvaluationFlags::values) for (unsigned int c = 0; c < n_components; ++c) { switch (dim) @@ -2473,7 +2307,7 @@ namespace internal values_dofs += 3 * size_deg; values_quad += n_q_points; } - else if (integrate_grad) + else if (integration_flag & EvaluationFlags::gradients) for (unsigned int c = 0; c < n_components; ++c) { switch (dim) @@ -2494,11 +2328,10 @@ namespace internal n_q_points_1d, n_q_points_1d, Number> - eval_grad( - AlignedVector(), - data.data.front().shape_gradients_collocation_eo, - AlignedVector()); - if (integrate_val) + eval_grad(AlignedVector(), + data.shape_gradients_collocation_eo, + AlignedVector()); + if (integration_flag & EvaluationFlags::values) eval_grad.template gradients<1, false, true>( gradients_quad + n_q_points, values_quad); else @@ -2513,7 +2346,7 @@ namespace internal } else { - if (integrate_val) + if (integration_flag & EvaluationFlags::values) { eval1.template values<1, false, false>(values_quad, scratch_data); @@ -2540,7 +2373,7 @@ namespace internal gradients_quad + n_q_points, values_dofs + size_deg); eval0.template gradients<0, false, false>(gradients_quad, values_dofs); - if (integrate_val == true) + if (integration_flag & EvaluationFlags::values) eval0.template values<0, false, true>(values_quad, values_dofs); break; @@ -2556,10 +2389,9 @@ namespace internal gradients_quad += dim * n_q_points; } - if (integrate_hessian) + if (integration_flag & EvaluationFlags::hessians) { - values_dofs = values_dofs_ptr; - constexpr unsigned int hdim = (dim * (dim + 1)) / 2; + values_dofs = values_dofs_ptr; for (unsigned int c = 0; c < n_components; ++c) { switch (dim) @@ -2568,7 +2400,8 @@ namespace internal // grad xx eval1.template values<1, false, false>(hessians_quad, scratch_data); - if (integrate_val || integrate_grad) + if (integration_flag & + (EvaluationFlags::values | EvaluationFlags::gradients)) eval0.template hessians<0, false, true>(scratch_data, values_dofs); else @@ -2601,7 +2434,7 @@ namespace internal eval1.template values<1, false, false>(hessians_quad + 4 * n_q_points, scratch_data); - if (integrate_grad) + if (integration_flag & EvaluationFlags::gradients) eval0.template gradients<0, false, true>(scratch_data, values_dofs + size_deg); @@ -2621,7 +2454,8 @@ namespace internal break; case 2: // grad xx - if (integrate_val || integrate_grad) + if (integration_flag & + (EvaluationFlags::values | EvaluationFlags::gradients)) eval0.template hessians<0, false, true>(hessians_quad, values_dofs); else @@ -2632,7 +2466,7 @@ namespace internal eval0.template values<0, false, false>( hessians_quad + n_q_points, values_dofs + 2 * size_deg); // grad xy - if (integrate_grad) + if (integration_flag & EvaluationFlags::gradients) eval0.template gradients<0, false, true>( hessians_quad + 2 * n_q_points, values_dofs + size_deg); else @@ -2641,16 +2475,16 @@ namespace internal break; case 1: values_dofs[2] = hessians_quad[0]; - if (!integrate_val) + if (!(integration_flag & EvaluationFlags::values)) values_dofs[0] = 0; - if (!integrate_grad) + if (!(integration_flag & EvaluationFlags::gradients)) values_dofs[1] = 0; break; default: AssertThrow(false, ExcNotImplemented()); } values_dofs += 3 * size_deg; - hessians_quad += hdim * n_q_points; + hessians_quad += dim * (dim + 1) / 2 * n_q_points; } } } @@ -2664,15 +2498,13 @@ namespace internal template static void interpolate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo &data, const Number * input, Number * output, - const bool do_gradients, - const bool do_hessians, const unsigned int face_no) { - Assert(static_cast(fe_degree) == - data.data.front().fe_degree || + Assert(static_cast(fe_degree) == data.data[0].fe_degree || fe_degree == -1, ExcInternalError()); @@ -2680,11 +2512,10 @@ namespace internal n_components, input, output, - do_gradients, - do_hessians, + flags, face_no, - data.data.front().fe_degree + 1, - data.data.front().shape_data_on_face, + data.data[0].fe_degree + 1, + data.data[0].shape_data_on_face, data.dofs_per_component_on_cell, 3 * data.dofs_per_component_on_face); } @@ -2694,16 +2525,15 @@ namespace internal */ template static void - interpolate_quadrature(const unsigned int n_components, + interpolate_quadrature(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo &data, const Number * input, Number * output, - const bool do_gradients, - const bool do_hessians, const unsigned int face_no) { Assert(static_cast(fe_degree + 1) == - data.data.front().quadrature.size() || + data.data[0].n_q_points_1d || fe_degree == -1, ExcInternalError()); @@ -2711,8 +2541,7 @@ namespace internal n_components, input, output, - do_gradients, - do_hessians, + flags, face_no, data.data.front().quadrature.size(), data.data.front().quadrature_data_on_face, @@ -2723,13 +2552,12 @@ namespace internal private: template static void - interpolate_generic(const unsigned int n_components, - const Number * input, - Number * output, - const bool do_gradients, - const bool do_hessians, - const unsigned int face_no, - const unsigned int n_points_1d, + 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) @@ -2756,13 +2584,13 @@ namespace internal for (unsigned int c = 0; c < n_components; ++c) { - if (do_hessians) + if (flag & EvaluationFlags::hessians) evalf.template apply_face(input, output); - else if (!do_hessians && do_gradients) + else if (flag & EvaluationFlags::gradients) evalf.template apply_face void - adjust_for_face_orientation(const unsigned int dim, - const unsigned int n_components, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map, - const bool integrate, - const bool values, - const bool gradients, - const bool hessians, - const unsigned int n_q_points, - Number * tmp_values, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad) + 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) { - Assert(face_orientation, ExcInternalError()); - const unsigned int *orientation = &orientation_map[face_orientation][0]; for (unsigned int c = 0; c < n_components; ++c) { - if (values == true) + if (flag & EvaluationFlags::values) { if (integrate) for (unsigned int q = 0; q < n_q_points; ++q) @@ -2935,7 +2757,7 @@ namespace internal for (unsigned int q = 0; q < n_q_points; ++q) values_quad[c * n_q_points + q] = tmp_values[q]; } - if (gradients == true) + if (flag & EvaluationFlags::gradients) for (unsigned int d = 0; d < dim; ++d) { if (integrate) @@ -2949,7 +2771,7 @@ namespace internal for (unsigned int q = 0; q < n_q_points; ++q) gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q]; } - if (hessians == true) + if (flag & EvaluationFlags::hessians) { const unsigned int hdim = (dim * (dim + 1)) / 2; for (unsigned int d = 0; d < hdim; ++d) @@ -2972,46 +2794,37 @@ namespace internal - template + template struct FEFaceEvaluationImplEvaluateSelector { template static bool - run(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - const VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map) + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + FEEvaluationBaseData + &eval) { + const auto & data = eval.get_shape_info(); + const auto & shape_data = data.data.front(); + const unsigned int face_no = eval.get_face_no(); + const unsigned int face_orientation = eval.get_face_orientation(); + if (data.element_type == MatrixFreeFunctions::tensor_none) { const unsigned int n_dofs = data.dofs_per_component_on_cell; const unsigned int n_q_points = data.n_q_points_faces[face_no]; - const auto shape_info = data.data.front(); - using Eval = EvaluatorTensorProduct; + using Eval = + EvaluatorTensorProduct; - if (evaluate_values) + if (evaluation_flag & EvaluationFlags::values) { const auto shape_values = - &shape_info.shape_values_face(face_no, face_orientation, 0); + &shape_data.shape_values_face(face_no, face_orientation, 0); - auto values_quad_ptr = values_quad; - auto values_dofs_actual_ptr = values_array; + auto values_quad_ptr = 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) @@ -3024,14 +2837,14 @@ namespace internal } } - if (evaluate_gradients) + if (evaluation_flag & EvaluationFlags::gradients) { - auto gradients_quad_ptr = gradients_quad; - auto values_dofs_actual_ptr = values_array; + auto gradients_quad_ptr = eval.begin_gradients(); + auto values_dofs_actual_ptr = values_dofs; - std::array shape_gradients; + std::array shape_gradients; for (unsigned int d = 0; d < dim; ++d) - shape_gradients[d] = &shape_info.shape_gradients_face( + shape_gradients[d] = &shape_data.shape_gradients_face( face_no, face_orientation, d, 0); for (unsigned int c = 0; c < n_components; ++c) @@ -3053,90 +2866,70 @@ namespace internal } } - if (evaluate_hessians) - { - AssertThrow(false, ExcNotImplemented()); - } + Assert(!(evaluation_flag & EvaluationFlags::hessians), + ExcNotImplemented()); return true; } - constexpr unsigned int static_dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; const unsigned int dofs_per_face = - fe_degree > -1 ? - static_dofs_per_face : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(shape_data.fe_degree + 1, dim - 1); - VectorizedArrayType *temp1 = scratch_data; + Number *temp = eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_face; - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - data, - values_array, - temp1, - evaluate_gradients, - evaluate_hessians, - face_no); + FEFaceNormalEvaluationImpl::template interpolate< + true, + false>(n_components, evaluation_flag, data, values_dofs, temp, face_no); - const unsigned int n_q_points_1d_actual = + const unsigned int subface_index = eval.get_subface_index(); + constexpr unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; + if (fe_degree > -1 && subface_index >= GeometryInfo::max_children_per_cell && data.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d_actual, - VectorizedArrayType>::evaluate_in_face(n_components, - data, - temp1, - values_quad, - gradients_quad, - hessians_quad, - scratch_data + 3 * - n_components * - dofs_per_face, - evaluate_values, - evaluate_gradients, - evaluate_hessians, - subface_index); + FEFaceEvaluationImpl::evaluate_in_face(n_components, + evaluation_flag, + shape_data, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); else - FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d_actual, - VectorizedArrayType>::evaluate_in_face(n_components, - data, - temp1, - values_quad, - gradients_quad, - hessians_quad, - scratch_data + 3 * - n_components * - dofs_per_face, - evaluate_values, - evaluate_gradients, - evaluate_hessians, - subface_index); + FEFaceEvaluationImpl::evaluate_in_face(n_components, + evaluation_flag, + shape_data, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); if (face_orientation) - adjust_for_face_orientation(dim, - n_components, - face_orientation, - orientation_map, - false, - evaluate_values, - evaluate_gradients, - evaluate_hessians, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad, - hessians_quad); + adjust_for_face_orientation( + dim, + n_components, + evaluation_flag, + &eval.get_orientation_map()(face_orientation, 0), + false, + data.n_q_points_face, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians()); return false; } @@ -3144,46 +2937,37 @@ namespace internal - template + template struct FEFaceEvaluationImplIntegrateSelector { template static bool - run(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const bool integrate_hessians, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map) + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + FEEvaluationBaseData + &eval) { + const auto & data = eval.get_shape_info(); + const auto & shape_data = data.data.front(); + const unsigned int face_no = eval.get_face_no(); + const unsigned int face_orientation = eval.get_face_orientation(); + if (data.element_type == MatrixFreeFunctions::tensor_none) { const unsigned int n_dofs = data.dofs_per_component_on_cell; const unsigned int n_q_points = data.n_q_points_faces[face_no]; - const auto shape_info = data.data.front(); - using Eval = EvaluatorTensorProduct; + using Eval = + EvaluatorTensorProduct; - if (integrate_values) + if (integration_flag & EvaluationFlags::values) { const auto shape_values = - &shape_info.shape_values_face(face_no, face_orientation, 0); + &shape_data.shape_values_face(face_no, face_orientation, 0); - auto values_quad_ptr = values_quad; - auto values_dofs_actual_ptr = values_array; + auto values_quad_ptr = 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) @@ -3196,14 +2980,14 @@ namespace internal } } - if (integrate_gradients) + if (integration_flag & EvaluationFlags::gradients) { - auto gradients_quad_ptr = gradients_quad; - auto values_dofs_actual_ptr = values_array; + auto gradients_quad_ptr = eval.begin_gradients(); + auto values_dofs_actual_ptr = values_dofs; - std::array shape_gradients; + std::array shape_gradients; for (unsigned int d = 0; d < dim; ++d) - shape_gradients[d] = &shape_info.shape_gradients_face( + shape_gradients[d] = &shape_data.shape_gradients_face( face_no, face_orientation, d, 0); for (unsigned int c = 0; c < n_components; ++c) @@ -3216,7 +3000,8 @@ namespace internal n_dofs, n_q_points); - if ((integrate_values == false) && d == 0) + if (!(integration_flag & EvaluationFlags::values) && + d == 0) eval.template gradients<0, false, false>( gradients_quad_ptr, values_dofs_actual_ptr); else @@ -3229,209 +3014,130 @@ namespace internal } } - if (integrate_hessians) - { - AssertThrow(false, ExcNotImplemented()); - } + Assert(!(integration_flag & EvaluationFlags::hessians), + ExcNotImplemented()); return true; } - if (face_orientation) - adjust_for_face_orientation(dim, - n_components, - face_orientation, - orientation_map, - true, - integrate_values, - integrate_gradients, - integrate_hessians, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad, - hessians_quad); - - constexpr unsigned int static_dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; const unsigned int dofs_per_face = - fe_degree > -1 ? - static_dofs_per_face : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(shape_data.fe_degree + 1, dim - 1); + + Number *temp = eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_face; - VectorizedArrayType *temp1 = scratch_data; + if (face_orientation) + adjust_for_face_orientation( + dim, + n_components, + integration_flag, + &eval.get_orientation_map()(face_orientation, 0), + true, + data.n_q_points_face, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians()); const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; + const unsigned int subface_index = eval.get_subface_index(); + if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && + eval.get_subface_index() >= + GeometryInfo::max_children_per_cell && data.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d_actual, - VectorizedArrayType>::integrate_in_face(n_components, - data, - temp1, - values_quad, - gradients_quad, - hessians_quad, - scratch_data + - 3 * n_components * - dofs_per_face, - integrate_values, - integrate_gradients, - integrate_hessians, - subface_index); + FEFaceEvaluationImpl::integrate_in_face(n_components, + integration_flag, + shape_data, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); else - FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d_actual, - VectorizedArrayType>::integrate_in_face(n_components, - data, - temp1, - values_quad, - gradients_quad, - hessians_quad, - scratch_data + - 3 * n_components * - dofs_per_face, - integrate_values, - integrate_gradients, - integrate_hessians, - subface_index); + FEFaceEvaluationImpl::integrate_in_face(n_components, + integration_flag, + shape_data, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - data, - temp1, - values_array, - integrate_gradients, - integrate_hessians, - face_no); + FEFaceNormalEvaluationImpl:: + template interpolate( + n_components, integration_flag, data, temp, values_dofs, face_no); return false; } }; - template - static bool - fe_face_evaluation_process_and_io(Processor &proc) + template + void + fe_face_evaluation_process_and_io( + const Processor & proc, + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + typename Processor::Number2_ * src_ptr, + const std::vector> *sm_ptr, + const EvaluationData & eval, + typename Processor::VectorizedArrayType_ * temp1) { - auto n_components = proc.n_components; - auto integrate = proc.integrate; - auto global_vector_ptr = proc.global_vector_ptr; - auto &sm_ptr = proc.sm_ptr; - auto &data = proc.data; - auto &dof_info = proc.dof_info; - auto values_quad = proc.values_quad; - auto gradients_quad = proc.gradients_quad; - auto hessians_quad = proc.hessians_quad; - auto scratch_data = proc.scratch_data; - auto do_values = proc.do_values; - auto do_gradients = proc.do_gradients; - auto do_hessians = proc.do_hessians; - auto active_fe_index = proc.active_fe_index; - auto first_selected_component = proc.first_selected_component; - auto cells = proc.cells; - auto face_nos = proc.face_nos; - auto subface_index = proc.subface_index; - auto dof_access_index = proc.dof_access_index; - auto face_orientations = proc.face_orientations; - auto &orientation_map = proc.orientation_map; - - static const int dim = Processor::dim_; - static const int fe_degree = Processor::fe_degree_; - using VectorizedArrayType = typename Processor::VectorizedArrayType_; + (void)proc; + (void)n_components; + (void)evaluation_flag; + (void)src_ptr; + (void)sm_ptr; + (void)eval; + (void)temp1; + + /* + constexpr int dim = Processor::dim_; + constexpr int fe_degree = Processor::fe_degree_; + using VectorizedArrayType = typename Processor::VectorizedArrayType_; using Number = typename Processor::Number_; using Number2_ = typename Processor::Number2_; - const unsigned int cell = cells[0]; - - // In the case of integration, we do not need to reshuffle the - // data at the quadrature points to adjust for the face - // orientation if the shape functions are nodal at the cell - // boundaries (and we only requested the integration of the - // values) or Hermite shape functions are used. These cases are - // handled later when the values are written back into the - // glrobal vector. - if (integrate && - (face_orientations[0] > 0 && - (subface_index < GeometryInfo::max_children_per_cell || - !(((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true && - fe_degree > 0) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) && - (dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_strided || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_mixed_strides || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - contiguous))))) - { - AssertDimension(n_face_orientations, 1); - adjust_for_face_orientation(dim, - n_components, - face_orientations[0], - orientation_map, - true, - do_values, - do_gradients, - do_hessians, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad, - hessians_quad); - } + const auto & shape_data = eval.get_shape_info().data[0]; + const auto element_type = eval.get_shape_info().element_type; + constexpr bool integrate = Processor::do_integrate; + const unsigned int face_no = eval.get_face_no(); // 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 = - (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 && - data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ? - data.data.front() - .shape_data_on_face[0][fe_degree + (integrate ? - (2 - (face_nos[0] % 2)) : - (1 + (face_nos[0] % 2)))] : - VectorizedArrayType(0.0 /*dummy*/); + shape_data + .shape_data_on_face[0][fe_degree + integrate ? (2 - face_no % 2) : + (1 + face_no % 2)]; constexpr unsigned int static_dofs_per_component = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) : - numbers::invalid_unsigned_int; - constexpr unsigned int static_dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; - const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); - - VectorizedArrayType *temp1 = scratch_data; - - const unsigned int dummy = 0; + Utilities::pow(fe_degree + 1, dim); + constexpr unsigned int dofs_per_face = + Utilities::pow(fe_degree + 1, dim - 1); // re-orientation std::array orientation = {}; if (n_face_orientations == 1) - orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ? - &data.face_orientations[face_orientations[0]][0] : - &dummy; + orientation[0] = + &eval.get_orientation_map()(eval.get_face_orientation(), 0); else { for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) @@ -3443,7 +3149,7 @@ namespace internal continue; orientation[v] = - (data.data.front().nodal_at_cell_boundaries == true) ? + (shape_data.nodal_at_cell_boundaries == true) ? &data.face_orientations[face_orientations[v]][0] : &dummy; } @@ -3517,541 +3223,471 @@ namespace internal for (unsigned int comp = 0; comp < n_components; ++comp) { - if (integrate) - proc.in_face_operation(temp1, comp); - // we can only use the fast functions if we know the polynomial degree // as a template parameter (fe_degree != -1), and it only makes sense // to use the functions for at least linear functions for values on // the faces and quadratic functions for gradients on the faces, so // include the switch here - if (((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true && - fe_degree > 0) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) && - do_hessians == false) + // 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) { - // 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(n_face_orientations, 1); + + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::size()); + Number2_ *vector_ptr = + global_vector_ptr + + dof_info + .dof_indices_contiguous[dof_access_index] + [cell * VectorizedArrayType::size()] + + (dof_info.component_dof_indices_offset[active_fe_index] + [first_selected_component] + + comp * static_dofs_per_component) * + VectorizedArrayType::size(); + + if (fe_degree > 1 && do_gradients == true) { - AssertDimension(n_face_orientations, 1); - - AssertDimension( - dof_info.n_vectorization_lanes_filled[dof_access_index][cell], - VectorizedArrayType::size()); - Number2_ *vector_ptr = - global_vector_ptr + - dof_info.dof_indices_contiguous[dof_access_index] - [cell * - VectorizedArrayType::size()] + - (dof_info - .component_dof_indices_offset[active_fe_index] - [first_selected_component] + - comp * static_dofs_per_component) * - VectorizedArrayType::size(); - - if (fe_degree > 1 && do_gradients == true) + for (unsigned int i = 0; i < dofs_per_face; ++i) { - for (unsigned int i = 0; i < dofs_per_face; ++i) + if (n_face_orientations == 1) { - if (n_face_orientations == 1) - { - const unsigned int ind1 = - index_array_hermite[0][2 * i]; - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1]; - AssertIndexRange(ind1, - data.dofs_per_component_on_cell); - AssertIndexRange(ind2, - data.dofs_per_component_on_cell); - const unsigned int i_ = reorientate(0, i); - proc.hermite_grad_vectorized( - temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1 * VectorizedArrayType::size(), - vector_ptr + ind2 * VectorizedArrayType::size(), - grad_weight); - } - else - { - Assert(false, ExcNotImplemented()); - } + const unsigned int ind1 = index_array_hermite[0][2 * i]; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1]; + AssertIndexRange(ind1, data.dofs_per_component_on_cell); + AssertIndexRange(ind2, data.dofs_per_component_on_cell); + const unsigned int i_ = reorientate(0, i); + proc.hermite_grad_vectorized( + temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1 * VectorizedArrayType::size(), + vector_ptr + ind2 * VectorizedArrayType::size(), + grad_weight); } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = index_array_nodal[0][i]; - proc.value_vectorized( - temp1[i_], - vector_ptr + ind * VectorizedArrayType::size()); - } - else - { - Assert(false, ExcNotImplemented()); - } + Assert(false, ExcNotImplemented()); } } } - - // 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) + else { - AssertDimension(n_face_orientations, 1); - - AssertDimension( - dof_info.n_vectorization_lanes_filled[dof_access_index][cell], - VectorizedArrayType::size()); - const unsigned int *indices = - &dof_info.dof_indices_contiguous[dof_access_index] - [cell * - VectorizedArrayType::size()]; - Number2_ *vector_ptr = - global_vector_ptr + - (comp * static_dofs_per_component + - dof_info - .component_dof_indices_offset[active_fe_index] - [first_selected_component]) * - VectorizedArrayType::size(); - if (fe_degree > 1 && do_gradients == true) + for (unsigned int i = 0; i < dofs_per_face; ++i) { - for (unsigned int i = 0; i < dofs_per_face; ++i) + if (n_face_orientations == 1) { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - const unsigned int ind1 = - index_array_hermite[0][2 * i] * - VectorizedArrayType::size(); - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1] * - VectorizedArrayType::size(); - proc.hermite_grad_vectorized_indexed( - temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1, - vector_ptr + ind2, - grad_weight, - indices, - indices); - } - else - { - Assert(false, 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 * + VectorizedArrayType::size()); } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = - index_array_nodal[0][i] * - VectorizedArrayType::size(); - proc.value_vectorized_indexed(temp1[i_], - vector_ptr + ind, - indices); - } - else - { - Assert(false, ExcNotImplemented()); - } + Assert(false, ExcNotImplemented()); } } } + } - // 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) + // 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(n_face_orientations, 1); + + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::size()); + const unsigned int *indices = + &dof_info + .dof_indices_contiguous[dof_access_index] + [cell * VectorizedArrayType::size()]; + Number2_ *vector_ptr = + global_vector_ptr + + (comp * static_dofs_per_component + + dof_info + .component_dof_indices_offset[active_fe_index] + [first_selected_component]) * + VectorizedArrayType::size(); + if (fe_degree > 1 && do_gradients == true) { - AssertDimension(n_face_orientations, 1); - - const unsigned int *strides = - &dof_info.dof_indices_interleave_strides - [dof_access_index][cell * VectorizedArrayType::size()]; - unsigned int indices[VectorizedArrayType::size()]; - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - indices[v] = - dof_info.dof_indices_contiguous - [dof_access_index] - [cell * VectorizedArrayType::size() + v] + - (dof_info - .component_dof_indices_offset[active_fe_index] - [first_selected_component] + - comp * static_dofs_per_component) * - strides[v]; - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - if (fe_degree > 1 && do_gradients == true) + for (unsigned int i = 0; i < dofs_per_face; ++i) { - if (n_filled_lanes == VectorizedArrayType::size()) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - unsigned int ind1[VectorizedArrayType::size()]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++v) - ind1[v] = - indices[v] + - index_array_hermite[0 /*TODO*/][2 * i] * - strides[v]; - unsigned int ind2[VectorizedArrayType::size()]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++v) - ind2[v] = - indices[v] + - index_array_hermite[0 /*TODO*/][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 - { - Assert(false, ExcNotImplemented()); - } - } - else + if (n_face_orientations == 1) { - 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]); - } + const unsigned int i_ = reorientate(0, i); + const unsigned int ind1 = + index_array_hermite[0][2 * i] * + VectorizedArrayType::size(); + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1] * + VectorizedArrayType::size(); + proc.hermite_grad_vectorized_indexed( + temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1, + vector_ptr + ind2, + grad_weight, + indices, + indices); } - } - else - { - if (n_filled_lanes == VectorizedArrayType::size()) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - if (n_face_orientations == 1) - { - unsigned int ind[VectorizedArrayType::size()]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++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 - { - Assert(false, ExcNotImplemented()); - } - } 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]]); + Assert(false, ExcNotImplemented()); } } } - - // 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)) + else { - const unsigned int *indices = - &dof_info.dof_indices_contiguous[dof_access_index] - [cell * - VectorizedArrayType::size()]; - Number2_ *vector_ptr = - global_vector_ptr + comp * static_dofs_per_component + - dof_info - .component_dof_indices_offset[active_fe_index] - [first_selected_component]; - - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - const bool vectorization_possible = - (n_face_orientations == 1) && - (n_filled_lanes == VectorizedArrayType::size()) && - (sm_ptr != nullptr); - - std::array - vector_ptrs = {}; - - if (vectorization_possible == false) + for (unsigned int i = 0; i < dofs_per_face; ++i) { if (n_face_orientations == 1) { - for (unsigned int v = 0; v < n_filled_lanes; ++v) - if (sm_ptr == nullptr) - { - vector_ptrs[v] = vector_ptr + indices[v]; - } - else - { - const auto &temp = - dof_info.dof_indices_contiguous_sm - [dof_access_index] - [cell * VectorizedArrayType::size() + v]; - vector_ptrs[v] = const_cast( - sm_ptr->operator[](temp.first).data() + - temp.second + comp * static_dofs_per_component + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component]); - } - } - else if (n_face_orientations == VectorizedArrayType::size()) - { - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++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 + - comp * static_dofs_per_component + - dof_info.component_dof_indices_offset - [active_fe_index] - [first_selected_component]); - } - } + const unsigned int i_ = reorientate(0, i); + const unsigned int ind = + index_array_nodal[0][i] * VectorizedArrayType::size(); + proc.value_vectorized_indexed(temp1[i_], + vector_ptr + ind, + indices); } else { Assert(false, ExcNotImplemented()); } } + } + } - if (do_gradients == true && - data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite) - { - if (vectorization_possible) - for (unsigned int i = 0; i < dofs_per_face; ++i) + // 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) + { + AssertDimension(n_face_orientations, 1); + + const unsigned int *strides = + &dof_info.dof_indices_interleave_strides + [dof_access_index][cell * VectorizedArrayType::size()]; + unsigned int indices[VectorizedArrayType::size()]; + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + indices[v] = + dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::size() + v] + + (dof_info + .component_dof_indices_offset[active_fe_index] + [first_selected_component] + + comp * static_dofs_per_component) * + strides[v]; + const unsigned int n_filled_lanes = + dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; + + if (fe_degree > 1 && do_gradients == true) + { + if (n_filled_lanes == VectorizedArrayType::size()) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + if (n_face_orientations == 1) { - 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); - + unsigned int ind1[VectorizedArrayType::size()]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++v) + ind1[v] = + indices[v] + + // TODO + index_array_hermite[0][2 * i] * strides[v]; + unsigned int ind2[VectorizedArrayType::size()]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++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], - vector_ptr + ind1, - vector_ptr + ind2, + global_vector_ptr, + global_vector_ptr, grad_weight, - indices, - indices); + ind1, + ind2); } - else if (n_face_orientations == 1) - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - 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 < VectorizedArrayType::size(); - ++v) - { - temp1[i_][v] = 0.0; - temp1[i_ + dofs_per_face][v] = 0.0; - } + Assert(false, ExcNotImplemented()); } - else - { - 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) + 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 ind = index_array_nodal[0][i]; - const unsigned int i_ = reorientate(0, 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 == VectorizedArrayType::size()) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + if (n_face_orientations == 1) + { + unsigned int ind[VectorizedArrayType::size()]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++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_], - vector_ptr + ind, - indices); + global_vector_ptr, + ind); } - else if (n_face_orientations == 1) - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - const unsigned int ind = index_array_nodal[0][i]; - const unsigned int i_ = reorientate(0, i); + Assert(false, ExcNotImplemented()); + } + } + 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) - proc.value(temp1[i_][v], vector_ptrs[v][ind]); + 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]]); + } + } + } - if (integrate == false) - for (unsigned int v = n_filled_lanes; - v < VectorizedArrayType::size(); - ++v) - temp1[i_][v] = 0.0; + // 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)) + { + const unsigned int *indices = + &dof_info + .dof_indices_contiguous[dof_access_index] + [cell * VectorizedArrayType::size()]; + Number2_ *vector_ptr = + global_vector_ptr + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset[active_fe_index] + [first_selected_component]; + + const unsigned int n_filled_lanes = + dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; + + const bool vectorization_possible = + (n_face_orientations == 1) && + (n_filled_lanes == VectorizedArrayType::size()) && + (sm_ptr != nullptr); + + std::array vector_ptrs = + {}; + + if (vectorization_possible == false) + { + if (n_face_orientations == 1) + { + for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (sm_ptr == nullptr) + { + vector_ptrs[v] = vector_ptr + indices[v]; } - else - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++v) - if (cells[v] != numbers::invalid_unsigned_int) - proc.value( - temp1[reorientate(v, i)][v], - vector_ptrs[v][index_array_nodal[v][i]]); - else if (integrate == false) - temp1[i][v] = 0.0; + const auto &temp = + dof_info.dof_indices_contiguous_sm + [dof_access_index] + [cell * VectorizedArrayType::size() + v]; + vector_ptrs[v] = const_cast( + sm_ptr->operator[](temp.first).data() + + temp.second + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component]); + } + } + else if (n_face_orientations == VectorizedArrayType::size()) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); + ++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 + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component]); + } } } + else + { + Assert(false, ExcNotImplemented()); + } + } + + if (do_gradients == true && + data.element_type == + MatrixFreeFunctions::tensor_symmetric_hermite) + { + 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, + indices, + 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 < VectorizedArrayType::size(); + ++v) + { + temp1[i_][v] = 0.0; + temp1[i_ + dofs_per_face][v] = 0.0; + } + } + else + { + 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 { - // case 5: default vector access - // for the integrate_scatter path (integrate == true), we - // need to only prepare the data in this function for all - // components to later call distribute_local_to_global(); - // for the gather_evaluate path (integrate == false), we - // instead want to leave early because we need to get the - // vector data from somewhere else - proc.default_operation(temp1, comp); - if (integrate) - accesses_global_vector = false; + 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, + 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 < VectorizedArrayType::size(); + ++v) + temp1[i_][v] = 0.0; + } else - return false; + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); + ++v) + if (cells[v] != numbers::invalid_unsigned_int) + proc.value(temp1[reorientate(v, i)][v], + vector_ptrs[v][index_array_nodal[v][i]]); + } } } else { - // case 5: default vector access - proc.default_operation(temp1, comp); - if (integrate) - accesses_global_vector = false; - else - return false; + // We should not end up here, this should be caught by + // FEFaceEvaluationImplGatherEvaluateSelector::supports() + Assert(false, ExcInternalError()); } - - if (!integrate) - proc.in_face_operation(temp1, comp); } - - if (!integrate && - (face_orientations[0] > 0 && - subface_index < GeometryInfo::max_children_per_cell)) - { - AssertDimension(n_face_orientations, 1); - adjust_for_face_orientation(dim, - n_components, - face_orientations[0], - orientation_map, - false, - do_values, - do_gradients, - do_hessians, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad, - hessians_quad); - } - - return accesses_global_vector; + */ } @@ -4065,121 +3701,118 @@ namespace internal template static bool run(const unsigned int n_components, - const unsigned int n_face_orientations, + const EvaluationFlags::EvaluationFlags evaluation_flag, const Number2 * src_ptr, const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map) + FEEvaluationBaseData &eval) { - if (src_ptr == nullptr) - return false; + Assert(fe_degree > -1, ExcInternalError()); + Assert(eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric, + ExcInternalError()); - if (data.element_type == MatrixFreeFunctions::tensor_none) - return false; + const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1); + + VectorizedArrayType *temp = eval.get_scratch_data().begin(); + VectorizedArrayType *scratch_data = + temp + 3 * n_components * dofs_per_face; + + Processor p; + + if (eval.get_dof_access_index() == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + eval.get_is_interior_face() == false) + fe_face_evaluation_process_and_io( + p, n_components, evaluation_flag, src_ptr, sm_ptr, eval, temp); + else + fe_face_evaluation_process_and_io<1>( + p, n_components, evaluation_flag, src_ptr, sm_ptr, eval, temp); + + const unsigned int subface_index = eval.get_subface_index(); + + if (subface_index >= GeometryInfo::max_children_per_cell) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d, + VectorizedArrayType>::evaluate_in_face(n_components, + evaluation_flag, + eval.get_shape_info().data[0], + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); + else + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d, + VectorizedArrayType>::evaluate_in_face(n_components, + evaluation_flag, + eval.get_shape_info().data[0], + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); + + if (eval.get_face_orientation()) + adjust_for_face_orientation( + dim, + n_components, + evaluation_flag, + &eval.get_orientation_map()(eval.get_face_orientation(), 0), + false, + eval.get_shape_info().n_q_points_face, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians()); + + return false; + } - (void)sm_ptr; - - Processor p(n_components, - false, - src_ptr, - sm_ptr, - data, - dof_info, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - evaluate_values, - evaluate_gradients, - evaluate_hessians, - active_fe_index, - first_selected_component, - cells, - face_nos, - subface_index, - dof_access_index, - face_orientations, - orientation_map); - - if (n_face_orientations == VectorizedArrayType::size()) - return fe_face_evaluation_process_and_io( - p); + static bool + supports( + const EvaluationFlags::EvaluationFlags evaluation_flag, + const internal::MatrixFreeFunctions::ShapeInfo + & shape_info, + const Number * vector_ptr, + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage) + { + const unsigned int fe_degree = shape_info.data[0].fe_degree; + if (fe_degree < 1 || !shape_info.data[0].nodal_at_cell_boundaries || + (evaluation_flag & EvaluationFlags::gradients && + (fe_degree < 2 || + shape_info.data[0].element_type != + MatrixFreeFunctions::tensor_symmetric_hermite)) || + (evaluation_flag & EvaluationFlags::hessians) || + vector_ptr == nullptr || + shape_info.data[0].element_type > + MatrixFreeFunctions::tensor_symmetric || + storage < + MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) + return false; else - return fe_face_evaluation_process_and_io<1>(p); + return false; } private: template struct Processor { - static const int dim_ = dim; - static const int fe_degree_ = fe_degree; - static const int n_q_points_1d_ = n_q_points_1d; - using VectorizedArrayType_ = VectorizedArrayType; - using Number_ = Number; - using Number2_ = const Number2; - - Processor( - const unsigned int n_components, - const bool integrate, - const Number2 * global_vector_ptr, - const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data, - const bool do_values, - const bool do_gradients, - const bool do_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map) - : n_components(n_components) - , integrate(integrate) - , global_vector_ptr(global_vector_ptr) - , sm_ptr(sm_ptr) - , data(data) - , dof_info(dof_info) - , values_quad(values_quad) - , gradients_quad(gradients_quad) - , hessians_quad(hessians_quad) - , scratch_data(scratch_data) - , do_values(do_values) - , do_gradients(do_gradients) - , do_hessians(do_hessians) - , active_fe_index(active_fe_index) - , first_selected_component(first_selected_component) - , cells(cells) - , face_nos(face_nos) - , subface_index(subface_index) - , dof_access_index(dof_access_index) - , face_orientations(face_orientations) - , orientation_map(orientation_map) - {} + static const bool do_integrate = false; + static const int dim_ = dim; + static const int fe_degree_ = fe_degree; + static const int n_q_points_1d_ = n_q_points_1d; + using VectorizedArrayType_ = VectorizedArrayType; + using Number_ = Number; + using Number2_ = const Number2; template void @@ -4243,86 +3876,6 @@ namespace internal // case 3b) temp = src_ptr; } - - template - void - default_operation(const T1 &, const unsigned int) - { - // case 5) - } - - template - void - in_face_operation(T1 &temp1, const unsigned int comp) - { - const unsigned int dofs_per_face = - fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); - const unsigned int n_q_points = - fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : - data.n_q_points_face; - constexpr unsigned int hdim = (dim * (dim + 1)) / 2; - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - data.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl:: - evaluate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + comp * dim * n_q_points, - hessians_quad + comp * hdim * n_q_points, - scratch_data + 3 * dofs_per_face, - do_values, - do_gradients, - do_hessians, - subface_index); - else - FEFaceEvaluationImpl:: - evaluate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + comp * dim * n_q_points, - hessians_quad + comp * hdim * n_q_points, - scratch_data + 3 * dofs_per_face, - do_values, - do_gradients, - do_hessians, - subface_index); - } - - const unsigned int n_components; - const bool integrate; - const Number2 * global_vector_ptr; - const std::vector> *sm_ptr; - const MatrixFreeFunctions::ShapeInfo &data; - const MatrixFreeFunctions::DoFInfo & dof_info; - VectorizedArrayType * values_quad; - VectorizedArrayType * gradients_quad; - VectorizedArrayType * hessians_quad; - VectorizedArrayType * scratch_data; - const bool do_values; - const bool do_gradients; - const bool do_hessians; - const unsigned int active_fe_index; - const unsigned int first_selected_component; - const std::array cells; - const std::array face_nos; - const unsigned int subface_index; - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index; - const std::array - face_orientations; - const Table<2, unsigned int> &orientation_map; }; }; @@ -4335,147 +3888,94 @@ namespace internal template static bool run(const unsigned int n_components, - const unsigned int n_face_orientations, + const EvaluationFlags::EvaluationFlags integration_flag, Number2 * dst_ptr, const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const bool integrate_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map) + FEEvaluationBaseData &eval) { - (void)sm_ptr; + Assert(fe_degree > -1, ExcInternalError()); + Assert(eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric, + ExcInternalError()); - if (dst_ptr == nullptr || - data.element_type == MatrixFreeFunctions::tensor_none) - { - AssertDimension(n_face_orientations, 1); - - // for block vectors simply integrate - FEFaceEvaluationImplIntegrateSelector:: - template run(n_components, - data, - values_array, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - integrate_values, - integrate_gradients, - integrate_hessians, - face_nos[0], - subface_index, - face_orientations[0], - orientation_map); - - // default vector access - return false; - } + const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1); + + VectorizedArrayType *temp = eval.get_scratch_data().begin(); + VectorizedArrayType *scratch_data = + temp + 3 * n_components * dofs_per_face; + + if (eval.get_face_orientation()) + adjust_for_face_orientation( + dim, + n_components, + integration_flag, + &eval.get_orientation_map()(eval.get_face_orientation(), 0), + true, + eval.get_shape_info().n_q_points_face, + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians()); + + const unsigned int subface_index = eval.get_subface_index(); + + if (fe_degree > -1 && eval.get_subface_index() >= + GeometryInfo::max_children_per_cell) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d, + VectorizedArrayType>::integrate_in_face(n_components, + integration_flag, + eval.get_shape_info().data[0], + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); + else + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d, + VectorizedArrayType>::integrate_in_face(n_components, + integration_flag, + eval.get_shape_info().data[0], + temp, + eval.begin_values(), + eval.begin_gradients(), + eval.begin_hessians(), + scratch_data, + subface_index); + Processor p; - Processor p(values_array, - n_components, - true, - dst_ptr, - sm_ptr, - data, - dof_info, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - integrate_values, - integrate_gradients, - integrate_hessians, - active_fe_index, - first_selected_component, - cells, - face_nos, - subface_index, - dof_access_index, - face_orientations, - orientation_map); - - if (n_face_orientations == VectorizedArrayType::size()) - return fe_face_evaluation_process_and_io( - p); + if (eval.get_dof_access_index() == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + eval.get_is_interior_face() == false) + fe_face_evaluation_process_and_io( + p, n_components, integration_flag, dst_ptr, sm_ptr, eval, temp); else - return fe_face_evaluation_process_and_io<1>(p); + fe_face_evaluation_process_and_io<1>( + p, n_components, integration_flag, dst_ptr, sm_ptr, eval, temp); + + return false; } private: template struct Processor { - static const int dim_ = dim; - static const int fe_degree_ = fe_degree; - static const int n_q_points_1d_ = n_q_points_1d; - using VectorizedArrayType_ = VectorizedArrayType; - using Number_ = Number; - using Number2_ = Number2; - - - Processor( - VectorizedArrayType * values_array, - const unsigned int n_components, - const bool integrate, - Number2 * global_vector_ptr, - const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data, - const bool do_values, - const bool do_gradients, - const bool do_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map) - : values_array(values_array) - , n_components(n_components) - , integrate(integrate) - , global_vector_ptr(global_vector_ptr) - , sm_ptr(sm_ptr) - , data(data) - , dof_info(dof_info) - , values_quad(values_quad) - , gradients_quad(gradients_quad) - , hessians_quad(hessians_quad) - , scratch_data(scratch_data) - , do_values(do_values) - , do_gradients(do_gradients) - , do_hessians(do_hessians) - , active_fe_index(active_fe_index) - , first_selected_component(first_selected_component) - , cells(cells) - , face_nos(face_nos) - , subface_index(subface_index) - , dof_access_index(dof_access_index) - , face_orientations(face_orientations) - , orientation_map(orientation_map) - {} + static const bool do_integrate = true; + static const int dim_ = dim; + static const int fe_degree_ = fe_degree; + static const int n_q_points_1d_ = n_q_points_1d; + using VectorizedArrayType_ = VectorizedArrayType; + using Number_ = Number; + using Number2_ = Number2; template void @@ -4547,101 +4047,6 @@ namespace internal // case 3b) dst_ptr += temp; } - - template - void - default_operation(const T0 &temp1, const unsigned int comp) - { - // case 5: default vector access, must be handled separately, just do - // the face-normal interpolation - - FEFaceNormalEvaluationImpl:: - template interpolate( - /* n_components */ 1, - data, - temp1, - values_array + comp * data.dofs_per_component_on_cell, - do_gradients, - do_hessians, - face_nos[0]); - } - - template - void - in_face_operation(T0 &temp1, const unsigned int comp) - { - const unsigned int dofs_per_face = - fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); - const unsigned int n_q_points = - fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : - data.n_q_points_face; - constexpr unsigned int hdim = (dim * (dim + 1)) / 2; - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - data.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl:: - integrate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + dim * comp * n_q_points, - hessians_quad + hdim * comp * n_q_points, - scratch_data + 3 * dofs_per_face, - do_values, - do_gradients, - do_hessians, - subface_index); - else - internal::FEFaceEvaluationImpl:: - integrate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + dim * comp * n_q_points, - hessians_quad + hdim * comp * n_q_points, - scratch_data + 3 * dofs_per_face, - do_values, - do_gradients, - do_hessians, - subface_index); - } - - VectorizedArrayType *values_array; - - - const unsigned int n_components; - const bool integrate; - Number2 * global_vector_ptr; - const std::vector> *sm_ptr; - const MatrixFreeFunctions::ShapeInfo &data; - const MatrixFreeFunctions::DoFInfo & dof_info; - VectorizedArrayType * values_quad; - VectorizedArrayType * gradients_quad; - VectorizedArrayType * hessians_quad; - VectorizedArrayType * scratch_data; - const bool do_values; - const bool do_gradients; - const bool do_hessians; - const unsigned int active_fe_index; - const unsigned int first_selected_component; - const std::array cells; - const std::array face_nos; - const unsigned int subface_index; - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index; - const std::array - face_orientations; - const Table<2, unsigned int> &orientation_map; }; }; diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h index 1c9d8b61ef..e461015786 100644 --- a/include/deal.II/matrix_free/evaluation_selector.h +++ b/include/deal.II/matrix_free/evaluation_selector.h @@ -35,6 +35,9 @@ DEAL_II_NAMESPACE_OPEN template struct SelectEvaluator { + using ScalarNumber = + typename internal::VectorizedArrayTrait::value_type; + /** * Chooses an appropriate evaluation strategy for the evaluate function, i.e. * this calls internal::FEEvaluationImpl::evaluate(), @@ -45,12 +48,8 @@ struct SelectEvaluator static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const internal::MatrixFreeFunctions::ShapeInfo &shape_info, - Number *values_dofs_actual, - Number *values_quad, - Number *gradients_quad, - Number *hessians_quad, - Number *scratch_data); + const Number * values_dofs, + FEEvaluationBaseData &eval); /** * Chooses an appropriate evaluation strategy for the integrate function, i.e. @@ -62,12 +61,8 @@ struct SelectEvaluator static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - const internal::MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, + Number * values_dofs, + FEEvaluationBaseData &eval, const bool sum_into_values_array = false); }; @@ -79,24 +74,14 @@ inline void SelectEvaluator::evaluate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const internal::MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data) + const Number * values_dofs, + FEEvaluationBaseData &eval) { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); - internal::FEEvaluationImplEvaluateSelector:: - template run(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + internal::FEEvaluationImplEvaluateSelector::template run< + fe_degree, + n_q_points_1d>(n_components, evaluation_flag, values_dofs, eval); } @@ -106,26 +91,15 @@ inline void SelectEvaluator::integrate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - const internal::MatrixFreeFunctions::ShapeInfo &shape_info, - Number * values_dofs_actual, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, + Number * values_dofs, + FEEvaluationBaseData &eval, const bool sum_into_values_array) { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); internal::FEEvaluationImplIntegrateSelector:: - template run(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - sum_into_values_array); + template run( + n_components, integration_flag, values_dofs, eval, sum_into_values_array); } #endif // DOXYGEN diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index 1dbeef01a3..e0152471f7 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -35,33 +35,23 @@ class FEEvaluationBaseData; namespace internal { - template > + template struct FEEvaluationFactory { static void evaluate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArrayType *values_dofs_actual, - VectorizedArrayType *values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data); + const VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval); static void integrate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArrayType *values_dofs_actual, - VectorizedArrayType *values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data, - const bool sum_into_values_array); + VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval, + const bool sum_into_values_array); static bool fast_evaluation_supported(const unsigned int given_degree, @@ -76,87 +66,34 @@ namespace internal struct FEFaceEvaluationFactory { static void - evaluate(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - const VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map); + evaluate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval); static void - integrate(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const bool integrate_hessians, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map); + integrate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval); - static bool + static void gather_evaluate( const unsigned int n_components, - const std::size_t n_face_orientations, + const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * src_ptr, const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map); + FEEvaluationBaseData &eval); - static bool + static void integrate_scatter( const unsigned int n_components, - const std::size_t n_face_orientations, + const EvaluationFlags::EvaluationFlags integration_flag, Number * dst_ptr, const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const bool integrate_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map); + FEEvaluationBaseData &eval); static bool fast_evaluation_supported(const unsigned int given_degree, diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index 3f6272b615..8dabe9598f 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -44,28 +44,20 @@ namespace internal template void FEEvaluationFactory::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArrayType *values_dofs_actual, - VectorizedArrayType *values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval) { instantiation_helper_run< 1, FEEvaluationImplEvaluateSelector>( - shape_info.data[0].fe_degree, - shape_info.data[0].n_q_points_1d, + eval.get_shape_info().data[0].fe_degree, + eval.get_shape_info().data[0].n_q_points_1d, n_components, evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); + values_dofs, + eval); } @@ -73,29 +65,21 @@ namespace internal template void FEEvaluationFactory::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArrayType *values_dofs_actual, - VectorizedArrayType *values_quad, - VectorizedArrayType *gradients_quad, - VectorizedArrayType *hessians_quad, - VectorizedArrayType *scratch_data, - const bool sum_into_values_array) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval, + const bool sum_into_values_array) { instantiation_helper_run< 1, FEEvaluationImplIntegrateSelector>( - shape_info.data[0].fe_degree, - shape_info.data[0].n_q_points_1d, + eval.get_shape_info().data[0].fe_degree, + eval.get_shape_info().data[0].n_q_points_1d, n_components, integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, + values_dofs, + eval, sum_into_values_array); } @@ -116,40 +100,20 @@ namespace internal template void FEFaceEvaluationFactory::evaluate( - const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - const VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval) { instantiation_helper_run< 1, FEFaceEvaluationImplEvaluateSelector>( - data.data[0].fe_degree, - data.data[0].n_q_points_1d, + eval.get_shape_info().data[0].fe_degree, + eval.get_shape_info().data[0].n_q_points_1d, n_components, - data, - values_array, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - evaluate_values, - evaluate_gradients, - evaluate_hessians, - face_no, - subface_index, - face_orientation, - orientation_map); + evaluation_flag, + values_dofs, + eval); } @@ -157,158 +121,70 @@ namespace internal template void FEFaceEvaluationFactory::integrate( - const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const bool integrate_hessians, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + VectorizedArrayType * values_dofs, + FEEvaluationBaseData &eval) { instantiation_helper_run< 1, FEFaceEvaluationImplIntegrateSelector>( - data.data[0].fe_degree, - data.data[0].n_q_points_1d, + eval.get_shape_info().data[0].fe_degree, + eval.get_shape_info().data[0].n_q_points_1d, n_components, - data, - values_array, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - integrate_values, - integrate_gradients, - integrate_hessians, - face_no, - subface_index, - face_orientation, - orientation_map); + integration_flag, + values_dofs, + eval); } template - bool + void FEFaceEvaluationFactory::gather_evaluate( const unsigned int n_components, - const std::size_t n_face_orientations, + const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * src_ptr, const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map) + FEEvaluationBaseData &eval) { - return instantiation_helper_run< + instantiation_helper_run< 1, FEFaceEvaluationImplGatherEvaluateSelector>( - data.data[0].fe_degree, - data.data[0].n_q_points_1d, + eval.get_shape_info().data[0].fe_degree, + eval.get_shape_info().data[0].n_q_points_1d, n_components, - n_face_orientations, + evaluation_flag, src_ptr, sm_ptr, - data, - dof_info, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - evaluate_values, - evaluate_gradients, - evaluate_hessians, - active_fe_index, - first_selected_component, - cells, - face_nos, - subface_index, - dof_access_index, - face_orientations, - orientation_map); + eval); } template - bool + void FEFaceEvaluationFactory::integrate_scatter( const unsigned int n_components, - const std::size_t n_face_orientations, + const EvaluationFlags::EvaluationFlags integration_flag, Number * dst_ptr, const std::vector> *sm_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const bool integrate_hessians, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array - face_orientations, - const Table<2, unsigned int> &orientation_map) + FEEvaluationBaseData &eval) { - return instantiation_helper_run< + instantiation_helper_run< 1, FEFaceEvaluationImplIntegrateScatterSelector>( - data.data[0].fe_degree, - data.data[0].n_q_points_1d, + eval.get_shape_info().data[0].fe_degree, + eval.get_shape_info().data[0].n_q_points_1d, n_components, - n_face_orientations, + integration_flag, dst_ptr, sm_ptr, - data, - dof_info, - values_array, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - integrate_values, - integrate_gradients, - integrate_hessians, - active_fe_index, - first_selected_component, - cells, - face_nos, - subface_index, - dof_access_index, - face_orientations, - orientation_map); + eval); } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index b171e4344b..81f6aaa3d4 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2309,7 +2309,6 @@ public: void integrate(const EvaluationFlags::EvaluationFlags integration_flag); - /** * @deprecated Please use the integrate() function with the EvaluationFlags argument. */ @@ -2329,7 +2328,8 @@ public: */ void integrate(const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType * values_array); + VectorizedArrayType * values_array, + const bool sum_into_values = false); /** * @deprecated Please use the integrate() function with the EvaluationFlags argument. @@ -2858,6 +2858,8 @@ namespace internal inline std::tuple< const internal::MatrixFreeFunctions::ShapeInfo *, const internal::MatrixFreeFunctions::DoFInfo *, + const internal::MatrixFreeFunctions:: + MappingInfoStorage *, unsigned int, unsigned int> get_shape_info_and_indices( @@ -2895,6 +2897,7 @@ namespace internal active_fe_index, active_quad_index), &dof_info, + &mapping_storage, active_fe_index, active_quad_index); } @@ -2932,10 +2935,6 @@ inline FEEvaluationBase::get( - data_in.get_mapping_info(), - quad_no), quad_no, is_interior_face, face_type) @@ -4992,8 +4991,7 @@ FEEvaluationBase:: for (unsigned int e = 0; e < dim; ++e) hessian_out[comp][d][d] += jac_grad[d][e] * - this->gradients_from_hessians_quad[(comp * dim + e) * nqp + - q_point]; + this->gradients_quad[(comp * dim + e) * nqp + q_point]; // add off-diagonal part of J' * grad(u) for (unsigned int d = 0, count = dim; d < dim; ++d) @@ -5001,8 +4999,7 @@ FEEvaluationBase:: for (unsigned int f = 0; f < dim; ++f) hessian_out[comp][d][e] += jac_grad[count][f] * - this->gradients_from_hessians_quad[(comp * dim + f) * nqp + - q_point]; + this->gradients_quad[(comp * dim + f) * nqp + q_point]; // take symmetric part for (unsigned int d = 0; d < dim; ++d) @@ -5099,8 +5096,7 @@ FEEvaluationBase:: for (unsigned int e = 0; e < dim; ++e) hessian_out[comp][d] += jac_grad[d][e] * - this->gradients_from_hessians_quad[(comp * dim + e) * nqp + - q_point]; + this->gradients_quad[(comp * dim + e) * nqp + q_point]; } } return hessian_out; @@ -7301,28 +7297,13 @@ FEEvaluation -1) SelectEvaluator:: - evaluate(n_components, - evaluation_flag_actual, - *this->data, - const_cast(values_array), - this->values_quad, - this->gradients_quad, - this->hessians_quad, - this->scratch_data.begin()); + evaluate(n_components, evaluation_flag_actual, values_array, *this); else internal::FEEvaluationFactory::evaluate( n_components, evaluation_flag_actual, - *this->data, const_cast(values_array), - this->values_quad, - this->gradients_quad, - this->hessians_quad, - this->scratch_data.begin()); - - if (hessians_on_general_cells) - for (unsigned int i = 0; i < n_components * dim * n_q_points; ++i) - this->gradients_from_hessians_quad[i] = this->gradients_quad[i]; + *this); # ifdef DEBUG if (evaluation_flag_actual & EvaluationFlags::values) @@ -7452,12 +7433,10 @@ namespace internal /** * Implementation for vectors that have the begin() methods. */ - template ::value && std::is_same().begin()), @@ -7465,24 +7444,18 @@ namespace internal VectorType>::type * = nullptr> bool try_integrate_scatter_inplace( + T phi, VectorType & destination, const unsigned int cell, - const unsigned int n_components, const unsigned int active_fe_index, const unsigned int first_selected_component, const internal::MatrixFreeFunctions::DoFInfo *dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * hessians_quad, - VectorizedArrayType * scratch_data, - const internal::MatrixFreeFunctions::ShapeInfo *data, - const EvaluationFlags::EvaluationFlags integration_flag) + const EvaluationFlags::EvaluationFlags evaluation_flag) { // If the index storage is interleaved and contiguous and the vector storage // has the correct alignment, we can directly pass the pointer into the - // vector to the integrate() call, without writing temporary results into a - // separate data field that will later be added into the vector. This saves - // some operations. + // vector to the evaluate() call, without reading the vector entries into a + // separate data field. This saves some operations. if (std::is_same::value && dof_info->index_storage_variants [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] == @@ -7505,28 +7478,8 @@ namespace internal dof_info->component_dof_indices_offset[active_fe_index] [first_selected_component] * VectorizedArrayType::size()); - if (fe_degree > -1) - SelectEvaluator:: - integrate(n_components, - integration_flag, - *data, - vec_values, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - true); - else - FEEvaluationFactory::integrate( - n_components, - integration_flag, - *data, - vec_values, - values_quad, - gradients_quad, - hessians_quad, - scratch_data, - true); + + phi->integrate(evaluation_flag, vec_values, true); return true; } @@ -7537,31 +7490,23 @@ namespace internal /** * Implementation for all other vectors like block vectors. */ - template ::value || !std::is_same().begin()), Number *>::value, VectorType>::type * = nullptr> bool - try_integrate_scatter_inplace( - VectorType &, - const unsigned int, - const unsigned int, - const unsigned int, - const unsigned int, - const internal::MatrixFreeFunctions::DoFInfo *, - const VectorizedArrayType *, - const VectorizedArrayType *, - const VectorizedArrayType *, - const VectorizedArrayType *, - const internal::MatrixFreeFunctions::ShapeInfo *, - const EvaluationFlags::EvaluationFlags) + try_integrate_scatter_inplace(T, + VectorType &, + const unsigned int, + const unsigned int, + const unsigned int, + const internal::MatrixFreeFunctions::DoFInfo *, + const EvaluationFlags::EvaluationFlags) { return false; } @@ -7689,7 +7634,8 @@ FEEvaluation:: integrate(const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType * values_array) + VectorizedArrayType * values_array, + const bool sum_into_values_array) { # ifdef DEBUG if (integration_flag & EvaluationFlags::values) @@ -7734,24 +7680,16 @@ FEEvaluation:: integrate(n_components, integration_flag_actual, - *this->data, values_array, - this->values_quad, - this->gradients_quad, - this->hessians_quad, - this->scratch_data.begin(), - false); + *this, + sum_into_values_array); else internal::FEEvaluationFactory::integrate( n_components, integration_flag_actual, - *this->data, values_array, - this->values_quad, - this->gradients_quad, - this->hessians_quad, - this->scratch_data.begin(), - false); + *this, + sum_into_values_array); # ifdef DEBUG this->dof_values_initialized = true; @@ -7805,22 +7743,13 @@ FEEvaluation( + if (internal::try_integrate_scatter_inplace( + this, destination, this->cell, - n_components, this->active_fe_index, this->first_selected_component, this->dof_info, - this->values_quad, - this->gradients_quad, - this->hessians_quad, - this->scratch_data.begin(), - this->data, integration_flag) == false) { integrate(integration_flag, this->begin_dof_values()); @@ -7930,9 +7859,7 @@ FEFaceEvaluationmapping_data != nullptr, ExcNotInitialized()); - const internal::MatrixFreeFunctions::FaceToCellTopology< - VectorizedArrayType::size()> &faces = - this->matrix_info->get_face_info(face_index); + if (face_index >= this->matrix_info->get_task_info().face_partition_data.back() && face_index < @@ -7943,21 +7870,7 @@ FEFaceEvaluationface_no = - (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no); - this->subface_index = this->is_interior_face == true ? - GeometryInfo::max_children_per_cell : - faces.subface_index; - - // First check if interior or exterior cell has non-standard orientation - // (i.e. the third bit is one or not). Then set zero if this cell has - // standard-orientation else copy the first three bits - // (which is equivalent to modulo 8). See also the documentation of - // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation. - this->face_orientation = - (this->is_interior_face == (faces.face_orientation >= 8)) ? - (faces.face_orientation % 8) : - 0; + this->reinit_face(this->matrix_info->get_face_info(face_index)); this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index]; const unsigned int offsets = @@ -8142,16 +8055,11 @@ FEFaceEvaluationdof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face == false) - { - const auto face_nos = this->compute_face_no_data(); - const auto face_orientations = this->compute_face_orientations(); - -# ifdef DEBUG - // currently on structured meshes are supported -> face numbers and - // orientations have to be the same for all filled lanes - for (unsigned int v = 1; v < VectorizedArrayType::size(); ++v) - { - if (face_nos[v] != numbers::invalid_unsigned_int) - AssertDimension(face_nos[0], face_nos[v]); - if (face_orientations[v] != numbers::invalid_unsigned_int) - AssertDimension(face_orientations[0], face_orientations[v]); - } -# endif - - internal::FEFaceEvaluationImplEvaluateSelector:: - template run( - n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag_actual & EvaluationFlags::values, - evaluation_flag_actual & EvaluationFlags::gradients, - evaluation_flag_actual & EvaluationFlags::hessians, - face_nos[0], - this->subface_index, - face_orientations[0], - this->descriptor->face_orientations); - } + if (fe_degree > -1) + internal::FEFaceEvaluationImplEvaluateSelector:: + template run(n_components, + evaluation_flag_actual, + values_array, + *this); else - { - if (fe_degree > -1) - internal::FEFaceEvaluationImplEvaluateSelector:: - template run( - n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag_actual & EvaluationFlags::values, - evaluation_flag_actual & EvaluationFlags::gradients, - evaluation_flag_actual & EvaluationFlags::hessians, - this->face_no, - this->subface_index, - this->face_orientation, - this->descriptor->face_orientations); - else - internal::FEFaceEvaluationFactory:: - evaluate(n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag_actual & EvaluationFlags::values, - evaluation_flag_actual & EvaluationFlags::gradients, - evaluation_flag_actual & EvaluationFlags::hessians, - this->face_no, - this->subface_index, - this->face_orientation, - this->descriptor->face_orientations); - } - - if (hessians_on_general_cells) - { - unsigned int size = n_components * dim * n_q_points; - for (unsigned int i = 0; i < size; ++i) - this->gradients_from_hessians_quad[i] = this->gradients_quad[i]; - } + internal::FEFaceEvaluationFactory:: + evaluate(n_components, evaluation_flag_actual, values_array, *this); # ifdef DEBUG if (evaluation_flag_actual & EvaluationFlags::values) @@ -8341,26 +8178,21 @@ FEFaceEvaluation:: - integrate(const EvaluationFlags::EvaluationFlags evaluation_flag, + integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType * values_array) { - Assert( - (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients | - EvaluationFlags::hessians)) == 0, - ExcMessage( - "Only EvaluationFlags::values, EvaluationFlags::gradients, and EvaluationFlags::hessians are supported.")); - - if (!(evaluation_flag & EvaluationFlags::values) && - !(evaluation_flag & EvaluationFlags::gradients) && - !(evaluation_flag & EvaluationFlags::hessians)) - return; + Assert((integration_flag & + ~(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians)) == 0, + ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, " + "and EvaluationFlags::hessians are supported.")); - EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag; - if (evaluation_flag & EvaluationFlags::hessians && + EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag; + if (integration_flag & EvaluationFlags::hessians && (this->cell_type > internal::MatrixFreeFunctions::affine)) { unsigned int size = n_components * dim * n_q_points; - if (evaluation_flag & EvaluationFlags::gradients) + if (integration_flag & EvaluationFlags::gradients) { for (unsigned int i = 0; i < size; ++i) this->gradients_quad[i] += this->gradients_from_hessians_quad[i]; @@ -8369,43 +8201,19 @@ FEFaceEvaluationgradients_quad[i] = this->gradients_from_hessians_quad[i]; - evaluation_flag_actual |= EvaluationFlags::gradients; + integration_flag_actual |= EvaluationFlags::gradients; } } if (fe_degree > -1) internal::FEFaceEvaluationImplIntegrateSelector:: - template run( - n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag_actual & EvaluationFlags::values, - evaluation_flag_actual & EvaluationFlags::gradients, - evaluation_flag_actual & EvaluationFlags::hessians, - this->face_no, - this->subface_index, - this->face_orientation, - this->descriptor->face_orientations); + template run(n_components, + integration_flag_actual, + values_array, + *this); else internal::FEFaceEvaluationFactory:: - integrate(n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag_actual & EvaluationFlags::values, - evaluation_flag_actual & EvaluationFlags::gradients, - evaluation_flag_actual & EvaluationFlags::hessians, - this->face_no, - this->subface_index, - this->face_orientation, - this->descriptor->face_orientations); + integrate(n_components, integration_flag_actual, values_array, *this); } @@ -8461,119 +8269,44 @@ FEFaceEvaluationdof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, - this->active_fe_index, - this->dof_info); - - if (this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face == false) - { - const auto cells = this->get_cell_or_face_ids(); - const auto face_nos = this->compute_face_no_data(); - const auto face_orientations = this->compute_face_orientations(); + const auto shared_vector_data = internal::get_shared_vector_data( + input_vector, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); - return internal::FEFaceEvaluationImplGatherEvaluateSelector< + if (fe_degree > 0 && + internal::FEFaceEvaluationImplGatherEvaluateSelector< + dim, + Number, + VectorizedArrayType>::supports(evaluation_flag, + *this->data, + internal::get_beginning( + input_vector), + this->dof_info->index_storage_variants + [this->dof_access_index][this->cell])) + { + if (fe_degree > -1) + internal::FEFaceEvaluationImplGatherEvaluateSelector< dim, Number, VectorizedArrayType>::template run( n_components, - VectorizedArrayType::size(), + evaluation_flag, internal::get_beginning(input_vector), shared_vector_data, - *this->data, - *this->dof_info, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - evaluation_flag & EvaluationFlags::hessians, - this->active_fe_index, - this->first_selected_component, - cells, - face_nos, - this->subface_index, - this->dof_access_index, - face_orientations, - this->descriptor->face_orientations); - } - else - { - // TODO: this copying should not be necessary once we have introduced - // an internal-data structure - std::array cells_ = {}; - std::array face_no_ = {}; - std::array - face_orientation_ = {}; - - cells_[0] = this->cell; - face_no_[0] = this->face_no; - face_orientation_[0] = this->face_orientation; - - if (fe_degree > -1) - { - return internal::FEFaceEvaluationImplGatherEvaluateSelector< - dim, - Number, - VectorizedArrayType>::template run( - n_components, - 1, - internal::get_beginning(input_vector), - shared_vector_data, - *this->data, - *this->dof_info, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - evaluation_flag & EvaluationFlags::hessians, - this->active_fe_index, - this->first_selected_component, - cells_, - face_no_, - this->subface_index, - this->dof_access_index, - face_orientation_, - this->descriptor->face_orientations); - } - else - return internal:: - FEFaceEvaluationFactory:: - gather_evaluate(n_components, - 1, - internal::get_beginning(input_vector), - shared_vector_data, - *this->data, - *this->dof_info, - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - evaluation_flag & EvaluationFlags::hessians, - this->active_fe_index, - this->first_selected_component, - cells_, - face_no_, - this->subface_index, - this->dof_access_index, - face_orientation_, - this->descriptor->face_orientations); - } - }; - - if (!fu()) + *this); + else + internal::FEFaceEvaluationFactory:: + gather_evaluate(n_components, + evaluation_flag, + internal::get_beginning(input_vector), + shared_vector_data, + *this); + } + else { this->read_dof_values(input_vector); this->evaluate(evaluation_flag); @@ -8633,7 +8366,7 @@ FEFaceEvaluation:: - integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag, + integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType & destination) { Assert((this->dof_access_index == @@ -8648,81 +8381,40 @@ FEFaceEvaluationactive_fe_index, this->dof_info); - // TODO: this copying should not be necessary once we have introduced - // an internal-data structure - std::array cells_ = {}; - std::array face_no_ = {}; - std::array face_orientation_ = {}; - - cells_[0] = this->cell; - face_no_[0] = this->face_no; - face_orientation_[0] = this->face_orientation; - - if (fe_degree > -1) + if (fe_degree > 0 && + internal::FEFaceEvaluationImplGatherEvaluateSelector< + dim, + Number, + VectorizedArrayType>::supports(integration_flag, + *this->data, + internal::get_beginning( + destination), + this->dof_info->index_storage_variants + [this->dof_access_index][this->cell])) { - if (!internal::FEFaceEvaluationImplIntegrateScatterSelector< - dim, - Number, - VectorizedArrayType>::template run( - n_components, - 1, - internal::get_beginning(destination), - shared_vector_data, - *this->data, - *this->dof_info, - this->begin_dof_values(), - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - evaluation_flag & EvaluationFlags::hessians, - this->active_fe_index, - this->first_selected_component, - cells_, - face_no_, - this->subface_index, - this->dof_access_index, - face_orientation_, - this->descriptor->face_orientations)) - { - // if we arrive here, writing into the destination vector did not - // succeed because some of the assumptions in integrate_scatter were - // not fulfilled (e.g. an element or degree that does not support - // direct writing), so we must do it here - this->distribute_local_to_global(destination); - } + if (fe_degree > -1) + internal::FEFaceEvaluationImplIntegrateScatterSelector< + dim, + Number, + VectorizedArrayType>::template run( + n_components, + integration_flag, + internal::get_beginning(destination), + shared_vector_data, + *this); + else + internal::FEFaceEvaluationFactory:: + integrate_scatter(n_components, + integration_flag, + internal::get_beginning(destination), + shared_vector_data, + *this); } else { - if (!internal::FEFaceEvaluationFactory:: - integrate_scatter(n_components, - 1, - internal::get_beginning(destination), - shared_vector_data, - *this->data, - *this->dof_info, - this->begin_dof_values(), - this->begin_values(), - this->begin_gradients(), - this->begin_hessians(), - this->scratch_data.begin(), - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - evaluation_flag & EvaluationFlags::hessians, - this->active_fe_index, - this->first_selected_component, - cells_, - face_no_, - this->subface_index, - this->dof_access_index, - face_orientation_, - this->descriptor->face_orientations)) - { - this->distribute_local_to_global(destination); - } + integrate(integration_flag); + this->distribute_local_to_global(destination); } } diff --git a/include/deal.II/matrix_free/fe_evaluation_base_data.h b/include/deal.II/matrix_free/fe_evaluation_base_data.h index ab1d057382..bf7457a1e5 100644 --- a/include/deal.II/matrix_free/fe_evaluation_base_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_base_data.h @@ -125,6 +125,22 @@ class FEEvaluationBaseData public: static constexpr unsigned int dimension = dim; + using NumberType = VectorizedArrayType; + + /** + * Constructor, taking data stored in MatrixFree. If applied to problems + * with more than one quadrature formula selected during construction of + * `matrix_free`, `quad_no` allows to select the appropriate formula. + */ + FEEvaluationBaseData(const std::tuple &info, + const unsigned int quad_no, + const bool is_interior_face, + const unsigned int face_type); + /** * Copy constructor. */ @@ -136,6 +152,23 @@ public: FEEvaluationBaseData & operator=(const FEEvaluationBaseData &other); + /** + * Sets the pointers for values, gradients, hessians to the central + * scratch_data_array inside the given scratch array, for a given number of + * components as provided by one of the derived classes. + */ + void + set_data_pointers(AlignedVector *scratch_data, + const unsigned int n_components); + + /** + * Initialize the indices related to the given face description. Used during + * the reinit() call of the FEFaceEvaluation class. + */ + void + reinit_face(const internal::MatrixFreeFunctions::FaceToCellTopology< + VectorizedArrayType::size()> &face); + /** * @name 1: Access to geometry data at quadrature points */ @@ -371,26 +404,78 @@ public: unsigned int get_active_quadrature_index() const; - //@} + /** + * If is_face is true, this function returns the face number within a cell for + * the face this object was initialized to. On cells where the face makes no + * sense, an exception is thrown. + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. + */ + unsigned int + get_face_no() const; -protected: /** - * Constructor. Made protected to prevent users from directly using this - * class. Takes all data stored in MatrixFree. If applied to problems with - * more than one quadrature formula selected during construction of - * `matrix_free`, `quad_no` allows to select the appropriate formula. + * If is_face is true, this function returns the index of a subface along a + * face in case this object was initialized to a face with hanging nodes. On + * faces On cells where the face makes no sense, an exception is thrown. + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. */ - FEEvaluationBaseData(const std::tuple &info, - const MappingInfoStorageType & mapping_data, - const unsigned int quad_no, - const bool is_interior_face, - const unsigned int face_type); + unsigned int + get_subface_index() const; + + /** + * If is_face is true, this function returns the orientation index within an + * array of orientations as returned by get_orientation_map(). + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. + */ + unsigned int + get_face_orientation() const; /** - * Constructor that comes with reduced functionality and works similar as + * If is_face is true, this function returns the map re-ordering from face + * quadrature points in standard orientation to a specific orientation. + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. + */ + const Table<2, unsigned int> & + get_orientation_map() const; + + /** + * Return the current index in the access to compressed indices. + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. + */ + internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex + get_dof_access_index() const; + + /** + * Return whether this object was set up to work on what is called + * "interior" faces in the evaluation context. + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. + */ + bool + get_is_interior_face() const; + + //@} + +protected: + /** + * Constructor that comes with reduced functionality and works similarly as * FEValues. */ FEEvaluationBaseData( @@ -398,15 +483,6 @@ protected: internal::MatrixFreeFunctions:: MappingDataOnTheFly> &mapping_data); - /** - * Sets the pointers for values, gradients, hessians to the central - * scratch_data_array inside the given scratch array, for a given number of - * components as provided by one of the derived classes. - */ - void - set_data_pointers(AlignedVector *scratch_data, - const unsigned int n_components); - /** * A pointer to the unit cell shape data, i.e., values, gradients and * Hessians in 1D at the quadrature points that constitute the tensor @@ -683,20 +759,20 @@ template inline FEEvaluationBaseData:: FEEvaluationBaseData(const std::tuple &shape_dof_info, - const MappingInfoStorageType & mapping_data, const unsigned int quad_no, const bool is_interior_face, const unsigned int face_type) : data(std::get<0>(shape_dof_info)) , dof_info(std::get<1>(shape_dof_info)) - , mapping_data(&mapping_data) + , mapping_data(std::get<2>(shape_dof_info)) , quad_no(quad_no) - , active_fe_index(std::get<2>(shape_dof_info)) - , active_quad_index(std::get<3>(shape_dof_info)) + , active_fe_index(std::get<3>(shape_dof_info)) + , active_quad_index(std::get<4>(shape_dof_info)) , descriptor( - &mapping_data.descriptor + &mapping_data->descriptor [is_face ? (active_quad_index * std::max(1, dim - 1) + (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) : @@ -1081,6 +1157,67 @@ FEEvaluationBaseData:: +template +inline unsigned int +FEEvaluationBaseData::get_face_no() + const +{ + return face_no; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_subface_index() const +{ + return subface_index; +} + + + +template +inline unsigned int +FEEvaluationBaseData:: + get_face_orientation() const +{ + return face_orientation; +} + + + +template +inline const Table<2, unsigned int> & +FEEvaluationBaseData:: + get_orientation_map() const +{ + Assert(descriptor != nullptr, ExcNotInitialized()); + return descriptor->face_orientations; +} + + + +template +inline internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex +FEEvaluationBaseData:: + get_dof_access_index() const +{ + return dof_access_index; +} + + + +template +inline bool +FEEvaluationBaseData:: + get_is_interior_face() const +{ + return is_interior_face; +} + + + template inline void FEEvaluationBaseData:: @@ -1122,6 +1259,29 @@ FEEvaluationBaseData:: +template +inline void +FEEvaluationBaseData::reinit_face( + const internal::MatrixFreeFunctions::FaceToCellTopology< + VectorizedArrayType::size()> &face) +{ + face_no = (is_interior_face ? face.interior_face_no : face.exterior_face_no); + subface_index = is_interior_face == true ? + GeometryInfo::max_children_per_cell : + face.subface_index; + + // First check if interior or exterior cell has non-standard orientation + // (i.e. the third bit is one or not). Then set zero if this cell has + // standard-orientation else copy the first three bits + // (which is equivalent to modulo 8). See also the documentation of + // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation. + face_orientation = (is_interior_face == (face.face_orientation >= 8)) ? + (face.face_orientation % 8) : + 0; +} + + + #endif // ifndef DOXYGEN diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index fc29f2a9d4..cf562d5740 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -1175,13 +1176,13 @@ namespace internal shape_info.dofs_per_component_on_cell; constexpr unsigned int hess_dim = dim * (dim + 1) / 2; - AlignedVector cell_points(dim * n_mapping_points); - AlignedVector cell_quads(dim * n_q_points); - AlignedVector cell_grads(dim * dim * n_q_points); - AlignedVector cell_grad_grads(dim * hess_dim * - n_q_points); - AlignedVector scratch_data( - dim * (2 * n_q_points + 3 * n_mapping_points)); + MappingInfoStorage temp_data; + temp_data.copy_descriptor(my_data); + FEEvaluationBaseData eval( + {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0); + + AlignedVector evaluation_data; + eval.set_data_pointers(&evaluation_data, dim); for (unsigned int cell = begin_cell; cell < end_cell; ++cell) for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d) @@ -1195,7 +1196,7 @@ namespace internal vectorized_load_and_transpose(n_mapping_points * dim, plain_quadrature_points.data(), start_indices, - cell_points.data()); + eval.begin_dof_values()); FEEvaluationFactory::evaluate( dim, @@ -1203,12 +1204,8 @@ namespace internal (update_flags_cells & update_jacobian_grads ? EvaluationFlags::hessians : EvaluationFlags::nothing), - shape_info, - cell_points.data(), - cell_quads.data(), - cell_grads.data(), - cell_grad_grads.data(), - scratch_data.data()); + eval.begin_dof_values(), + eval); } if (update_flags_cells & update_quadrature_points) { @@ -1225,9 +1222,10 @@ namespace internal else for (unsigned int d = 0; d < dim; ++d) for (unsigned int q = 0; q < n_q_points; ++q) - store_vectorized_array(cell_quads[q + d * n_q_points], - vv, - quadrature_points[q][d]); + store_vectorized_array( + eval.begin_values()[q + d * n_q_points], + vv, + quadrature_points[q][d]); } const unsigned int n_points = @@ -1240,7 +1238,9 @@ namespace internal Tensor<2, dim, VectorizedDouble> jac; for (unsigned int d = 0; d < dim; ++d) for (unsigned int e = 0; e < dim; ++e) - jac[d][e] = cell_grads[q + (d * dim + e) * n_q_points]; + jac[d][e] = + eval + .begin_gradients()[q + (d * dim + e) * n_q_points]; // eliminate roundoff errors if (cell_type[cell] == cartesian) @@ -1288,13 +1288,13 @@ namespace internal { for (unsigned int e = 0; e < dim; ++e) jac_grad[d][e][e] = - cell_grad_grads[q + (d * hess_dim + e) * - n_q_points]; + eval.begin_hessians()[q + (d * hess_dim + e) * + n_q_points]; for (unsigned int c = dim, e = 0; e < dim; ++e) for (unsigned int f = e + 1; f < dim; ++f, ++c) jac_grad[d][e][f] = jac_grad[d][f][e] = - cell_grad_grads[q + (d * hess_dim + c) * - n_q_points]; + eval.begin_hessians()[q + (d * hess_dim + c) * + n_q_points]; const auto inv_jac_grad = process_jacobian_gradient(inv_jac, inv_jac, @@ -2088,21 +2088,30 @@ namespace internal shape_info.dofs_per_component_on_cell; constexpr unsigned int hess_dim = dim * (dim + 1) / 2; - AlignedVector cell_points(dim * n_mapping_points); - AlignedVector face_quads(dim * n_q_points); - AlignedVector face_grads(dim * dim * n_q_points); - AlignedVector face_grad_grads(dim * hess_dim * - n_q_points); + MappingInfoStorage temp_data; + temp_data.copy_descriptor(my_data); + FEEvaluationBaseData eval_int( + {&shape_info, nullptr, &temp_data, 0, 0}, 0, true, 0); + FEEvaluationBaseData eval_ext( + {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0); - AlignedVector scratch_data( - dim * (2 * n_q_points + 3 * n_mapping_points)); + AlignedVector evaluation_data, evaluation_data_ext; + eval_int.set_data_pointers(&evaluation_data, dim); + eval_ext.set_data_pointers(&evaluation_data_ext, dim); for (unsigned int face = begin_face; face < end_face; ++face) for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d) { + internal::MatrixFreeFunctions::FaceToCellTopology< + VectorizedDouble::size()> + face_double; + face_double.interior_face_no = faces[face].interior_face_no; + face_double.exterior_face_no = faces[face].exterior_face_no; + face_double.face_orientation = faces[face].face_orientation; + face_double.subface_index = faces[face].subface_index; + // load the geometry field for all SIMD lanes - unsigned int start_indices[n_lanes_d]; - const unsigned int face_no = faces[face].interior_face_no; + unsigned int start_indices[n_lanes_d]; for (unsigned int v = 0; v < n_lanes_d; ++v) if (faces[face].cells_interior[vv + v] != numbers::invalid_unsigned_int) @@ -2111,37 +2120,29 @@ namespace internal else start_indices[v] = faces[face].cells_interior[0] * n_mapping_points * dim; + + eval_int.reinit_face(face_double); vectorized_load_and_transpose(n_mapping_points * dim, plain_quadrature_points.data(), start_indices, - cell_points.data()); + eval_int.begin_dof_values()); // now let the matrix-free evaluators provide us with the // data on faces FEFaceEvaluationFactory::evaluate( dim, - shape_info, - cell_points.data(), - face_quads.data(), - face_grads.data(), - face_grad_grads.data(), - scratch_data.data(), - true, - true, - update_flags_faces & update_jacobian_grads, - face_no, - GeometryInfo::max_children_per_cell, - faces[face].face_orientation > 8 ? - faces[face].face_orientation - 8 : - 0, - my_data.descriptor[0].face_orientations); - + EvaluationFlags::values | EvaluationFlags::gradients | + (update_flags_faces & update_jacobian_grads ? + EvaluationFlags::hessians : + EvaluationFlags::nothing), + eval_int.begin_dof_values(), + eval_int); if (update_flags_faces & update_quadrature_points) for (unsigned int q = 0; q < n_q_points; ++q) for (unsigned int d = 0; d < dim; ++d) store_vectorized_array( - face_quads[d * n_q_points + q], + eval_int.begin_values()[d * n_q_points + q], vv, my_data.quadrature_points [my_data.quadrature_point_offsets[face] + q][d]); @@ -2155,7 +2156,9 @@ namespace internal [&](unsigned int face_no, int is_exterior, unsigned int q, - Tensor<2, dim, VectorizedDouble> inv_jac) { + Tensor<2, dim, VectorizedDouble> inv_jac, + FEEvaluationBaseData + &eval) { Tensor<2, dim, VectorizedDouble> inv_transp_jac_permut; for (unsigned int d = 0; d < dim; ++d) for (unsigned int e = 0; e < dim; ++e) @@ -2168,7 +2171,8 @@ namespace internal Tensor<2, dim, VectorizedDouble> jacobi; for (unsigned int e = 0; e < dim; ++e) for (unsigned int d = 0; d < dim; ++d) - jacobi[d][e] = face_grads[(d * dim + e) * n_q_points + q]; + jacobi[d][e] = + eval.begin_gradients()[(d * dim + e) * n_q_points + q]; Tensor<2, dim, VectorizedDouble> inv_transp_jac = transpose(invert(jacobi)); Tensor<3, dim, VectorizedDouble> jac_grad; @@ -2176,12 +2180,13 @@ namespace internal { for (unsigned int e = 0; e < dim; ++e) jac_grad[d][e][e] = - face_grad_grads[q + (d * hess_dim + e) * n_q_points]; + eval.begin_gradients()[q + (d * hess_dim + e) * + n_q_points]; for (unsigned int c = dim, e = 0; e < dim; ++e) for (unsigned int f = e + 1; f < dim; ++f, ++c) jac_grad[d][e][f] = jac_grad[d][f][e] = - face_grad_grads[q + - (d * hess_dim + c) * n_q_points]; + eval.begin_gradients()[q + (d * hess_dim + c) * + n_q_points]; const auto inv_jac_grad = process_jacobian_gradient(inv_transp_jac_permut, inv_transp_jac, @@ -2203,6 +2208,7 @@ namespace internal face_type[face] <= affine ? 1 : n_q_points; for (unsigned int q = 0; q < n_points_compute; ++q) { + const unsigned int face_no = faces[face].interior_face_no; Tensor<2, dim, VectorizedDouble> jac; for (unsigned int e = 0; e < dim; ++e) { @@ -2210,7 +2216,9 @@ namespace internal ExtractFaceHelper::reorder_face_derivative_indices( face_no, e); for (unsigned int d = 0; d < dim; ++d) - jac[d][ee] = face_grads[(d * dim + e) * n_q_points + q]; + jac[d][ee] = + eval_int + .begin_gradients()[(d * dim + e) * n_q_points + q]; } Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac); for (unsigned int e = 0; e < dim; ++e) @@ -2227,7 +2235,7 @@ namespace internal if (update_flags_faces & update_jacobian_grads) { - compute_jacobian_grad(face_no, 0, q, inv_jac); + compute_jacobian_grad(face_no, 0, q, inv_jac, eval_int); } std::array, dim - 1> @@ -2285,47 +2293,43 @@ namespace internal start_indices[v] = faces[face].cells_exterior[0] * n_mapping_points * dim; + eval_ext.reinit_face(face_double); vectorized_load_and_transpose(n_mapping_points * dim, plain_quadrature_points.data(), start_indices, - cell_points.data()); + eval_ext.begin_dof_values()); FEFaceEvaluationFactory:: evaluate(dim, - shape_info, - cell_points.data(), - face_quads.data(), - face_grads.data(), - face_grad_grads.data(), - scratch_data.data(), - false, - true, - update_flags_faces & update_jacobian_grads, - faces[face].exterior_face_no, - faces[face].subface_index, - faces[face].face_orientation < 8 ? - faces[face].face_orientation : - 0, - my_data.descriptor[0].face_orientations); + EvaluationFlags::values | + EvaluationFlags::gradients | + (update_flags_faces & update_jacobian_grads ? + EvaluationFlags::hessians : + EvaluationFlags::nothing), + eval_ext.begin_dof_values(), + eval_ext); for (unsigned int q = 0; q < n_points_compute; ++q) { + const unsigned int face_no = faces[face].exterior_face_no; Tensor<2, dim, VectorizedDouble> jac; for (unsigned int e = 0; e < dim; ++e) { const unsigned int ee = ExtractFaceHelper::reorder_face_derivative_indices< - dim>(faces[face].exterior_face_no, e); + dim>(face_no, e); for (unsigned int d = 0; d < dim; ++d) jac[d][ee] = - face_grads[(d * dim + e) * n_q_points + q]; + eval_ext + .begin_gradients()[(d * dim + e) * n_q_points + + q]; } Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac); for (unsigned int e = 0; e < dim; ++e) { const unsigned int ee = ExtractFaceHelper::reorder_face_derivative_indices< - dim>(faces[face].exterior_face_no, e); + dim>(face_no, e); for (unsigned int d = 0; d < dim; ++d) store_vectorized_array( inv_jac[ee][d], @@ -2335,10 +2339,8 @@ namespace internal if (update_flags_faces & update_jacobian_grads) { - compute_jacobian_grad(faces[face].exterior_face_no, - 1, - q, - inv_jac); + compute_jacobian_grad( + face_no, 1, q, inv_jac, eval_ext); } my_data.normals_times_jacobians[1][offset + q] = diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h index d121e96a06..4457046663 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.h +++ b/include/deal.II/matrix_free/mapping_info_storage.h @@ -304,6 +304,17 @@ namespace internal unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const; + /** + * Get the descriptor from another instance of QuadratureDescriptor, + * clearing all other data fields. + */ + template + void + copy_descriptor(const MappingInfoStorage &other); + /** * Prints a detailed summary of memory consumption in the different * structures of this class to the given output stream. diff --git a/include/deal.II/matrix_free/mapping_info_storage.templates.h b/include/deal.II/matrix_free/mapping_info_storage.templates.h index 69509d1024..aff20a206e 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.templates.h +++ b/include/deal.II/matrix_free/mapping_info_storage.templates.h @@ -178,6 +178,47 @@ namespace internal + template + template + void + MappingInfoStorage:: + copy_descriptor(const MappingInfoStorage &other) + { + clear_data_fields(); + descriptor.clear(); + descriptor.resize(other.descriptor.size()); + for (unsigned int i = 0; i < descriptor.size(); ++i) + { + descriptor[i].n_q_points = other.descriptor[i].n_q_points; + descriptor[i].quadrature_1d = other.descriptor[i].quadrature_1d; + descriptor[i].quadrature = other.descriptor[i].quadrature; + for (unsigned int d = 0; d < structdim; ++d) + { + descriptor[i].tensor_quadrature_weights[d].resize( + other.descriptor[i].tensor_quadrature_weights[d].size()); + std::copy( + other.descriptor[i].tensor_quadrature_weights[d].begin(), + other.descriptor[i].tensor_quadrature_weights[d].end(), + descriptor[i].tensor_quadrature_weights[d].begin()); + } + descriptor[i].quadrature_weights.resize( + other.descriptor[i].quadrature_weights.size()); + std::copy(other.descriptor[i].quadrature_weights.begin(), + other.descriptor[i].quadrature_weights.end(), + descriptor[i].quadrature_weights.begin()); + descriptor[i].face_orientations = + other.descriptor[i].face_orientations; + } + } + + + template