]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Choose appropriate tensorized evaluation template instantiation at run-time
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 4 Aug 2017 22:11:42 +0000 (00:11 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 18 Aug 2017 09:41:46 +0000 (11:41 +0200)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_selector.h [new file with mode: 0644]
include/deal.II/matrix_free/fe_evaluation.h

index 2da6db43163316ca25929c691d46d285bc7c6c29..105f47ed1ffdacc212c17237f1c1e12785546ea6 100644 (file)
@@ -376,7 +376,8 @@ namespace internal
                                    : (Eval::dofs_per_cell > Eval::n_q_points ?
                                       Eval::dofs_per_cell : Eval::n_q_points);
     const unsigned int max_stack_size = 100;
-    VectorizedArray<Number>  temp_data[(temp_size > 0 && temp_size < max_stack_size) ? 2*temp_size : 1];
+    AlignedVector<VectorizedArray<Number>>  temp_data;
+    temp_data.resize_fast((temp_size > 0 && temp_size < max_stack_size) ? 2*temp_size : 1);
     VectorizedArray<Number> *temp1;
     VectorizedArray<Number> *temp2;
     if (temp_size == 0)
diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h
new file mode 100644 (file)
index 0000000..680be1b
--- /dev/null
@@ -0,0 +1,574 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii__matrix_free_evaluation_selector_h
+#define dealii__matrix_free_evaluation_selector_h
+
+#include <deal.II/matrix_free/evaluation_kernels.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#ifndef DOXYGEN
+namespace
+{
+// The following classes serve the purpose of choosing the correct template
+// specialization of the FEEvaluationImpl* classes in case fe_degree
+// and n_q_points_1d are only given as runtime parameters.
+// The logic is the following:
+// 1. Start with fe_degree=0, n_q_points_1d=0 and DEPTH=0.
+// 2. If the current assumption on fe_degree doesn't match the runtime
+//    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.
+// 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.
+
+  /**
+   * This class serves as a fallback in case we don't have the appropriate template
+   * specialization for the run time and template parameters given.
+   */
+  template <int dim, int n_components, typename Number>
+  struct Default
+  {
+    static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+                                 VectorizedArray<Number> *values_dofs_actual[],
+                                 VectorizedArray<Number> *values_quad[],
+                                 VectorizedArray<Number> *gradients_quad[][dim],
+                                 VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                                 VectorizedArray<Number> *scratch_data,
+                                 const bool               evaluate_values,
+                                 const bool               evaluate_gradients,
+                                 const bool               evaluate_hessians)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+               dim, -1, 0, n_components, Number>
+               ::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)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+               dim, -1, 0, n_components, Number>
+               ::integrate(shape_info, values_dofs_actual, values_quad,
+                           gradients_quad, scratch_data,
+                           integrate_values, integrate_gradients);
+    }
+  };
+
+
+  /**
+   * This class implements the actual choice of the template specialization.
+   */
+  template<int dim, int n_components, typename Number,
+           int DEPTH=0, int degree=0, int n_q_points_1d=0, class Enable = void>
+  struct Factory : Default<dim, n_components, Number> {};
+
+  /**
+   * This specialization sets the maximal fe_degree for
+   * which we want to determine the correct template parameters based at runtime.
+   */
+  template<int n_q_points_1d, int dim, int n_components, typename Number>
+  struct Factory<dim, n_components, Number, 0, 10, n_q_points_1d> : Default<dim, n_components, Number> {};
+
+  /**
+   * This specialization sets the maximal number of n_q_points_1d for
+   * which we want to determine the correct template parameters based at runtime.
+   */
+  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> {};
+
+  /**
+   * This class chooses the correct template degree.
+   */
+  template<int degree, int n_q_points_1d, int dim, int n_components, typename Number>
+  struct Factory<dim, n_components, Number, 0, degree, n_q_points_1d>
+  {
+    static inline void evaluate (
+      const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+      VectorizedArray<Number> *values_dofs_actual[],
+      VectorizedArray<Number> *values_quad[],
+      VectorizedArray<Number> *gradients_quad[][dim],
+      VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+      VectorizedArray<Number> *scratch_data,
+      const bool               evaluate_values,
+      const bool               evaluate_gradients,
+      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;
+      if (runtime_degree == degree)
+        Factory<dim, n_components, Number, 1, degree, nonnegative_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
+        Factory<dim, n_components, Number, 0, degree+1, n_q_points_1d>::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_degree = shape_info.fe_degree;
+      constexpr unsigned int nonnegative_start_n_q_points = (degree>0)?degree-1:0;
+      if (runtime_degree == degree)
+        Factory<dim, n_components, Number, 1, degree, nonnegative_start_n_q_points>::integrate
+        (shape_info, values_dofs_actual, values_quad, gradients_quad,
+         scratch_data, integrate_values, integrate_gradients);
+      else
+        Factory<dim, n_components, Number, 0, degree+1, n_q_points_1d>::integrate
+        (shape_info, values_dofs_actual, values_quad, gradients_quad,
+         scratch_data, integrate_values, integrate_gradients);
+    }
+  };
+
+  /**
+   * 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>
+  {
+    static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+                                 VectorizedArray<Number> *values_dofs_actual[],
+                                 VectorizedArray<Number> *values_quad[],
+                                 VectorizedArray<Number> *gradients_quad[][dim],
+                                 VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                                 VectorizedArray<Number> *scratch_data,
+                                 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);
+    }
+
+    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);
+    }
+  };
+
+
+
+  /**
+   * This is the entry point for choosing the correct runtime parameters
+   * for the 'evaluate' function.
+   */
+  template<int dim, int n_components, typename Number>
+  void symmetric_selector_evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+                                    VectorizedArray<Number> *values_dofs_actual[],
+                                    VectorizedArray<Number> *values_quad[],
+                                    VectorizedArray<Number> *gradients_quad[][dim],
+                                    VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                                    VectorizedArray<Number> *scratch_data,
+                                    const bool               evaluate_values,
+                                    const bool               evaluate_gradients,
+                                    const bool               evaluate_hessians)
+  {
+    Assert(shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric||
+           shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation,
+           ExcInternalError());
+    Factory<dim, n_components, Number>::evaluate
+    (shape_info, values_dofs_actual, values_quad, gradients_quad, hessians_quad,
+     scratch_data, evaluate_values, evaluate_gradients, evaluate_hessians);
+  }
+
+
+
+  /**
+   * This is the entry point for choosing the correct runtime parameters
+   * for the 'integrate' function.
+   */
+  template<int dim, int n_components, typename Number>
+  void symmetric_selector_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)
+  {
+    Assert(shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric||
+           shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation,
+           ExcInternalError());
+    Factory<dim, n_components, Number>::integrate
+    (shape_info, values_dofs_actual, values_quad, gradients_quad,
+     scratch_data, integrate_values, integrate_gradients);
+  }
+}
+#endif
+
+
+/**
+ * This class chooses an appropriate evaluation strategy based on the
+ * template parameters and the shape_info variable which contains runtime
+ * parameters. In case the template parameters fe_degree and n_q_points_1d
+ * contain valid information (i.e. fe_degree>-1 and n_q_points_1d>0), we simply
+ * pass these values to the respective template specializations.
+ * Otherwise, we perform a runtime matching of the runtime parameters to find
+ * the correct specialization. This matching currently supports
+ * $0\leq fe\_degree$ and $degree-1\leq n_q_points_1d\leq degree+1$.
+ */
+template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+struct SelectEvaluator
+{
+  static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+                       VectorizedArray<Number> *values_dofs_actual[],
+                       VectorizedArray<Number> *values_quad[],
+                       VectorizedArray<Number> *gradients_quad[][dim],
+                       VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                       VectorizedArray<Number> *scratch_data,
+                       const bool               evaluate_values,
+                       const bool               evaluate_gradients,
+                       const bool               evaluate_hessians);
+
+  static 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);
+};
+
+/**
+ * Specialization for invalid template parameters, i.e. fe_degree==-1.
+ */
+template <int dim, int n_q_points_1d, int n_components, typename Number>
+struct SelectEvaluator<dim, -1, n_q_points_1d, n_components, Number>
+{
+  static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+                       VectorizedArray<Number> *values_dofs_actual[],
+                       VectorizedArray<Number> *values_quad[],
+                       VectorizedArray<Number> *gradients_quad[][dim],
+                       VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+                       VectorizedArray<Number> *scratch_data,
+                       const bool               evaluate_values,
+                       const bool               evaluate_gradients,
+                       const bool               evaluate_hessians);
+
+  static 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);
+};
+
+//----------------------Implementation for SelectEvaluator---------------------
+#ifndef DOXYGEN
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+inline
+void
+SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate
+(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+ VectorizedArray<Number> *values_dofs_actual[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+ VectorizedArray<Number> *scratch_data,
+ const bool               evaluate_values,
+ const bool               evaluate_gradients,
+ const bool               evaluate_hessians)
+{
+  Assert(fe_degree>=0  && n_q_points_1d>0, ExcInternalError());
+
+  if (fe_degree+1 == n_q_points_1d &&
+      shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+    {
+      internal::FEEvaluationImplCollocation<dim, fe_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 if (fe_degree+1 == n_q_points_1d &&
+           shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
+    {
+      internal::FEEvaluationImplTransformToCollocation<dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
+               dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+               dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+               dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+               dim, fe_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
+    AssertThrow(false, ExcNotImplemented());
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+inline
+void
+SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::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)
+{
+  Assert(fe_degree>=0  && n_q_points_1d>0, ExcInternalError());
+
+  if (fe_degree+1 == n_q_points_1d &&
+      shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+    {
+      internal::FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
+      ::integrate(shape_info, values_dofs_actual, values_quad,
+                  gradients_quad, scratch_data,
+                  integrate_values, integrate_gradients);
+    }
+  else if (fe_degree+1 == n_q_points_1d &&
+           shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
+    {
+      internal::FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
+      ::integrate(shape_info, values_dofs_actual, values_quad,
+                  gradients_quad, scratch_data,
+                  integrate_values, integrate_gradients);
+    }
+  else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
+               dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+               dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+               dim, fe_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 if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+               dim, fe_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
+    AssertThrow(false, ExcNotImplemented());
+}
+
+
+
+template <int dim, int dummy, int n_components, typename Number>
+inline
+void
+SelectEvaluator<dim, -1, dummy, n_components, Number>::evaluate
+(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
+ VectorizedArray<Number> *values_dofs_actual[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+ VectorizedArray<Number> *scratch_data,
+ const bool               evaluate_values,
+ const bool               evaluate_gradients,
+ const bool               evaluate_hessians)
+{
+  if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+               dim, -1, 0, n_components, Number>
+               ::evaluate(shape_info, values_dofs_actual, values_quad,
+                          gradients_quad, hessians_quad, scratch_data,
+                          evaluate_values, evaluate_gradients, evaluate_hessians);
+    }
+  else if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+               dim, -1, 0, n_components, Number>
+               ::evaluate(shape_info, values_dofs_actual, values_quad,
+                          gradients_quad, hessians_quad, scratch_data,
+                          evaluate_values, evaluate_gradients, evaluate_hessians);
+    }
+  else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general)
+    internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+             dim, -1, 0, n_components, Number>
+             ::evaluate(shape_info, values_dofs_actual, values_quad,
+                        gradients_quad, hessians_quad, scratch_data,
+                        evaluate_values, evaluate_gradients, evaluate_hessians);
+  else
+    symmetric_selector_evaluate<dim, n_components, Number>
+    (shape_info, values_dofs_actual, values_quad,
+     gradients_quad, hessians_quad, scratch_data,
+     evaluate_values, evaluate_gradients, evaluate_hessians);
+}
+
+
+
+template <int dim, int dummy, int n_components, typename Number>
+inline
+void
+SelectEvaluator<dim, -1, dummy, n_components, Number>::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)
+{
+  if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
+               dim, -1, 0, n_components, Number>
+               ::integrate(shape_info, values_dofs_actual, values_quad,
+                           gradients_quad, scratch_data,
+                           integrate_values, integrate_gradients);
+    }
+  else if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor)
+    {
+      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
+               dim, -1, 0, n_components, Number>
+               ::integrate(shape_info, values_dofs_actual, values_quad,
+                           gradients_quad, scratch_data,
+                           integrate_values, integrate_gradients);
+    }
+  else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general)
+    internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
+             dim, -1, 0, n_components, Number>
+             ::integrate(shape_info, values_dofs_actual, values_quad,
+                         gradients_quad, scratch_data,
+                         integrate_values, integrate_gradients);
+  else
+    symmetric_selector_integrate<dim, n_components, Number>
+    (shape_info, values_dofs_actual, values_quad,
+     gradients_quad, scratch_data,
+     integrate_values, integrate_gradients);
+}
+#endif //DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index bfd8d3b42a97451d890c2abc56d05b5c047b627b..eca7223d3a73bf875ba4545980cd828aa0987e87 100644 (file)
@@ -30,6 +30,7 @@
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/evaluation_kernels.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/evaluation_selector.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -5338,66 +5339,10 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   Assert(this->matrix_info != nullptr ||
          this->mapped_geometry->is_initialized(), ExcNotInitialized());
 
-  // Run time selection of algorithm matching the element type. Since some
-  // template combinations do not work for fe_degree==-1, create safe template
-  // values for the calls below.
-  const unsigned int fe_degree_nonnegative = fe_degree != -1 ? fe_degree : 0;
-  const unsigned int n_q_points_1d_positive = n_q_points_1d > 0 ? n_q_points_1d : 1;
-  const unsigned int n_q_points_1d_adjusted = fe_degree == -1 ? 0 : n_q_points_1d;
-  if (fe_degree >= 0 &&
-      fe_degree+1 == n_q_points_1d &&
-      this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
-    {
-      internal::FEEvaluationImplCollocation<dim, fe_degree_nonnegative, n_components_, Number>
-      ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
-                 this->gradients_quad, this->hessians_quad, this->scratch_data,
-                 evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
-  else if (fe_degree >= 0 &&
-           fe_degree+1 == n_q_points_1d &&
-           this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
-    {
-      internal::FEEvaluationImplTransformToCollocation<dim, fe_degree_nonnegative, n_components_, Number>
-      ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
-                 this->gradients_quad, this->hessians_quad, this->scratch_data,
-                 evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
-  else if (fe_degree >= 0 &&
-           this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
-               dim, fe_degree_nonnegative, n_q_points_1d_positive, n_components_, Number>
-               ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
-                          this->gradients_quad, this->hessians_quad, this->scratch_data,
-                          evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
-  else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-               dim, fe_degree, n_q_points_1d_adjusted, n_components_, Number>
-               ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
-                          this->gradients_quad, this->hessians_quad, this->scratch_data,
-                          evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
-  else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
-               dim, fe_degree, n_q_points_1d_adjusted, n_components_, Number>
-               ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
-                          this->gradients_quad, this->hessians_quad, this->scratch_data,
-                          evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
-  else if (fe_degree == -1 ||
-           this->data->element_type == internal::MatrixFreeFunctions::tensor_general)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
-               dim, fe_degree, n_q_points_1d_adjusted, n_components_, Number>
-               ::evaluate(*this->data, &this->values_dofs[0], this->values_quad,
-                          this->gradients_quad, this->hessians_quad, this->scratch_data,
-                          evaluate_values, evaluate_gradients, evaluate_hessians);
-    }
-  else
-    AssertThrow(false, ExcNotImplemented());
+  SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>
+  ::evaluate (*this->data, &this->values_dofs[0], this->values_quad,
+              this->gradients_quad, this->hessians_quad, this->scratch_data,
+              evaluate_values, evaluate_gradients, evaluate_hessians);
 
 #ifdef DEBUG
   if (evaluate_values == true)
@@ -5428,66 +5373,10 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   Assert(this->matrix_info != nullptr ||
          this->mapped_geometry->is_initialized(), ExcNotInitialized());
 
-  // Run time selection of algorithm matching the element type. Since some
-  // template combinations do not work for fe_degree==-1, create safe template
-  // values for the calls below.
-  const unsigned int fe_degree_nonnegative = fe_degree != -1 ? fe_degree : 0;
-  const unsigned int n_q_points_1d_positive = n_q_points_1d > 0 ? n_q_points_1d : 1;
-  const unsigned int n_q_points_1d_adjusted = fe_degree == -1 ? 0 : n_q_points_1d;
-  if (fe_degree >= 0 &&
-      fe_degree+1 == n_q_points_1d &&
-      this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
-    {
-      internal::FEEvaluationImplCollocation<dim, fe_degree_nonnegative, n_components_, Number>
-      ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
-                  this->gradients_quad, this->scratch_data,
-                  integrate_values, integrate_gradients);
-    }
-  else if (fe_degree >= 0 &&
-           fe_degree+1 == n_q_points_1d &&
-           this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
-    {
-      internal::FEEvaluationImplTransformToCollocation<dim, fe_degree_nonnegative, n_components_, Number>
-      ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
-                  this->gradients_quad, this->scratch_data,
-                  integrate_values, integrate_gradients);
-    }
-  else if (fe_degree >= 0 &&
-           this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric,
-               dim, fe_degree_nonnegative, n_q_points_1d_positive, n_components_, Number>
-               ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
-                           this->gradients_quad, this->scratch_data,
-                           integrate_values, integrate_gradients);
-    }
-  else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0,
-               dim, fe_degree, n_q_points_1d_adjusted, n_components_, Number>
-               ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
-                           this->gradients_quad, this->scratch_data,
-                           integrate_values, integrate_gradients);
-    }
-  else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::truncated_tensor,
-               dim, fe_degree, n_q_points_1d_adjusted, n_components_, Number>
-               ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
-                           this->gradients_quad, this->scratch_data,
-                           integrate_values, integrate_gradients);
-    }
-  else if (fe_degree == -1 ||
-           this->data->element_type == internal::MatrixFreeFunctions::tensor_general)
-    {
-      internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
-               dim, fe_degree, n_q_points_1d_adjusted, n_components_, Number>
-               ::integrate(*this->data, &this->values_dofs[0], this->values_quad,
-                           this->gradients_quad, this->scratch_data,
-                           integrate_values, integrate_gradients);
-    }
-  else
-    AssertThrow(false, ExcNotImplemented());
+  SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>
+  ::integrate (*this->data, &this->values_dofs[0], this->values_quad,
+               this->gradients_quad, this->scratch_data,
+               integrate_values, integrate_gradients);
 
 #ifdef DEBUG
   this->dof_values_initialized = true;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.