cosserat-continuum.cc: Read input parameters as Python file

Rather than as our homegrown .parset.  That way, after more changes
this will allow to have the entire problem description in a single
file, including the boundary data.

problems/cosserat-continuum-cantilever.parset → problems/cosserat-continuum-cantilever.py

+import math
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+
 #############################################
 #  Grid parameters
 #############################################
 
-structuredGrid = cube
+parameterSet.structuredGrid = "cube"
 
 # bounding box
-lower = 0 0
-upper = 100 10
+parameterSet.lower = "0 0"
+parameterSet.upper = "100 10"
 
-elements = 10 1
+parameterSet.elements = "10 1"
 
 # Number of grid levels
-numLevels = 1
+parameterSet.numLevels = 1
 
 #############################################
 #  Solver parameters
 #############################################
 
 # Number of homotopy steps for the Dirichlet boundary conditions
-numHomotopySteps = 1
+parameterSet.numHomotopySteps = 1
 
-# Solver type: trustRegion or proximalNewton
-solvertype = trustRegion
+# Solver type: "trustRegion" or "proximalNewton"
+parameterSet.solvertype = "trustRegion"
 
 # Tolerance of the trust region solver
-tolerance = 1e-3
+parameterSet.tolerance = 1e-3
 
 # Max number of steps of the trust region solver
-maxSolverSteps = 200
+parameterSet.maxSolverSteps = 200
 
-solverScaling = 1 1 1 0.01 0.01 0.01
+parameterSet.solverScaling = "1 1 1 0.01 0.01 0.01"
 
 # Initial trust-region radius
-initialTrustRegionRadius = 3.125
+parameterSet.initialTrustRegionRadius = 3.125
 
 # Number of multigrid iterations per trust-region step
-numIt = 400
+parameterSet.numIt = 400
 
 # Number of presmoothing steps
-nu1 = 3
+parameterSet.nu1 = 3
 
 # Number of postsmoothing steps
-nu2 = 3
+parameterSet.nu2 = 3
 
 # Number of coarse grid corrections
-mu = 1
+parameterSet.mu = 1
 
 # Number of base solver iterations
-baseIt = 1
+parameterSet.baseIt = 1
 
 # Tolerance of the multigrid solver
-mgTolerance = 1e-5
+parameterSet.mgTolerance = 1e-5
 
 # Tolerance of the base grid solver
-baseTolerance = 1e-8
+parameterSet.baseTolerance = 1e-8
 
 # Measure convergence
-instrumented = 0
+parameterSet.instrumented = 0
 
 ############################
 #   Material parameters
 ############################
 
-
-## For the Wriggers L-shape example
-[materialParameters]
+parameterSet.materialParameters = ParameterSet()
 
 # shell thickness
-thickness = 0.6
+parameterSet.materialParameters.thickness = 0.6
 
 # Lame parameters
 # corresponds to E = 71240 N/mm^2, nu=0.31
 # However, we use units N/m^2
-mu = 2.7191e+4
-lambda = 4.4364e+4
+parameterSet.materialParameters.mu = 2.7191e+4
+parameterSet.materialParameters["lambda"] = 4.4364e+4
 
 # Cosserat couple modulus
-mu_c = 0
+parameterSet.materialParameters.mu_c = 0
 
 # Length scale parameter
-L_c = 0.6
+parameterSet.materialParameters.L_c = 0.6
 
 # Curvature exponent
-q = 2
+parameterSet.materialParameters.q = 2
 
 # Shear correction factor
-kappa = 1
+parameterSet.materialParameters.kappa = 1
+
+# TODO: These three parameters are not actually used,
+# but the current implementation requires them nevertheless.
+parameterSet.materialParameters.b1 = 1
+parameterSet.materialParameters.b2 = 1
+parameterSet.materialParameters.b3 = 1
 
-[]
 
 #############################################
 #  Boundary values
 #############################################
 
-problem = cantilever
+parameterSet.problem = "cantilever"
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
-dirichletVerticesPredicate = "x[0] < 0.01"
+parameterSet.dirichletVerticesPredicate = "[x[0] < 0.01, x[0] < 0.01, x[0] < 0.01]"
+parameterSet.dirichletRotationVerticesPredicate = "x[0] < 0.01"
 
 ###  Python predicate specifying all Neumann grid vertices
 # x is the vertex coordinate
