localnonlinearity.hh

#ifndef DUNE_TECTONIC_LOCAL_NONLINEARITY_HH
#define DUNE_TECTONIC_LOCAL_NONLINEARITY_HH

#include <cmath>
#include <limits>

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

#include <dune/fufem/interval.hh>

#include "nicefunction.hh"

namespace Dune {
template <int dimension> class LocalNonlinearity {
public:
  typedef FieldVector<double, dimension> VectorType;
  typedef FieldMatrix<double, dimension, dimension> MatrixType;

  LocalNonlinearity(shared_ptr<NiceFunction const> func) : func_(func) {}

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

  double regularity(VectorType const &x) const {
    if (x.two_norm() < 1e-8) // TODO
      return std::numeric_limits<double>::infinity();

    return func_->regularity(x.two_norm());
  }

  // directional subdifferential: at u on the line u + t*v
  // u and v are assumed to be non-zero
  void directionalSubDiff(VectorType const &u, VectorType const &v,
                          Interval<double> &D) const {
    if (u.two_norm() == 0) {
      D[0] = D[1] = func_->rightDifferential(0) * v.two_norm();
      return;
    }
    double const un = u.two_norm();
    double const ndotp = (u * v) / un;
    // Our coordinate system is now such that v is a unit vector!
    if (ndotp > 0) {
      D[1] = ndotp * func_->rightDifferential(un);
      D[0] = ndotp * func_->leftDifferential(un);
    } else {
      D[1] = ndotp * func_->leftDifferential(un);
      D[0] = ndotp * func_->rightDifferential(un);
    }
  }

  /*
    H''(|x|) x \otimes x / |x|^2
    + H'(|x|) [|x|^2 id - x \otimes x] / |x|^3

    For i == j, this is (writing k and using the notation from below)

    h2 * x[k] * x[k] / normX2 + h1 [normX2 - x[k]^2]/normX^3
    = h2 * x[k] * x[k] / normX2 + h1ox - h1 x[k]^2/normX^3
    = (h2-h1ox) * x[k] * x[k] / normX2 + h1ox

    For i != j, this is (again, using the notation from below)

    h2 * x[i] * x[j] / normX2 - h1 * x[i] * x[j]/normX^3
    = (h2 - h1ox) * x[i] * x[j] / normX2
  */
  void addHessian(VectorType const &x, MatrixType &A) const {
    if (x == VectorType(0))
      return;

    double const normX2 = x.two_norm2();
    double const normX = sqrt(normX2);
    double const h1ox = func_->rightDifferential(normX) / normX;
    double const h2 = func_->second_deriv(normX);

    for (int i = 0; i < dimension; ++i)
      for (int j = 0; j < i; ++j) {
        double const entry = (h2 - h1ox) * x[i] * x[j] / normX2;
        A[i][j] += entry;
        A[j][i] += entry;
      }

    for (int k = 0; k < dimension; ++k)
      A[k][k] += (h2 - h1ox) * x[k] * x[k] / normX2 + h1ox;
  }

  void addGradient(VectorType const &x, VectorType &y) const {
    double const xnorm = x.two_norm();
    if (xnorm < 1e-8 or std::isinf(func_->regularity(xnorm))) // TODO
      return;

    // left and right differential coincide in this case
    y.axpy(func_->rightDifferential(xnorm) / xnorm, x);
  }

  void upperGradient(VectorType const &x, VectorType &ret) const {
    double const xnorm = x.two_norm();
    assert(xnorm != 0);

    ret = x;
    ret *= func_->rightDifferential(xnorm) / xnorm;
  }

  void lowerGradient(VectorType const &x, VectorType &ret) const {
    double const xnorm = x.two_norm();
    assert(xnorm != 0);

    ret = x;
    ret *= func_->leftDifferential(xnorm) / xnorm;
  }

  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<NiceFunction const> const func_;
};
}
#endif