Do not use bisection for the Ruina update

src/compute_state_ruina.cc

 #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 <cmath>
 
 #include "compute_state_ruina.hh"
 
-// Ruina's nonlinearity
-class Ruina {
-public:
-  typedef Dune::FieldVector<double, 1> VectorType;
-  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
-
-  Ruina(double h, double uol) : h(h), uol(uol) {}
-
-  void directionalSubDiff(VectorType const &u, VectorType const &v,
-                          Interval<double> &D) const {
-    if (uol == 0)
-      D[0] = D[1] = v[0] * uol * (u[0]);
-    else
-      D[0] = D[1] = v[0] * uol * (u[0] + std::log(uol / h));
-  }
-
-  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 const uol;
-};
-
-double state_update_ruina_bisection(double h, double uol, double old_state) {
-  Ruina::VectorType const start(0);
-  Ruina::VectorType const direction(1);
-  Ruina const phi(h, uol);
-  MyDirectionalConvexFunction<Ruina> const J(1.0, old_state, 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_ruina(double h, double uol, // unorm/L
                           double old_state) {
-  double const uol_log_uol = (uol == 0) ? 0.0 : uol * std::log(uol);
-
-  double ret1 = state_update_ruina_bisection(h, uol, old_state);
-  assert(std::abs(ret1 + uol * ret1 - old_state + uol_log_uol -
-                  uol * std::log(h)) < 1e-14);
+  if (uol == 0)
+    return old_state;
 
-  return ret1;
+  double const ret = old_state - uol * std::log(uol / h);
+  return ret / (1 + uol);
 }