evaluation of friction energy

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);

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 {