frictionpotential.hh

#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