Introduce a dedicated class LocalQuickAndDirtyFEFunction to produce initial...

Introduce a dedicated class LocalQuickAndDirtyFEFunction to produce initial iterates for the GFE minimization problems

Previously, a LocalProjectedFEFunction was used for this.  This turned out to be a bad idea for two reasons:
- Evaluating such a function can actually pretty expensive.  For example,
  to project onto SO(3), the polar decomposition needs to be evaluated by
  an iterative method.
- Screwing up the LocalProjectedFEFunction implementation can actually screw up
  the GFE implementation as well in very subtle ways.  Specifically, if the
  LocalProjectedFEFunction is not ADOL-C-differentiable, then the GFE function
  will not be differentiable, either.  This fact cost me a few hours searching...

dune/gfe/localgeodesicfefunction.hh

@@ -9,7 +9,7 @@
 
 #include <dune/gfe/averagedistanceassembler.hh>
 #include <dune/gfe/targetspacertrsolver.hh>
-#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/localquickanddirtyfefunction.hh>
 #include <dune/gfe/rigidbodymotion.hh>
 
 #include <dune/gfe/tensor3.hh>
@@ -186,7 +186,7 @@ evaluate(const Dune::FieldVector<ctype, dim>& local) const
     AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
 
     // Create a reasonable initial iterate for the iterative solver
-    Dune::GFE::LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> localProjectedFEFunction(localFiniteElement_, coefficients_);
+    Dune::GFE::LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> localProjectedFEFunction(localFiniteElement_, coefficients_);
     TargetSpace initialIterate = localProjectedFEFunction.evaluate(local);
 
     // Iteratively solve the GFE minimization problem

dune/gfe/localquickanddirtyfefunction.hh

+#ifndef DUNE_GFE_LOCALQUICKANDDIRTYFEFUNCTION_HH
+#define DUNE_GFE_LOCALQUICKANDDIRTYFEFUNCTION_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/type.hh>
+
+#include <dune/gfe/rotation.hh>
+#include <dune/gfe/linearalgebra.hh>
+
+namespace Dune {
+
+  namespace GFE {
+
+    /** \brief Interpolate on a manifold, as fast as we can
+     *
+     * This class implements interpolation of values on a manifold in a 'quick-and-dirty' way.
+     * No particular interpolation rule is used consistenly for all manifolds.  Rather, it is
+     * used whatever is quickest and 'reasonable' for any given space.  The reason to have this
+     * is to provide initial iterates for the iterative solvers used to solve the local
+     * geodesic-FE minimization problems.
+     *
+     * \note In particular, we currently do not implement the derivative of the interpolation,
+     *   because it is not needed for producing initial iterates!
+     *
+     * \tparam dim Dimension of the reference element
+     * \tparam ctype Type used for coordinates on the reference element
+     * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
+     * \tparam TargetSpace The manifold that the function takes its values in
+     */
+    template <int dim, class ctype, class LocalFiniteElement, class TS>
+    class LocalQuickAndDirtyFEFunction
+    {
+    public:
+      using TargetSpace=TS;
+    private:
+      typedef typename TargetSpace::ctype RT;
+
+      typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+      static const int embeddedDim = EmbeddedTangentVector::dimension;
+
+      static const int spaceDim = TargetSpace::TangentVector::dimension;
+
+    public:
+
+      /** \brief The type used for derivatives */
+      typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
+
+      /** \brief Constructor
+       * \param localFiniteElement A Lagrangian finite element that provides the interpolation points
+       * \param coefficients Values of the function at the Lagrange points
+       */
+      LocalQuickAndDirtyFEFunction(const LocalFiniteElement& localFiniteElement,
+                               const std::vector<TargetSpace>& coefficients)
+      : localFiniteElement_(localFiniteElement),
+      coefficients_(coefficients)
+      {
+        assert(localFiniteElement_.localBasis().size() == coefficients_.size());
+      }
+
+      /** \brief The number of Lagrange points */
+      unsigned int size() const
+      {
+        return localFiniteElement_.localBasis().size();
+      }
+
+      /** \brief The type of the reference element */
+      Dune::GeometryType type() const
+      {
+        return localFiniteElement_.type();
+      }
+
+      /** \brief Evaluate the function */
+      TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const;
+#if 0
+      /** \brief Evaluate the derivative of the function */
+      DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const;
+
+      /** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
+       *        \param local Local coordinates in the reference element where to evaluate the derivative
+       *        \param q Value of the local gfe function at 'local'.  If you provide something wrong here the result will be wrong, too!
+       */
+      DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
+                                        const TargetSpace& q) const;
+#endif
+      /** \brief Get the i'th base coefficient. */
+      TargetSpace coefficient(int i) const
+      {
+        return coefficients_[i];
+      }
+    private:
+
+      /** \brief The scalar local finite element, which provides the weighting factors
+       *        \todo We really only need the local basis
+       */
+      const LocalFiniteElement& localFiniteElement_;
+
+      /** \brief The coefficient vector */
+      std::vector<TargetSpace> coefficients_;
+
+    };
+
+    template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+    TargetSpace LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
+    evaluate(const Dune::FieldVector<ctype, dim>& local) const
+    {
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<Dune::FieldVector<ctype,1> > w;
+      localFiniteElement_.localBasis().evaluateFunction(local,w);
+
+      typename TargetSpace::CoordinateType c(0);
+      for (size_t i=0; i<coefficients_.size(); i++)
+        c.axpy(w[i][0], coefficients_[i].globalCoordinates());
+
+      return TargetSpace(c);
+    }
+#if 0
+    template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+    typename LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
+    LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
+    evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
+    {
+      // the function value at the point where we are evaluating the derivative
+      TargetSpace q = evaluate(local);
+
+      // Actually compute the derivative
+      return evaluateDerivative(local,q);
+    }
+
+    template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+    typename LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
+    LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
+    evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace& q) const
+    {
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<Dune::FieldVector<ctype,1> > w;
+      localFiniteElement_.localBasis().evaluateFunction(local,w);
+
+      std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
+      localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
+
+      typename TargetSpace::CoordinateType embeddedInterpolation(0);
+      for (size_t i=0; i<coefficients_.size(); i++)
+        embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());
+
+      Dune::FieldMatrix<RT,embeddedDim,dim> derivative(0);
+      for (size_t i=0; i<embeddedDim; i++)
+        for (size_t j=0; j<dim; j++)
+          for (size_t k=0; k<coefficients_.size(); k++)
+            derivative[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];
+
+      auto derivativeOfProjection = TargetSpace::derivativeOfProjection(embeddedInterpolation);
+
+      typename LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType result;
+
+      for (size_t i=0; i<result.N(); i++)
+        for (size_t j=0; j<result.M(); j++)
+        {
+          result[i][j] = 0;
+          for (size_t k=0; k<derivativeOfProjection.M(); k++)
+            result[i][j] += derivativeOfProjection[i][k]*derivative[k][j];
+        }
+
+      return result;
+    }
+#endif
+  }   // namespace GFE
+
+}   // namespace Dune
+#endif