u4 -> u, s4 -> alpha, b -> rhs

src/one-body-sample.cc

@@ -189,20 +189,20 @@ int main(int argc, char *argv[]) {
             leafView, p1Basis, frictionalNodes);
 
     // {{{ Initialise vectors
-    VectorType u4(finestSize);
-    u4 = 0.0; // Has to be zero!
+    VectorType u(finestSize);
+    u = 0.0; // Has to be zero!
 
-    SingletonVectorType s4_old(finestSize);
-    s4_old = parset.get<double>("boundary.friction.state.initial");
+    SingletonVectorType alpha_old(finestSize);
+    alpha_old = parset.get<double>("boundary.friction.state.initial");
 
-    VectorType u4_diff(finestSize);
-    u4_diff = 0.0; // Has to be zero!
+    VectorType u_diff(finestSize);
+    u_diff = 0.0; // Has to be zero!
 
-    auto s4_new = Dune::make_shared<SingletonVectorType>(finestSize);
-    *s4_new = s4_old;
+    auto alpha_new = Dune::make_shared<SingletonVectorType>(finestSize);
+    *alpha_new = alpha_old;
 
     SingletonVectorType vonMisesStress;
-    VectorType b4;
+    VectorType rhs;
     // }}}
 
     typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
@@ -253,13 +253,13 @@ int main(int argc, char *argv[]) {
 
       {
         assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
-            leafView, p1Basis, neumannNodes, b4, neumannFunction, time);
-        stiffnessMatrix.mmv(u4, b4);
+            leafView, p1Basis, neumannNodes, rhs, neumannFunction, time);
+        stiffnessMatrix.mmv(u, rhs);
         // Apply Dirichlet condition
         for (size_t i = 0; i < finestSize; ++i)
           if (ignoreNodes[i].count() == dim) {
-            dirichletFunction.evaluate(time, u4_diff[i][0]);
-            u4_diff[i][0] /= refinement_factor;
+            dirichletFunction.evaluate(time, u_diff[i][0]);
+            u_diff[i][0] /= refinement_factor;
           }
 
         for (size_t state_fpi = 0;
@@ -268,15 +268,16 @@ int main(int argc, char *argv[]) {
              ++state_fpi) {
           auto myGlobalNonlinearity =
               assemble_nonlinearity<MatrixType, VectorType>(
-                  parset.sub("boundary.friction"), nodalIntegrals, s4_new, h);
+                  parset.sub("boundary.friction"), nodalIntegrals, alpha_new,
+                  h);
           MyConvexProblemType const myConvexProblem(stiffnessMatrix,
-                                                    *myGlobalNonlinearity, b4);
+                                                    *myGlobalNonlinearity, rhs);
 
           MyBlockProblemType myBlockProblem(parset, myConvexProblem);
           auto multigridStep = mySolver.getSolver();
-          multigridStep->setProblem(u4_diff, myBlockProblem);
+          multigridStep->setProblem(u_diff, myBlockProblem);
 
-          VectorType const u4_diff_old = u4_diff;
+          VectorType const u_diff_old = u_diff;
           LoopSolver<VectorType> overallSolver(
               multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
               solver_tolerance, &energyNorm, verbosity,
@@ -285,7 +286,7 @@ int main(int argc, char *argv[]) {
 
           for (size_t i = 0; i < frictionalNodes.size(); ++i) {
             if (frictionalNodes[i][0]) {
-              double const unorm = u4_diff[i].two_norm();
+              double const unorm = u_diff[i].two_norm();
 
               // // the (logarithmic) steady state corresponding to the
               // // current velocity
@@ -294,11 +295,12 @@ int main(int argc, char *argv[]) {
               switch (parset.get<Config::state_model>(
                   "boundary.friction.state.model")) {
                 case Config::Dieterich:
-                  (*s4_new)[i] =
-                      state_update_dieterich(h, unorm / L, s4_old[i]);
+                  (*alpha_new)[i] =
+                      state_update_dieterich(h, unorm / L, alpha_old[i]);
                   break;
                 case Config::Ruina:
-                  (*s4_new)[i] = state_update_ruina(h, unorm / L, s4_old[i]);
+                  (*alpha_new)[i] =
+                      state_update_ruina(h, unorm / L, alpha_old[i]);
                   break;
               }
             }
@@ -307,7 +309,7 @@ int main(int argc, char *argv[]) {
             std::cout << '.';
             std::cout.flush();
           }
-          if (energyNorm.diff(u4_diff_old, u4_diff) <
+          if (energyNorm.diff(u_diff_old, u_diff) <
               parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
             break;
         }
@@ -317,8 +319,8 @@ int main(int argc, char *argv[]) {
         if (parset.get<bool>("writeEvolution")) {
           double out;
           neumannFunction.evaluate(time, out);
-          octave_writer << (*s4_new)[first_frictional_node][0] << " "
-                        << u4[first_frictional_node][0] * 1e6 << " " << out
+          octave_writer << (*alpha_new)[first_frictional_node][0] << " "
+                        << u[first_frictional_node][0] * 1e6 << " " << out
                         << std::endl;
         }
 
@@ -326,9 +328,9 @@ int main(int argc, char *argv[]) {
         // with the Ruina state evolution law.
         // Jumps at 120 and 360 timesteps; v1 = .6 * v2;
         if (parset.get<bool>("printVelocitySteppingComparison")) {
-          double const v = u4_diff[first_frictional_node].two_norm() / h / L;
+          double const v = u_diff[first_frictional_node].two_norm() / h / L;
 
-          double const euler = (*s4_new)[first_frictional_node];
+          double const euler = (*alpha_new)[first_frictional_node];
           double direct;
           if (run < 120) {
             direct = euler;
@@ -350,28 +352,28 @@ int main(int argc, char *argv[]) {
 
         // Record the coefficient of friction at a fixed node
         if (parset.get<bool>("printCoefficient")) {
-          double const V = u4_diff[first_frictional_node].two_norm();
-          double const state = (*s4_new)[first_frictional_node];
+          double const V = u_diff[first_frictional_node].two_norm();
+          double const state = (*alpha_new)[first_frictional_node];
 
           coefficient_writer << (mu + a * std::log(V * eta) +
                                  b * (state - std::log(eta * L))) << std::endl;
         }
       }
 
-      u4 += u4_diff;
-      s4_old = *s4_new;
+      u += u_diff;
+      alpha_old = *alpha_new;
 
       // Compute von Mises stress and write everything to a file
       if (parset.get<bool>("writeVTK")) {
         VonMisesStressAssembler<GridType> localStressAssembler(
             E, nu,
             Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>(
-                p1Basis, u4));
+                p1Basis, u));
         FunctionalAssembler<P0Basis>(p0Basis)
             .assemble(localStressAssembler, vonMisesStress, true);
 
         writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>(
-            p1Basis, u4, *s4_new, p0Basis, vonMisesStress, leafView,
+            p1Basis, u, *alpha_new, p0Basis, vonMisesStress, leafView,
             (boost::format("obs%d") % run).str());
       }
     }
@@ -381,10 +383,10 @@ int main(int argc, char *argv[]) {
 
     if (parset.get<bool>("printFrictionalBoundary")) {
       // Print displacement on frictional boundary
-      boost::format const formatter("u4[%02d] = %+3e");
+      boost::format const formatter("u[%02d] = %+3e");
       for (size_t i = 0; i < frictionalNodes.size(); ++i)
         if (frictionalNodes[i][0])
-          std::cout << (boost::format(formatter) % i % u4[i]) << std::endl;
+          std::cout << (boost::format(formatter) % i % u[i]) << std::endl;
     }
 
     octave_writer.close();