diff --git a/dune/gfe/Makefile.am b/dune/gfe/Makefile.am
index 97b5f27c82e35ced58b0c02085639a558b2594c9..4a7d5d05a5ee69280196b8fce723a99b523464e5 100644
--- a/dune/gfe/Makefile.am
+++ b/dune/gfe/Makefile.am
@@ -20,6 +20,7 @@ srcinclude_HEADERS = adolcnamespaceinjections.hh \
linearalgebra.hh \
localgeodesicfefunction.hh \
localgeodesicfestiffness.hh \
+ localgeodesicfefdstiffness.hh \
localgfetestfunctionbasis.hh \
localprojectedfefunction.hh \
maxnormtrustregion.hh \
diff --git a/dune/gfe/localgeodesicfefdstiffness.hh b/dune/gfe/localgeodesicfefdstiffness.hh
new file mode 100644
index 0000000000000000000000000000000000000000..178c758369dc6f03c2f3f422b911353859b928b2
--- /dev/null
+++ b/dune/gfe/localgeodesicfefdstiffness.hh
@@ -0,0 +1,231 @@
+#ifndef DUNE_GFE_LOCAL_GEODESIC_FE_FD_STIFFNESS_HH
+#define DUNE_GFE_LOCAL_GEODESIC_FE_FD_STIFFNESS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dune/gfe/localgeodesicfestiffness.hh>
+
+/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
+ */
+template<class GridView, class LocalFiniteElement, class TargetSpace>
+class LocalGeodesicFEFDStiffness
+ : public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,TargetSpace>
+{
+ // grid types
+ typedef typename GridView::Grid::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+ typedef typename GridView::template Codim<0>::Entity Entity;
+
+ // some other sizes
+ enum {gridDim=GridView::dimension};
+
+public:
+
+ //! Dimension of a tangent space
+ enum { blocksize = TargetSpace::TangentVector::dimension };
+
+ //! Dimension of the embedding space
+ enum { embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension };
+
+ LocalGeodesicFEFDStiffness(const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, TargetSpace>* energy)
+ : localEnergy_(energy)
+ {}
+
+ /** \brief Compute the energy at the current configuration */
+ virtual RT energy (const Entity& element,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution) const
+ {
+ return localEnergy_->energy(element,localFiniteElement,localSolution);
+ }
+
+ /** \brief Assemble the element gradient of the energy functional
+
+ The default implementation in this class uses a finite difference approximation */
+ virtual void assembleGradient(const Entity& element,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& solution,
+ std::vector<typename TargetSpace::TangentVector>& gradient) const;
+
+ /** \brief Assemble the local stiffness matrix at the current position
+
+ This default implementation used finite-difference approximations to compute the second derivatives
+
+ The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
+ 'Optimization algorithms on matrix manifolds', page 107. There it says that
+ \f[
+ \langle Hess f(x)[\xi], \eta \rangle
+ = \frac 12 \frac{d^2}{dt^2} \Big(f(\exp_x(t(\xi + \eta))) - f(\exp_x(t\xi)) - f(\exp_x(t\eta))\Big)\Big|_{t=0}.
+ \f]
+ We compute that using a finite difference approximation.
+
+ */
+ virtual void assembleGradientAndHessian(const Entity& e,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient);
+
+
+ const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, TargetSpace>* localEnergy_;
+
+};
+
+template <class GridView, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEFDStiffness<GridView, LocalFiniteElement, TargetSpace>::
+assembleGradient(const Entity& element,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient) const
+{
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+
+ double eps = 1e-6;
+
+ localGradient.resize(localSolution.size());
+
+ std::vector<TargetSpace> forwardSolution = localSolution;
+ std::vector<TargetSpace> backwardSolution = localSolution;
+
+ for (size_t i=0; i<localSolution.size(); i++) {
+
+ // basis vectors of the tangent space of the i-th entry of localSolution
+ const Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
+
+ for (int j=0; j<blocksize; j++) {
+
+ typename TargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
+ forwardCorrection *= eps;
+
+ typename TargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
+ backwardCorrection *= -eps;
+
+ forwardSolution[i] = TargetSpace::exp(localSolution[i], forwardCorrection);
+ backwardSolution[i] = TargetSpace::exp(localSolution[i], backwardCorrection);
+
+ localGradient[i][j] = (energy(element,localFiniteElement,forwardSolution) - energy(element,localFiniteElement, backwardSolution)) / (2*eps);
+
+ }
+
+ forwardSolution[i] = localSolution[i];
+ backwardSolution[i] = localSolution[i];
+
+ }
+
+}
+
+
+// ///////////////////////////////////////////////////////////
+// Compute gradient by finite-difference approximation
+// ///////////////////////////////////////////////////////////
+template <class GridType, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEFDStiffness<GridType, LocalFiniteElement, TargetSpace>::
+assembleGradientAndHessian(const Entity& element,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient)
+{
+ // Number of degrees of freedom for this element
+ size_t nDofs = localSolution.size();
+
+ // Clear assemble data
+ this->A_.setSize(nDofs, nDofs);
+
+ this->A_ = 0;
+
+ const double eps = 5e-4;
+
+ std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
+ for (size_t i=0; i<B.size(); i++)
+ B[i] = localSolution[i].orthonormalFrame();
+
+ // Precompute negative energy at the current configuration
+ // (negative because that is how we need it as part of the 2nd-order fd formula)
+ RT centerValue = -energy(element, localFiniteElement, localSolution);
+
+ // Precompute energy infinitesimal corrections in the directions of the local basis vectors
+ std::vector<Dune::array<RT,blocksize> > forwardEnergy(nDofs);
+ std::vector<Dune::array<RT,blocksize> > backwardEnergy(nDofs);
+
+ //#pragma omp parallel for schedule (dynamic)
+ for (size_t i=0; i<localSolution.size(); i++) {
+ for (size_t i2=0; i2<blocksize; i2++) {
+ typename TargetSpace::EmbeddedTangentVector epsXi = B[i][i2];
+ epsXi *= eps;
+ typename TargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
+ minusEpsXi *= -1;
+
+ std::vector<TargetSpace> forwardSolution = localSolution;
+ std::vector<TargetSpace> backwardSolution = localSolution;
+
+ forwardSolution[i] = TargetSpace::exp(localSolution[i],epsXi);
+ backwardSolution[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
+
+ forwardEnergy[i][i2] = energy(element, localFiniteElement, forwardSolution);
+ backwardEnergy[i][i2] = energy(element, localFiniteElement, backwardSolution);
+
+ }
+
+ }
+
+ //////////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ //////////////////////////////////////////////////////////////
+
+ localGradient.resize(localSolution.size());
+
+ for (size_t i=0; i<localSolution.size(); i++)
+ for (int j=0; j<blocksize; j++)
+ localGradient[i][j] = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
+
+
+ ///////////////////////////////////////////////////////////////////////////
+ // Compute Riemannian Hesse matrix by finite-difference approximation.
+ // We loop over the lower left triangular half of the matrix.
+ // The other half follows from symmetry.
+ ///////////////////////////////////////////////////////////////////////////
+ //#pragma omp parallel for schedule (dynamic)
+ for (size_t i=0; i<localSolution.size(); i++) {
+ for (size_t i2=0; i2<blocksize; i2++) {
+ for (size_t j=0; j<=i; j++) {
+ for (size_t j2=0; j2<((i==j) ? i2+1 : blocksize); j2++) {
+
+ std::vector<TargetSpace> forwardSolutionXiEta = localSolution;
+ std::vector<TargetSpace> backwardSolutionXiEta = localSolution;
+
+ typename TargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; epsXi *= eps;
+ typename TargetSpace::EmbeddedTangentVector epsEta = B[j][j2]; epsEta *= eps;
+
+ typename TargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1;
+ typename TargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1;
+
+ if (i==j)
+ forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
+ else {
+ forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi);
+ forwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],epsEta);
+ }
+
+ if (i==j)
+ backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi+minusEpsEta);
+ else {
+ backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
+ backwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],minusEpsEta);
+ }
+
+ RT forwardValue = energy(element, localFiniteElement, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
+ RT backwardValue = energy(element, localFiniteElement, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
+
+ this->A_[i][j][i2][j2] = this->A_[j][i][j2][i2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
+
+ }
+ }
+ }
+ }
+}
+
+#endif
+