From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Fri, 18 Aug 2017 10:35:18 +0000 (+0200)
Subject: Consider n_q_points_1d in [degree+1, degree+3) instead
X-Git-Tag: v9.0.0-rc1~1180^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c374c5cd81b17cbab0edb139ff81cf9f687be698;p=dealii.git

Consider n_q_points_1d in [degree+1, degree+3) instead
---

diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h
index 680be1b330..0c74cfa717 100644
--- a/include/deal.II/matrix_free/evaluation_selector.h
+++ b/include/deal.II/matrix_free/evaluation_selector.h
@@ -33,10 +33,10 @@ namespace
 //    parameter, increase fe_degree  by one and try again.
 //    If fe_degree==10 use the class Default which serves as a fallback.
 // 3. After fixing the fe_degree, DEPTH is increased (DEPTH=1) and we start with
-//    n_q_points=fe_degree-1.
+//    n_q_points=fe_degree+1.
 // 4. If the current assumption on n_q_points_1d doesn't match the runtime
 //    parameter, increase n_q_points_1d by one and try again.
-//    If n_q_points_1d==degree+2 use the class Default which serves as a fallback.
+//    If n_q_points_1d==degree+3 use the class Default which serves as a fallback.
 
   /**
    * This class serves as a fallback in case we don't have the appropriate template
@@ -99,7 +99,7 @@ namespace
    */
   template<int degree, int n_q_points_1d, int dim, int n_components, typename Number>
   struct Factory<dim, n_components, Number, 1, degree, n_q_points_1d,
-    typename std::enable_if<n_q_points_1d==degree+2>::type> : Default<dim, n_components, Number> {};
+    typename std::enable_if<n_q_points_1d==degree+3>::type> : Default<dim, n_components, Number> {};
 
   /**
    * This class chooses the correct template degree.
@@ -119,9 +119,9 @@ namespace
       const bool               evaluate_hessians)
     {
       const unsigned int runtime_degree = shape_info.fe_degree;
-      constexpr unsigned int nonnegative_start_n_q_points = (degree>0)?degree-1:0;
+      constexpr unsigned int start_n_q_points = degree+1;
       if (runtime_degree == degree)
-        Factory<dim, n_components, Number, 1, degree, nonnegative_start_n_q_points>::evaluate
+        Factory<dim, n_components, Number, 1, degree, start_n_q_points>::evaluate
         (shape_info, values_dofs_actual, values_quad, gradients_quad, hessians_quad,
          scratch_data, evaluate_values, evaluate_gradients, evaluate_hessians);
       else
@@ -140,9 +140,9 @@ namespace
       const bool               integrate_gradients)
     {
       const int runtime_degree = shape_info.fe_degree;
-      constexpr unsigned int nonnegative_start_n_q_points = (degree>0)?degree-1:0;
+      constexpr unsigned int start_n_q_points = degree+1;
       if (runtime_degree == degree)
-        Factory<dim, n_components, Number, 1, degree, nonnegative_start_n_q_points>::integrate
+        Factory<dim, n_components, Number, 1, degree, start_n_q_points>::integrate
         (shape_info, values_dofs_actual, values_quad, gradients_quad,
          scratch_data, integrate_values, integrate_gradients);
       else
@@ -156,8 +156,7 @@ namespace
    * This class chooses the correct template n_q_points_1d after degree was chosen.
    */
   template<int degree, int n_q_points_1d, int dim, int n_components, typename Number>
-  struct Factory<dim, n_components, Number, 1, degree, n_q_points_1d,
-    typename std::enable_if<n_q_points_1d<degree+2>::type>
+  struct Factory<dim, n_components, Number, 1, degree, n_q_points_1d, typename std::enable_if<(n_q_points_1d<degree+3)>::type>
   {
     static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
                                  VectorizedArray<Number> *values_dofs_actual[],
@@ -168,69 +167,69 @@ namespace
                                  const bool               evaluate_values,
                                  const bool               evaluate_gradients,
                                  const bool               evaluate_hessians)
