From 5c2b8eea9376ff622a85c07d85cda2602227a3ab Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 6 Feb 2022 07:52:04 +0100
Subject: [PATCH] Make FEEval more consistent with FEValues

---
 doc/news/changes/minor/20220206Munch          |  6 +++
 include/deal.II/matrix_free/fe_evaluation.h   | 54 +++++++++++++++++++
 .../deal.II/matrix_free/fe_evaluation_data.h  | 17 ++++++
 tests/matrix_free/ecl_02.cc                   | 12 ++---
 4 files changed, 83 insertions(+), 6 deletions(-)
 create mode 100644 doc/news/changes/minor/20220206Munch

diff --git a/doc/news/changes/minor/20220206Munch b/doc/news/changes/minor/20220206Munch
new file mode 100644
index 0000000000..e78837613b
--- /dev/null
+++ b/doc/news/changes/minor/20220206Munch
@@ -0,0 +1,6 @@
+Improved: To be consistent with the FEValues classes, one can now 
+loop over all DoF and quadrature-point indices in a 
+range-based loop with the help of the new functions dof_indices() and
+quadrature_point_indices().
+<br>
+(Peter Munch, 2022/02/06)
diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index b4c4ab72ff..16e236948b 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -2302,6 +2302,14 @@ public:
                     const bool  integrate_gradients,
                     VectorType &output_vector);
 
+  /**
+   * Return an object that can be thought of as an array containing all indices
+   * from zero to @p dofs_per_cell. This allows to write code using
+   * range-based for loops.
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  dof_indices() const;
+
   /**
    * The number of degrees of freedom of a single component on the cell for
    * the underlying evaluation object. Usually close to
@@ -2700,6 +2708,14 @@ public:
                     const bool  integrate_gradients,
                     VectorType &output_vector);
 
+  /**
+   * Return an object that can be thought of as an array containing all indices
+   * from zero to @p dofs_per_cell. This allows to write code using
+   * range-based for loops.
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  dof_indices() const;
+
   /**
    * The number of degrees of freedom of a single component on the cell for
    * the underlying evaluation object. Usually close to
@@ -7552,6 +7568,25 @@ FEEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+FEEvaluation<dim,
+             fe_degree,
+             n_q_points_1d,
+             n_components_,
+             Number,
+             VectorizedArrayType>::dof_indices() const
+{
+  return {0U, dofs_per_cell};
+}
+
+
+
 /*-------------------------- FEFaceEvaluation ---------------------------*/
 
 
@@ -8351,6 +8386,25 @@ FEFaceEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+FEFaceEvaluation<dim,
+                 fe_degree,
+                 n_q_points_1d,
+                 n_components_,
+                 Number,
+                 VectorizedArrayType>::dof_indices() const
+{
+  return {0U, dofs_per_cell};
+}
+
+
+
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h
index 57ff0ac602..26df97bf9e 100644
--- a/include/deal.II/matrix_free/fe_evaluation_data.h
+++ b/include/deal.II/matrix_free/fe_evaluation_data.h
@@ -533,6 +533,14 @@ public:
       return face_ids;
   }
 
+  /**
+   * Return an object that can be thought of as an array containing all indices
+   * from zero to @p n_quadrature_points. This allows to write code using
+   * range-based for loops.
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  quadrature_point_indices() const;
+
   //@}
 
   /**
@@ -1582,6 +1590,15 @@ FEEvaluationData<dim, Number, is_face>::is_interior_face() const
 
 
 
+template <int dim, typename Number, bool is_face>
+inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+FEEvaluationData<dim, Number, is_face>::quadrature_point_indices() const
+{
+  return {0U, n_quadrature_points};
+}
+
+
+
 namespace internal
 {
   template <std::size_t N,
diff --git a/tests/matrix_free/ecl_02.cc b/tests/matrix_free/ecl_02.cc
index 04a699695d..2785699537 100644
--- a/tests/matrix_free/ecl_02.cc
+++ b/tests/matrix_free/ecl_02.cc
@@ -125,24 +125,24 @@ test(const unsigned int n_refinements = 1)
 
               phi_m.read_dof_values(src);
 
-              for (unsigned int i = 0; i < phi_m.static_dofs_per_component; ++i)
+              for (const auto i : phi_m.dof_indices())
                 deallog << static_cast<int>(phi_m.begin_dof_values()[i][0])
                         << ' ';
               deallog << std::endl;
               phi_m.gather_evaluate(src, EvaluationFlags::values);
-              for (unsigned int i = 0; i < phi_m.static_n_q_points; ++i)
-                deallog << static_cast<int>(phi_m.begin_values()[i][0]) << ' ';
+              for (const auto q : phi_p.quadrature_point_indices())
+                deallog << static_cast<int>(phi_m.begin_values()[q][0]) << " ";
               deallog << std::endl;
 
               phi_p.read_dof_values(src);
-              for (unsigned int i = 0; i < phi_p.static_dofs_per_component; ++i)
+              for (const auto i : phi_p.dof_indices())
                 deallog << static_cast<int>(phi_p.begin_dof_values()[i][0])
                         << ' ';
               deallog << std::endl;
 
               phi_p.gather_evaluate(src, EvaluationFlags::values);
-              for (unsigned int i = 0; i < phi_p.static_n_q_points; ++i)
-                deallog << static_cast<int>(phi_p.begin_values()[i][0]) << ' ';
+              for (const auto q : phi_p.quadrature_point_indices())
+                deallog << static_cast<int>(phi_p.begin_values()[q][0]) << " ";
               deallog << std::endl << std::endl;
             }
         }
-- 
2.39.5