#include <cassert>

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

#include <dune/tnnmg/problem-classes/bisection.hh>

#include <dune/tectonic/mydirectionalconvexfunction.hh>

#include "compute_state.hh"

// The nonlinearity exp(-x) [ Note: the Dieterich law has an additional linear
// term! ]
class DecayingExponential {
public:
  typedef Dune::FieldVector<double, 1> VectorType;
  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;

  DecayingExponential(double h) : h(h) {}

  void directionalSubDiff(VectorType const &u, VectorType const &v,
                          Interval<double> &D) const {
    D[0] = D[1] = v[0] * (-h * std::exp(-u[0]));
  }

  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:
  double const h;
};

double state_update_dieterich_bisection(double h, double uol,
                                        double old_state) {
  DecayingExponential::VectorType const start(0);
  DecayingExponential::VectorType const direction(1);
  DecayingExponential const phi(h);
  MyDirectionalConvexFunction<DecayingExponential> const J(
      1.0, old_state - uol, phi, start, direction);

  int bisectionsteps = 0;
  Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO
  return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
}

double state_update_dieterich(double h, double uol, // unorm / L
                              double old_state) {
  double const ret = state_update_dieterich_bisection(h, uol, old_state);
  /* We have

     ret - old_state + uol = h*exp(-ret)

     or

     log((ret - old_state + uol)/h) = -ret
  */
  assert(std::abs(std::log((ret - old_state + uol) / h) + ret) < 1e-13);
  return ret;
}