diff --git a/dune/gfe/hyperbolichalfspacepoint.hh b/dune/gfe/hyperbolichalfspacepoint.hh
index 393f776a291ccf794b556b941680bd77308e8615..c3919be70f72a4981c60d09fb6f2e2af183233b1 100644
--- a/dune/gfe/hyperbolichalfspacepoint.hh
+++ b/dune/gfe/hyperbolichalfspacepoint.hh
@@ -48,6 +48,21 @@ class HyperbolicHalfspacePoint
         }
     }
 
+    /** \brief Compute derivative of $F(p,q) = 1 + ||p-q||^2 / 2p_nq_n with respect to p 
+     \param[in] diffNormSquared Expected to be ||p-q||^2, taken from the caller for efficiency reasons
+     */
+    static Dune::FieldVector<T,N> computeDFdp(const HyperbolicHalfspacePoint& p, const HyperbolicHalfspacePoint& q, const T& diffNormSquared)
+    {
+        Dune::FieldVector<T,N> result;
+        
+        for (size_t i=0; i<N-1; i++)
+            result[i] = ( q.data_[i] - p.data_[i] ) / (q.data_[N-1] * p.data_[N-1]);
+        
+        result[N-1] = - diffNormSquared / (2*p.data_[N-1]*p.data_[N-1]*q.data_[N-1]) + (p.data_[N-1] - q.data_[N-1]) / (p.data_[N-1]*q.data_[N-1]);
+        
+        return result;
+    }
+    
     /** \brief Compute derivative of $F(p,q) = 1 + ||p-q||^2 / 2p_nq_n with respect to q 
      \param[in] diffNormSquared Expected to be ||p-q||^2, taken from the caller for efficiency reasons
      */
@@ -240,12 +255,7 @@ public:
         T diffNormSquared = (p-q).two_norm2();
 
         // Compute first derivatives of F with respect to p and q
-        Dune::FieldVector<T,N> dFdp;
-        for (size_t i=0; i<N-1; i++)
-            dFdp[i] = ( p[i] - q[i] ) / (p[N-1] * q[N-1]);
-        
-        dFdp[N-1] = - diffNormSquared / (2*p[N-1]*q[N-1]*q[N-1]) - (p[N-1] - q[N-1]) / (p[N-1]*q[N-1]);
-
+        Dune::FieldVector<T,N> dFdp = computeDFdp(a,b,diffNormSquared);
         Dune::FieldVector<T,N> dFdq = computeDFdq(a,b,diffNormSquared);
 
         // Compute second derivatives of F
@@ -420,12 +430,7 @@ public:
         T diffNormSquared = (p-q).two_norm2();
 
         // Compute first derivatives of F with respect to p and q
-        Dune::FieldVector<T,N> dFdp;
-        for (size_t i=0; i<N-1; i++)
-            dFdp[i] = ( p[i] - q[i] ) / (p[N-1] * q[N-1]);
-        
-        dFdp[N-1] = - diffNormSquared / (2*p[N-1]*q[N-1]*q[N-1]) - (p[N-1] - q[N-1]) / (p[N-1]*q[N-1]);
-
+        Dune::FieldVector<T,N> dFdp = computeDFdp(a,b,diffNormSquared);
         Dune::FieldVector<T,N> dFdq = computeDFdq(a,b,diffNormSquared);
 
         // Compute second derivatives of F