From d4f163494bc0381c9bf7a84eacf6d91680f165a2 Mon Sep 17 00:00:00 2001 From: Elias Pipping <elias.pipping@fu-berlin.de> Date: Mon, 27 Aug 2012 10:07:34 +0200 Subject: [PATCH] Rename: h -> tau --- src/assemblers.cc | 4 ++-- src/assemblers_tmpl.cc | 2 +- src/compute_state.cc | 22 +++++++++++----------- src/compute_state_ruina.cc | 4 ++-- src/one-body-sample.cc | 27 ++++++++++++++------------- 5 files changed, 30 insertions(+), 29 deletions(-) diff --git a/src/assemblers.cc b/src/assemblers.cc index 1c5d5ed0..cf215ea3 100644 --- a/src/assemblers.cc +++ b/src/assemblers.cc @@ -57,7 +57,7 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const> assemble_nonlinearity( Dune::ParameterTree const &parset, Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals, - Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h) { + Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double tau) { auto const size = nodalIntegrals.size(); typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType; @@ -84,7 +84,7 @@ assemble_nonlinearity( return Dune::make_shared< Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>( - nodalIntegrals, a, mu, eta, normalStress, b, state, L, h); + nodalIntegrals, a, mu, eta, normalStress, b, state, L, tau); } case Config::Laursen: assert(false); diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc index 388e60cf..d4e815c9 100644 --- a/src/assemblers_tmpl.cc +++ b/src/assemblers_tmpl.cc @@ -34,4 +34,4 @@ template Dune::shared_ptr< assemble_nonlinearity<MatrixType, VectorType>( Dune::ParameterTree const &parset, Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals, - Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h); + Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double tau); diff --git a/src/compute_state.cc b/src/compute_state.cc index 2ff511c2..b928618b 100644 --- a/src/compute_state.cc +++ b/src/compute_state.cc @@ -17,11 +17,11 @@ class DecayingExponential { typedef Dune::FieldVector<double, 1> VectorType; typedef Dune::FieldMatrix<double, 1, 1> MatrixType; - DecayingExponential(double h) : h(h) {} + DecayingExponential(double tau) : tau(tau) {} void directionalSubDiff(VectorType const &u, VectorType const &v, Interval<double> &D) const { - D[0] = D[1] = v[0] * (-h * std::exp(-u[0])); + D[0] = D[1] = v[0] * (-tau * std::exp(-u[0])); } void directionalDomain(VectorType const &, VectorType const &, @@ -31,14 +31,14 @@ class DecayingExponential { } private: - double const h; + double const tau; }; -double state_update_dieterich_bisection(double h, double uol, +double state_update_dieterich_bisection(double tau, double uol, double old_state) { DecayingExponential::VectorType const start(0); DecayingExponential::VectorType const direction(1); - DecayingExponential const phi(h); + DecayingExponential const phi(tau); MyDirectionalConvexFunction<DecayingExponential> const J( 1.0, old_state - uol, phi, start, direction); @@ -47,18 +47,18 @@ double state_update_dieterich_bisection(double h, double uol, return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO } -double state_update_dieterich(double h, double uol, double old_state) { - double const ret = state_update_dieterich_bisection(h, uol, old_state); +double state_update_dieterich(double tau, double uol, double old_state) { + double const ret = state_update_dieterich_bisection(tau, uol, old_state); /* We have - ret - old_state + uol = h*exp(-ret) + ret - old_state + uol = tau*exp(-ret) or - log((ret - old_state + uol)/h) = -ret + log((ret - old_state + uol)/tau) = -ret */ - assert(std::min(std::abs(ret - old_state + uol - h * std::exp(-ret)), - std::abs(std::log((ret - old_state + uol) / h) + ret)) < + assert(std::min(std::abs(ret - old_state + uol - tau * std::exp(-ret)), + std::abs(std::log((ret - old_state + uol) / tau) + ret)) < 1e-12); return ret; } diff --git a/src/compute_state_ruina.cc b/src/compute_state_ruina.cc index de844ce2..077650fc 100644 --- a/src/compute_state_ruina.cc +++ b/src/compute_state_ruina.cc @@ -4,10 +4,10 @@ #include "compute_state_ruina.hh" -double state_update_ruina(double h, double uol, double old_state) { +double state_update_ruina(double tau, double uol, double old_state) { if (uol == 0) return old_state; - double const ret = old_state - uol * std::log(uol / h); + double const ret = old_state - uol * std::log(uol / tau); return ret / (1 + uol); } diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 248c1d91..f0820d1f 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -222,7 +222,7 @@ int main(int argc, char *argv[]) { auto const eta = parset.get<double>("boundary.friction.eta"); auto const mu = parset.get<double>("boundary.friction.mu"); auto const timesteps = parset.get<size_t>("timesteps"); - double const h = 1.0 / timesteps; + double const tau = 1.0 / timesteps; octave_writer << "# name: A" << std::endl << "# type: matrix" << std::endl << "# rows: " << timesteps << std::endl << "# columns: 3" @@ -242,7 +242,7 @@ int main(int argc, char *argv[]) { ++first_frictional_node; for (size_t run = 1; run <= timesteps; ++run) { - double const time = h * run; + double const time = tau * run; if (parset.get<bool>("printProgress")) { std::cout << '*'; std::cout.flush(); @@ -263,7 +263,7 @@ int main(int argc, char *argv[]) { ++state_fpi) { auto myGlobalNonlinearity = assemble_nonlinearity<MatrixType, VectorType>( - parset.sub("boundary.friction"), *nodalIntegrals, alpha, h); + parset.sub("boundary.friction"), *nodalIntegrals, alpha, tau); MyConvexProblemType const myConvexProblem(stiffnessMatrix, *myGlobalNonlinearity, rhs); @@ -285,15 +285,16 @@ int main(int argc, char *argv[]) { // // the (logarithmic) steady state corresponding to the // // current velocity - // std::cout << std::log(L/unorm * h) << std::endl; + // std::cout << std::log(L/unorm * tau) << std::endl; switch (parset.get<Config::state_model>( "boundary.friction.state.model")) { case Config::Dieterich: - alpha[i] = state_update_dieterich(h, unorm / L, alpha_old[i]); + alpha[i] = + state_update_dieterich(tau, unorm / L, alpha_old[i]); break; case Config::Ruina: - alpha[i] = state_update_ruina(h, unorm / L, alpha_old[i]); + alpha[i] = state_update_ruina(tau, unorm / L, alpha_old[i]); break; } } @@ -321,7 +322,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() / h / L; + double const v = u_diff[first_frictional_node].two_norm() / tau / L; double const euler = alpha[first_frictional_node]; double direct; @@ -330,15 +331,15 @@ int main(int argc, char *argv[]) { } 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))); + direct = std::log( + 1.0 / v2 * + std::pow((v2 / v1), std::exp(-v2 * (run - 120) * tau))); } 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))); + direct = std::log( + 1.0 / v1 * + std::pow((v1 / v2), std::exp(-v1 * (run - 360) * tau))); } velocity_stepping_writer << euler << " " << direct << std::endl; } -- GitLab