From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Tue, 28 Apr 2020 21:55:21 +0000 (+0200)
Subject: Use FEValuesBase::quadrature_point_indices() in fe tests
X-Git-Tag: v9.2.0-rc1~155^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d7ebc8b20b883297b84ed35565047b8f2aba78d;p=dealii.git

Use FEValuesBase::quadrature_point_indices() in fe tests
---

diff --git a/tests/fe/abf_01.cc b/tests/fe/abf_01.cc
index 764f5ae108..8c3210b5fc 100644
--- a/tests/fe/abf_01.cc
+++ b/tests/fe/abf_01.cc
@@ -99,7 +99,7 @@ void EvaluateDerivative(DoFHandler<2> *dof_handler, Vector<double> &solution)
         n_q_points, std::vector<Tensor<1, 2>>(n_components));
       fe_values.get_function_gradients(solution, grads_here);
 
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         {
           //    double u0 = this_value[q_point](0);
           // double v0 = this_value[q_point](1);
@@ -186,7 +186,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             {
               coefficient->value_list(fe_values.get_quadrature_points(),
                                       coefficient_values);
-              for (unsigned int point = 0; point < n_q_points; ++point)
+              for (const auto point : fe_values.quadrature_point_indices())
                 {
                   const double weight = fe_values.JxW(point);
                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -209,7 +209,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             {
               coefficient->vector_value_list(fe_values.get_quadrature_points(),
                                              coefficient_vector_values);
-              for (unsigned int point = 0; point < n_q_points; ++point)
+              for (const auto point : fe_values.quadrature_point_indices())
                 {
                   const double weight = fe_values.JxW(point);
                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -260,7 +260,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             }
           */
 
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (const auto point : fe_values.quadrature_point_indices())
             {
               const double weight = fe_values.JxW(point);
               //      const double weight = q.weight(point);
@@ -368,7 +368,7 @@ create_right_hand_side(const Mapping<dim> &   mapping,
                                   rhs_values);
 
           cell_vector = 0;
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (const auto point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               cell_vector(i) += rhs_values[point] *
                                 fe_values.shape_value(i, point) *
@@ -394,7 +394,7 @@ create_right_hand_side(const Mapping<dim> &   mapping,
                                          rhs_values);
 
           cell_vector = 0;
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (const auto point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               for (unsigned int comp_i = 0; comp_i < fe.n_components();
                    ++comp_i)
diff --git a/tests/fe/abf_02.cc b/tests/fe/abf_02.cc
index b25ee66938..afe2f4ad1e 100644
--- a/tests/fe/abf_02.cc
+++ b/tests/fe/abf_02.cc
@@ -98,7 +98,7 @@ void EvaluateDerivative(DoFHandler<3> &dof_handler, Vector<double> &solution)
         n_q_points, std::vector<Tensor<1, 3>>(n_components));
       fe_values.get_function_gradients(solution, grads_here);
 
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         {
           double u0 = 0;
           double v0 = 0;
@@ -197,7 +197,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             {
               coefficient->value_list(fe_values.get_quadrature_points(),
                                       coefficient_values);
-              for (unsigned int point = 0; point < n_q_points; ++point)
+              for (const auto point : fe_values.quadrature_point_indices())
                 {
                   const double weight = fe_values.JxW(point);
                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -220,7 +220,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             {
               coefficient->vector_value_list(fe_values.get_quadrature_points(),
                                              coefficient_vector_values);
-              for (unsigned int point = 0; point < n_q_points; ++point)
+              for (const auto point : fe_values.quadrature_point_indices())
                 {
                   const double weight = fe_values.JxW(point);
                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -250,7 +250,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
           const unsigned int        dofs_per_face = fe.dofs_per_face;
           std::vector<unsigned int> face_dof_indices(dofs_per_face);
 
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (const auto point : fe_values.quadrature_point_indices())
             {
               const double weight = fe_values.JxW(point);
               //      const double weight = q.weight(point);
@@ -350,7 +350,7 @@ create_right_hand_side(const Mapping<dim> &   mapping,
                                   rhs_values);
 
           cell_vector = 0;
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (const auto point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               cell_vector(i) += rhs_values[point] *
                                 fe_values.shape_value(i, point) *
@@ -376,7 +376,7 @@ create_right_hand_side(const Mapping<dim> &   mapping,
                                          rhs_values);
 
           cell_vector = 0;
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (const auto point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               for (unsigned int comp_i = 0; comp_i < fe.n_components();
                    ++comp_i)
diff --git a/tests/fe/br_approximation_01.cc b/tests/fe/br_approximation_01.cc
index 0e4758d858..365bb002e3 100644
--- a/tests/fe/br_approximation_01.cc
+++ b/tests/fe/br_approximation_01.cc
@@ -376,7 +376,7 @@ double TestProjection(Mapping<2> &mapping, DoFHandler<2> *dof_handler)
                                                  Vector<double>(n_components));
           fe_values.get_function_values(solution, this_value);
 
-          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (const auto q_point : fe_values.quadrature_point_indices())
             {
               double   u = this_value[q_point](0), v = this_value[q_point](1);
               Point<2> p = fe_values.quadrature_point(q_point);
diff --git a/tests/fe/cell_similarity_01.cc b/tests/fe/cell_similarity_01.cc
index bc076c02a7..81f60f915a 100644
--- a/tests/fe/cell_similarity_01.cc
+++ b/tests/fe/cell_similarity_01.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_02.cc b/tests/fe/cell_similarity_02.cc
index fa58ded7ef..3f81faed71 100644
--- a/tests/fe/cell_similarity_02.cc
+++ b/tests/fe/cell_similarity_02.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_03.cc b/tests/fe/cell_similarity_03.cc
index a90412198d..85c3caa9e9 100644
--- a/tests/fe/cell_similarity_03.cc
+++ b/tests/fe/cell_similarity_03.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_04.cc b/tests/fe/cell_similarity_04.cc
index f561ed6144..ec3e226664 100644
--- a/tests/fe/cell_similarity_04.cc
+++ b/tests/fe/cell_similarity_04.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_05.cc b/tests/fe/cell_similarity_05.cc
index ed0b67819e..4bb9d30bef 100644
--- a/tests/fe/cell_similarity_05.cc
+++ b/tests/fe/cell_similarity_05.cc
@@ -99,7 +99,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_06.cc b/tests/fe/cell_similarity_06.cc
index 2beeba42eb..69dcb5942e 100644
--- a/tests/fe/cell_similarity_06.cc
+++ b/tests/fe/cell_similarity_06.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_07.cc b/tests/fe/cell_similarity_07.cc
index 8c460603c4..3b0a4fbcdd 100644
--- a/tests/fe/cell_similarity_07.cc
+++ b/tests/fe/cell_similarity_07.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_08.cc b/tests/fe/cell_similarity_08.cc
index d188ccbb62..1ed45d7dec 100644
--- a/tests/fe/cell_similarity_08.cc
+++ b/tests/fe/cell_similarity_08.cc
@@ -99,7 +99,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_09.cc b/tests/fe/cell_similarity_09.cc
index 4546c54813..4a4320d510 100644
--- a/tests/fe/cell_similarity_09.cc
+++ b/tests/fe/cell_similarity_09.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_10.cc b/tests/fe/cell_similarity_10.cc
index 622a48a317..c689b873b8 100644
--- a/tests/fe/cell_similarity_10.cc
+++ b/tests/fe/cell_similarity_10.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_11.cc b/tests/fe/cell_similarity_11.cc
index 59af03ce22..1995cab4cc 100644
--- a/tests/fe/cell_similarity_11.cc
+++ b/tests/fe/cell_similarity_11.cc
@@ -78,7 +78,7 @@ test(const Triangulation<dim> &tr)
       fe_values.reinit(cell);
 
       deallog << "Jacobians: ";
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << "[ ";
           for (unsigned int d = 0; d < dim; ++d)
@@ -88,7 +88,7 @@ test(const Triangulation<dim> &tr)
         }
       deallog << std::endl;
       deallog << "Derivatives of shape function: ";
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << "[ ";
           for (unsigned int d = 0; d < dim; ++d)
diff --git a/tests/fe/cell_similarity_12.cc b/tests/fe/cell_similarity_12.cc
index fe0e482eb4..1135bdf409 100644
--- a/tests/fe/cell_similarity_12.cc
+++ b/tests/fe/cell_similarity_12.cc
@@ -79,7 +79,7 @@ test(const unsigned int degree)
         fe_values.reinit(cell);
 
         deallog << "Jacobians: ";
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             deallog << "[ ";
             for (unsigned int d = 0; d < dim; ++d)
@@ -89,7 +89,7 @@ test(const unsigned int degree)
           }
         deallog << std::endl;
         deallog << "Derivatives of shape function: ";
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             deallog << "[ ";
             for (unsigned int d = 0; d < dim; ++d)
diff --git a/tests/fe/cell_similarity_13.cc b/tests/fe/cell_similarity_13.cc
index 6e47fb7dfc..1312451d01 100644
--- a/tests/fe/cell_similarity_13.cc
+++ b/tests/fe/cell_similarity_13.cc
@@ -72,7 +72,7 @@ test()
       fe_values.reinit(cell);
 
       deallog << "Jacobians: ";
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << "[ ";
           for (unsigned int d = 0; d < dim; ++d)
@@ -82,7 +82,7 @@ test()
         }
       deallog << std::endl;
       deallog << "Derivatives of shape function: ";
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << "[ ";
           for (unsigned int d = 0; d < dim; ++d)
diff --git a/tests/fe/cell_similarity_dgp_monomial_01.cc b/tests/fe/cell_similarity_dgp_monomial_01.cc
index 37716797cb..d9452be4c0 100644
--- a/tests/fe/cell_similarity_dgp_monomial_01.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_01.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_02.cc b/tests/fe/cell_similarity_dgp_monomial_02.cc
index 0f7b51e11f..b40ade2469 100644
--- a/tests/fe/cell_similarity_dgp_monomial_02.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_02.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_03.cc b/tests/fe/cell_similarity_dgp_monomial_03.cc
index ccac23ea31..78a75cfcdc 100644
--- a/tests/fe/cell_similarity_dgp_monomial_03.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_03.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_04.cc b/tests/fe/cell_similarity_dgp_monomial_04.cc
index 0cab031577..7e55929402 100644
--- a/tests/fe/cell_similarity_dgp_monomial_04.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_04.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_05.cc b/tests/fe/cell_similarity_dgp_monomial_05.cc
index 13dd9c027d..51ae79c985 100644
--- a/tests/fe/cell_similarity_dgp_monomial_05.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_05.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_06.cc b/tests/fe/cell_similarity_dgp_monomial_06.cc
index 4a479c5c9f..5161000f2c 100644
--- a/tests/fe/cell_similarity_dgp_monomial_06.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_06.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_07.cc b/tests/fe/cell_similarity_dgp_monomial_07.cc
index 5acc0e4db3..09729810f8 100644
--- a/tests/fe/cell_similarity_dgp_monomial_07.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_07.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_08.cc b/tests/fe/cell_similarity_dgp_monomial_08.cc
index 36f4b7e3d5..d722d680e1 100644
--- a/tests/fe/cell_similarity_dgp_monomial_08.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_08.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_09.cc b/tests/fe/cell_similarity_dgp_monomial_09.cc
index bd8724a47f..84678074f0 100644
--- a/tests/fe/cell_similarity_dgp_monomial_09.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_09.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_monomial_10.cc b/tests/fe/cell_similarity_dgp_monomial_10.cc
index 9b06eca60c..2d9df34f06 100644
--- a/tests/fe/cell_similarity_dgp_monomial_10.cc
+++ b/tests/fe/cell_similarity_dgp_monomial_10.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_01.cc b/tests/fe/cell_similarity_dgp_nonparametric_01.cc
index fe4536309c..9fd1b1fa29 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_01.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_01.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_02.cc b/tests/fe/cell_similarity_dgp_nonparametric_02.cc
index 098d36c2ff..3bdba170d6 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_02.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_02.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_03.cc b/tests/fe/cell_similarity_dgp_nonparametric_03.cc
index a54bac82ba..cefbbc02b6 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_03.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_03.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_04.cc b/tests/fe/cell_similarity_dgp_nonparametric_04.cc
index 8cd206e03f..3b3e9f5314 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_04.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_04.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_05.cc b/tests/fe/cell_similarity_dgp_nonparametric_05.cc
index 01b84e7d51..31d79e2d19 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_05.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_05.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_06.cc b/tests/fe/cell_similarity_dgp_nonparametric_06.cc
index f4e73f97db..6951cb5ed5 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_06.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_06.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_07.cc b/tests/fe/cell_similarity_dgp_nonparametric_07.cc
index 81becf6d93..b7279e242a 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_07.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_07.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_08.cc b/tests/fe/cell_similarity_dgp_nonparametric_08.cc
index cfb9ccc434..b6ace3c40d 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_08.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_08.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_09.cc b/tests/fe/cell_similarity_dgp_nonparametric_09.cc
index a38b42f408..8a5e8dc404 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_09.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_09.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/cell_similarity_dgp_nonparametric_10.cc b/tests/fe/cell_similarity_dgp_nonparametric_10.cc
index e4d71656e2..683ae7f92d 100644
--- a/tests/fe/cell_similarity_dgp_nonparametric_10.cc
+++ b/tests/fe/cell_similarity_dgp_nonparametric_10.cc
@@ -98,7 +98,7 @@ test(const Triangulation<dim> &tr)
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
         for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               mass_matrix[cell->index()](i, j) += fe_values.shape_value(i, q) *
                                                   fe_values.shape_value(j, q) *
diff --git a/tests/fe/curl_curl_01.cc b/tests/fe/curl_curl_01.cc
index 6ae2aee7e6..95568c013d 100644
--- a/tests/fe/curl_curl_01.cc
+++ b/tests/fe/curl_curl_01.cc
@@ -348,7 +348,7 @@ MaxwellProblem<dim>::assemble_system()
       fe_values.reinit(cell);
       right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
                                         rhs_values);
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         {
           for (unsigned int i = 0; i < dofs_per_cell; ++i)
             {
diff --git a/tests/fe/deformed_projection.h b/tests/fe/deformed_projection.h
index 81f67dfb35..f957134a5d 100644
--- a/tests/fe/deformed_projection.h
+++ b/tests/fe/deformed_projection.h
@@ -151,7 +151,7 @@ void EvaluateDerivative(DoFHandler<2> *dof_handler, Vector<double> &solution)
         n_q_points, std::vector<Tensor<1, 2>>(n_components));
       fe_values.get_function_gradients(solution, grads_here);
 
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         {
           //   double u0 = this_value[q_point](0);
           // double v0 = this_value[q_point](1);
@@ -243,7 +243,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             {
               coefficient->value_list(fe_values.get_quadrature_points(),
                                       coefficient_values);
-              for (unsigned int point = 0; point < n_q_points; ++point)
+              for (auto point : fe_values.quadrature_point_indices())
                 {
                   const double weight = fe_values.JxW(point);
                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -266,7 +266,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
             {
               coefficient->vector_value_list(fe_values.get_quadrature_points(),
                                              coefficient_vector_values);
-              for (unsigned int point = 0; point < n_q_points; ++point)
+              for (auto point : fe_values.quadrature_point_indices())
                 {
                   const double weight = fe_values.JxW(point);
                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -316,7 +316,7 @@ create_mass_matrix(const Mapping<dim> &       mapping,
                 }
             }
 
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (auto point : fe_values.quadrature_point_indices())
             {
               const double weight = fe_values.JxW(point);
               //     const double weight = q.weight(point);
@@ -425,7 +425,7 @@ create_right_hand_side(const Mapping<dim> &   mapping,
                                   rhs_values);
 
           cell_vector = 0;
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (auto point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               cell_vector(i) += rhs_values[point] *
                                 fe_values.shape_value(i, point) *
@@ -451,7 +451,7 @@ create_right_hand_side(const Mapping<dim> &   mapping,
                                          rhs_values);
 
           cell_vector = 0;
-          for (unsigned int point = 0; point < n_q_points; ++point)
+          for (auto point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               for (unsigned int comp_i = 0; comp_i < fe.n_components();
                    ++comp_i)
diff --git a/tests/fe/dgp_01.cc b/tests/fe/dgp_01.cc
index 896990794a..6017e4da4f 100644
--- a/tests/fe/dgp_01.cc
+++ b/tests/fe/dgp_01.cc
@@ -57,7 +57,7 @@ test(const unsigned int degree)
        ++cell)
     {
       fe_values.reinit(cell);
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         AssertThrow(fe_values.shape_value(0, q) == 1, ExcInternalError());
     }
 
diff --git a/tests/fe/fe_abf_gradient_divergence_theorem.cc b/tests/fe/fe_abf_gradient_divergence_theorem.cc
index ad5ebc2d0a..1d7c7ae1f8 100644
--- a/tests/fe/fe_abf_gradient_divergence_theorem.cc
+++ b/tests/fe/fe_abf_gradient_divergence_theorem.cc
@@ -104,7 +104,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
@@ -114,9 +114,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       boundary_integral +=
                         fe_face_values[single_component].value(i, q) *
diff --git a/tests/fe/fe_br.cc b/tests/fe/fe_br.cc
index 45636bf8d7..2e55487363 100644
--- a/tests/fe/fe_br.cc
+++ b/tests/fe/fe_br.cc
@@ -93,7 +93,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
diff --git a/tests/fe/fe_dgp_3rd_derivative_divergence_theorem.cc b/tests/fe/fe_dgp_3rd_derivative_divergence_theorem.cc
index 4add01e1d2..d7b6cb6afb 100644
--- a/tests/fe/fe_dgp_3rd_derivative_divergence_theorem.cc
+++ b/tests/fe/fe_dgp_3rd_derivative_divergence_theorem.cc
@@ -102,7 +102,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<3, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral +=
                     fe_values[single_component].third_derivative(i, q) *
@@ -116,9 +116,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<2, dim> hessian =
                         fe_face_values[single_component].hessian(i, q);
diff --git a/tests/fe/fe_dgq_3rd_derivative_divergence_theorem.cc b/tests/fe/fe_dgq_3rd_derivative_divergence_theorem.cc
index 60719e0c5a..c2fed59e2b 100644
--- a/tests/fe/fe_dgq_3rd_derivative_divergence_theorem.cc
+++ b/tests/fe/fe_dgq_3rd_derivative_divergence_theorem.cc
@@ -101,7 +101,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<3, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral +=
                     fe_values[single_component].third_derivative(i, q) *
@@ -115,9 +115,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<2, dim> hessian =
                         fe_face_values[single_component].hessian(i, q);
diff --git a/tests/fe/fe_dgq_gradient_divergence_theorem.cc b/tests/fe/fe_dgq_gradient_divergence_theorem.cc
index d1b4eacc53..d988e963d2 100644
--- a/tests/fe/fe_dgq_gradient_divergence_theorem.cc
+++ b/tests/fe/fe_dgq_gradient_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       boundary_integral +=
                         fe_face_values[single_component].value(i, q) *
diff --git a/tests/fe/fe_dgq_hessian_divergence_theorem.cc b/tests/fe/fe_dgq_hessian_divergence_theorem.cc
index f8b2e29add..292043c293 100644
--- a/tests/fe/fe_dgq_hessian_divergence_theorem.cc
+++ b/tests/fe/fe_dgq_hessian_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<2, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].hessian(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<1, dim> gradient =
                         fe_face_values[single_component].gradient(i, q);
diff --git a/tests/fe/fe_enriched_02.cc b/tests/fe/fe_enriched_02.cc
index 5a1649369d..90e5178aac 100644
--- a/tests/fe/fe_enriched_02.cc
+++ b/tests/fe/fe_enriched_02.cc
@@ -111,7 +111,7 @@ test1()
         fe_values.get_quadrature_points();
 
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           deallog << "dof=" << i << " qp=" << q_points[q_point]
                   << " f(qp)=" << function.value(q_points[q_point])
                   << " N(qp)=" << fe_values.shape_value(i, q_point)
diff --git a/tests/fe/fe_enriched_03.cc b/tests/fe/fe_enriched_03.cc
index 4b7a3249b4..b37de778f6 100644
--- a/tests/fe/fe_enriched_03.cc
+++ b/tests/fe/fe_enriched_03.cc
@@ -105,7 +105,7 @@ test2()
           fe_face_values.get_quadrature_points();
 
         for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (const auto q_point : fe_face_values.quadrature_point_indices())
             deallog << "dof=" << i << " qp=" << q_points[q_point]
                     << " f(qp)=" << function.value(q_points[q_point])
                     << " N(qp)=" << fe_face_values.shape_value(i, q_point)
diff --git a/tests/fe/fe_enriched_04.cc b/tests/fe/fe_enriched_04.cc
index b83b6e8124..53624c2b4f 100644
--- a/tests/fe/fe_enriched_04.cc
+++ b/tests/fe/fe_enriched_04.cc
@@ -130,7 +130,7 @@ test3()
       fe_values.get_function_values(solution_pou, solution_values_pou);
       fe_values.get_function_values(solution, solution_values);
 
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         deallog << " qp=" << q_points[q_point]
                 << " f(qp)=" << function.value(q_points[q_point])
                 << " U_fe(qp)=" << solution_values_fe[q_point]
diff --git a/tests/fe/fe_enriched_06.cc b/tests/fe/fe_enriched_06.cc
index d75e79a890..e635250d0e 100644
--- a/tests/fe/fe_enriched_06.cc
+++ b/tests/fe/fe_enriched_06.cc
@@ -144,7 +144,7 @@ test5()
         fe_values.get_function_values(solution, solution_values);
 
         deallog << " cell=" << cell->center() << std::endl;
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           {
             // find non-zero shape_value
             deallog << " qp=" << q_points[q_point]
diff --git a/tests/fe/fe_enriched_color_07.cc b/tests/fe/fe_enriched_color_07.cc
index 89e4f4f55c..dfdba4fbba 100644
--- a/tests/fe/fe_enriched_color_07.cc
+++ b/tests/fe/fe_enriched_color_07.cc
@@ -752,7 +752,7 @@ EstimateEnrichmentFunction::assemble_system()
       rhs_values.resize(n_q_points);
       rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
 
-      for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
+      for (const auto q_index : fe_values.quadrature_point_indices())
         {
           double radius = center.distance(fe_values.quadrature_point(q_index));
 
@@ -1521,7 +1521,7 @@ LaplaceProblem<dim>::assemble_system()
                                   tmp_rhs_value);
 
             // add tmp to the total one at quadrature points
-            for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            for (const auto q_point : fe_values.quadrature_point_indices())
               {
                 rhs_value[q_point] += tmp_rhs_value[q_point];
               }
@@ -1534,7 +1534,7 @@ LaplaceProblem<dim>::assemble_system()
         cell_system_matrix = 0;
         cell_rhs           = 0;
 
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           for (unsigned int i = 0; i < dofs_per_cell; ++i)
             {
               for (unsigned int j = i; j < dofs_per_cell; ++j)
diff --git a/tests/fe/fe_enriched_compare_to_fe_system.cc b/tests/fe/fe_enriched_compare_to_fe_system.cc
index e03ad8e44a..f2f1ca9dd7 100644
--- a/tests/fe/fe_enriched_compare_to_fe_system.cc
+++ b/tests/fe/fe_enriched_compare_to_fe_system.cc
@@ -229,7 +229,7 @@ test(const FiniteElement<dim> & fe1,
 
       // check shape functions
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values_system.quadrature_point_indices())
           check_consistency(
             q_points[q_point],
             function,
@@ -252,7 +252,8 @@ test(const FiniteElement<dim> & fe1,
             fe_face_values_system.get_quadrature_points();
 
           for (unsigned int i = 0; i < dofs_per_cell; ++i)
-            for (unsigned int q_point = 0; q_point < n_q_points_face; ++q_point)
+            for (const auto q_point :
+                 fe_face_values_system.quadrature_point_indices())
               check_consistency(
                 q_points[q_point],
                 function,
diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.cc b/tests/fe/fe_enriched_compare_to_fe_system_2.cc
index ab3550f94e..21dc8784a5 100644
--- a/tests/fe/fe_enriched_compare_to_fe_system_2.cc
+++ b/tests/fe/fe_enriched_compare_to_fe_system_2.cc
@@ -308,7 +308,7 @@ test(const FiniteElement<dim> & fe_base,
 
       // check shape functions
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values_system.quadrature_point_indices())
           check_consistency(
             q_points[q_point],
             fun1,
@@ -339,7 +339,8 @@ test(const FiniteElement<dim> & fe_base,
             fe_face_values_system.get_quadrature_points();
 
           for (unsigned int i = 0; i < dofs_per_cell; ++i)
-            for (unsigned int q_point = 0; q_point < n_q_points_face; ++q_point)
+            for (const auto q_point :
+                 fe_face_values_system.quadrature_point_indices())
               check_consistency(
                 q_points[q_point],
                 fun1,
diff --git a/tests/fe/fe_enriched_step-36.cc b/tests/fe/fe_enriched_step-36.cc
index 82cd66b790..0add928260 100644
--- a/tests/fe/fe_enriched_step-36.cc
+++ b/tests/fe/fe_enriched_step-36.cc
@@ -410,7 +410,7 @@ namespace Step36
           potential.value_list(fe_values.get_quadrature_points(),
                                potential_values);
 
-          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (const auto q_point : fe_values.quadrature_point_indices())
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               for (unsigned int j = i; j < dofs_per_cell; ++j)
                 {
diff --git a/tests/fe/fe_enriched_step-36b.cc b/tests/fe/fe_enriched_step-36b.cc
index 8a6187bf81..38c79aa84e 100644
--- a/tests/fe/fe_enriched_step-36b.cc
+++ b/tests/fe/fe_enriched_step-36b.cc
@@ -480,7 +480,7 @@ namespace Step36
 
           if (cell->active_fe_index() == 0) // plain FE
             {
-              for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+              for (const auto q_point : fe_values.quadrature_point_indices())
                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
                   for (unsigned int j = i; j < dofs_per_cell; ++j)
                     {
@@ -518,8 +518,8 @@ namespace Step36
 
                       if (i_group == fe_group && j_group == fe_group) // fe - fe
                         {
-                          for (unsigned int q_point = 0; q_point < n_q_points;
-                               ++q_point)
+                          for (const auto q_point :
+                               fe_values.quadrature_point_indices())
                             {
                               cell_stiffness_matrix(i, j) +=
                                 (fe_values[fe_extractor].gradient(i, q_point) *
@@ -540,8 +540,8 @@ namespace Step36
                       else if (i_group == fe_group &&
                                j_group == pou_group) // fe - pou
                         {
-                          for (unsigned int q_point = 0; q_point < n_q_points;
-                               ++q_point)
+                          for (const auto q_point :
+                               fe_values.quadrature_point_indices())
                             {
                               cell_stiffness_matrix(i, j) +=
                                 (fe_values[fe_extractor].gradient(i, q_point) *
@@ -567,8 +567,8 @@ namespace Step36
                       else if (i_group == pou_group &&
                                j_group == fe_group) // pou - fe
                         {
-                          for (unsigned int q_point = 0; q_point < n_q_points;
-                               ++q_point)
+                          for (const auto q_point :
+                               fe_values.quadrature_point_indices())
                             {
                               cell_stiffness_matrix(i, j) +=
                                 ((fe_values[pou_extractor].gradient(i,
@@ -597,8 +597,8 @@ namespace Step36
                           Assert(i_group == pou_group && j_group == pou_group,
                                  ExcInternalError());
 
-                          for (unsigned int q_point = 0; q_point < n_q_points;
-                               ++q_point)
+                          for (const auto q_point :
+                               fe_values.quadrature_point_indices())
                             {
                               cell_stiffness_matrix(i, j) +=
                                 ((fe_values[pou_extractor].gradient(i,
diff --git a/tests/fe/fe_face_orientation_nedelec.h b/tests/fe/fe_face_orientation_nedelec.h
index 8d8c0c4fb1..b0ac7a98f4 100644
--- a/tests/fe/fe_face_orientation_nedelec.h
+++ b/tests/fe/fe_face_orientation_nedelec.h
@@ -212,7 +212,7 @@ evaluate(const FE_Nedelec<3> & fe,
       fe_field_function.vector_value_list(fe_values.get_quadrature_points(),
                                           values);
 
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < 3; ++d)
             deallog << values_ref[q_point][d] - values[q_point](d) << "  ";
diff --git a/tests/fe/fe_face_orientation_nedelec_0.cc b/tests/fe/fe_face_orientation_nedelec_0.cc
index c5e6deaeee..f47f8d7de3 100644
--- a/tests/fe/fe_face_orientation_nedelec_0.cc
+++ b/tests/fe/fe_face_orientation_nedelec_0.cc
@@ -129,7 +129,7 @@ evaluate(const FE_Nedelec<3> & fe,
                                           values);
       std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
       cell->get_dof_indices(dof_indices);
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (const auto q_point : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < 3; ++d)
             deallog << values_ref[q_point][d] - values[q_point](d) << "  ";
diff --git a/tests/fe/fe_nedelec_gradient_divergence_theorem.cc b/tests/fe/fe_nedelec_gradient_divergence_theorem.cc
index d8818c4705..f5091d50c4 100644
--- a/tests/fe/fe_nedelec_gradient_divergence_theorem.cc
+++ b/tests/fe/fe_nedelec_gradient_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       boundary_integral +=
                         fe_face_values[single_component].value(i, q) *
diff --git a/tests/fe/fe_nedelec_hessian_divergence_theorem.cc b/tests/fe/fe_nedelec_hessian_divergence_theorem.cc
index 629f04d0e0..f669ee5310 100644
--- a/tests/fe/fe_nedelec_hessian_divergence_theorem.cc
+++ b/tests/fe/fe_nedelec_hessian_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<2, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].hessian(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<1, dim> gradient =
                         fe_face_values[single_component].gradient(i, q);
diff --git a/tests/fe/fe_nedelec_singularity_01.cc b/tests/fe/fe_nedelec_singularity_01.cc
index 52f3b4428f..2456920458 100644
--- a/tests/fe/fe_nedelec_singularity_01.cc
+++ b/tests/fe/fe_nedelec_singularity_01.cc
@@ -448,7 +448,7 @@ namespace nedelec_singularity
             right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
                                               rhs_values);
 
-            for (unsigned int q = 0; q < n_q_points; ++q)
+            for (const auto q : fe_values.quadrature_point_indices())
               {
                 Tensor<1, dim, std::complex<double>> rhs;
 
diff --git a/tests/fe/fe_nedelec_singularity_02.cc b/tests/fe/fe_nedelec_singularity_02.cc
index 3987aebf31..107a153f2e 100644
--- a/tests/fe/fe_nedelec_singularity_02.cc
+++ b/tests/fe/fe_nedelec_singularity_02.cc
@@ -437,7 +437,7 @@ namespace nedelec_singularity
             right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
                                               rhs_values);
 
-            for (unsigned int q = 0; q < n_q_points; ++q)
+            for (const auto q : fe_values.quadrature_point_indices())
               {
                 Tensor<1, dim> rhs;
 
diff --git a/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc
index ffc09e9ec9..55ed223edc 100644
--- a/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc
+++ b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc
@@ -102,7 +102,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
@@ -112,9 +112,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       boundary_integral +=
                         fe_face_values[single_component].value(i, q) *
diff --git a/tests/fe/fe_nedelec_sz_non_rect_2d.cc b/tests/fe/fe_nedelec_sz_non_rect_2d.cc
index b06269d40d..2e1e97010e 100644
--- a/tests/fe/fe_nedelec_sz_non_rect_2d.cc
+++ b/tests/fe/fe_nedelec_sz_non_rect_2d.cc
@@ -252,7 +252,7 @@ namespace polytest
 
         right_hand_side.rhs_value_list(fe_values.get_quadrature_points(),
                                        rhs_value_list);
-        for (unsigned int q = 0; q < n_q_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             Tensor<1, dim> rhs_value;
             for (unsigned int d = 0; d < dim; ++d)
diff --git a/tests/fe/fe_nedelec_sz_non_rect_face.cc b/tests/fe/fe_nedelec_sz_non_rect_face.cc
index 2546263d0f..13d9f25431 100644
--- a/tests/fe/fe_nedelec_sz_non_rect_face.cc
+++ b/tests/fe/fe_nedelec_sz_non_rect_face.cc
@@ -254,7 +254,7 @@ namespace Maxwell
         // Calc values of curlE from fe solution:
         cell->get_dof_indices(local_dof_indices);
         // Loop over quad points to calculate solution:
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           {
             // Split exact solution into real/imaginary parts:
             for (unsigned int component = 0; component < dim; component++)
@@ -371,7 +371,7 @@ namespace Maxwell
                                          rhs_value_list);
 
         // Loop over all element quad points:
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           {
             // store rhs value at this q point & turn into tensor
             for (unsigned int component = 0; component < dim; component++)
diff --git a/tests/fe/fe_project_2d.cc b/tests/fe/fe_project_2d.cc
index c52d5828ef..c0ebd4c984 100644
--- a/tests/fe/fe_project_2d.cc
+++ b/tests/fe/fe_project_2d.cc
@@ -228,7 +228,7 @@ test(const FiniteElement<dim> &fe,
           const std::vector<double> &JxW_values = fe_values.get_JxW_values();
           fe_values[vec].get_function_divergences(v, div_v);
           fe_values[vec].get_function_curls(v, curl_v);
-          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (const auto q_point : fe_values.quadrature_point_indices())
             {
               total_div += JxW_values[q_point] * div_v[q_point];
               total_curl += JxW_values[q_point] * curl_v[q_point];
diff --git a/tests/fe/fe_project_3d.cc b/tests/fe/fe_project_3d.cc
index fb94bb407a..9c64fd11e6 100644
--- a/tests/fe/fe_project_3d.cc
+++ b/tests/fe/fe_project_3d.cc
@@ -311,7 +311,7 @@ test(const FiniteElement<dim> &fe,
           fe_values[vec].get_function_divergences(v, div_v);
           fe_values[vec].get_function_curls(v, curl_v);
           fe_values[vec].get_function_hessians(v, hessians);
-          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (const auto q_point : fe_values.quadrature_point_indices())
             {
               total_div += JxW_values[q_point] * div_v[q_point];
               total_curl += JxW_values[q_point] * curl_v[q_point];
diff --git a/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc b/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc
index a6e826c6c0..5b40e594e9 100644
--- a/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc
+++ b/tests/fe/fe_q_3rd_derivative_divergence_theorem.cc
@@ -101,7 +101,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<3, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral +=
                     fe_values[single_component].third_derivative(i, q) *
@@ -115,9 +115,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<2, dim> hessian =
                         fe_face_values[single_component].hessian(i, q);
diff --git a/tests/fe/fe_q_bubbles.cc b/tests/fe/fe_q_bubbles.cc
index 2cb779ee82..aa0df74496 100644
--- a/tests/fe/fe_q_bubbles.cc
+++ b/tests/fe/fe_q_bubbles.cc
@@ -222,7 +222,7 @@ Step3<dim>::assemble_system(unsigned int i)
       cell_rhs    = 0;
 
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           {
             for (unsigned int j = 0; j < dofs_per_cell; ++j)
               cell_matrix(i, j) +=
diff --git a/tests/fe/fe_q_dg0.cc b/tests/fe/fe_q_dg0.cc
index 618224d198..b216d956bf 100644
--- a/tests/fe/fe_q_dg0.cc
+++ b/tests/fe/fe_q_dg0.cc
@@ -556,7 +556,7 @@ namespace Step22
 
         right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
                                           rhs_values);
-        for (unsigned int q = 0; q < n_q_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             for (unsigned int k = 0; k < dofs_per_cell; ++k)
               {
@@ -604,7 +604,7 @@ namespace Step22
                     const std::vector<Point<dim>> &quad_points =
                       fe_v_face.get_quadrature_points();
 
-                    for (unsigned int q = 0; q < n_q_face; ++q)
+                    for (const auto q : fe_v_face.quadrature_point_indices())
                       {
                         double jump =
                           jumpfunction.jump(quad_points[q], normals[q]);
@@ -932,7 +932,7 @@ namespace Step22
         fe_v.reinit(cell);
         cell->get_dof_indices(local_dof_indices);
 
-        for (unsigned int q = 0; q < n_q_points; ++q)
+        for (const auto q : fe_v.quadrature_point_indices())
           {
             double div = 0;
             for (unsigned int i = 0; i < dofs_per_cell - 1; ++i)
diff --git a/tests/fe/fe_q_gradient_divergence_theorem.cc b/tests/fe/fe_q_gradient_divergence_theorem.cc
index 85dd24de4b..3f221ee4c7 100644
--- a/tests/fe/fe_q_gradient_divergence_theorem.cc
+++ b/tests/fe/fe_q_gradient_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       boundary_integral +=
                         fe_face_values[single_component].value(i, q) *
diff --git a/tests/fe/fe_q_hessian_divergence_theorem.cc b/tests/fe/fe_q_hessian_divergence_theorem.cc
index 42f83c22a5..fa016d9b0a 100644
--- a/tests/fe/fe_q_hessian_divergence_theorem.cc
+++ b/tests/fe/fe_q_hessian_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<2, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].hessian(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<1, dim> gradient =
                         fe_face_values[single_component].gradient(i, q);
diff --git a/tests/fe/fe_q_hierarchical_3rd_derivative_divergence_theorem.cc b/tests/fe/fe_q_hierarchical_3rd_derivative_divergence_theorem.cc
index 5b60be2f45..0ff447428f 100644
--- a/tests/fe/fe_q_hierarchical_3rd_derivative_divergence_theorem.cc
+++ b/tests/fe/fe_q_hierarchical_3rd_derivative_divergence_theorem.cc
@@ -102,7 +102,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<3, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral +=
                     fe_values[single_component].third_derivative(i, q) *
@@ -116,9 +116,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<2, dim> hessian =
                         fe_face_values[single_component].hessian(i, q);
diff --git a/tests/fe/fe_rt_gradient_divergence_theorem.cc b/tests/fe/fe_rt_gradient_divergence_theorem.cc
index f9e3615777..c9931a78ef 100644
--- a/tests/fe/fe_rt_gradient_divergence_theorem.cc
+++ b/tests/fe/fe_rt_gradient_divergence_theorem.cc
@@ -103,7 +103,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<1, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].gradient(i, q) *
                                    fe_values.JxW(q);
@@ -113,9 +113,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       boundary_integral +=
                         fe_face_values[single_component].value(i, q) *
diff --git a/tests/fe/fe_rt_hessian_divergence_theorem.cc b/tests/fe/fe_rt_hessian_divergence_theorem.cc
index 7e1b90f2b2..59370e7993 100644
--- a/tests/fe/fe_rt_hessian_divergence_theorem.cc
+++ b/tests/fe/fe_rt_hessian_divergence_theorem.cc
@@ -97,7 +97,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<2, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral += fe_values[single_component].hessian(i, q) *
                                    fe_values.JxW(q);
@@ -107,9 +107,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<1, dim> gradient =
                         fe_face_values[single_component].gradient(i, q);
diff --git a/tests/fe/fe_system_3rd_derivative_divergence_theorem.cc b/tests/fe/fe_system_3rd_derivative_divergence_theorem.cc
index 7b08f3b088..6ba744632a 100644
--- a/tests/fe/fe_system_3rd_derivative_divergence_theorem.cc
+++ b/tests/fe/fe_system_3rd_derivative_divergence_theorem.cc
@@ -101,7 +101,7 @@ test(const Triangulation<dim> &tr,
               ss << "component=" << c << ", dof=" << i << std::endl;
 
               Tensor<3, dim> bulk_integral;
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   bulk_integral +=
                     fe_values[single_component].third_derivative(i, q) *
@@ -115,9 +115,7 @@ test(const Triangulation<dim> &tr,
               for (const unsigned int face : GeometryInfo<dim>::face_indices())
                 {
                   fe_face_values.reinit(cell, face);
-                  for (unsigned int q = 0;
-                       q < fe_face_values.n_quadrature_points;
-                       ++q)
+                  for (const auto q : fe_face_values.quadrature_point_indices())
                     {
                       Tensor<2, dim> hessian =
                         fe_face_values[single_component].hessian(i, q);
diff --git a/tests/fe/fe_values_function_manifold.cc b/tests/fe/fe_values_function_manifold.cc
index 2bdad7b9e8..58b706a69d 100644
--- a/tests/fe/fe_values_function_manifold.cc
+++ b/tests/fe/fe_values_function_manifold.cc
@@ -320,9 +320,7 @@ JxWError<dim>::setup_matrices()
       cell_rhs    = 0.0;
       fe_values.reinit(cell);
 
-      for (unsigned int q_point_n = 0;
-           q_point_n < fe_values.n_quadrature_points;
-           ++q_point_n)
+      for (const auto q_point_n : fe_values.quadrature_point_indices())
         {
           const double point_forcing =
             manufactured_forcing->value(fe_values.quadrature_point(q_point_n));
diff --git a/tests/fe/fe_values_view_01.cc b/tests/fe/fe_values_view_01.cc
index 9c91670ba6..7682ae65c3 100644
--- a/tests/fe/fe_values_view_01.cc
+++ b/tests/fe/fe_values_view_01.cc
@@ -58,7 +58,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
       FEValuesExtractors::Scalar single_component(c);
 
       for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             deallog << "i=" << i << ", q=" << q << std::endl;
             deallog << "   " << fe_values[single_component].value(i, q) << ' ';
diff --git a/tests/fe/fe_values_view_02.cc b/tests/fe/fe_values_view_02.cc
index 3876bec6bd..dac397d4c2 100644
--- a/tests/fe/fe_values_view_02.cc
+++ b/tests/fe/fe_values_view_02.cc
@@ -64,7 +64,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
         FEValuesExtractors::Vector vec_components(c);
 
         for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               deallog << "i=" << i << ", q=" << q << std::endl;
               deallog << "   " << fe_values[vec_components].value(i, q) << ' ';
diff --git a/tests/fe/fe_values_view_03.cc b/tests/fe/fe_values_view_03.cc
index c9fddb271f..1db642e4a9 100644
--- a/tests/fe/fe_values_view_03.cc
+++ b/tests/fe/fe_values_view_03.cc
@@ -64,7 +64,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           FEValuesExtractors::Scalar single_component(c);
 
           for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
-            for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+            for (const auto q : fe_values.quadrature_point_indices())
               {
                 deallog << "i=" << i << ", q=" << q << std::endl;
                 deallog << "   " << fe_values[single_component].value(i, q)
diff --git a/tests/fe/fe_values_view_04.cc b/tests/fe/fe_values_view_04.cc
index 71349450cb..bbfee129cf 100644
--- a/tests/fe/fe_values_view_04.cc
+++ b/tests/fe/fe_values_view_04.cc
@@ -69,7 +69,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
             FEValuesExtractors::Vector vec_components(c);
 
             for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i)
-              for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   deallog << "i=" << i << ", q=" << q << std::endl;
                   deallog << "   ";
diff --git a/tests/fe/fe_values_view_05.cc b/tests/fe/fe_values_view_05.cc
index ea7e346438..aae9fcd18b 100644
--- a/tests/fe/fe_values_view_05.cc
+++ b/tests/fe/fe_values_view_05.cc
@@ -70,7 +70,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                       scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << scalar_values[q] << std::endl;
           Assert(std::fabs(scalar_values[q] - vector_values[q](c)) <=
diff --git a/tests/fe/fe_values_view_06.cc b/tests/fe/fe_values_view_06.cc
index 7a94c8ed66..8830080b62 100644
--- a/tests/fe/fe_values_view_06.cc
+++ b/tests/fe/fe_values_view_06.cc
@@ -70,7 +70,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                          scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << scalar_values[q] << std::endl;
           Assert((scalar_values[q] - vector_values[q][c]).norm() <=
diff --git a/tests/fe/fe_values_view_07.cc b/tests/fe/fe_values_view_07.cc
index 3edbbb0b2c..d49bfd85c2 100644
--- a/tests/fe/fe_values_view_07.cc
+++ b/tests/fe/fe_values_view_07.cc
@@ -70,7 +70,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                         scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << scalar_values[q] << std::endl;
           Assert((scalar_values[q] - vector_values[q][c]).norm() <=
diff --git a/tests/fe/fe_values_view_08.cc b/tests/fe/fe_values_view_08.cc
index c6672347a5..bfe9fe2d5f 100644
--- a/tests/fe/fe_values_view_08.cc
+++ b/tests/fe/fe_values_view_08.cc
@@ -72,7 +72,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                       scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           deallog << scalar_values[q] << std::endl;
           Assert(std::fabs(scalar_values[q] - vector_values[q](c)) <=
diff --git a/tests/fe/fe_values_view_09.cc b/tests/fe/fe_values_view_09.cc
index fba93b0b80..272c05792b 100644
--- a/tests/fe/fe_values_view_09.cc
+++ b/tests/fe/fe_values_view_09.cc
@@ -72,7 +72,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                          scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < dim; ++d)
             deallog << scalar_values[q][d] << (d < dim - 1 ? " " : "");
diff --git a/tests/fe/fe_values_view_10.cc b/tests/fe/fe_values_view_10.cc
index d458788d2f..202e4b3bd0 100644
--- a/tests/fe/fe_values_view_10.cc
+++ b/tests/fe/fe_values_view_10.cc
@@ -72,7 +72,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                         scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_10_single_01.cc b/tests/fe/fe_values_view_10_single_01.cc
index 3f62b6aca0..6223c54825 100644
--- a/tests/fe/fe_values_view_10_single_01.cc
+++ b/tests/fe/fe_values_view_10_single_01.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                         scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_10_single_02.cc b/tests/fe/fe_values_view_10_single_02.cc
index d12c5388b9..16eef00a8a 100644
--- a/tests/fe/fe_values_view_10_single_02.cc
+++ b/tests/fe/fe_values_view_10_single_02.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                         scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_10_single_03.cc b/tests/fe/fe_values_view_10_single_03.cc
index 4d211b6f10..33f582cfe9 100644
--- a/tests/fe/fe_values_view_10_single_03.cc
+++ b/tests/fe/fe_values_view_10_single_03.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
                                                         scalar_values);
       deallog << "component=" << c << std::endl;
 
-      for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+      for (const auto q : fe_values.quadrature_point_indices())
         {
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_11.cc b/tests/fe/fe_values_view_11.cc
index 5468337e94..9290cebf57 100644
--- a/tests/fe/fe_values_view_11.cc
+++ b/tests/fe/fe_values_view_11.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               deallog << selected_vector_values[q][d] << std::endl;
diff --git a/tests/fe/fe_values_view_12.cc b/tests/fe/fe_values_view_12.cc
index e71e558ec9..8a983d0596 100644
--- a/tests/fe/fe_values_view_12.cc
+++ b/tests/fe/fe_values_view_12.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               deallog << selected_vector_values[q][d] << std::endl;
diff --git a/tests/fe/fe_values_view_13.cc b/tests/fe/fe_values_view_13.cc
index 4f11d8bc15..8766887332 100644
--- a/tests/fe/fe_values_view_13.cc
+++ b/tests/fe/fe_values_view_13.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               deallog << selected_vector_values[q][d] << std::endl;
diff --git a/tests/fe/fe_values_view_14.cc b/tests/fe/fe_values_view_14.cc
index 418d0411a3..573438f541 100644
--- a/tests/fe/fe_values_view_14.cc
+++ b/tests/fe/fe_values_view_14.cc
@@ -77,7 +77,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               deallog << selected_vector_values[q][d] << std::endl;
diff --git a/tests/fe/fe_values_view_15.cc b/tests/fe/fe_values_view_15.cc
index 828608e575..f7e8271019 100644
--- a/tests/fe/fe_values_view_15.cc
+++ b/tests/fe/fe_values_view_15.cc
@@ -77,7 +77,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_16.cc b/tests/fe/fe_values_view_16.cc
index 3e1acfa267..4b16042cf2 100644
--- a/tests/fe/fe_values_view_16.cc
+++ b/tests/fe/fe_values_view_16.cc
@@ -77,7 +77,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_17.cc b/tests/fe/fe_values_view_17.cc
index 4f179c88a2..b834b6c6ed 100644
--- a/tests/fe/fe_values_view_17.cc
+++ b/tests/fe/fe_values_view_17.cc
@@ -77,7 +77,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             deallog << selected_vector_values[q] << std::endl;
 
diff --git a/tests/fe/fe_values_view_18.cc b/tests/fe/fe_values_view_18.cc
index 2ccbc247dd..52ca56c42a 100644
--- a/tests/fe/fe_values_view_18.cc
+++ b/tests/fe/fe_values_view_18.cc
@@ -79,7 +79,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             for (unsigned int d = 0; d < dim; ++d)
               for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/fe_values_view_19.cc b/tests/fe/fe_values_view_19.cc
index d39fb0c8fc..52d0b58207 100644
--- a/tests/fe/fe_values_view_19.cc
+++ b/tests/fe/fe_values_view_19.cc
@@ -75,7 +75,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             deallog << selected_vector_values[q] << std::endl;
 
diff --git a/tests/fe/fe_values_view_20.cc b/tests/fe/fe_values_view_20.cc
index 8f22e12114..ebf84a2014 100644
--- a/tests/fe/fe_values_view_20.cc
+++ b/tests/fe/fe_values_view_20.cc
@@ -77,7 +77,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           fe_function, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             deallog << selected_vector_values[q] << std::endl;
 
diff --git a/tests/fe/fe_values_view_30.cc b/tests/fe/fe_values_view_30.cc
index 3796c427f1..bc993f9efe 100644
--- a/tests/fe/fe_values_view_30.cc
+++ b/tests/fe/fe_values_view_30.cc
@@ -155,7 +155,7 @@ test(const Triangulation<dim> &tr, const FiniteElement<dim> &fe)
           function_vals, selected_vector_values);
         deallog << "component=" << c << std::endl;
 
-        for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           for (unsigned int d = 0; d < dim; ++d)
             {
               for (unsigned int e = 0; e < dim; ++e)
diff --git a/tests/fe/nedelec_non_rect_2d.cc b/tests/fe/nedelec_non_rect_2d.cc
index 61e3d6a53f..8b94c637da 100644
--- a/tests/fe/nedelec_non_rect_2d.cc
+++ b/tests/fe/nedelec_non_rect_2d.cc
@@ -250,7 +250,7 @@ namespace polytest
 
         right_hand_side.rhs_value_list(fe_values.get_quadrature_points(),
                                        rhs_value_list);
-        for (unsigned int q = 0; q < n_q_points; ++q)
+        for (const auto q : fe_values.quadrature_point_indices())
           {
             Tensor<1, dim> rhs_value;
             for (unsigned int d = 0; d < dim; ++d)
diff --git a/tests/fe/nedelec_non_rect_face.cc b/tests/fe/nedelec_non_rect_face.cc
index 87294f7fd5..e90ab77adf 100644
--- a/tests/fe/nedelec_non_rect_face.cc
+++ b/tests/fe/nedelec_non_rect_face.cc
@@ -255,7 +255,7 @@ namespace Maxwell
         // Calc values of curlE from fe solution:
         cell->get_dof_indices(local_dof_indices);
         // Loop over quad points to calculate solution:
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           {
             // Split exact solution into real/imaginary parts:
             for (unsigned int component = 0; component < dim; component++)
@@ -372,7 +372,7 @@ namespace Maxwell
                                          rhs_value_list);
 
         // Loop over all element quad points:
-        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (const auto q_point : fe_values.quadrature_point_indices())
           {
             // store rhs value at this q point & turn into tensor
             for (unsigned int component = 0; component < dim; component++)
diff --git a/tests/fe/non_primitive_1.cc b/tests/fe/non_primitive_1.cc
index 71cfda8a1e..1910e32a51 100644
--- a/tests/fe/non_primitive_1.cc
+++ b/tests/fe/non_primitive_1.cc
@@ -79,7 +79,7 @@ create_stokes_matrix_1(const DoFHandler<dim> &dof_handler,
 
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
         for (unsigned int j = 0; j < dofs_per_cell; ++j)
-          for (unsigned int q = 0; q < n_q_points; ++q)
+          for (const auto q : fe_values.quadrature_point_indices())
             {
               const unsigned int comp_i = fe.system_to_component_index(i).first,
                                  comp_j = fe.system_to_component_index(j).first;
@@ -150,7 +150,7 @@ create_stokes_matrix_2(const DoFHandler<dim> &dof_handler,
         for (unsigned int comp_i = 0; comp_i < fe.n_components(); ++comp_i)
           for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
             for (unsigned int comp_j = 0; comp_j < fe.n_components(); ++comp_j)
-              for (unsigned int q = 0; q < n_q_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   // velocity-velocity coupling?
                   if ((comp_i < dim) && (comp_j < dim))
@@ -227,7 +227,7 @@ create_stokes_matrix_3(const DoFHandler<dim> &dof_handler,
               for (unsigned int comp_j = 0; comp_j < fe.n_components();
                    ++comp_j)
                 if (fe.get_nonzero_components(j)[comp_j] == true)
-                  for (unsigned int q = 0; q < n_q_points; ++q)
+                  for (const auto q : fe_values.quadrature_point_indices())
                     {
                       // velocity-velocity coupling?
                       if ((comp_i < dim) && (comp_j < dim))
diff --git a/tests/fe/non_primitive_2.cc b/tests/fe/non_primitive_2.cc
index 90a6f00229..8f7e5f0d65 100644
--- a/tests/fe/non_primitive_2.cc
+++ b/tests/fe/non_primitive_2.cc
@@ -83,7 +83,7 @@ create_stokes_matrix_2(const DoFHandler<dim> &dof_handler,
         for (unsigned int comp_i = 0; comp_i < fe.n_components(); ++comp_i)
           for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
             for (unsigned int comp_j = 0; comp_j < fe.n_components(); ++comp_j)
-              for (unsigned int q = 0; q < n_q_points; ++q)
+              for (const auto q : fe_values.quadrature_point_indices())
                 {
                   // velocity-velocity coupling?
                   if ((comp_i < dim) && (comp_j < dim))
@@ -160,7 +160,7 @@ create_stokes_matrix_3(const DoFHandler<dim> &dof_handler,
               for (unsigned int comp_j = 0; comp_j < fe.n_components();
                    ++comp_j)
                 if (fe.get_nonzero_components(j)[comp_j] == true)
-                  for (unsigned int q = 0; q < n_q_points; ++q)
+                  for (const auto q : fe_values.quadrature_point_indices())
                     {
                       // velocity-velocity coupling?
                       if ((comp_i < dim) && (comp_j < dim))
diff --git a/tests/fe/rt_approximation_01.cc b/tests/fe/rt_approximation_01.cc
index 465f8fb7e8..2b1d535685 100644
--- a/tests/fe/rt_approximation_01.cc
+++ b/tests/fe/rt_approximation_01.cc
@@ -388,7 +388,7 @@ double TestProjection(Mapping<2> &mapping, DoFHandler<2> *dof_handler)
                                                  Vector<double>(n_components));
           fe_values.get_function_values(solution, this_value);
 
-          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (const auto q_point : fe_values.quadrature_point_indices())
             {
               double   u = this_value[q_point](0), v = this_value[q_point](1);
               Point<2> p = fe_values.quadrature_point(q_point);
diff --git a/tests/fe/rt_normal_02.cc b/tests/fe/rt_normal_02.cc
index 37b6c43a66..dbf4e82f48 100644
--- a/tests/fe/rt_normal_02.cc
+++ b/tests/fe/rt_normal_02.cc
@@ -136,7 +136,7 @@ void EvaluateNormal2(DoFHandler<2> *dof_handler, Vector<double> &solution)
                 n_q_proj, Vector<double>(n_components));
               fe_v_n.get_function_values(solution, this_value_n);
 
-              for (unsigned int q_point = 0; q_point < n_q_face; ++q_point)
+              for (const auto q_point : fe_v_face.quadrature_point_indices())
                 {
                   Tensor<1, 2> vn = fe_v_face.normal_vector(q_point);
                   double       nx = vn[0];
@@ -218,7 +218,7 @@ void EvaluateNormal(DoFHandler<2> *dof_handler, Vector<double> &solution)
                 n_q_face, Vector<double>(n_components));
               fe_v_face_n.get_function_values(solution, this_value_n);
 
-              for (unsigned int q_point = 0; q_point < n_q_face; ++q_point)
+              for (const auto q_point : fe_v_face.quadrature_point_indices())
                 {
                   Tensor<1, 2> vn = fe_v_face.normal_vector(q_point);
                   double       nx = vn[0];