diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index 9a13a6826ccc98a99033e830c5c4b618d0a9ac2a..3e072407e0d5bb42abd24d555fb1c8a4a915db3f 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -73,6 +73,25 @@ public:
             matrix.axpy(weights_[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
     }
 
+    void assembleHessian(const TargetSpace& x,
+                         Dune::FieldMatrix<double,size,size>& matrix) const
+    {
+        Dune::FieldMatrix<double,embeddedSize,embeddedSize> embeddedHessian;
+        assembleEmbeddedHessian(x,embeddedHessian);
+        
+        Dune::FieldMatrix<double,size,embeddedSize> frame = x.orthonormalFrame();
+        
+        // this is frame * embeddedHessian * frame^T
+        for (int i=0; i<size; i++)
+            for (int j=0; j<size; j++) {
+                matrix[i][j] = 0;
+                for (int k=0; k<embeddedSize; k++)
+                    for (int l=0; l<embeddedSize; l++)
+                        matrix[i][j] += frame[i][k]*embeddedHessian[k][l]*frame[j][l];
+            }
+        
+    }
+
     const std::vector<TargetSpace> coefficients_;
 
     const std::vector<double> weights_;