-neumannVerticesPredicate = "x[0] > 99.99"
+parameterSet.neumannVerticesPredicate = "x[0] > 99.99"
 
 ###  Neumann values
-neumannValues = 0 0 3
+parameterSet.neumannValues = "0 0 3"
 
 # Initial deformation
-initialDeformation = "[x[0], x[1], 0]"
+parameterSet.initialDeformation = "[x[0], x[1], 0]"
 
-#startFromFile = yes
-#initialIterateFilename = initial_iterate.vtu
+#parameterSet.startFromFile = yes
+#parameterSet.initialIterateFilename = initial_iterate.vtu

problems/cosserat-continuum-twisted-strip.parset → problems/cosserat-continuum-twisted-strip.py

+import math
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+
 #############################################
 #  Grid parameters
 #############################################
 
-structuredGrid = cube
-lower = 0  -0.005
-upper = 0.1 0.005
-elements = 10 1
+parameterSet.structuredGrid = "cube"
+parameterSet.lower = "0  -0.005"
+parameterSet.upper = "0.1 0.005"
+parameterSet.elements = "10 1"
 
 # Number of grid levels
-numLevels = 1
+parameterSet.numLevels = 1
 
 #############################################
 #  Solver parameters
 #############################################
 
 # Number of homotopy steps for the Dirichlet boundary conditions
-numHomotopySteps = 24
+parameterSet.numHomotopySteps = 24
 
 # Tolerance of the trust region solver
-tolerance = 1e-8
+parameterSet.tolerance = 1e-8
 
 # Max number of steps of the trust region solver
-maxSolverSteps = 200
+parameterSet.maxSolverSteps = 200
 
-solverScaling = 1 1 1 1 1 1
+parameterSet.solverScaling = "1 1 1 1 1 1"
 
 # Initial trust-region radius
-initialTrustRegionRadius = 0.1
+parameterSet.initialTrustRegionRadius = 0.1
 
 # Number of multigrid iterations per trust-region step
-numIt = 200
+parameterSet.numIt = 200
 
 # Number of presmoothing steps
-nu1 = 3
+parameterSet.nu1 = 3
 
 # Number of postsmoothing steps
-nu2 = 3
+parameterSet.nu2 = 3
 
 # Number of coarse grid corrections
-mu = 1
+parameterSet.mu = 1
 
 # Number of base solver iterations
-baseIt = 100
+parameterSet.baseIt = 100
 
 # Tolerance of the multigrid solver
-mgTolerance = 1e-7
+parameterSet.mgTolerance = 1e-7
 
 # Tolerance of the base grid solver
-baseTolerance = 1e-8
+parameterSet.baseTolerance = 1e-8
 
 # Measure convergence
-instrumented = 0
+parameterSet.instrumented = 0
 
 ############################
 #   Material parameters
@@ -58,48 +67,50 @@ instrumented = 0
 
 
 # Parameters for the twisted 0.1x0.01 strip
-[materialParameters]
+parameterSet.materialParameters = ParameterSet()
 
 # shell thickness
-thickness = 0.002
+parameterSet.materialParameters.thickness = 0.002
 
 # Lame parameters
 # corresponds to E = 3.5GPa, nu=0.31
-mu = 5.6452e+09
-lambda = 2.1796e+09
+# (Notice that lambda is a key word in python such that
+# we cannot use attribute syntax to set it but have to
+# resort to dictionary syntax.)
+parameterSet.materialParameters.mu = 5.6452e+09
+parameterSet.materialParameters['lambda'] = 2.1796e+09
 
 # Cosserat couple modulus
-mu_c = 0
+parameterSet.materialParameters.mu_c = 0
 
 # Length scale parameter
-L_c = 2e-6
+parameterSet.materialParameters.L_c = 2e-6
 
 # Curvature exponent
-q = 2
+parameterSet.materialParameters.q = 2
 
 # Shear correction factor
-kappa = 1
+parameterSet.materialParameters.kappa = 1
 
 # TODO: These three parameters are not actually used,
 # but the current implementation requires them nevertheless.
