Change

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index ffd07d08a80030651e08940165b7b47cb6dcdf30..5c4010eea4efbec5350457eb8160fc05197fe123 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -648,10 +648,11 @@ namespace internal
    * This class allows for dimension-independent application of the operation,
    * implemented by template recursion. It has been tested up to 6D.
    */
-  template <EvaluatorVariant variant,
-            int              dim,
-            int              basis_size_1,
-            int              basis_size_2,
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            int               dim,
+            int               basis_size_1,
+            int               basis_size_2,
             typename Number,
             typename Number2>
   struct FEEvaluationImplBasisChange
@@ -696,6 +697,8 @@ namespace internal
         basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
         ExcMessage("The second dimension must not be smaller than the first"));
 
+      Assert(quantity == EvaluatorQuantity::value, ExcInternalError());
+
       // we do recursion until dim==1 or dim==2 and we have
       // basis_size_1==basis_size_2. The latter optimization increases
       // optimization possibilities for the compiler but does only work for
@@ -739,6 +742,7 @@ namespace internal
             for (unsigned int q = np_1; q != 0; --q)
               FEEvaluationImplBasisChange<
                 variant,
+                quantity,
                 next_dim,
                 basis_size_1,
                 basis_size_2,
@@ -821,6 +825,10 @@ namespace internal
                "Input and output cannot alias with each other when "
                "adding the result of the basis change to existing data"));
 
+      Assert(quantity == EvaluatorQuantity::value ||
+               quantity == EvaluatorQuantity::hessian,
+             ExcInternalError());
+
       constexpr int next_dim =
         (dim > 2 ||
          ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
@@ -833,8 +841,8 @@ namespace internal
                              Number,
                              Number2>
                          eval_val(transformation_matrix,
-                 AlignedVector<Number2>(),
-                 AlignedVector<Number2>(),
+                 transformation_matrix,
+                 transformation_matrix,
                  basis_size_1_variable,
                  basis_size_2_variable);
       const unsigned int np_1 =
@@ -850,27 +858,65 @@ namespace internal
         {
           if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
             {
-              eval_val.template values<1, false, false>(values_in, values_in);
+              if (quantity == EvaluatorQuantity::value)
+                eval_val.template values<1, false, false>(values_in, values_in);
+              else
+                eval_val.template hessians<1, false, false>(values_in,
+                                                            values_in);
+
               if (add_into_result)
-                eval_val.template values<0, false, true>(values_in, values_out);
+                {
+                  if (quantity == EvaluatorQuantity::value)
+                    eval_val.template values<0, false, true>(values_in,
+                                                             values_out);
+                  else
+                    eval_val.template hessians<0, false, true>(values_in,
+                                                               values_out);
+                }
               else
-                eval_val.template values<0, false, false>(values_in,
-                                                          values_out);
+                {
+                  if (quantity == EvaluatorQuantity::value)
+                    eval_val.template values<0, false, false>(values_in,
+                                                              values_out);
+                  else
+                    eval_val.template hessians<0, false, false>(values_in,
+                                                                values_out);
+                }
             }
           else
             {
               if (dim == 1 && add_into_result)
-                eval_val.template values<0, false, true>(values_in, values_out);
+                {
+                  if (quantity == EvaluatorQuantity::value)
+                    eval_val.template values<0, false, true>(values_in,
+                                                             values_out);
+                  else
+                    eval_val.template hessians<0, false, true>(values_in,
+                                                               values_out);
+                }
               else if (dim == 1)
-                eval_val.template values<0, false, false>(values_in,
-                                                          values_out);
+                {
+                  if (quantity == EvaluatorQuantity::value)
+                    eval_val.template values<0, false, false>(values_in,
+                                                              values_out);
+                  else
+                    eval_val.template hessians<0, false, false>(values_in,
+                                                                values_out);
+                }
               else
-                eval_val.template values<dim - 1, false, false>(values_in,
-                                                                values_in);
+                {
+                  if (quantity == EvaluatorQuantity::value)
+                    eval_val.template values<dim - 1, false, false>(values_in,
+                                                                    values_in);
+                  else
+                    eval_val.template hessians<dim - 1, false, false>(
+                      values_in, values_in);
+                }
             }
           if (next_dim < dim)
             for (unsigned int q = 0; q < np_1; ++q)
               FEEvaluationImplBasisChange<variant,
+                                          quantity,
                                           next_dim,
                                           basis_size_1,
                                           basis_size_2,
@@ -891,107 +937,6 @@ namespace internal
         }
     }
 
