Do not use transformation to collocation for large n_q_points.

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 7cfa435b6ad3dd7aaa58d0aeacfe1cc58248101d..3f47fbe06eb12f496af7acb45c72798c5bd3d63f 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -1337,6 +1337,16 @@ namespace internal
             typename Number>
   struct FEFaceEvaluationImpl
   {
+    // We enable a transformation to collocation for derivatives if it gives
+    // correct results (first two conditions), if it is the most efficient
+    // choice in terms of operation counts (third condition) and if we were
+    // able to initialize the fields in shape_info.templates.h from the
+    // polynomials (fourth condition).
+    static constexpr bool use_collocation =
+      symmetric_evaluate &&
+      n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
+      n_q_points_1d < 200;
+
     static void
     evaluate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
                      Number *                                      values_dofs,
@@ -1431,7 +1441,7 @@ namespace internal
             switch (dim)
               {
                 case 3:
-                  if (symmetric_evaluate && n_q_points_1d > fe_degree)
+                  if (use_collocation)
                     {
                       eval1.template values<0, true, false>(values_dofs,
                                                             values_quad);
@@ -1590,7 +1600,7 @@ namespace internal
                                                            2 * n_q_points);
                   eval1.template values<0, false, false>(
                     gradients_quad + 2 * n_q_points, values_dofs + size_deg);
-                  if (symmetric_evaluate && n_q_points_1d > fe_degree)
+                  if (use_collocation)
                     {
                       internal::EvaluatorTensorProduct<
                         internal::evaluate_evenodd,

diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h
index 2fd382aa06478a2789915cd9400d08d35263d4ca..30115e481be149d72a3710d3d2beef81aef89f23 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_selector.h
+++ b/include/deal.II/matrix_free/evaluation_selector.h
@@ -251,6 +251,17 @@ namespace internal
                    n_q_points_1d,
                    typename std::enable_if<(n_q_points_1d < degree + 3)>::type>
     {
+      /**
+       * We enable a transformation to collocation for derivatives if it gives
+       * correct results (first condition), if it is the most efficient choice
+       * in terms of operation counts (second condition) and if we were able
+       * to initialize the fields in shape_info.templates.h from the
+       * polynomials (third condition).
+       */
+      static constexpr bool      use_collocation =
+        n_q_points_1d > degree &&n_q_points_1d <= 3 * degree / 2 + 1 &&
+        n_q_points_1d < 200;
+
       static inline void
       evaluate(
         const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
@@ -280,7 +291,7 @@ namespace internal
                            evaluate_values,
                            evaluate_gradients,
                            evaluate_hessians);
-            else if (degree < n_q_points_1d)
+            else if (use_collocation)
               internal::FEEvaluationImplTransformToCollocation<
                 dim,
                 degree,
@@ -352,7 +363,7 @@ namespace internal
                             integrate_values,
                             integrate_gradients,
                             sum_into_values_array);
-            else if (degree < n_q_points_1d)
+            else if (use_collocation)
               internal::FEEvaluationImplTransformToCollocation<
                 dim,
                 degree,
@@ -481,6 +492,17 @@ template <int dim,
           typename Number>
 struct SelectEvaluator
 {
+  /**
+   * We enable a transformation to collocation for derivatives if it gives
+   * correct results (first condition), if it is the most efficient choice in
+   * terms of operation counts (second condition) and if we were able to
+   * initialize the fields in shape_info.templates.h from the polynomials
+   * (third condition).
+   */
+  static constexpr bool         use_collocation =
+    n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
+    n_q_points_1d < 200;
+
   /**
    * Chooses an appropriate evaluation strategy for the evaluate function, i.e.
    * this calls internal::FEEvaluationImpl::evaluate(),
@@ -606,9 +628,10 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate(
                    evaluate_gradients,
                    evaluate_hessians);
     }
-  else if (fe_degree < n_q_points_1d &&
-           shape_info.element_type <=
-             internal::MatrixFreeFunctions::tensor_symmetric)
+  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
+  // shape_info.h for more details
+  else if (use_collocation && shape_info.element_type <=
+                                internal::MatrixFreeFunctions::tensor_symmetric)
     {
       internal::FEEvaluationImplTransformToCollocation<
         dim,
@@ -625,7 +648,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate(
                           evaluate_gradients,
                           evaluate_hessians);
     }
-  else if (shape_info.element_type ==
+  else if (shape_info.element_type <=
            internal::MatrixFreeFunctions::tensor_symmetric)
     {
       internal::FEEvaluationImpl<
@@ -739,9 +762,10 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                     integrate_gradients,
                     sum_into_values_array);
     }
-  else if (fe_degree < n_q_points_1d &&
-           shape_info.element_type <=
-             internal::MatrixFreeFunctions::tensor_symmetric)
+  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
+  // shape_info.h for more details
+  else if (use_collocation && shape_info.element_type <=
+                                internal::MatrixFreeFunctions::tensor_symmetric)
     {
       internal::FEEvaluationImplTransformToCollocation<
         dim,
@@ -757,7 +781,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                            integrate_gradients,
                            sum_into_values_array);
     }
-  else if (shape_info.element_type ==
+  else if (shape_info.element_type <=
            internal::MatrixFreeFunctions::tensor_symmetric)
     {
       internal::FEEvaluationImpl<