Do not split the Dieterich nonlinearity

It makes the code harder to read and the computational benefit is
questionable

src/compute_state_dieterich.cc

@@ -10,18 +10,17 @@
 
 #include "compute_state_dieterich.hh"
 
-// The nonlinearity exp(-x)
-// NOTE: the Dieterich law has an additional linear term!
-class DecayingExponential {
+// psi(beta) = V/L beta + e^(-beta)
+class DieterichNonlinearity {
 public:
   typedef Dune::FieldVector<double, 1> VectorType;
   typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
 
-  DecayingExponential(double tau) : tau(tau) {}
+  DieterichNonlinearity(double _VoL) : VoL(_VoL) {}
 
   void directionalSubDiff(VectorType const &u, VectorType const &v,
                           Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-tau * std::exp(-u[0]));
+    D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
@@ -31,16 +30,16 @@ class DecayingExponential {
   }
 
 private:
-  double const tau;
+  double const VoL;
 };
 
 double state_update_dieterich_bisection(double tau, double VoL,
                                         double old_state) {
-  DecayingExponential::VectorType const start(0);
-  DecayingExponential::VectorType const direction(1);
-  DecayingExponential const phi(tau);
-  MyDirectionalConvexFunction<DecayingExponential> const J(
-      1.0, old_state - tau * VoL, phi, start, direction);
+  DieterichNonlinearity::VectorType const start(0);
+  DieterichNonlinearity::VectorType const direction(1);
+  DieterichNonlinearity const phi(VoL);
+  MyDirectionalConvexFunction<DieterichNonlinearity> const J(
+      1.0 / tau, old_state / tau, phi, start, direction);
 
   int bisectionsteps = 0;
   Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO