diff --git a/src/unitvector.hh b/src/unitvector.hh
index f1b8186f19c99b5d4f5a125e61ae2730e51c23ff..ca0dee985bf9c58090ff142d69eead2a7a0e4a01 100644
--- a/src/unitvector.hh
+++ b/src/unitvector.hh
@@ -55,7 +55,7 @@ public:
         }
 
         // Project gradient onto the tangent plane at b in order to obtain the surface gradient
-        result.axpy(-1*(b.data_*result), b.data_);
+        result = b.projectOntoTangentSpace(result);
 
         // Gradient must be a tangent vector at b, in other words, orthogonal to it
         assert( std::abs(b.data_ * result) < 1e-7);
@@ -63,6 +63,13 @@ public:
         return result;
     }
 
+    /** \brief Project tangent vector of R^n onto the tangent space */
+    EmbeddedTangentVector projectOntoTangentSpace(const EmbeddedTangentVector& v) const {
+        EmbeddedTangentVector result = v;
+        result.axpy(-1*(data_*result), data_);
+        return result;
+    }
+
     /** \brief The global coordinates, if you really want them */
     const Dune::FieldVector<double,dim>& globalCoordinates() const {
         return data_;