Add helper class RemotePointEvaluation::CellData

diff --git a/doc/news/changes/minor/20230917SchreterMunch b/doc/news/changes/minor/20230917SchreterMunch
new file mode 100644 (file)
index 0000000..3cadc71
--- /dev/null
+++ b/doc/news/changes/minor/20230917SchreterMunch
@@ -0,0 +1,4 @@
+New: Added helper class RemotePointEvaluation::CellData to store
+and access data of cell-specific points.
+<br>
+(Magdalena Schreter, Peter Munch, 2023/09/17)

diff --git a/include/deal.II/base/mpi_remote_point_evaluation.h b/include/deal.II/base/mpi_remote_point_evaluation.h
index af7de524372016392f49bbbe54d9e0f28ced6c36..0c441c2854a2b5b92b210c07ac7f127a9ed4bbde 100644 (file)
--- a/include/deal.II/base/mpi_remote_point_evaluation.h
+++ b/include/deal.II/base/mpi_remote_point_evaluation.h
@@ -111,17 +111,74 @@ namespace Utilities
              const Mapping<dim, spacedim> &mapping);
 
       /**
-       * Data of points positioned in a cell.
+       * Helper class to store and to access data of points positioned in
+       * processed cells.
        */
-      struct CellData
+      class CellData
       {
+      public:
         /**
-         * Level and index of cells.
+         * Constructor.
+         *
+         * @param triangulation Triangulation of the domain.
+         */
+        CellData(const Triangulation<dim, spacedim> &triangulation);
+
+        /**
+         * Return an object that can be thought of as an array containing all
+         * indices from zero (inclusive) to the total number of cells
+         * where points reside given by `cells.size()` (exclusive). This allows
+         * one to write code using range-based `for` loops of the following
+         * kind:
+         *
+         * @code
+         * for (const auto cell : cell_data.cell_indices())
+         *   {
+         *     const auto cell_dofs =
+         *       cell_data.get_active_cell_iterator(cell)->as_dof_handler_iterator(
+         *         dof_handler);
+         *
+         *     const auto unit_points = cell_data.get_unit_points(cell);
+         *     const auto local_value = cell_data.get_data_view(cell, values);
+         *
+         *     // user code: not shown
+         *   }
+         * @endcode
+         *
+         * Here, we are looping over all cells where points reside and
+         * use the index to call the functions get_active_cell_iterator(),
+         * get_unit_points(), and get_data_view().
+         */
+        std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+        cell_indices() const;
+
+        /**
+         * Return active cell iterator of the processed cell @p cell.
+         */
+        typename Triangulation<dim, spacedim>::active_cell_iterator
+        get_active_cell_iterator(const unsigned int cell) const;
+
+        /**
+         * Return unit points of the processed cell @p cell.
+         */
+        ArrayView<const Point<dim>>
+        get_unit_points(const unsigned int cell) const;
+
+        /**
+         * Return local view of the processed cell @p cell for the vector @p values.
+         */
+        template <typename T>
+        ArrayView<T>
+        get_data_view(const unsigned int  cell,
+                      const ArrayView<T> &values) const;
+
+        /**
+         * Level and index of processed cells.
          */
         std::vector<std::pair<int, int>> cells;
 
         /**
-         * Pointers to beginning and ending of the (reference) points
+         * Pointers to the start and end of the (reference) points
          * associated to the cell.
          */
         std::vector<unsigned int> reference_point_ptrs;
@@ -130,6 +187,13 @@ namespace Utilities
          * Reference points in the interval [0,1]^dim.
          */
         std::vector<Point<dim>> reference_point_values;
+
+      private:
+        /**
+         * Reference to the underlying triangulation needed for
+         * get_active_cell_iterator().
+         */
+        const Triangulation<dim, spacedim> &triangulation;
       };
 
       /**
@@ -327,9 +391,9 @@ namespace Utilities
 
       /**
        * Point data sorted according to cells so that evaluation (incl. reading
-       * of degrees of freedoms) needs to performed only once per cell.
+       * of degrees of freedoms) needs to be performed only once per cell.
        */