-    {
-      const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
-      if (runtime_n_q_points_1d == n_q_points_1d)
-        {
-          if (n_q_points_1d == degree+1)
-            {
-              if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
-                internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
-                ::evaluate(shape_info, values_dofs_actual, values_quad,
-                           gradients_quad, hessians_quad, scratch_data,
-                           evaluate_values, evaluate_gradients, evaluate_hessians);
-              else
-                internal::FEEvaluationImplTransformToCollocation<dim, degree, n_components, Number>
-                ::evaluate(shape_info, values_dofs_actual, values_quad,
-                           gradients_quad, hessians_quad, scratch_data,
-                           evaluate_values, evaluate_gradients, evaluate_hessians);
-            }
-          else
-            internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric, dim, degree, n_q_points_1d, n_components, Number>
-            ::evaluate(shape_info, values_dofs_actual, values_quad,
-                       gradients_quad, hessians_quad, scratch_data,
-                       evaluate_values, evaluate_gradients, evaluate_hessians);
-        }
-      else
-        Factory<dim, n_components, Number, 1, degree, n_q_points_1d+1>::evaluate (shape_info, values_dofs_actual, values_quad,
-            gradients_quad, hessians_quad, scratch_data,
-            evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
+  {
+    const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
+    if (runtime_n_q_points_1d == n_q_points_1d)
+      {
+        if (n_q_points_1d == degree+1)
+          {
+            if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+              internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
+              ::evaluate(shape_info, values_dofs_actual, values_quad,
+                         gradients_quad, hessians_quad, scratch_data,
+                         evaluate_values, evaluate_gradients, evaluate_hessians);
+            else
+              internal::FEEvaluationImplTransformToCollocation<dim, degree, n_components, Number>
+              ::evaluate(shape_info, values_dofs_actual, values_quad,
+                         gradients_quad, hessians_quad, scratch_data,
+                         evaluate_values, evaluate_gradients, evaluate_hessians);
+          }
+        else
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric, dim, degree, n_q_points_1d, n_components, Number>
+          ::evaluate(shape_info, values_dofs_actual, values_quad,
+                     gradients_quad, hessians_quad, scratch_data,
+                     evaluate_values, evaluate_gradients, evaluate_hessians);
+      }
+    else
+      Factory<dim, n_components, Number, 1, degree, n_q_points_1d+1>::evaluate (shape_info, values_dofs_actual, values_quad,
+          gradients_quad, hessians_quad, scratch_data,
+          evaluate_values, evaluate_gradients, evaluate_hessians);
+  }
 
-    static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
-                                  VectorizedArray<Number> *values_dofs_actual[],
-                                  VectorizedArray<Number> *values_quad[],
-                                  VectorizedArray<Number> *gradients_quad[][dim],
-                                  VectorizedArray<Number> *scratch_data,
-                                  const bool               integrate_values,
-                                  const bool               integrate_gradients)
-    {
-      const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
-      if (runtime_n_q_points_1d == n_q_points_1d)
-        {
-          if (n_q_points_1d == degree+1)
-            {
-              if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
-                internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
-                ::integrate(shape_info, values_dofs_actual, values_quad,
-                            gradients_quad, scratch_data,
-                            integrate_values, integrate_gradients);
-              else
-                internal::FEEvaluationImplTransformToCollocation<dim, degree, n_components, Number>
-                ::integrate(shape_info, values_dofs_actual, values_quad,
-                            gradients_quad, scratch_data,
-                            integrate_values, integrate_gradients);
-            }
-          else
-            internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric, dim, degree, n_q_points_1d, n_components, Number>
-            ::integrate(shape_info, values_dofs_actual, values_quad, gradients_quad,
-                        scratch_data, integrate_values, integrate_gradients);
-        }
-      else
-        Factory<dim, n_components, Number, 1, degree, n_q_points_1d+1>
-        ::integrate (shape_info, values_dofs_actual, values_quad, gradients_quad,
-                     scratch_data, integrate_values, integrate_gradients);
-    }
+  static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+                                VectorizedArray<Number> *values_dofs_actual[],
+                                VectorizedArray<Number> *values_quad[],
+                                VectorizedArray<Number> *gradients_quad[][dim],
+                                VectorizedArray<Number> *scratch_data,
+                                const bool               integrate_values,
+                                const bool               integrate_gradients)
+  {
+    const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
+    if (runtime_n_q_points_1d == n_q_points_1d)
+      {
+        if (n_q_points_1d == degree+1)
+          {
+            if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+              internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
+              ::integrate(shape_info, values_dofs_actual, values_quad,
+                          gradients_quad, scratch_data,
+                          integrate_values, integrate_gradients);
+            else
+              internal::FEEvaluationImplTransformToCollocation<dim, degree, n_components, Number>
+              ::integrate(shape_info, values_dofs_actual, values_quad,
+                          gradients_quad, scratch_data,
+                          integrate_values, integrate_gradients);
+          }
+        else
+          internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric, dim, degree, n_q_points_1d, n_components, Number>
+          ::integrate(shape_info, values_dofs_actual, values_quad, gradients_quad,
+                      scratch_data, integrate_values, integrate_gradients);
+      }
+    else
+      Factory<dim, n_components, Number, 1, degree, n_q_points_1d+1>
+      ::integrate (shape_info, values_dofs_actual, values_quad, gradients_quad,
+                   scratch_data, integrate_values, integrate_gradients);
+  }
   };