From dcebc8a1e1f390f63be4b8326bf8ad14724c66ee Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 16 May 2012 14:41:38 +0200
Subject: [PATCH] Clean up velocity stepping test
---
src/one-body-sample.cc | 42 ++++++++++++++++++++++++------------------
1 file changed, 24 insertions(+), 18 deletions(-)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 45574e24..0f7b7c38 100644
--- a/src/one-body-sample.cc
+++ b/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();
}
--
GitLab