The Dirichlet values classes are in separate Python files now

[[Imported from SVN: r9888]]

cosserat-continuum.cc

@@ -50,89 +50,6 @@ const int blocksize = TargetSpace::TangentVector::dimension;
 
 using namespace Dune;
 
-class WongPellegrinoDirichletValuesPythonWriter
-{
-  FieldVector<double,2> upper_;
-  double homotopy_;
-public:
-
-  WongPellegrinoDirichletValuesPythonWriter(FieldVector<double,2> upper, double homotopy)
-  : upper_(upper), homotopy_(homotopy)
-  {}
-
-  void writeDeformation()
-  {
-    Python::runStream()
-        << std::endl << "def deformationDirichletValues(x):"
-        << std::endl << "    out = [x[0], x[1], 0]"
-        << std::endl << "    if x[1] > " << upper_[1]-1e-4 << ":"
-        << std::endl << "        out[0] += 0.003 * " << homotopy_
-        << std::endl << "    return out";
-  }
-
-  void writeOrientation()
-  {
-    Python::runStream()
-        << std::endl << "def orientationDirichletValues(x):"
-        << std::endl << "    rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]"
-        << std::endl << "    return rotation";
-  }
-
-};
-
-
-class TwistedStripDirichletValuesPythonWriter
-{
-  FieldVector<double,2> upper_;
-  double homotopy_;
-  double totalAngle_;
-public:
-
-  TwistedStripDirichletValuesPythonWriter(FieldVector<double,2> upper, double homotopy)
-  : upper_(upper), homotopy_(homotopy)
-  {
-    totalAngle_ = 6*M_PI;
-  }
-
-  void writeDeformation()
-  {
-    Python::runStream()
-        << std::endl << "def deformationDirichletValues(x):"
-        << std::endl << "    upper = [" << upper_[0] << ", " << upper_[1] << "]"
-        << std::endl << "    homotopy = " << homotopy_
-        << std::endl << "    angle = " << totalAngle_ << " * x[0]/upper[0]"
-        << std::endl << "    angle *= homotopy"
-
-        // center of rotation
-        << std::endl << "    center = [0, 0, 0]"
-        << std::endl << "    center[1] = upper[1]/2.0"
-
-        << std::endl << "    rotation = [[1,0,0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]"
-
-        << std::endl << "    inEmbedded = [x[0]-center[0], x[1]-center[1], 0-center[2]]"
-
-        // Matrix-vector product
-        << std::endl << "    out = [rotation[0][0]*inEmbedded[0]+rotation[0][1]*inEmbedded[1], rotation[1][0]*inEmbedded[0]+rotation[1][1]*inEmbedded[1], rotation[2][0]*inEmbedded[0]+rotation[2][1]*inEmbedded[1]]"
-
-        << std::endl << "    out = [out[0]+center[0], out[1]+center[1], out[2]+center[2]]"
-        << std::endl << "    return out";
-  }
-
-  void writeOrientation()
-  {
-    Python::runStream()
-        << std::endl << "def orientationDirichletValues(x):"
-        << std::endl << "    upper = [" << upper_[0] << ", " << upper_[1] << "]"
-        << std::endl << "    angle = " << totalAngle_ << " * x[0]/upper[0];"
-        << std::endl << "    angle *= " << homotopy_
-
-        << std::endl << "    rotation = [[1,0,0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]"
-        << std::endl << "    return rotation";
-  }
-
-};
-
-
 /** \brief A constant vector-valued function, for simple Neumann boundary values */
 struct NeumannFunction
     : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
@@ -164,6 +81,10 @@ int main (int argc, char *argv[]) try
     Python::run("import math");
 
     //feenableexcept(FE_INVALID);
+    Python::runStream()
+        << std::endl << "import sys"
+        << std::endl << "sys.path.append('/home/sander/dune/dune-gfe/')"
+        << std::endl;
 
     typedef std::vector<TargetSpace> SolutionType;
 
@@ -361,26 +282,16 @@ int main (int argc, char *argv[]) try
         //   Set Dirichlet values
         ////////////////////////////////////////////////////////
 
-        if (parameterSet.get<std::string>("problem") == "twisted-strip")
-        {
-          TwistedStripDirichletValuesPythonWriter twistedStripDirichletValuesPythonWriter(upper, homotopyParameter);
-          twistedStripDirichletValuesPythonWriter.writeDeformation();
-          twistedStripDirichletValuesPythonWriter.writeOrientation();
-
-        } else if (parameterSet.get<std::string>("problem") == "wong-pellegrino")
-        {
-          WongPellegrinoDirichletValuesPythonWriter wongPellegrinoDirichletValuesPythonWriter(upper, homotopyParameter);
-          wongPellegrinoDirichletValuesPythonWriter.writeDeformation();
-          wongPellegrinoDirichletValuesPythonWriter.writeOrientation();
+        Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
 
-        } else if (parameterSet.get<std::string>("problem") == "wriggers-L-shape")
-        {
+        Python::Callable C = dirichletValuesClass.get("DirichletValues");
 
-        } else
-          DUNE_THROW(Exception, "Unknown problem type");
+        // Call a constructor.
+        Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
-        PythonFunction<FieldVector<double,dim>, FieldVector<double,3> >   deformationDirichletValues(main.get("deformationDirichletValues"));
-        PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(main.get("orientationDirichletValues"));
+        // Extract object member functions as Dune functions
+        PythonFunction<FieldVector<double,dim>, FieldVector<double,3> >   deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
+        PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
 
         std::vector<FieldVector<double,3> > ddV;
         Functions::interpolate(feBasis, ddV, deformationDirichletValues, dirichletDofs);

twisted-strip-dirichlet-values.py

+import math
+
+class DirichletValues:
+    def __init__(self, homotopyParameter):
+        self.homotopyParameter = homotopyParameter
+        self.upper = [0.1, 0.01]
+        self.totalAngle = 6*math.pi
+
+    def deformation(self, x):
+        angle = self.totalAngle * x[0]/self.upper[0]
+        angle *= self.homotopyParameter
+
+        # Center of rotation
+        center = [0, 0, 0]
+        center[1] = self.upper[1]/2.0
+
+        rotation = [[1,0,0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]
+
+        inEmbedded = [x[0]-center[0], x[1]-center[1], 0-center[2]]
+
+        # Matrix-vector product
+        out = [rotation[0][0]*inEmbedded[0]+rotation[0][1]*inEmbedded[1], rotation[1][0]*inEmbedded[0]+rotation[1][1]*inEmbedded[1], rotation[2][0]*inEmbedded[0]+rotation[2][1]*inEmbedded[1]]
+
+        out = [out[0]+center[0], out[1]+center[1], out[2]+center[2]]
+        return out
+
+
+    def orientation(self, x):
+        angle = self.totalAngle * x[0]/self.upper[0]
+        angle *= self.homotopyParameter
+
+        rotation = [[1,0,0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]
+        return rotation
+

wong-pellegrino-dirichlet-values.py

+class DirichletValues:
+    def __init__(self, homotopyParameter):
+        self.homotopyParameter = homotopyParameter
+
+    def deformation(self, x):
+        out = [x[0], x[1], 0]
+        if x[1] >  0.128-1e-4 :
+            out[0] += 0.003 * self.homotopyParameter
+        return out
+
+    def orientation(self, x):
+        rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]
+        return rotation