-    /**
-     * Same as above but with hessians.
-     */
-#ifndef DEBUG
-    DEAL_II_ALWAYS_INLINE
-#endif
-    static void
-    do_backward_hessians(
-      const unsigned int                    n_components,
-      const dealii::AlignedVector<Number2> &transformation_matrix,
-      const bool                            add_into_result,
-      Number *                              values_in,
-      Number *                              values_out,
-      const unsigned int                    basis_size_1_variable =
-        dealii::numbers::invalid_unsigned_int,
-      const unsigned int basis_size_2_variable =
-        dealii::numbers::invalid_unsigned_int)
-    {
-      Assert(basis_size_1 != 0 ||
-               basis_size_1_variable <= basis_size_2_variable,
-             dealii::ExcMessage(
-               "The second dimension must not be smaller than the first"));
-      Assert(add_into_result == false || values_in != values_out,
-             dealii::ExcMessage(
-               "Input and output cannot alias with each other when "
-               "adding the result of the basis change to existing data"));
-
-      constexpr int next_dim =
-        (dim > 2 ||
-         ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
-          dim - 1 :
-          dim;
-      dealii::internal::EvaluatorTensorProduct<
-        variant,
-        dim,
-        basis_size_1,
-        (basis_size_1 == 0 ? 0 : basis_size_2),
-        Number,
-        Number2>
-                         eval_val(dealii::AlignedVector<Number2>(),
-                 dealii::AlignedVector<Number2>(),
-                 transformation_matrix,
-                 basis_size_1_variable,
-                 basis_size_2_variable);
-      const unsigned int np_1 =
-        basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
-      const unsigned int np_2 =
-        basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
-      Assert(np_1 > 0 && np_1 != dealii::numbers::invalid_unsigned_int,
-             dealii::ExcMessage("Cannot transform with 0-point basis"));
-      Assert(np_2 > 0 && np_2 != dealii::numbers::invalid_unsigned_int,
-             dealii::ExcMessage("Cannot transform with 0-point basis"));
-
-      for (unsigned int c = 0; c < n_components; ++c)
-        {
-          if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
-            {
-              eval_val.template hessians<1, false, false>(values_in, values_in);
-              if (add_into_result)
-                eval_val.template hessians<0, false, true>(values_in,
-                                                           values_out);
-              else
-                eval_val.template hessians<0, false, false>(values_in,
-                                                            values_out);
-            }
-          else
-            {
-              if (dim == 1 && add_into_result)
-                eval_val.template hessians<0, false, true>(values_in,
-                                                           values_out);
-              else if (dim == 1)
-                eval_val.template hessians<0, false, false>(values_in,
-                                                            values_out);
-              else
-                eval_val.template hessians<dim - 1, false, false>(values_in,
-                                                                  values_in);
-            }
-          if (next_dim < dim)
-            for (unsigned int q = 0; q < np_1; ++q)
-              FEEvaluationImplBasisChange<variant,
-                                          next_dim,
-                                          basis_size_1,
-                                          basis_size_2,
-                                          Number,
-                                          Number2>::
-                do_backward_hessians(
-                  1,
-                  transformation_matrix,
-                  add_into_result,
-                  values_in +
-                    q * dealii::Utilities::fixed_power<next_dim>(np_2),
-                  values_out +
-                    q * dealii::Utilities::fixed_power<next_dim>(np_1),
-                  basis_size_1_variable,
-                  basis_size_2_variable);
-
-          values_in += dealii::Utilities::fixed_power<dim>(np_2);
-          values_out += dealii::Utilities::fixed_power<dim>(np_1);
-        }
-    }
-
     /**
      * This operation applies a mass-matrix-like operation, consisting of a
      * do_forward() operation, multiplication by the coefficients in the
@@ -1033,6 +978,7 @@ namespace internal
       for (unsigned int q = basis_size_1; q != 0; --q)
         FEEvaluationImplBasisChange<
           variant,
+          EvaluatorQuantity::value,
           next_dim,
           basis_size_1,
           basis_size_2,
@@ -1073,6 +1019,7 @@ namespace internal
       for (unsigned int q = 0; q < basis_size_1; ++q)
         FEEvaluationImplBasisChange<
           variant,
+          EvaluatorQuantity::value,
           next_dim,
           basis_size_1,
           basis_size_2,
@@ -1330,6 +1277,7 @@ namespace internal
       {
         FEEvaluationImplBasisChange<
           evaluate_evenodd,
+          EvaluatorQuantity::value,
           dim,
           (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
           n_q_points_1d,
@@ -1405,6 +1353,7 @@ namespace internal
         // transform back to the original space
         FEEvaluationImplBasisChange<
           evaluate_evenodd,
+          EvaluatorQuantity::value,
           dim,
           (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
           n_q_points_1d,

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index 5d62238c8310b5cced8c7837029e2a1b01219288..278f486a8d69ea63be25b4d428f65f4debedb122 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -67,6 +67,27 @@ namespace internal
 
 
 
+  /**
+   * Determine which quantity should be computed via the tensor product kernels.
+   */
+  enum class EvaluatorQuantity
+  {
+    /**
+     * Evaluate/integrate by shape functions.
+     */
+    value,
+    /**
+     * Evaluate/integrate by gradients of the shape functions.
+     */
+    gradient,
+    /**
+     * Evaluate/integrate by hessians of the shape functions.
+     */
+    hessian
+  };
+
+
+
   /**
    * Generic evaluator framework that valuates the given shape data in general
    * dimensions using the tensor product form. Depending on the particular

diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc
index 153f0d28b694724d14ae24e6df57535c3e139976..f68c62a40bff0eb98c9fc5c9d5d1337c542b15a6 100644 (file)
--- a/source/multigrid/mg_transfer_matrix_free.cc
+++ b/source/multigrid/mg_transfer_matrix_free.cc
@@ -458,19 +458,22 @@ MGTransferMatrixFree<dim, Number>::do_prolongate_add(
           // must go through the components backwards because we want to write
           // the output to the same array as the input
           for (int c = n_components - 1; c >= 0; --c)
-            internal::FEEvaluationImplBasisChange<internal::evaluate_general,
-                                                  dim,
-                                                  degree + 1,
-                                                  2 * degree + 1,
-                                                  VectorizedArray<Number>,
-                                                  VectorizedArray<Number>>::
-              do_forward(1,
-                         prolongation_matrix_1d,
-                         evaluation_data.begin() +
-                           c * Utilities::fixed_power<dim>(degree_size),
-                         evaluation_data.begin() + c * n_scalar_cell_dofs,
-                         fe_degree + 1,
-                         2 * fe_degree + 1);
+            internal::FEEvaluationImplBasisChange<
+              internal::evaluate_general,
+              internal::EvaluatorQuantity::value,
+              dim,
+              degree + 1,
+              2 * degree + 1,
+              VectorizedArray<Number>,
+              VectorizedArray<Number>>::do_forward(1,
+                                                   prolongation_matrix_1d,
+                                                   evaluation_data.begin() +
+                                                     c * Utilities::fixed_power<
+                                                           dim>(degree_size),
+                                                   evaluation_data.begin() +
+                                                     c * n_scalar_cell_dofs,
+                                                   fe_degree + 1,
+                                                   2 * fe_degree + 1);
           weight_dofs_on_child<dim, degree, Number>(
             &weights_on_refined[to_level - 1][(cell / vec_size) * three_to_dim],
             n_components,
@@ -480,19 +483,22 @@ MGTransferMatrixFree<dim, Number>::do_prolongate_add(
       else
         {
           for (int c = n_components - 1; c >= 0; --c)
-            internal::FEEvaluationImplBasisChange<internal::evaluate_general,
-                                                  dim,
-                                                  degree + 1,
-                                                  2 * degree + 2,
-                                                  VectorizedArray<Number>,
-                                                  VectorizedArray<Number>>::
-              do_forward(1,
-                         prolongation_matrix_1d,
-                         evaluation_data.begin() +
-                           c * Utilities::fixed_power<dim>(degree_size),
-                         evaluation_data.begin() + c * n_scalar_cell_dofs,
-                         fe_degree + 1,
-                         2 * fe_degree + 2);
+            internal::FEEvaluationImplBasisChange<
+              internal::evaluate_general,
+              internal::EvaluatorQuantity::value,
+              dim,
+              degree + 1,
+              2 * degree + 2,
+              VectorizedArray<Number>,
+              VectorizedArray<Number>>::do_forward(1,
+                                                   prolongation_matrix_1d,
+                                                   evaluation_data.begin() +
+                                                     c * Utilities::fixed_power<
+                                                           dim>(degree_size),
+                                                   evaluation_data.begin() +
+                                                     c * n_scalar_cell_dofs,
+                                                   fe_degree + 1,
+                                                   2 * fe_degree + 2);
         }
 
       // write into dst vector
@@ -556,38 +562,46 @@ MGTransferMatrixFree<dim, Number>::do_restrict_add(
             fe_degree,
             evaluation_data.data());
           for (unsigned int c = 0; c < n_components; ++c)
-            internal::FEEvaluationImplBasisChange<internal::evaluate_general,
-                                                  dim,
-                                                  degree + 1,
-                                                  2 * degree + 1,
-                                                  VectorizedArray<Number>,
-                                                  VectorizedArray<Number>>::
-              do_backward(1,
-                          prolongation_matrix_1d,
-                          false,
-                          evaluation_data.begin() + c * n_scalar_cell_dofs,
-                          evaluation_data.begin() +
-                            c * Utilities::fixed_power<dim>(degree_size),
-                          fe_degree + 1,
-                          2 * fe_degree + 1);
+            internal::FEEvaluationImplBasisChange<
+              internal::evaluate_general,
+              internal::EvaluatorQuantity::value,
+              dim,
+              degree + 1,
+              2 * degree + 1,
+              VectorizedArray<Number>,
+              VectorizedArray<Number>>::do_backward(1,
+                                                    prolongation_matrix_1d,
+                                                    false,
+                                                    evaluation_data.begin() +
+                                                      c * n_scalar_cell_dofs,
+                                                    evaluation_data.begin() +
+                                                      c *
+                                                        Utilities::fixed_power<
+                                                          dim>(degree_size),
+                                                    fe_degree + 1,
+                                                    2 * fe_degree + 1);
         }
       else
         {
           for (unsigned int c = 0; c < n_components; ++c)
-            internal::FEEvaluationImplBasisChange<internal::evaluate_general,
-                                                  dim,
-                                                  degree + 1,
-                                                  2 * degree + 2,
-                                                  VectorizedArray<Number>,
-                                                  VectorizedArray<Number>>::
-              do_backward(1,
-                          prolongation_matrix_1d,
-                          false,
-                          evaluation_data.begin() + c * n_scalar_cell_dofs,
-                          evaluation_data.begin() +
-                            c * Utilities::fixed_power<dim>(degree_size),
-                          fe_degree + 1,
-                          2 * fe_degree + 2);
+            internal::FEEvaluationImplBasisChange<
+              internal::evaluate_general,
+              internal::EvaluatorQuantity::value,
+              dim,
+              degree + 1,
+              2 * degree + 2,
+              VectorizedArray<Number>,
+              VectorizedArray<Number>>::do_backward(1,
+                                                    prolongation_matrix_1d,
+                                                    false,
+                                                    evaluation_data.begin() +
+                                                      c * n_scalar_cell_dofs,
+                                                    evaluation_data.begin() +
+                                                      c *
+                                                        Utilities::fixed_power<
+                                                          dim>(degree_size),
+                                                    fe_degree + 1,
+                                                    2 * fe_degree + 2);
         }
 
       // write into dst vector