diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 6ff986da031af0e10360a25e5a5bf1c38825d82b..1552f595fefad0f6710ed2a108530fe023af1011 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/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
diff --git a/dune/gfe/localquickanddirtyfefunction.hh b/dune/gfe/localquickanddirtyfefunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b77022b098b68df2174bb58695075e23a61c4bfd
--- /dev/null
+++ b/dune/gfe/localquickanddirtyfefunction.hh
@@ -0,0 +1,172 @@
+#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