From 91b858f7dc306a27bc910b128a32fd59a365c211 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 11 Mar 2012 14:16:57 +0100
Subject: [PATCH] Externalise computation of state
---
src/Makefile.am | 2 +-
src/compute_state.cc | 45 ++++++++++++++++++++++++++++++++++++++++++
src/compute_state.hh | 9 +++++++++
src/one-body-sample.cc | 38 +----------------------------------
4 files changed, 56 insertions(+), 38 deletions(-)
create mode 100644 src/compute_state.cc
create mode 100644 src/compute_state.hh
diff --git a/src/Makefile.am b/src/Makefile.am
index cb534a6c..85f5ef83 100644
--- a/src/Makefile.am
+++ b/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)\"
diff --git a/src/compute_state.cc b/src/compute_state.cc
new file mode 100644
index 00000000..15c13a11
--- /dev/null
+++ b/src/compute_state.cc
@@ -0,0 +1,45 @@
+#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;
+}
diff --git a/src/compute_state.hh b/src/compute_state.hh
new file mode 100644
index 00000000..4ed5680e
--- /dev/null
+++ b/src/compute_state.hh
@@ -0,0 +1,9 @@
+#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
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 2c8be4f7..124786b2 100644
--- a/src/one-body-sample.cc
+++ b/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>>
--
GitLab