From: Martin Kronbichler Date: Fri, 18 Sep 2020 15:10:38 +0000 (+0200) Subject: Add inverse mass matrix with degree -1 X-Git-Tag: v9.3.0-rc1~1087^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d4697e22dce73690da44523ff57974d2c83a3b04;p=dealii.git Add inverse mass matrix with degree -1 --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index cbe9fe3b4f..cd10369e7b 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3529,7 +3529,8 @@ namespace internal false, Number> &fe_eval, const Number * in_array, - Number * out_array) + Number * out_array, + typename std::enable_if::type * = nullptr) { constexpr unsigned int dofs_per_component = Utilities::pow(fe_degree + 1, dim); @@ -3578,6 +3579,61 @@ namespace internal } return false; } + + template + static bool + run(const unsigned int n_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 unsigned int dofs_per_component = + fe_eval.get_shape_info().dofs_per_component_on_cell; + + Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); + + internal:: + EvaluatorTensorProduct + evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values, + AlignedVector(), + AlignedVector(), + fe_eval.get_shape_info().data.front().fe_degree + 1, + fe_eval.get_shape_info().data.front().fe_degree + 1); + + 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 values<0, true, false>(in, out); + if (dim > 1) + evaluator.template values<1, true, false>(out, out); + if (dim > 2) + evaluator.template values<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 values<2, false, false>(out, out); + if (dim > 1) + evaluator.template values<1, false, false>(out, out); + evaluator.template values<0, false, false>(out, out); + } + return false; + } }; @@ -3595,7 +3651,8 @@ namespace internal const AlignedVector &inverse_shape, const AlignedVector &inverse_coefficients, const Number * in_array, - Number * out_array) + Number * out_array, + typename std::enable_if::type * = nullptr) { constexpr unsigned int dofs_per_component = Utilities::pow(fe_degree + 1, dim); @@ -3647,6 +3704,23 @@ namespace internal } return false; } + + /** + * Version for degree = -1 + */ + template + static bool + run(const unsigned int, + const AlignedVector &, + const AlignedVector &, + const Number *, + Number *, + typename std::enable_if::type * = nullptr) + { + static_assert(fe_degree == -1, "Only usable for degree -1"); + Assert(false, ExcNotImplemented()); + return false; + } }; @@ -3663,7 +3737,8 @@ namespace internal run(const unsigned int n_desired_components, const AlignedVector &inverse_shape, const Number * in_array, - Number * out_array) + Number * out_array, + typename std::enable_if::type * = nullptr) { constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim); internal::EvaluatorTensorProduct + static bool + run(const unsigned int, + const AlignedVector &, + const Number *, + Number *, + typename std::enable_if::type * = nullptr) + { + static_assert(fe_degree == -1, "Only usable for degree -1"); + Assert(false, ExcNotImplemented()); + return false; + } }; } // end of namespace internal