From 17bdd8cda8a8e79cfa83c995cc7a65fa80b6182c Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 15 Mar 2012 00:28:44 +0100
Subject: [PATCH] Rescale functional; Stricter assertions

---
 src/compute_state.cc | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/compute_state.cc b/src/compute_state.cc
index 24e420a8..599c0de4 100644
--- a/src/compute_state.cc
+++ b/src/compute_state.cc
@@ -15,9 +15,11 @@ class DecayingExponential {
   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] * (-std::exp(-u[0]));
+    D[0] = D[1] = v[0] * (-h * std::exp(-u[0]));
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
@@ -25,12 +27,15 @@ class DecayingExponential {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
+
+private:
+  double const h;
 };
 
 double compute_state_update_bisection(double h, double unorm, double L,
                                       double old_state) {
   MyDirectionalConvexFunction<DecayingExponential> const J(
-      1.0 / h, (old_state - unorm / L) / h, DecayingExponential(), 0, 1);
+      1.0, old_state - unorm / L, DecayingExponential(h), 0, 1);
   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
@@ -39,7 +44,7 @@ double compute_state_update_bisection(double h, double unorm, double L,
 double compute_state_update(double h, double unorm, double L,
                             double old_state) {
   double ret1 = compute_state_update_bisection(h, unorm, L, old_state);
-  assert(std::abs(ret1 - old_state + unorm / L - h * std::exp(-ret1)) < 1e-10);
+  assert(std::abs(ret1 - old_state + unorm / L - h * std::exp(-ret1)) < 1e-12);
 
   return ret1;
 }
-- 
GitLab