From b8421de607b395cabbdf7294b86309307a9cd60f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 8 Dec 2019 20:09:19 +0100 Subject: [PATCH] Use inverse shape values for inverse mass operation --- .../deal.II/matrix_free/evaluation_kernels.h | 27 +++++----- include/deal.II/matrix_free/operators.h | 53 ++++++------------- 2 files changed, 30 insertions(+), 50 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index e77ae68e11..4ab979c2aa 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3037,6 +3037,9 @@ namespace internal + /** + * This struct implements the action of the inverse mass matrix operation + */ template struct CellwiseInverseMassMatrixImpl { @@ -3078,30 +3081,30 @@ namespace internal Number * out = out_array + d * dofs_per_component; // Need to select 'apply' method with hessian slot because values // assume symmetries that do not exist in the inverse shapes - evaluator.template hessians<0, false, false>(in, out); + evaluator.template hessians<0, true, false>(in, out); if (dim > 1) { - evaluator.template hessians<1, false, false>(out, out); + evaluator.template hessians<1, true, false>(out, out); if (dim == 3) { - evaluator.template hessians<2, false, false>(out, out); + evaluator.template hessians<2, true, false>(out, out); for (unsigned int q = 0; q < dofs_per_component; ++q) out[q] *= inv_coefficient[q]; - evaluator.template hessians<2, true, false>(out, out); + evaluator.template hessians<2, false, false>(out, out); } else if (dim == 2) for (unsigned int q = 0; q < dofs_per_component; ++q) out[q] *= inv_coefficient[q]; - evaluator.template hessians<1, true, false>(out, out); + evaluator.template hessians<1, false, false>(out, out); } else { for (unsigned int q = 0; q < dofs_per_component; ++q) out[q] *= inv_coefficient[q]; } - evaluator.template hessians<0, true, false>(out, out); + evaluator.template hessians<0, false, false>(out, out); inv_coefficient += shift_coefficient; } @@ -3130,17 +3133,17 @@ namespace internal if (dim == 3) { - evaluator.template hessians<2, true, false>(in, out); - evaluator.template hessians<1, true, false>(out, out); - evaluator.template hessians<0, true, false>(out, out); + 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, true, false>(in, out); - evaluator.template hessians<0, true, false>(out, out); + evaluator.template hessians<1, false, false>(in, out); + evaluator.template hessians<0, false, false>(out, out); } if (dim == 1) - evaluator.template hessians<0, true, false>(in, out); + evaluator.template hessians<0, false, false>(in, out); } } }; diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 2920c3837b..abdc4cf9a2 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -709,11 +709,6 @@ namespace MatrixFreeOperators Number, false, VectorizedArrayType> &fe_eval; - - /** - * A structure to hold inverse shape functions - */ - AlignedVector inverse_shape; }; @@ -958,26 +953,7 @@ namespace MatrixFreeOperators false, VectorizedArrayType> &fe_eval) : fe_eval(fe_eval) - { - FullMatrix shapes_1d(fe_degree + 1, fe_degree + 1); - for (unsigned int i = 0, c = 0; i < shapes_1d.m(); ++i) - for (unsigned int j = 0; j < shapes_1d.n(); ++j, ++c) - shapes_1d(i, j) = fe_eval.get_shape_info().shape_values[c][0]; - shapes_1d.gauss_jordan(); - const unsigned int stride = (fe_degree + 2) / 2; - inverse_shape.resize(stride * (fe_degree + 1)); - for (unsigned int i = 0; i < stride; ++i) - for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q) - { - inverse_shape[i * stride + q] = - 0.5 * (shapes_1d(i, q) + shapes_1d(i, fe_degree - q)); - inverse_shape[(fe_degree - i) * stride + q] = - 0.5 * (shapes_1d(i, q) - shapes_1d(i, fe_degree - q)); - } - if (fe_degree % 2 == 0) - for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q) - inverse_shape[fe_degree / 2 * stride + q] = shapes_1d(fe_degree / 2, q); - } + {} @@ -1030,15 +1006,15 @@ namespace MatrixFreeOperators const VectorizedArrayType * in_array, VectorizedArrayType * out_array) const { - internal::CellwiseInverseMassMatrixImpl< - dim, - fe_degree, - n_components, - VectorizedArrayType>::apply(inverse_shape, - inverse_coefficients, - n_actual_components, - in_array, - out_array); + internal::CellwiseInverseMassMatrixImpl:: + apply(fe_eval.get_shape_info().inverse_shape_values_eo, + inverse_coefficients, + n_actual_components, + in_array, + out_array); } @@ -1062,10 +1038,11 @@ namespace MatrixFreeOperators fe_degree, n_components, VectorizedArrayType>:: - transform_from_q_points_to_basis(inverse_shape, - n_actual_components, - in_array, - out_array); + transform_from_q_points_to_basis( + fe_eval.get_shape_info().inverse_shape_values_eo, + n_actual_components, + in_array, + out_array); } -- 2.39.5