From 2c3f735d68f19f1076852555873d307ec7cb5eee Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 2 Feb 2022 15:52:54 -0700
Subject: [PATCH] Use prefix ++ over postfix.

---
 CeresFE/src/ceres.cc                          | 36 ++++-----
 CeresFE/support_code/config_in.h              | 24 +++---
 CeresFE/support_code/ellipsoid_fit.h          |  6 +-
 Distributed_LDG_Method/LDGPoisson.cc          | 44 +++++------
 .../nonlinear-poro-viscoelasticity.cc         |  6 +-
 .../elastoplastic.cc                          | 78 +++++++++----------
 two_phase_flow/LevelSetSolver.cc              | 40 +++++-----
 7 files changed, 117 insertions(+), 117 deletions(-)

diff --git a/CeresFE/src/ceres.cc b/CeresFE/src/ceres.cc
index 9a9045d..88f020a 100644
--- a/CeresFE/src/ceres.cc
+++ b/CeresFE/src/ceres.cc
@@ -484,7 +484,7 @@ namespace Step22
   {
     vector<double> etas;
 
-    for (unsigned int i=0; i < system_parameters::sizeof_depths_eta; i++)
+    for (unsigned int i=0; i < system_parameters::sizeof_depths_eta; ++i)
       {
         etas.push_back(system_parameters::depths_eta[i]);
         etas.push_back(system_parameters::eta_kinks[i]);
@@ -551,7 +551,7 @@ namespace Step22
 
     //find the correct interval to do the interpolation in
     int n_minus_one = -1;
-    for (unsigned int n = 1; n <= n_visc_kinks; n++)
+    for (unsigned int n = 1; n <= n_visc_kinks; ++n)
       {
         unsigned int ndeep = 2 * n - 2;
         unsigned int nshallow = 2 * n;
@@ -693,7 +693,7 @@ namespace Step22
     grav_parameters.push_back(system_parameters::rho[1]);
 
     std::cout << "Body parameters are: " ;
-    for (int i=0; i<6; i++)
+    for (int i=0; i<6; ++i)
       std::cout << grav_parameters[i] << " ";
     std::cout << endl;
 
@@ -1201,7 +1201,7 @@ namespace Step22
             cell_area += fe_values.JxW(q);
             velocities[q] *= fe_values.JxW(q);
             current_cell_velocity += velocities[q];
-            for (unsigned int i = 0; i < (dim+1); i++)
+            for (unsigned int i = 0; i < (dim+1); ++i)
               {
                 velocity_grads[q][i] *= fe_values.JxW(q);
                 current_cell_grads[i] += velocity_grads[q][i];
@@ -1210,7 +1210,7 @@ namespace Step22
             current_cell_old_phiphi_stress += local_quadrature_points_history[q].old_phiphi_stress * fe_values.JxW(q);
           }
         current_cell_velocity /= cell_area;
-        for (unsigned int i = 0; i < (dim+1); i++)
+        for (unsigned int i = 0; i < (dim+1); ++i)
           current_cell_grads[i] /= cell_area;
         current_cell_old_stress /= cell_area;
         current_cell_old_phiphi_stress /= cell_area;
@@ -1232,7 +1232,7 @@ namespace Step22
     if (plastic_iteration == 0)
       cell_viscosities.resize(0);
     //loop across all the cells to find and adjust eta of failing cells
-    for (unsigned int i = 0; i < triangulation.n_active_cells(); i++)
+    for (unsigned int i = 0; i < triangulation.n_active_cells(); ++i)
       {
         double current_cell_viscosity = 0;
 
@@ -1395,7 +1395,7 @@ namespace Step22
     if (total_fails <= 100 || decrease_in_plasticity <= 0.2)
       {
         system_parameters::continue_plastic_iterations = false;
-        for (unsigned int j=0; j < triangulation.n_active_cells(); j++)
+        for (unsigned int j=0; j < triangulation.n_active_cells(); ++j)
           {
             // Writes to file the undisturbed cell viscosities
             std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
@@ -1434,7 +1434,7 @@ namespace Step22
                   local_quadrature_points_history[q].new_eta = system_parameters::eta_floor;
 
                 quad_viscosities[cell_no][q].reinit(dim+1);
-                for (unsigned int ii=0; ii<dim; ii++)
+                for (unsigned int ii=0; ii<dim; ++ii)
                   quad_viscosities[cell_no][q](ii) = fe_values.quadrature_point(q)[ii];
                 quad_viscosities[cell_no][q](dim) = local_quadrature_points_history[q].new_eta;
               }
@@ -1506,7 +1506,7 @@ namespace Step22
             while (find_more_cells)
               {
                 new_cells_found = 0;
-                for (unsigned int i = start_cell; i<neighbor_cells.size(); i++)
+                for (unsigned int i = start_cell; i<neighbor_cells.size(); ++i)
                   {
                     for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
                       {
@@ -1539,11 +1539,11 @@ namespace Step22
             for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
               {
                 std::vector<double> nearby_etas_q;
-                for (unsigned int i = 0; i<neighbor_indices.size(); i++)
-                  for (unsigned int j=0; j<quadrature_formula.size(); j++)
+                for (unsigned int i = 0; i<neighbor_indices.size(); ++i)
+                  for (unsigned int j=0; j<quadrature_formula.size(); ++j)
                     {
                       Point<dim> test_q;
-                      for (unsigned int d=0; d<dim; d++)
+                      for (unsigned int d=0; d<dim; ++d)
                         test_q(d) = quad_viscosities[neighbor_indices[i]][j][d];
                       double qq_distance = fe_values.quadrature_point(q).distance(test_q);
                       if (qq_distance < system_parameters::smoothing_radius)
@@ -1551,7 +1551,7 @@ namespace Step22
                     }
                 // Write smoothed viscosities to quadrature_points_history; simple boxcar function is the smoothing kernel
                 double mean_eta = 0;
-                for (unsigned int l = 0; l<nearby_etas_q.size(); l++)
+                for (unsigned int l = 0; l<nearby_etas_q.size(); ++l)
                   {
                     mean_eta += nearby_etas_q[l];
                   }
@@ -1668,8 +1668,8 @@ namespace Step22
         if (cell->at_boundary())
           {
             int zero_faces = 0;
-            for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; f++)
-              for (unsigned int i=0; i<dim; i++)
+            for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+              for (unsigned int i=0; i<dim; ++i)
                 if (fabs(cell->face(f)->center()[i]) < zero_tolerance)
                   zero_faces++;
             if (zero_faces==0)
@@ -1758,7 +1758,7 @@ namespace Step22
     std::vector<double> ellipse_axes(0);
     // compute fit to boundary 0, 1, 2 ...
     std::cout << endl;
-    for (unsigned int i = 0; i<system_parameters::sizeof_material_id; i++)
+    for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
       {
         ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
         system_parameters::q_axes.push_back(ellipse_axes[0]);
@@ -1905,7 +1905,7 @@ namespace Step22
                 fout_boundaries.close();
 
                 how_many = 0;
-                for (unsigned int i=0; i<dim; i++)
+                for (unsigned int i=0; i<dim; ++i)
                   if (fabs(cell->face(f)->center()[i]) > zero_tolerance)
                     how_many++;
                 if (how_many==dim)
@@ -1925,7 +1925,7 @@ namespace Step22
     std::vector<double> ellipse_axes(0);
     // compute fit to boundary 0, 1, 2 ...
     std::cout << endl;
-    for (unsigned int i = 0; i<system_parameters::sizeof_material_id; i++)
+    for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
       {
         ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
         system_parameters::q_axes.push_back(ellipse_axes[0]);
diff --git a/CeresFE/support_code/config_in.h b/CeresFE/support_code/config_in.h
index 0ee79b9..1bafc21 100644
--- a/CeresFE/support_code/config_in.h
+++ b/CeresFE/support_code/config_in.h
@@ -122,22 +122,22 @@ namespace Step22
 
     // rheology parameters
 
-    for (unsigned int i=0; i<system_parameters::sizeof_depths_eta; i++)
+    for (unsigned int i=0; i<system_parameters::sizeof_depths_eta; ++i)
       fout_config << "depths_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
 
-    for (unsigned int i=0; i<system_parameters::sizeof_eta_kinks; i++)
+    for (unsigned int i=0; i<system_parameters::sizeof_eta_kinks; ++i)
       fout_config << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
 
-    for (unsigned int i=0; i<system_parameters::sizeof_depths_rho; i++)
+    for (unsigned int i=0; i<system_parameters::sizeof_depths_rho; ++i)
       fout_config << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
 
-    for (unsigned int i=0; i<system_parameters::sizeof_rho; i++)
+    for (unsigned int i=0; i<system_parameters::sizeof_rho; ++i)
       fout_config << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
 
-    for (unsigned int i=0; i<system_parameters::sizeof_material_id; i++)
+    for (unsigned int i=0; i<system_parameters::sizeof_material_id; ++i)
       fout_config << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
 
-    for (unsigned int i=0; i<system_parameters::sizeof_G; i++)
+    for (unsigned int i=0; i<system_parameters::sizeof_G; ++i)
       fout_config << "G[" << i << "] = " << system_parameters::G[i] << endl;
 
     fout_config << "eta_ceiling = " << system_parameters::eta_ceiling << endl;
@@ -262,7 +262,7 @@ namespace Step22
         unsigned int ndepths_eta = set_depths_eta.getLength();
         system_parameters::sizeof_depths_eta = ndepths_eta;
 
-        for (unsigned int i=0; i<ndepths_eta; i++)
+        for (unsigned int i=0; i<ndepths_eta; ++i)
           {
             system_parameters::depths_eta.push_back(set_depths_eta[i]);
             cout << "depth_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
@@ -276,7 +276,7 @@ namespace Step22
 
         //        cout << "Number of depth = " << ndepths << endl;
 
-        for (unsigned int i=0; i<neta_kinks; i++)
+        for (unsigned int i=0; i<neta_kinks; ++i)
           {
             system_parameters::eta_kinks.push_back(set_eta_kinks[i]);
             cout << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
@@ -290,7 +290,7 @@ namespace Step22
 
         //        cout << "Number of depth = " << ndepths << endl;
 
-        for (unsigned int i=0; i<ndepths_rho; i++)
+        for (unsigned int i=0; i<ndepths_rho; ++i)
           {
             system_parameters::depths_rho.push_back(set_depths_rho[i]);
             cout << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
@@ -304,7 +304,7 @@ namespace Step22
 
         //        cout << "Number of depth = " << ndepths << endl;
 
-        for (unsigned int i=0; i<nrho; i++)
+        for (unsigned int i=0; i<nrho; ++i)
           {
             system_parameters::rho.push_back(set_rho[i]);
             cout << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
@@ -318,7 +318,7 @@ namespace Step22
 
         //        cout << "Number of depth = " << ndepths << endl;
 
-        for (unsigned int i=0; i<nmaterial_id; i++)
+        for (unsigned int i=0; i<nmaterial_id; ++i)
           {
             system_parameters::material_id.push_back(set_material_id[i]);
             cout << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
@@ -332,7 +332,7 @@ namespace Step22
 
         //        cout << "Number of depth = " << ndepths << endl;
 
-        for (unsigned int i=0; i<nG; i++)
+        for (unsigned int i=0; i<nG; ++i)
           {
             system_parameters::G.push_back(set_G[i]);
             cout << "G[" << i << "] = " << system_parameters::G[i] << endl;
diff --git a/CeresFE/support_code/ellipsoid_fit.h b/CeresFE/support_code/ellipsoid_fit.h
index f935a23..5379916 100644
--- a/CeresFE/support_code/ellipsoid_fit.h
+++ b/CeresFE/support_code/ellipsoid_fit.h
@@ -155,10 +155,10 @@ void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char bou
       FullMatrix<double> A_out(ind_bnry_row.size(),dim);
       Vector<double>     b_out(ind_bnry_row.size());
 
-      for (unsigned int i=0; i<dim; i++)
+      for (unsigned int i=0; i<dim; ++i)
         ind_bnry_col.push_back(i);
 
-      for (unsigned int i=0; i<ind_bnry_row.size(); i++)
+      for (unsigned int i=0; i<ind_bnry_row.size(); ++i)
         b_out(i) = 1;
 
       A_out.extract_submatrix_from(A, ind_bnry_row, ind_bnry_col);
@@ -171,7 +171,7 @@ void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char bou
       solver.solve (AtA, x, Atb, PreconditionIdentity());
 
       // find ellipsoidal axes
-      for (unsigned int i=0; i<dim; i++)
+      for (unsigned int i=0; i<dim; ++i)
         ell.push_back(sqrt(1.0/x[i]));
     }
   else
diff --git a/Distributed_LDG_Method/LDGPoisson.cc b/Distributed_LDG_Method/LDGPoisson.cc
index 542c7e3..5db57a2 100644
--- a/Distributed_LDG_Method/LDGPoisson.cc
+++ b/Distributed_LDG_Method/LDGPoisson.cc
@@ -230,7 +230,7 @@ make_grid()
   triangulation.refine_global(n_refine);
 
   unsigned int local_refine = 2;
-  for (unsigned int i =0; i <local_refine; i++)
+  for (unsigned int i =0; i <local_refine; ++i)
     {
       typename Triangulation<dim>::active_cell_iterator
       cell = triangulation.begin_active(),
@@ -246,7 +246,7 @@ make_grid()
       // more realistic refinement stategies are
       // discussed elsewhere in the deal.ii
       // documentation.
-      for (; cell != endc; cell++)
+      for (; cell != endc; ++cell)
         {
           if ((cell->center()[1]) > 0.9 )
             {
@@ -273,7 +273,7 @@ make_grid()
   typename Triangulation<dim>::cell_iterator
   cell = triangulation.begin(),
   endc = triangulation.end();
-  for (; cell != endc; cell++)
+  for (; cell != endc; ++cell)
     {
       for (unsigned int face_no=0;
            face_no < GeometryInfo<dim>::faces_per_cell;
@@ -514,7 +514,7 @@ assemble_system()
   cell = dof_handler.begin_active(),
   endc = dof_handler.end();
 
-  for (; cell!=endc; cell++)
+  for (; cell!=endc; ++cell)
     {
       // Now, since we are working in a distributed setting,
       // we can only work on cells and write to dofs in the
@@ -856,16 +856,16 @@ assemble_cell_terms(
   // Now, we loop over the quadrature points in the
   // cell and then loop over the degrees of freedom and perform
   // quadrature to approximate the integrals.
-  for (unsigned int q=0; q<n_q_points; q++)
+  for (unsigned int q=0; q<n_q_points; ++q)
     {
-      for (unsigned int i=0; i<dofs_per_cell; i++)
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
         {
           const Tensor<1, dim> psi_i_field          = cell_fe[VectorField].value(i,q);
           const double         div_psi_i_field      = cell_fe[VectorField].divergence(i,q);
           const double         psi_i_potential      = cell_fe[Potential].value(i,q);
           const Tensor<1, dim> grad_psi_i_potential = cell_fe[Potential].gradient(i,q);
 
-          for (unsigned int j=0; j<dofs_per_cell; j++)
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
             {
               const Tensor<1, dim> psi_j_field        = cell_fe[VectorField].value(j,q);
               const double         psi_j_potential    = cell_fe[Potential].value(j,q);
@@ -926,14 +926,14 @@ assemble_Dirichlet_boundary_terms(
   Dirichlet_bc_function.value_list(face_fe.get_quadrature_points(),
                                    Dirichlet_bc_values);
 
-  for (unsigned int q=0; q<n_q_points; q++)
+  for (unsigned int q=0; q<n_q_points; ++q)
     {
-      for (unsigned int i=0; i<dofs_per_cell; i++)
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
         {
           const Tensor<1, dim> psi_i_field     = face_fe[VectorField].value(i,q);
           const double         psi_i_potential = face_fe[Potential].value(i,q);
 
-          for (unsigned int j=0; j<dofs_per_cell; j++)
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
             {
               const Tensor<1, dim> psi_j_field    = face_fe[VectorField].value(j,q);
               const double         psi_j_potential = face_fe[Potential].value(j,q);
@@ -994,14 +994,14 @@ assemble_Neumann_boundary_terms(
   // element faces.
   std::vector<double >    Neumann_bc_values(n_q_points);
 
-  for (unsigned int q=0; q<n_q_points; q++)
+  for (unsigned int q=0; q<n_q_points; ++q)
     {
-      for (unsigned int i=0; i<dofs_per_cell; i++)
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
         {
           const Tensor<1, dim> psi_i_field     = face_fe[VectorField].value(i,q);
           const double         psi_i_potential = face_fe[Potential].value(i,q);
 
-          for (unsigned int j=0; j<dofs_per_cell; j++)
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
             {
 
               const double    psi_j_potential = face_fe[Potential].value(j,q);
@@ -1057,23 +1057,23 @@ assemble_flux_terms(
   // the unit vector $\boldsymbol \beta$ that is used in defining
   // the LDG/ALternating fluxes.
   Point<dim> beta;
-  for (int i=0; i<dim; i++)
+  for (int i=0; i<dim; ++i)
     beta(i) = 1.0;
   beta /= sqrt(beta.square() );
 
   // Now we loop over all the quadrature points on the element face
   // and loop over all the degrees of freedom and approximate
   // the following flux integrals.
-  for (unsigned int q=0; q<n_face_points; q++)
+  for (unsigned int q=0; q<n_face_points; ++q)
     {
-      for (unsigned int i=0; i<dofs_this_cell; i++)
+      for (unsigned int i=0; i<dofs_this_cell; ++i)
         {
           const Tensor<1,dim>  psi_i_field_minus  =
             fe_face_values[VectorField].value(i,q);
           const double psi_i_potential_minus  =
             fe_face_values[Potential].value(i,q);
 
-          for (unsigned int j=0; j<dofs_this_cell; j++)
+          for (unsigned int j=0; j<dofs_this_cell; ++j)
             {
               const Tensor<1,dim> psi_j_field_minus   =
                 fe_face_values[VectorField].value(j,q);
@@ -1115,7 +1115,7 @@ assemble_flux_terms(
                                      fe_face_values.JxW(q);
             }
 
-          for (unsigned int j=0; j<dofs_neighbor_cell; j++)
+          for (unsigned int j=0; j<dofs_neighbor_cell; ++j)
             {
               const Tensor<1,dim> psi_j_field_plus    =
                 fe_neighbor_face_values[VectorField].value(j,q);
@@ -1159,14 +1159,14 @@ assemble_flux_terms(
             }
         }
 
-      for (unsigned int i=0; i<dofs_neighbor_cell; i++)
+      for (unsigned int i=0; i<dofs_neighbor_cell; ++i)
         {
           const Tensor<1,dim>  psi_i_field_plus =
             fe_neighbor_face_values[VectorField].value(i,q);
           const double         psi_i_potential_plus =
             fe_neighbor_face_values[Potential].value(i,q);
 
-          for (unsigned int j=0; j<dofs_this_cell; j++)
+          for (unsigned int j=0; j<dofs_this_cell; ++j)
             {
               const Tensor<1,dim> psi_j_field_minus               =
                 fe_face_values[VectorField].value(j,q);
@@ -1210,7 +1210,7 @@ assemble_flux_terms(
                                     fe_face_values.JxW(q);
             }
 
-          for (unsigned int j=0; j<dofs_neighbor_cell; j++)
+          for (unsigned int j=0; j<dofs_neighbor_cell; ++j)
             {
               const Tensor<1,dim> psi_j_field_plus =
                 fe_neighbor_face_values[VectorField].value(j,q);
@@ -1420,7 +1420,7 @@ output_results()    const
 
   Vector<float>   subdomain(triangulation.n_active_cells());
 
-  for (unsigned int i=0; i<subdomain.size(); i++)
+  for (unsigned int i=0; i<subdomain.size(); ++i)
     subdomain(i) = triangulation.locally_owned_subdomain();
 
   data_out.add_data_vector(subdomain,"subdomain");
diff --git a/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc b/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc
index c96f2b1..999f6a1 100644
--- a/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc
+++ b/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc
@@ -2556,9 +2556,9 @@ namespace NonLinearPoroViscoElasticity
                   {
                     const Tensor<2, dim> Grad_Nx_u =
                            scratch.fe_values_ref[u_fe].gradient(k, q);
-                    for (unsigned int dd = 0; dd < dim; dd++)
+                    for (unsigned int dd = 0; dd < dim; ++dd)
                       {
-                        for (unsigned int ee = 0; ee < dim; ee++)
+                        for (unsigned int ee = 0; ee < dim; ++ee)
                           {
                             scratch.solution_grads_u_total[q][dd][ee]
                              += scratch.local_dof_values[k] * Grad_Nx_u[dd][ee];
@@ -2573,7 +2573,7 @@ namespace NonLinearPoroViscoElasticity
 
                     scratch.solution_values_p_fluid_total[q]
                              += scratch.local_dof_values[k] * Nx_p;
-                    for (unsigned int dd = 0; dd < dim; dd++)
+                    for (unsigned int dd = 0; dd < dim; ++dd)
                       {
                         scratch.solution_grads_p_fluid_total[q][dd]
                             += scratch.local_dof_values[k] * Grad_Nx_p[dd];
diff --git a/goal_oriented_elastoplasticity/elastoplastic.cc b/goal_oriented_elastoplasticity/elastoplastic.cc
index b687c7f..cc3c5eb 100644
--- a/goal_oriented_elastoplasticity/elastoplastic.cc
+++ b/goal_oriented_elastoplasticity/elastoplastic.cc
@@ -4828,7 +4828,7 @@ namespace ElastoPlastic
           }
         else
           {
-            for (unsigned int i = 0; i < 5; i++)
+            for (unsigned int i = 0; i < 5; ++i)
               {
                 distributed_solution = incremental_displacement;
 
@@ -5469,8 +5469,8 @@ namespace ElastoPlastic
                          local_history_strain_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                          local_history_strain_fe_values (dim, std::vector< Vector<double> >(dim));
 
-    for (unsigned int i=0; i<dim; i++)
-      for (unsigned int j=0; j<dim; j++)
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
         {
           history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
           local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
@@ -5501,8 +5501,8 @@ namespace ElastoPlastic
           Assert (local_quadrature_points_history <
                   &quadrature_point_history.back(),
                   ExcInternalError());
-          for (unsigned int i=0; i<dim; i++)
-            for (unsigned int j=0; j<dim; j++)
+          for (unsigned int i=0; i<dim; ++i)
+            for (unsigned int j=0; j<dim; ++j)
               {
                 for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                   {
@@ -5614,8 +5614,8 @@ namespace ElastoPlastic
     // stress
     std::vector< std::vector< Vector<double> > >
     distributed_history_stress_field (dim, std::vector< Vector<double> >(dim));
-    for (unsigned int i=0; i<dim; i++)
-      for (unsigned int j=0; j<dim; j++)
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
         {
           distributed_history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
         }
@@ -5635,8 +5635,8 @@ namespace ElastoPlastic
     // strain
     std::vector< std::vector< Vector<double> > >
     distributed_history_strain_field (dim, std::vector< Vector<double> >(dim));
-    for (unsigned int i=0; i<dim; i++)
-      for (unsigned int j=0; j<dim; j++)
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
         {
           distributed_history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
         }
@@ -5679,8 +5679,8 @@ namespace ElastoPlastic
           Assert (local_quadrature_points_history <
                   &quadrature_point_history.back(),
                   ExcInternalError());
-          for (unsigned int i=0; i<dim; i++)
-            for (unsigned int j=0; j<dim; j++)
+          for (unsigned int i=0; i<dim; ++i)
+            for (unsigned int j=0; j<dim; ++j)
               {
                 dg_cell->get_dof_values (history_stress_field[i][j],
                                          local_history_stress_fe_values[i][j]);
@@ -6070,8 +6070,8 @@ namespace ElastoPlastic
     history_stress_field (dim, std::vector< Vector<double> >(dim)),
                          local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                          local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
-    for (unsigned int i=0; i<dim; i++)
-      for (unsigned int j=0; j<dim; j++)
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
         {
           history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
           local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
@@ -6113,8 +6113,8 @@ namespace ElastoPlastic
             {
               stress_at_qpoint = local_quadrature_points_history[q].old_stress;
 
-              for (unsigned int i=0; i<dim; i++)
-                for (unsigned int j=i; j<dim; j++)
+              for (unsigned int i=0; i<dim; ++i)
+                for (unsigned int j=i; j<dim; ++j)
                   {
                     local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
                   }
@@ -6124,8 +6124,8 @@ namespace ElastoPlastic
             }
 
 
-          for (unsigned int i=0; i<dim; i++)
-            for (unsigned int j=i; j<dim; j++)
+          for (unsigned int i=0; i<dim; ++i)
+            for (unsigned int j=i; j<dim; ++j)
               {
                 qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
                                             local_history_stress_values_at_qpoints[i][j]);
@@ -6154,8 +6154,8 @@ namespace ElastoPlastic
 
     std::vector< std::vector< Vector<double> > >
     history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
-    for (unsigned int i=0; i<dim; i++)
-      for (unsigned int j=0; j<dim; j++)
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
         {
           history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
         }
@@ -6175,8 +6175,8 @@ namespace ElastoPlastic
           dg_cell->get_dof_values (VM_stress_field,
                                    local_VM_stress_fe_values);
 
-          for (unsigned int i=0; i<dim; i++)
-            for (unsigned int j=0; j<dim; j++)
+          for (unsigned int i=0; i<dim; ++i)
+            for (unsigned int j=0; j<dim; ++j)
               {
                 dg_cell->get_dof_values (history_stress_field[i][j],
                                          local_history_stress_fe_values[i][j]);
@@ -6197,8 +6197,8 @@ namespace ElastoPlastic
 
               VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v);
 
-              for (unsigned int i=0; i<dim; i++)
-                for (unsigned int j=0; j<dim; j++)
+              for (unsigned int i=0; i<dim; ++i)
+                for (unsigned int j=0; j<dim; ++j)
                   {
                     history_stress_on_vertices[i][j](dof_1_vertex) +=
                       local_history_stress_fe_values[i][j](v);
@@ -6211,8 +6211,8 @@ namespace ElastoPlastic
       {
         VM_stress_on_vertices(id) /= counter_on_vertices(id);
 
-        for (unsigned int i=0; i<dim; i++)
-          for (unsigned int j=0; j<dim; j++)
+        for (unsigned int i=0; i<dim; ++i)
+          for (unsigned int j=0; j<dim; ++j)
             {
               history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
             }
@@ -6421,8 +6421,8 @@ namespace ElastoPlastic
         history_stress_field (dim, std::vector< Vector<double> >(dim)),
                              local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                              local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
-        for (unsigned int i=0; i<dim; i++)
-          for (unsigned int j=0; j<dim; j++)
+        for (unsigned int i=0; i<dim; ++i)
+          for (unsigned int j=0; j<dim; ++j)
             {
               history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
               local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
@@ -6473,8 +6473,8 @@ namespace ElastoPlastic
                   // stress in polar coordinate
                   stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix);
 
-                  for (unsigned int i=0; i<dim; i++)
-                    for (unsigned int j=i; j<dim; j++)
+                  for (unsigned int i=0; i<dim; ++i)
+                    for (unsigned int j=i; j<dim; ++j)
                       {
                         local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
                       }
@@ -6482,8 +6482,8 @@ namespace ElastoPlastic
                 }
 
 
-              for (unsigned int i=0; i<dim; i++)
-                for (unsigned int j=i; j<dim; j++)
+              for (unsigned int i=0; i<dim; ++i)
+                for (unsigned int j=i; j<dim; ++j)
                   {
                     qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
                                                 local_history_stress_values_at_qpoints[i][j]);
@@ -6545,8 +6545,8 @@ namespace ElastoPlastic
 
         std::vector< std::vector< Vector<double> > >
         history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
-        for (unsigned int i=0; i<dim; i++)
-          for (unsigned int j=0; j<dim; j++)
+        for (unsigned int i=0; i<dim; ++i)
+          for (unsigned int j=0; j<dim; ++j)
             {
               history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
             }
@@ -6562,8 +6562,8 @@ namespace ElastoPlastic
           if (cell->is_locally_owned())
             {
 
-              for (unsigned int i=0; i<dim; i++)
-                for (unsigned int j=0; j<dim; j++)
+              for (unsigned int i=0; i<dim; ++i)
+                for (unsigned int j=0; j<dim; ++j)
                   {
                     dg_cell->get_dof_values (history_stress_field[i][j],
                                              local_history_stress_fe_values[i][j]);
@@ -6582,8 +6582,8 @@ namespace ElastoPlastic
 
                   counter_on_vertices (dof_1_vertex) += 1;
 
-                  for (unsigned int i=0; i<dim; i++)
-                    for (unsigned int j=0; j<dim; j++)
+                  for (unsigned int i=0; i<dim; ++i)
+                    for (unsigned int j=0; j<dim; ++j)
                       {
                         history_stress_on_vertices[i][j](dof_1_vertex) +=
                           local_history_stress_fe_values[i][j](v);
@@ -6594,8 +6594,8 @@ namespace ElastoPlastic
 
         for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
           {
-            for (unsigned int i=0; i<dim; i++)
-              for (unsigned int j=0; j<dim; j++)
+            for (unsigned int i=0; i<dim; ++i)
+              for (unsigned int j=0; j<dim; ++j)
                 {
                   history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
                 }
diff --git a/two_phase_flow/LevelSetSolver.cc b/two_phase_flow/LevelSetSolver.cc
index d2c445e..6f36032 100644
--- a/two_phase_flow/LevelSetSolver.cc
+++ b/two_phase_flow/LevelSetSolver.cc
@@ -629,7 +629,7 @@ void LevelSetSolver<dim>::invert_ML()
 {
   // loop on locally owned i-DOFs (rows)
   IndexSet::ElementIterator idofs_iter = locally_owned_dofs_LS.begin();
-  for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+  for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
     {
       int gi = *idofs_iter;
       inverse_ML_vector(gi) = 1./ML_vector(gi);
@@ -747,7 +747,7 @@ void LevelSetSolver<dim>::assemble_C_Matrix ()
               }
 
             for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
-              for (unsigned int j=0; j < dofs_per_cell_LS; j++)
+              for (unsigned int j=0; j < dofs_per_cell_LS; ++j)
                 {
                   cell_Cij_x(i,j) += (shape_grads_LS[j][0])*shape_values_LS[i]*JxW;
                   cell_Cij_y(i,j) += (shape_grads_LS[j][1])*shape_values_LS[i]*JxW;
@@ -874,7 +874,7 @@ void LevelSetSolver<dim>::assemble_K_DL_DH_times_vector
   // loop on locally owned i-DOFs (rows)
   IndexSet::ElementIterator idofs_iter = locally_owned_dofs_LS.begin();
 
-  for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+  for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
     {
       PetscInt gi = *idofs_iter;
       //double ith_K_times_solution = 0;
@@ -915,7 +915,7 @@ void LevelSetSolver<dim>::assemble_K_DL_DH_times_vector
       std::vector<double> dLi(ncolumns), dCi(ncolumns);
       double dLii = 0, dCii = 0;
       // loop on sparsity pattern of i-th DOF
-      for (int j =0; j < ncolumns; j++)
+      for (int j =0; j < ncolumns; ++j)
         {
           C[0] = Cxi[j];
           C[1] = Cyi[j];
@@ -1073,7 +1073,7 @@ void LevelSetSolver<dim>::assemble_EntRes_Matrix ()
               }
 
             for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
-              for (unsigned int j=0; j < dofs_per_cell_LS; j++)
+              for (unsigned int j=0; j < dofs_per_cell_LS; ++j)
                 {
                   cell_EntRes (i,j) += Rk*shape_values_LS[i]*shape_values_LS[j]*JxW;
                   cell_volume (i,j) += JxW;
@@ -1113,7 +1113,7 @@ void LevelSetSolver<dim>::compute_bounds(PETScWrappers::MPI::Vector &un_solution
   umax_vector = 0;
   // loop on locally owned i-DOFs (rows)
   IndexSet::ElementIterator idofs_iter = locally_owned_dofs_LS.begin();
-  for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+  for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
     {
       int gi = *idofs_iter;
 
@@ -1123,7 +1123,7 @@ void LevelSetSolver<dim>::compute_bounds(PETScWrappers::MPI::Vector &un_solution
       get_vector_values(un_solution,gj_indices,soln);
       // compute bounds, ith row of flux matrix, P vectors
       double mini=1E10, maxi=-1E10;
-      for (unsigned int j =0; j < gj_indices.size(); j++)
+      for (unsigned int j =0; j < gj_indices.size(); ++j)
         {
           // bounds
           mini = std::min(mini,soln[j]);
@@ -1152,7 +1152,7 @@ void LevelSetSolver<dim>::check_max_principle(PETScWrappers::MPI::Vector &unp1_s
     if (cell_LS->is_locally_owned() && !cell_LS->at_boundary())
       {
         cell_LS->get_dof_indices(local_dof_indices);
-        for (unsigned int i=0; i<dofs_per_cell; i++)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
           if (locally_owned_dofs_LS.is_element(local_dof_indices[i]))
             {
               double solni = unp1_solution(local_dof_indices[i]);
@@ -1212,7 +1212,7 @@ void LevelSetSolver<dim>::compute_MPP_uH
   const PetscScalar *MCi, *dLi, *dCi;
   double solni, mi, solLi, solHi;
 
-  for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+  for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
     {
       int gi = *idofs_iter;
       // read vectors at i-th DOF
@@ -1238,7 +1238,7 @@ void LevelSetSolver<dim>::compute_MPP_uH
       // compute bounds, ith row of flux matrix, P vectors
       double mini=1E10, maxi=-1E10;
       double Pposi=0 ,Pnegi=0;
-      for (int j =0; j < ncolumns; j++)
+      for (int j =0; j < ncolumns; ++j)
         {
           // bounds
           mini = std::min(mini,soln[j]);
@@ -1281,7 +1281,7 @@ void LevelSetSolver<dim>::compute_MPP_uH
   const double *Ai;
   double Rposi, Rnegi;
   idofs_iter=locally_owned_dofs_LS.begin();
-  for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+  for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
     {
       int gi = *idofs_iter;
       Rposi = R_pos_vector(gi);
@@ -1300,7 +1300,7 @@ void LevelSetSolver<dim>::compute_MPP_uH
       // Array for i-th row of A_times_L matrix
       std::vector<double> LxAi(ncolumns);
       // loop in sparsity pattern of i-th DOF
-      for (int j =0; j < ncolumns; j++)
+      for (int j =0; j < ncolumns; ++j)
         LxAi[j] = Ai[j] * ((Ai[j]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]));
 
       // save i-th row of LxA
@@ -1334,14 +1334,14 @@ void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
       const PetscInt *gj;
       const PetscScalar *Akp1i;
       double mi;
-      for (int iter=0; iter<NUM_ITER; iter++)
+      for (int iter=0; iter<NUM_ITER; ++iter)
         {
           MPP_uLkp1_solution_ghosted = MPP_uH_solution;
           Akp1_matrix.add(-1.0, LxAkp1_matrix); //new matrix to limit: A-LxA
 
           // loop on locally owned i-DOFs (rows)
           IndexSet::ElementIterator idofs_iter = locally_owned_dofs_LS.begin();
-          for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+          for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
             {
               int gi = *idofs_iter;
 
@@ -1359,7 +1359,7 @@ void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
               // compute bounds, ith row of flux matrix, P vectors
               double mini=1E10, maxi=-1E10;
               double Pposi=0 ,Pnegi=0;
-              for (int j =0; j < ncolumns; j++)
+              for (int j =0; j < ncolumns; ++j)
                 {
                   // bounds
                   mini = std::min(mini,soln[j]);
@@ -1390,7 +1390,7 @@ void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
           // compute limiters. NOTE: this is a different loop due to need of i- and j-th entries of R vectors
           double Rposi, Rnegi;
           idofs_iter=locally_owned_dofs_LS.begin();
-          for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+          for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
             {
               int gi = *idofs_iter;
               Rposi = R_pos_vector(gi);
@@ -1408,7 +1408,7 @@ void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
 
               // Array for i-th row of LxAkp1 matrix
               std::vector<double> LxAkp1i(ncolumns);
-              for (int j =0; j < ncolumns; j++)
+              for (int j =0; j < ncolumns; ++j)
                 LxAkp1i[j] = Akp1i[j] * ((Akp1i[j]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]));
 
               // save i-th row of LxA
@@ -1496,7 +1496,7 @@ void LevelSetSolver<dim>::get_sparsity_pattern()
   const PetscInt *gj;
   const PetscScalar *MCi;
 
-  for (; idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++)
+  for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
     {
       PetscInt gi = *idofs_iter;
       // get i-th row of mass matrix (dummy, I just need the indices gj)
@@ -1595,7 +1595,7 @@ void LevelSetSolver<dim>::get_vector_values (PETScWrappers::VectorBase &vector,
   PetscScalar *soln;
   VecGetArray(solution_in_local_form, &soln);
 
-  for (i = 0; i < n_idx; i++)
+  for (i = 0; i < n_idx; ++i)
     {
       int index = indices[i];
       if (index >= begin && index < end)
@@ -1641,7 +1641,7 @@ void LevelSetSolver<dim>::get_vector_values (PETScWrappers::VectorBase &vector,
   PetscScalar *soln;
   VecGetArray(solution_in_local_form, &soln);
 
-  for (i = 0; i < n_idx; i++)
+  for (i = 0; i < n_idx; ++i)
     {
       int index = map_from_Q1_to_Q2[indices[i]];
       if (index >= begin && index < end)
-- 
2.39.5