From ef695ed05c84593dd101caebb59db70523f1fa75 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 5 Sep 2011 17:40:13 +0200
Subject: [PATCH] Populate with sample code
---
src/Makefile.am | 41 ++++---
src/bisection-simpler-example.cc | 53 +++++++++
src/bisection-simpler-example2-gradient.cc | 125 +++++++++++++++++++++
src/bisection-simpler-example2.cc | 58 ++++++++++
4 files changed, 256 insertions(+), 21 deletions(-)
create mode 100644 src/bisection-simpler-example.cc
create mode 100644 src/bisection-simpler-example2-gradient.cc
create mode 100644 src/bisection-simpler-example2.cc
diff --git a/src/Makefile.am b/src/Makefile.am
index 963b653b..a57492b1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,33 +1,32 @@
SUBDIRS =
-noinst_PROGRAMS = dune_tectonic
-
-dune_tectonic_SOURCES = dune_tectonic.cc
-
-dune_tectonic_CPPFLAGS = $(AM_CPPFLAGS) \
- $(DUNEMPICPPFLAGS) \
- $(UG_CPPFLAGS) \
- $(AMIRAMESH_CPPFLAGS) \
- $(ALBERTA_CPPFLAGS) \
+noinst_PROGRAMS = \
+ bisection-simpler-example \
+ bisection-simpler-example2 \
+ bisection-simpler-example2-gradient
+
+bisection_simpler_example_SOURCES = bisection-simpler-example.cc
+bisection_simpler_example2_SOURCES = bisection-simpler-example2.cc
+bisection_simpler_example2_gradient_SOURCES = bisection-simpler-example2-gradient.cc
+
+check-am:
+ ./bisection-simpler-example
+ ./bisection-simpler-example2
+ ./bisection-simpler-example2-gradient
+
+AM_CPPFLAGS = \
+ $(DUNE_CPPFLAGS) \
$(ALUGRID_CPPFLAGS)
+
# The libraries have to be given in reverse order (most basic libraries
# last). Also, due to some misunderstanding, a lot of libraries include the
# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
# here as well.
-dune_tectonic_LDADD = \
+LDADD = \
$(DUNE_LDFLAGS) $(DUNE_LIBS) \
- $(ALUGRID_LDFLAGS) $(ALUGRID_LIBS) \
- $(ALBERTA_LDFLAGS) $(ALBERTA_LIBS) \
- $(AMIRAMESH_LDFLAGS) $(AMIRAMESH_LIBS) \
- $(UG_LDFLAGS) $(UG_LIBS) \
- $(DUNEMPILIBS) \
- $(LDADD)
-dune_tectonic_LDFLAGS = $(AM_LDFLAGS) \
- $(DUNEMPILDFLAGS) \
- $(UG_LDFLAGS) \
- $(AMIRAMESH_LDFLAGS) \
- $(ALBERTA_LDFLAGS) \
+ $(ALUGRID_LDFLAGS) $(ALUGRID_LIBS)
+AM_LDFLAGS = \
$(ALUGRID_LDFLAGS) \
$(DUNE_LDFLAGS)
diff --git a/src/bisection-simpler-example.cc b/src/bisection-simpler-example.cc
new file mode 100644
index 00000000..f3ca7273
--- /dev/null
+++ b/src/bisection-simpler-example.cc
@@ -0,0 +1,53 @@
+#include <dune/fufem/interval.hh>
+
+#include <dune/tnnmg/nonlinearities/smallfunctional.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include <cassert>
+#include <limits>
+
+class SampleFunctional : public SmallFunctional<2> {
+public:
+ SampleFunctional(SmallMatrix A, SmallVector b) : A_(A), b_(b) {}
+
+ double operator()(const SmallVector v) const {
+ SmallVector y;
+ A_.mv(v, y); // y = Av
+ y /= 2; // y = 1/2 Av
+ y -= b_; // y = 1/2 Av - b
+ return y * v; // <1/2 Av - b,v>
+ }
+
+ SmallVector d(const SmallVector v) const {
+ SmallVector y;
+ A_.mv(v, y); // y = Av
+ y -= b_; // y = Av - b
+ return y;
+ }
+
+ SmallMatrix d2(const SmallVector v) const { return A_; }
+
+private:
+ SmallMatrix A_;
+ SmallVector b_;
+};
+
+int main() {
+ int dim = 2;
+ SampleFunctional::SmallMatrix A;
+ A[0][0] = 3;
+ A[0][1] = 0;
+ A[1][0] = 0;
+ A[1][1] = 3;
+ SampleFunctional::SmallVector b;
+ b[0] = 1;
+ b[1] = 2;
+
+ SampleFunctional J(A, b);
+
+ assert(J(b) == 2.5);
+
+ assert(b.two_norm() == sqrt(5));
+ assert(b.two_norm2() == 5);
+ return 0;
+}
diff --git a/src/bisection-simpler-example2-gradient.cc b/src/bisection-simpler-example2-gradient.cc
new file mode 100644
index 00000000..2637a48f
--- /dev/null
+++ b/src/bisection-simpler-example2-gradient.cc
@@ -0,0 +1,125 @@
+// FIXME: is there a bette way to do this?
+//#define DUNE_MINIMAL_DEBUG_LEVEL 2
+#include <dune/common/stdstreams.hh>
+
+#include <dune/fufem/interval.hh>
+
+#include <dune/tnnmg/nonlinearities/smallfunctional.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include <cassert>
+#include <limits>
+
+template <int dimension>
+class SampleFunctional : public SmallFunctional<dimension> {
+public:
+ typedef typename SmallFunctional<dimension>::SmallMatrix SmallMatrix;
+ typedef typename SmallFunctional<dimension>::SmallVector SmallVector;
+
+ SampleFunctional(SmallMatrix A, SmallVector b) : A_(A), b_(b) {}
+
+ double operator()(const SmallVector v) const {
+ SmallVector y;
+ A_.mv(v, y); // y = Av
+ y /= 2; // y = 1/2 Av
+ y -= b_; // y = 1/2 Av - b
+ return y * v; // <1/2 Av - b,v>
+ }
+
+ SmallVector d(const SmallVector v) const {
+ SmallVector y;
+ A_.mv(v, y); // y = Av
+ y -= b_; // y = Av - b
+ return y;
+ }
+
+ SmallMatrix d2(const SmallVector v) const { return A_; }
+
+ // extending the interface here
+ double directionalDerivative(const SmallVector x,
+ const SmallVector dir) const {
+ SmallVector grad = d(x);
+ return grad * dir;
+ }
+
+ SmallVector minimise(const SmallVector x, int iterations) const {
+ SmallVector descDir = d(x);
+ descDir *= -1; // The negative gradient
+
+ Dune::dverb << "Starting at x with J(x) = " << operator()(x) << std::endl;
+ Dune::dverb << "Minimizing in direction w with dJ(x,w) = "
+ << directionalDerivative(x, descDir) << std::endl;
+
+ double l = 0;
+ double r = 1;
+ SmallVector tmp = descDir;
+ tmp *= r;
+ while (directionalDerivative(x + tmp, descDir) < 0) {
+ l = r;
+ r *= 2;
+ tmp *= 2;
+ Dune::dverb << "Widened interval!" << std::endl;
+ }
+ Dune::dverb << "Interval now [" << l << "," << r << "]" << std::endl;
+
+ // Debugging
+ {
+ SmallVector tmpl = tmp;
+ tmpl *= l;
+ SmallVector tmpr = tmp;
+ tmpr *= r;
+ assert(directionalDerivative(x + tmpl, descDir) < 0);
+ assert(directionalDerivative(x + tmpr, descDir) > 0);
+ }
+
+ double m = l / 2 + r / 2;
+ SmallVector middle;
+ for (size_t count = 0; count < iterations; ++count) {
+ Dune::dverb << "now at m = " << m << std::endl;
+ Dune::dverb << "Value of J here: " << operator()(x + middle) << std::endl;
+
+ middle = descDir;
+ middle *= m;
+
+ double derivative = directionalDerivative(x + middle, descDir);
+
+ if (derivative < 0) {
+ l = m;
+ m = (m + r) / 2;
+ } else if (derivative > 0) {
+ r = m;
+ m = (l + m) / 2;
+ } else
+ break;
+ }
+ return middle;
+ }
+
+private:
+ SmallMatrix A_;
+ SmallVector b_;
+};
+
+int main() {
+ int const dim = 2;
+ typedef SampleFunctional<dim> SampleFunctional;
+
+ SampleFunctional::SmallMatrix A;
+ A[0][0] = 3;
+ A[0][1] = 0;
+ A[1][0] = 0;
+ A[1][1] = 3;
+ SampleFunctional::SmallVector b;
+ b[0] = 1;
+ b[1] = 2;
+
+ SampleFunctional J(A, b);
+
+ SampleFunctional::SmallVector start = b;
+ start *= 17;
+ SampleFunctional::SmallVector correction = J.minimise(start, 20);
+ assert(J(start + correction) <= J(start));
+ std::cout << J(start + correction) << std::endl;
+
+ return 0;
+}
diff --git a/src/bisection-simpler-example2.cc b/src/bisection-simpler-example2.cc
new file mode 100644
index 00000000..dde6268c
--- /dev/null
+++ b/src/bisection-simpler-example2.cc
@@ -0,0 +1,58 @@
+#include <dune/fufem/interval.hh>
+
+#include <dune/tnnmg/nonlinearities/smallfunctional.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include <cassert>
+#include <limits>
+
+template <int dimension>
+class SampleFunctional : public SmallFunctional<dimension> {
+public:
+ typedef typename SmallFunctional<dimension>::SmallMatrix SmallMatrix;
+ typedef typename SmallFunctional<dimension>::SmallVector SmallVector;
+
+ SampleFunctional(SmallMatrix A, SmallVector b) : A_(A), b_(b) {}
+
+ double operator()(const SmallVector v) const {
+ SmallVector y;
+ A_.mv(v, y); // y = Av
+ y /= 2; // y = 1/2 Av
+ y -= b_; // y = 1/2 Av - b
+ return y * v; // <1/2 Av - b,v>
+ }
+
+ SmallVector d(const SmallVector v) const {
+ SmallVector y;
+ A_.mv(v, y); // y = Av
+ y -= b_; // y = Av - b
+ return y;
+ }
+
+ SmallMatrix d2(const SmallVector v) const { return A_; }
+
+private:
+ SmallMatrix A_;
+ SmallVector b_;
+};
+
+int main() {
+ int const dim = 2;
+
+ SampleFunctional<dim>::SmallMatrix A;
+ A[0][0] = 3;
+ A[0][1] = 0;
+ A[1][0] = 0;
+ A[1][1] = 3;
+ SampleFunctional<dim>::SmallVector b;
+ b[0] = 1;
+ b[1] = 2;
+
+ SampleFunctional<dim> J(A, b);
+
+ assert(J(b) == 2.5);
+
+ assert(b.two_norm() == sqrt(5));
+ assert(b.two_norm2() == 5);
+ return 0;
+}
--
GitLab