From 3afb70dfc68c38559197d8691eb0ffd225e22d7a Mon Sep 17 00:00:00 2001
From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Wed, 24 Nov 2021 14:28:11 -0500
Subject: [PATCH] Separate mapping related computations from FEPointEvaluation

---
 .../deal.II/matrix_free/fe_point_evaluation.h | 105 +++++++++---------
 source/matrix_free/fe_point_evaluation.cc     |  21 ++++
 .../matrix_free/fe_point_evaluation.inst.in   |   2 +
 3 files changed, 75 insertions(+), 53 deletions(-)

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index 0ab23f57f8..ddd9007e0e 100644
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -357,13 +357,24 @@ namespace internal
     is_fast_path_supported(const FiniteElement<dim, spacedim> &fe,
                            const unsigned int base_element_number);
 
+    template <int dim, int spacedim>
+    bool
+    is_fast_path_supported(const Mapping<dim, spacedim> &mapping);
+
     template <int dim, int spacedim>
     std::vector<Polynomials::Polynomial<double>>
     get_polynomial_space(const FiniteElement<dim, spacedim> &fe);
   } // namespace FEPointEvaluation
 } // namespace internal
 
-
+template <int dim, int spacedim>
+void
+compute_mapping_data(
+  const Mapping<dim> &                                        mapping,
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const ArrayView<const Point<dim>> &                         unit_points,
+  dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &mapping_data);
 
 /**
  * This class provides an interface to the evaluation of interpolated solution
@@ -615,18 +626,6 @@ private:
    */
   SmartPointer<const Mapping<dim, spacedim>> mapping;
 
-  /**
-   * Pointer to the function of the mapping that computes the necessary
-   * mapping quantities like quadrature points and Jacobians.
-   */
-  std::function<void(
-    const typename Triangulation<dim, spacedim>::cell_iterator & /*cell*/,
-    const ArrayView<const Point<dim>> & /*unit_points*/,
-    const UpdateFlags /*update_flags*/,
-    dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-      & /*output_data*/)>
-    fill_mapping_data_for_generic_points;
-
   /**
    * Pointer to the FiniteElement object passed to the constructor.
    */
@@ -747,41 +746,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
   AssertIndexRange(first_selected_component + n_components,
                    fe.n_components() + 1);
 
-  if (const MappingQ<dim, spacedim> *mapping_q =
-        dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
-    {
-      fill_mapping_data_for_generic_points =
-        [mapping_q](
-          const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-          const ArrayView<const Point<dim>> &unit_points,
-          const UpdateFlags                  update_flags,
-          dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                       spacedim>
-            &output_data) -> void {
-        mapping_q->fill_mapping_data_for_generic_points(cell,
-                                                        unit_points,
-                                                        update_flags,
-                                                        output_data);
-      };
-    }
-  else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
-             dynamic_cast<const MappingCartesian<dim, spacedim> *>(&mapping))
-    {
-      fill_mapping_data_for_generic_points =
-        [mapping_cartesian](
-          const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-          const ArrayView<const Point<dim>> &unit_points,
-          const UpdateFlags                  update_flags,
-          dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                       spacedim>
-            &output_data) -> void {
-        mapping_cartesian->fill_mapping_data_for_generic_points(cell,
-                                                                unit_points,
-                                                                update_flags,
-                                                                output_data);
-      };
-    }
-
   bool         same_base_element   = true;
   unsigned int base_element_number = 0;
   component_in_base_element        = 0;
@@ -799,7 +763,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
     else
       component += fe.element_multiplicity(base_element_number);
 
-  if (fill_mapping_data_for_generic_points &&
+  if (internal::FEPointEvaluation::is_fast_path_supported(mapping) &&
       internal::FEPointEvaluation::is_fast_path_supported(
         fe, base_element_number) &&
       same_base_element)
@@ -847,6 +811,43 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
 }
 
 
+template <int dim, int spacedim>
+void
+compute_mapping_data_for_generic_points(
+  const Mapping<dim> &                                        mapping,
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const ArrayView<const Point<dim>> &                         unit_points,
+  const UpdateFlags                                           update_flags,
+  dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &mapping_data)
+{
+  UpdateFlags update_flags_mapping = update_default;
+  // translate update flags
+  if (update_flags & update_jacobians)
+    update_flags_mapping |= update_jacobians;
+  if (update_flags & update_gradients ||
+      update_flags & update_inverse_jacobians)
+    update_flags_mapping |= update_inverse_jacobians;
+  if (update_flags & update_quadrature_points)
+    update_flags_mapping |= update_quadrature_points;
+
+  if (const MappingQ<dim, spacedim> *mapping_q =
+        dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
+    {
+      mapping_q->fill_mapping_data_for_generic_points(cell,
+                                                      unit_points,
+                                                      update_flags_mapping,
+                                                      mapping_data);
+    }
+  else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
+             dynamic_cast<const MappingCartesian<dim, spacedim> *>(&mapping))
+    {
+      mapping_cartesian->fill_mapping_data_for_generic_points(
+        cell, unit_points, update_flags_mapping, mapping_data);
+    }
+}
+
+
 
 template <int n_components, int dim, int spacedim, typename Number>
 void
@@ -856,10 +857,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
 {
   // If using the fast path, we need to precompute the mapping data.
   if (!poly.empty())
-    fill_mapping_data_for_generic_points(cell,
-                                         unit_points,
-                                         update_flags_mapping,
-                                         mapping_data);
+    compute_mapping_data_for_generic_points(
+      *mapping, cell, unit_points, update_flags_mapping, mapping_data);
 
   // then call the other version of this function with the precomputed data
   reinit(cell, unit_points, mapping_data);
diff --git a/source/matrix_free/fe_point_evaluation.cc b/source/matrix_free/fe_point_evaluation.cc
index 7ecde1907a..3770f18435 100644
--- a/source/matrix_free/fe_point_evaluation.cc
+++ b/source/matrix_free/fe_point_evaluation.cc
@@ -78,6 +78,27 @@ namespace internal
     }
 
 
+
+    template <int dim, int spacedim>
+    bool
+    is_fast_path_supported(const Mapping<dim, spacedim> &mapping)
+    {
+      if (const MappingQ<dim, spacedim> *mapping_q =
+            dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
+        {
+          return true;
+        }
+      else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
+                 dynamic_cast<const MappingCartesian<dim, spacedim> *>(
+                   &mapping))
+        {
+          return true;
+        }
+      return false;
+    }
+
+
+
     template <int dim, int spacedim>
     std::vector<Polynomials::Polynomial<double>>
     get_polynomial_space(const FiniteElement<dim, spacedim> &fe)
diff --git a/source/matrix_free/fe_point_evaluation.inst.in b/source/matrix_free/fe_point_evaluation.inst.in
index 74736be9ff..6a50a3fc68 100644
--- a/source/matrix_free/fe_point_evaluation.inst.in
+++ b/source/matrix_free/fe_point_evaluation.inst.in
@@ -21,6 +21,8 @@ for (deal_II_dimension, deal_II_space_dimension : DIMENSIONS)
     template bool internal::FEPointEvaluation::is_fast_path_supported(
       const FiniteElement<deal_II_dimension, deal_II_space_dimension> &,
       const unsigned int);
+    template bool internal::FEPointEvaluation::is_fast_path_supported(
+      const Mapping<deal_II_dimension, deal_II_space_dimension> &);
     template std::vector<Polynomials::Polynomial<double>>
     internal::FEPointEvaluation::get_polynomial_space(
       const FiniteElement<deal_II_dimension, deal_II_space_dimension> &);
-- 
2.39.5