-      CellData cell_data;
+      std::unique_ptr<CellData> cell_data;
 
       /**
        * Permutation index within a send buffer.
@@ -349,6 +413,21 @@ namespace Utilities
     };
 
 
+
+    template <int dim, int spacedim>
+    template <typename T>
+    ArrayView<T>
+    RemotePointEvaluation<dim, spacedim>::CellData::get_data_view(
+      const unsigned int  cell,
+      const ArrayView<T> &values) const
+    {
+      AssertIndexRange(cell, cells.size());
+      return {values.data() + reference_point_ptrs[cell],
+              reference_point_ptrs[cell + 1] - reference_point_ptrs[cell]};
+    }
+
+
+
     template <int dim, int spacedim>
     template <typename T>
     void
@@ -383,7 +462,7 @@ namespace Utilities
                                send_permutation.size());
 
       // evaluate functions at points
-      evaluation_function(buffer_eval, cell_data);
+      evaluation_function(buffer_eval, *cell_data);
 
       // sort for communication
       unsigned int my_rank_local_recv = numbers::invalid_unsigned_int;
@@ -694,7 +773,7 @@ namespace Utilities
         }
 
       // evaluate function at points
-      evaluation_function(buffer_eval, cell_data);
+      evaluation_function(buffer_eval, *cell_data);
 #endif
     }
 

diff --git a/source/base/mpi_remote_point_evaluation.cc b/source/base/mpi_remote_point_evaluation.cc
index 07f67dfe540e1fe5d5868da6d7d7c5f59e786722..e0c52cceffd30f3caf72fd83298bba4e6fab4d4d 100644 (file)
--- a/source/base/mpi_remote_point_evaluation.cc
+++ b/source/base/mpi_remote_point_evaluation.cc
@@ -188,7 +188,7 @@ namespace Utilities
       Assert(enforce_unique_mapping == false || unique_mapping,
              ExcInternalError());
 
-      cell_data        = {};
+      cell_data        = std::make_unique<CellData>(tria);
       send_permutation = {};
 
       std::pair<int, int> dummy{-1, -1};
@@ -197,28 +197,68 @@ namespace Utilities
           if (dummy != std::get<0>(i))
             {
               dummy = std::get<0>(i);
-              cell_data.cells.emplace_back(dummy);
-              cell_data.reference_point_ptrs.emplace_back(
-                cell_data.reference_point_values.size());
+              cell_data->cells.emplace_back(dummy);
+              cell_data->reference_point_ptrs.emplace_back(
+                cell_data->reference_point_values.size());
             }
 
-          cell_data.reference_point_values.emplace_back(std::get<3>(i));
+          cell_data->reference_point_values.emplace_back(std::get<3>(i));
           send_permutation.emplace_back(std::get<5>(i));
         }
 
-      cell_data.reference_point_ptrs.emplace_back(
-        cell_data.reference_point_values.size());
+      cell_data->reference_point_ptrs.emplace_back(
+        cell_data->reference_point_values.size());
 
       this->ready_flag = true;
     }
 
 
 
+    template <int dim, int spacedim>
+    RemotePointEvaluation<dim, spacedim>::CellData::CellData(
+      const Triangulation<dim, spacedim> &triangulation)
+      : triangulation(triangulation)
+    {}
+
+
+
+    template <int dim, int spacedim>
+    std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+    RemotePointEvaluation<dim, spacedim>::CellData::cell_indices() const
+    {
+      return {0, static_cast<unsigned int>(cells.size())};
+    }
+
+
+
+    template <int dim, int spacedim>
+    typename Triangulation<dim, spacedim>::active_cell_iterator
+    RemotePointEvaluation<dim, spacedim>::CellData::get_active_cell_iterator(
+      const unsigned int cell) const
+    {
+      AssertIndexRange(cell, cells.size());
+      return {&triangulation, cells[cell].first, cells[cell].second};
+    }
+
+
+
+    template <int dim, int spacedim>
+    ArrayView<const Point<dim>>
+    RemotePointEvaluation<dim, spacedim>::CellData::get_unit_points(
+      const unsigned int cell) const
+    {
+      AssertIndexRange(cell, cells.size());
+      return {reference_point_values.data() + reference_point_ptrs[cell],
+              reference_point_ptrs[cell + 1] - reference_point_ptrs[cell]};
+    }
+
+
+
     template <int dim, int spacedim>
     const typename RemotePointEvaluation<dim, spacedim>::CellData &
     RemotePointEvaluation<dim, spacedim>::get_cell_data() const
     {
-      return cell_data;
+      return *cell_data;
     }