From d10cd7e0dbc95372f0b06237679594655ec4e9e0 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 10 Dec 2021 15:13:32 +0100 Subject: [PATCH] Remove second template argument in FEEvaluationBaseData Rename FEEvaluationBaseData to FEEvaluationData. --- include/deal.II/fe/mapping_q_internal.h | 10 +- .../deal.II/matrix_free/evaluation_kernels.h | 310 +++++++------- .../deal.II/matrix_free/evaluation_selector.h | 25 +- .../matrix_free/evaluation_template_factory.h | 116 +++-- .../evaluation_template_factory.templates.h | 155 ++++--- include/deal.II/matrix_free/fe_evaluation.h | 140 +++--- ...ation_base_data.h => fe_evaluation_data.h} | 397 ++++++++---------- .../matrix_free/mapping_data_on_the_fly.h | 77 ++-- include/deal.II/matrix_free/mapping_info.h | 26 +- .../matrix_free/mapping_info.templates.h | 145 +++---- .../matrix_free/mapping_info_storage.h | 64 +-- .../mapping_info_storage.templates.h | 129 +++--- include/deal.II/matrix_free/operators.h | 20 +- .../evaluation_template_factory.inst.in | 17 +- source/matrix_free/mapping_info.inst.in | 2 - 15 files changed, 746 insertions(+), 887 deletions(-) rename include/deal.II/matrix_free/{fe_evaluation_base_data.h => fe_evaluation_data.h} (76%) diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 1fed26230c..49c4268f4f 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -35,7 +35,7 @@ #include #include -#include +#include #include #include #include @@ -1131,11 +1131,11 @@ namespace internal } internal::MatrixFreeFunctions:: - MappingInfoStorage + MappingInfoStorage temp_data; temp_data.descriptor.resize(1); temp_data.descriptor[0].n_q_points = n_q_points; - FEEvaluationBaseData eval( + FEEvaluationData eval( {&data.shape_info, nullptr, &temp_data, 0, 0}); // prepare arrays @@ -1156,8 +1156,8 @@ namespace internal } // do the actual tensorized evaluation - internal::FEEvaluationFactory:: - evaluate(n_comp, evaluation_flag, eval.begin_dof_values(), eval); + internal::FEEvaluationFactory::evaluate( + n_comp, evaluation_flag, eval.begin_dof_values(), eval); } // do the postprocessing diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 576187616b..f2b01daf1d 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -26,11 +26,10 @@ #include #include -#include +#include #include #include #include -#include DEAL_II_NAMESPACE_OPEN @@ -121,20 +120,17 @@ namespace internal n_q_points_1d, Number>; - using EvaluationData = - FEEvaluationBaseData; - static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs_actual, - EvaluationData & eval); + FEEvaluationData & eval); static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs_actual, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array); static Eval @@ -170,20 +166,17 @@ namespace internal n_q_points_1d, Number> { - using EvaluationData = - FEEvaluationBaseData; - static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs_actual, - EvaluationData & eval); + FEEvaluationData & eval); static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs_actual, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array); }; @@ -199,7 +192,7 @@ namespace internal const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs_actual, - EvaluationData & eval) + FEEvaluationData & eval) { if (evaluation_flag == EvaluationFlags::nothing) return; @@ -249,10 +242,10 @@ namespace internal const Number *values_dofs = values_dofs_actual; if (type == MatrixFreeFunctions::truncated_tensor) { - const std::size_t stored_dofs_per_comp = + const std::size_t n_dofs_per_comp = eval.get_shape_info().dofs_per_component_on_cell; Number *values_dofs_tmp = - temp1 + 2 * (std::max(stored_dofs_per_comp, n_q_points)); + temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points)); const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree; for (unsigned int c = 0; c < n_components; ++c) @@ -265,7 +258,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 * stored_dofs_per_comp + count_p]; + values_dofs_actual[c * n_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(); @@ -472,7 +465,7 @@ namespace internal const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs_actual, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array) { std::array *, 3> @@ -509,7 +502,7 @@ namespace internal temp2 = temp1 + temp_size; } - const unsigned int n_q_points = temp_size == 0 ? + const std::size_t n_q_points = temp_size == 0 ? eval.get_shape_info().n_q_points : Eval::n_columns_of_product; const unsigned int dofs_per_comp = @@ -519,8 +512,9 @@ namespace internal // expand dof_values to tensor product for truncated tensor products Number *values_dofs = (type == MatrixFreeFunctions::truncated_tensor) ? - temp1 + 2 * (std::max(eval.get_shape_info().dofs_per_component_on_cell, - 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(); @@ -750,7 +744,7 @@ namespace internal if (type == MatrixFreeFunctions::truncated_tensor) { - const std::size_t stored_dofs_per_comp = + const std::size_t n_dofs_per_comp = eval.get_shape_info().dofs_per_component_on_cell; values_dofs -= dofs_per_comp * n_components; const int degree = @@ -764,7 +758,7 @@ namespace internal { for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q) - values_dofs_actual[c * stored_dofs_per_comp + count_p] = + values_dofs_actual[c * n_dofs_per_comp + count_p] = values_dofs[c * dofs_per_comp + count_q]; count_q += j + i; } @@ -783,13 +777,12 @@ namespace internal n_q_points_1d, Number>::evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags - evaluation_flag, - const Number * values_dofs_actual, - EvaluationData &eval) + evaluation_flag, + const Number *values_dofs_actual, + FEEvaluationData &eval) { - 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 std::size_t n_dofs = eval.get_shape_info().dofs_per_component_on_cell; + const std::size_t n_q_points = eval.get_shape_info().n_q_points; const auto &shape_data = eval.get_shape_info().data; @@ -846,24 +839,23 @@ namespace internal template inline void - FEEvaluationImpl::integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags - integration_flag, - Number * values_dofs_actual, - EvaluationData &eval, - const bool add_into_values_array) + FEEvaluationImpl< + MatrixFreeFunctions::tensor_none, + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs_actual, + FEEvaluationData &eval, + const bool add_into_values_array) { // TODO: implement hessians AssertThrow(!(integration_flag & EvaluationFlags::hessians), ExcNotImplemented()); - 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 std::size_t n_dofs = eval.get_shape_info().dofs_per_component_on_cell; + const std::size_t n_q_points = eval.get_shape_info().n_q_points; const auto &shape_data = eval.get_shape_info().data; @@ -1340,14 +1332,11 @@ 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); + FEEvaluationData & eval); static void do_evaluate(const MatrixFreeFunctions::UnivariateShapeData &shape, @@ -1360,7 +1349,7 @@ namespace internal integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array); static void @@ -1380,7 +1369,7 @@ namespace internal const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, - EvaluationData & eval) + FEEvaluationData & eval) { constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); @@ -1468,7 +1457,7 @@ namespace internal const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array) { constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); @@ -1617,20 +1606,17 @@ namespace internal template struct FEEvaluationImplTransformToCollocation { - using EvaluationData = - FEEvaluationBaseData; - static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, - EvaluationData & eval); + FEEvaluationData & eval); static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number * values_dofs, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array); }; @@ -1645,7 +1631,7 @@ namespace internal Number>::evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, - EvaluationData & eval) + FEEvaluationData & eval) { const auto &shape_data = eval.get_shape_info().data[0]; @@ -1654,6 +1640,7 @@ namespace internal "of lower degree, so the evaluation results would be " "wrong. Thus, this class does not permit the desired " "operation.")); + constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim); constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim); for (unsigned int c = 0; c < n_components; ++c) @@ -1667,8 +1654,7 @@ namespace internal Number, Number>::do_forward(1, shape_data.shape_values_eo, - values_dofs + - c * Utilities::pow(fe_degree + 1, dim), + values_dofs + c * n_dofs, eval.begin_values() + c * n_q_points); // apply derivatives in the collocation space @@ -1696,7 +1682,7 @@ namespace internal Number>::integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, - EvaluationData & eval, + FEEvaluationData & eval, const bool add_into_values_array) { const auto &shape_data = eval.get_shape_info().data[0]; @@ -1706,7 +1692,7 @@ namespace internal "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) { @@ -1743,6 +1729,22 @@ namespace internal + /** + * Helper function to specify whether transformation to collocation should + * be used: It should give correct results (first condition), we need to be + * able to initialize the fields in shape_info.templates.h from the + * polynomials (second condition), and it should be the most efficient + * choice in terms of operation counts (third condition). + */ + constexpr bool + use_collocation_evaluation(const unsigned int fe_degree, + const unsigned int n_q_points_1d) + { + return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) && + (n_q_points_1d <= 3 * fe_degree / 2 + 1); + } + + /** * This class chooses an appropriate evaluation strategy based on the * template parameters and the shape_info variable which contains runtime @@ -1766,21 +1768,11 @@ namespace internal run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, - FEEvaluationBaseData - &eval) + FEEvaluationData & 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). - 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), @@ -1794,7 +1786,8 @@ namespace internal } // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see // shape_info.h for more details - else if (fe_degree >= 0 && use_collocation && + else if (fe_degree >= 0 && + use_collocation_evaluation(fe_degree, n_q_points_1d) && element_type <= ElementType::tensor_symmetric) { internal::FEEvaluationImplTransformToCollocation< @@ -1888,22 +1881,12 @@ namespace internal run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, - FEEvaluationBaseData - & eval, - const bool sum_into_values_array) + FEEvaluationData & 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 < 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), @@ -1921,7 +1904,8 @@ namespace internal } // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see // shape_info.h for more details - else if (fe_degree >= 0 && use_collocation && + else if (fe_degree >= 0 && + use_collocation_evaluation(fe_degree, n_q_points_1d) && element_type <= ElementType::tensor_symmetric) { internal::FEEvaluationImplTransformToCollocation< @@ -2013,10 +1997,6 @@ namespace internal // choice in terms of operation counts (third condition) and if we were // able to initialize the fields in shape_info.templates.h from the // polynomials (fourth condition). - static constexpr bool use_collocation = - 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, @@ -2073,9 +2050,9 @@ namespace internal 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) : - Utilities::pow(data.fe_degree + 1, dim - 1); + const std::size_t n_dofs = fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(data.fe_degree + 1, dim - 1); const std::size_t n_q_points = fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : Utilities::pow(data.n_q_points_1d, dim - 1); @@ -2105,7 +2082,9 @@ namespace internal default: Assert(false, ExcNotImplemented()); } - values_dofs += 3 * size_deg; + // Note: we always keep storage of values, 1st and 2nd derivatives + // in an array + values_dofs += 3 * n_dofs; values_quad += n_q_points; } else if (evaluation_flag & EvaluationFlags::gradients) @@ -2114,7 +2093,8 @@ namespace internal switch (dim) { case 3: - if (use_collocation) + if (symmetric_evaluate && + use_collocation_evaluation(fe_degree, n_q_points_1d)) { eval0.template values<0, true, false>(values_dofs, values_quad); @@ -2154,14 +2134,14 @@ namespace internal values_quad); } // grad z - eval0.template values<0, true, false>(values_dofs + size_deg, + eval0.template values<0, true, false>(values_dofs + n_dofs, scratch_data); eval1.template values<1, true, false>( scratch_data, gradients_quad + (dim - 1) * n_q_points); break; case 2: - eval0.template values<0, true, false>(values_dofs + size_deg, + eval0.template values<0, true, false>(values_dofs + n_dofs, gradients_quad + n_q_points); eval0.template gradients<0, true, false>(values_dofs, @@ -2177,7 +2157,7 @@ namespace internal default: AssertThrow(false, ExcNotImplemented()); } - values_dofs += 3 * size_deg; + values_dofs += 3 * n_dofs; values_quad += n_q_points; gradients_quad += dim * n_q_points; } @@ -2205,7 +2185,7 @@ namespace internal // grad zz eval0.template values<0, true, false>(values_dofs + - 2 * size_deg, + 2 * n_dofs, scratch_data); eval1.template values<1, true, false>(scratch_data, hessians_quad + @@ -2220,15 +2200,14 @@ namespace internal // grad xz eval0.template gradients<0, true, false>(values_dofs + - size_deg, + n_dofs, scratch_data); eval1.template values<1, true, false>(scratch_data, hessians_quad + 4 * n_q_points); // grad yz - eval0.template values<0, true, false>(values_dofs + - size_deg, + eval0.template values<0, true, false>(values_dofs + n_dofs, scratch_data); eval1.template gradients<1, true, false>(scratch_data, hessians_quad + @@ -2241,10 +2220,10 @@ namespace internal hessians_quad); // grad yy eval0.template values<0, true, false>( - values_dofs + 2 * size_deg, hessians_quad + n_q_points); + values_dofs + 2 * n_dofs, hessians_quad + n_q_points); // grad xy eval0.template gradients<0, true, false>( - values_dofs + size_deg, hessians_quad + 2 * n_q_points); + values_dofs + n_dofs, hessians_quad + 2 * n_q_points); break; case 1: hessians_quad[0] = values_dofs[2]; @@ -2252,7 +2231,7 @@ namespace internal default: AssertThrow(false, ExcNotImplemented()); } - values_dofs += 3 * size_deg; + values_dofs += 3 * n_dofs; hessians_quad += dim * (dim + 1) / 2 * n_q_points; } } @@ -2273,7 +2252,7 @@ namespace internal 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 = + const std::size_t n_dofs = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : (dim > 1 ? Utilities::fixed_power(data.fe_degree + 1) : 1); @@ -2306,7 +2285,7 @@ namespace internal default: Assert(false, ExcNotImplemented()); } - values_dofs += 3 * size_deg; + values_dofs += 3 * n_dofs; values_quad += n_q_points; } else if (integration_flag & EvaluationFlags::gradients) @@ -2320,9 +2299,11 @@ namespace internal 2 * n_q_points, gradients_quad + 2 * n_q_points); - eval0.template values<0, false, false>( - gradients_quad + 2 * n_q_points, values_dofs + size_deg); - if (use_collocation) + eval0.template values<0, false, false>(gradients_quad + + 2 * n_q_points, + values_dofs + n_dofs); + if (symmetric_evaluate && + use_collocation_evaluation(fe_degree, n_q_points_1d)) { internal::EvaluatorTensorProduct< internal::evaluate_evenodd, @@ -2371,8 +2352,9 @@ namespace internal } break; case 2: - eval0.template values<0, false, false>( - gradients_quad + n_q_points, values_dofs + size_deg); + eval0.template values<0, false, false>(gradients_quad + + n_q_points, + values_dofs + n_dofs); eval0.template gradients<0, false, false>(gradients_quad, values_dofs); if (integration_flag & EvaluationFlags::values) @@ -2386,7 +2368,7 @@ namespace internal default: AssertThrow(false, ExcNotImplemented()); } - values_dofs += 3 * size_deg; + values_dofs += 3 * n_dofs; values_quad += n_q_points; gradients_quad += dim * n_q_points; } @@ -2423,7 +2405,7 @@ namespace internal scratch_data); eval0.template values<0, false, false>(scratch_data, values_dofs + - 2 * size_deg); + 2 * n_dofs); // grad xy eval1.template gradients<1, false, false>(hessians_quad + @@ -2439,19 +2421,18 @@ namespace internal if (integration_flag & EvaluationFlags::gradients) eval0.template gradients<0, false, true>(scratch_data, values_dofs + - size_deg); + n_dofs); else eval0.template gradients<0, false, false>(scratch_data, values_dofs + - size_deg); + n_dofs); // grad yz eval1.template gradients<1, false, false>(hessians_quad + 5 * n_q_points, scratch_data); eval0.template values<0, false, true>(scratch_data, - values_dofs + - size_deg); + values_dofs + n_dofs); break; case 2: @@ -2466,14 +2447,14 @@ namespace internal // grad yy eval0.template values<0, false, false>( - hessians_quad + n_q_points, values_dofs + 2 * size_deg); + hessians_quad + n_q_points, values_dofs + 2 * n_dofs); // grad xy if (integration_flag & EvaluationFlags::gradients) eval0.template gradients<0, false, true>( - hessians_quad + 2 * n_q_points, values_dofs + size_deg); + hessians_quad + 2 * n_q_points, values_dofs + n_dofs); else eval0.template gradients<0, false, false>( - hessians_quad + 2 * n_q_points, values_dofs + size_deg); + hessians_quad + 2 * n_q_points, values_dofs + n_dofs); break; case 1: values_dofs[2] = hessians_quad[0]; @@ -2485,7 +2466,7 @@ namespace internal default: AssertThrow(false, ExcNotImplemented()); } - values_dofs += 3 * size_deg; + values_dofs += 3 * n_dofs; hessians_quad += dim * (dim + 1) / 2 * n_q_points; } } @@ -2553,7 +2534,7 @@ namespace internal private: template - static DEAL_II_ALWAYS_INLINE void + static void interpolate_generic(const unsigned int n_components, const Number * input, Number * output, @@ -2804,8 +2785,7 @@ namespace internal run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, - FEEvaluationBaseData - &eval) + FEEvaluationData & eval) { const auto & data = eval.get_shape_info(); const auto & shape_data = data.data.front(); @@ -2814,8 +2794,8 @@ namespace internal 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 std::size_t n_dofs = data.dofs_per_component_on_cell; + const std::size_t n_q_points = data.n_q_points_faces[face_no]; using Eval = EvaluatorTensorProduct; @@ -2878,6 +2858,8 @@ namespace internal fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : Utilities::pow(shape_data.fe_degree + 1, dim - 1); + // Note: we always keep storage of values, 1st and 2nd derivatives in an + // array, so reserve space for all three here Number *temp = eval.get_scratch_data().begin(); Number *scratch_data = temp + 3 * n_components * dofs_per_face; @@ -2947,8 +2929,7 @@ namespace internal run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, - FEEvaluationBaseData - &eval) + FEEvaluationData & eval) { const auto & data = eval.get_shape_info(); const auto & shape_data = data.data.front(); @@ -2957,8 +2938,8 @@ namespace internal 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 std::size_t n_dofs = data.dofs_per_component_on_cell; + const std::size_t n_q_points = data.n_q_points_faces[face_no]; using Eval = EvaluatorTensorProduct; @@ -3716,11 +3697,11 @@ namespace internal { template static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number2 * src_ptr, - const std::vector> *sm_ptr, - FEEvaluationBaseData &eval) + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number2 * src_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &eval) { Assert(fe_degree > -1, ExcInternalError()); Assert(eval.get_shape_info().element_type <= @@ -3903,11 +3884,11 @@ namespace internal { template static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number2 * dst_ptr, - const std::vector> *sm_ptr, - FEEvaluationBaseData &eval) + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number2 * dst_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &eval) { Assert(fe_degree > -1, ExcInternalError()); Assert(eval.get_shape_info().element_type <= @@ -4070,20 +4051,17 @@ namespace internal /** * This struct implements the action of the inverse mass matrix operation - * using an FEEvaluationBaseData argument. + * using an FEEvaluationData argument. */ template struct CellwiseInverseMassMatrixImplBasic { template static bool - run(const unsigned int n_components, - const FEEvaluationBaseData &fe_eval, - const Number * in_array, - Number * out_array, + run(const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array, typename std::enable_if::type * = nullptr) { constexpr unsigned int dofs_per_component = @@ -4136,13 +4114,10 @@ namespace internal template static bool - run(const unsigned int n_components, - const FEEvaluationBaseData &fe_eval, - const Number * in_array, - Number * out_array, + run(const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array, typename std::enable_if::type * = nullptr) { static_assert(fe_degree == -1, "Only usable for degree -1"); @@ -4194,7 +4169,7 @@ namespace internal /** * This struct implements the action of the inverse mass matrix operation - * using an FEEvaluationBaseData argument. + * using an FEEvaluationData argument. */ template struct CellwiseInverseMassMatrixImplFlexible @@ -4281,20 +4256,17 @@ namespace internal /** * This struct implements the action of the inverse mass matrix operation - * using an FEEvaluationBaseData argument. + * using an FEEvaluationData argument. */ template struct CellwiseInverseMassMatrixImplTransformFromQPoints { template static bool - run(const unsigned int n_desired_components, - const FEEvaluationBaseData &fe_eval, - const Number * in_array, - Number * out_array) + run(const unsigned int n_desired_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array) { static const bool do_inplace = fe_degree > -1 && (fe_degree + 1 == n_q_points_1d); @@ -4308,12 +4280,12 @@ namespace internal fe_eval.get_shape_info().data.front().inverse_shape_values_eo : fe_eval.get_shape_info().data.front().inverse_shape_values; - const unsigned int dofs_per_component = + const std::size_t dofs_per_component = do_inplace ? Utilities::pow(fe_degree + 1, dim) : fe_eval.get_shape_info().dofs_per_component_on_cell; - const unsigned int n_q_points = do_inplace ? - Utilities::pow(fe_degree + 1, dim) : - fe_eval.get_shape_info().n_q_points; + const std::size_t n_q_points = do_inplace ? + Utilities::pow(fe_degree + 1, dim) : + fe_eval.get_shape_info().n_q_points; internal::EvaluatorTensorProduct 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(), @@ -49,7 +46,7 @@ struct SelectEvaluator evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, - FEEvaluationBaseData &eval); + FEEvaluationData & eval); /** * Chooses an appropriate evaluation strategy for the integrate function, i.e. @@ -62,7 +59,7 @@ struct SelectEvaluator integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, - FEEvaluationBaseData &eval, + FEEvaluationData & eval, const bool sum_into_values_array = false); }; @@ -72,10 +69,10 @@ struct SelectEvaluator template inline void SelectEvaluator::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number * values_dofs, - FEEvaluationBaseData &eval) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + FEEvaluationData & eval) { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); @@ -89,11 +86,11 @@ SelectEvaluator::evaluate( template inline void SelectEvaluator::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number * values_dofs, - FEEvaluationBaseData &eval, - const bool sum_into_values_array) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + FEEvaluationData & eval, + const bool sum_into_values_array) { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index e0152471f7..22c3591a66 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -29,29 +29,27 @@ DEAL_II_NAMESPACE_OPEN -template -class FEEvaluationBaseData; +template +class FEEvaluationData; namespace internal { - template + template struct FEEvaluationFactory { static void - evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval); + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + FEEvaluationData & eval); static void - integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval, - const bool sum_into_values_array); + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + FEEvaluationData & eval, + const bool sum_into_values_array); static bool fast_evaluation_supported(const unsigned int given_degree, @@ -60,75 +58,69 @@ namespace internal - template > + template struct FEFaceEvaluationFactory { static void - evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval); + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + FEEvaluationData & eval); static void - integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval); - - static void - gather_evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number * src_ptr, - const std::vector> *sm_ptr, - FEEvaluationBaseData &eval); - - static void - integrate_scatter( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number * dst_ptr, - const std::vector> *sm_ptr, - FEEvaluationBaseData &eval); + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + FEEvaluationData & eval); static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d); }; + template + struct FEFaceEvaluationGatherFactory + { + static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * src_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &eval); + + static void + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * dst_ptr, + const std::vector> *sm_ptr, + FEEvaluationData &eval); + }; + - template > + template struct CellwiseInverseMassFactory { static void - apply(const unsigned int n_components, - const FEEvaluationBaseData - & fe_eval, - const VectorizedArrayType *in_array, - VectorizedArrayType * out_array); + apply(const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array); static void - apply(const unsigned int n_components, - const unsigned int fe_degree, - const AlignedVector &inverse_shape, - const AlignedVector &inverse_coefficients, - const VectorizedArrayType * in_array, - VectorizedArrayType * out_array); + apply(const unsigned int n_components, + const unsigned int fe_degree, + const AlignedVector &inverse_shape, + const AlignedVector &inverse_coefficients, + const Number * in_array, + Number * out_array); static void transform_from_q_points_to_basis( - const unsigned int n_components, - const FEEvaluationBaseData - & fe_eval, - const VectorizedArrayType *in_array, - VectorizedArrayType * out_array); + const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array); }; template + template void - FEEvaluationFactory::evaluate( + FEEvaluationFactory::evaluate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval) + const Number * values_dofs, + FEEvaluationData & eval) { - instantiation_helper_run< - 1, - FEEvaluationImplEvaluateSelector>( + instantiation_helper_run<1, FEEvaluationImplEvaluateSelector>( eval.get_shape_info().data[0].fe_degree, eval.get_shape_info().data[0].n_q_points_1d, n_components, @@ -62,18 +60,16 @@ namespace internal - template + template void - FEEvaluationFactory::integrate( + FEEvaluationFactory::integrate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval, - const bool sum_into_values_array) + Number * values_dofs, + FEEvaluationData & eval, + const bool sum_into_values_array) { - instantiation_helper_run< - 1, - FEEvaluationImplIntegrateSelector>( + instantiation_helper_run<1, FEEvaluationImplIntegrateSelector>( eval.get_shape_info().data[0].fe_degree, eval.get_shape_info().data[0].n_q_points_1d, n_components, @@ -85,11 +81,11 @@ namespace internal - template + template bool - FEEvaluationFactory:: - fast_evaluation_supported(const unsigned int given_degree, - const unsigned int n_q_points_1d) + FEEvaluationFactory::fast_evaluation_supported( + const unsigned int given_degree, + const unsigned int n_q_points_1d) { return instantiation_helper_run<1, FastEvaluationSupported>(given_degree, n_q_points_1d); @@ -97,17 +93,16 @@ namespace internal - template + template void - FEFaceEvaluationFactory::evaluate( + FEFaceEvaluationFactory::evaluate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval) + const Number * values_dofs, + FEEvaluationData & eval) { - instantiation_helper_run< - 1, - FEFaceEvaluationImplEvaluateSelector>( + instantiation_helper_run<1, + FEFaceEvaluationImplEvaluateSelector>( eval.get_shape_info().data[0].fe_degree, eval.get_shape_info().data[0].n_q_points_1d, n_components, @@ -118,17 +113,17 @@ namespace internal - template + template void - FEFaceEvaluationFactory::integrate( + FEFaceEvaluationFactory::integrate( const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType * values_dofs, - FEEvaluationBaseData &eval) + Number * values_dofs, + FEEvaluationData & eval) { instantiation_helper_run< 1, - FEFaceEvaluationImplIntegrateSelector>( + FEFaceEvaluationImplIntegrateSelector>( eval.get_shape_info().data[0].fe_degree, eval.get_shape_info().data[0].n_q_points_1d, n_components, @@ -141,12 +136,12 @@ namespace internal template void - FEFaceEvaluationFactory::gather_evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number * src_ptr, - const std::vector> *sm_ptr, - FEEvaluationBaseData &eval) + FEFaceEvaluationGatherFactory::evaluate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * src_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &eval) { instantiation_helper_run< 1, @@ -166,12 +161,12 @@ namespace internal template void - FEFaceEvaluationFactory::integrate_scatter( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number * dst_ptr, - const std::vector> *sm_ptr, - FEEvaluationBaseData &eval) + FEFaceEvaluationGatherFactory::integrate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * dst_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &eval) { instantiation_helper_run< 1, @@ -189,11 +184,11 @@ namespace internal - template + template bool - FEFaceEvaluationFactory:: - fast_evaluation_supported(const unsigned int given_degree, - const unsigned int n_q_points_1d) + FEFaceEvaluationFactory::fast_evaluation_supported( + const unsigned int given_degree, + const unsigned int n_q_points_1d) { return instantiation_helper_run<1, FastEvaluationSupported>(given_degree, n_q_points_1d); @@ -201,65 +196,59 @@ namespace internal - template + template void - CellwiseInverseMassFactory::apply( - const unsigned int n_components, - const FEEvaluationBaseData - & fe_eval, - const VectorizedArrayType *in_array, - VectorizedArrayType * out_array) + CellwiseInverseMassFactory::apply( + const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array) { const unsigned int fe_degree = fe_eval.get_shape_info().data[0].fe_degree; - instantiation_helper_run< - 1, - CellwiseInverseMassMatrixImplBasic>( + instantiation_helper_run<1, + CellwiseInverseMassMatrixImplBasic>( fe_degree, fe_degree + 1, n_components, fe_eval, in_array, out_array); } - template + template void - CellwiseInverseMassFactory::apply( - const unsigned int n_components, - const unsigned int fe_degree, - const AlignedVector &inverse_shape, - const AlignedVector &inverse_coefficients, - const VectorizedArrayType * in_array, - VectorizedArrayType * out_array) + CellwiseInverseMassFactory::apply( + const unsigned int n_components, + const unsigned int fe_degree, + const AlignedVector &inverse_shape, + const AlignedVector &inverse_coefficients, + const Number * in_array, + Number * out_array) { instantiation_helper_run< 1, - CellwiseInverseMassMatrixImplFlexible>( - fe_degree, - fe_degree + 1, - n_components, - inverse_shape, - inverse_coefficients, - in_array, - out_array); + CellwiseInverseMassMatrixImplFlexible>(fe_degree, + fe_degree + 1, + n_components, + inverse_shape, + inverse_coefficients, + in_array, + out_array); } - template + template void - CellwiseInverseMassFactory:: - transform_from_q_points_to_basis( - const unsigned int n_components, - const FEEvaluationBaseData - & fe_eval, - const VectorizedArrayType *in_array, - VectorizedArrayType * out_array) + CellwiseInverseMassFactory::transform_from_q_points_to_basis( + const unsigned int n_components, + const FEEvaluationData &fe_eval, + const Number * in_array, + Number * out_array) { const unsigned int fe_degree = fe_eval.get_shape_info().data[0].fe_degree; const unsigned int n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d; instantiation_helper_run< 1, - CellwiseInverseMassMatrixImplTransformFromQPoints>( + CellwiseInverseMassMatrixImplTransformFromQPoints>( fe_degree, n_q_points_1d, n_components, fe_eval, in_array, out_array); } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a3192566af..d6f88674ad 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -33,7 +33,7 @@ #include #include #include -#include +#include #include #include #include @@ -59,7 +59,7 @@ DEAL_II_NAMESPACE_OPEN * FEEvaluationBase::distribute_local_to_global() functions, as well as * methods to access values and gradients of finite element functions. It also * inherits the geometry access functions provided by the class - * FEEvaluationBaseData. + * FEEvaluationData. * * This class has five template arguments: * @@ -91,7 +91,7 @@ template class FEEvaluationBase - : public FEEvaluationBaseData + : public FEEvaluationData { public: using number_type = Number; @@ -671,7 +671,7 @@ protected: * possibly within the element if the evaluate/integrate routines are * combined inside user code (e.g. for computing cell matrices). * - * The optional FEEvaluationBaseData object allows several + * The optional FEEvaluationData object allows several * FEEvaluation objects to share the geometry evaluation, i.e., the * underlying mapping and quadrature points do only need to be evaluated * once. This only works if the quadrature formulas are the same. Otherwise, @@ -686,8 +686,7 @@ protected: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other); + const FEEvaluationData *other); /** * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe, @@ -845,8 +844,7 @@ protected: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other); + const FEEvaluationData *other); /** * Copy constructor @@ -1003,8 +1001,7 @@ protected: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other); + const FEEvaluationData *other); /** * Copy constructor @@ -1169,8 +1166,7 @@ protected: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other); + const FEEvaluationData *other); /** * Copy constructor @@ -1351,7 +1347,7 @@ protected: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData<1, Number, is_face, VectorizedArrayType> *other); + const FEEvaluationData<1, VectorizedArrayType, is_face> *other); /** * Copy constructor @@ -2122,10 +2118,9 @@ public: * use the FEEvaluation object in %parallel to the given one because * otherwise the intended sharing may create race conditions. */ - FEEvaluation( - const FiniteElement & fe, - const FEEvaluationBaseData &other, - const unsigned int first_selected_component = 0); + FEEvaluation(const FiniteElement & fe, + const FEEvaluationData &other, + const unsigned int first_selected_component = 0); /** * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe, @@ -2815,7 +2810,7 @@ namespace internal const internal::MatrixFreeFunctions::ShapeInfo *, const internal::MatrixFreeFunctions::DoFInfo *, const internal::MatrixFreeFunctions:: - MappingInfoStorage *, + MappingInfoStorage *, unsigned int, unsigned int> get_shape_info_and_indices( @@ -2882,7 +2877,7 @@ inline FEEvaluationBase( + : FEEvaluationData( internal::get_shape_info_and_indices(data_in, dof_no, first_selected_component, @@ -2946,15 +2941,13 @@ inline FEEvaluationBase & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other) - : FEEvaluationBaseData( + const FEEvaluationData *other) + : FEEvaluationData( other != nullptr && other->mapped_geometry->get_quadrature() == quadrature ? other->mapped_geometry : - std::make_shared< - internal::MatrixFreeFunctions:: - MappingDataOnTheFly>( + std::make_shared>( mapping, quadrature, update_flags), @@ -3000,7 +2993,7 @@ inline FEEvaluationBase &other) - : FEEvaluationBaseData(other) + : FEEvaluationData(other) , scratch_data_array(other.matrix_info == nullptr ? new AlignedVector() : other.matrix_info->acquire_scratch_data()) @@ -3014,12 +3007,12 @@ inline FEEvaluationBasemapped_geometry = std::make_shared< - internal::MatrixFreeFunctions:: - MappingDataOnTheFly>( - other.mapped_geometry->get_fe_values().get_mapping(), - other.mapped_geometry->get_quadrature(), - other.mapped_geometry->get_fe_values().get_update_flags()); + this->mapped_geometry = + std::make_shared>( + other.mapped_geometry->get_fe_values().get_mapping(), + other.mapped_geometry->get_quadrature(), + other.mapped_geometry->get_fe_values().get_update_flags()); this->mapping_data = &this->mapped_geometry->get_data_storage(); this->cell = 0; @@ -3062,8 +3055,7 @@ operator=(const FEEvaluationBaserelease_scratch_data(scratch_data_array); } - this->FEEvaluationBaseData:: - operator=(other); + this->FEEvaluationData::operator=(other); matrix_info = other.matrix_info; @@ -3076,12 +3068,12 @@ operator=(const FEEvaluationBase(); // Create deep copy of mapped geometry for use in parallel - this->mapped_geometry = std::make_shared< - internal::MatrixFreeFunctions:: - MappingDataOnTheFly>( - other.mapped_geometry->get_fe_values().get_mapping(), - other.mapped_geometry->get_quadrature(), - other.mapped_geometry->get_fe_values().get_update_flags()); + this->mapped_geometry = + std::make_shared>( + other.mapped_geometry->get_fe_values().get_mapping(), + other.mapped_geometry->get_quadrature(), + other.mapped_geometry->get_fe_values().get_update_flags()); this->cell = 0; this->mapping_data = &this->mapped_geometry->get_data_storage(); this->jacobian = @@ -5332,8 +5324,7 @@ inline FEEvaluationAccess & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other) + const FEEvaluationData *other) : FEEvaluationBase( mapping, fe, @@ -5432,8 +5423,7 @@ inline FEEvaluationAccess:: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other) + const FEEvaluationData *other) : FEEvaluationBase( mapping, fe, @@ -5759,8 +5749,7 @@ inline FEEvaluationAccess:: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData - *other) + const FEEvaluationData *other) : FEEvaluationBase( mapping, fe, @@ -6170,7 +6159,7 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: const Quadrature<1> & quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBaseData<1, Number, is_face, VectorizedArrayType> *other) + const FEEvaluationData<1, VectorizedArrayType, is_face> *other) : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>( mapping, fe, @@ -6591,10 +6580,9 @@ inline FEEvaluation:: - FEEvaluation( - const FiniteElement & fe, - const FEEvaluationBaseData &other, - const unsigned int first_selected_component) + FEEvaluation(const FiniteElement & fe, + const FEEvaluationData &other, + const unsigned int first_selected_component) : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(), fe, other.mapped_geometry->get_quadrature(), @@ -7098,7 +7086,7 @@ FEEvaluation:: evaluate(n_components, evaluation_flag_actual, values_array, *this); else - internal::FEEvaluationFactory::evaluate( + internal::FEEvaluationFactory::evaluate( n_components, evaluation_flag_actual, const_cast(values_array), @@ -7267,10 +7255,10 @@ namespace internal VectorizedArrayType *vec_values = reinterpret_cast( destination.begin() + - dof_info->dof_indices_contiguous + dof_info.dof_indices_contiguous [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] [cell * VectorizedArrayType::size()] + - dof_info->component_dof_indices_offset + dof_info.component_dof_indices_offset [phi.get_active_fe_index()] [phi.get_first_selected_component()] * VectorizedArrayType::size()); @@ -7470,7 +7458,7 @@ FEEvaluation::integrate( + internal::FEEvaluationFactory::integrate( n_components, integration_flag_actual, values_array, @@ -7925,8 +7913,8 @@ FEFaceEvaluation:: - evaluate(n_components, evaluation_flag_actual, values_array, *this); + internal::FEFaceEvaluationFactory::evaluate( + n_components, evaluation_flag_actual, values_array, *this); # ifdef DEBUG if (evaluation_flag_actual & EvaluationFlags::values) @@ -8062,8 +8050,8 @@ FEFaceEvaluation:: - integrate(n_components, integration_flag_actual, values_array, *this); + internal::FEFaceEvaluationFactory::integrate( + n_components, integration_flag_actual, values_array, *this); } @@ -8149,12 +8137,15 @@ FEFaceEvaluation:: - gather_evaluate(n_components, - evaluation_flag, - internal::get_beginning(input_vector), - shared_vector_data, - *this); + internal::FEFaceEvaluationGatherFactory< + dim, + Number, + VectorizedArrayType>::evaluate(n_components, + evaluation_flag, + internal::get_beginning( + input_vector), + shared_vector_data, + *this); } else { @@ -8254,12 +8245,15 @@ FEFaceEvaluation:: - integrate_scatter(n_components, - integration_flag, - internal::get_beginning(destination), - shared_vector_data, - *this); + internal::FEFaceEvaluationGatherFactory< + dim, + Number, + VectorizedArrayType>::integrate(n_components, + integration_flag, + internal::get_beginning( + destination), + shared_vector_data, + *this); } else { @@ -8337,7 +8331,7 @@ FEEvaluation:: + internal::FEEvaluationFactory:: fast_evaluation_supported(given_degree, give_n_q_points_1d) : true; } @@ -8361,7 +8355,7 @@ FEFaceEvaluation:: + internal::FEFaceEvaluationFactory:: fast_evaluation_supported(given_degree, give_n_q_points_1d) : true; } diff --git a/include/deal.II/matrix_free/fe_evaluation_base_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h similarity index 76% rename from include/deal.II/matrix_free/fe_evaluation_base_data.h rename to include/deal.II/matrix_free/fe_evaluation_data.h index 447297c1c8..020dd83dc1 100644 --- a/include/deal.II/matrix_free/fe_evaluation_base_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -14,8 +14,8 @@ // --------------------------------------------------------------------- -#ifndef dealii_matrix_free_fe_evaluation_base_data_h -#define dealii_matrix_free_fe_evaluation_base_data_h +#ifndef dealii_matrix_free_fe_evaluation_data_h +#define dealii_matrix_free_fe_evaluation_data_h #include @@ -74,7 +74,7 @@ namespace internal { namespace MatrixFreeFunctions { - template + template class MappingDataOnTheFly; } } // namespace internal @@ -85,7 +85,7 @@ namespace internal * mapping-related information independent of the degrees of freedom and * finite element in use. This class provides access functionality for user * code. The main usage is through the class FEEvaluation. However, there is - * functionality to construct an object of type FEEvaluationBaseData given + * functionality to construct an object of type FEEvaluationData given * suitable data pointers to internal data, which allows usage in some * scenarios where no full MatrixFree object is available. * @@ -108,52 +108,49 @@ namespace internal * * @ingroup matrixfree */ -template -class FEEvaluationBaseData +template +class FEEvaluationData { - static_assert( - std::is_same::value, - "Type of Number and of VectorizedArrayType do not match."); - - using ShapeInfoType = - internal::MatrixFreeFunctions::ShapeInfo; - using MappingInfoStorageType = - internal::MatrixFreeFunctions::MappingInfoStorage<(is_face ? dim - 1 : dim), - dim, - Number, - VectorizedArrayType>; + using ShapeInfoType = internal::MatrixFreeFunctions::ShapeInfo; + using MappingInfoStorageType = internal::MatrixFreeFunctions:: + MappingInfoStorage<(is_face ? dim - 1 : dim), dim, Number>; using DoFInfo = internal::MatrixFreeFunctions::DoFInfo; public: static constexpr unsigned int dimension = dim; - using NumberType = VectorizedArrayType; + using NumberType = Number; + + using ScalarNumber = + typename internal::VectorizedArrayTrait::value_type; + + static constexpr unsigned int n_lanes = sizeof(Number) / sizeof(ScalarNumber); /** * 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 bool is_interior_face = true, - const unsigned int quad_no = 0, - const unsigned int face_type = 0, - const unsigned int first_selected_component = 0); + FEEvaluationData(const std::tuple &info, + const bool is_interior_face = true, + const unsigned int quad_no = 0, + const unsigned int face_type = 0, + const unsigned int first_selected_component = 0); /** * Copy constructor. */ - FEEvaluationBaseData(const FEEvaluationBaseData &other) = default; + FEEvaluationData(const FEEvaluationData &other) = default; /** * Copy assignment operator. */ - FEEvaluationBaseData & - operator=(const FEEvaluationBaseData &other); + FEEvaluationData & + operator=(const FEEvaluationData &other); /** * Sets the pointers for values, gradients, hessians to the central @@ -161,16 +158,16 @@ public: * components as provided by one of the derived classes. */ void - set_data_pointers(AlignedVector *scratch_data, - const unsigned int n_components); + 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); + reinit_face( + const internal::MatrixFreeFunctions::FaceToCellTopology &face); /** * @name 1: Access to geometry data at quadrature points @@ -183,14 +180,14 @@ public: * do not assume it to be stable over those calls. The maximum size you can * write into is 3*dofs_per_cell+2*n_q_points. */ - ArrayView + ArrayView get_scratch_data() const; /** * Return the determinant of the Jacobian from the unit to the real cell * times the quadrature weight. */ - VectorizedArrayType + Number JxW(const unsigned int q_point) const; /** @@ -204,7 +201,7 @@ public: * gradients to gradients on the real cell by a multiplication $J^{-\mathrm * T} \hat{\nabla} u_h$. */ - Tensor<2, dim, VectorizedArrayType> + Tensor<2, dim, Number> inverse_jacobian(const unsigned int q_point) const; /** @@ -219,7 +216,7 @@ public: * * @note Only implemented in case `is_face == true`. */ - Tensor<1, dim, VectorizedArrayType> + Tensor<1, dim, Number> get_normal_vector(const unsigned int q_point) const; //@} @@ -236,7 +233,7 @@ public: * used in this class. In general, it is safer to use the get_dof_value() * function instead. */ - const VectorizedArrayType * + const Number * begin_dof_values() const; /** @@ -247,7 +244,7 @@ public: * structures used in this class. In general, it is safer to use the * get_dof_value() function instead. */ - VectorizedArrayType * + Number * begin_dof_values(); /** @@ -260,7 +257,7 @@ public: * applied manually. In general, it is safer to use the get_value() function * instead, which does all the transformation internally. */ - const VectorizedArrayType * + const Number * begin_values() const; /** @@ -273,7 +270,7 @@ public: * applied manually. In general, it is safer to use the get_value() function * instead, which does all the transformation internally. */ - VectorizedArrayType * + Number * begin_values(); /** @@ -287,7 +284,7 @@ public: * manually. In general, it is safer to use the get_gradient() function * instead, which does all the transformation internally. */ - const VectorizedArrayType * + const Number * begin_gradients() const; /** @@ -301,7 +298,7 @@ public: * manually. In general, it is safer to use the get_gradient() function * instead, which does all the transformation internally. */ - VectorizedArrayType * + Number * begin_gradients(); /** @@ -316,7 +313,7 @@ public: * general, it is safer to use the get_laplacian() or get_hessian() * functions instead, which does all the transformation internally. */ - const VectorizedArrayType * + const Number * begin_hessians() const; /** @@ -331,7 +328,7 @@ public: * general, it is safer to use the get_laplacian() or get_hessian() * functions instead, which does all the transformation internally. */ - VectorizedArrayType * + Number * begin_hessians(); //@} @@ -484,14 +481,14 @@ public: * Return the id of the cells this FEEvaluation or FEFaceEvaluation is * associated with. */ - const std::array & + const std::array & get_cell_ids() const; /** * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is * associated with. */ - const std::array & + const std::array & get_cell_or_face_ids() const; /** @@ -506,7 +503,7 @@ public: * is not stable between releases of deal.II, and is hence mostly for * internal use. */ - const std::array & + const std::array & get_all_face_numbers() const; /** @@ -517,7 +514,7 @@ public: * @note Only available for `dof_access_index == dof_access_cell` and * `is_interior_face == false`. */ - const std::array & + const std::array & get_all_face_orientations() const; //@} @@ -527,11 +524,11 @@ protected: * Constructor that comes with reduced functionality and works similarly as * FEValues. */ - FEEvaluationBaseData( + FEEvaluationData( const std::shared_ptr< - internal::MatrixFreeFunctions:: - MappingDataOnTheFly> &mapping_data, - const unsigned int n_fe_components, + internal::MatrixFreeFunctions::MappingDataOnTheFly> + & mapping_data, + const unsigned int n_fe_components, const unsigned int first_selected_component); /** @@ -603,7 +600,7 @@ protected: * A pointer to the Jacobian information of the present cell. Only set to a * useful value if on a non-Cartesian cell. */ - const Tensor<2, dim, VectorizedArrayType> *jacobian; + const Tensor<2, dim, Number> *jacobian; /** * A pointer to the Jacobian determinant of the present cell. If on a @@ -611,22 +608,22 @@ protected: * Jacobian determinant, otherwise the Jacobian determinant times the * quadrature weight. */ - const VectorizedArrayType *J_value; + const Number *J_value; /** * A pointer to the normal vectors at faces. */ - const Tensor<1, dim, VectorizedArrayType> *normal_vectors; + const Tensor<1, dim, Number> *normal_vectors; /** * A pointer to the normal vectors times the jacobian at faces. */ - const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian; + const Tensor<1, dim, Number> *normal_x_jacobian; /** * A pointer to the quadrature weights of the underlying quadrature formula. */ - const Number *quadrature_weights; + const ScalarNumber *quadrature_weights; /** * This is the user-visible part of FEEvaluationBase::scratch_data_array, @@ -634,7 +631,7 @@ protected: * during the allocation of the internal data structures in * FEEvaluationBase. */ - mutable ArrayView scratch_data; + mutable ArrayView scratch_data; /** * This field stores the values for local degrees of freedom (e.g. after @@ -646,7 +643,7 @@ protected: * @p scratch_data_array. Due to its access as a thread local memory, the * memory can get reused between different calls. */ - VectorizedArrayType *values_dofs; + Number *values_dofs; /** * This field stores the values of the finite element function on quadrature @@ -657,7 +654,7 @@ protected: * @p scratch_data_array. Due to its access as a thread local memory, the * memory can get reused between different calls. */ - VectorizedArrayType *values_quad; + Number *values_quad; /** * This field stores the gradients of the finite element function on @@ -670,7 +667,7 @@ protected: * @p scratch_data_array. Due to its access as a thread local memory, the * memory can get reused between different calls. */ - VectorizedArrayType *gradients_quad; + Number *gradients_quad; /** * This field stores the gradients of the finite element function on @@ -683,7 +680,7 @@ protected: * @p scratch_data_array. Due to its access as a thread local memory, the * memory can get reused between different calls. */ - VectorizedArrayType *gradients_from_hessians_quad; + Number *gradients_from_hessians_quad; /** * This field stores the Hessians of the finite element function on @@ -694,7 +691,7 @@ protected: * @p scratch_data_array. Due to its access as a thread local memory, the * memory can get reused between different calls. */ - VectorizedArrayType *hessians_quad; + Number *hessians_quad; /** * Debug information to track whether dof values have been initialized @@ -798,13 +795,13 @@ protected: * Stores the (non-vectorized) id of the cells or faces this object is * initialized to. Relevant for ECL. */ - std::array cell_ids; + std::array cell_ids; /** * Stores the (non-vectorized) id of the cells or faces this object is * initialized to. Relevant for ECL. */ - std::array cell_or_face_ids; + std::array cell_or_face_ids; /** * Stores the (non-vectorized) number of faces within cells in case of ECL @@ -814,7 +811,7 @@ protected: * @note Only available for `dof_access_index == dof_access_cell` and * `is_interior_face == false`. */ - std::array all_face_numbers; + std::array all_face_numbers; /** * Store the orientation of the neighbor's faces with respect to the current @@ -824,15 +821,15 @@ protected: * @note Only available for `dof_access_index == dof_access_cell` and * `is_interior_face == false`. */ - std::array all_face_orientations; + std::array all_face_orientations; /** * Geometry data that can be generated FEValues on the fly with the * respective constructor, as an alternative to the entry point with * MatrixFree. */ - std::shared_ptr> + std::shared_ptr< + internal::MatrixFreeFunctions::MappingDataOnTheFly> mapped_geometry; // Make FEEvaluation and FEEvaluationBase objects friends for access to @@ -850,17 +847,17 @@ protected: #ifndef DOXYGEN -template -inline FEEvaluationBaseData:: - FEEvaluationBaseData(const std::tuple &shape_dof_info, - const bool is_interior_face, - const unsigned int quad_no, - const unsigned int face_type, - const unsigned int first_selected_component) +template +inline FEEvaluationData::FEEvaluationData( + const std::tuple &shape_dof_info, + const bool is_interior_face, + const unsigned int quad_no, + const unsigned int face_type, + const unsigned int first_selected_component) : data(std::get<0>(shape_dof_info)) , dof_info(std::get<1>(shape_dof_info)) , mapping_data(std::get<2>(shape_dof_info)) @@ -903,14 +900,13 @@ inline FEEvaluationBaseData:: -template -inline FEEvaluationBaseData:: - FEEvaluationBaseData( - const std::shared_ptr< - internal::MatrixFreeFunctions:: - MappingDataOnTheFly> &mapped_geometry, - const unsigned int n_fe_components, - const unsigned int first_selected_component) +template +inline FEEvaluationData::FEEvaluationData( + const std::shared_ptr< + internal::MatrixFreeFunctions::MappingDataOnTheFly> + & mapped_geometry, + const unsigned int n_fe_components, + const unsigned int first_selected_component) : data(nullptr) , dof_info(nullptr) , mapping_data(nullptr) @@ -940,10 +936,9 @@ inline FEEvaluationBaseData:: -template -inline FEEvaluationBaseData & -FEEvaluationBaseData::operator=( - const FEEvaluationBaseData &other) +template +inline FEEvaluationData & +FEEvaluationData::operator=(const FEEvaluationData &other) { AssertDimension(quad_no, other.quad_no); AssertDimension(n_fe_components, other.n_fe_components); @@ -987,11 +982,11 @@ FEEvaluationBaseData::operator=( -template +template inline void -FEEvaluationBaseData:: - set_data_pointers(AlignedVector *scratch_data_array, - const unsigned int n_components) +FEEvaluationData::set_data_pointers( + AlignedVector *scratch_data_array, + const unsigned int n_components) { Assert(scratch_data_array != nullptr, ExcInternalError()); @@ -1028,11 +1023,10 @@ FEEvaluationBaseData:: -template +template inline void -FEEvaluationBaseData::reinit_face( - const internal::MatrixFreeFunctions::FaceToCellTopology< - VectorizedArrayType::size()> &face) +FEEvaluationData::reinit_face( + const internal::MatrixFreeFunctions::FaceToCellTopology &face) { Assert(is_face == true, ExcMessage("Faces can only be set if the is_face template parameter " @@ -1059,10 +1053,10 @@ FEEvaluationBaseData::reinit_face( -template -inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> -FEEvaluationBaseData:: - get_normal_vector(const unsigned int q_point) const +template +inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number> +FEEvaluationData::get_normal_vector( + const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); Assert(normal_vectors != nullptr, @@ -1076,10 +1070,9 @@ FEEvaluationBaseData:: -template -inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -FEEvaluationBaseData::JxW( - const unsigned int q_point) const +template +inline DEAL_II_ALWAYS_INLINE Number +FEEvaluationData::JxW(const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); Assert(J_value != nullptr, @@ -1096,10 +1089,10 @@ FEEvaluationBaseData::JxW( -template -inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> -FEEvaluationBaseData:: - inverse_jacobian(const unsigned int q_point) const +template +inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number> +FEEvaluationData::inverse_jacobian( + const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); Assert(jacobian != nullptr, @@ -1113,20 +1106,18 @@ FEEvaluationBaseData:: -template -inline const VectorizedArrayType * -FEEvaluationBaseData:: - begin_dof_values() const +template +inline const Number * +FEEvaluationData::begin_dof_values() const { return values_dofs; } -template -inline VectorizedArrayType * -FEEvaluationBaseData:: - begin_dof_values() +template +inline Number * +FEEvaluationData::begin_dof_values() { # ifdef DEBUG dof_values_initialized = true; @@ -1136,10 +1127,9 @@ FEEvaluationBaseData:: -template -inline const VectorizedArrayType * -FEEvaluationBaseData::begin_values() - const +template +inline const Number * +FEEvaluationData::begin_values() const { # ifdef DEBUG Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized()); @@ -1149,9 +1139,9 @@ FEEvaluationBaseData::begin_values() -template -inline VectorizedArrayType * -FEEvaluationBaseData::begin_values() +template +inline Number * +FEEvaluationData::begin_values() { # ifdef DEBUG values_quad_initialized = true; @@ -1162,10 +1152,9 @@ FEEvaluationBaseData::begin_values() -template -inline const VectorizedArrayType * -FEEvaluationBaseData:: - begin_gradients() const +template +inline const Number * +FEEvaluationData::begin_gradients() const { # ifdef DEBUG Assert(gradients_quad_initialized || gradients_quad_submitted, @@ -1176,10 +1165,9 @@ FEEvaluationBaseData:: -template -inline VectorizedArrayType * -FEEvaluationBaseData:: - begin_gradients() +template +inline Number * +FEEvaluationData::begin_gradients() { # ifdef DEBUG gradients_quad_submitted = true; @@ -1190,10 +1178,9 @@ FEEvaluationBaseData:: -template -inline const VectorizedArrayType * -FEEvaluationBaseData:: - begin_hessians() const +template +inline const Number * +FEEvaluationData::begin_hessians() const { # ifdef DEBUG Assert(hessians_quad_initialized, ExcNotInitialized()); @@ -1203,10 +1190,9 @@ FEEvaluationBaseData:: -template -inline VectorizedArrayType * -FEEvaluationBaseData:: - begin_hessians() +template +inline Number * +FEEvaluationData::begin_hessians() { # ifdef DEBUG hessians_quad_initialized = true; @@ -1216,10 +1202,9 @@ FEEvaluationBaseData:: -template +template inline unsigned int -FEEvaluationBaseData:: - get_mapping_data_index_offset() const +FEEvaluationData::get_mapping_data_index_offset() const { Assert(mapping_data != nullptr, ExcInternalError()); @@ -1234,10 +1219,9 @@ FEEvaluationBaseData:: -template +template inline internal::MatrixFreeFunctions::GeometryType -FEEvaluationBaseData::get_cell_type() - const +FEEvaluationData::get_cell_type() const { Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_type; @@ -1245,10 +1229,9 @@ FEEvaluationBaseData::get_cell_type() -template -inline const internal::MatrixFreeFunctions::ShapeInfo & -FEEvaluationBaseData:: - get_shape_info() const +template +inline const internal::MatrixFreeFunctions::ShapeInfo & +FEEvaluationData::get_shape_info() const { Assert(data != nullptr, ExcInternalError()); return *data; @@ -1256,10 +1239,9 @@ FEEvaluationBaseData:: -template +template inline const internal::MatrixFreeFunctions::DoFInfo & -FEEvaluationBaseData::get_dof_info() - const +FEEvaluationData::get_dof_info() const { Assert(dof_info != nullptr, ExcMessage( @@ -1269,30 +1251,27 @@ FEEvaluationBaseData::get_dof_info() -template +template inline const std::vector & -FEEvaluationBaseData:: - get_internal_dof_numbering() const +FEEvaluationData::get_internal_dof_numbering() const { return data->lexicographic_numbering; } -template +template inline unsigned int -FEEvaluationBaseData:: - get_quadrature_index() const +FEEvaluationData::get_quadrature_index() const { return quad_no; } -template +template inline unsigned int -FEEvaluationBaseData:: - get_current_cell_index() const +FEEvaluationData::get_current_cell_index() const { if (is_face && dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) @@ -1303,80 +1282,72 @@ FEEvaluationBaseData:: -template +template inline unsigned int -FEEvaluationBaseData:: - get_active_fe_index() const +FEEvaluationData::get_active_fe_index() const { return active_fe_index; } -template +template inline unsigned int -FEEvaluationBaseData:: - get_active_quadrature_index() const +FEEvaluationData::get_active_quadrature_index() const { return active_quad_index; } -template +template inline unsigned int -FEEvaluationBaseData:: - get_first_selected_component() const +FEEvaluationData::get_first_selected_component() const { return first_selected_component; } -template -inline ArrayView -FEEvaluationBaseData:: - get_scratch_data() const +template +inline ArrayView +FEEvaluationData::get_scratch_data() const { return scratch_data; } -template +template inline unsigned int -FEEvaluationBaseData::get_face_no() - const +FEEvaluationData::get_face_no() const { return face_no; } -template +template inline unsigned int -FEEvaluationBaseData:: - get_subface_index() const +FEEvaluationData::get_subface_index() const { return subface_index; } -template +template inline unsigned int -FEEvaluationBaseData:: - get_face_orientation() const +FEEvaluationData::get_face_orientation() const { return face_orientation; } -template +template inline const Table<2, unsigned int> & -FEEvaluationBaseData:: - get_orientation_map() const +FEEvaluationData::get_orientation_map() const { Assert(descriptor != nullptr, ExcNotInitialized()); return descriptor->face_orientations; @@ -1384,30 +1355,28 @@ FEEvaluationBaseData:: -template +template inline internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex -FEEvaluationBaseData:: - get_dof_access_index() const +FEEvaluationData::get_dof_access_index() const { return dof_access_index; } -template +template inline bool -FEEvaluationBaseData:: - get_is_interior_face() const +FEEvaluationData::get_is_interior_face() const { return is_interior_face; } -template -inline const std::array & -FEEvaluationBaseData::get_cell_ids() - const +template +inline const std::array::n_lanes> & +FEEvaluationData::get_cell_ids() const { Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_ids; @@ -1415,10 +1384,10 @@ FEEvaluationBaseData::get_cell_ids() -template -inline const std::array & -FEEvaluationBaseData:: - get_cell_or_face_ids() const +template +inline const std::array::n_lanes> & +FEEvaluationData::get_cell_or_face_ids() const { Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); if (!is_face || dof_access_index == @@ -1430,10 +1399,10 @@ FEEvaluationBaseData:: -template -inline const std::array & -FEEvaluationBaseData:: - get_all_face_numbers() const +template +inline const std::array::n_lanes> & +FEEvaluationData::get_all_face_numbers() const { Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); Assert(is_face && @@ -1448,10 +1417,10 @@ FEEvaluationBaseData:: -template -inline const std::array & -FEEvaluationBaseData:: - get_all_face_orientations() const +template +inline const std::array::n_lanes> & +FEEvaluationData::get_all_face_orientations() const { Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); Assert(is_face && diff --git a/include/deal.II/matrix_free/mapping_data_on_the_fly.h b/include/deal.II/matrix_free/mapping_data_on_the_fly.h index b1c7390090..822a4c66f2 100644 --- a/include/deal.II/matrix_free/mapping_data_on_the_fly.h +++ b/include/deal.II/matrix_free/mapping_data_on_the_fly.h @@ -45,10 +45,10 @@ namespace internal * deal.II information in a form that FEEvaluation and friends can use for * vectorized access. Since no vectorization over cells is available with * the DoFHandler/Triangulation cell iterators, the interface to - * FEEvaluation's vectorization model is to use @p - * VectorizedArray::n_array_element copies of the same element. This - * interface is thus primarily useful for evaluating several operators on - * the same cell, e.g., when assembling cell matrices. + * FEEvaluation's vectorization model is to use VectorizedArray::size() + * copies of the same element. This interface is thus primarily useful for + * evaluating several operators on the same cell, e.g., when assembling + * cell matrices. * * As opposed to the Mapping classes in deal.II, this class does not * actually provide a boundary description that can be used to evaluate @@ -56,13 +56,9 @@ namespace internal * given deal.II mapping (as passed to the constructor of this class) in a * form accessible to FEEvaluation. */ - template + template class MappingDataOnTheFly { - static_assert( - std::is_same::value, - "Type of Number and of VectorizedArrayType do not match."); - public: /** * Constructor, similar to FEValues. Since this class only evaluates the @@ -116,7 +112,7 @@ namespace internal * MappingInfo. This ensures compatibility with the precomputed data * fields in the MappingInfo class. */ - const MappingInfoStorage & + const MappingInfoStorage & get_data_storage() const; /** @@ -152,24 +148,22 @@ namespace internal * The storage part created for a single cell and held in analogy to * MappingInfo. */ - MappingInfoStorage - mapping_info_storage; + MappingInfoStorage mapping_info_storage; }; /*-------------------------- Inline functions ---------------------------*/ - template - inline MappingDataOnTheFly:: - MappingDataOnTheFly(const Mapping & mapping, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags) - : fe_values( - mapping, - fe_dummy, - Quadrature(quadrature), - MappingInfo::compute_update_flags( - update_flags)) + template + inline MappingDataOnTheFly::MappingDataOnTheFly( + const Mapping & mapping, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags) + : fe_values(mapping, + fe_dummy, + Quadrature(quadrature), + MappingInfoStorage::compute_update_flags( + update_flags)) , quadrature_1d(quadrature) { mapping_info_storage.descriptor.resize(1); @@ -196,10 +190,10 @@ namespace internal - template - inline MappingDataOnTheFly:: - MappingDataOnTheFly(const Quadrature<1> &quadrature, - const UpdateFlags update_flags) + template + inline MappingDataOnTheFly::MappingDataOnTheFly( + const Quadrature<1> &quadrature, + const UpdateFlags update_flags) : MappingDataOnTheFly(::dealii::StaticMappingQ1::mapping, quadrature, update_flags) @@ -207,9 +201,9 @@ namespace internal - template + template inline void - MappingDataOnTheFly::reinit( + MappingDataOnTheFly::reinit( typename dealii::Triangulation::cell_iterator cell) { if (present_cell == cell) @@ -246,10 +240,9 @@ namespace internal - template + template inline bool - MappingDataOnTheFly::is_initialized() - const + MappingDataOnTheFly::is_initialized() const { return present_cell != typename dealii::Triangulation::cell_iterator(); @@ -257,38 +250,36 @@ namespace internal - template + template inline typename dealii::Triangulation::cell_iterator - MappingDataOnTheFly::get_cell() const + MappingDataOnTheFly::get_cell() const { return fe_values.get_cell(); } - template + template inline const dealii::FEValues & - MappingDataOnTheFly::get_fe_values() const + MappingDataOnTheFly::get_fe_values() const { return fe_values; } - template - inline const MappingInfoStorage & - MappingDataOnTheFly::get_data_storage() - const + template + inline const MappingInfoStorage & + MappingDataOnTheFly::get_data_storage() const { return mapping_info_storage; } - template + template inline const Quadrature<1> & - MappingDataOnTheFly::get_quadrature() - const + MappingDataOnTheFly::get_quadrature() const { return quadrature_1d; } diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index 4fc62e48be..7087e5fda3 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -160,20 +160,19 @@ namespace internal /** * The data cache for the cells. */ - std::vector> - cell_data; + std::vector> cell_data; /** * The data cache for the faces. */ - std::vector> + std::vector> face_data; /** * The data cache for the face-associated-with-cell topology, following * the @p cell_type variable for the cell types. */ - std::vector> + std::vector> face_data_by_cells; /** @@ -253,16 +252,6 @@ namespace internal const dealii::Triangulation & tria, const std::vector> &cells, const dealii::hp::MappingCollection & mapping); - - /** - * Helper function to determine which update flags must be set in the - * internal functions to initialize all data as requested by the user. - */ - static UpdateFlags - compute_update_flags( - const UpdateFlags update_flags, - const std::vector> &quad = - std::vector>()); }; @@ -277,7 +266,7 @@ namespace internal template struct MappingInfoCellsOrFaces { - static const MappingInfoStorage & + static const MappingInfoStorage & get(const MappingInfo &mapping_info, const unsigned int quad_no) { @@ -289,10 +278,9 @@ namespace internal template struct MappingInfoCellsOrFaces { - static const MappingInfoStorage - & - get(const MappingInfo &mapping_info, - const unsigned int quad_no) + static const MappingInfoStorage & + get(const MappingInfo &mapping_info, + const unsigned int quad_no) { AssertIndexRange(quad_no, mapping_info.face_data.size()); return mapping_info.face_data[quad_no]; diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index a434fe9b01..8371df6190 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include #include @@ -56,53 +56,6 @@ namespace internal - template - UpdateFlags - MappingInfo::compute_update_flags( - const UpdateFlags update_flags, - const std::vector> &quad) - { - // this class is build around the evaluation of jacobians, so compute - // them in any case. The Jacobians will be inverted manually. Since we - // always do support integration, we also include the JxW values - UpdateFlags new_flags = update_jacobians | update_JxW_values; - - // for Hessian information, need inverse Jacobians and the derivative of - // Jacobians (these two together will give use the gradients of the - // inverse Jacobians, which is what we need) - if (update_flags & update_hessians || - update_flags & update_jacobian_grads) - new_flags |= update_jacobian_grads; - - if (update_flags & update_quadrature_points) - new_flags |= update_quadrature_points; - - // there is one more thing: if we have a quadrature formula with only - // one quadrature point on the first component, but more points on later - // components, we need to have Jacobian gradients anyway in order to - // determine whether the Jacobian is constant throughout a cell - if (quad.empty() == false) - { - bool formula_with_one_point = false; - for (unsigned int i = 0; i < quad[0].size(); ++i) - if (quad[0][i].size() == 1) - { - formula_with_one_point = true; - break; - } - if (formula_with_one_point == true) - for (unsigned int comp = 1; comp < quad.size(); ++comp) - for (unsigned int i = 0; i < quad[comp].size(); ++i) - if (quad[comp][i].size() > 1) - { - new_flags |= update_jacobian_grads; - } - } - return new_flags; - } - - - template void MappingInfo::initialize( @@ -130,7 +83,9 @@ namespace internal // dummy FE that is used to set up an FEValues object. Do not need the // actual finite element because we will only evaluate quantities for // the mapping that are independent of the FE - this->update_flags_cells = compute_update_flags(update_flags_cells, quad); + this->update_flags_cells = + MappingInfoStorage::compute_update_flags( + update_flags_cells, quad); this->update_flags_boundary_faces = ((update_flags_inner_faces | update_flags_boundary_faces) & @@ -646,9 +601,9 @@ namespace internal const std::vector & active_fe_index, const dealii::hp::MappingCollection & mapping, MappingInfo &mapping_info, - std::pair>, - CompressedCellData> &data) + std::pair< + std::vector>, + CompressedCellData> &data) { FE_Nothing dummy_fe; @@ -923,14 +878,14 @@ namespace internal typename Number, typename VectorizedArrayType> void - copy_data(const unsigned int first_cell, - const std::array &data_shift, - const std::vector & indices_compressed, - const std::vector & cell_type, - MappingInfoStorage - &data_cells_local, - MappingInfoStorage - &data_cells) + copy_data( + const unsigned int first_cell, + const std::array &data_shift, + const std::vector & indices_compressed, + const std::vector & cell_type, + MappingInfoStorage + &data_cells_local, + MappingInfoStorage &data_cells) { // Copy the index offsets and shift by the appropriate value for (unsigned int lcell = 0; @@ -1166,7 +1121,7 @@ namespace internal const UpdateFlags update_flags_cells, const AlignedVector & plain_quadrature_points, const ShapeInfo &shape_info, - MappingInfoStorage &my_data) + MappingInfoStorage &my_data) { constexpr unsigned int n_lanes = VectorizedArrayType::size(); constexpr unsigned int n_lanes_d = VectorizedDouble::size(); @@ -1176,9 +1131,9 @@ namespace internal shape_info.dofs_per_component_on_cell; constexpr unsigned int hess_dim = dim * (dim + 1) / 2; - MappingInfoStorage temp_data; + MappingInfoStorage temp_data; temp_data.copy_descriptor(my_data); - FEEvaluationBaseData eval( + FEEvaluationData eval( {&shape_info, nullptr, &temp_data, 0, 0}); AlignedVector evaluation_data; @@ -1198,7 +1153,7 @@ namespace internal start_indices, eval.begin_dof_values()); - FEEvaluationFactory::evaluate( + FEEvaluationFactory::evaluate( dim, EvaluationFlags::values | EvaluationFlags::gradients | (update_flags_cells & update_jacobian_grads ? @@ -1339,7 +1294,7 @@ namespace internal MultithreadInfo::n_threads()); std::vector>, + std::vector>, ExtractCellHelper:: CompressedCellData>> data_cells_local; @@ -1353,8 +1308,7 @@ namespace internal while (cell_range.first < n_cell_batches) { data_cells_local.push_back(std::make_pair( - std::vector< - MappingInfoStorage>( + std::vector>( cell_data.size()), ExtractCellHelper:: CompressedCellData( @@ -1569,8 +1523,7 @@ namespace internal const dealii::hp::MappingCollection & mapping_in, MappingInfo &mapping_info, std::pair< - std::vector< - MappingInfoStorage>, + std::vector>, CompressedFaceData> &data) { std::vector>>> dummy_fe( @@ -2029,9 +1982,8 @@ namespace internal const unsigned int last_face, const std::vector &face_type, const std::vector> - &faces, - MappingInfoStorage - &data_faces) + & faces, + MappingInfoStorage &data_faces) { for (unsigned int face = first_face; face < last_face; ++face) { @@ -2078,7 +2030,7 @@ namespace internal const UpdateFlags update_flags_faces, const AlignedVector & plain_quadrature_points, const ShapeInfo &shape_info, - MappingInfoStorage &my_data) + MappingInfoStorage &my_data) { constexpr unsigned int n_lanes = VectorizedArrayType::size(); constexpr unsigned int n_lanes_d = VectorizedDouble::size(); @@ -2088,11 +2040,11 @@ namespace internal shape_info.dofs_per_component_on_cell; constexpr unsigned int hess_dim = dim * (dim + 1) / 2; - MappingInfoStorage temp_data; + MappingInfoStorage temp_data; temp_data.copy_descriptor(my_data); - FEEvaluationBaseData eval_int( + FEEvaluationData eval_int( {&shape_info, nullptr, &temp_data, 0, 0}, true); - FEEvaluationBaseData eval_ext( + FEEvaluationData eval_ext( {&shape_info, nullptr, &temp_data, 0, 0}, false); AlignedVector evaluation_data, evaluation_data_ext; @@ -2129,7 +2081,7 @@ namespace internal // now let the matrix-free evaluators provide us with the // data on faces - FEFaceEvaluationFactory::evaluate( + FEFaceEvaluationFactory::evaluate( dim, EvaluationFlags::values | EvaluationFlags::gradients | (update_flags_faces & update_jacobian_grads ? @@ -2153,12 +2105,11 @@ namespace internal const unsigned int offset = my_data.data_index_offsets[face]; const auto compute_jacobian_grad = - [&](unsigned int face_no, - int is_exterior, - unsigned int q, - Tensor<2, dim, VectorizedDouble> inv_jac, - FEEvaluationBaseData - &eval) { + [&](unsigned int face_no, + int is_exterior, + unsigned int q, + Tensor<2, dim, VectorizedDouble> inv_jac, + FEEvaluationData &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) @@ -2299,15 +2250,14 @@ namespace internal start_indices, eval_ext.begin_dof_values()); - FEFaceEvaluationFactory:: - evaluate(dim, - EvaluationFlags::values | - EvaluationFlags::gradients | - (update_flags_faces & update_jacobian_grads ? - EvaluationFlags::hessians : - EvaluationFlags::nothing), - eval_ext.begin_dof_values(), - eval_ext); + FEFaceEvaluationFactory::evaluate( + dim, + 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) { @@ -2377,8 +2327,7 @@ namespace internal MultithreadInfo::n_threads()); std::vector>, + std::vector>, ExtractFaceHelper:: CompressedFaceData>> data_faces_local; @@ -2393,7 +2342,7 @@ namespace internal { data_faces_local.push_back(std::make_pair( std::vector< - MappingInfoStorage>( + MappingInfoStorage>( face_data.size()), ExtractFaceHelper:: CompressedFaceData( @@ -2678,7 +2627,7 @@ namespace internal // points, filling up all SIMD lanes as appropriate for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q) { - MappingInfoStorage &my_data = + MappingInfoStorage &my_data = cell_data[my_q]; // step 4a: set the index offsets, find out how much to allocate, @@ -2800,8 +2749,8 @@ namespace internal // points, filling up all SIMD lanes as appropriate for (unsigned int my_q = 0; my_q < face_data.size(); ++my_q) { - MappingInfoStorage - &my_data = face_data[my_q]; + MappingInfoStorage &my_data = + face_data[my_q]; // step 6a: set the index offsets, find out how much to allocate, // and allocate the memory diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h index 4457046663..9dd239683c 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.h +++ b/include/deal.II/matrix_free/mapping_info_storage.h @@ -109,18 +109,18 @@ namespace internal * * @ingroup matrixfree */ - template + template struct MappingInfoStorage { - static_assert( - std::is_same::value, - "Type of Number and of VectorizedArrayType do not match."); - struct QuadratureDescriptor { + /** + * In case this class is instantiated for VectorizedArray types, this + * indicates the underlying scalar type for data which is the same on + * all lanes like the quadrature weights. + */ + using ScalarNumber = typename VectorizedArrayTrait::value_type; + /** * Constructor. Does nothing. */ @@ -167,14 +167,15 @@ namespace internal * Quadrature weights separated by dimension for use in specific * situations. */ - std::array, structdim> tensor_quadrature_weights; + std::array, structdim> + tensor_quadrature_weights; /** * A cached vector of quadrature weights in the given number format * (non-vectorized, as it is cheap to broadcast the value to all lanes * when it is used in a vectorized context). */ - AlignedVector quadrature_weights; + AlignedVector quadrature_weights; /** * For quadrature on faces, the evaluation of basis functions is not @@ -217,14 +218,14 @@ namespace internal * * Indexed by @p data_index_offsets. */ - AlignedVector JxW_values; + AlignedVector JxW_values; /** * Stores the normal vectors. * * Indexed by @p data_index_offsets. */ - AlignedVector> normal_vectors; + AlignedVector> normal_vectors; /** * The storage of covariant transformation on quadrature points, i.e., @@ -237,8 +238,7 @@ namespace internal * but the default case (cell integrals or boundary integrals) only * fills the zeroth component and ignores the first one. */ - std::array>, 2> - jacobians; + std::array>, 2> jacobians; /** * The storage of the gradients of the inverse Jacobian @@ -256,9 +256,8 @@ namespace internal * fills the zeroth component and ignores the first one. */ std::array< - AlignedVector>>, + AlignedVector< + Tensor<1, spacedim *(spacedim + 1) / 2, Tensor<1, spacedim, Number>>>, 2> jacobian_gradients; @@ -269,7 +268,7 @@ namespace internal * * Indexed by @p data_index_offsets. */ - std::array>, 2> + std::array>, 2> normals_times_jacobians; /** @@ -287,7 +286,7 @@ namespace internal * * Indexed by @p quadrature_point_offsets. */ - AlignedVector> quadrature_points; + AlignedVector> quadrature_points; /** * Clears all data fields except the descriptor vector. @@ -308,12 +307,20 @@ namespace internal * Get the descriptor from another instance of QuadratureDescriptor, * clearing all other data fields. */ - template + template void - copy_descriptor(const MappingInfoStorage &other); + copy_descriptor( + const MappingInfoStorage &other); + + /** + * Helper function to determine which update flags must be set in the + * internal functions to initialize all data as requested by the user. + */ + static UpdateFlags + compute_update_flags( + const UpdateFlags update_flags, + const std::vector> &quad = + std::vector>()); /** * Prints a detailed summary of memory consumption in the different @@ -335,13 +342,10 @@ namespace internal /* ------------------- inline functions ----------------------------- */ - template + template inline unsigned int - MappingInfoStorage:: - quad_index_from_n_q_points(const unsigned int n_q_points) const + MappingInfoStorage::quad_index_from_n_q_points( + const unsigned int n_q_points) const { for (unsigned int i = 0; i < descriptor.size(); ++i) if (n_q_points == descriptor[i].n_q_points) 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 aff20a206e..d4c15ef09a 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.templates.h +++ b/include/deal.II/matrix_free/mapping_info_storage.templates.h @@ -40,27 +40,20 @@ namespace internal { namespace MatrixFreeFunctions { - template - MappingInfoStorage:: - QuadratureDescriptor::QuadratureDescriptor() + template + MappingInfoStorage::QuadratureDescriptor:: + QuadratureDescriptor() : n_q_points(numbers::invalid_unsigned_int) {} - template + template template void - MappingInfoStorage:: - QuadratureDescriptor::initialize( - const Quadrature &quadrature, - const UpdateFlags update_flags_inner_faces) + MappingInfoStorage::QuadratureDescriptor:: + initialize(const Quadrature &quadrature, + const UpdateFlags update_flags_inner_faces) { Assert(structdim + 1 <= spacedim || update_flags_inner_faces == update_default, @@ -79,15 +72,11 @@ namespace internal - template + template void - MappingInfoStorage:: - QuadratureDescriptor::initialize( - const Quadrature<1> &quadrature_1d, - const UpdateFlags update_flags_inner_faces) + MappingInfoStorage::QuadratureDescriptor:: + initialize(const Quadrature<1> &quadrature_1d, + const UpdateFlags update_flags_inner_faces) { Assert(structdim + 1 <= spacedim || update_flags_inner_faces == update_default, @@ -137,13 +126,10 @@ namespace internal - template + template std::size_t - MappingInfoStorage:: - QuadratureDescriptor::memory_consumption() const + MappingInfoStorage::QuadratureDescriptor:: + memory_consumption() const { std::size_t memory = sizeof(this) + quadrature.memory_consumption() + quadrature_weights.memory_consumption() + @@ -155,13 +141,9 @@ namespace internal - template + template void - MappingInfoStorage:: - clear_data_fields() + MappingInfoStorage::clear_data_fields() { data_index_offsets.clear(); JxW_values.clear(); @@ -178,17 +160,11 @@ namespace internal - template - template + template + template void - MappingInfoStorage:: - copy_descriptor(const MappingInfoStorage &other) + MappingInfoStorage::copy_descriptor( + const MappingInfoStorage &other) { clear_data_fields(); descriptor.clear(); @@ -219,13 +195,56 @@ namespace internal - template + template + UpdateFlags + MappingInfoStorage::compute_update_flags( + const UpdateFlags update_flags, + const std::vector> &quad) + { + // this class is build around the evaluation of jacobians, so compute + // them in any case. The Jacobians will be inverted manually. Since we + // always do support integration, we also include the JxW values + UpdateFlags new_flags = update_jacobians | update_JxW_values; + + // for Hessian information, need inverse Jacobians and the derivative of + // Jacobians (these two together will give use the gradients of the + // inverse Jacobians, which is what we need) + if (update_flags & update_hessians || + update_flags & update_jacobian_grads) + new_flags |= update_jacobian_grads; + + if (update_flags & update_quadrature_points) + new_flags |= update_quadrature_points; + + // there is one more thing: if we have a quadrature formula with only + // one quadrature point on the first component, but more points on later + // components, we need to have Jacobian gradients anyway in order to + // determine whether the Jacobian is constant throughout a cell + if (quad.empty() == false) + { + bool formula_with_one_point = false; + for (unsigned int i = 0; i < quad[0].size(); ++i) + if (quad[0][i].size() == 1) + { + formula_with_one_point = true; + break; + } + if (formula_with_one_point == true) + for (unsigned int comp = 1; comp < quad.size(); ++comp) + for (unsigned int i = 0; i < quad[comp].size(); ++i) + if (quad[comp][i].size() > 1) + { + new_flags |= update_jacobian_grads; + } + } + return new_flags; + } + + + + template std::size_t - MappingInfoStorage:: - memory_consumption() const + MappingInfoStorage::memory_consumption() const { return MemoryConsumption::memory_consumption(descriptor) + MemoryConsumption::memory_consumption(data_index_offsets) + @@ -243,14 +262,12 @@ namespace internal - template + template template void - MappingInfoStorage:: - print_memory_consumption(StreamType &out, const TaskInfo &task_info) const + MappingInfoStorage::print_memory_consumption( + StreamType & out, + const TaskInfo &task_info) const { // print_memory_statistics involves global communication, so we can // disable the check here only if no processor has any such data diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 319e35ccb5..64f18fae96 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -1067,8 +1067,8 @@ namespace MatrixFreeOperators internal::CellwiseInverseMassMatrixImplBasic:: template run(n_components, fe_eval, in_array, out_array); else - internal::CellwiseInverseMassFactory:: - apply(n_components, fe_eval, in_array, out_array); + internal::CellwiseInverseMassFactory::apply( + n_components, fe_eval, in_array, out_array); } @@ -1101,13 +1101,13 @@ namespace MatrixFreeOperators in_array, out_array); else - internal::CellwiseInverseMassFactory:: - apply(n_actual_components, - given_degree, - fe_eval.get_shape_info().data.front().inverse_shape_values_eo, - inverse_coefficients, - in_array, - out_array); + internal::CellwiseInverseMassFactory::apply( + n_actual_components, + given_degree, + fe_eval.get_shape_info().data.front().inverse_shape_values_eo, + inverse_coefficients, + in_array, + out_array); } @@ -1138,7 +1138,7 @@ namespace MatrixFreeOperators in_array, out_array); else - internal::CellwiseInverseMassFactory:: + internal::CellwiseInverseMassFactory:: transform_from_q_points_to_basis(n_actual_components, fe_eval, in_array, diff --git a/source/matrix_free/evaluation_template_factory.inst.in b/source/matrix_free/evaluation_template_factory.inst.in index 468f6f3724..6e99a7faff 100644 --- a/source/matrix_free/evaluation_template_factory.inst.in +++ b/source/matrix_free/evaluation_template_factory.inst.in @@ -17,19 +17,18 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) { - template struct dealii::internal::FEEvaluationFactory< - deal_II_dimension, - deal_II_scalar_vectorized::value_type, - deal_II_scalar_vectorized>; + template struct dealii::internal:: + FEEvaluationFactory; + + template struct dealii::internal:: + FEFaceEvaluationFactory; - template struct dealii::internal::FEFaceEvaluationFactory< + template struct dealii::internal::FEFaceEvaluationGatherFactory< deal_II_dimension, deal_II_scalar_vectorized::value_type, deal_II_scalar_vectorized>; // inverse mass - template struct dealii::internal::CellwiseInverseMassFactory< - deal_II_dimension, - deal_II_scalar_vectorized::value_type, - deal_II_scalar_vectorized>; + template struct dealii::internal:: + CellwiseInverseMassFactory; } diff --git a/source/matrix_free/mapping_info.inst.in b/source/matrix_free/mapping_info.inst.in index 7b3996bf0b..1ebb0afa32 100644 --- a/source/matrix_free/mapping_info.inst.in +++ b/source/matrix_free/mapping_info.inst.in @@ -25,14 +25,12 @@ for (deal_II_dimension : DIMENSIONS; template struct internal::MatrixFreeFunctions::MappingInfoStorage< deal_II_dimension, deal_II_dimension, - deal_II_scalar_vectorized::value_type, deal_II_scalar_vectorized>; #if deal_II_dimension > 1 template struct internal::MatrixFreeFunctions::MappingInfoStorage< deal_II_dimension - 1, deal_II_dimension, - deal_II_scalar_vectorized::value_type, deal_II_scalar_vectorized>; #endif -- 2.39.5