From 7f3161902689a7a6734522252484c5a6ed58c4d8 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Tue, 16 Jul 2019 17:54:09 +0200
Subject: [PATCH] evaluation of friction energy

---
 dune/tectonic/frictionpotential.hh | 11 +++++++++++
 dune/tectonic/localfriction.hh     |  5 +++++
 2 files changed, 16 insertions(+)

diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index fcfcc5be..3bf6a11d 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -58,6 +58,13 @@ class TruncatedRateState : public FrictionPotential {
     return std::abs(second_deriv(V));
   }
 
+  double evaluate(double V) const override {
+      if (V <= Vmin)
+        return 0.0;
+
+      return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
+  }
+
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
@@ -92,6 +99,10 @@ class RegularisedRateState : public FrictionPotential {
     return std::abs(second_deriv(V));
   }
 
+  double evaluate(double V) const override {
+      return weight * fd.C * V - weightedNormalStress * 4 * fd.a * Vmin * std::pow(std::asinh(0.25 * V / Vmin), 2.0);
+  }
+
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 4bdd054e..227b8ff5 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -19,6 +19,7 @@ template <size_t dimension> class LocalFriction {
   using VectorType = Dune::FieldVector<double, dimension>;
   using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
 
+  double virtual operator()(VectorType const &x) const = 0;
   void virtual updateAlpha(double alpha) = 0;
   double virtual regularity(VectorType const &x) const = 0;
   double virtual coefficientOfFriction(VectorType const &x) const = 0;
@@ -44,6 +45,10 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
   WrappedScalarFriction(Args... args)
       : func_(args...) {}
 
+  double operator()(VectorType const &x) const override {
+    return func_.evaluate(x.two_norm());
+  }
+
   void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
 
   double regularity(VectorType const &x) const override {
-- 
GitLab