mynonlinearity.hh

/* -*- mode:c++; mode:semantic -*- */

#ifndef MYNONLINEARITY_HH
#define MYNONLINEARITY_HH

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

#include <dune/fufem/interval.hh>

#include <limits>

#include "nicefunction.hh"

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

  MyNonlinearity(NiceFunction const &func) : func_(func) {}

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

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

  void upperGradient(VectorType const x, VectorType &ret) const {
    ret = x;
    ret *= func_.rightDifferential(x.two_norm()) / x.two_norm();
  }

  void lowerGradient(VectorType const x, VectorType &ret) const {
    ret = x;
    ret *= func_.leftDifferential(x.two_norm()) / x.two_norm();
  }

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