From 60e975803d017e007bc1ffe747a764e3832dae54 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sat, 20 Jul 2013 11:38:54 +0200
Subject: [PATCH] [Problem] Allow mu to be cut off (from below)
---
dune/tectonic/frictiondata.hh | 2 ++
dune/tectonic/frictionpotential.hh | 30 ++++++++++++------------------
src/one-body-sample.parset | 1 +
3 files changed, 15 insertions(+), 18 deletions(-)
diff --git a/dune/tectonic/frictiondata.hh b/dune/tectonic/frictiondata.hh
index 8cc72469..656683f9 100644
--- a/dune/tectonic/frictiondata.hh
+++ b/dune/tectonic/frictiondata.hh
@@ -9,6 +9,7 @@ struct FrictionData {
a(parset.get<double>("a")),
b(parset.get<double>("b")),
mu0(parset.get<double>("mu0")),
+ mumin(parset.get<double>("mumin")),
normalStress(normalStress) {}
double const L;
@@ -16,6 +17,7 @@ struct FrictionData {
double const a;
double const b;
double const mu0;
+ double const mumin;
double const normalStress;
};
#endif
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index bc2b2273..73d17dba 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -34,23 +34,15 @@ class FrictionPotential : public FrictionPotentialWrapper {
FrictionPotential(double coefficient, FrictionData const &fd)
: fd(fd), weightTimesNormalStress(coefficient * fd.normalStress) {}
- double evaluate(double V) const {
- assert(V >= 0);
- if (V <= V_m)
- return 0;
-
- // V log(V/V_m) - V + V_m
- return weightTimesNormalStress * fd.a * (V * std::log(V / V_m) - V + V_m);
- }
-
// log(V/V_m) if V >= V_0
// 0 otherwise
double differential(double V) const {
- assert(V >= 0);
- if (V <= V_m)
- return 0;
+ assert(V >= 0.0);
+ if (V <= V_cutoff) {
+ return fd.mumin;
+ }
- return weightTimesNormalStress * fd.a * std::log(V / V_m);
+ return weightTimesNormalStress * fd.a * (std::log(V) - logV_m);
}
// 1/V if V > V_0
@@ -58,7 +50,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
// 0 if V < V_0
double second_deriv(double V) const {
assert(V >= 0);
- if (V <= V_m)
+ if (V <= V_cutoff)
return 0;
return weightTimesNormalStress * fd.a / V;
@@ -67,7 +59,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
double regularity(double V) const {
assert(V >= 0);
// TODO: Make this controllable
- if (std::abs(V - V_m) < 1e-14)
+ if (std::abs(V - V_cutoff) < 1e-14)
return std::numeric_limits<double>::infinity();
return std::abs(second_deriv(V));
@@ -77,14 +69,16 @@ class FrictionPotential : public FrictionPotentialWrapper {
// with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ]
void updateState(double state) {
// state is assumed to be logarithmic
- V_m = fd.V0 *
- std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
+ logV_m = std::log(fd.V0) +
+ (-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
+ V_cutoff = std::exp(logV_m + fd.mumin / fd.a);
}
private:
FrictionData const &fd;
double const weightTimesNormalStress;
- double V_m;
+ double logV_m;
+ double V_cutoff; // velocity at which mu falls short of mumin
};
class TrivialFunction : public FrictionPotentialWrapper {
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 1f819b24..d3eecdd8 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -51,6 +51,7 @@ requiredResidual = 1e-12
[boundary.friction]
mu0 = 0.6
+mumin = 0.0
a = 0.010
b = 0.015
V0 = 1e-6
--
GitLab