Clean up velocity stepping test

src/one-body-sample.cc

@@ -215,6 +215,8 @@ int main(int argc, char *argv[]) {
 
     std::fstream octave_writer("data", std::fstream::out);
     std::fstream coefficient_writer("coefficient", std::fstream::out);
+    std::fstream velocity_stepping_writer("velocity_stepping",
+                                          std::fstream::out);
     timer.reset();
 
     auto const L = parset.get<double>("boundary.friction.ruina.L");
@@ -229,6 +231,10 @@ int main(int argc, char *argv[]) {
                   << "# rows: " << timesteps << std::endl << "# columns: 3"
                   << std::endl;
 
+    velocity_stepping_writer << "# name: B" << std::endl << "# type: matrix"
+                             << std::endl << "# rows: " << timesteps
+                             << std::endl << "# columns: 2" << std::endl;
+
     auto const &dirichletFunction = functions.get("dirichletCondition");
     auto const &neumannFunction = functions.get("neumannCondition");
 
@@ -329,25 +335,24 @@ int main(int argc, char *argv[]) {
         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;
+          double const euler = (*s4_new)[first_frictional_node];
+          double direct;
+          if (run < 120) {
+            direct = euler;
+          } else 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)));
           }
+          velocity_stepping_writer << euler << " " << direct << std::endl;
         }
 
         // Record the coefficient of friction at a fixed node
@@ -392,6 +397,7 @@ int main(int argc, char *argv[]) {
 
     octave_writer.close();
     coefficient_writer.close();
+    velocity_stepping_writer.close();
 
     Python::stop();
   }