Apply requests from the review

- Remove explicit unrolling of tensors
- Use ArrayView in internal functions instead of passing around AlignedVector
- Make vmult private and inline
- Documentation

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index f999a7eaa60075cd0342160eb5076ff0a28f414f..b0145fa8e472d85664d708c4cc66c930746b8ff2 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -5707,7 +5707,7 @@ namespace internal
     static bool
     run(const unsigned int                          n_desired_components,
         const FEEvaluationData<dim, Number, false> &fe_eval,
-        const AlignedVector<Number> &               inverse_coefficients,
+        const ArrayView<const Number> &             inverse_coefficients,
         const bool                                  dyadic_coefficients,
         const Number *                              in_array,
         Number *                                    out_array)
@@ -5728,7 +5728,7 @@ namespace internal
         {
           if (inverse_coefficients.size() != dofs_per_component)
             AssertDimension(n_desired_components * dofs_per_component,
-                            inverse_coefficients.size())
+                            inverse_coefficients.size());
         }
       else
         {
@@ -5821,8 +5821,9 @@ namespace internal
       return false;
     }
 
+  private:
     template <int n_components>
-    static void
+    static inline void
     vmult(const Number *     inverse_coefficients,
           const Number *     src,
           Number *           dst,

diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h
index 616386cb8b52cc4911bb3d70f43e409957695244..a179402e54712d5e908fbb446fd04714be2c008f 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_template_factory.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.h
@@ -107,7 +107,7 @@ namespace internal
     static void
     apply(const unsigned int                          n_components,
           const FEEvaluationData<dim, Number, false> &fe_eval,
-          const AlignedVector<Number> &               inverse_coefficients,
+          const ArrayView<const Number> &             inverse_coefficients,
           const bool                                  dyadic_coefficients,
           const Number *                              in_array,
           Number *                                    out_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 95b93c6f29cc933b439a095a0cdbb503cd5eccbf..67eba41712baf560a0520af0890922020eca0e1c 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -114,7 +114,7 @@ namespace internal
   CellwiseInverseMassFactory<dim, Number>::apply(
     const unsigned int                          n_components,
     const FEEvaluationData<dim, Number, false> &fe_eval,
-    const AlignedVector<Number> &               inverse_coefficients,
+    const ArrayView<const Number> &             inverse_coefficients,
     const bool                                  dyadic_coefficients,
     const Number *                              in_array,
     Number *                                    out_array)
@@ -122,15 +122,14 @@ namespace internal
     const unsigned int fe_degree = fe_eval.get_shape_info().data[0].fe_degree;
     instantiation_helper_run<
       1,
-      CellwiseInverseMassMatrixImplFlexible<dim, Number>>(
-      fe_degree,
-      fe_degree + 1,
-      n_components,
-      fe_eval,
-      inverse_coefficients,
-      dyadic_coefficients,
-      in_array,
-      out_array);
+      CellwiseInverseMassMatrixImplFlexible<dim, Number>>(fe_degree,
+                                                          fe_degree + 1,
+                                                          n_components,
+                                                          fe_eval,
+                                                          inverse_coefficients,
+                                                          dyadic_coefficients,
+                                                          in_array,
+                                                          out_array);
   }
 
 

diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 480aafb288bf6c4e83b0b41359d530ec7fa4805b..e2c5334900e6c531cd05f30e3caa612a1fa26cc9 100644 (file)
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -611,7 +611,8 @@ namespace MatrixFreeOperators
    * The equation may contain variable coefficients, so the user is required
    * to provide an array for the inverse of the local coefficient (this class
    * provide a helper method 'fill_inverse_JxW_values' to get the inverse of a
-   * constant-coefficient operator).
+   * constant-coefficient operator). The local coefficient can either be scalar
+   * in each component, or dyadic, i.e. couple between components.
    */
   template <int dim,
             int fe_degree,
@@ -642,14 +643,14 @@ namespace MatrixFreeOperators
      * 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.
+     * coefficient is allowed. The coefficients are interpreted as scalar
+     * in each component.
      */
     void
     apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
           const unsigned int                        n_actual_components,
           const VectorizedArrayType *               in_array,
