Work with velocities rather than differences

src/one-body-sample.cc

@@ -194,10 +194,10 @@ int main(int argc, char *argv[]) {
     VectorType u_old(finestSize);
     u_old = 0.0; // Has to be zero!
     VectorType u_old_old;
-    VectorType u_diff(finestSize);
-    u_diff = 0.0; // Has to be zero!
 
     VectorType u(finestSize);
+    VectorType ud(finestSize);
+    ud = 0.0;
 
     SingletonVectorType alpha_old(finestSize);
     alpha_old = parset.get<double>("boundary.friction.state.initial");
@@ -298,12 +298,12 @@ int main(int argc, char *argv[]) {
           overallSolver.solve();
 
           timeSteppingScheme->extractSolution(problem_iterate, u);
-          timeSteppingScheme->extractDifference(problem_iterate, u_diff);
+          timeSteppingScheme->extractVelocity(problem_iterate, ud);
 
           // Update the state
           for (size_t i = 0; i < frictionalNodes.size(); ++i) {
             if (frictionalNodes[i][0]) {
-              double const unorm = u_diff[i].two_norm();
+              double const unorm = ud[i].two_norm() * tau;
 
               // // the (logarithmic) steady state corresponding to the
               // // current velocity
@@ -353,7 +353,7 @@ 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 = u_diff[first_frictional_node].two_norm() / tau / L;
+          double const v = ud[first_frictional_node].two_norm() / L;
 
           double const euler = alpha[first_frictional_node];
           double direct;
@@ -377,7 +377,7 @@ int main(int argc, char *argv[]) {
 
         // Record the coefficient of friction at a fixed node
         if (parset.get<bool>("printCoefficient")) {
-          double const V = u_diff[first_frictional_node].two_norm();
+          double const V = ud[first_frictional_node].two_norm();
           double const state = alpha[first_frictional_node];
 
           coefficient_writer << (mu + a * std::log(V * eta) +
@@ -388,7 +388,7 @@ int main(int argc, char *argv[]) {
       if (parset.get<bool>("printFrictionalVelocity")) {
         for (size_t i = 0; i < frictionalNodes.size(); ++i)
           if (frictionalNodes[i][0])
-            std::cout << u_diff[i][0] * timesteps << " ";
+            std::cout << ud[i][0] << " ";
         std::cout << std::endl;
       }
       u_old_old = u_old;

src/timestepping.cc

@@ -23,8 +23,8 @@ class TimeSteppingScheme {
   void virtual extractSolution(VectorType const &problem_iterate,
                                VectorType &solution) const = 0;
 
-  void virtual extractDifference(VectorType const &problem_iterate,
-                                 VectorType &difference) const = 0;
+  void virtual extractVelocity(VectorType const &problem_iterate,
+                               VectorType &velocity) const = 0;
 
 protected:
   VectorType const ℓ
@@ -84,10 +84,9 @@ class ImplicitEuler
     solution.axpy(this->tau, problem_iterate);
   }
 
-  void virtual extractDifference(VectorType const &problem_iterate,
-                                 VectorType &difference) const {
-    difference = problem_iterate;
-    difference *= this->tau;
+  void virtual extractVelocity(VectorType const &problem_iterate,
+                               VectorType &velocity) const {
+    velocity = problem_iterate;
   }
 };
 
@@ -151,9 +150,8 @@ class ImplicitTwoStep
     }
   }
 
-  void virtual extractDifference(VectorType const &problem_iterate,
-                                 VectorType &difference) const {
-    difference = problem_iterate;
-    difference *= this->tau;
+  void virtual extractVelocity(VectorType const &problem_iterate,
+                               VectorType &velocity) const {
+    velocity = problem_iterate;
   }
 };