frictionpotential.hh

#ifndef DUNE_TECTONIC_FRICTIONPOTENTIAL_HH
#define DUNE_TECTONIC_FRICTIONPOTENTIAL_HH

#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>

#include <dune/common/exceptions.hh>
#include <dune/common/function.hh>

#include "frictiondata.hh"

class FrictionPotential {
public:
  virtual ~FrictionPotential() {}

  double virtual differential(double s) const = 0;
  double virtual second_deriv(double x) const = 0;
  double virtual regularity(double s) const = 0;
  double virtual coefficientOfFriction(double s) const = 0;

  double virtual evaluate(double x) const {
    DUNE_THROW(Dune::NotImplemented, "evaluation not implemented");
  }

  void virtual updateAlpha(double) = 0;
};

class TruncatedRateState : public FrictionPotential {
public:
  TruncatedRateState(double _weight, double _weightedNormalStress,
                     FrictionData _fd)
      : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}

  double coefficientOfFriction(double V) const override {
    if (V <= Vmin)
      return 0.0;

    //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;

    auto res = fd.a * std::log(V / Vmin);
    if (std::isinf(res)) {
        //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
    }

    return res;
  }

  double differential(double V) const override {
    //std::cout << "differential: " <<  weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
    //if (V <= Vmin or regularity(V)>10e8)
    //  return 0.0;

    return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
  }

  double second_deriv(double V) const override {
    //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;

    if (V <= Vmin)
      return 0.0;

    return - weightedNormalStress * (fd.a / V);
  }

  double regularity(double V) const override {
    //std::cout << "regularity: " << ((std::abs(V - Vmin) < 1e-14) ? std::numeric_limits<double>::infinity() : std::abs(second_deriv(V))) << std::endl;

    if (std::abs(V - Vmin) < 1e-14) // TODO
      return std::numeric_limits<double>::infinity();

    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);
    //Vmin = 0.0;
  }

private:
  FrictionData const fd;
  double const weight;
  double const weightedNormalStress;
  double Vmin;
};

class RegularisedRateState : public FrictionPotential {
public:
  RegularisedRateState(double _weight, double _weightedNormalStress,
                       FrictionData _fd)
      : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}

  double coefficientOfFriction(double V) const override {
    return fd.a * std::asinh(0.5 * V / Vmin);
  }

  double differential(double V) const override {
    return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
  }

  double second_deriv(double V) const override {
    return -weightedNormalStress * fd.a / std::hypot(2.0 * Vmin, V);
  }

  double regularity(double V) const override {
    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);
  }

private:
  FrictionData const fd;
  double const weight;
  double const weightedNormalStress;
  double Vmin;
};

class Tresca : public FrictionPotential {
public:
  Tresca(double _weight, double _weightedNormalStress, FrictionData _fd)
      : fd(_fd), weightedNormalStress(_weightedNormalStress) {}

  double coefficientOfFriction(double V) const override {
    return fd.mu0;
  }

  double differential(double V) const override {
    //if (std::isinf(regularity(V)))
    //    return 0.0;
    return - weightedNormalStress * coefficientOfFriction(V);
  }

  double second_deriv(double V) const override {
    return 0;
  }

  double regularity(double V) const override {
    if (std::abs(V) < 1e-7) // TODO
        return std::numeric_limits<double>::infinity();

    return std::abs(second_deriv(V));
  }

  double evaluate(double V) const override {
      return - weightedNormalStress * fd.mu0 * V;
  }

  void updateAlpha(double alpha) override {}

private:
  FrictionData const fd;
  double const weightedNormalStress;
};

class ZeroFunction : public FrictionPotential {
public:
  template <typename... Args>
  ZeroFunction(Args... args) {}

  double evaluate(double) const override { return 0; }

  double coefficientOfFriction(double s) const override { return 0; }

  double differential(double) const override { return 0; }

  double second_deriv(double) const override { return 0; }

  double regularity(double) const override { return 0; }

  void updateAlpha(double) override {}
};
#endif