diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 1d9bb430fc1dfdc7df2662a53c1453134059bdf8..c9f11e4b865b031f5e22aa5352f1b16faf6f8b71 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -892,6 +892,27 @@ public:
 
     }
 
+    /** \brief Derivative of the map from unit quaternions to orthogonal matrices
+     */
+    static Tensor3<T,3,3,4> derivativeOfQuaternionToMatrix(const Dune::FieldVector<T,4>& p)
+    {
+      Tensor3<T,3,3,4> result;
+
+      result[0][0] = { 2*p[0], -2*p[1], -2*p[2],  2*p[3]};
+      result[0][1] = { 2*p[1],  2*p[0], -2*p[3], -2*p[2]};
+      result[0][2] = { 2*p[2],  2*p[3],  2*p[0],  2*p[1]};
+
+      result[1][0] = { 2*p[1],  2*p[0],  2*p[3],  2*p[2]};
+      result[1][1] = {-2*p[0],  2*p[1], -2*p[2],  2*p[3]};
+      result[1][2] = {-2*p[3],  2*p[2],  2*p[1], -2*p[0]};
+
+      result[2][0] = { 2*p[2], -2*p[3],  2*p[0], -2*p[1]};
+      result[2][1] = { 2*p[3],  2*p[2],  2*p[1],  2*p[0]};
+      result[2][2] = {-2*p[0], -2*p[1],  2*p[2],  2*p[3]};
+
+      return result;
+    }
+
     /** \brief Set rotation from orthogonal matrix
 
     We tacitly assume that the matrix really is orthogonal */