From 61c1e0f93b52bbee0686db884aaa7796c3ffdee9 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 5 Aug 2023 14:44:42 +0200
Subject: [PATCH] FEPointEvaluation/MappingInfo for simplices

---
 .../deal.II/matrix_free/fe_point_evaluation.h | 26 ++++++++
 include/deal.II/non_matching/mapping_info.h   | 59 ++++++++++++++++++-
 2 files changed, 84 insertions(+), 1 deletion(-)

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index 4f71313c46..841adb4725 100644
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -1603,6 +1603,32 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
   current_face_number = numbers::invalid_unsigned_int;
 
   do_reinit();
+
+  if (!fast_path)
+    {
+      const auto unit_points_vectorized = mapping_info->get_unit_point(
+        mapping_info->compute_unit_point_index_offset(current_cell_index,
+                                                      current_face_number));
+      const unsigned int n_q_points_unvectorized =
+        mapping_info->get_n_q_points_unvectorized(current_cell_index,
+                                                  current_face_number);
+
+      std::vector<Point<dim>> unit_points(n_q_points_unvectorized);
+
+      for (unsigned int v = 0; v < n_q_points_unvectorized; ++v)
+        for (unsigned int d = 0; d < dim; ++d)
+          unit_points[v][d] = unit_points_vectorized[v / n_lanes_internal][d]
+                                                    [v % n_lanes_internal];
+
+      fe_values = std::make_shared<FEValues<dim, spacedim>>(
+        *mapping,
+        *fe,
+        Quadrature<dim>(
+          std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+        update_flags);
+
+      fe_values->reinit(mapping_info->get_cell_iterator(current_cell_index));
+    }
 }
 
 
diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h
index 532dbe635b..b8c571fd95 100644
--- a/include/deal.II/non_matching/mapping_info.h
+++ b/include/deal.II/non_matching/mapping_info.h
@@ -216,8 +216,10 @@ namespace NonMatching
       /**
        * Constructor which sets the default arguments.
        */
-      AdditionalData(const bool use_global_weights = false)
+      AdditionalData(const bool use_global_weights = false,
+                     const bool store_cells        = false)
         : use_global_weights(use_global_weights)
+        , store_cells(store_cells)
       {}
 
       /**
@@ -226,6 +228,15 @@ namespace NonMatching
        * QSimplex::compute_affine_transformation().
        */
       bool use_global_weights;
+
+      /**
+       * During the reinit() function calls, cells are passed as
+       * argument. In the default case, the cell is not stored,
+       * since all relevant mapping related information is precomputed.
+       * Hoever, this flag enables that the cells are stored so that
+       * they can be accessed later on.
+       */
+      bool store_cells;
     };
 
     /**
@@ -422,6 +433,14 @@ namespace NonMatching
     get_n_q_points_unvectorized(const unsigned int cell_index,
                                 const unsigned int face_number) const;
 
+    /**
+     * Return cell iterator.
+     *
+     * @note This call is only possible if AdditionalData::store_cells is enabled.
+     */
+    typename Triangulation<dim, spacedim>::cell_iterator
+    get_cell_iterator(const unsigned int cell_index) const;
+
   private:
     using MappingData =
       dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
@@ -633,6 +652,19 @@ namespace NonMatching
      * dependent objects know that they need to reinitialize as well.
      */
     boost::signals2::signal<void()> is_reinitialized;
+
+    /**
+     * Reference to the triangulation passed via the cells to the
+     * reinit functions. This field is only set if
+     * AdditionalData::store_cells is enabled.
+     */
+    SmartPointer<const Triangulation<dim, spacedim>> triangulation;
+
+    /**
+     * Level and indices of cells passed to the reinit functions. This
+     * vector is only filled if AdditionalData::store_cells is enabled.
+     */
+    std::vector<std::pair<int, int>> cell_level_and_indices;
   };
 
   // ----------------------- template functions ----------------------
@@ -807,6 +839,9 @@ namespace NonMatching
 
     n_q_points_unvectorized.reserve(n_cells);
 
+    if (additional_data.store_cells)
+      cell_level_and_indices.resize(n_cells);
+
     // fill unit points index offset vector
     unit_points_index.reserve(n_cells + 1);
     unit_points_index.push_back(0);
@@ -844,6 +879,12 @@ namespace NonMatching
     unsigned int cell_index = 0;
     for (const auto &cell : cell_iterator_range)
       {
+        if (additional_data.store_cells)
+          {
+            this->triangulation                = &cell->get_triangulation();
+            cell_level_and_indices[cell_index] = {cell->level(), cell->index()};
+          }
+
         // store unit points
         const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
           n_q_points_unvectorized[cell_index]);
@@ -1164,6 +1205,22 @@ namespace NonMatching
 
 
 
+  template <int dim, int spacedim, typename Number>
+  typename Triangulation<dim, spacedim>::cell_iterator
+  MappingInfo<dim, spacedim, Number>::get_cell_iterator(
+    const unsigned int cell_index) const
+  {
+    Assert(
+      additional_data.store_cells,
+      ExcMessage(
+        "Cells have been not stored. You can enable this by Additional::store_cells."));
+    return {triangulation.get(),
+            cell_level_and_indices[cell_index].first,
+            cell_level_and_indices[cell_index].second};
+  }
+
+
+
   template <int dim, int spacedim, typename Number>
   template <typename NumberType>
   unsigned int
-- 
2.39.5