diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh index fcfcc5be7abcff22d90d0063ec6d6d8a1cf4beb7..3bf6a11debb3173bbd97c5ff53e0795a2dff76ec 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 4bdd054e17438d0c5bcb89dfc4d1c0778b9018e5..227b8ff5c3af92f90a9c1c83d6bad6a163c9a079 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 {