From df123c4b1d08f815ea1fb63129964972095ada8c Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 21 Jan 2013 15:49:57 +0100
Subject: [PATCH] Clean up RuinaFunction
Drop members a, V_0, K
Move a into coefficientProduct
Move V_0 into V_m
Add comments that show what we're computing
---
dune/tectonic/nicefunction.hh | 49 +++++++++++++++++++----------------
1 file changed, 26 insertions(+), 23 deletions(-)
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index eb938803..f7623586 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -38,56 +38,62 @@ class NiceFunction {
}
};
+// phi(V) = V log(V/V_m) - V + V_m if V >= V_m
+// = 0 otherwise
+// with V_m = V_0 exp(-K/a),
+// i.e. K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha
class RuinaFunction : public NiceFunction {
public:
RuinaFunction(double coefficient, double a, double mu_0, double V_0,
double normalStress, double b, double state, double L)
: NiceFunction(),
- a(a),
- V_0(V_0),
- coefficientProduct(coefficient * normalStress),
+ coefficientProduct(coefficient * a * normalStress),
// state is assumed to be logarithmic
- K(mu_0 + b * (state + std::log(V_0 / L))),
- V_m(std::exp(-K / a)) {}
+ V_m(V_0 * std::exp(-(mu_0 + b * (state + std::log(V_0 / L))) / a))
+ // We could also compute V_m as
+ // V_0 * std::exp(-(mu_0 + b * state)/a)
+ // * std::pow(V_0 / L, -b/a)
+ // which would avoid the std::exp(std::log())
+ {}
double virtual evaluate(double V) const {
- double const arg = V / V_0;
- if (arg <= V_m)
+ assert(V >= 0);
+ if (V <= V_m)
return 0;
- return coefficientProduct *
- (+a * arg * (std::log(arg) - 1) + K * arg + a * V_m);
+ // V log(V/V_m) - V + V_m
+ return coefficientProduct * (V * std::log(V / V_m) - V + V_m);
}
+ // log(V/V_m) if V >= V_0
+ // 0 otherwise
double virtual leftDifferential(double V) const {
- double const arg = V / V_0;
- if (arg <= V_m)
+ assert(V >= 0);
+ if (V <= V_m)
return 0;
- double const ret = a * std::log(arg) + K;
- return coefficientProduct * ret;
+ return coefficientProduct * std::log(V / V_m);
}
double virtual rightDifferential(double V) const {
return leftDifferential(V);
}
+ // 1/V if V > V_0
+ // undefined if V == V_0
+ // 0 if V < V_0
double virtual second_deriv(double V) const {
assert(V >= 0);
- assert(V_0 > 0);
- double const arg = V / V_0;
- if (arg <= V_m)
+ if (V <= V_m)
return 0;
- return coefficientProduct * a / V;
+ return coefficientProduct / V;
}
double virtual regularity(double V) const {
assert(V >= 0);
- assert(V_0 > 0);
- double const arg = V / V_0;
// TODO: Make this controllable
- if (std::abs(arg - V_m) < 1e-14)
+ if (std::abs(V - V_m) < 1e-14)
return std::numeric_limits<double>::infinity();
return std::abs(second_deriv(V));
@@ -96,10 +102,7 @@ class RuinaFunction : public NiceFunction {
bool virtual smoothesNorm() const { return true; }
private:
- double const a;
- double const V_0;
double const coefficientProduct;
- double const K;
double const V_m;
};
--
GitLab