diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index ce516ace831adcd4cfe47103c98ff78e49a521e2..101c37b9986b73c96746603e2e8014d9bf130d7d 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -210,17 +210,18 @@ public:
   /** \brief Constructor with a set of material parameters
    * \param parameters The material parameters
    */
-  SurfaceCosseratEnergy(const Dune::ParameterTree& parameters, const std::vector<UnitVector<double,3> >& vertexNormals, const BoundaryPatch<GridView>* shellBoundary)
+  SurfaceCosseratEnergy(const Dune::ParameterTree& parameters,
+    const std::vector<UnitVector<double,3> >& vertexNormals,
+    const BoundaryPatch<GridView>* shellBoundary,
+    const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dim-1, dim>>& geometriesOnShellBoundary,
+    const std::function<double(Dune::FieldVector<double,dim>)> thicknessF,
+    const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF)
   : shellBoundary_(shellBoundary),
-    vertexNormals_(vertexNormals)
+    vertexNormals_(vertexNormals),
+    geometriesOnShellBoundary_(geometriesOnShellBoundary),
+    thicknessF_(thicknessF),
+    lameF_(lameF)
   {
-    // The shell thickness
-    thickness_ = parameters.template get<double>("thickness");
-
-    // Lame constants
-    mu_ = parameters.template get<double>("mu_cosserat");
-    lambda_ = parameters.template get<double>("lambda_cosserat");
-
     // Cosserat couple modulus
     mu_c_ = parameters.template get<double>("mu_c");
 
@@ -243,7 +244,6 @@ RT energy(const typename Basis::LocalView& localView,
 {
   // The element geometry
   auto element = localView.element();
-  auto geometry = element.geometry();
 
   // The set of shape functions on this element
   const auto& localFiniteElement = localView.tree().finiteElement();
@@ -273,10 +273,14 @@ RT energy(const typename Basis::LocalView& localView,
 
   RT energy = 0;
 
+  auto& idSet = gridView.grid().globalIdSet();
+
   for (auto&& it : intersections(shellBoundary_->gridView(), element)) {
     if (not shellBoundary_->contains(it))
       continue;
     
+    auto id = idSet.subId(it.inside(), it.indexInInside(), 1);
+    auto boundaryGeometry = geometriesOnShellBoundary_.at(id);
     auto quadOrder = (it.type().isSimplex()) ? localFiniteElement.localBasis().order()
                                                   : localFiniteElement.localBasis().order() * gridDim;
 
@@ -286,7 +290,14 @@ RT energy(const typename Basis::LocalView& localView,
       // Local position of the quadrature point
       const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());;
 
-      const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
+      // Global position of the quadrature point
+      auto quadPosGlobal = it.geometry().global(quad[pt].position());
+      double thickness = thicknessF_(quadPosGlobal);
+      auto lameConstants = lameF_(quadPosGlobal);
+      double mu = lameConstants[0];
+      double lambda = lameConstants[1];
+
+      const DT integrationElement = boundaryGeometry.integrationElement(quad[pt].position());
 
       // The value of the local function
       RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
@@ -321,7 +332,7 @@ RT energy(const typename Basis::LocalView& localView,
 
       // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
       // of the element.  If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
-      const auto jacobianTransposed = it.geometry().jacobianTransposed(quad[pt].position());
+      const auto jacobianTransposed = boundaryGeometry.jacobianTransposed(quad[pt].position());
       // auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
 
       for (int i=0; i<2; i++)
@@ -409,15 +420,15 @@ RT energy(const typename Basis::LocalView& localView,
       //////////////////////////////////////////////////////////
 
       // Add the membrane energy density
-      auto energyDensity = (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_m(Ee);
-      energyDensity += (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_m(Ee*b + c*Ke);
-      energyDensity += Dune::Power<3>::eval(thickness_) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke);
-      energyDensity += Dune::Power<5>::eval(thickness_) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
+      auto energyDensity = (thickness - K*Dune::Power<3>::eval(thickness) / 12.0) * W_m(Ee, mu, lambda);
+      energyDensity += (Dune::Power<3>::eval(thickness) / 12.0 - K * Dune::Power<5>::eval(thickness) / 80.0)*W_m(Ee*b + c*Ke, mu, lambda);
+      energyDensity += Dune::Power<3>::eval(thickness) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke, mu, lambda);
+      energyDensity += Dune::Power<5>::eval(thickness) / 80.0 * W_mp( (Ee*b + c*Ke)*b, mu, lambda);
 
       // Add the bending energy density
-      energyDensity += (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_curv(Ke)
-                     + (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_curv(Ke*b)
-                     + Dune::Power<5>::eval(thickness_) / 80.0 * W_curv(Ke*b*b);
+      energyDensity += (thickness - K*Dune::Power<3>::eval(thickness) / 12.0) * W_curv(Ke, mu)
+                     + (Dune::Power<3>::eval(thickness) / 12.0 - K * Dune::Power<5>::eval(thickness) / 80.0)*W_curv(Ke*b, mu)
+                     + Dune::Power<5>::eval(thickness) / 80.0 * W_curv(Ke*b*b, mu);
 
       // Add energy density
       energy += quad[pt].weight() * integrationElement * energyDensity;
@@ -427,37 +438,43 @@ RT energy(const typename Basis::LocalView& localView,
   return energy;
 }
 
-  RT W_m(const Dune::FieldMatrix<field_type,3,3>& S) const
+  RT W_m(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
   {
-    return W_mixt(S,S);
+    return W_mixt(S,S, mu, lambda);
   }
 
-  RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T) const
+  RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T, double mu, double lambda) const
   {
-    return mu_ * frobeniusProduct(sym(S), sym(T))
+    return mu * frobeniusProduct(sym(S), sym(T))
          + mu_c_ * frobeniusProduct(skew(S), skew(T))
-         + lambda_ * mu_ / (lambda_ + 2*mu_) * trace(S) * trace(T);
+         + lambda * mu / (lambda + 2*mu) * trace(S) * trace(T);
   }
 
-  RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S) const
+  RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
   {
-    return mu_ * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda_ * 0.5 * traceSquared(S);
+    return mu * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda * 0.5 * traceSquared(S);
   }
 
-  RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S) const
+  RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S, double mu) const
   {
-    return mu_ * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
+    return mu * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
   }
 
 private:
-  /** \brief The Neumann boundary */
+  /** \brief The shell boundary */
   const BoundaryPatch<GridView>* shellBoundary_;
 
-  /** \brief The normal vectors at the grid vertices.  This are used to compute the reference surface curvature. */
+  /** \brief Stress-free geometries of the shell elements*/
+  const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dim-1, dim>> geometriesOnShellBoundary_;
+
+  /** \brief The normal vectors at the grid vertices. They are used to compute the reference surface curvature. */
   std::vector<UnitVector<double,3> > vertexNormals_;
 
-  /** \brief The shell thickness */
-  double thickness_;
+  /** \brief The shell thickness as a function*/
+  std::function<double(Dune::FieldVector<double,dim>)> thicknessF_;
+
+  /** \brief The Lamé-parameters as a function*/
+  std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF_;
 
   /** \brief Lame constants */
   double mu_, lambda_;
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index c88e018ced274672058f297fe8efff1e5d2232d6..471ea478564f912ead8f141f3acc0c03caf597d9 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -68,6 +68,9 @@
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
+#include <iostream>
+#include <fstream>
+
 // grid dimension
 #ifndef WORLD_DIM
 #  define WORLD_DIM 3
@@ -149,7 +152,7 @@ int main (int argc, char *argv[]) try
   ParameterTreeParser::readOptions(argc, argv, parameterSet);
 
   // read solver settings
-  int numLevels                   = parameterSet.get<int>("numLevels");
+  int numLevels                         = parameterSet.get<int>("numLevels");
   int numHomotopySteps                  = parameterSet.get<int>("numHomotopySteps");
   const double tolerance                = parameterSet.get<double>("tolerance");
   const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
@@ -162,7 +165,8 @@ int main (int argc, char *argv[]) try
   const double mgTolerance              = parameterSet.get<double>("mgTolerance");
   const double baseTolerance            = parameterSet.get<double>("baseTolerance");
   const bool instrumented               = parameterSet.get<bool>("instrumented");
-  std::string resultPath                = parameterSet.get("resultPath", "");
+  const bool startFromFile              = parameterSet.get<bool>("startFromFile");
+  const std::string resultPath          = parameterSet.get("resultPath", "");
 
   // ///////////////////////////////////////
   //    Create the grid
@@ -198,7 +202,7 @@ int main (int argc, char *argv[]) try
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
   PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
 
-  // Same for the boundary that will get wrinkled
+  // Same for the Surface Shell Boundary
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("surfaceShellVerticesPredicate", "0") + std::string(")");
   PythonFunction<FieldVector<double,dim>, bool> pythonSurfaceShellVertices(Python::evaluate(lambda));
 
@@ -242,7 +246,6 @@ int main (int argc, char *argv[]) try
 
   const GridView::IndexSet& indexSet = gridView.indexSet();
 
-
   for (auto&& v : vertices(gridView))
   {
     bool isDirichlet = pythonDirichletVertices(v.geometry().corner(0));
@@ -306,14 +309,14 @@ int main (int argc, char *argv[]) try
   
   SolutionTypeCosserat x(feBasis.size());
 
-  //Solution in 3D, without the Cosserat directors on the boundary
   BlockVector<FieldVector<double,3> > v;
   
+  //Initial deformation of the underlying substrate
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
   PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
   ::Functions::interpolate(fufemFEBasis, v, pythonInitialDeformation);
 
-  //Copy over the parts of the boundary
+  //Copy over the initial deformation
   for (size_t i=0; i<x.size(); i++) {
       x[i].r = v[i];
   }
@@ -349,9 +352,95 @@ int main (int argc, char *argv[]) try
   std::cout << "Neumann values: " << neumannValues << std::endl;
 
   //  We need to subsample, because VTK cannot natively display real second-order functions
-  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order));
+  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
   vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
-  vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") +  "_" + std::to_string(neumannValues[0]) + "_0");
+  vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_0");
+
+  typedef MultiLinearGeometry<double, dim-1, dim> ML;
+  std::unordered_map<GridType::GlobalIdSet::IdType, ML> geometriesOnShellBoundary;
+  
+  auto& idSet = grid->globalIdSet();
+
+  // Read in the grid deformation
+  if (startFromFile) {
+    const std::string pathToGridDeformationFile = parameterSet.get("pathToGridDeformationFile", "");
+
+    // for this, we create a basis of order 1 in order to deform the geometries on the surface shell boundary
+    typedef Dune::Functions::LagrangeBasis<GridView, 1> FEBasisOrder1;
+    FEBasisOrder1 feBasisOrder1(gridView);
+  
+    // Read grid deformation information from the file specified in the parameter set via gridDeformationFile
+    BlockVector<FieldVector<double,3> > gridDeformationFromFile(feBasisOrder1.size());
+    std::string line;
+    std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile"));
+    if (file.is_open()) {
+      size_t i = 0;
+      while (std::getline(file, line)) {
+        size_t j = 0;
+        std::stringstream entries(line);
+        std::string entry;
+        FieldVector<double,3> coord(0);
+        while(entries >> entry) {
+          coord[j++] = std::stod(entry);
+        }
+        gridDeformationFromFile[i++] = coord;
+      }
+      if (i != feBasisOrder1.size())
+        DUNE_THROW(Exception, "Error: Grid and deformation vector do not match!");
+      file.close();
+    } else {
+      DUNE_THROW(Exception, "Error: Could not open the file containing the deformation vector!");
+    }
+
+    auto gridDeformationFromFileFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasisOrder1, gridDeformationFromFile);
+    auto localGridDeformationFromFileFunction = localFunction(gridDeformationFromFileFunction);
+
+    //Write out the stress-free geometries that were read in
+    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+    vtkWriter.addVertexData(localGridDeformationFromFileFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+    vtkWriter.write("stress-free-geometries");
+
+    //Iterate over boundary, each facet on the boundary has an element (boundaryElement.inside()) with a unique global id (idSet.subId);
+    //we store the new geometry in the map with this id as reference
+    for (auto boundaryElement : surfaceShellBoundary) {
+      localGridDeformationFromFileFunction.bind(boundaryElement.inside());
+      std::vector<Dune::FieldVector<double,dim>> corners;
+      for (int i = 0; i < boundaryElement.geometry().corners(); i++) {
+        auto corner = boundaryElement.geometry().corner(i);
+        corner += localGridDeformationFromFileFunction(boundaryElement.inside().geometry().local(boundaryElement.geometry().corner(i)));
+        corners.push_back(corner);
+      }
+      localGridDeformationFromFileFunction.unbind();
+      GridType::GlobalIdSet::IdType id = idSet.subId(boundaryElement.inside(), boundaryElement.indexInInside(), 1);
+      ML boundaryGeometry(boundaryElement.geometry().type(), corners);
+      geometriesOnShellBoundary.insert({id, boundaryGeometry});
+    }
+  } else {
+    // Read grid deformation from deformation function
+    auto gridDeformationLambda = std::string("lambda x: (") + parameterSet.get<std::string>("gridDeformation") + std::string(")");
+    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > gridDeformation(Python::evaluate(gridDeformationLambda));
+
+    //Iterate over boundary, each facet on the boundary has an element (boundaryElement.inside()) with a unique global id (idSet.subId);
+    //we store the new geometry in the map with this id as reference
+    for (auto boundaryElement : surfaceShellBoundary) {
+      std::vector<Dune::FieldVector<double,dim>> corners;
+      for (int i = 0; i < boundaryElement.geometry().corners(); i++) {
+        auto corner = gridDeformation(boundaryElement.geometry().corner(i));
+        corners.push_back(corner);
+      }
+      GridType::GlobalIdSet::IdType id = idSet.subId(boundaryElement.inside(), boundaryElement.indexInInside(), 1);
+      ML boundaryGeometry(boundaryElement.geometry().type(), corners);
+      geometriesOnShellBoundary.insert({id, boundaryGeometry});
+    }
+  }
+
+  const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+  Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
+  Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
+
+  Python::Reference pythonObject = surfaceShellCallable();
+  PythonFunction<Dune::FieldVector<double, dim>, double> fThickness(pythonObject.get("thickness"));
+  PythonFunction<Dune::FieldVector<double, dim>, Dune::FieldVector<double, 2>> fLame(pythonObject.get("lame"));
 
   for (int i=0; i<numHomotopySteps; i++)
   {
@@ -361,7 +450,6 @@ int main (int argc, char *argv[]) try
     //   Create an assembler for the energy functional
     // ////////////////////////////////////////////////////////////
 
-    const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
 
 #if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
     std::shared_ptr<NeumannFunction> neumannFunction;
@@ -465,7 +553,8 @@ int main (int argc, char *argv[]) try
       UnitVector vertexNormal(vertexNormalRaw);
       vertexNormals[i] = vertexNormal;
     }
-    surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary);
+
+    surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
 
     std::shared_ptr<LocalEnergyBase> totalEnergy;
     totalEnergy = std::make_shared<GFE::SumCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble>> (elasticAndNeumann, surfaceCosseratEnergy);
@@ -553,7 +642,7 @@ int main (int argc, char *argv[]) try
     auto localDisplacementFunction = localFunction(displacementFunction);
 
     //  We need to subsample, because VTK cannot natively display real second-order functions
-    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order));
+    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
     vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
     vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_" + std::to_string(i+1));
   }
