Define all possible variables const in FV local operator

Const all the things!

dune/dorie/model/richards/local_operator_FV.hh

 #ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
 #define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
 
-#include <vector>
-
 #include <dune/geometry/referenceelements.hh>
 
 #include <dune/pdelab/common/quadraturerules.hh>
@@ -114,19 +112,19 @@ public:
     const auto& entity_o = ig.outside();
 
     // Get geometries
-    auto geo_f = entity_f.geometry();
-    auto geo_i = entity_i.geometry();
-    auto geo_o = entity_o.geometry();
+    const auto geo_f = entity_f.geometry();
+    const auto geo_i = entity_i.geometry();
+    const auto geo_o = entity_o.geometry();
 
     // Get volume of entities
     const auto volume_f = geo_f.volume();
 
     // Get normal
-    auto normal_f = entity_f.centerUnitOuterNormal();
+    const auto normal_f = entity_f.centerUnitOuterNormal();
 
     // Entity centers in global coordinates
-    auto center_position_i_g = geo_i.center();
-    auto center_position_o_g = geo_o.center();
+    const auto center_position_i_g = geo_i.center();
+    const auto center_position_o_g = geo_o.center();
 
     // Inside/outside solute value
     const RangeU u_i = x_i(lfsu_i,0);
@@ -150,8 +148,8 @@ public:
                                 + 1.0/(cond_o + 1E-30) );
 
     // Distance between the two entity centers
-    center_position_i_g -= center_position_o_g;
-    auto distance = center_position_i_g.two_norm();
+    const auto center_position_diff = center_position_i_g - center_position_o_g;
+    const auto distance = center_position_diff.two_norm();
 
     // Finite difference of u between the two entities
     RangeU dudn = (u_o-u_i)/distance;
@@ -210,9 +208,17 @@ public:
     const auto& entity_i = ig.inside();
 
     // Get geometries
-    auto geo_f = entity_f.geometry();
+    const auto geo_f = entity_f.geometry();
     const auto& ref_el_f = referenceElement(geo_f);
     const auto& center_position_f = ref_el_f.position(0,0);
+
+    // Face volume for integration
+    const auto volume_f = geo_f.volume();
+
+    // Normal vector of intersection
+    const auto normal_f = ig.centerUnitOuterNormal();
+
+    // Retrieve boundary condition
     auto bc = _boundary->bc(entity_f,center_position_f,_time);
 
     // bind parameterization and retrieve functions
@@ -222,27 +228,21 @@ public:
 
     if (BoundaryCondition::isDirichlet(bc))
     {
-      auto geo_i = entity_i.geometry();
-      auto center_position_i_g = geo_i.center();
-
-      // Get normal
-      auto normal_f = entity_f.centerUnitOuterNormal();
+      const auto geo_i = entity_i.geometry();
+      const auto center_position_i_g = geo_i.center();
 
       // Face center in global coordinates
-      auto center_position_f_g = geo_f.center();
+      const auto center_position_f_g = geo_f.center();
 
       // Compute distance of these two points
-      center_position_i_g -= center_position_f_g;
-      auto distance = center_position_i_g.two_norm();
+      const auto position_diff = center_position_i_g - center_position_f_g;
+      const auto distance = position_diff.two_norm();
 
       // Evaluate Dirichlet condition
-      auto g = _boundary->g(entity_f,center_position_f,_time);
-
-      // Face volume for integration
-      auto volume_f = geo_f.volume();
+      const auto g = _boundary->g(entity_f,center_position_f,_time);
 
       // Inside unknown value
-      RangeU u_i = x_i(lfsu_i,0);
+      const RangeU u_i = x_i(lfsu_i,0);
 
       // Compute the conductivity
       const RangeU cond_i = conductivity_f_i(saturation_f_i(u_i));
@@ -254,21 +254,18 @@ public:
       dudn -= _param->gravity()*normal_f;
 
       // Water flux in normal direction w.r.t the intersection
-      auto water_flux_n = - cond_i*dudn;
+      const auto water_flux_n = - cond_i*dudn;
 
       // Contribution to residual from Dirichlet boundary
       r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f);
     }
     else if (BoundaryCondition::isNeumann(bc))
     {
-      // Face volume for integration
-      auto volume_f = geo_f.volume();
-
       // Contribution to residual from Neumann boundary
       RangeFieldU water_flux_f = _boundary->j(entity_f,center_position_f,_time);
 
       // Convert water flux into units of matric head flux
-      water_flux_f *= std::abs( ig.centerUnitOuterNormal()*_param->gravity() );
+      water_flux_f *= std::abs( normal_f * _param->gravity() );
       // water_flux_f *= _param->gravity().two_norm();
 
       r_i.accumulate(lfsv_i, 0, water_flux_f*volume_f );
@@ -362,11 +359,10 @@ public:
                         Traits::LocalBasisType::Traits;
     // Get range for trial space
     using RangeU = typename LocalBasisTraitsU::RangeType;
-
-    RangeU u = x(lfsu,0);
+    const RangeU u = x(lfsu,0);
 
     // entity geometry
-    auto geo = eg.geometry();
+    const auto geo = eg.geometry();
 
     // bind parameterization and retrieve functions
     _param->bind(eg.entity());