diff --git a/src/harmonicenergystiffness.hh b/src/harmonicenergystiffness.hh
index 42ea57fdb24aa3062e9d8a847b846d093700c1e5..c0306a84785cfbdd17fc4c68cf22fb04967a9cd6 100644
--- a/src/harmonicenergystiffness.hh
+++ b/src/harmonicenergystiffness.hh
@@ -69,7 +69,7 @@ energy(const Entity& element,
         double weight = quad[pt].weight() * integrationElement;
 
         // The derivative of the local function defined on the reference element
-        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
+        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivativeFD(quadPos);
 
         // The derivative of the function defined on the actual element
         Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> derivative(0);
@@ -77,6 +77,25 @@ energy(const Entity& element,
         for (int comp=0; comp<4; comp++)
             jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
 
+        TargetSpace value = localGeodesicFEFunction.evaluate(quadPos);
+
+        for (int i=0; i<gridDim; i++) {
+            double dotproduct = 0;
+            for (int j=0; j<4; j++)
+                dotproduct += value[j]*derivative[j][i];
+            assert(std::fabs(dotproduct) < 1e-6);
+        }
+
+#if 0
+        double derivativenorm0 = 0;
+        double derivativenorm1 = 0;
+        for (int i=0; i<4; i++) {
+            derivativenorm0 += derivative[i][0]*derivative[i][0];
+            derivativenorm1 += derivative[i][1]*derivative[i][1];
+        }
+        std::cout << "Derivative norm: " << derivativenorm0 << ",  " << derivativenorm1 << std::endl;
+#endif
+
         for (int comp=0; comp<4; comp++) {
 
             for (int dir=0; dir<gridDim; dir++)
@@ -85,7 +104,7 @@ energy(const Entity& element,
         }
 
     }
-    
+
     return energy;
 }