diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index df4749eb9be221a20381c136f34bfa3218e99013..1aca4d04ea179e7efe604f6bc2e942c13c58435c 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -3,10 +3,16 @@
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/parametertree.hh>
+#include <dune/common/version.hh>
+
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/fufem/boundarypatch.hh>
 
+#if DUNE_VERSION_LT(DUNE_COMMON, 2, 9)
+#include <dune/matrix-vector/transpose.hh>
+#endif
+
 #include <dune/gfe/localenergy.hh>
 #include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
 #ifdef PROJECTED_INTERPOLATION
@@ -267,9 +273,17 @@ public:
 
         Dune::FieldMatrix<field_type,3,3> wryness(0);
 
+#if DUNE_VERSION_LT(DUNE_COMMON, 2, 9)
+        Dune::FieldMatrix<field_type,dim,dim> transposeR;
+        Dune::MatrixVector::transpose(R,transposeR);
+        auto axialVectorx1 = SkewMatrix<field_type,3>(transposeR*dRx1).axial();
+        auto axialVectorx2 = SkewMatrix<field_type,3>(transposeR*dRx2).axial();
+        auto axialVectorx3 = SkewMatrix<field_type,3>(transposeR*dRx3).axial();
+#else
         auto axialVectorx1 = SkewMatrix<field_type,3>(transpose(R)*dRx1).axial();
         auto axialVectorx2 = SkewMatrix<field_type,3>(transpose(R)*dRx2).axial();
         auto axialVectorx3 = SkewMatrix<field_type,3>(transpose(R)*dRx3).axial();
+#endif
         for (int i=0; i<3; i++) {
             wryness[i][0] = axialVectorx1[i];
             wryness[i][1] = axialVectorx2[i];