From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 4 Jun 2021 17:22:48 +0000 (+0200)
Subject: CellwiseInverseMassMatrix::transform_from_q_points_to_basis: for degree = -1 and... 
X-Git-Tag: v9.4.0-rc1~1277^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=80ca15ae1fe0bc0e5871b0dfa85f3f58a1d0892b;p=dealii.git

CellwiseInverseMassMatrix::transform_from_q_points_to_basis: for degree = -1 and degree+1!=q_points_1d
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index b06955fa61..4901a4a2df 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -4351,12 +4351,18 @@ namespace internal
   {
     template <int fe_degree, int = 0>
     static bool
-    run(const unsigned int           n_desired_components,
-        const AlignedVector<Number> &inverse_shape,
-        const Number *               in_array,
-        Number *                     out_array,
+    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)
     {
+      const AlignedVector<Number> &inverse_shape =
+        fe_eval.get_shape_info().data.front().inverse_shape_values_eo;
+
       constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
       internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
                                        dim,
@@ -4391,14 +4397,54 @@ namespace internal
 
     template <int fe_degree, int = 0>
     static bool
-    run(const unsigned int,
-        const AlignedVector<Number> &,
-        const Number *,
-        Number *,
+    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");
-      Assert(false, ExcNotImplemented());
+
+      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);
+
+      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;
+
+          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);
+            }
+          if (dim == 2)
+            {
+              evaluator.template values<1, false, false>(in, temp_1);
+              evaluator.template values<0, false, false>(temp_1, out);
+            }
+          if (dim == 1)
+            evaluator.template values<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 2523f20d32..5b596874b3 100644
--- a/include/deal.II/matrix_free/evaluation_template_factory.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.h
@@ -178,11 +178,12 @@ namespace internal
 
     static void
     transform_from_q_points_to_basis(
-      const unsigned int                        n_components,
-      const unsigned int                        fe_degree,
-      const AlignedVector<VectorizedArrayType> &inverse_shape,
-      const VectorizedArrayType *               in_array,
-      VectorizedArrayType *                     out_array);
+      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);
   };
 
 } // end of namespace internal
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 6256b0bca5..bbd5bf4443 100644
--- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -387,22 +387,18 @@ namespace internal
   void
   CellwiseInverseMassFactory<dim, Number, VectorizedArrayType>::
     transform_from_q_points_to_basis(
-      const unsigned int                        n_components,
-      const unsigned int                        fe_degree,
-      const AlignedVector<VectorizedArrayType> &inverse_shape,
-      const VectorizedArrayType *               in_array,
-      VectorizedArrayType *                     out_array)
+      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)
   {
     instantiation_helper_run<
       1,
       CellwiseInverseMassMatrixImplTransformFromQPoints<dim,
                                                         VectorizedArrayType>>(
-      fe_degree,
-      fe_degree + 1,
-      n_components,
-      inverse_shape,
-      in_array,
-      out_array);
+      fe_degree, fe_degree + 1, 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 0a36f49c3b..971b5d9e05 100644
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -1067,9 +1067,7 @@ namespace MatrixFreeOperators
     internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
       dim,
       VectorizedArrayType>::template run<fe_degree>(n_actual_components,
-                                                    fe_eval.get_shape_info()
-                                                      .data.front()
-                                                      .inverse_shape_values_eo,
+                                                    fe_eval,
                                                     in_array,
                                                     out_array);
   }