diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh index eb9388030a0965439d604b3422ce77b51e084419..f76235862b2d0d35278b6c76099123474b99141f 100644 --- a/dune/tectonic/nicefunction.hh +++ b/dune/tectonic/nicefunction.hh @@ -38,56 +38,62 @@ class NiceFunction { } }; +// phi(V) = V log(V/V_m) - V + V_m if V >= V_m +// = 0 otherwise +// with V_m = V_0 exp(-K/a), +// i.e. K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha class RuinaFunction : public NiceFunction { public: RuinaFunction(double coefficient, double a, double mu_0, double V_0, double normalStress, double b, double state, double L) : NiceFunction(), - a(a), - V_0(V_0), - coefficientProduct(coefficient * normalStress), + coefficientProduct(coefficient * a * normalStress), // state is assumed to be logarithmic - K(mu_0 + b * (state + std::log(V_0 / L))), - V_m(std::exp(-K / a)) {} + V_m(V_0 * std::exp(-(mu_0 + b * (state + std::log(V_0 / L))) / a)) + // We could also compute V_m as + // V_0 * std::exp(-(mu_0 + b * state)/a) + // * std::pow(V_0 / L, -b/a) + // which would avoid the std::exp(std::log()) + {} double virtual evaluate(double V) const { - double const arg = V / V_0; - if (arg <= V_m) + assert(V >= 0); + if (V <= V_m) return 0; - return coefficientProduct * - (+a * arg * (std::log(arg) - 1) + K * arg + a * V_m); + // V log(V/V_m) - V + V_m + return coefficientProduct * (V * std::log(V / V_m) - V + V_m); } + // log(V/V_m) if V >= V_0 + // 0 otherwise double virtual leftDifferential(double V) const { - double const arg = V / V_0; - if (arg <= V_m) + assert(V >= 0); + if (V <= V_m) return 0; - double const ret = a * std::log(arg) + K; - return coefficientProduct * ret; + return coefficientProduct * std::log(V / V_m); } double virtual rightDifferential(double V) const { return leftDifferential(V); } + // 1/V if V > V_0 + // undefined if V == V_0 + // 0 if V < V_0 double virtual second_deriv(double V) const { assert(V >= 0); - assert(V_0 > 0); - double const arg = V / V_0; - if (arg <= V_m) + if (V <= V_m) return 0; - return coefficientProduct * a / V; + return coefficientProduct / V; } double virtual regularity(double V) const { assert(V >= 0); - assert(V_0 > 0); - double const arg = V / V_0; // TODO: Make this controllable - if (std::abs(arg - V_m) < 1e-14) + if (std::abs(V - V_m) < 1e-14) return std::numeric_limits<double>::infinity(); return std::abs(second_deriv(V)); @@ -96,10 +102,7 @@ class RuinaFunction : public NiceFunction { bool virtual smoothesNorm() const { return true; } private: - double const a; - double const V_0; double const coefficientProduct; - double const K; double const V_m; };