#ifndef NICE_FUNCTION_HH #define NICE_FUNCTION_HH #include <algorithm> #include <cassert> #include <cmath> #include <limits> #include <dune/common/exceptions.hh> #include <dune/common/function.hh> #include "frictiondata.hh" class FrictionPotentialWrapper { public: virtual ~FrictionPotentialWrapper() {} double virtual differential(double s) const = 0; double virtual second_deriv(double x) const = 0; double virtual regularity(double s) const = 0; double virtual evaluate(double x) const { DUNE_THROW(Dune::NotImplemented, "evaluation not implemented"); } void virtual updateLogState(double) = 0; }; class FrictionPotential : public FrictionPotentialWrapper { public: FrictionPotential(double coefficient, double _normalStress, FrictionData const &_fd) : fd(_fd), weight(coefficient), normalStress(_normalStress) {} double differential(double V) const { assert(V >= 0.0); if (V <= V_cutoff) return 0.0; return weight * (-normalStress) * fd.a * (std::log(V) - logV_m); } double second_deriv(double V) const { assert(V >= 0); if (V <= V_cutoff) return 0; return weight * (-normalStress) * (fd.a / V); } double regularity(double V) const { assert(V >= 0); if (std::abs(V - V_cutoff) < 1e-14) // TODO return std::numeric_limits<double>::infinity(); return std::abs(second_deriv(V)); } void updateLogState(double logState) { double const tmp = (fd.mu0 + fd.b * (logState + std::log(fd.V0 / fd.L))) / fd.a; logV_m = std::log(fd.V0) - tmp; V_cutoff = fd.V0 / std::exp(tmp); } private: FrictionData const &fd; double const weight; double const normalStress; double logV_m; double V_cutoff; // velocity at which mu falls short of mumin }; class TrivialFunction : public FrictionPotentialWrapper { public: double evaluate(double) const { return 0; } double differential(double) const { return 0; } double second_deriv(double) const { return 0; } double regularity(double) const { return 0; } void updateLogState(double) {} }; #endif