From 927d8b7201caed18dc28384f8ddd3140cfa85f18 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 5 Jan 2020 18:20:49 +0100 Subject: [PATCH] Implement simpler CellwiseInverseMassMatrix::apply() function --- .../deal.II/matrix_free/evaluation_kernels.h | 80 ++++++++++++++----- include/deal.II/matrix_free/operators.h | 40 +++++++++- 2 files changed, 100 insertions(+), 20 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 4ab979c2aa..79623b3e8a 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3043,6 +3043,58 @@ namespace internal template struct CellwiseInverseMassMatrixImpl { + template + static void + apply(const FEEvaluationType &fe_eval, + const Number * in_array, + Number * out_array) + { + constexpr unsigned int dofs_per_component = + Utilities::pow(fe_degree + 1, dim); + + Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); + Assert(fe_eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric, + ExcNotImplemented()); + + internal::EvaluatorTensorProduct + evaluator(AlignedVector(), + AlignedVector(), + fe_eval.get_shape_info().inverse_shape_values_eo); + + for (unsigned int d = 0; d < n_components; ++d) + { + const Number *in = in_array + d * dofs_per_component; + 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, true, false>(in, out); + if (dim > 1) + evaluator.template hessians<1, true, false>(out, out); + if (dim > 2) + evaluator.template hessians<2, true, false>(out, out); + } + for (unsigned int q = 0; q < dofs_per_component; ++q) + { + const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q); + for (unsigned int d = 0; d < n_components; ++d) + out_array[q + d * dofs_per_component] *= inverse_JxW_q; + } + for (unsigned int d = 0; d < n_components; ++d) + { + Number *out = out_array + d * dofs_per_component; + if (dim > 2) + evaluator.template hessians<2, false, false>(out, out); + if (dim > 1) + evaluator.template hessians<1, false, false>(out, out); + evaluator.template hessians<0, false, false>(out, out); + } + } + static void apply(const AlignedVector &inverse_shape, const AlignedVector &inverse_coefficients, @@ -3083,27 +3135,17 @@ namespace internal // assume symmetries that do not exist in the inverse shapes evaluator.template hessians<0, true, false>(in, out); if (dim > 1) - { - evaluator.template hessians<1, true, false>(out, out); + evaluator.template hessians<1, true, false>(out, out); + if (dim > 2) + evaluator.template hessians<2, true, false>(out, out); - if (dim == 3) - { - 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, false, false>(out, out); - } - else if (dim == 2) - for (unsigned int q = 0; q < dofs_per_component; ++q) - out[q] *= inv_coefficient[q]; + for (unsigned int q = 0; q < dofs_per_component; ++q) + out[q] *= inv_coefficient[q]; - evaluator.template hessians<1, false, false>(out, out); - } - else - { - for (unsigned int q = 0; q < dofs_per_component; ++q) - out[q] *= inv_coefficient[q]; - } + if (dim > 2) + evaluator.template hessians<2, false, false>(out, out); + if (dim > 1) + evaluator.template hessians<1, false, false>(out, out); evaluator.template hessians<0, false, false>(out, out); inv_coefficient += shift_coefficient; diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index abdc4cf9a2..9d71deab6d 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -642,7 +642,7 @@ namespace MatrixFreeOperators /** * Applies the inverse mass matrix operation on an input array. It is * assumed that the passed input and output arrays are of correct size, - * namely FEEval::dofs_per_cell * n_components long. The inverse of the + * namely FEEvaluation::dofs_per_cell long. The inverse of the * local coefficient (also containing the inverse JxW values) must be * passed as first argument. Passing more than one component in the * coefficient is allowed. @@ -653,6 +653,21 @@ namespace MatrixFreeOperators const VectorizedArrayType * in_array, VectorizedArrayType * out_array) const; + /** + * Applies the inverse mass matrix operation on an input array, using the + * inverse of the JxW values provided by the `fe_eval` argument passed to + * the constructor of this class. Note that the user code must call + * FEEvaluation::reinit() on the underlying evaluator to make the + * FEEvaluationBase::JxW() method return the information of the correct + * cell. It is assumed that the pointers of the input and output arrays + * are valid over the length FEEvaluation::dofs_per_cell, which is the + * number of entries processed by this function. The `in_array` and + * `out_array` arguments may point to the same memory position. + */ + void + apply(const VectorizedArrayType *in_array, + VectorizedArrayType * out_array) const; + /** * This operation performs a projection from the data given in quadrature * points to the actual basis underlying this object. This projection can @@ -990,6 +1005,29 @@ namespace MatrixFreeOperators + template + inline void + CellwiseInverseMassMatrix< + dim, + fe_degree, + n_components, + Number, + VectorizedArrayType>::apply(const VectorizedArrayType *in_array, + VectorizedArrayType * out_array) const + { + internal::CellwiseInverseMassMatrixImpl< + dim, + fe_degree, + n_components, + VectorizedArrayType>::apply(fe_eval, in_array, out_array); + } + + + template