diff --git a/cosserat-continuum-twisted-strip.parset b/cosserat-continuum-twisted-strip.parset
index 2aa1856601a018e80ad7c6e9dc308da442ae3a61..d2cdfd399e11c2e54fc2b7e9e10432035ce508c3 100644
--- a/cosserat-continuum-twisted-strip.parset
+++ b/cosserat-continuum-twisted-strip.parset
@@ -88,4 +88,7 @@ kappa = 1
 # x is the vertex coordinate
 dirichletVerticesPredicate = "x[0] < 0.001 or x[0] > 0.0999"
 
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0]"
+
 
diff --git a/cosserat-continuum-wong-pellegrino-wrinkling.parset b/cosserat-continuum-wong-pellegrino-wrinkling.parset
index 3d748989532e8cca885c7ff941dd0c6a7940debf..fba98912ec211bba521e40234d81850aca1127e2 100644
--- a/cosserat-continuum-wong-pellegrino-wrinkling.parset
+++ b/cosserat-continuum-wong-pellegrino-wrinkling.parset
@@ -137,3 +137,6 @@ neumannValues =  -1e3 0 0
 structuredGrid = false
 path = /home/sander/data/shells/wriggers_L_shape/
 gridFile = wriggers-L-shape_99.msh
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0.002*cos(1e4*x[0])]"
diff --git a/cosserat-continuum-wriggers-l-shape.parset b/cosserat-continuum-wriggers-l-shape.parset
index 89066ff02c225a43862835086714b3382d4e25d7..8e809176e8c4275c50c544dee5a7fb8ed8552c97 100644
--- a/cosserat-continuum-wriggers-l-shape.parset
+++ b/cosserat-continuum-wriggers-l-shape.parset
@@ -95,3 +95,5 @@ neumannVerticesPredicate = "x[1] < -0.239"
 ###  Neumann values, if needed
 neumannValues =  -1e3 0 0
 
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0]"
diff --git a/cosserat-continuum.cc b/cosserat-continuum.cc
index d507b7c449e7e5d45b932b500136735e588dab1c..c3ada0cb3e9f7b27773e1c0975a700b0355b5e0b 100644
--- a/cosserat-continuum.cc
+++ b/cosserat-continuum.cc
@@ -51,21 +51,6 @@ const int blocksize = TargetSpace::TangentVector::dimension;
 
 using namespace Dune;
 
-class Identity
-: public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3>>
-{
-public:
-  void evaluate(const FieldVector<double,dim>& x, FieldVector<double,3>& y) const
-  {
-    y = 0;
-    for (int i=0; i<dim; i++)
-      y[i] = x[i];
-
-    //y[2] = 0.002*std::cos(1e4*x[0]);
-  }
-};
-
-
 // Dirichlet boundary data for the shear/wrinkling example
 class WrinklingDirichletValues
 : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
@@ -312,9 +297,11 @@ int main (int argc, char *argv[]) try
 
     SolutionType x(feBasis.size());
 
-    Identity identity;
+    lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+    PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+
     std::vector<FieldVector<double,3> > v;
-    Functions::interpolate(feBasis, v, identity);
+    Functions::interpolate(feBasis, v, pythonInitialDeformation);
 
     for (size_t i=0; i<x.size(); i++)
       x[i].r = v[i];