diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index 3c23d23285b334c81971459c86a6cbd50b3af045..be41bba0819d9223c4c5b586db2ddbc9da4b176f 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -81,20 +81,15 @@ public:
     void assembleHessian(const TargetSpace& x,
                          Dune::FieldMatrix<ctype,size,size>& matrix) const
     {
-        Dune::FieldMatrix<ctype,embeddedSize,embeddedSize> embeddedHessian;
+        Dune::SymmetricMatrix<ctype,embeddedSize> embeddedHessian;
         assembleEmbeddedHessian(x,embeddedHessian);
 
         Dune::FieldMatrix<ctype,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];
-            }
-
+            for (int j=0; j<size; j++)
+                matrix[i][j] = embeddedHessian.energyScalarProduct(frame[i], frame[j]);
     }
 
     const std::vector<TargetSpace> coefficients_;