Control Neumann boundary conditions through python

src/Makefile.am

@@ -51,6 +51,7 @@ AM_CXXFLAGS = \
 
 AM_CPPFLAGS = \
 	$(DUNE_CPPFLAGS) \
+	$(PYTHON_CPPFLAGS) \
 	-I$(top_srcdir)
 
 # The libraries have to be given in reverse order (most basic libraries
@@ -58,9 +59,11 @@ AM_CPPFLAGS = \
 # -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
 # here as well.
 LDADD = \
-	$(DUNE_LDFLAGS) $(DUNE_LIBS)
+	$(DUNE_LDFLAGS) $(DUNE_LIBS) \
+	$(PYTHON_LIBS)
 AM_LDFLAGS = \
-	$(DUNE_LDFLAGS)
+	$(DUNE_LDFLAGS) \
+	$(PYTHON_LDFLAGS)
 
 # don't follow the full GNU-standard
 # we need automake 1.5

src/one-body-sample.cc

@@ -10,11 +10,17 @@
 #undef HAVE_IPOPT
 #endif
 
+#ifndef HAVE_PYTHON
+#error Python is required
+#endif
+
 #include <exception>
 #include <iostream>
 
 #include <boost/format.hpp>
 
+#include <Python.h>
+
 #include <cmath>
 
 #include "LambertW.h"
@@ -42,11 +48,13 @@
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/transferoperatorassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/functions/vtkbasisgridfunction.hh>
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/sharedpointermap.hh>
 #include <dune/solvers/common/numproc.hh> // Solver::FULL
 #include <dune/solvers/iterationsteps/blockgsstep.hh>
 #include <dune/solvers/iterationsteps/multigridstep.hh>
@@ -108,12 +116,12 @@ template <class GridType, class GridView, class LocalVectorType, class FEBasis>
 void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
                       Dune::BitSetVector<1> const &neumannNodes,
                       Dune::BlockVector<LocalVectorType> &f,
+                      Dune::VirtualFunction<double, double> const &neumann,
                       double time) { // constant sample function on neumann
                                      // boundary
   BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
   LocalVectorType SampleVector(0);
-  // FIXME: random values (time-dependent)
-  SampleVector[0] = (time <= 0.5) ? sin(time * 2 * M_PI) : (time - 0.5) * 2;
+  neumann.evaluate(time, SampleVector[0]);
   SampleVector[1] = 0;
   ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
       SampleVector);
@@ -225,6 +233,20 @@ double compute_state_update_lambert(double h, double unorm, double L,
 
 int main(int argc, char *argv[]) {
   try {
+    typedef SharedPointerMap<std::string, Dune::VirtualFunction<double, double>>
+    FunctionMap;
+    FunctionMap functions;
+    {
+      Python::start();
+
+      Python::run("import sys");
+      Python::run("sys.path.append('" srcdir "')");
+
+      Python::import("one-body-sample_neumann")
+          .get("Functions")
+          .toC<FunctionMap::Base>(functions);
+    }
+
     Dune::ParameterTree parset;
     Dune::ParameterTreeParser::readINITree(srcdir "/one-body-sample.parset",
                                            parset); // FIXME
@@ -380,7 +402,8 @@ int main(int argc, char *argv[]) {
 
       // b = neumann
       assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
-          leafView, p1Basis, neumannNodes, b1, h * run);
+          leafView, p1Basis, neumannNodes, b1, functions.get("sampleFunction"),
+          h * run);
       b2 = b1;
       b3 = b1;
       b4 = b1;
@@ -472,10 +495,10 @@ int main(int argc, char *argv[]) {
                                       "%|40t|u[%03d] = %+3e");
         std::cout << boost::format(formatter) % run % (*s4_new)[i] % run % u4[i]
                   << std::endl;
-        // FIXME: (Duplicate!) Magic values
-        octave_writer << (*s4_new)[i] << " " << u4[i][0] * 1e6 << " "
-                      << ((h * run <= 0.5) ? sin(h * run * 2 * M_PI)
-                                           : (h * run - 0.5) * 2) << std::endl;
+        double out;
+        functions.get("sampleFunction").evaluate(h * run, out);
+        octave_writer << (*s4_new)[i] << " " << u4[i][0] * 1e6 << " " << out
+                      << std::endl;
         break;
       }
 
@@ -573,6 +596,7 @@ int main(int argc, char *argv[]) {
           std::cout << boost::format(formatter) % i % u1[i] % i % u2[i] % i %
                            u3[i] % i % u4[i] << std::endl;
     }
+    Python::stop();
   }
   catch (Dune::Exception &e) {
     Dune::derr << "Dune reported error: " << e << std::endl;

src/one-body-sample_neumann.py

+class sampleFunction:
+    def __call__(self, x):
+        fst = 0.3
+        snd = 0.1
+        trd = 0.3
+        if x < 1.0/3:
+            return fst * x
+        elif x < 2.0/3:
+            return snd * (x - 1.0/3) + fst * 1.0/3
+        else:
+            return trd * (x - 2.0/3) + (fst + snd) * 1.0/3
+
+Functions = {'sampleFunction' : sampleFunction()}