diff --git a/src/film-on-substrate.parset b/src/film-on-substrate.parset
index 3ddb62b72f50eb14df411569550eb5f88bc64c95..645c408bf9efca144da5cdc2c0dfe1b101901442 100644
--- a/src/film-on-substrate.parset
+++ b/src/film-on-substrate.parset
@@ -3,16 +3,49 @@
 #############################################
 
 structuredGrid = true
-
-# bounding box
+# 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:
 lower = 0 0 0
-upper = 10 10 0.5
+upper = 200 100 200
 
-elements = 10 10 5
+elements = 10 5 5
 
-# Number of grid levels
+# Number of grid levels, all elements containing surfaceshell grid vertices will get adaptively refined
 numLevels = 2
 
+# When starting from a file, the stress-free configuration of the surfaceShell is read from a file, this file needs to match the *finest* grid level!
+startFromFile = true
+pathToGridDeformationFile = ./
+gridDeformationFile = deformation
+
+# When not starting from a file, deformation of the surface shell part can be given here using the gridDeformation function
+gridDeformation="[1.3*x[0], x[1], x[2]]"
+
+#############################################
+#  Boundary values
+#############################################
+
+problem = cantilever
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+###  Python predicate specifying all surfaceshell grid vertices, elements conataining these vertices will get adaptively refined
+# x is the vertex coordinate
+surfaceShellVerticesPredicate = "x[2] > 199.99"
+
+###  Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 199.99"
+
+##  Neumann values
+neumannValues = 0 0 0 
+
+# Initial deformation
+initialDeformation = "x"
+
 #############################################
 #  Solver parameters
 #############################################
