diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh index aa65fb9ebb5d8c4432c04db2fd6d1b526eb1336b..e284e95a83dc3286987d0410b18aa314f155bd1e 100644 --- a/dune/tectonic/nicefunction.hh +++ b/dune/tectonic/nicefunction.hh @@ -34,43 +34,40 @@ class RuinaFunction : public NiceFunction { public: RuinaFunction(double coefficient, double a, double mu, double eta, double normalStress, double b, double state, double L, double h) - : coefficient(coefficient), - a(a), - mu(mu), - eta(eta), - normalStress(normalStress), - compound_state(b * (state - std::log(eta * L))), - h(h), - rho(exp(-(mu + compound_state) / a)) {} - - // TODO: untested + : a(a), + eta_h(eta / h), + coefficientProduct(coefficient * normalStress), + c(mu + (a - b) * std::log(eta) - b * std::log(L)), + K(b * state) // state is assumed to be logarithmic + {} + void virtual evaluate(double const &x, double &y) const { - assert(false); - double const arg = eta * x / h; - double const expstar = (arg == 0) ? 0 : arg * std::log(arg) - arg; - double const gamma = - normalStress * (mu + a * std::log(eta) + compound_state); - y = a * normalStress * expstar + gamma * arg; - y *= coefficient * h * normalStress / eta; + double const arg = eta_h * x; + if (arg == 0) { // TODO: Make this controllable + y = 0; + return; + } + + double const expstar = arg * std::log(arg) - arg; + y = a * expstar + (c + K) * arg; + y *= coefficientProduct / eta_h; } - double virtual leftDifferential(double s) const { - double const arg = eta * s / h; - if (arg == 0) + double virtual leftDifferential(double x) const { + double const arg = eta_h * x; + if (arg == 0) // TODO: Make this controllable return 0; - double const gamma = - normalStress * (mu + a * std::log(eta) + compound_state); - return coefficient * (a * normalStress * std::log(arg) + gamma); + double const ret = (a * std::log(arg) + c + K); + return coefficientProduct * ret; } - /* see above */ double virtual rightDifferential(double s) const { return leftDifferential(s); } double virtual second_deriv(double s) const { - return coefficient * a * normalStress / s; + return coefficientProduct * a / s; } double virtual regularity(double s) const { @@ -81,15 +78,11 @@ class RuinaFunction : public NiceFunction { } private: - double const coefficient; double const a; - double const mu; - double const eta; - double const normalStress; - double const compound_state; - double const h; - - double const rho; + double const eta_h; + double const coefficientProduct; + double const c; + double const K; }; class LinearFunction : public NiceFunction {