From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Sun, 8 Dec 2019 19:09:19 +0000 (+0100)
Subject: Use inverse shape values for inverse mass operation
X-Git-Tag: v9.2.0-rc1~767^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b8421de607b395cabbdf7294b86309307a9cd60f;p=dealii.git

Use inverse shape values for inverse mass operation
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index e77ae68e11..4ab979c2aa 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -3037,6 +3037,9 @@ namespace internal
 
 
 
+  /**
+   * This struct implements the action of the inverse mass matrix operation
+   */
   template <int dim, int fe_degree, int n_components, typename Number>
   struct CellwiseInverseMassMatrixImpl
   {
@@ -3078,30 +3081,30 @@ namespace internal
           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, false, false>(in, out);
+          evaluator.template hessians<0, true, false>(in, out);
           if (dim > 1)
             {
-              evaluator.template hessians<1, false, false>(out, out);
+              evaluator.template hessians<1, true, false>(out, out);
 
               if (dim == 3)
                 {
-                  evaluator.template hessians<2, false, false>(out, out);
+                  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, true, false>(out, out);
+                  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];
 
-              evaluator.template hessians<1, true, false>(out, out);
+              evaluator.template hessians<1, false, false>(out, out);
             }
           else
             {
               for (unsigned int q = 0; q < dofs_per_component; ++q)
                 out[q] *= inv_coefficient[q];
             }
-          evaluator.template hessians<0, true, false>(out, out);
+          evaluator.template hessians<0, false, false>(out, out);
 
           inv_coefficient += shift_coefficient;
         }
@@ -3130,17 +3133,17 @@ namespace internal
 
           if (dim == 3)
             {
-              evaluator.template hessians<2, true, false>(in, out);
-              evaluator.template hessians<1, true, false>(out, out);
-              evaluator.template hessians<0, true, false>(out, out);
+              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, true, false>(in, out);
-              evaluator.template hessians<0, true, false>(out, out);
+              evaluator.template hessians<1, false, false>(in, out);
+              evaluator.template hessians<0, false, false>(out, out);
             }
           if (dim == 1)
-            evaluator.template hessians<0, true, false>(in, out);
+            evaluator.template hessians<0, false, false>(in, out);
         }
     }
   };
diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 2920c3837b..abdc4cf9a2 100644
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -709,11 +709,6 @@ namespace MatrixFreeOperators
                            Number,
                            false,
                            VectorizedArrayType> &fe_eval;
-
-    /**
-     * A structure to hold inverse shape functions
-     */
-    AlignedVector<VectorizedArrayType> inverse_shape;
   };
 
 
@@ -958,26 +953,7 @@ namespace MatrixFreeOperators
                              false,
                              VectorizedArrayType> &fe_eval)
     : fe_eval(fe_eval)
-  {
-    FullMatrix<double> shapes_1d(fe_degree + 1, fe_degree + 1);
-    for (unsigned int i = 0, c = 0; i < shapes_1d.m(); ++i)
-      for (unsigned int j = 0; j < shapes_1d.n(); ++j, ++c)
-        shapes_1d(i, j) = fe_eval.get_shape_info().shape_values[c][0];
-    shapes_1d.gauss_jordan();
-    const unsigned int stride = (fe_degree + 2) / 2;
-    inverse_shape.resize(stride * (fe_degree + 1));
-    for (unsigned int i = 0; i < stride; ++i)
-      for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
-        {
-          inverse_shape[i * stride + q] =
-            0.5 * (shapes_1d(i, q) + shapes_1d(i, fe_degree - q));
-          inverse_shape[(fe_degree - i) * stride + q] =
-            0.5 * (shapes_1d(i, q) - shapes_1d(i, fe_degree - q));
-        }
-    if (fe_degree % 2 == 0)
-      for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
-        inverse_shape[fe_degree / 2 * stride + q] = shapes_1d(fe_degree / 2, q);
-  }
+  {}
 
 
 
@@ -1030,15 +1006,15 @@ namespace MatrixFreeOperators
           const VectorizedArrayType *               in_array,
           VectorizedArrayType *                     out_array) const
   {
-    internal::CellwiseInverseMassMatrixImpl<
-      dim,
-      fe_degree,
-      n_components,
-      VectorizedArrayType>::apply(inverse_shape,
-                                  inverse_coefficients,
-                                  n_actual_components,
-                                  in_array,
-                                  out_array);
+    internal::CellwiseInverseMassMatrixImpl<dim,
+                                            fe_degree,
+                                            n_components,
+                                            VectorizedArrayType>::
+      apply(fe_eval.get_shape_info().inverse_shape_values_eo,
+            inverse_coefficients,
+            n_actual_components,
+            in_array,
+            out_array);
   }
 
 
@@ -1062,10 +1038,11 @@ namespace MatrixFreeOperators
                                             fe_degree,
                                             n_components,
                                             VectorizedArrayType>::
-      transform_from_q_points_to_basis(inverse_shape,
-                                       n_actual_components,
-                                       in_array,
-                                       out_array);
+      transform_from_q_points_to_basis(
+        fe_eval.get_shape_info().inverse_shape_values_eo,
+        n_actual_components,
+        in_array,
+        out_array);
   }