Implement simpler CellwiseInverseMassMatrix::apply() function

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 4ab979c2aafa3cbb966bc20a3557c033df72a757..79623b3e8ad805d799d1b50f7c884644b3de3ef6 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -3043,6 +3043,58 @@ namespace internal
   template <int dim, int fe_degree, int n_components, typename Number>
   struct CellwiseInverseMassMatrixImpl
   {
+    template <typename FEEvaluationType>
+    static void
+    apply(const FEEvaluationType &fe_eval,
+          const Number *          in_array,
+          Number *                out_array)
+    {
+      constexpr unsigned int dofs_per_component =
+        Utilities::pow(fe_degree + 1, dim);
+
+      Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
+      Assert(fe_eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric,
+             ExcNotImplemented());
+
+      internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
+                                       dim,
+                                       fe_degree + 1,
+                                       fe_degree + 1,
+                                       Number>
+        evaluator(AlignedVector<Number>(),
+                  AlignedVector<Number>(),
+                  fe_eval.get_shape_info().inverse_shape_values_eo);
+
+      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 hessians<0, true, false>(in, out);
+          if (dim > 1)
+            evaluator.template hessians<1, true, false>(out, out);
+          if (dim > 2)
+            evaluator.template hessians<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 hessians<2, false, false>(out, out);
+          if (dim > 1)
+            evaluator.template hessians<1, false, false>(out, out);
+          evaluator.template hessians<0, false, false>(out, out);
+        }
+    }
+
     static void
     apply(const AlignedVector<Number> &inverse_shape,
           const AlignedVector<Number> &inverse_coefficients,
@@ -3083,27 +3135,17 @@ namespace internal
           // assume symmetries that do not exist in the inverse shapes
           evaluator.template hessians<0, true, false>(in, out);
           if (dim > 1)
-            {
-              evaluator.template hessians<1, true, false>(out, out);
+            evaluator.template hessians<1, true, false>(out, out);
+          if (dim > 2)
+            evaluator.template hessians<2, true, false>(out, out);
 
-              if (dim == 3)
-                {
-                  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, false, false>(out, out);
-                }
-              else if (dim == 2)
-                for (unsigned int q = 0; q < dofs_per_component; ++q)
-                  out[q] *= inv_coefficient[q];
+          for (unsigned int q = 0; q < dofs_per_component; ++q)
+            out[q] *= inv_coefficient[q];
 
-              evaluator.template hessians<1, false, false>(out, out);
-            }
-          else
-            {
-              for (unsigned int q = 0; q < dofs_per_component; ++q)
-                out[q] *= inv_coefficient[q];
-            }
+          if (dim > 2)
+            evaluator.template hessians<2, false, false>(out, out);
+          if (dim > 1)
+            evaluator.template hessians<1, false, false>(out, out);
           evaluator.template hessians<0, false, false>(out, out);
 
           inv_coefficient += shift_coefficient;

diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index abdc4cf9a2c5332aba59971d76c609e200c6c7a6..9d71deab6d537b2774967058116a546308e79d08 100644 (file)
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -642,7 +642,7 @@ namespace MatrixFreeOperators
     /**
      * Applies the inverse mass matrix operation on an input array. It is
      * assumed that the passed input and output arrays are of correct size,
-     * namely FEEval::dofs_per_cell * n_components long. The inverse of the
+     * namely FEEvaluation::dofs_per_cell long. The inverse of the
      * local coefficient (also containing the inverse JxW values) must be
      * passed as first argument. Passing more than one component in the
      * coefficient is allowed.
@@ -653,6 +653,21 @@ namespace MatrixFreeOperators
           const VectorizedArrayType *               in_array,
           VectorizedArrayType *                     out_array) const;
 
+    /**
+     * Applies the inverse mass matrix operation on an input array, using the
+     * inverse of the JxW values provided by the `fe_eval` argument passed to
+     * the constructor of this class. Note that the user code must call
+     * FEEvaluation::reinit() on the underlying evaluator to make the
+     * FEEvaluationBase::JxW() method return the information of the correct
+     * cell. It is assumed that the pointers of the input and output arrays
+     * are valid over the length FEEvaluation::dofs_per_cell, which is the
+     * number of entries processed by this function. The `in_array` and
+     * `out_array` arguments may point to the same memory position.
+     */
+    void
+    apply(const VectorizedArrayType *in_array,
+          VectorizedArrayType *      out_array) const;
+
     /**
      * This operation performs a projection from the data given in quadrature
      * points to the actual basis underlying this object. This projection can
@@ -990,6 +1005,29 @@ namespace MatrixFreeOperators
 
 
 
+  template <int dim,
+            int fe_degree,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType>
+  inline void
+  CellwiseInverseMassMatrix<
+    dim,
+    fe_degree,
+    n_components,
+    Number,
+    VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
+                                VectorizedArrayType *      out_array) const
+  {
+    internal::CellwiseInverseMassMatrixImpl<
+      dim,
+      fe_degree,
+      n_components,
+      VectorizedArrayType>::apply(fe_eval, in_array, out_array);
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_components,