From 5566be307498b3b9d9cf6435584c824b8e516cb9 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@it.uu.se>
Date: Mon, 29 Nov 2021 11:18:02 +0100
Subject: [PATCH] Use pre-compiled code for CellwiseInverseMassMatrix

---
 .../changes/minor/20211129MartinKronbichler   |  6 ++++
 .../matrix_free/evaluation_template_factory.h |  3 --
 .../evaluation_template_factory.templates.h   |  7 ++--
 include/deal.II/matrix_free/operators.h       | 36 +++++++++++++------
 4 files changed, 35 insertions(+), 17 deletions(-)
 create mode 100644 doc/news/changes/minor/20211129MartinKronbichler

diff --git a/doc/news/changes/minor/20211129MartinKronbichler b/doc/news/changes/minor/20211129MartinKronbichler
new file mode 100644
index 0000000000..eea4c28eb2
--- /dev/null
+++ b/doc/news/changes/minor/20211129MartinKronbichler
@@ -0,0 +1,6 @@
+Fixed: MatrixFreeOperators::CellwiseInverseMassMatrix with template argument
+`-1` for the polynomial degree now chooses the pre-compiled code with fast
+tensor-product algorithms if available, rather than a slower general-purpose
+code.
+<br>
+(Martin Kronbichler, 2021/11/29)
diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h
index 8afb37d092..50ac727c87 100644
--- a/include/deal.II/matrix_free/evaluation_template_factory.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.h
@@ -172,7 +172,6 @@ namespace internal
   {
     static void
     apply(const unsigned int n_components,
-          const unsigned int fe_degree,
           const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
             &                        fe_eval,
           const VectorizedArrayType *in_array,
@@ -189,8 +188,6 @@ namespace internal
     static void
     transform_from_q_points_to_basis(
       const unsigned int n_components,
-      const unsigned int fe_degree,
-      const unsigned int n_q_points_1d,
       const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
         &                        fe_eval,
       const VectorizedArrayType *in_array,
diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
index b246ee4e44..5a98c3781b 100644
--- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -330,12 +330,12 @@ namespace internal
   void
   CellwiseInverseMassFactory<dim, Number, VectorizedArrayType>::apply(
     const unsigned int n_components,
-    const unsigned int fe_degree,
     const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
       &                        fe_eval,
     const VectorizedArrayType *in_array,
     VectorizedArrayType *      out_array)
   {
+    const unsigned int fe_degree = fe_eval.get_shape_info().data[0].fe_degree;
     instantiation_helper_run<
       1,
       CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>>(
@@ -373,13 +373,14 @@ namespace internal
   CellwiseInverseMassFactory<dim, Number, VectorizedArrayType>::
     transform_from_q_points_to_basis(
       const unsigned int n_components,
-      const unsigned int fe_degree,
-      const unsigned int n_q_points_1d,
       const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
         &                        fe_eval,
       const VectorizedArrayType *in_array,
       VectorizedArrayType *      out_array)
   {
+    const unsigned int fe_degree = fe_eval.get_shape_info().data[0].fe_degree;
+    const unsigned int n_q_points_1d =
+      fe_eval.get_shape_info().data[0].n_q_points_1d;
     instantiation_helper_run<
       1,
       CellwiseInverseMassMatrixImplTransformFromQPoints<dim,
diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 017af8daa4..319e35ccb5 100644
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -1063,8 +1063,12 @@ namespace MatrixFreeOperators
     VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
                                 VectorizedArrayType *      out_array) const
   {
-    internal::CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>::
-      template run<fe_degree>(n_components, fe_eval, in_array, out_array);
+    if (fe_degree > -1)
+      internal::CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>::
+        template run<fe_degree>(n_components, fe_eval, in_array, out_array);
+    else
+      internal::CellwiseInverseMassFactory<dim, Number, VectorizedArrayType>::
+        apply(n_components, fe_eval, in_array, out_array);
   }
 
 
@@ -1085,13 +1089,25 @@ namespace MatrixFreeOperators
           const VectorizedArrayType *               in_array,
           VectorizedArrayType *                     out_array) const
   {
-    internal::CellwiseInverseMassMatrixImplFlexible<dim, VectorizedArrayType>::
-      template run<fe_degree>(
-        n_actual_components,
-        fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
-        inverse_coefficients,
-        in_array,
-        out_array);
+    const unsigned int given_degree =
+      fe_eval.get_shape_info().data[0].fe_degree;
+    if (fe_degree > -1)
+      internal::CellwiseInverseMassMatrixImplFlexible<dim,
+                                                      VectorizedArrayType>::
+        template run<fe_degree>(
+          n_actual_components,
+          fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
+          inverse_coefficients,
+          in_array,
+          out_array);
+    else
+      internal::CellwiseInverseMassFactory<dim, Number, VectorizedArrayType>::
+        apply(n_actual_components,
+              given_degree,
+              fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
+              inverse_coefficients,
+              in_array,
+              out_array);
   }
 
 
@@ -1124,8 +1140,6 @@ namespace MatrixFreeOperators
     else
       internal::CellwiseInverseMassFactory<dim, Number, VectorizedArrayType>::
         transform_from_q_points_to_basis(n_actual_components,
-                                         fe_degree,
-                                         n_q_points_1d,
                                          fe_eval,
                                          in_array,
                                          out_array);
-- 
2.39.5