From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 5 Jun 2021 15:07:04 +0000 (+0200)
Subject: Refactor CellwiseInverseMassMatrixImplTransformFromQPoints and use CellwiseInverseMas... 
X-Git-Tag: v9.4.0-rc1~1274^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=80a756c902c1e75abf1b0867764ac8c4f15d3cbb;p=dealii.git

Refactor CellwiseInverseMassMatrixImplTransformFromQPoints and use CellwiseInverseMassFactory
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 4901a4a2df..7e8bb94e4f 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -4349,7 +4349,7 @@ namespace internal
   template <int dim, typename Number>
   struct CellwiseInverseMassMatrixImplTransformFromQPoints
   {
-    template <int fe_degree, int = 0>
+    template <int fe_degree, int n_q_points_1d>
     static bool
     run(const unsigned int                  n_desired_components,
         const FEEvaluationBaseData<dim,
@@ -4357,93 +4357,62 @@ namespace internal
                                    false,
                                    Number> &fe_eval,
         const Number *                      in_array,
-        Number *                            out_array,
-        typename std::enable_if<fe_degree != -1>::type * = nullptr)
+        Number *                            out_array)
     {
-      const AlignedVector<Number> &inverse_shape =
-        fe_eval.get_shape_info().data.front().inverse_shape_values_eo;
+      static const bool do_inplace =
+        fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
 
-      constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
-      internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
+      Assert(fe_eval.get_shape_info().element_type !=
+               MatrixFreeFunctions::ElementType::tensor_none,
+             ExcNotImplemented());
+
+      const auto &inverse_shape =
+        do_inplace ?
+          fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
+          fe_eval.get_shape_info().data.front().inverse_shape_values;
+
+      const unsigned int dofs_per_component =
+        do_inplace ? Utilities::pow(fe_degree + 1, dim) :
+                     fe_eval.get_shape_info().dofs_per_component_on_cell;
+      const unsigned int n_q_points = do_inplace ?
+                                        Utilities::pow(fe_degree + 1, dim) :
+                                        fe_eval.get_shape_info().n_q_points;
+
+      internal::EvaluatorTensorProduct<do_inplace ? internal::evaluate_evenodd :
+                                                    internal::evaluate_general,
                                        dim,
                                        fe_degree + 1,
-                                       fe_degree + 1,
+                                       n_q_points_1d,
                                        Number>
         evaluator(AlignedVector<Number>(),
                   AlignedVector<Number>(),
-                  inverse_shape);
-
-      for (unsigned int d = 0; d < n_desired_components; ++d)
-        {
-          const Number *in  = in_array + d * dofs_per_cell;
-          Number *      out = out_array + d * dofs_per_cell;
-
-          if (dim == 3)
-            {
-              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, false, false>(in, out);
-              evaluator.template hessians<0, false, false>(out, out);
-            }
-          if (dim == 1)
-            evaluator.template hessians<0, false, false>(in, out);
-        }
-      return false;
-    }
-
-    template <int fe_degree, int = 0>
-    static bool
-    run(const unsigned int                  n_desired_components,
-        const FEEvaluationBaseData<dim,
-                                   typename Number::value_type,
-                                   false,
-                                   Number> &fe_eval,
-        const Number *                      in_array,
-        Number *                            out_array,
-        typename std::enable_if<fe_degree == -1>::type * = nullptr)
-    {
-      static_assert(fe_degree == -1, "Only usable for degree -1");
-
-      const AlignedVector<Number> &inverse_shape =
-        fe_eval.get_shape_info().data.front().inverse_shape_values;
-
-      const unsigned int dofs_per_component =
-        fe_eval.get_shape_info().dofs_per_component_on_cell;
-      const unsigned int n_q_points = fe_eval.get_shape_info().n_q_points;
-
-      internal::
-        EvaluatorTensorProduct<internal::evaluate_general, dim, 0, 0, Number>
-          evaluator(inverse_shape,
-                    AlignedVector<Number>(),
-                    AlignedVector<Number>(),
-                    fe_eval.get_shape_info().data.front().fe_degree + 1,
-                    fe_eval.get_shape_info().data.front().n_q_points_1d);
-
-      auto temp_1 = fe_eval.get_scratch_data().begin();
-      auto temp_2 = temp_1 + std::max(n_q_points, dofs_per_component);
+                  inverse_shape,
+                  fe_eval.get_shape_info().data.front().fe_degree + 1,
+                  fe_eval.get_shape_info().data.front().n_q_points_1d);
 
       for (unsigned int d = 0; d < n_desired_components; ++d)
         {
           const Number *in  = in_array + d * n_q_points;
           Number *      out = out_array + d * dofs_per_component;
 
+          auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
+          auto temp_2 = do_inplace ?
+                          out :
+                          (temp_1 + std::max(n_q_points, dofs_per_component));
+
           if (dim == 3)
             {
-              evaluator.template values<2, false, false>(in, temp_1);
-              evaluator.template values<1, false, false>(temp_1, temp_2);
-              evaluator.template values<0, false, false>(temp_2, out);
+              evaluator.template hessians<2, false, false>(in, temp_1);
+              evaluator.template hessians<1, false, false>(temp_1, temp_2);
+              evaluator.template hessians<0, false, false>(temp_2, out);
             }
           if (dim == 2)
             {
-              evaluator.template values<1, false, false>(in, temp_1);
-              evaluator.template values<0, false, false>(temp_1, out);
+              evaluator.template hessians<1, false, false>(in, temp_1);
+              evaluator.template hessians<0, false, false>(temp_1, out);
             }
           if (dim == 1)
-            evaluator.template values<0, false, false>(in, out);
+            evaluator.template hessians<0, false, false>(in, out);
         }
       return false;
     }
diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h
index 5b596874b3..9299bd4921 100644
--- a/include/deal.II/matrix_free/evaluation_template_factory.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.h
@@ -180,6 +180,7 @@ namespace internal
     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 bbd5bf4443..db988b6bc4 100644
--- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -389,6 +389,7 @@ namespace internal
     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,
@@ -398,7 +399,7 @@ namespace internal
       1,
       CellwiseInverseMassMatrixImplTransformFromQPoints<dim,
                                                         VectorizedArrayType>>(
-      fe_degree, fe_degree + 1, n_components, fe_eval, in_array, out_array);
+      fe_degree, n_q_points_1d, n_components, fe_eval, in_array, out_array);
   }
 
 } // end of namespace internal
diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 971b5d9e05..8465ea3415 100644
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -1064,12 +1064,24 @@ namespace MatrixFreeOperators
                                      const VectorizedArrayType *in_array,
                                      VectorizedArrayType *      out_array) const
   {
-    internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
-      dim,
-      VectorizedArrayType>::template run<fe_degree>(n_actual_components,
-                                                    fe_eval,
-                                                    in_array,
-                                                    out_array);
+    const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
+
+    if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
+      internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
+        dim,
+        VectorizedArrayType>::template run<fe_degree,
+                                           fe_degree + 1>(n_actual_components,
+                                                          fe_eval,
+                                                          in_array,
+                                                          out_array);
+    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);
   }