New GFE assembler using finite differences for gradient and Hesse matrix

This used to be in the base class as default implementations, but I think
it is nicer to have it in a derived class.  In particular, this will
eventually allow to do the FD approximation with high-precision number types.

[[Imported from SVN: r9842]]

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 \

dune/gfe/localgeodesicfefdstiffness.hh

+#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
+