From b75bb0903d8b4c0840d7c29748c3f143e26a7816 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 11 Mar 2012 11:57:35 +0100
Subject: [PATCH] Fix up RuinaFunction
---
dune/tectonic/nicefunction.hh | 52 ++++++++---------------------------
1 file changed, 12 insertions(+), 40 deletions(-)
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index aea4dbd7..a779b440 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -41,44 +41,25 @@ class RuinaFunction : public NiceFunction {
h(h),
rho(exp(-(mu + compound_state) / a)) {}
- /*
- If mu and sigma_n denote the coefficient of friction and the
- normal stress, respectively, this function is given by
-
- 1/eta sigma_n r(max(eta id, rho)) + a rho + mu
-
- with the constants a and b from Ruina's model and
-
- r(beta) = beta (a (log beta - 1) + mu)
-
- as well as
-
- rho = exp(-mu/a)
- */
+ // TODO: untested
void virtual evaluate(double const &x, double &y) const {
+ assert(false);
double const arg = eta * x / h;
- if (arg <= rho) {
- y = -std::log(rho); // TODO: We can write this out explicitly
- return;
- }
- y = arg * (std::log(arg) - 1) - arg * std::log(rho) + rho - std::log(rho);
+ 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;
}
- /*
- (leaving some terms aside): with s > rho
-
- 1/eta d/dx [ a * (s log s - s) + mu s ] where s = eta x
- = 1/eta [ a * (log (eta x) * eta) + eta mu ]
- = a * log(eta x) + mu
- */
double virtual leftDifferential(double s) const {
double const arg = eta * s / h;
- if (arg <= rho)
+ if (arg == 0)
return 0;
- return coefficient * normalStress *
- (a * std::log(arg) + mu + compound_state);
+ double const gamma =
+ normalStress * (mu + a * std::log(eta) + compound_state);
+ return coefficient * (a * normalStress * std::log(arg) + gamma);
}
/* see above */
@@ -86,21 +67,12 @@ class RuinaFunction : public NiceFunction {
return leftDifferential(s);
}
- /*
- d/dx a * log(eta x) + mu
- = a * 1/(eta x) * eta
- = a/x
- */
double virtual second_deriv(double s) const {
- // includes the case eta * s = rho for which there is no second derivative
- if (eta * s / h <= rho)
- return 0;
- else
- return coefficient * normalStress * a / s;
+ return coefficient * a * normalStress / s;
}
double virtual regularity(double s) const {
- if (eta * s == rho)
+ if (s == 0)
return std::numeric_limits<double>::infinity();
return std::abs(second_deriv(s));
--
GitLab