nicefunction.hh

#ifndef NICE_FUNCTION_HH
#define NICE_FUNCTION_HH

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

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

namespace Dune {
class NiceFunction {
protected:
  NiceFunction(std::vector<double> const &_kinks) : kinks(_kinks) {}

private:
  std::vector<double> const kinks;

public:
  std::vector<double> const &get_kinks() const { return kinks; }

  NiceFunction() : kinks() {}

  virtual ~NiceFunction() {}

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

  double virtual second_deriv(double x) const = 0;

  double virtual regularity(double s) const = 0;

  // Whether H(|.|) is smooth at zero
  bool virtual smoothesNorm() const { return false; }

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

class RuinaFunction : public NiceFunction {
public:
  RuinaFunction(double coefficient, double a, double mu, double V0,
                double normalStress, double b, double state, double L)
      : NiceFunction(),
        a(a),
        V0(V0),
        coefficientProduct(coefficient * normalStress),
        K(mu +
          b * (state + std::log(V0 / L))), // state is assumed to be logarithmic
        rho(std::exp(-K / a)) {}

  double virtual evaluate(double x) const {
    double const arg = x / V0;
    if (arg <= rho)
      return 0;

    return coefficientProduct *
           (+a * arg * (std::log(arg) - 1) + K * arg + a * rho);
  }

  double virtual leftDifferential(double x) const {
    double const arg = x / V0;
    if (arg <= rho)
      return 0;

    double const ret = a * std::log(arg) + K;
    return coefficientProduct * ret;
  }

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

  double virtual second_deriv(double x) const {
    assert(x >= 0);
    assert(V0 > 0);
    double const arg = x / V0;
    if (arg <= rho)
      return 0;

    return coefficientProduct * a / x;
  }

  double virtual regularity(double x) const {
    assert(x >= 0);
    assert(V0 > 0);
    double const arg = x / V0;
    // TODO: Make this controllable
    if (std::abs(arg - rho) < 1e-14)
      return std::numeric_limits<double>::infinity();

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

  bool virtual smoothesNorm() const { return true; }

private:
  double const a;
  double const V0;
  double const coefficientProduct;
  double const K;
  double const rho;
};

class TrivialFunction : public NiceFunction {
public:
  TrivialFunction() : NiceFunction() {}

  double virtual evaluate(double) const { return 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; }

  bool virtual smoothesNorm() const { return true; }
};
}
#endif