localfriction.hh

#ifndef DUNE_TECTONIC_LOCAL_FRICTION_HH
#define DUNE_TECTONIC_LOCAL_FRICTION_HH

#include <cmath>
#include <limits>

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/shared_ptr.hh>

#include <dune/fufem/arithmetic.hh>
#include <dune/fufem/interval.hh>

#include "frictionpotential.hh"

namespace Dune {
template <int dimension> class LocalFriction {
public:
  using VectorType = FieldVector<double, dimension>;
  using MatrixType = FieldMatrix<double, dimension, dimension>;

  explicit LocalFriction(shared_ptr<FrictionPotentialWrapper> func)
      : func(func), smp(func->smallestPositivePoint()) {}

  double operator()(VectorType const &x) const {
    return func->evaluate(x.two_norm());
  }

  void updateState(double state) {
    func->updateState(state);
    smp = func->smallestPositivePoint();
  }

  double regularity(VectorType const &x) const {
    double const xnorm = x.two_norm();
    if (xnorm < 1e-14 && xnorm >= smp)
      return std::numeric_limits<double>::infinity();

    return func->regularity(xnorm);
  }

  // directional subdifferential: at u on the line u + t*v
  // u and v are assumed to be non-zero
  void directionalSubDiff(VectorType const &x, VectorType const &v,
                          Interval<double> &D) const {
    double const xnorm = x.two_norm();
    if (xnorm < smp)
      D[0] = D[1] = 0;
    else
      D[0] = D[1] = func->differential(xnorm) / xnorm * (x * v);
  }

  /** Formula for the derivative:

      \f{align*}{
      \frac {d^2}{dz^2} H(|z|)
      &= \frac d{dz} \left( H'(|z|) \otimes \frac z{|z|} \right)\\
      &= H''(|z|) \frac z{|z|} \otimes \frac z{|z|}
      + H'(|z|) \left( \frac {|z| \operatorname{id} - z \otimes z/|z|}{|z|^2}
     \right)\\
      &= \frac {H''(|z|)}{|z|^2} z \otimes z
      + \frac {H'(|z|)}{|z|} \operatorname{id}
      - \frac {H'(|z|)}{|z|^3} z \otimes z\\
      &= \left( \frac {H''(|z|)}{|z|^2} - \frac {H'(|z|)}{|z|^3} \right) z
     \otimes z
      + \frac {H'(|z|)}{|z|} \operatorname{id}
      \f}
  */
  void addHessian(VectorType const &x, MatrixType &A) const {
    double const xnorm2 = x.two_norm2();
    double const xnorm = std::sqrt(xnorm2);
    if (xnorm < smp)
      return;
    double const xnorm3 = xnorm * xnorm2;

    double const H1 = func->differential(xnorm);
    double const H2 = func->second_deriv(xnorm);

    // TODO: potential optimisation: factor out (H1 / xnorm), get rid of xnorm3
    double const weight1 = H2 / xnorm2;
    double const weight2 = -H1 / xnorm3;
    double const weight3 = H1 / xnorm;

    // {{{ In what follows, we handle the case 0 * (1/x) = 0 with x
    // close to 0.
    double tensorweight = 0;

    if (std::isnan(weight1))
      assert(H2 == 0);
    else
      tensorweight += weight1;

    if (std::isnan(weight2))
      assert(H1 == 0);
    else
      tensorweight += weight2;

    double idweight = 0;
    if (std::isnan(weight3))
      assert(H1 == 0);
    else
      idweight += weight3;
    // }}}

    for (int i = 0; i < dimension; ++i)
      for (int j = 0; j < i; ++j) {
        double const entry = tensorweight * x[i] * x[j];
        A[i][j] += entry;
        A[j][i] += entry;
      }

    for (int k = 0; k < dimension; ++k) {
      double const entry = tensorweight * x[k] * x[k];
      A[k][k] += entry + idweight;
    }
  }

  void addGradient(VectorType const &x, VectorType &y) const {
    double const xnorm = x.two_norm();
    if (std::isinf(func->regularity(xnorm)))
      return;

    if (xnorm >= smp)
      Arithmetic::addProduct(y, func->differential(xnorm) / xnorm, x);
  }

  void directionalDomain(VectorType const &, VectorType const &,
                         Interval<double> &dom) const {
    dom[0] = -std::numeric_limits<double>::max();
    dom[1] = std::numeric_limits<double>::max();
  }

private:
  shared_ptr<FrictionPotentialWrapper> const func;
  double smp;
};
}
#endif