Externalise computation of state

src/Makefile.am

@@ -23,7 +23,7 @@ run-one-body-sample-gdb: one-body-sample
 	libtool --mode execute gdb ./one-body-sample
 
 one_body_sample_SOURCES = \
-	one-body-sample.cc assemblers.cc LambertW.cc
+	one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc
 one_body_sample_CPPFLAGS = \
 	$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\"
 

src/compute_state.cc

+#include "LambertW.h"
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/fufem/interval.hh>
+
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "compute_state.hh"
+
+// The nonlinearity exp(-x)
+class DecayingExponential {
+public:
+  typedef Dune::FieldVector<double, 1> VectorType;
+  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
+
+  double operator()(VectorType const &u) const { return exp(-u[0]); }
+
+  void directionalSubDiff(VectorType const &u, VectorType const &v,
+                          Interval<double> &D) const {
+    D[0] = D[1] = v[0] * (-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();
+  }
+};
+
+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);
+  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 compute_state_update_lambert(double h, double unorm, double L,
+                                    double old_state) {
+  double const rhs = unorm / L - old_state;
+  return LambertW(0, h * exp(rhs)) - rhs;
+}

src/compute_state.hh

+#ifndef COMPUTE_STATE_HH
+#define COMPUTE_STATE_HH
+
+double compute_state_update_bisection(double h, double unorm, double L,
+                                      double old_state);
+double compute_state_update_lambert(double h, double unorm, double L,
+                                    double old_state);
+
+#endif

src/one-body-sample.cc

@@ -23,8 +23,6 @@
 
 #include <cmath>
 
-#include "LambertW.h"
-
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/exceptions.hh>
 #include <dune/common/fmatrix.hh>
@@ -67,6 +65,7 @@
 #include <dune/tectonic/mydirectionalconvexfunction.hh>
 
 #include "assemblers.hh"
+#include "compute_state.hh"
 
 int const dim = 2;
 
@@ -108,41 +107,6 @@ void setup_boundary(GridView const &gridView,
   std::cout << "Number of frictional nodes: " << dirichlet_nodes << std::endl;
 }
 
-// The nonlinearity exp(-x)
-class DecayingExponential {
-public:
-  typedef Dune::FieldVector<double, 1> VectorType;
-  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
-
-  double operator()(VectorType const &u) const { return exp(-u[0]); }
-
-  void directionalSubDiff(VectorType const &u, VectorType const &v,
-                          Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-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();
-  }
-};
-
-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);
-  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 compute_state_update_lambert(double h, double unorm, double L,
-                                    double old_state) {
-  double const rhs = unorm / L - old_state;
-  return LambertW(0, h * exp(rhs)) - rhs;
-}
-
 int main(int argc, char *argv[]) {
   try {
     typedef SharedPointerMap<std::string, Dune::VirtualFunction<double, double>>