-          VectorizedArrayType *                     out_array,
-          const bool dyadic_coefficients = false) const;
+          VectorizedArrayType *                     out_array) const;
 
     /**
      * Applies the inverse @ref GlossMassMatrix "mass matrix" operation on an input array, using the
@@ -672,11 +673,11 @@ namespace MatrixFreeOperators
      * The second-rank tensor at each quadrature point defines a linear operator
      * on a vector holding the dof components. It is assumed that the passed
      * input and output arrays are of correct size, namely
-     * FEEvaluation::dofs_per_cell long.
+     * FEEvaluation::dofs_per_cell long. The `in_array` and `out_array`
+     * arguments may point to the same memory position.
      * `inverse_dyadic_coefficients` must be dofs_per_component long, and every
      * element must be a second-rank tensor of dimension `n_components`. All
-     * entries should also contain the inverse JxW values. The `in_array` and
-     * `out_array` arguments may point to the same memory position.
+     * entries should also contain the inverse JxW values.
      */
     void
     apply(const AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
@@ -1112,24 +1113,26 @@ namespace MatrixFreeOperators
     apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
           const unsigned int                        n_actual_components,
           const VectorizedArrayType *               in_array,
-          VectorizedArrayType *                     out_array,
-          const bool                                dyadic_coefficients) const
+          VectorizedArrayType *                     out_array) const
   {
     if (fe_degree > -1)
-      internal::CellwiseInverseMassMatrixImplFlexible<
-        dim,
-        VectorizedArrayType>::template run<fe_degree>(n_actual_components,
-                                                      fe_eval,
-                                                      inverse_coefficients,
-                                                      dyadic_coefficients,
-                                                      in_array,
-                                                      out_array);
+      internal::CellwiseInverseMassMatrixImplFlexible<dim,
+                                                      VectorizedArrayType>::
+        template run<fe_degree>(
+          n_actual_components,
+          fe_eval,
+          ArrayView<const VectorizedArrayType>(inverse_coefficients.data(),
+                                               inverse_coefficients.size()),
+          false,
+          in_array,
+          out_array);
     else
       internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
         n_actual_components,
         fe_eval,
-        inverse_coefficients,
-        dyadic_coefficients,
+        ArrayView<const VectorizedArrayType>(inverse_coefficients.data(),
+                                             inverse_coefficients.size()),
+        false,
         in_array,
         out_array);
   }
@@ -1150,24 +1153,29 @@ namespace MatrixFreeOperators
           const VectorizedArrayType *in_array,
           VectorizedArrayType *      out_array) const
   {
-    const unsigned int dofs_per_component = inverse_dyadic_coefficients.size();
-    constexpr unsigned int n_tensor_components = n_components * n_components;
-
-    AlignedVector<VectorizedArrayType> inverse_coefficients(
-      dofs_per_component * n_tensor_components);
-
-    // Flatten the inverse dyadic coefficients into `inverse_coefficients`
-    {
-      auto begin = inverse_coefficients.begin();
-      for (unsigned int q = 0; q < dofs_per_component; ++q)
-        {
-          const auto end = std::next(begin, n_tensor_components);
-          inverse_dyadic_coefficients[q].unroll(begin, end);
-          begin = end;
-        }
-    }
+    const unsigned int unrolled_size =
+      inverse_dyadic_coefficients.size() * (n_components * n_components);
 
-    apply(inverse_coefficients, n_components, in_array, out_array, true);
+    if (fe_degree > -1)
+      internal::CellwiseInverseMassMatrixImplFlexible<dim,
+                                                      VectorizedArrayType>::
+        template run<fe_degree>(n_components,
+                                fe_eval,
+                                ArrayView<const VectorizedArrayType>(
+                                  &inverse_dyadic_coefficients[0][0][0],
+                                  unrolled_size),
+                                true,
+                                in_array,
+                                out_array);
+    else
+      internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
+        n_components,
+        fe_eval,
+        ArrayView<const VectorizedArrayType>(
+          &inverse_dyadic_coefficients[0][0][0], unrolled_size),
+        true,
+        in_array,
+        out_array);
   }