Implement InverseCellwiseMassMatrix with tensor var. coeff.

based on the `Flexible` implementation with its own new interface (vector of rank-2 tensors), see #14843

- Unify the flexible implementation as well as fill inverse_JxW_values() for all `fe_degree`s.
- Add boolean template parameter `dyadic_coefficients` to the flexible implementation.
- Flexible implementation no longer takes the shape values, but fe_eval.
- Matrix-vector product for dyadic coefficients implemented in vmult.

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index da52686722be4ab040f02c29571ee2277d3d96ed..f999a7eaa60075cd0342160eb5076ff0a28f414f 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -5705,24 +5705,42 @@ namespace internal
   {
     template <int fe_degree, int = 0>
     static bool
-    run(const unsigned int           n_desired_components,
-        const AlignedVector<Number> &inverse_shape,
-        const AlignedVector<Number> &inverse_coefficients,
-        const Number *               in_array,
-        Number *                     out_array,
-        std::enable_if_t<fe_degree != -1> * = nullptr)
+    run(const unsigned int                          n_desired_components,
+        const FEEvaluationData<dim, Number, false> &fe_eval,
+        const AlignedVector<Number> &               inverse_coefficients,
+        const bool                                  dyadic_coefficients,
+        const Number *                              in_array,
+        Number *                                    out_array)
     {
-      constexpr unsigned int dofs_per_component =
-        Utilities::pow(fe_degree + 1, dim);
+      const unsigned int given_degree =
+        (fe_degree > -1) ? fe_degree :
+                           fe_eval.get_shape_info().data.front().fe_degree;
+
+      const unsigned int dofs_per_component =
+        Utilities::pow(given_degree + 1, dim);
+
       Assert(inverse_coefficients.size() > 0 &&
                inverse_coefficients.size() % dofs_per_component == 0,
              ExcMessage(
                "Expected diagonal to be a multiple of scalar dof per cells"));
-      if (inverse_coefficients.size() != dofs_per_component)
-        AssertDimension(n_desired_components * dofs_per_component,
-                        inverse_coefficients.size());
+
+      if (!dyadic_coefficients)
+        {
+          if (inverse_coefficients.size() != dofs_per_component)
+            AssertDimension(n_desired_components * dofs_per_component,
+                            inverse_coefficients.size())
+        }
+      else
+        {
+          AssertDimension(n_desired_components * n_desired_components *
+                            dofs_per_component,
+                          inverse_coefficients.size());
+        }
 
       Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
+      Assert(fe_eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
+             ExcNotImplemented());
 
       EvaluatorTensorProduct<evaluate_evenodd,
                              dim,
@@ -5731,53 +5749,106 @@ namespace internal
                              Number>
         evaluator(AlignedVector<Number>(),
                   AlignedVector<Number>(),
-                  inverse_shape);
+                  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
+                  given_degree + 1,
+                  given_degree + 1);
+
+      const Number *in  = in_array;
+      Number *      out = out_array;
+
+      const Number *inv_coefficient = inverse_coefficients.data();
 
       const unsigned int shift_coefficient =
         inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
                                                            0;
-      const Number *inv_coefficient = inverse_coefficients.data();
-      for (unsigned int d = 0; d < n_desired_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)
-            out[q] *= inv_coefficient[q];
+      const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
+      const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
 
-          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);
+      for (unsigned int d = 0; d < n_comp_outer; ++d)
+        {
+          for (unsigned int di = 0; di < n_comp_inner; ++di)
+            {
+              const Number *in_  = in + di * dofs_per_component;
+              Number *      out_ = out + di * dofs_per_component;
+              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_);
+            }
+          if (dyadic_coefficients)
+            {
+              const auto n_coeff_components =
+                n_desired_components * n_desired_components;
+              if (n_desired_components == dim)
+                {
+                  for (unsigned int q = 0; q < dofs_per_component; ++q)
+                    vmult<dim>(&inv_coefficient[q * n_coeff_components],
+                               &in[q],
+                               &out[q],
+                               dofs_per_component);
+                }
+              else
+                {
+                  for (unsigned int q = 0; q < dofs_per_component; ++q)
+                    vmult<-1>(&inv_coefficient[q * n_coeff_components],
+                              &in[q],
+                              &out[q],
+                              dofs_per_component,
+                              n_desired_components);
+                }
+            }
+          else
+            for (unsigned int q = 0; q < dofs_per_component; ++q)
+              out[q] *= inv_coefficient[q];
+
+          for (unsigned int di = 0; di < n_comp_inner; ++di)
+            {
+              Number *out_ = out + di * 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_);
+            }
 
+          in += dofs_per_component;
+          out += dofs_per_component;
           inv_coefficient += shift_coefficient;
         }
+
       return false;
     }
 
-    /**
-     * Version for degree = -1
-     */
-    template <int fe_degree, int = 0>
-    static bool
-    run(const unsigned int,
-        const AlignedVector<Number> &,
-        const AlignedVector<Number> &,
-        const Number *,
-        Number *,
-        std::enable_if_t<fe_degree == -1> * = nullptr)
+    template <int n_components>
+    static void
+    vmult(const Number *     inverse_coefficients,
+          const Number *     src,
+          Number *           dst,
+          const unsigned int dofs_per_component,
+          const unsigned int n_given_components = 0)
     {
-      static_assert(fe_degree == -1, "Only usable for degree -1");
-      Assert(false, ExcNotImplemented());
-      return false;
+      const unsigned int n_desired_components =
+        (n_components > -1) ? n_components : n_given_components;
+
+      std::array<Number, dim + 2> tmp;
+      Assert(n_desired_components <= dim + 2,
+             ExcMessage(
+               "Number of components larger than dim+2 not supported."));
+
+      for (unsigned int d = 0; d < n_desired_components; ++d)
+        tmp[d] = src[d * dofs_per_component];
+
+      for (unsigned int d1 = 0; d1 < n_desired_components; ++d1)
+        {
+          const Number *inv_coeff_row =
+            &inverse_coefficients[d1 * n_desired_components];
+          Number sum = inv_coeff_row[0] * tmp[0];
+          for (unsigned int d2 = 1; d2 < n_desired_components; ++d2)
+            sum += inv_coeff_row[d2] * tmp[d2];
+          dst[d1 * dofs_per_component] = sum;
+        }
     }
   };
 

diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h
index 64aaae2b4a251f576455b1a03bf93fde8cb85bdc..616386cb8b52cc4911bb3d70f43e409957695244 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_template_factory.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.h
@@ -105,12 +105,12 @@ namespace internal
           Number *                                    out_array);
 
     static void
-    apply(const unsigned int           n_components,
-          const unsigned int           fe_degree,
-          const AlignedVector<Number> &inverse_shape,
-          const AlignedVector<Number> &inverse_coefficients,
-          const Number *               in_array,
-          Number *                     out_array);
+    apply(const unsigned int                          n_components,
+          const FEEvaluationData<dim, Number, false> &fe_eval,
+          const AlignedVector<Number> &               inverse_coefficients,
+          const bool                                  dyadic_coefficients,
+          const Number *                              in_array,
+          Number *                                    out_array);
 
     static void
     transform_from_q_points_to_basis(

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 f02a7383e458a3d200e4fbbfbe41f07d419214e9..95b93c6f29cc933b439a095a0cdbb503cd5eccbf 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h
+++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -112,22 +112,25 @@ namespace internal
   template <int dim, typename Number>
   void
   CellwiseInverseMassFactory<dim, Number>::apply(
-    const unsigned int           n_components,
-    const unsigned int           fe_degree,
-    const AlignedVector<Number> &inverse_shape,
-    const AlignedVector<Number> &inverse_coefficients,
-    const Number *               in_array,
-    Number *                     out_array)
+    const unsigned int                          n_components,
+    const FEEvaluationData<dim, Number, false> &fe_eval,
+    const AlignedVector<Number> &               inverse_coefficients,
+    const bool                                  dyadic_coefficients,
+    const Number *                              in_array,
+    Number *                                    out_array)
   {
+    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,
-                                                          inverse_shape,
-                                                          inverse_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 f6bfcfe0d9bf63a56b982f6217f73c807f2c8239..480aafb288bf6c4e83b0b41359d530ec7fa4805b 100644 (file)
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -648,7 +648,8 @@ namespace MatrixFreeOperators
     apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
           const unsigned int                        n_actual_components,
           const VectorizedArrayType *               in_array,
-          VectorizedArrayType *                     out_array) const;
+          VectorizedArrayType *                     out_array,
+          const bool dyadic_coefficients = false) const;
 
     /**
      * Applies the inverse @ref GlossMassMatrix "mass matrix" operation on an input array, using the
@@ -665,6 +666,24 @@ namespace MatrixFreeOperators
     apply(const VectorizedArrayType *in_array,
           VectorizedArrayType *      out_array) const;
 
+    /**
+     * This operation applies the inverse @ref GlossMassMatrix "mass matrix"
+     * operation on an input array with local dyadic-valued coefficients.
+     * 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.
+     * `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.
+     */
+    void
+    apply(const AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
+            &                        inverse_dyadic_coefficients,
+          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
@@ -1033,8 +1052,12 @@ namespace MatrixFreeOperators
     fill_inverse_JxW_values(
       AlignedVector<VectorizedArrayType> &inverse_jxw) const
   {
-    constexpr unsigned int dofs_per_component_on_cell =
-      Utilities::pow(fe_degree + 1, dim);
+    const unsigned int dofs_per_component_on_cell =
+      (fe_degree > -1) ?
+        Utilities::pow(fe_degree + 1, dim) :
+        Utilities::pow(fe_eval.get_shape_info().data.front().fe_degree + 1,
+                       dim - 1);
+
     Assert(inverse_jxw.size() > 0 &&
              inverse_jxw.size() % dofs_per_component_on_cell == 0,
            ExcMessage(
@@ -1089,29 +1112,64 @@ namespace MatrixFreeOperators
     apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
           const unsigned int                        n_actual_components,
           const VectorizedArrayType *               in_array,
-          VectorizedArrayType *                     out_array) const
+          VectorizedArrayType *                     out_array,
+          const bool                                dyadic_coefficients) const
   {
-    const unsigned int given_degree =
-      fe_eval.get_shape_info().data[0].fe_degree;
     if (fe_degree > -1)
-      internal::CellwiseInverseMassMatrixImplFlexible<dim,
-                                                      VectorizedArrayType>::
-        template run<fe_degree>(
-          n_actual_components,
-          fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
-          inverse_coefficients,
-          in_array,
-          out_array);
+      internal::CellwiseInverseMassMatrixImplFlexible<
+        dim,
+        VectorizedArrayType>::template run<fe_degree>(n_actual_components,
+                                                      fe_eval,
+                                                      inverse_coefficients,
+                                                      dyadic_coefficients,
+                                                      in_array,
+                                                      out_array);
     else
       internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
         n_actual_components,
-        given_degree,
-        fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
+        fe_eval,
         inverse_coefficients,
+        dyadic_coefficients,
         in_array,
         out_array);
   }
 
+  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 AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
+            &                        inverse_dyadic_coefficients,
+          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;
+        }
+    }
+
+    apply(inverse_coefficients, n_components, in_array, out_array, true);
+  }
+
 
 
   template <int dim,