From 3a625ee19bf869e4d165af6ec7e4551b188db272 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Sat, 6 Mar 2021 13:30:35 +0100
Subject: [PATCH] fix Tresca friction
---
.../friction/frictionpotential.hh | 27 ++++++++++---------
1 file changed, 15 insertions(+), 12 deletions(-)
diff --git a/dune/tectonic/data-structures/friction/frictionpotential.hh b/dune/tectonic/data-structures/friction/frictionpotential.hh
index becfb81d..81b47741 100644
--- a/dune/tectonic/data-structures/friction/frictionpotential.hh
+++ b/dune/tectonic/data-structures/friction/frictionpotential.hh
@@ -34,14 +34,14 @@ class TruncatedRateState : public FrictionPotential {
: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
double coefficientOfFriction(double V) const override {
- if (V <= Vmin or Vmin == 0.0)
+ if (V <= Vmin or regularity(V)>10e8)
return 0.0;
//std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
auto res = fd.a * std::log(V / Vmin);
if (std::isinf(res)) {
- std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
+ //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
}
return res;
@@ -49,16 +49,19 @@ class TruncatedRateState : public FrictionPotential {
double differential(double V) const override {
//std::cout << "differential: " << weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
+ if (V <= Vmin or regularity(V)>10e8)
+ return 0.0;
+
return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
}
double second_deriv(double V) const override {
//std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;
- if (V <= Vmin or Vmin == 0.0)
+ if (V <= Vmin)
return 0.0;
- return -weightedNormalStress * (fd.a / V);
+ return - weightedNormalStress * (fd.a / V);
}
double regularity(double V) const override {
@@ -71,10 +74,7 @@ class TruncatedRateState : public FrictionPotential {
}
double evaluate(double V) const override {
-
- //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1)) << std::endl;
-
- if (V <= Vmin or Vmin == 0.0)
+ if (V <= Vmin or regularity(V)>10e8)
return 0.0;
return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
@@ -83,9 +83,7 @@ class TruncatedRateState : public FrictionPotential {
void updateAlpha(double alpha) override {
double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
Vmin = fd.V0 / std::exp(logrest);
-
- if (Vmin == 0.0)
- std::cout << " logrest: " << logrest << " alpha: " << alpha << std::endl;
+ //Vmin = 0.0;
}
private:
@@ -143,6 +141,8 @@ class Tresca : public FrictionPotential {
}
double differential(double V) const override {
+ //if (std::isinf(regularity(V)))
+ // return 0.0;
return - weightedNormalStress * coefficientOfFriction(V);
}
@@ -151,11 +151,14 @@ class Tresca : public FrictionPotential {
}
double regularity(double V) const override {
+ if (std::abs(V) < 1e-14) // TODO
+ return std::numeric_limits<double>::infinity();
+
return std::abs(second_deriv(V));
}
double evaluate(double V) const override {
- return 0;
+ return - weightedNormalStress * fd.mu0 * V;
}
void updateAlpha(double alpha) override {}
--
GitLab