From d18a6e20dd6cf7bd803a983d51c6f328757a9239 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 4 Jul 2018 15:32:31 +0200 Subject: [PATCH] Enable 1D case in inverse mass matrix operator. --- .../deal.II/matrix_free/evaluation_kernels.h | 12 ++++- include/deal.II/matrix_free/operators.h | 51 +++++++++++-------- 2 files changed, 42 insertions(+), 21 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 1173c86f95..e6262b4822 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -917,6 +917,14 @@ namespace internal constexpr int next_dim = dim > 1 ? dim - 1 : dim; Number * my_scratch = basis_size_1 != basis_size_2 ? scratch_data : values_out; + + const unsigned int size_per_component = Utilities::pow(basis_size_2, dim); + Assert(coefficients.size() == size_per_component || + coefficients.size() == n_components * size_per_component, + ExcDimensionMismatch(coefficients.size(), size_per_component)); + const unsigned int stride = + coefficients.size() == size_per_component ? 0 : 1; + for (unsigned int q = basis_size_1; q != 0; --q) FEEvaluationImplBasisChange< variant, @@ -950,7 +958,9 @@ namespace internal my_scratch + i, my_scratch + i); for (unsigned int q = 0; q < basis_size_2; ++q) for (unsigned int i = ii; i < ii + n_inner_blocks; ++i) - my_scratch[i + q * n_blocks] *= coefficients[i + q * n_blocks]; + my_scratch[i + q * n_blocks + c * size_per_component] *= + coefficients[i + q * n_blocks + + c * stride * size_per_component]; for (unsigned int i = ii; i < ii + n_inner_blocks; ++i) eval_val.template values_one_line( my_scratch + i, my_scratch + i); diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 29f4b75b23..9f5d20f4a8 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -939,16 +939,17 @@ namespace MatrixFreeOperators const VectorizedArray * in_array, VectorizedArray * out_array) const { - constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim); + constexpr unsigned int dofs_per_component = + Utilities::pow(fe_degree + 1, dim); Assert(inverse_coefficients.size() > 0 && - inverse_coefficients.size() % dofs_per_cell == 0, + inverse_coefficients.size() % dofs_per_component == 0, ExcMessage( "Expected diagonal to be a multiple of scalar dof per cells")); - if (inverse_coefficients.size() != dofs_per_cell) - AssertDimension(n_actual_components * dofs_per_cell, + if (inverse_coefficients.size() != dofs_per_component) + AssertDimension(n_actual_components * dofs_per_component, inverse_coefficients.size()); - Assert(dim == 2 || dim == 3, ExcNotImplemented()); + Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); internal::EvaluatorTensorProduct dofs_per_cell ? dofs_per_cell : 0; + inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0; const VectorizedArray *inv_coefficient = &inverse_coefficients[0]; - VectorizedArray temp_data_field[dofs_per_cell]; + VectorizedArray temp_data_field[dofs_per_component]; for (unsigned int d = 0; d < n_actual_components; ++d) { - const VectorizedArray *in = in_array + d * dofs_per_cell; - VectorizedArray * out = out_array + d * dofs_per_cell; + const VectorizedArray *in = in_array + d * dofs_per_component; + VectorizedArray * 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, temp_data_field); - evaluator.template hessians<1, false, false>(temp_data_field, out); - - if (dim == 3) + if (dim > 1) { - evaluator.template hessians<2, false, false>(out, temp_data_field); - for (unsigned int q = 0; q < dofs_per_cell; ++q) + evaluator.template hessians<1, false, false>(temp_data_field, out); + + if (dim == 3) + { + evaluator.template hessians<2, false, false>(out, + temp_data_field); + for (unsigned int q = 0; q < dofs_per_component; ++q) + temp_data_field[q] *= inv_coefficient[q]; + evaluator.template hessians<2, true, false>(temp_data_field, + 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, temp_data_field); + } + else + { + for (unsigned int q = 0; q < dofs_per_component; ++q) temp_data_field[q] *= inv_coefficient[q]; - evaluator.template hessians<2, true, false>(temp_data_field, out); } - else if (dim == 2) - for (unsigned int q = 0; q < dofs_per_cell; ++q) - out[q] *= inv_coefficient[q]; - - evaluator.template hessians<1, true, false>(out, temp_data_field); evaluator.template hessians<0, true, false>(temp_data_field, out); inv_coefficient += shift_coefficient; -- 2.39.5