Populate with sample code

src/Makefile.am

 
 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)
 

src/bisection-simpler-example.cc

+#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;
+}

src/bisection-simpler-example2-gradient.cc

+// 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;
+}

src/bisection-simpler-example2.cc

+#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;
+}