Velocity stepping tests

src/one-body-sample.cc

@@ -323,6 +323,33 @@ int main(int argc, char *argv[]) {
                         << std::endl;
         }
 
+        // Comparison with the analytic solution of a velocity stepping test
+        // 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;
+
+          if (run >= 120) {
+            double const euler = (*s4_new)[first_frictional_node];
+            double direct;
+            if (run < 360) {
+              double const v2 = v;
+              double const v1 = 0.6 * v2;
+              direct = std::log(
+                  1.0 / v2 *
+                  std::pow((v2 / v1), std::exp(-v2 * (run - 120) * h)));
+            } else {
+              double const v1 = v;
+              double const v2 = v1 / 0.6;
+              direct = std::log(
+                  1.0 / v1 *
+                  std::pow((v1 / v2), std::exp(-v1 * (run - 360) * h)));
+            }
+            std::cout << "[" << run << "]: " << euler << " " << direct << " "
+                      << (euler - direct) / direct << std::endl;
+          }
+        }
+
         // Record the coefficient of friction at a fixed node
         if (parset.get<bool>("printCoefficient")) {
           double const V = u4_diff[first_frictional_node].two_norm();

src/one-body-sample.parset

@@ -9,6 +9,8 @@ printProgress = true
 writeEvolution = false
 writeVTK = false
 
+printVelocitySteppingComparison = false
+
 [grid]
 refinements = 4