diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index 3bbbcdd82bb6cda4cc504c8656aad5244c9f7b7e..819ac2f29863cd5b0957538b9efde6eb4dbe9ba4 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/localgeodesicfeadolcstiffness.hh
@@ -244,9 +244,9 @@ assembleGradientAndHessian(const Entity& element,
     for (size_t i=0; i<nDofs; i++)
         orthonormalFrame[i] = localSolution[i].orthonormalFrame();
 
-#ifndef ADOLC_VECTOR_MODE
     Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
 
+#ifndef ADOLC_VECTOR_MODE
     std::vector<double> v(nDoubles);
     std::vector<double> w(nDoubles);
 
@@ -270,42 +270,9 @@ assembleGradientAndHessian(const Entity& element,
         // Make v the null vector again
         std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
       }
-
-    // From this, compute the Hessian with respect to the manifold (which we assume here is embedded
-    // isometrically in a Euclidean space.
-    // For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
-    // at the Riemannian Hessian".
-
-    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
-    typedef typename TargetSpace::TangentVector         TangentVector;
-
-    this->A_.setSize(nDofs,nDofs);
-
-    for (size_t col=0; col<nDofs; col++) {
-
-        for (size_t subCol=0; subCol<blocksize; subCol++) {
-
-            EmbeddedTangentVector z = orthonormalFrame[col][subCol];
-
-            // P_x \partial^2 f z
-            for (size_t row=0; row<nDofs; row++) {
-                TangentVector semiEmbeddedProduct;
-                embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
-
-                for (int subRow=0; subRow<blocksize; subRow++)
-                    this->A_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
-            }
-
-        }
-
-    }
-
 #else
-    Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
-
     int n = nDoubles;
-    //std::cout << "n: " << n << std::endl;
-    int nDirections = nDoubles;
+    int nDirections = nDofs * blocksize;
     double* tangent[nDoubles];
     for(size_t i=0; i<nDoubles; i++)
         tangent[i] = (double*)malloc(nDirections*sizeof(double));
@@ -314,26 +281,27 @@ assembleGradientAndHessian(const Entity& element,
     for(size_t i=0; i<nDoubles; i++)
         rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
 
-    //std::cout << "nDirections: " << nDirections << std::endl;
-    for (int i=0; i<n; i++)
-      for (int j=0; j<nDirections; j++)
-        tangent[i][j] = (i==j);
+    for (int j=0; j<nDirections; j++)
+    {
+      for (int i=0; i<n; i++)
+        tangent[i][j] = 0.0;
 
-    hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+      for (int i=0; i<embeddedBlocksize; i++)
+        tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
+    }
 
-    //std::cout << "Hallo Welt!" << std::endl;
+    hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
 
     // Copy Hessian into Dune data type
     for(size_t i=0; i<nDoubles; i++)
-      for (size_t j=0; j<nDoubles; j++)
-        embeddedHessian[i/embeddedBlocksize][j/embeddedBlocksize][i%embeddedBlocksize][j%embeddedBlocksize] = rawHessian[i][j];
-
+      for (size_t j=0; j<nDirections; j++)
+        embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
 
     for(size_t i=0; i<nDoubles; i++) {
         free(rawHessian[i]);
         free(tangent[i]);
     }
-    //printmatrix(std::cout, embeddedHessian, "hessian", "--");
+#endif
 
     // From this, compute the Hessian with respect to the manifold (which we assume here is embedded
     // isometrically in a Euclidean space.
@@ -352,22 +320,17 @@ assembleGradientAndHessian(const Entity& element,
             EmbeddedTangentVector z = orthonormalFrame[col][subCol];
 
             // P_x \partial^2 f z
-            std::vector<EmbeddedTangentVector> semiEmbeddedProduct(nDofs);
-
             for (size_t row=0; row<nDofs; row++) {
-                embeddedHessian[row][col].mv(z,semiEmbeddedProduct[row]);
-                //tmp1 = localSolution[row].projectOntoTangentSpace(tmp1);
-                TangentVector tmp2;
-                orthonormalFrame[row].mv(semiEmbeddedProduct[row],tmp2);
+                TangentVector semiEmbeddedProduct;
+                embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
 
                 for (int subRow=0; subRow<blocksize; subRow++)
-                    this->A_[row][col][subRow][subCol] = tmp2[subRow];
+                    this->A_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
             }
 
         }
 
     }
-#endif
 
     //////////////////////////////////////////////////////////////////////////////////////
     //  Further correction due to non-planar configuration space