-b1 = 1
-b2 = 1
-b3 = 1
+parameterSet.materialParameters.b1 = 1
+parameterSet.materialParameters.b2 = 1
+parameterSet.materialParameters.b3 = 1
 
-[]
 
 #############################################
 #  Boundary values
 #############################################
 
-problem = twisted-strip
+parameterSet.problem = "twisted-strip"
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
-dirichletVerticesPredicate = "[x[0] < 0.001 or x[0] > 0.0999, x[0] < 0.001 or x[0] > 0.0999, x[0] < 0.001 or x[0] > 0.0999]"
-dirichletRotationVerticesPredicate = "x[0] < 0.001 or x[0] > 0.0999"
+parameterSet.dirichletVerticesPredicate = "[x[0] < 0.001 or x[0] > 0.0999, x[0] < 0.001 or x[0] > 0.0999, x[0] < 0.001 or x[0] > 0.0999]"
+parameterSet.dirichletRotationVerticesPredicate = "x[0] < 0.001 or x[0] > 0.0999"
 
 # Initial deformation
-initialDeformation = "[x[0], x[1], 0]"
+parameterSet.initialDeformation = "[x[0], x[1], 0]"
 
 

problems/cosserat-continuum-wong-pellegrino-wrinkling.parset → problems/cosserat-continuum-wong-pellegrino-wrinkling.py

+import math
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+
 #############################################
 #  Grid parameters
 #############################################
 
-structuredGrid = true
-lower = 0 0
-upper = 0.38 0.128
-elements = 15 5
+parameterSet.structuredGrid = "cube"
+parameterSet.lower = "0 0"
+parameterSet.upper = "0.38 0.128"
+parameterSet.elements = "15 5"
 
 # Number of grid levels
-numLevels = 1
+parameterSet.numLevels = 1
 
 #############################################
 #  Solver parameters
 #############################################
 
 # Number of homotopy steps for the Dirichlet boundary conditions
-numHomotopySteps = 1
+parameterSet.numHomotopySteps = 1
 
 # Tolerance of the trust region solver
-tolerance = 1e-8
+parameterSet.tolerance = 1e-8
 
 # Max number of steps of the trust region solver
-maxTrustRegionSteps = 500
+parameterSet.maxSolverSteps = 500
 
-solverScaling = 1 1 1 0.01 0.01 0.01
+parameterSet.solverScaling = "1 1 1 0.01 0.01 0.01"
 
 # Initial trust-region radius
-initialTrustRegionRadius = 0.001
+parameterSet.initialTrustRegionRadius = 0.001
 
 # Number of multigrid iterations per trust-region step
-numIt = 200
+parameterSet.numIt = 200
 
 # Number of presmoothing steps
-nu1 = 3
+parameterSet.nu1 = 3
 
 # Number of postsmoothing steps
-nu2 = 3
+parameterSet.nu2 = 3
 
 # Number of coarse grid corrections
-mu = 1
+parameterSet.mu = 1
 
 # Number of base solver iterations
-baseIt = 100
+parameterSet.baseIt = 100
 
 # Tolerance of the multigrid solver
-mgTolerance = 1e-7
+parameterSet.mgTolerance = 1e-7
 
 # Tolerance of the base grid solver
-baseTolerance = 1e-8
+parameterSet.baseTolerance = 1e-8
 
 # Measure convergence
-instrumented = 0
+parameterSet.instrumented = 0
+
 
 ############################
 #   Material parameters
@@ -59,41 +69,47 @@ instrumented = 0
 
 # Parameters for the shearing/wrinkling example from Wong/Pellegrino 2006
 # We use 'meters' as the length unit
-[materialParameters]
+parameterSet.materialParameters = ParameterSet()
 
 # shell thickness
-thickness = 2.5e-5
+parameterSet.materialParameters.thickness = 2.5e-5
 
 # Lame parameters
 # corresponds to E = 3.5GPa, nu=0.31
-mu = 5.6452e+09
-lambda = 2.1796e+09
+parameterSet.materialParameters.mu = 5.6452e+09
+parameterSet.materialParameters["lambda"] = 2.1796e+09
 
 # Cosserat couple modulus