@@ -64,16 +97,12 @@ energy = stvenantkirchhoff
 ## For the Wriggers L-shape example
 [materialParameters]
 
-# shell thickness
-thickness = 0.6
+surfaceShellParameters = surface-shell-parameters-1-3
 
 ## Lame parameters for stvenantkirchhoff, E = mu(3*lambda + 2*mu)/(lambda + mu)
 # mu = 2.7191e+4
 # lambda = 4.4364e+4
 
-#Lame parameters for the cosserat material
-mu_cosserat = 7.64e+6 #Parameters for E = 20.4 MPa, nu=0.335
-lambda_cosserat = 1.55e+7
 
 # Cosserat couple modulus
 mu_c = 0
@@ -111,31 +140,4 @@ mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute th
 # 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 = cantilever
-
-###  Python predicate specifying all Dirichlet grid vertices
-# x is the vertex coordinate
-dirichletVerticesPredicate = "x[0] < 0.01"
-
-###  Python predicate specifying all Neumann grid vertices
-# x is the vertex coordinate
-surfaceShellVerticesPredicate = "x[2] > 0.49"
-
-###  Python predicate specifying all Neumann grid vertices
-# x is the vertex coordinate
-neumannVerticesPredicate = "x[0] > 9.99"
-
-##  Neumann values
-neumannValues = 1e5 0 0 
-
-# Initial deformation
-initialDeformation = "x"
-
-#startFromFile = yes
-#initialIterateFilename = initial_iterate.vtu
+[]
\ No newline at end of file
diff --git a/src/surface-shell-parameters-1-3.py b/src/surface-shell-parameters-1-3.py
new file mode 100644
index 0000000000000000000000000000000000000000..17ebcee453e862d3fcfd9fd9401f0404dd68e4a0
--- /dev/null
+++ b/src/surface-shell-parameters-1-3.py
@@ -0,0 +1,18 @@
+import math
+
+class SurfaceShellParameters:
+    def thickness(self, x):
+        if (x[1] > 50.0):
+            return 1.0
+        else:
+            return 3.0
+
+    #Lame parameters for the surface shell material:
+    #Parameters for E = 20.4 MPa, nu=0.335
+    #mu_cosserat = 7.64e+6
+    #lambda_cosserat = 1.55e+7
+    def lame(self, x):
+        if (x[1] > 50.0):
+            return [7.64e+6,1.55e+7]
+        else:
+            return [7.64e+6,1.55e+7]