From a9a0e442a0306c6f8a7674b0d72d9988e18bada0 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 3 Sep 2020 20:39:03 +0200 Subject: [PATCH] Change --- .../deal.II/matrix_free/evaluation_kernels.h | 181 +++++++----------- .../matrix_free/tensor_product_kernels.h | 21 ++ source/multigrid/mg_transfer_matrix_free.cc | 122 ++++++------ 3 files changed, 154 insertions(+), 170 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index ffd07d08a8..5c4010eea4 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -648,10 +648,11 @@ namespace internal * This class allows for dimension-independent application of the operation, * implemented by template recursion. It has been tested up to 6D. */ - template struct FEEvaluationImplBasisChange @@ -696,6 +697,8 @@ namespace internal basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable, ExcMessage("The second dimension must not be smaller than the first")); + Assert(quantity == EvaluatorQuantity::value, ExcInternalError()); + // we do recursion until dim==1 or dim==2 and we have // basis_size_1==basis_size_2. The latter optimization increases // optimization possibilities for the compiler but does only work for @@ -739,6 +742,7 @@ namespace internal for (unsigned int q = np_1; q != 0; --q) FEEvaluationImplBasisChange< variant, + quantity, next_dim, basis_size_1, basis_size_2, @@ -821,6 +825,10 @@ namespace internal "Input and output cannot alias with each other when " "adding the result of the basis change to existing data")); + Assert(quantity == EvaluatorQuantity::value || + quantity == EvaluatorQuantity::hessian, + ExcInternalError()); + constexpr int next_dim = (dim > 2 || ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ? @@ -833,8 +841,8 @@ namespace internal Number, Number2> eval_val(transformation_matrix, - AlignedVector(), - AlignedVector(), + transformation_matrix, + transformation_matrix, basis_size_1_variable, basis_size_2_variable); const unsigned int np_1 = @@ -850,27 +858,65 @@ namespace internal { if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2) { - eval_val.template values<1, false, false>(values_in, values_in); + if (quantity == EvaluatorQuantity::value) + eval_val.template values<1, false, false>(values_in, values_in); + else + eval_val.template hessians<1, false, false>(values_in, + values_in); + if (add_into_result) - eval_val.template values<0, false, true>(values_in, values_out); + { + if (quantity == EvaluatorQuantity::value) + eval_val.template values<0, false, true>(values_in, + values_out); + else + eval_val.template hessians<0, false, true>(values_in, + values_out); + } else - eval_val.template values<0, false, false>(values_in, - values_out); + { + if (quantity == EvaluatorQuantity::value) + eval_val.template values<0, false, false>(values_in, + values_out); + else + eval_val.template hessians<0, false, false>(values_in, + values_out); + } } else { if (dim == 1 && add_into_result) - eval_val.template values<0, false, true>(values_in, values_out); + { + if (quantity == EvaluatorQuantity::value) + eval_val.template values<0, false, true>(values_in, + values_out); + else + eval_val.template hessians<0, false, true>(values_in, + values_out); + } else if (dim == 1) - eval_val.template values<0, false, false>(values_in, - values_out); + { + if (quantity == EvaluatorQuantity::value) + eval_val.template values<0, false, false>(values_in, + values_out); + else + eval_val.template hessians<0, false, false>(values_in, + values_out); + } else - eval_val.template values(values_in, - values_in); + { + if (quantity == EvaluatorQuantity::value) + eval_val.template values(values_in, + values_in); + else + eval_val.template hessians( + values_in, values_in); + } } if (next_dim < dim) for (unsigned int q = 0; q < np_1; ++q) FEEvaluationImplBasisChange &transformation_matrix, - const bool add_into_result, - Number * values_in, - Number * values_out, - const unsigned int basis_size_1_variable = - dealii::numbers::invalid_unsigned_int, - const unsigned int basis_size_2_variable = - dealii::numbers::invalid_unsigned_int) - { - Assert(basis_size_1 != 0 || - basis_size_1_variable <= basis_size_2_variable, - dealii::ExcMessage( - "The second dimension must not be smaller than the first")); - Assert(add_into_result == false || values_in != values_out, - dealii::ExcMessage( - "Input and output cannot alias with each other when " - "adding the result of the basis change to existing data")); - - constexpr int next_dim = - (dim > 2 || - ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ? - dim - 1 : - dim; - dealii::internal::EvaluatorTensorProduct< - variant, - dim, - basis_size_1, - (basis_size_1 == 0 ? 0 : basis_size_2), - Number, - Number2> - eval_val(dealii::AlignedVector(), - dealii::AlignedVector(), - transformation_matrix, - basis_size_1_variable, - basis_size_2_variable); - const unsigned int np_1 = - basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable; - const unsigned int np_2 = - basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable; - Assert(np_1 > 0 && np_1 != dealii::numbers::invalid_unsigned_int, - dealii::ExcMessage("Cannot transform with 0-point basis")); - Assert(np_2 > 0 && np_2 != dealii::numbers::invalid_unsigned_int, - dealii::ExcMessage("Cannot transform with 0-point basis")); - - for (unsigned int c = 0; c < n_components; ++c) - { - if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2) - { - eval_val.template hessians<1, false, false>(values_in, values_in); - if (add_into_result) - eval_val.template hessians<0, false, true>(values_in, - values_out); - else - eval_val.template hessians<0, false, false>(values_in, - values_out); - } - else - { - if (dim == 1 && add_into_result) - eval_val.template hessians<0, false, true>(values_in, - values_out); - else if (dim == 1) - eval_val.template hessians<0, false, false>(values_in, - values_out); - else - eval_val.template hessians(values_in, - values_in); - } - if (next_dim < dim) - for (unsigned int q = 0; q < np_1; ++q) - FEEvaluationImplBasisChange:: - do_backward_hessians( - 1, - transformation_matrix, - add_into_result, - values_in + - q * dealii::Utilities::fixed_power(np_2), - values_out + - q * dealii::Utilities::fixed_power(np_1), - basis_size_1_variable, - basis_size_2_variable); - - values_in += dealii::Utilities::fixed_power(np_2); - values_out += dealii::Utilities::fixed_power(np_1); - } - } - /** * This operation applies a mass-matrix-like operation, consisting of a * do_forward() operation, multiplication by the coefficients in the @@ -1033,6 +978,7 @@ namespace internal for (unsigned int q = basis_size_1; q != 0; --q) FEEvaluationImplBasisChange< variant, + EvaluatorQuantity::value, next_dim, basis_size_1, basis_size_2, @@ -1073,6 +1019,7 @@ namespace internal for (unsigned int q = 0; q < basis_size_1; ++q) FEEvaluationImplBasisChange< variant, + EvaluatorQuantity::value, next_dim, basis_size_1, basis_size_2, @@ -1330,6 +1277,7 @@ namespace internal { FEEvaluationImplBasisChange< evaluate_evenodd, + EvaluatorQuantity::value, dim, (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), n_q_points_1d, @@ -1405,6 +1353,7 @@ namespace internal // transform back to the original space FEEvaluationImplBasisChange< evaluate_evenodd, + EvaluatorQuantity::value, dim, (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), n_q_points_1d, diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 5d62238c83..278f486a8d 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -67,6 +67,27 @@ namespace internal + /** + * Determine which quantity should be computed via the tensor product kernels. + */ + enum class EvaluatorQuantity + { + /** + * Evaluate/integrate by shape functions. + */ + value, + /** + * Evaluate/integrate by gradients of the shape functions. + */ + gradient, + /** + * Evaluate/integrate by hessians of the shape functions. + */ + hessian + }; + + + /** * Generic evaluator framework that valuates the given shape data in general * dimensions using the tensor product form. Depending on the particular diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index 153f0d28b6..f68c62a40b 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -458,19 +458,22 @@ MGTransferMatrixFree::do_prolongate_add( // must go through the components backwards because we want to write // the output to the same array as the input for (int c = n_components - 1; c >= 0; --c) - internal::FEEvaluationImplBasisChange, - VectorizedArray>:: - do_forward(1, - prolongation_matrix_1d, - evaluation_data.begin() + - c * Utilities::fixed_power(degree_size), - evaluation_data.begin() + c * n_scalar_cell_dofs, - fe_degree + 1, - 2 * fe_degree + 1); + internal::FEEvaluationImplBasisChange< + internal::evaluate_general, + internal::EvaluatorQuantity::value, + dim, + degree + 1, + 2 * degree + 1, + VectorizedArray, + VectorizedArray>::do_forward(1, + prolongation_matrix_1d, + evaluation_data.begin() + + c * Utilities::fixed_power< + dim>(degree_size), + evaluation_data.begin() + + c * n_scalar_cell_dofs, + fe_degree + 1, + 2 * fe_degree + 1); weight_dofs_on_child( &weights_on_refined[to_level - 1][(cell / vec_size) * three_to_dim], n_components, @@ -480,19 +483,22 @@ MGTransferMatrixFree::do_prolongate_add( else { for (int c = n_components - 1; c >= 0; --c) - internal::FEEvaluationImplBasisChange, - VectorizedArray>:: - do_forward(1, - prolongation_matrix_1d, - evaluation_data.begin() + - c * Utilities::fixed_power(degree_size), - evaluation_data.begin() + c * n_scalar_cell_dofs, - fe_degree + 1, - 2 * fe_degree + 2); + internal::FEEvaluationImplBasisChange< + internal::evaluate_general, + internal::EvaluatorQuantity::value, + dim, + degree + 1, + 2 * degree + 2, + VectorizedArray, + VectorizedArray>::do_forward(1, + prolongation_matrix_1d, + evaluation_data.begin() + + c * Utilities::fixed_power< + dim>(degree_size), + evaluation_data.begin() + + c * n_scalar_cell_dofs, + fe_degree + 1, + 2 * fe_degree + 2); } // write into dst vector @@ -556,38 +562,46 @@ MGTransferMatrixFree::do_restrict_add( fe_degree, evaluation_data.data()); for (unsigned int c = 0; c < n_components; ++c) - internal::FEEvaluationImplBasisChange, - VectorizedArray>:: - do_backward(1, - prolongation_matrix_1d, - false, - evaluation_data.begin() + c * n_scalar_cell_dofs, - evaluation_data.begin() + - c * Utilities::fixed_power(degree_size), - fe_degree + 1, - 2 * fe_degree + 1); + internal::FEEvaluationImplBasisChange< + internal::evaluate_general, + internal::EvaluatorQuantity::value, + dim, + degree + 1, + 2 * degree + 1, + VectorizedArray, + VectorizedArray>::do_backward(1, + prolongation_matrix_1d, + false, + evaluation_data.begin() + + c * n_scalar_cell_dofs, + evaluation_data.begin() + + c * + Utilities::fixed_power< + dim>(degree_size), + fe_degree + 1, + 2 * fe_degree + 1); } else { for (unsigned int c = 0; c < n_components; ++c) - internal::FEEvaluationImplBasisChange, - VectorizedArray>:: - do_backward(1, - prolongation_matrix_1d, - false, - evaluation_data.begin() + c * n_scalar_cell_dofs, - evaluation_data.begin() + - c * Utilities::fixed_power(degree_size), - fe_degree + 1, - 2 * fe_degree + 2); + internal::FEEvaluationImplBasisChange< + internal::evaluate_general, + internal::EvaluatorQuantity::value, + dim, + degree + 1, + 2 * degree + 2, + VectorizedArray, + VectorizedArray>::do_backward(1, + prolongation_matrix_1d, + false, + evaluation_data.begin() + + c * n_scalar_cell_dofs, + evaluation_data.begin() + + c * + Utilities::fixed_power< + dim>(degree_size), + fe_degree + 1, + 2 * fe_degree + 2); } // write into dst vector -- 2.39.5