diff --git a/src/gsMaterialMatrixBase.h b/src/gsMaterialMatrixBase.h
index 6f2147e..cef3f66 100644
--- a/src/gsMaterialMatrixBase.h
+++ b/src/gsMaterialMatrixBase.h
@@ -196,6 +196,7 @@ class gsMaterialMatrixBase
      * @param[in]  u       The in-plane shell coordinates to be eveluated on
      * @param      result  The result
      */
+    GISMO_DEPRECATED
     virtual void    pstretch_into(const index_t /*patch*/, const gsMatrix<T>& /*u*/, gsMatrix<T>& /*result*/) const
     { GISMO_NO_IMPLEMENTATION; }
     /**
@@ -205,6 +206,7 @@ class gsMaterialMatrixBase
      * @param[in]  u       The in-plane shell coordinates to be eveluated on
      * @param      result  The result
      */
+    GISMO_DEPRECATED
     virtual void    pstretchDir_into(const index_t /*patch*/, const gsMatrix<T>& /*u*/, gsMatrix<T>& /*result*/) const
     { GISMO_NO_IMPLEMENTATION; }
         /**
@@ -214,6 +216,7 @@ class gsMaterialMatrixBase
      * @param[in]  u       The in-plane shell coordinates to be eveluated on
      * @param      result  The result
      */
+    GISMO_DEPRECATED
     virtual void    pstress_into(const index_t /*patch*/, const gsMatrix<T>& /*u*/, gsMatrix<T>& /*result*/) const
     { GISMO_NO_IMPLEMENTATION; }
     /**
@@ -223,6 +226,7 @@ class gsMaterialMatrixBase
      * @param[in]  u       The in-plane shell coordinates to be eveluated on
      * @param      result  The result
      */
+    GISMO_DEPRECATED
     virtual void    pstressDir_into(const index_t /*patch*/, const gsMatrix<T>& /*u*/, gsMatrix<T>& /*result*/) const
     { GISMO_NO_IMPLEMENTATION; }
     /**
@@ -242,6 +246,7 @@ class gsMaterialMatrixBase
      * @param[in]  u       The in-plane shell coordinates to be eveluated on
      * @param      result  The result
      */
+    GISMO_DEPRECATED
     virtual void  parameters_into(const index_t /*patch*/, const gsMatrix<T>& /*u*/, gsMatrix<T>& /*result*/) const
     { GISMO_NO_IMPLEMENTATION; }
 
@@ -252,6 +257,7 @@ class gsMaterialMatrixBase
      * @param[in]  u       The in-plane shell coordinates to be eveluated on
      * @param      result  The result
      */
+    GISMO_DEPRECATED
     virtual void  transform_into(const index_t /*patch*/, const gsMatrix<T>& /*u*/, gsMatrix<T>& /*result*/) const
     { GISMO_NO_IMPLEMENTATION; }
 
diff --git a/src/gsMaterialMatrixBaseDim.hpp b/src/gsMaterialMatrixBaseDim.hpp
index eb8d3cb..8ca1bfc 100644
--- a/src/gsMaterialMatrixBaseDim.hpp
+++ b/src/gsMaterialMatrixBaseDim.hpp
@@ -178,7 +178,6 @@ gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_deformation(const index_t pat
     this->_computePoints(patch,u);
 
     gsMatrix<T> result(9, u.cols() * z.rows());
-    std::pair<gsVector<T>,gsMatrix<T>> res;
     index_t colIdx;
     for (index_t k=0; k!=u.cols(); k++)
     {
diff --git a/src/gsMaterialMatrixNonlinear.h b/src/gsMaterialMatrixNonlinear.h
index c57569e..41d3238 100644
--- a/src/gsMaterialMatrixNonlinear.h
+++ b/src/gsMaterialMatrixNonlinear.h
@@ -159,6 +159,9 @@ class gsMaterialMatrixNonlinear :    public gsMaterialMatrixBaseDim<dim,T>
     /// See \ref gsMaterialMatrixBase for details
     gsMatrix<T> eval3D_CauchyVector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput out = MaterialOutput::Generic) const override;
 
+    /// See \ref gsMaterialMatrixBase for details
+    gsMatrix<T> eval3D_pstretch(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const override;
+
     /// See \ref gsMaterialMatrixBase for details
     gsMatrix<T> eval3D_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput out = MaterialOutput::Generic) const override;
 
@@ -569,6 +572,14 @@ class gsMaterialMatrixNonlinear :    public gsMaterialMatrixBaseDim<dim,T>
     template<enum Material _mat, bool _com>
     typename std::enable_if<!_com && !(_mat==Material::SvK), gsMatrix<T>>::type _eval3D_pstress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const;
 
+    // Stretch overrides
+    template<bool _com>
+    typename std::enable_if<_com, gsMatrix<T>>::type
+    eval3D_pstretch_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const;
+    template<bool _com>
+    typename std::enable_if<!_com, gsMatrix<T>>::type
+    eval3D_pstretch_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const;
+
     template<enum Material _mat>
     typename std::enable_if<_mat==Material::MR, void>::type _setRatio_impl(const gsFunctionSet<T> & Ratio);
     template<enum Material _mat>
diff --git a/src/gsMaterialMatrixNonlinear.hpp b/src/gsMaterialMatrixNonlinear.hpp
index 34ba047..488a24c 100644
--- a/src/gsMaterialMatrixNonlinear.hpp
+++ b/src/gsMaterialMatrixNonlinear.hpp
@@ -563,6 +563,56 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_CauchyVe
     return this->eval3D_CauchyStress(patch,u,z,out);
 }
 
+template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
+gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstretch(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
+{
+    return eval3D_pstretch_impl<comp>(patch,u,z);
+}
+
+template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
+template <bool _comp>
+typename std::enable_if< _comp, gsMatrix<T>>::type
+gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstretch_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
+{
+    this->_computePoints(patch,u);
+
+    gsMatrix<T> result(3, u.cols() * z.rows());
+    std::pair<gsVector<T>,gsMatrix<T>> res;
+    index_t colIdx;
+
+    gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,z);
+    T C33;
+
+    gsMatrix<T,3,3> c;
+    for (index_t k=0; k!= u.cols(); k++)
+    {
+        for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
+        {
+            colIdx = j*u.cols()+k;
+            this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
+
+            C33 = C33s(0,k);
+
+            // Compute c
+            c.setZero();
+            c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
+            c(2,2) = C33; // c33
+
+            res = this->_evalStretch(c,m_data.mine().m_gcon_ori);
+            result.col(colIdx) = res.first;
+        }
+    }
+    return result;
+}
+
+template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
+template <bool _comp>
+typename std::enable_if<!_comp, gsMatrix<T>>::type
+gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstretch_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
+{
+    return Base::eval3D_pstretch(patch,u,z);
+}
+
 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
 gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
 {
diff --git a/src/gsThinShellUtils.h b/src/gsThinShellUtils.h
index 57ede7f..eec8b42 100644
--- a/src/gsThinShellUtils.h
+++ b/src/gsThinShellUtils.h
@@ -26,7 +26,7 @@
 #include <gsUtils/gsPointGrid.h>
 
 #include <gsAssembler/gsAssemblerOptions.h>
-#include <gsAssembler/gsExpressions.h>
+#include <gsExpressions/gsExpressions.h>
 
 
 #  define MatExprType  auto