diff --git a/.gitignore b/.gitignore
index d312bf402f10cd055fbe6cbd1b52dc0ef359e545..e14662c1ab4521f3e1e41ef3c10078a7dc3399b7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,6 @@
 # Default cmake build directory
 /build-cmake
 
+# ignore Python files
+*.pyc
+
diff --git a/problems/film-on-substrate-deformation.parset b/problems/film-on-substrate-deformation.parset
new file mode 100644
index 0000000000000000000000000000000000000000..0a1c371b1c6de45d6f86369a35ff43a68ef8dfeb
--- /dev/null
+++ b/problems/film-on-substrate-deformation.parset
@@ -0,0 +1,133 @@
+#############################################
+#  Grid parameters
+#############################################
+
+structuredGrid = true
+lower = 0 0 0
+
+# whole experiment: 45 mm x 10 mm x 2 mm, scaling with 10^7 such that the thickness, which is around 100 nm, so 100x10^-9 = 10^-7 is equal to 1.
+# upper = 45e4 10e4 2e4
+# using only a section of the whole experiment as deformed grid to start with for dune-gfe: use much smaller dimensions!
+upper = 600 200 200
+
+elements = 15 5 5
+
+# Number of grid levels
+numLevels = 3
+
+adaptiveRefinement = true
+
+#############################################
+#  Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 500
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.1
+
+# Number of multigrid iterations per trust-region step
+numIt = 1000
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Damping for the smoothers of the multigrid solver
+damping = 1.0
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+#   Material parameters
+############################
+
+energy = mooneyrivlin # stvenantkirchhoff, neohooke, hencky, exphencky or mooneyrivlin
+
+[materialParameters]
+
+## Lame parameters for stvenantkirchhoff, E = mu(3*lambda + 2*mu)/(lambda + mu)
+#mu = 2.7191e+4
+#lambda = 4.4364e+4
+
+#mooneyrivlin_a = 2.7191e+6
+#mooneyrivlin_b = 2.7191e+6
+#mooneyrivlin_c = 2.7191e+6
+
+#mooneyrivlin_10 = -7.28e+5 #182 20:1
+#mooneyrivlin_01 = 9.17e+5
+#mooneyrivlin_20 = 1.23e+5
+#mooneyrivlin_02 = 8.15e+5
+#mooneyrivlin_11 = -5.14e+5
+
+#mooneyrivlin_10 = -3.01e+6 #182 2:1
+#mooneyrivlin_01 = 3.36e+6
+#mooneyrivlin_20 = 5.03e+6
+#mooneyrivlin_02 = 13.1e+6
+#mooneyrivlin_11 = -15.2e+6
+
+mooneyrivlin_10 = -1.67e+6 #184 2:1
+mooneyrivlin_01 = 1.94e+6
+mooneyrivlin_20 = 2.42e+6
+mooneyrivlin_02 = 6.52e+6
+mooneyrivlin_11 = -7.34e+6
+
+mooneyrivlin_30 = 0
+mooneyrivlin_21 = 0
+mooneyrivlin_12 = 0
+mooneyrivlin_03 = 0
+
+# volume-preserving parameter 
+mooneyrivlin_k = 75e+6
+
+# How to choose the volume-preserving parameter?
+# We need a stretch of 30%  (45e4 10e4 2e4 in x-direction, so a stretch of 45e4*0.3 = 13.5e4)
+# 184 2:1, mooneyrivlin_k = .. and mooneyrivlin_energy = square, neumannValues = 27e4 0 0
+# 184 2:1, mooneyrivlin_k = 57e+6 and mooneyrivlin_energy = log, neumannValues = 27e4 0 0 
+
+mooneyrivlin_energy = square # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
+
+# ciarlet: Fomula from "Ciarlet: Three-Dimensional Elasticity", here no penalty term is 
+# log: Generalized Rivlin model or polynomial hyperelastic model, using  0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
+# square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
+
+[]
+
+#############################################
+#  Boundary values
+#############################################
+
+problem = wriggers-l-shape
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+###  Python predicate specifying all neumannValues grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 599.99"
+
+###  Neumann values
+neumannValues = 27e4 0 0
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], x[2]]"
diff --git a/problems/wriggers-l-shape-dirichlet-values.py b/problems/wriggers-l-shape-dirichlet-values.py
new file mode 100644
index 0000000000000000000000000000000000000000..eac737a0f584c379922fdf950115dd5654daac28
--- /dev/null
+++ b/problems/wriggers-l-shape-dirichlet-values.py
@@ -0,0 +1,18 @@
+class DirichletValues:
+    def __init__(self, homotopyParameter):
+        self.homotopyParameter = homotopyParameter
+
+    # Deformation of 3d classical materials
+    def dirichletValues(self, x):
+        # Clamp the L-shape in its reference configuration
+        return [x[0], x[1], x[2]]
+
+    # Deformation of Cosserat shells
+    def deformation(self, x):
+        # Clamp the L-shape in its reference configuration
+        return [x[0], x[1], 0]
+
+    # Orientation of Cosserat materials
+    def orientation(self, x):
+        rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]
+        return rotation
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index cbc5a51d386f5f871020937dfaf1e289dfd08c4f..3138ee9decb5236604d851e2dafe65c77a9db36f 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -46,7 +46,7 @@
 // grid dimension
 const int dim = 3;
 
-const int order = 1;
+const int order = 2;
 
 using namespace Dune;
 
@@ -55,6 +55,8 @@ int main (int argc, char *argv[]) try
   // initialize MPI, finalize is done automatically on exit
   Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
 
+  if (mpiHelper.rank()==0)
+    std::cout << "ORDER = " << order << std::endl;
   // Start Python interpreter
   Python::start();
   Python::Reference main = Python::import("__main__");