-mu_c = 0
+parameterSet.materialParameters.mu_c = 0
 
 # Length scale parameter
-L_c = 2.5e-8
+parameterSet.materialParameters.L_c = 2.5e-8
 
 # Curvature exponent
-q = 2
+parameterSet.materialParameters.q = 2
 
 # Shear correction factor
-kappa = 1
+parameterSet.materialParameters.kappa = 1
+
+# TODO: These three parameters are not actually used,
+# but the current implementation requires them nevertheless.
+parameterSet.materialParameters.b1 = 1
+parameterSet.materialParameters.b2 = 1
+parameterSet.materialParameters.b3 = 1
 
-[]
 
 #############################################
 #  Boundary values
 #############################################
 
-problem = wong-pellegrino
+parameterSet.problem = "wong-pellegrino"
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
-dirichletVerticesPredicate = "x[1] < 0.0001 or x[1] > 0.128 - 0.0001"
+parameterSet.dirichletVerticesPredicate = "[x[1] < 0.0001 or x[1] > 0.128 - 0.0001, x[1] < 0.0001 or x[1] > 0.128 - 0.0001, x[1] < 0.0001 or x[1] > 0.128 - 0.0001]"
+parameterSet.dirichletRotationVerticesPredicate = "x[1] < 0.0001 or x[1] > 0.128 - 0.0001"
 
 # Initial deformation
-#startFromFile = true
-initialIterateFilename = cosserat_iterate_2.vtu
-initialDeformation = "[x[0] + 0.003*x[1] / 0.128, x[1], 0.002*math.cos(1e4*x[0])]"
+#parameterSet.startFromFile = True
+parameterSet.initialIterateFilename = "cosserat_iterate_2.vtu"
+parameterSet.initialDeformation = "[x[0] + 0.003*x[1] / 0.128, x[1], 0.002*math.cos(1e4*x[0])]"

problems/cosserat-continuum-wriggers-l-shape.parset → problems/cosserat-continuum-wriggers-l-shape.py

+import math
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+
 #############################################
 #  Grid parameters
 #############################################
 
-structuredGrid = false
-path = /home/sander/data/shells/wriggers_L_shape/
-gridFile = wriggers-L-shape_99_mm.msh
+parameterSet.structuredGrid = "false"
+parameterSet.path = "/home/sander/data/shells/wriggers_L_shape/"
+parameterSet.gridFile = "wriggers-L-shape_99_mm.msh"
 
 # Number of grid levels
-numLevels = 2
+parameterSet.numLevels = 2
 
 #############################################
 #  Solver parameters
 #############################################
 
 # Number of homotopy steps for the Dirichlet boundary conditions
-numHomotopySteps = 1
+parameterSet.numHomotopySteps = 1
 
 # Tolerance of the trust region solver
-tolerance = 1e-8
+parameterSet.tolerance = 1e-8
 
 # Max number of steps of the trust region solver
-maxSolverSteps = 1000
+parameterSet.maxSolverSteps = 1000
 
-solverScaling = 1 1 1 0.01 0.01 0.01
+parameterSet.solverScaling = "1 1 1 0.01 0.01 0.01"
 
 # Initial trust-region radius
-initialTrustRegionRadius = 1
+parameterSet.initialTrustRegionRadius = 1
 
 # Number of multigrid iterations per trust-region step
-numIt = 200
+parameterSet.numIt = 200
 
 # Number of presmoothing steps
-nu1 = 3
+parameterSet.nu1 = 3
 
 # Number of postsmoothing steps
-nu2 = 3
+parameterSet.nu2 = 3
 
 # Number of coarse grid corrections
-mu = 1
+parameterSet.mu = 1
 
 # Number of base solver iterations
-baseIt = 1
+parameterSet.baseIt = 1
 
 # Tolerance of the multigrid solver
-mgTolerance = 1e-7
+parameterSet.mgTolerance = 1e-7
 
 # Tolerance of the base grid solver
-baseTolerance = 1e-8
+parameterSet.baseTolerance = 1e-8
 
 # Measure convergence
-instrumented = 0
+parameterSet.instrumented = 0
 
 ############################
 #   Material parameters
 ############################
 
-
-## For the Wriggers L-shape example
-[materialParameters]
+parameterSet.materialParameters = ParameterSet()
 
 # shell thickness
