From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 17 Sep 2020 14:20:16 +0000 (+0200)
Subject: Use new factory class within library
X-Git-Tag: v9.3.0-rc1~1087^2~7
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=42e7433ca9a2f8c9b24e718f68f3afa9a5bc56e0;p=dealii.git

Use new factory class within library
---

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 99b0349a47..c7561253db 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -28,8 +28,7 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q_generic.h>
 
-#include <deal.II/matrix_free/evaluation_kernels.h>
-#include <deal.II/matrix_free/evaluation_selector.h>
+#include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/mapping_info.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -1302,7 +1301,7 @@ namespace internal
                                                 start_indices,
                                                 cell_points.data());
 
-                  SelectEvaluator<dim, -1, 0, VectorizedDouble>::evaluate(
+                  FEEvaluationFactory<dim, double, VectorizedDouble>::evaluate(
                     dim,
                     EvaluationFlags::values | EvaluationFlags::gradients |
                       (update_flags_cells & update_jacobian_grads ?
@@ -2126,21 +2125,21 @@ namespace internal
 
               // now let the matrix-free evaluators provide us with the
               // data on faces
-              FEFaceEvaluationImplEvaluateSelector<dim, VectorizedDouble>::
-                template run<-1, 0>(dim,
-                                    shape_info,
-                                    cell_points.data(),
-                                    face_quads.data(),
-                                    face_grads.data(),
-                                    scratch_data.data(),
-                                    true,
-                                    true,
-                                    face_no,
-                                    GeometryInfo<dim>::max_children_per_cell,
-                                    faces[face].face_orientation > 8 ?
-                                      faces[face].face_orientation - 8 :
-                                      0,
-                                    my_data.descriptor[0].face_orientations);
+              FEFaceEvaluationFactory<dim, double, VectorizedDouble>::evaluate(
+                dim,
+                shape_info,
+                cell_points.data(),
+                face_quads.data(),
+                face_grads.data(),
+                scratch_data.data(),
+                true,
+                true,
+                face_no,
+                GeometryInfo<dim>::max_children_per_cell,
+                faces[face].face_orientation > 8 ?
+                  faces[face].face_orientation - 8 :
+                  0,
+                my_data.descriptor[0].face_orientations);
 
 
               if (update_flags_faces & update_quadrature_points)
@@ -2243,22 +2242,21 @@ namespace internal
                                                 start_indices,
                                                 cell_points.data());
 
-                  FEFaceEvaluationImplEvaluateSelector<dim, VectorizedDouble>::
-                    template run<-1, 0>(
-                      dim,
-                      shape_info,
-                      cell_points.data(),
-                      face_quads.data(),
-                      face_grads.data(),
-                      scratch_data.data(),
-                      false,
-                      true,
-                      faces[face].exterior_face_no,
-                      faces[face].subface_index,
-                      faces[face].face_orientation < 8 ?
-                        faces[face].face_orientation :
-                        0,
-                      my_data.descriptor[0].face_orientations);
+                  FEFaceEvaluationFactory<dim, double, VectorizedDouble>::
+                    evaluate(dim,
+                             shape_info,
+                             cell_points.data(),
+                             face_quads.data(),
+                             face_grads.data(),
+                             scratch_data.data(),
+                             false,
+                             true,
+                             faces[face].exterior_face_no,
+                             faces[face].subface_index,
+                             faces[face].face_orientation < 8 ?
+                               faces[face].face_orientation :
+                               0,
+                             my_data.descriptor[0].face_orientations);
 
                   for (unsigned int q = 0; q < n_points_compute; ++q)
                     {
diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h
index 6db682e4e1..45d366dd92 100644
--- a/include/deal.II/numerics/vector_tools_project.templates.h
+++ b/include/deal.II/numerics/vector_tools_project.templates.h
@@ -149,13 +149,9 @@ namespace VectorTools
 
     /*
      * MatrixFree implementation of project() for an arbitrary number of
-     * components and arbitrary degree of the FiniteElement.
+     * components of the FiniteElement.
      */
-    template <int components,
-              int fe_degree,
-              int dim,
-              typename Number,
-              int spacedim>
+    template <int components, int dim, typename Number, int spacedim>
     void
     project_matrix_free(
       const Mapping<dim, spacedim> &   mapping,
@@ -180,9 +176,6 @@ namespace VectorTools
       Assert(dof.get_fe(0).n_components() == function.n_components,
              ExcDimensionMismatch(dof.get_fe(0).n_components(),
                                   function.n_components));
-      Assert(fe_degree == -1 ||
-               dof.get_fe().degree == static_cast<unsigned int>(fe_degree),
-             ExcDimensionMismatch(fe_degree, dof.get_fe().degree));
       Assert(dof.get_fe(0).n_components() == components,
              ExcDimensionMismatch(components, dof.get_fe(0).n_components()));
 
@@ -201,8 +194,8 @@ namespace VectorTools
                           additional_data);
       using MatrixType = MatrixFreeOperators::MassOperator<
         dim,
-        fe_degree,
-        fe_degree + 2,
+        -1,
+        0,
         components,
         LinearAlgebra::distributed::Vector<Number>>;
       MatrixType mass_matrix;
@@ -222,8 +215,7 @@ namespace VectorTools
 
         // account for inhomogeneous constraints
         inhomogeneities.update_ghost_values();
-        FEEvaluation<dim, fe_degree, fe_degree + 2, components, Number> phi(
-          *matrix_free);
+        FEEvaluation<dim, -1, 0, components, Number> phi(*matrix_free);
         for (unsigned int cell = 0; cell < matrix_free->n_macro_cells(); ++cell)
           {
             phi.reinit(cell);
@@ -255,81 +247,6 @@ namespace VectorTools
 
 
 
-    /**
-     * Helper interface. After figuring out the number of components in
-     * project_matrix_free_component, we determine the degree of the
-     * FiniteElement and call project_matrix_free with the appropriate
-     * template arguments.
-     */
-    template <int components, int dim, typename Number, int spacedim>
-    void
-    project_matrix_free_degree(
-      const Mapping<dim, spacedim> &   mapping,
-      const DoFHandler<dim, spacedim> &dof,
-      const AffineConstraints<Number> &constraints,
-      const Quadrature<dim> &          quadrature,
-      const Function<
-        spacedim,
-        typename LinearAlgebra::distributed::Vector<Number>::value_type>
-        &                                         function,
-      LinearAlgebra::distributed::Vector<Number> &work_result,
-      const bool                                  enforce_zero_boundary,
-      const Quadrature<dim - 1> &                 q_boundary,
-      const bool                                  project_to_boundary_first)
-    {
-      switch (dof.get_fe().degree)
-        {
-          case 1:
-            project_matrix_free<components, 1>(mapping,
-                                               dof,
-                                               constraints,
-                                               quadrature,
-                                               function,
-                                               work_result,
-                                               enforce_zero_boundary,
-                                               q_boundary,
-                                               project_to_boundary_first);
-            break;
-
-          case 2:
-            project_matrix_free<components, 2>(mapping,
-                                               dof,
-                                               constraints,
-                                               quadrature,
-                                               function,
-                                               work_result,
-                                               enforce_zero_boundary,
-                                               q_boundary,
-                                               project_to_boundary_first);
-            break;
-
-          case 3:
-            project_matrix_free<components, 3>(mapping,
-                                               dof,
-                                               constraints,
-                                               quadrature,
-                                               function,
-                                               work_result,
-                                               enforce_zero_boundary,
-                                               q_boundary,
-                                               project_to_boundary_first);
-            break;
-
-          default:
-            project_matrix_free<components, -1>(mapping,
-                                                dof,
-                                                constraints,
-                                                quadrature,
-                                                function,
-                                                work_result,
-                                                enforce_zero_boundary,
-                                                q_boundary,
-                                                project_to_boundary_first);
-        }
-    }
-
-
-
     // Helper interface for the matrix-free implementation of project().
     // Used to determine the number of components.
     template <int dim, typename Number, int spacedim>
@@ -351,51 +268,51 @@ namespace VectorTools
       switch (dof.get_fe(0).n_components())
         {
           case 1:
-            project_matrix_free_degree<1>(mapping,
-                                          dof,
-                                          constraints,
-                                          quadrature,
-                                          function,
-                                          work_result,
-                                          enforce_zero_boundary,
-                                          q_boundary,
-                                          project_to_boundary_first);
+            project_matrix_free<1>(mapping,
+                                   dof,
+                                   constraints,
+                                   quadrature,
+                                   function,
+                                   work_result,
+                                   enforce_zero_boundary,
+                                   q_boundary,
+                                   project_to_boundary_first);
             break;
 
           case 2:
-            project_matrix_free_degree<2>(mapping,
-                                          dof,
-                                          constraints,
-                                          quadrature,
-                                          function,
-                                          work_result,
-                                          enforce_zero_boundary,
-                                          q_boundary,
-                                          project_to_boundary_first);
+            project_matrix_free<2>(mapping,
+                                   dof,
+                                   constraints,
+                                   quadrature,
+                                   function,
+                                   work_result,
+                                   enforce_zero_boundary,
+                                   q_boundary,
+                                   project_to_boundary_first);
             break;
 
           case 3:
-            project_matrix_free_degree<3>(mapping,
-                                          dof,
-                                          constraints,
-                                          quadrature,
-                                          function,
-                                          work_result,
-                                          enforce_zero_boundary,
-                                          q_boundary,
-                                          project_to_boundary_first);
+            project_matrix_free<3>(mapping,
+                                   dof,
+                                   constraints,
+                                   quadrature,
+                                   function,
+                                   work_result,
+                                   enforce_zero_boundary,
+                                   q_boundary,
+                                   project_to_boundary_first);
             break;
 
           case 4:
-            project_matrix_free_degree<4>(mapping,
-                                          dof,
-                                          constraints,
-                                          quadrature,
-                                          function,
-                                          work_result,
-                                          enforce_zero_boundary,
-                                          q_boundary,
-                                          project_to_boundary_first);
+            project_matrix_free<4>(mapping,
+                                   dof,
+                                   constraints,
+                                   quadrature,
+                                   function,
+                                   work_result,
+                                   enforce_zero_boundary,
+                                   q_boundary,
+                                   project_to_boundary_first);
             break;
 
           default:
@@ -581,7 +498,7 @@ namespace VectorTools
                                                            vec_result);
     }
 
-    template <int dim, typename VectorType, int spacedim, int fe_degree>
+    template <int dim, typename VectorType, int spacedim>
     void
     project_parallel(
       const Mapping<dim, spacedim> &                            mapping,
@@ -598,9 +515,6 @@ namespace VectorTools
              ExcDimensionMismatch(dof.get_fe(0).n_components(), 1));
       Assert(vec_result.size() == dof.n_dofs(),
              ExcDimensionMismatch(vec_result.size(), dof.n_dofs()));
-      Assert(fe_degree == -1 ||
-               dof.get_fe().degree == static_cast<unsigned int>(fe_degree),
-             ExcDimensionMismatch(fe_degree, dof.get_fe().degree));
 
       // set up mass matrix and right hand side
       typename MatrixFree<dim, Number>::AdditionalData additional_data;
@@ -615,12 +529,8 @@ namespace VectorTools
                           constraints,
                           QGauss<1>(dof.get_fe().degree + 2),
                           additional_data);
-      using MatrixType = MatrixFreeOperators::MassOperator<
-        dim,
-        fe_degree,
-        fe_degree + 2,
-        1,
-        LinearAlgebra::distributed::Vector<Number>>;
+      using MatrixType = MatrixFreeOperators::
+        MassOperator<dim, -1, 0, 1, LinearAlgebra::distributed::Vector<Number>>;
       MatrixType mass_matrix;
       mass_matrix.initialize(matrix_free);
       mass_matrix.compute_diagonal();
@@ -696,11 +606,7 @@ namespace VectorTools
 
 
 
-    template <int dim,
-              typename VectorType,
-              int spacedim,
-              int fe_degree,
-              int n_q_points_1d>
+    template <int dim, typename VectorType, int spacedim>
     void
     project_parallel(
       std::shared_ptr<const MatrixFree<dim, typename VectorType::value_type>>
@@ -720,16 +626,9 @@ namespace VectorTools
              ExcDimensionMismatch(dof.get_fe(0).n_components(), 1));
       Assert(vec_result.size() == dof.n_dofs(),
              ExcDimensionMismatch(vec_result.size(), dof.n_dofs()));
-      Assert(fe_degree == -1 ||
-               dof.get_fe().degree == static_cast<unsigned int>(fe_degree),
-             ExcDimensionMismatch(fe_degree, dof.get_fe().degree));
 
-      using MatrixType = MatrixFreeOperators::MassOperator<
-        dim,
-        fe_degree,
-        n_q_points_1d,
-        1,
-        LinearAlgebra::distributed::Vector<Number>>;
+      using MatrixType = MatrixFreeOperators::
+        MassOperator<dim, -1, 0, 1, LinearAlgebra::distributed::Vector<Number>>;
       MatrixType mass_matrix;
       mass_matrix.initialize(matrix_free, {fe_component});
       mass_matrix.compute_diagonal();
@@ -744,9 +643,8 @@ namespace VectorTools
 
       // assemble right hand side:
       {
-        FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number> fe_eval(
-          *matrix_free, fe_component);
-        const unsigned int n_cells    = matrix_free->n_macro_cells();
+        FEEvaluation<dim, -1, 0, 1, Number> fe_eval(*matrix_free, fe_component);
+        const unsigned int n_cells    = matrix_free->n_cell_batches();
         const unsigned int n_q_points = fe_eval.n_q_points;
 
         for (unsigned int cell = 0; cell < n_cells; ++cell)
@@ -857,24 +755,8 @@ namespace VectorTools
             const unsigned int)> &                                  func,
           VectorType &                                              vec_result)
   {
-    switch (dof.get_fe().degree)
-      {
-        case 1:
-          internal::project_parallel<dim, VectorType, spacedim, 1>(
-            mapping, dof, constraints, quadrature, func, vec_result);
-          break;
-        case 2:
-          internal::project_parallel<dim, VectorType, spacedim, 2>(
-            mapping, dof, constraints, quadrature, func, vec_result);
-          break;
-        case 3:
-          internal::project_parallel<dim, VectorType, spacedim, 3>(
-            mapping, dof, constraints, quadrature, func, vec_result);
-          break;
-        default:
-          internal::project_parallel<dim, VectorType, spacedim, -1>(
-            mapping, dof, constraints, quadrature, func, vec_result);
-      }
+    internal::project_parallel<dim, VectorType, spacedim>(
+      mapping, dof, constraints, quadrature, func, vec_result);
   }
 
 
@@ -886,38 +768,15 @@ namespace VectorTools
             typename VectorType::value_type,
             VectorizedArray<typename VectorType::value_type>>>      matrix_free,
           const AffineConstraints<typename VectorType::value_type> &constraints,
-          const unsigned int      n_q_points_1d,
+          const unsigned int,
           const std::function<VectorizedArray<typename VectorType::value_type>(
             const unsigned int,
             const unsigned int)> &func,
           VectorType &            vec_result,
           const unsigned int      fe_component)
   {
-    const unsigned int fe_degree =
-      matrix_free->get_dof_handler(fe_component).get_fe().degree;
-
-    if (fe_degree + 1 == n_q_points_1d)
-      switch (fe_degree)
-        {
-          case 1:
-            internal::project_parallel<dim, VectorType, dim, 1, 2>(
-              matrix_free, constraints, func, vec_result, fe_component);
-            break;
-          case 2:
-            internal::project_parallel<dim, VectorType, dim, 2, 3>(
-              matrix_free, constraints, func, vec_result, fe_component);
-            break;
-          case 3:
-            internal::project_parallel<dim, VectorType, dim, 3, 4>(
-              matrix_free, constraints, func, vec_result, fe_component);
-            break;
-          default:
-            internal::project_parallel<dim, VectorType, dim, -1, 0>(
-              matrix_free, constraints, func, vec_result, fe_component);
-        }
-    else
-      internal::project_parallel<dim, VectorType, dim, -1, 0>(
-        matrix_free, constraints, func, vec_result, fe_component);
+    internal::project_parallel<dim, VectorType, dim>(
+      matrix_free, constraints, func, vec_result, fe_component);
   }
 
 
diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index fa08b3466a..3d878bbc04 100644
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -36,8 +36,7 @@
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/tensor_product_matrix.h>
 
-#include <deal.II/matrix_free/evaluation_kernels.h>
-#include <deal.II/matrix_free/evaluation_selector.h>
+#include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 
@@ -1541,15 +1540,17 @@ namespace internal
                 }
 
             // do the actual tensorized evaluation
-            SelectEvaluator<dim, -1, 0, VectorizedArray<double>>::evaluate(
-              n_comp,
-              evaluation_flag,
-              data.shape_info,
-              data.values_dofs.begin(),
-              data.values_quad.begin(),
-              data.gradients_quad.begin(),
-              data.hessians_quad.begin(),
-              data.scratch.begin());
+            internal::FEEvaluationFactory<
+              dim,
+              double,
+              VectorizedArray<double>>::evaluate(n_comp,
+                                                 evaluation_flag,
+                                                 data.shape_info,
+                                                 data.values_dofs.begin(),
+                                                 data.values_quad.begin(),
+                                                 data.gradients_quad.begin(),
+                                                 data.hessians_quad.begin(),
+                                                 data.scratch.begin());
           }
 
         // do the postprocessing