diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index 3f31724ca165fdfb387c9119ce4ebee94b94fd95..ce516ace831adcd4cfe47103c98ff78e49a521e2 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -238,124 +238,6 @@ public:
     b3_ = parameters.template get<double>("b3");
   }
 
-    /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
-   * the first equation of (4.4) in Neff's paper
-   */
-  RT quadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,dim>& U) const
-  {
-      Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
-      for (int i=0; i<dim; i++)
-          UMinus1[i][i] -= 1;
-
-      return mu_ * sym(UMinus1).frobenius_norm2()
-              + mu_c_ * skew(UMinus1).frobenius_norm2()
-              + (mu_*lambda_)/(2*mu_ + lambda_) * traceSquared(sym(UMinus1));
-  }
-
-  /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
-   * the second equation of (4.4) in Neff's paper
-   */
-  RT longQuadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,dim>& U) const
-  {
-      RT result = 0;
-
-      // shear-stretch energy
-      Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
-      for (int i=0; i<dim-1; i++)
-          for (int j=0; j<dim-1; j++)
-              sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
-
-      result += mu_ * sym2x2.frobenius_norm2();
-
-      // first order drill energy
-      Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
-      for (int i=0; i<dim-1; i++)
-          for (int j=0; j<dim-1; j++)
-              skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
-
-      result += mu_c_ * skew2x2.frobenius_norm2();
-
-
-      // classical transverse shear energy
-      result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
-
-      // elongational stretch energy
-      result += mu_*lambda_ / (2*mu_ + lambda_) * traceSquared(sym2x2);
-
-      return result;
-  }
-
-  /** \brief Energy for large-deformation problems (private communication by Patrizio Neff)
-   */
-  RT nonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,dim>& U) const
-  {
-      Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
-      for (int i=0; i<dim; i++)
-          UMinus1[i][i] -= 1;
-
-      RT detU = U.determinant();
-
-      return mu_ * sym(UMinus1).frobenius_norm2() + mu_c_ * skew(UMinus1).frobenius_norm2()
-              + (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
-  }
-
-  /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
-   * the second equation of (4.4) in Neff's paper
-   */
-  RT longNonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,dim,dim-1>& U) const
-  {
-      RT result = 0;
-
-      // shear-stretch energy
-      Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
-      for (int i=0; i<dim-1; i++)
-          for (int j=0; j<dim-1; j++)
-              sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
-
-      result += mu_ * sym2x2.frobenius_norm2();
-
-      // first order drill energy
-      Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
-      for (int i=0; i<dim-1; i++)
-          for (int j=0; j<dim-1; j++)
-              skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
-
-      result += mu_c_ * skew2x2.frobenius_norm2();
-
-
-      // classical transverse shear energy
-      result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
-
-      // elongational stretch energy
-      RT detU = U.determinant();
-      result += (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
-
-      return result;
-  }
-
-  RT curvatureEnergy(const Tensor3<field_type,3,3,dim-1>& DR) const
-  {
-#ifdef DONT_USE_CURL
-      return mu_ * std::pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0);
-#else
-      return mu_ * std::pow(L_c_ * L_c_ * curl(DR).frobenius_norm2(),q_/2.0);
-#endif
-  }
-
-  RT bendingEnergy(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,dim-1>& DR) const
-  {
-      // left-multiply the derivative of the third director (in DR[][2][]) with R^T
-      Dune::FieldMatrix<field_type,3,3> RT_DR3(0);
-      for (int i=0; i<3; i++)
-          for (int j=0; j<dim; j++)
-              for (int k=0; k<3; k++)
-                  RT_DR3[i][j] += R[k][i] * DR[k][2][j];
-
-      return mu_ * sym(RT_DR3).frobenius_norm2()
-             + mu_c_ * skew(RT_DR3).frobenius_norm2()
-             + mu_*lambda_/(2*mu_+lambda_) * traceSquared(RT_DR3);
-  }
-
 RT energy(const typename Basis::LocalView& localView,
        const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
 {