-thickness = 0.6
+parameterSet.materialParameters.thickness = 0.6
 
 # Lame parameters
 # corresponds to E = 71240 N/mm^2, nu=0.31
 # However, we use units N/m^2
-mu = 2.7191e+4
-lambda = 4.4364e+4
+parameterSet.materialParameters.mu = 2.7191e+4
+parameterSet.materialParameters["lambda"] = 4.4364e+4
 
 # Cosserat couple modulus
-mu_c = 0
+parameterSet.materialParameters.mu_c = 0
 
 # Length scale parameter
-L_c = 0.6e-3
+parameterSet.materialParameters.L_c = 0.6e-3
 
 # Curvature exponent
-q = 2
+parameterSet.materialParameters.q = 2
 
 # Shear correction factor
-kappa = 1
+parameterSet.materialParameters.kappa = 1
+
+# TODO: These three parameters are not actually used,
+# but the current implementation requires them nevertheless.
+parameterSet.materialParameters.b1 = 1
+parameterSet.materialParameters.b2 = 1
+parameterSet.materialParameters.b3 = 1
 
-[]
 
 #############################################
 #  Boundary values
 #############################################
 
-problem = wriggers-l-shape
+parameterSet.problem = "wriggers-l-shape"
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
-dirichletVerticesPredicate = "x[0] < 1"
+parameterSet.dirichletVerticesPredicate = "[x[0] < 1, x[0] < 1, x[0] < 1]"
+parameterSet.dirichletRotationVerticesPredicate = "x[0] < 1"
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
-neumannVerticesPredicate = "x[1] < -239"
+parameterSet.neumannVerticesPredicate = "x[1] < -239"
 
 ###  Neumann values, if needed
-neumannValues =  0.09 0 0
+parameterSet.neumannValues =  "0.09 0 0"
 
 # Initial deformation
-#initialDeformation = "[x[0], x[1], 0 if (x[0] < 225 or x[1] < -15) else 0.001*(x[0]-225)*(x[1]+15)]"
-##initialDeformation = "[x[0], x[1], 0]"
+#parameterSet.initialDeformation = "[x[0], x[1], 0 if (x[0] < 225 or x[1] < -15) else 0.001*(x[0]-225)*(x[1]+15)]"
+##parameterSet.initialDeformation = "[x[0], x[1], 0]"
 
-startFromFile = yes
-initialIterateGridFilename = wriggers-L-shape_99_mm.msh
-initialIterateFilename = initial-wriggers-l-shape.vtu
+parameterSet.startFromFile = True
+parameterSet.initialIterateGridFilename = "wriggers-L-shape_99_mm.msh"
+parameterSet.initialIterateFilename = "initial-wriggers-l-shape.vtu"

src/cosserat-continuum.cc

@@ -127,7 +127,7 @@ int main (int argc, char *argv[]) try
 
   // Start Python interpreter
   Python::start();
-  Python::Reference main = Python::import("__main__");
+  auto pyMain = Python::main();
   Python::run("import math");
 
   //feenableexcept(FE_INVALID);
@@ -136,13 +136,12 @@ int main (int argc, char *argv[]) try
     << std::endl << "sys.path.append('" << argv[1] << "')"
     << std::endl;
 
-  using namespace TypeTree::Indices;
-  using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
-      std::vector<Rotation<double,3> > >;
-
   // parse data file
+  auto pyModule = pyMain.import(argv[2]);
+
+  // Get main parameter set
   ParameterTree parameterSet;
-  ParameterTreeParser::readINITree(argv[2], parameterSet);
+  pyModule.get("parameterSet").toC(parameterSet);
 
   ParameterTreeParser::readOptions(argc, argv, parameterSet);
 
@@ -164,6 +163,10 @@ int main (int argc, char *argv[]) try
   const bool adolcScalarMode            = parameterSet.get<bool>("adolcScalarMode", false);
   const std::string resultPath          = parameterSet.get("resultPath", "");
 
+  using namespace Dune::Indices;
+  using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
+      std::vector<Rotation<double,3> > >;
+
   // ///////////////////////////////////////
   //    Create the grid
   // ///////////////////////////////////////