From: Peter Munch Date: Sat, 5 Jun 2021 15:07:04 +0000 (+0200) Subject: Refactor CellwiseInverseMassMatrixImplTransformFromQPoints and use CellwiseInverseMas... X-Git-Tag: v9.4.0-rc1~1274^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12392%2Fhead;p=dealii.git Refactor CellwiseInverseMassMatrixImplTransformFromQPoints and use CellwiseInverseMassFactory --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 4901a4a2df..7e8bb94e4f 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -4349,7 +4349,7 @@ namespace internal template struct CellwiseInverseMassMatrixImplTransformFromQPoints { - template + template static bool run(const unsigned int n_desired_components, const FEEvaluationBaseData &fe_eval, const Number * in_array, - Number * out_array, - typename std::enable_if::type * = nullptr) + Number * out_array) { - const AlignedVector &inverse_shape = - fe_eval.get_shape_info().data.front().inverse_shape_values_eo; + static const bool do_inplace = + fe_degree > -1 && (fe_degree + 1 == n_q_points_1d); - constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim); - internal::EvaluatorTensorProduct evaluator(AlignedVector(), AlignedVector(), - inverse_shape); - - for (unsigned int d = 0; d < n_desired_components; ++d) - { - const Number *in = in_array + d * dofs_per_cell; - Number * out = out_array + d * dofs_per_cell; - - if (dim == 3) - { - evaluator.template hessians<2, false, false>(in, out); - evaluator.template hessians<1, false, false>(out, out); - evaluator.template hessians<0, false, false>(out, out); - } - if (dim == 2) - { - evaluator.template hessians<1, false, false>(in, out); - evaluator.template hessians<0, false, false>(out, out); - } - if (dim == 1) - evaluator.template hessians<0, false, false>(in, out); - } - return false; - } - - template - static bool - run(const unsigned int n_desired_components, - const FEEvaluationBaseData &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"); - - const AlignedVector &inverse_shape = - fe_eval.get_shape_info().data.front().inverse_shape_values; - - const unsigned int dofs_per_component = - fe_eval.get_shape_info().dofs_per_component_on_cell; - const unsigned int n_q_points = fe_eval.get_shape_info().n_q_points; - - internal:: - EvaluatorTensorProduct - evaluator(inverse_shape, - AlignedVector(), - AlignedVector(), - fe_eval.get_shape_info().data.front().fe_degree + 1, - fe_eval.get_shape_info().data.front().n_q_points_1d); - - auto temp_1 = fe_eval.get_scratch_data().begin(); - auto temp_2 = temp_1 + std::max(n_q_points, dofs_per_component); + inverse_shape, + fe_eval.get_shape_info().data.front().fe_degree + 1, + fe_eval.get_shape_info().data.front().n_q_points_1d); for (unsigned int d = 0; d < n_desired_components; ++d) { const Number *in = in_array + d * n_q_points; Number * out = out_array + d * dofs_per_component; + auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin(); + auto temp_2 = do_inplace ? + out : + (temp_1 + std::max(n_q_points, dofs_per_component)); + if (dim == 3) { - evaluator.template values<2, false, false>(in, temp_1); - evaluator.template values<1, false, false>(temp_1, temp_2); - evaluator.template values<0, false, false>(temp_2, out); + evaluator.template hessians<2, false, false>(in, temp_1); + evaluator.template hessians<1, false, false>(temp_1, temp_2); + evaluator.template hessians<0, false, false>(temp_2, out); } if (dim == 2) { - evaluator.template values<1, false, false>(in, temp_1); - evaluator.template values<0, false, false>(temp_1, out); + evaluator.template hessians<1, false, false>(in, temp_1); + evaluator.template hessians<0, false, false>(temp_1, out); } if (dim == 1) - evaluator.template values<0, false, false>(in, out); + evaluator.template hessians<0, false, false>(in, out); } return false; } diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index 5b596874b3..9299bd4921 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -180,6 +180,7 @@ namespace internal transform_from_q_points_to_basis( const unsigned int n_components, const unsigned int fe_degree, + const unsigned int n_q_points_1d, const FEEvaluationBaseData & fe_eval, const VectorizedArrayType *in_array, diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index bbd5bf4443..db988b6bc4 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -389,6 +389,7 @@ namespace internal transform_from_q_points_to_basis( const unsigned int n_components, const unsigned int fe_degree, + const unsigned int n_q_points_1d, const FEEvaluationBaseData & fe_eval, const VectorizedArrayType *in_array, @@ -398,7 +399,7 @@ namespace internal 1, CellwiseInverseMassMatrixImplTransformFromQPoints>( - fe_degree, fe_degree + 1, n_components, fe_eval, in_array, out_array); + fe_degree, n_q_points_1d, n_components, fe_eval, in_array, out_array); } } // end of namespace internal diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 971b5d9e05..8465ea3415 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -1064,12 +1064,24 @@ namespace MatrixFreeOperators const VectorizedArrayType *in_array, VectorizedArrayType * out_array) const { - internal::CellwiseInverseMassMatrixImplTransformFromQPoints< - dim, - VectorizedArrayType>::template run(n_actual_components, - fe_eval, - in_array, - out_array); + const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d; + + if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d)) + internal::CellwiseInverseMassMatrixImplTransformFromQPoints< + dim, + VectorizedArrayType>::template run(n_actual_components, + fe_eval, + in_array, + out_array); + else + internal::CellwiseInverseMassFactory:: + transform_from_q_points_to_basis(n_actual_components, + fe_degree, + n_q_points_1d, + fe_eval, + in_array, + out_array); }