diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index ab5ab42d0c89907bb20cc64de4927f4cec892956..c491b6fd9c5c8d625a7c9adde21880e4aa2989b2 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -434,8 +434,6 @@ assemble(const Entity& element,
this->A[i][i][j][j] = 1;
#else
-#if 1
-
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
@@ -477,152 +475,6 @@ assemble(const Entity& element,
for (int k=0; k<blocksize; k++)
this->A[i][j][k] = localSolution[j].projectOntoTangentSpace(this->A[i][j][k]);
-
-
-
-
-#else
-
-
- double eps = 1e-4;
-
- typedef typename Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> >::row_type::iterator ColumnIterator;
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- std::vector<TargetSpace> forwardSolution = localSolution;
- std::vector<TargetSpace> backwardSolution = localSolution;
-
- std::vector<TargetSpace> forwardForwardSolution = localSolution;
- std::vector<TargetSpace> forwardBackwardSolution = localSolution;
- std::vector<TargetSpace> backwardForwardSolution = localSolution;
- std::vector<TargetSpace> backwardBackwardSolution = localSolution;
-
- // ///////////////////////////////////////////////////////////////
- // Loop over all blocks of the element matrix
- // ///////////////////////////////////////////////////////////////
- for (int i=0; i<this->A.N(); i++) {
-
- ColumnIterator cIt = this->A[i].begin();
- ColumnIterator cEndIt = this->A[i].end();
-
- for (; cIt!=cEndIt; ++cIt) {
-
- // compute only the upper right triangular matrix
- if (cIt.index() < i)
- continue;
-
- // ////////////////////////////////////////////////////////////////////////////
- // Compute a finite-difference approximation of this hessian matrix block
- // ////////////////////////////////////////////////////////////////////////////
-
- for (int j=0; j<blocksize; j++) {
-
- for (int k=0; k<blocksize; k++) {
-
- // compute only the upper right triangular matrix
- if (i==cIt.index() && k<j)
- continue;
-
- // Diagonal entries
- if (i==cIt.index() && j==k) {
-
- infinitesimalVariation(forwardSolution[i], eps, j);
- infinitesimalVariation(backwardSolution[i], -eps, j);
-
- double forwardEnergy = energy(element, forwardSolution);
-
- double solutionEnergy = energy(element, localSolution);
-
- double backwardEnergy = energy(element, backwardSolution);
-
- // Second derivative
- (*cIt)[j][k] = (forwardEnergy - 2*solutionEnergy + backwardEnergy) / (eps*eps);
-
- forwardSolution[i] = localSolution[i];
- backwardSolution[i] = localSolution[i];
-
- } else {
-
- // Off-diagonal entries
- infinitesimalVariation(forwardForwardSolution[i], eps, j);
- infinitesimalVariation(forwardForwardSolution[cIt.index()], eps, k);
- infinitesimalVariation(forwardBackwardSolution[i], eps, j);
- infinitesimalVariation(forwardBackwardSolution[cIt.index()], -eps, k);
- infinitesimalVariation(backwardForwardSolution[i], -eps, j);
- infinitesimalVariation(backwardForwardSolution[cIt.index()], eps, k);
- infinitesimalVariation(backwardBackwardSolution[i], -eps, j);
- infinitesimalVariation(backwardBackwardSolution[cIt.index()],-eps, k);
-
- double forwardForwardEnergy = energy(element, forwardForwardSolution);
-
- double forwardBackwardEnergy = energy(element, forwardBackwardSolution);
-
- double backwardForwardEnergy = energy(element, backwardForwardSolution);
-
- double backwardBackwardEnergy = energy(element, backwardBackwardSolution);
-
- (*cIt)[j][k] = (forwardForwardEnergy + backwardBackwardEnergy
- - forwardBackwardEnergy - backwardForwardEnergy) / (4*eps*eps);
-
- forwardForwardSolution[i] = localSolution[i];
- forwardForwardSolution[cIt.index()] = localSolution[cIt.index()];
- forwardBackwardSolution[i] = localSolution[i];
- forwardBackwardSolution[cIt.index()] = localSolution[cIt.index()];
- backwardForwardSolution[i] = localSolution[i];
- backwardForwardSolution[cIt.index()] = localSolution[cIt.index()];
- backwardBackwardSolution[i] = localSolution[i];
- backwardBackwardSolution[cIt.index()] = localSolution[cIt.index()];
-
- }
-
- }
-
- }
-
- }
-
- }
-
- // ///////////////////////////////////////////////////////////////
- // Symmetrize the matrix
- // This is possible expensive, but I want to be absolute sure
- // that the matrix is symmetric.
- // ///////////////////////////////////////////////////////////////
- for (int i=0; i<this->A.N(); i++) {
-
- ColumnIterator cIt = this->A[i].begin();
- ColumnIterator cEndIt = this->A[i].end();
-
- for (; cIt!=cEndIt; ++cIt) {
-
- if (cIt.index()>i)
- continue;
-
-
- if (cIt.index()==i) {
-
- for (int j=1; j<blocksize; j++)
- for (int k=0; k<j; k++)
- (*cIt)[j][k] = (*cIt)[k][j];
-
- } else {
-
- const Dune::FieldMatrix<double,blocksize,blocksize>& other = this->A[cIt.index()][i];
-
- for (int j=0; j<blocksize; j++)
- for (int k=0; k<blocksize; k++)
- (*cIt)[j][k] = other[k][j];
-
-
- }
-
-
- }
-
- }
-#endif
#endif
}