diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index 4901f8d5f6d1aaf08437c28ce1ca716d36496fac..bb019ef9bb8b48a811f7495d30694c50ca9e0175 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -32,8 +32,7 @@
  * \tparam field_type                  The coordinate type of the TargetSpace
  * \tparam StressFreeStateGridFunction Type of the GridFunction representing the Cosserat shell in a stress free state
  */
-template<class Basis, int dim, class field_type=double, class StressFreeStateGridFunction = 
-  Dune::Functions::DiscreteGlobalBasisFunction<Basis,std::vector<Dune::FieldVector<double, Basis::GridView::dimensionworld>> > >
+template<class Basis, int dim, class field_type, class StressFreeStateGridFunction>
 class NonplanarCosseratShellEnergy
   : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >,
     public MixedLocalGeodesicFEStiffness<Basis,
@@ -161,7 +160,7 @@ public:
              Alternative derivation of the higher-order constitudtive model for six-parameter elastic shells, equations (119) and (126). */
   bool useAlternativeEnergyWCoss_;
 
-  /** \brief The geometry used for assembling */
+  /** \brief The geometry of the reference deformation used for assembling */
   const StressFreeStateGridFunction* stressFreeStateGridFunction_;
 
   /** \brief The Neumann boundary */
@@ -189,23 +188,17 @@ energy(const typename Basis::LocalView& localView,
   const auto& localFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
 
 #if HAVE_DUNE_CURVEDGEOMETRY
-  // Construct a curved geometry of this element of the Cosserat shell in stress-free state
-  // When using element.geometry(), then the curvatures on the element are zero, when using a curved geometry, they are not
-  // If a parametrization representing the Cosserat shell in a stress-free state is given,
-  // this is used for the curved geometry approximation.
+  // Construct a curved geometry of this element of the Cosserat shell in its stress-free state
   // The variable local holds the local coordinates in the reference element
   // and localGeometry.global maps them to the world coordinates
-  auto curvedGeometryGridFunctionOrder = localFiniteElement.localBasis().order();
   Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> geometry(referenceElement(element),
     [this,element](const auto& local) {
-      if (not stressFreeStateGridFunction_) {
-        return element.geometry().global(local);
-      }
       auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
       localGridFunction.bind(element);
       return localGridFunction(local);
-    }, curvedGeometryGridFunctionOrder);
+    }, stressFreeStateGridFunction_->order());
 #else
+  // When using element.geometry(), the geometry of the element is flat
   auto geometry = element.geometry();
 #endif
 
@@ -436,23 +429,17 @@ energy(const typename Basis::LocalView& localView,
   const auto& orientationLocalFiniteElement = LocalFiniteElementFactory<Basis,1>::get(localView,_1);
 
 #if HAVE_DUNE_CURVEDGEOMETRY
-  // Construct a curved geometry of this element of the Cosserat shell in stress-free state
-  // When using element.geometry(), then the curvatures on the element are zero, when using a curved geometry, they are not
-  // If a parametrization representing the Cosserat shell in a stress-free state is given,
-  // this is used for the curved geometry approximation.
+  // Construct a curved geometry of this element of the Cosserat shell in its stress-free state
   // The variable local holds the local coordinates in the reference element
   // and localGeometry.global maps them to the world coordinates
-  auto curvedGeometryGridFunctionOrder = deformationLocalFiniteElement.localBasis().order();
   Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> geometry(referenceElement(element),
     [this,element](const auto& local) {
-      if (not stressFreeStateGridFunction_) {
-        return element.geometry().global(local);
-      }
       auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
       localGridFunction.bind(element);
       return localGridFunction(local);
-    }, curvedGeometryGridFunctionOrder);
+    }, stressFreeStateGridFunction_->order());
 #else
+  // When using element.geometry(), the geometry of the element is flat
   auto geometry = element.geometry();
 #endif
 
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 13b000bfc85c1c03d20c7aab4f4cf3a80fad90ce..8c8a94e5fbcdfa626ac75fb46edfb4e2e16ddfab 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -580,16 +580,15 @@ int main (int argc, char *argv[]) try
             using StiffnessType = MixedLocalGFEADOLCStiffness<CompositeBasis, RealTuple<double,3>, Rotation<double,3>>;
             std::shared_ptr<StiffnessType> localGFEStiffness;
 
-            NonplanarCosseratShellEnergy< CompositeBasis, 3, adouble,
-            Dune::Functions::DiscreteGlobalBasisFunction< DeformationFEBasis,std::vector<Dune::FieldVector<double, dimworld>> >>
-                                                                localCosseratEnergyPlanar(materialParameters,
-                                                                                                    nullptr,
-                                                                                                    &neumannBoundary,
-                                                                                                    neumannFunction,
-                                                                                                    volumeLoad);
-
-              localGFEStiffness = std::make_shared<StiffnessType>(&localCosseratEnergyPlanar, adolcScalarMode);
-
+#if HAVE_DUNE_CURVEDGEOMETRY && WORLD_DIM == 3 && GRID_DIM == 2
+            NonplanarCosseratShellEnergy<CompositeBasis, 3, adouble, decltype(creator)> localCosseratEnergy(materialParameters,
+                                                                                        &creator,
+                                                                                        &neumannBoundary,
+                                                                                        neumannFunction,
+                                                                                        volumeLoad);
+
+            localGFEStiffness = std::make_shared<StiffnessType>(&localCosseratEnergy, adolcScalarMode);
+#endif
             MixedGFEAssembler<CompositeBasis,
                       RealTuple<double,3>,
                       Rotation<double,3> > mixedAssembler(compositeBasis, localGFEStiffness.get());