nicefunction.hh

#ifndef NICE_FUNCTION_HH
#define NICE_FUNCTION_HH

#include <algorithm>
#include <cmath>

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

namespace Dune {
class NiceFunction : public VirtualFunction<double, double> {
public:
  virtual ~NiceFunction() {}

  double virtual leftDifferential(double s) const = 0;
  double virtual rightDifferential(double s) const = 0;

  double virtual second_deriv(double x) const {
    DUNE_THROW(NotImplemented, "second derivative not implemented");
  }

  double virtual regularity(double s) const {
    DUNE_THROW(NotImplemented, "regularity not implemented");
  }
};

class RuinaFunction : public NiceFunction {
public:
  RuinaFunction(double coefficient, double a, double mu, double eta,
                double normalStress)
      : coefficient(coefficient),
        a(a),
        mu(mu),
        rho(exp(-mu / a)),
        eta(eta),
        normalStress(normalStress) {}

  /*
    If mu and sigma_n denote the coefficient of friction and the
    normal stress, respectively, this function is given by

    1/eta sigma_n h(max(eta id, rho)) + a rho + mu

    with the constants a and b from Ruina's model and

    h(beta) = beta (a (log beta - 1) + mu)

    as well as

    rho = exp(-mu/a)
  */
  void virtual evaluate(double const &x, double &y) const {
    double const arg = std::max(eta * x, rho);
    double const h = arg * (a * (std::log(arg) - 1) + mu);
    y = 1 / eta * normalStress * h + a * rho + mu;
    y *= coefficient;
  }

  /*
    (leaving some terms aside): with s > rho

    1/eta d/dx [ a * (s log s - s) + mu s ] where s = eta x
    = 1/eta [ a * (log (eta x) * eta) + eta mu ]
    = a * log(eta x) + mu
  */
  double virtual leftDifferential(double s) const {
    if (eta * s <= rho)
      return 0;

    return coefficient * normalStress * (a * std::log(eta * s) + mu);
  }

  /* see above */
  double virtual rightDifferential(double s) const {
    if (eta * s <= rho)
      return 0;

    return coefficient * normalStress * (a * std::log(eta * s) + mu);
  }

  /*
    d/dx a * log(eta x) + mu
    = a * 1/(eta x) * eta
    = a/x
  */
  double virtual second_deriv(double s) const {
    // includes the case eta * s = rho for which there is no second derivative
    if (eta * s <= rho)
      return 0;
    else
      return coefficient * normalStress * (a / s);
  }

  double virtual regularity(double s) const {
    if (eta * s == rho)
      return std::numeric_limits<double>::infinity();

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

private:
  double coefficient;
  double a;
  double mu;
  double eta;
  double normalStress;

  double rho;
};

class LinearFunction : public NiceFunction {
public:
  LinearFunction(double a) : coefficient(a) {}

  void virtual evaluate(double const &x, double &y) const {
    y = coefficient * x;
  }

  double virtual leftDifferential(double s) const { return coefficient; }

  double virtual rightDifferential(double s) const { return coefficient; }

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

  double virtual regularity(double s) const { return 0; }

private:
  double coefficient;
};

template <int slope> class SampleFunction : public NiceFunction {
public:
  void virtual evaluate(double const &x, double &y) const {
    y = (x < 1) ? x : (slope * (x - 1) + 1);
  }

  double virtual leftDifferential(double s) const {
    return (s <= 1) ? 1 : slope;
  }

  double virtual rightDifferential(double s) const {
    return (s < 1) ? 1 : slope;
  }
};

class SteepFunction : public NiceFunction {
public:
  void virtual evaluate(double const &x, double &y) const { y = 100 * x; }

  double virtual leftDifferential(double s) const { return 100; }

  double virtual rightDifferential(double s) const { return 100; }
};

class TrivialFunction : public NiceFunction {
public:
  void virtual evaluate(double const &x, double &y) const { y = 0; }

  double virtual leftDifferential(double) const { return 0; }

  double virtual rightDifferential(double) const { return 0; }

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

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

// slope in [n-1,n] is n
class HorribleFunction : public NiceFunction {
public:
  void virtual evaluate(double const &x, double &y) const {
    double const fl = floor(x);
    double const sum = fl * (fl + 1) / 2;
    y = sum + (fl + 1) * (x - fl);
  }

  double virtual leftDifferential(double x) const {
    double const fl = floor(x);
    if (x - fl < 1e-14)
      return fl;
    else
      return fl + 1;
  }

  double virtual rightDifferential(double x) const {
    double const c = ceil(x);
    if (c - x < 1e-14)
      return c + 1;
    else
      return c;
  }
};

// slope in [n-1,n] is log(n+1)
class HorribleFunctionLogarithmic : public NiceFunction {
public:
  void virtual evaluate(double const &x, double &y) const {
    y = 0;
    size_t const fl = floor(x);
    for (size_t i = 1; i <= fl;)
      y += log(++i); // factorials grow to fast so we compute this incrementally

    y += log(fl + 2) * (x - fl);
  }

  double virtual leftDifferential(double x) const {
    double const fl = floor(x);
    if (x - fl < 1e-14)
      return log(fl + 1);
    else
      return log(fl + 2);
  }

  double virtual rightDifferential(double x) const {
    double const c = ceil(x);
    if (c - x < 1e-14)
      return log(c + 2);
    else
      return log(c + 1);
  }
};
}
#endif