dune/gfe/densities/localdensity.hh

 #ifndef DUNE_GFE_DENSITIES_LOCALDENSITY_HH
 #define DUNE_GFE_DENSITIES_LOCALDENSITY_HH
 
+#include <optional>
+
 #include <dune/common/fmatrix.hh>
 
 #include <adolc/adalloc.h>
@@ -8,19 +10,25 @@
 
 #include <dune/istl/matrix.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   /** \brief A base class for energy densities to be evaluated in an integral energy
    *
-   * \tparam Position The evaluation point in the integration domain
+   * \tparam ElementOrIntersection The domain of the density.
+   *   Can be either a grid element (i.e., a codimension-0 Entity)
+   *   or an Intersection.
    * \tparam TargetSpace Type for the function value
    */
-  template<class Position, class TargetSpace>
+  template<class ElementOrIntersection, class TargetSpace>
   class LocalDensity
   {
+    using Geometry = typename ElementOrIntersection::Geometry;
+    using LocalCoordinate = typename Geometry::LocalCoordinate;
+
     using field_type = typename TargetSpace::field_type;
     using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
-    using DerivativeType = FieldMatrix<field_type,TargetSpace::EmbeddedTangentVector::dimension,Position::size()>;
+    using DerivativeType = FieldMatrix<field_type,TargetSpace::EmbeddedTangentVector::dimension,LocalCoordinate::size ()>;
 
     // Number of independent variables
     static constexpr auto m = TargetSpace::EmbeddedTangentVector::dimension + DerivativeType::rows*DerivativeType::cols;
@@ -35,6 +43,10 @@ namespace Dune::GFE {
     double*** Yppp_;   /* results of hov_wk_forward  */
     double*** Zppp_;   /* result of hos_ov_reverse */
 
+  protected:
+    /** \brief The element or intersection that this density is defined on */
+    std::optional<ElementOrIntersection> elementOrIntersection_ = std::nullopt;
+
   public:
 
     LocalDensity()
@@ -51,20 +63,29 @@ namespace Dune::GFE {
 
     }
 
-    ~LocalDensity()
+    virtual ~LocalDensity()
     {
       myfree3(densityTangent_);
       myfree3(Yppp_);
       myfree3(Zppp_);
     }
 
+    /** \brief Bind the density to a `ElementOrIntersection` object
+     *
+     * Stores a copy of the `ElementOrIntersection`'s geometry.
+     **/
+    virtual void bind(const ElementOrIntersection& elementOrIntersection)
+    {
+      elementOrIntersection_.emplace(elementOrIntersection);
+    }
+
     /** \brief Evaluate the density for a given value and first derivative
      *
      * \param x The current position
      * \param value The value of the integrand at x
      * \param derivative The derivative of the integrand at x
      */
-    virtual field_type operator() (const Position& x,
+    virtual field_type operator() (const LocalCoordinate& x,
                                    const typename TargetSpace::CoordinateType& value,
                                    const DerivativeType& derivative) const = 0;
 
@@ -72,7 +93,7 @@ namespace Dune::GFE {
      *
      * The default implementation here uses ADOL-C for this.
      */
-    virtual void derivatives(const Position& x,
+    virtual void derivatives(const LocalCoordinate& x,
                              const TargetSpace& value,
                              const DerivativeType& derivative,
                              field_type& densityValue,
@@ -208,7 +229,7 @@ namespace Dune::GFE {
     }
 
     // Construct a copy of this density but using 'adouble' as the number type
-    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const = 0;
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const = 0;
 
     /** \brief Whether the density depends on the 'value' parameter
      *
@@ -216,7 +237,8 @@ namespace Dune::GFE {
      * for one factor space only.
      *
      * \param factor The factor space that is being asked about.
-     *   The default value -1 means: Does any of the factors depend on the value?
+     *   The default value -1 means: Does the density depend on the value
+     *   of any of the factors?
      */
     virtual bool dependsOnValue(int factor=-1) const = 0;
 
@@ -226,7 +248,8 @@ namespace Dune::GFE {
      * for one factor space only.
      *
      * \param factor The factor space that is being asked about
-     *   The default value -1 means: Does any of the factors depend on the derivative?
+     *   The default value -1 means: Does the density depend on the derivative
+     *   of any of the factors?
      */
     virtual bool dependsOnDerivative(int factor=-1) const = 0;
   };

dune/gfe/densities/planarcosseratshelldensity.hh

@@ -21,14 +21,15 @@
 
 namespace Dune::GFE
 {
-  template<class Position, class field_type>
+  template<class ElementOrIntersection, class field_type>
   class PlanarCosseratShellDensity final
-    : public GFE::LocalDensity<Position, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
+    : public GFE::LocalDensity<ElementOrIntersection, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
   {
   private:
+    using LocalCoordinate = typename ElementOrIntersection::Geometry::LocalCoordinate;
     // PlanarCosseratShellDensity only works for 2d grids
-    static_assert(Position::size()==2);
-    static constexpr int gridDim = Position::size();
+    static_assert(LocalCoordinate::size()==2);
+    static constexpr int gridDim = LocalCoordinate::size();
 
     // The target space with 'adouble' as the number type
     using ATargetSpace = GFE::ProductManifold<RealTuple<adouble,3>,Rotation<adouble,3> >;
@@ -210,7 +211,7 @@ namespace Dune::GFE
 
     /** \brief Evaluate the density
      */
-    virtual field_type operator() (const Position& x,
+    virtual field_type operator() (const LocalCoordinate& x,
                                    const typename GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >::CoordinateType& value,
                                    const FieldMatrix<field_type,7,gridDim>& derivative) const override
     {
@@ -253,9 +254,13 @@ namespace Dune::GFE
     }
 
     // Construct a copy of this density but using 'adouble' as the number type
-    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const override
+    virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const override
     {
-      return std::make_unique<PlanarCosseratShellDensity<Position,adouble> >(thickness_, mu_, lambda_, mu_c_, L_c_, q_, kappa_);
+      auto result = std::make_unique<PlanarCosseratShellDensity<ElementOrIntersection,adouble> >(thickness_, mu_, lambda_, mu_c_, L_c_, q_, kappa_);
+
+      if (this->elementOrIntersection_)
+        result->bind(*this->elementOrIntersection_);
+      return result;
     }
 
     /** \brief The density depends on the microrotation value, but not on the deformation
@@ -291,6 +296,6 @@ namespace Dune::GFE
     double kappa_;
   };
 
-}
+}  // namespace Dune::GFE
 
 #endif   //#ifndef DUNE_GFE_DENSITIES_PLANARCOSSERATSHELLDENSITY_HH

dune/gfe/filereader.hh

@@ -8,36 +8,37 @@
 // This is a custom file format parser which reads in a grid deformation file
 // The file must contain lines of the format <grid vertex>:<displacement or rotation>
 // !!! THIS IS A TEMPORARY SOLUTION AND WILL BE REPLACED BY THE Dune::VtkReader in dune-vtk/dune/vtk/vtkreader.hh !!!
-namespace Dune {
-  namespace GFE {
-    // Convert the pairs {grid vertex, vector of dimension d} in the given file to a map
-    template <int d>
-    static std::unordered_map<std::string, FieldVector<double,d> > transformFileToMap(std::string pathToFile) {
-
-      std::unordered_map<std::string, FieldVector<double,d> > map;
-
-      std::string line, displacement, entry;
-
-      std::ifstream file(pathToFile, std::ios::in);
-
-      if (file.is_open()) {
-        while (std::getline(file, line)) {
-          size_t j = 0;
-          size_t pos = line.find(":");
-          displacement = line.substr(pos + 1);
-          line.erase(pos);
-          std::stringstream entries(displacement);
-          FieldVector<double,d> vector(0);
-          while(entries >> entry)
-            vector[j++] = std::stod(entry);
-          map.insert({line,vector});
-        }
-        file.close();
-      } else {
-        DUNE_THROW(Exception, "Error: Could not open the file " + pathToFile + " !");
+namespace Dune::GFE
+{
+
+  // Convert the pairs {grid vertex, vector of dimension d} in the given file to a map
+  template <int d>
+  static std::unordered_map<std::string, FieldVector<double,d> > transformFileToMap(std::string pathToFile) {
+
+    std::unordered_map<std::string, FieldVector<double,d> > map;
+
+    std::string line, displacement, entry;
+
+    std::ifstream file(pathToFile, std::ios::in);
+
+    if (file.is_open()) {
+      while (std::getline(file, line)) {
+        size_t j = 0;
+        size_t pos = line.find(":");
+        displacement = line.substr(pos + 1);
+        line.erase(pos);
+        std::stringstream entries(displacement);
+        FieldVector<double,d> vector(0);
+        while(entries >> entry)
+          vector[j++] = std::stod(entry);
+        map.insert({line,vector});
       }
-      return map;
+      file.close();
+    } else {
+      DUNE_THROW(Exception, "Error: Could not open the file " + pathToFile + " !");
     }
+    return map;
   }
+
 }
 #endif

dune/gfe/functions/CMakeLists.txt

+install(FILES
+        embeddedglobalgfefunction.hh
+        globalgfefunction.hh
+        interpolationderivatives.hh
+        localgeodesicfefunction.hh
+        localprojectedfefunction.hh
+        localquickanddirtyfefunction.hh
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gfe/functions)

dune/gfe/functions/discretekirchhoffbendingisometry.hh

+#ifndef DUNE_GFE_FUNCTIONS_DISCRETEKIRCHHOFFBENDINGISOMETRY_HH
+#define DUNE_GFE_FUNCTIONS_DISCRETEKIRCHHOFFBENDINGISOMETRY_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+#include <dune/gfe/functions/localisometrycomponentfunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+#include <dune/gfe/tensor3.hh>
+
+#include <dune/functions/functionspacebases/cubichermitebasis.hh>
+
+namespace Dune::GFE
+{
+  /** \brief Interpolate a Discrete Kirchhoff finite element function from a set of RealTuple x Rotation coefficients
+   *
+   * The Discrete Kirchhoff finite element is a Hermite-type element.
+   * The Rotation component of the coefficient set represents the deformation gradient.
+   *
+   * The class stores both a global coefficient vector 'coefficients_'
+   * as well as a local coefficient vector 'localHermiteCoefficients_'
+   * that is used for local evaluation after binding to an element.
+   *
+   *
+   * \tparam DiscreteKirchhoffBasis The basis used to compute function values
+   * \tparam CoefficientBasis The basis used to index the coefficient vector.
+   * \tparam Coefficients The container of global coefficients
+   */
+  template <class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients>
+  class DiscreteKirchhoffBendingIsometry
+  {
+    using GridView = typename DiscreteKirchhoffBasis::GridView;
+    using ctype = typename GridView::ctype;
+
+  public:
+    constexpr static int gridDim = GridView::dimension;
+    using Coefficient = typename Coefficients::value_type;
+    using RT = typename Coefficient::ctype;
+
+    static constexpr int components_ = 3;
+
+    typedef GFE::LocalProjectedFEFunction<CoefficientBasis, GFE::ProductManifold<RealTuple<RT,components_>, Rotation<RT,components_> > > LocalInterpolationRule;
+
+    /** \brief The type used for derivatives */
+    typedef Dune::FieldMatrix<RT, components_, gridDim> DerivativeType;
+
+    /** \brief The type used for discrete Hessian */
+    typedef Tensor3<RT, components_, gridDim, gridDim> discreteHessianType;
+
+
+    /** \brief Constructor
+     * \param basis An object of Hermite basis type
+     * \param coefficientBasis An object of type LagrangeBasis<GridView,1>
+     * \param globalIsometryCoefficients Values and derivatives of the function at the Lagrange points
+     */
+    DiscreteKirchhoffBendingIsometry(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
+                                     const CoefficientBasis& coefficientBasis,
+                                     Coefficients& globalIsometryCoefficients)
+      : coefficientBasis_(coefficientBasis),
+      basis_(discreteKirchhoffBasis),
+      localView_(basis_),
+      localViewCoefficient_(coefficientBasis_),
+      globalIsometryCoefficients_(globalIsometryCoefficients)
+    {}
+
+    void bind(const typename GridView::template Codim<0>::Entity& element)
+    {
+      localView_.bind(element);
+
+      /** extract the local coefficients from the global coefficient vector */
+      std::vector<Coefficient> localIsometryCoefficients;
+      getLocalIsometryCoefficients(localIsometryCoefficients,element);
+      updateLocalHermiteCoefficients(localIsometryCoefficients,element);
+
+      /**
+       * @brief  Get Coefficients of the discrete Gradient
+       *
+       * The discrete Gradient is a special linear combination represented in a [P2]^2 - (Lagrange) space
+       * The coefficients of this linear combination correspond to certain linear combinations of the Gradients of
+       * the deformation (hermite) basis. The coefficients are stored in the form [Basisfunctions x components x gridDim]
+       * in a BlockVector<FieldMatrix<RT, 3, gridDim>>.
+       */
+      Impl::discreteGradientCoefficients(discreteJacobianCoefficients_,rotationLFE_,*this,element);
+    }
+
+    /** bind to an element and update the coefficient member variables for a set of new local coefficients */
+    void bind(const typename GridView::template Codim<0>::Entity& element,const Coefficients& newLocalCoefficients)
+    {
+      this->bind(element);
+
+      // Update the global isometry coefficients.
+      for (unsigned int vertex=0; vertex<3; ++vertex)
+      {
+        size_t localIdx = localViewCoefficient_.tree().localIndex(vertex);
+        size_t globalIdx = localViewCoefficient_.index(localIdx);
+        globalIsometryCoefficients_[globalIdx] = newLocalCoefficients[vertex];
+      }
+
+      // Update the local hermite coefficients.
+      updateLocalHermiteCoefficients(newLocalCoefficients,element);
+
+      // After setting a new local configurarion, the coefficients of the discrete Gradient need to be updated as well.
+      Impl::discreteGradientCoefficients(discreteJacobianCoefficients_,rotationLFE_,*this,element);
+    }
+
+    /** \brief The type of the element we are bound to */
+    Dune::GeometryType type() const
+    {
+      return localView_.element().type();
+    }
+
+
+    /** \brief Evaluate the function */
+    auto evaluate(const Dune::FieldVector<ctype, gridDim>& local) const
+    {
+      const auto &localFiniteElement = localView_.tree().child(0).finiteElement();
+
+      // Evaluate the shape functions
+      std::vector<FieldVector<double, 1> > values;
+      localFiniteElement.localBasis().evaluateFunction(local, values);
+
+      FieldVector<RT,components_> result(0);
+
+      for(size_t i=0; i<localFiniteElement.size(); i++)
+        result.axpy(values[i][0], localHermiteCoefficients_[i]);
+
+      return (RealTuple<RT, components_>)result;
+    }
+
+    /** \brief Evaluate the derivative of the function */
+    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, gridDim>& local) const
+    {
+      const auto &localFiniteElement = localView_.tree().child(0).finiteElement();
+
+      const auto jacobianInverse = localView_.element().geometry().jacobianInverse(local);
+      std::vector<FieldMatrix<double,1,gridDim> > jacobianValues;
+      localFiniteElement.localBasis().evaluateJacobian(local, jacobianValues);
+      for (size_t i = 0; i<jacobianValues.size(); i++)
+        jacobianValues[i] = jacobianValues[i] * jacobianInverse;
+
+      DerivativeType result(0);
+
+      for(size_t i=0; i<localFiniteElement.size(); i++)
+        for (unsigned int k = 0; k<components_; k++)
+          for (unsigned int j = 0; j<gridDim; ++j)
+            result[k][j] += (jacobianValues[i][0][j] * localHermiteCoefficients_[i][k]);
+
+      return result;
+    }
+
+    /** \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, gridDim>& local,
+                                      const Coefficient& q) const
+    {
+      return evaluateDerivative(local);
+    }
+
+
+
+    /** \brief Evaluate the discrete gradient operator (This quantity lives in the P2-rotation space) */
+    DerivativeType evaluateDiscreteGradient(const Dune::FieldVector<ctype, gridDim>& local) const
+    {
+      //Get values of the P2-Basis functions on local point
+      std::vector<FieldVector<double,1> > basisValues;
+      rotationLFE_.localBasis().evaluateFunction(local, basisValues);
+
+      DerivativeType discreteGradient(0);
+
+      for (int k=0; k<components_; k++)
+        for (int l=0; l<gridDim; l++)
+          for (std::size_t i=0; i<rotationLFE_.size(); i++)
+            discreteGradient[k][l]  += discreteJacobianCoefficients_[i][k][l]*basisValues[i];
+
+      return discreteGradient;
+    }
+
+
+    /** \brief Evaluate the discrete hessian operator (This quantity lives in the P2-rotation space) */
+    discreteHessianType evaluateDiscreteHessian(const Dune::FieldVector<ctype, gridDim>& local) const
+    {
+      // Get gradient values of the P2-Basis functions on local point
+      const auto geometryJacobianInverse = localView_.element().geometry().jacobianInverse(local);
+      std::vector<FieldMatrix<double,1,gridDim> > gradients;
+      rotationLFE_.localBasis().evaluateJacobian(local, gradients);
+
+      for (size_t i=0; i< gradients.size(); i++)
+        gradients[i] = gradients[i] * geometryJacobianInverse;
+
+      discreteHessianType discreteHessian(0);
+
+      for (int k=0; k<components_; k++)
+        for (int l=0; l<gridDim; l++)
+          for (std::size_t i=0; i<rotationLFE_.size(); i++)
+          {
+            discreteHessian[k][l][0]  += discreteJacobianCoefficients_[i][k][l]*gradients[i][0][0];
+            discreteHessian[k][l][1]  += discreteJacobianCoefficients_[i][k][l]*gradients[i][0][1];
+          }
+
+      return discreteHessian;
+    }
+
+
+    //! Obtain the grid view that the basis is defined on
+    const GridView &gridView() const
+    {
+      return basis_.gridView();
+    }
+
+    void updateglobalIsometryCoefficients(const Coefficients& newCoefficients)
+    {
+      globalIsometryCoefficients_ = newCoefficients;
+    }
+
+    auto getDiscreteJacobianCoefficients()
+    {
+      return discreteJacobianCoefficients_;
+    }
+
+    /**
+        Update the local hermite basis coefficient vector
+     */
+    void updateLocalHermiteCoefficients(const Coefficients& localIsometryCoefficients,const typename GridView::template Codim<0>::Entity& element)
+    {
+      localViewCoefficient_.bind(element);
+
+      //Create a LocalProjected-Finite element from the local coefficients used for interpolation.
+      LocalInterpolationRule localPBfunction{coefficientBasis_};
+      localPBfunction.bind(element,localIsometryCoefficients);
+
+      /**
+          Interpolate into the local hermite space for each component.
+       */
+      const auto &localHermiteFiniteElement = localView_.tree().child(0).finiteElement();
+
+      for(size_t k=0; k<components_; k++)
+      {
+        std::vector<RT> hermiteComponentCoefficients;
+        Dune::GFE::Impl::LocalIsometryComponentFunction<double,LocalInterpolationRule> localIsometryComponentFunction(localPBfunction,k);
+
+        localHermiteFiniteElement.localInterpolation().interpolate(localIsometryComponentFunction,hermiteComponentCoefficients);
+
+        for(size_t i=0; i<localHermiteFiniteElement.size(); i++)
+          localHermiteCoefficients_[i][k] = hermiteComponentCoefficients[i];
+      }
+    }
+
+    /** Extract the local isometry coefficients. */
+    void getLocalIsometryCoefficients(std::vector<Coefficient>& in,const typename GridView::template Codim<0>::Entity& element)
+    {
+      localViewCoefficient_.bind(element);
+
+      in.resize(3);
+      for (unsigned int vertex = 0; vertex<3; ++vertex)
+      {
+        size_t localIdx = localViewCoefficient_.tree().localIndex(vertex);
+        size_t globalIdx = localViewCoefficient_.index(localIdx);
+        in[vertex] = globalIsometryCoefficients_[globalIdx];
+      }
+    }
+
+    /** \brief Return the representation LocalFiniteElement for the rotation space */
+    auto const &getRotationFiniteElement() const
+    {
+      return rotationLFE_;
+    }
+
+    /** \brief Return the discrete Jacobian coefficients */
+    BlockVector<FieldMatrix<RT, components_, gridDim> > getDiscreteJacobianCoefficients() const
+    {
+      return discreteJacobianCoefficients_;
+    }
+
+    /** \brief Get the i'th base coefficient. */
+    Coefficients coefficient(int i) const
+    {
+      return globalIsometryCoefficients_[i];
+    }
+
+  private:
+
+    /** \brief The Lagrange basis used to assign manifold-valued degrees of freedom */
+    const CoefficientBasis& coefficientBasis_;
+
+    /** \brief The hermite basis used to represent deformations */
+    const DiscreteKirchhoffBasis& basis_;
+
+    /** \brief LocalFiniteElement (P2-Lagrange) used to represent the discrete Jacobian (rotation) on the current element. */
+    Dune::LagrangeSimplexLocalFiniteElement<ctype, double, gridDim, 2> rotationLFE_;
+
+    /** \brief The coefficients of the discrete Gradient/Hessian operator */
+    BlockVector<FieldMatrix<RT, components_, gridDim> > discreteJacobianCoefficients_;
+
+    mutable typename DiscreteKirchhoffBasis::LocalView localView_;
+    mutable typename CoefficientBasis::LocalView localViewCoefficient_;
+
+    /** \brief The global coefficient vector */
+    Coefficients& globalIsometryCoefficients_;
+
+    static constexpr int localHermiteBasisSize_ = Dune::Functions::Impl::CubicHermiteLocalCoefficients<gridDim,true>::size();
+
+    /** \brief The local coefficient vector of the hermite basis. */
+    std::array<Dune::FieldVector<RT,components_>,localHermiteBasisSize_> localHermiteCoefficients_;
+
+  };
+
+}  // end namespace Dune::GFE
+#endif   // DUNE_GFE_FUNCTIONS_DISCRETEKIRCHHOFFBENDINGISOMETRY_HH

dune/gfe/embeddedglobalgfefunction.hh → dune/gfe/functions/embeddedglobalgfefunction.hh

 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 // vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_GFE_EMBEDDEDGLOBALGFEFUNCTION_HH
-#define DUNE_GFE_EMBEDDEDGLOBALGFEFUNCTION_HH
+#ifndef DUNE_GFE_FUNCTIONS_EMBEDDEDGLOBALGFEFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_EMBEDDEDGLOBALGFEFUNCTION_HH
 
 #include <memory>
 #include <optional>
@@ -16,10 +16,10 @@
 #include <dune/functions/gridfunctions/gridfunction.hh>
 #include <dune/functions/backends/concepts.hh>
 
-#include <dune/gfe/globalgfefunction.hh>
-
-namespace Dune::GFE {
+#include <dune/gfe/functions/globalgfefunction.hh>
 
+namespace Dune::GFE
+{
   template<typename EGGF>
   class EmbeddedGlobalGFEFunctionDerivative;
 
@@ -64,7 +64,7 @@ namespace Dune::GFE {
     using Domain = typename Base::Domain;
     using Range = typename TargetSpace::CoordinateType;
 
-    using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+    using Traits = Functions::Imp::GridFunctionTraits<Range (Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
 
     class LocalFunction
       : public Base::LocalFunctionBase
@@ -97,7 +97,7 @@ namespace Dune::GFE {
        */
       Range operator()(const Domain& x) const
       {
-        return this->localInterpolationRule_->evaluate(x).globalCoordinates();
+        return this->localInterpolationRule_.evaluate(x).globalCoordinates();
       }
 
       //! Local function of the derivative
@@ -177,7 +177,8 @@ namespace Dune::GFE {
         localCoeff[i] = this->dofs()[localView.index(i)];
 
       // create local gfe function
-      LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
+      LocalInterpolationRule localInterpolationRule(localView.globalBasis());
+      localInterpolationRule.bind(element,localCoeff);
       return localInterpolationRule.evaluate(local).globalCoordinates();
     }
 
@@ -201,7 +202,8 @@ namespace Dune::GFE {
         localCoeff[i] = this->dofs()[localView.index(i)];
 
       // create local gfe function
-      LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
+      LocalInterpolationRule localInterpolationRule(localView.globalBasis());
+      localInterpolationRule.bind(element,localCoeff);
 
       // use it to evaluate the derivative
       auto refJac = localInterpolationRule.evaluateDerivative(local);
@@ -256,7 +258,7 @@ namespace Dune::GFE {
     using Domain = typename Base::Domain;
     using Range = typename Functions::SignatureTraits<typename EmbeddedGlobalGFEFunction::Traits::DerivativeInterface>::Range;
 
-    using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+    using Traits = Functions::Imp::GridFunctionTraits<Range (Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
 
     /**
      * \brief local function evaluating the derivative in reference coordinates
@@ -319,7 +321,7 @@ namespace Dune::GFE {
       Range operator()(const Domain& x) const
       {
         // Jacobian with respect to local coordinates
-        auto refJac = this->localInterpolationRule_->evaluateDerivative(x);
+        auto refJac = this->localInterpolationRule_.evaluateDerivative(x);
 
         // Transform to world coordinates
         return refJac * geometry_->jacobianInverse(x);
@@ -386,4 +388,4 @@ namespace Dune::GFE {
 
 } // namespace Dune::GFE
 
-#endif // DUNE_GFE_EMBEDDEDGLOBALGFEFUNCTION_HH
+#endif // DUNE_GFE_FUNCTIONS_EMBEDDEDGLOBALGFEFUNCTION_HH

dune/gfe/globalgfefunction.hh → dune/gfe/functions/globalgfefunction.hh

 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 // vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_GFE_GLOBALGFEFUNCTION_HH
-#define DUNE_GFE_GLOBALGFEFUNCTION_HH
+#ifndef DUNE_GFE_FUNCTIONS_GLOBALGFEFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_GLOBALGFEFUNCTION_HH
 
 #include <memory>
 #include <optional>
@@ -20,8 +20,8 @@
 #include <dune/fufem/functions/virtualgridfunction.hh>
 #endif
 
-namespace Dune::GFE {
-
+namespace Dune::GFE
+{
   namespace Impl {
 
 #if DUNE_VERSION_LTE(DUNE_FUFEM, 2, 9)
@@ -102,6 +102,7 @@ namespace Dune::GFE {
 #endif
           : data_(data)
           , localView_(data_->basis->localView())
+          , localInterpolationRule_(*data_->basis)
         {
           localDoFs_.reserve(localView_.maxSize());
         }
@@ -115,6 +116,7 @@ namespace Dune::GFE {
         LocalFunctionBase(const LocalFunctionBase& other)
           : data_(other.data_)
           , localView_(other.localView_)
+          , localInterpolationRule_(other.localInterpolationRule_)
         {
           localDoFs_.reserve(localView_.maxSize());
           if (bound())
@@ -160,8 +162,7 @@ namespace Dune::GFE {
           }
 
           // create local GFE function
-          // TODO Store this object by value
-          localInterpolationRule_ = std::make_unique<LocalInterpolationRule>(this->localView_.tree().finiteElement(),localDoFs_);
+          localInterpolationRule_.bind(element,localDoFs_);
         }
 
         //! Unbind the local-function.
@@ -191,7 +192,7 @@ namespace Dune::GFE {
 #endif
         LocalView localView_;
         std::vector<Coefficient> localDoFs_;
-        std::unique_ptr<LocalInterpolationRule> localInterpolationRule_;
+        LocalInterpolationRule localInterpolationRule_;
       };
 
     protected:
@@ -275,7 +276,7 @@ namespace Dune::GFE {
     using Domain = typename Base::Domain;
     using Range = typename TargetSpace::CoordinateType;
 
-    using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+    using Traits = Functions::Imp::GridFunctionTraits<Range (Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
 
     class LocalFunction
       : public Base::LocalFunctionBase
@@ -308,7 +309,7 @@ namespace Dune::GFE {
        */
       Range operator()(const Domain& x) const
       {
-        return this->localInterpolationRule_->evaluate(x).globalCoordinates();
+        return this->localInterpolationRule_.evaluate(x).globalCoordinates();
       }
 
       //! Local function of the derivative
@@ -414,7 +415,7 @@ namespace Dune::GFE {
     using Domain = typename Base::Domain;
     using Range = typename Functions::SignatureTraits<typename GlobalGFEFunction::Traits::DerivativeInterface>::Range;
 
-    using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+    using Traits = Functions::Imp::GridFunctionTraits<Range (Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
 
     /**
      * \brief local function evaluating the derivative in reference coordinates
@@ -477,7 +478,7 @@ namespace Dune::GFE {
       Range operator()(const Domain& x) const
       {
         // Jacobian with respect to local coordinates
-        auto refJac = this->localInterpolationRule_->evaluateDerivative(x);
+        auto refJac = this->localInterpolationRule_.evaluateDerivative(x);
 
         // Transform to world coordinates
         return refJac * geometry_->jacobianInverse(x);
@@ -544,4 +545,4 @@ namespace Dune::GFE {
 
 } // namespace Dune::GFE
 
-#endif // DUNE_GFE_GLOBALGFEFUNCTION_HH
+#endif // DUNE_GFE_FUNCTIONS_GLOBALGFEFUNCTION_HH

dune/gfe/interpolationderivatives.hh → dune/gfe/functions/interpolationderivatives.hh

-#ifndef DUNE_GFE_INTERPOLATIONDERIVATIVES_HH
-#define DUNE_GFE_INTERPOLATIONDERIVATIVES_HH
+#ifndef DUNE_GFE_FUNCTIONS_INTERPOLATIONDERIVATIVES_HH
+#define DUNE_GFE_FUNCTIONS_INTERPOLATIONDERIVATIVES_HH
 
 // Includes for the ADOL-C automatic differentiation library
 #include <adolc/adolc.h>
@@ -8,8 +8,8 @@
 
 #include <dune/fufem/utilities/adolcnamespaceinjections.hh>
 
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/functions/localgeodesicfefunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
@@ -223,8 +223,8 @@ namespace Dune::GFE
       }
 
       // Create the functions, we want to tape the function evaluation and the evaluation of the derivatives
-      const auto& scalarFiniteElement = localInterpolationRule_.localFiniteElement();
-      ALocalInterpolationRule localGFEFunction(scalarFiniteElement,localAConfiguration);
+      ALocalInterpolationRule localGFEFunction(localInterpolationRule_.globalBasis());
+      localGFEFunction.bind(element,localAConfiguration);
 
       if (doValue_)
       {
@@ -427,10 +427,20 @@ namespace Dune::GFE
    * standard FE interpolation, and the derivatives with respect to the coefficients can be
    * computed much simpler and faster than for the general case.
    */
-  template <int gridDim, typename field_type, typename LocalFiniteElement,int dim>
-  class InterpolationDerivatives<LocalGeodesicFEFunction<gridDim, field_type, LocalFiniteElement, RealTuple<field_type,dim> > >
+  template <typename Basis, typename field_type,int dim>
+  class InterpolationDerivatives<LocalGeodesicFEFunction<Basis, RealTuple<field_type,dim> > >
   {
-    using LocalInterpolationRule = LocalGeodesicFEFunction<gridDim, field_type, LocalFiniteElement, RealTuple<field_type,dim> >;
+    using LocalInterpolationRule = LocalGeodesicFEFunction<Basis, RealTuple<field_type,dim> >;
+
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using WeightType = typename LocalBasis::Traits::RangeType;
+    using WeightJacobianType = typename LocalBasis::Traits::JacobianType;
+
     using TargetSpace = typename LocalInterpolationRule::TargetSpace;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
@@ -449,11 +459,10 @@ namespace Dune::GFE
     const bool doDerivative_;
 
     // Values of all scalar shape functions at the point we are bound to
-    std::vector<FieldVector<double,1> > shapeFunctionValues_;
+    std::vector<WeightType> shapeFunctionValues_;
 
     // Gradients of all scalar shape functions at the point we are bound to
-    // TODO: The second dimension must be WorldDim
-    std::vector<FieldMatrix<double,1,gridDim> > shapeFunctionGradients_;
+    std::vector<WeightJacobianType> shapeFunctionGradients_;
 
 
   public:
@@ -479,10 +488,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -546,7 +554,7 @@ namespace Dune::GFE
       // First derivatives of the function gradient wrt to the FE coefficients
       for (size_t i=0; i<nDofs; ++i)
         for (int j=0; j<blocksize; ++j)
-          for (int k=0; k<gridDim; ++k)
+          for (std::size_t k=0; k<gridDim; ++k)
             firstDerivative[blocksize + j*gridDim + k][i*blocksize+j] = shapeFunctionGradients_[i][0][k];
 
       // For RealTuple, firstDerivative and embeddedFirstDerivative coincide
@@ -568,11 +576,16 @@ namespace Dune::GFE
    * standard FE interpolation, and the derivatives with respect to the coefficients can be
    * computed much simpler and faster than for the general case.
    */
-  template <int gridDim, typename field_type, typename LocalFiniteElement,int dim>
-  class InterpolationDerivatives<LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, RealTuple<field_type,dim> > >
+  template <typename Basis, typename field_type, int dim>
+  class InterpolationDerivatives<LocalProjectedFEFunction<Basis, RealTuple<field_type,dim> > >
   {
     // TODO: The implementation here would be identical to the geodesic FE case
-    using LocalInterpolationRule = LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, RealTuple<field_type,dim> >;
+    using LocalInterpolationRule = LocalProjectedFEFunction<Basis, RealTuple<field_type,dim> >;
+
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
     using TargetSpace = typename LocalInterpolationRule::TargetSpace;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
@@ -618,10 +631,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -684,7 +696,7 @@ namespace Dune::GFE
       // First derivatives of the function gradient wrt to the FE coefficients
       for (size_t i=0; i<nDofs; ++i)
         for (int j=0; j<blocksize; ++j)
-          for (int k=0; k<gridDim; ++k)
+          for (std::size_t k=0; k<gridDim; ++k)
             firstDerivative[blocksize + j*gridDim + k][i*blocksize+j] = shapeFunctionGradients_[i][0][k];
 
       // For RealTuple, firstDerivative and embeddedFirstDerivative coincide
@@ -704,10 +716,19 @@ namespace Dune::GFE
    * interpolation.  No matter what the target space is, the interpolation is always
    * Euclidean in the surrounding space.
    */
-  template <int gridDim, typename ctype, typename LocalFiniteElement, typename TargetSpace>
-  class InterpolationDerivatives<LocalProjectedFEFunction<gridDim, ctype, LocalFiniteElement, TargetSpace, false> >
+  template <typename Basis, typename TargetSpace>
+  class InterpolationDerivatives<LocalProjectedFEFunction<Basis, TargetSpace, false> >
   {
-    using LocalInterpolationRule = LocalProjectedFEFunction<gridDim, ctype, LocalFiniteElement, TargetSpace, false>;
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using FERangeType = typename LocalBasis::Traits::RangeType;
+    using FEJacobianType = typename LocalBasis::Traits::JacobianType;
+
+    using LocalInterpolationRule = LocalProjectedFEFunction<Basis, TargetSpace, false>;
     using CoordinateType = typename TargetSpace::CoordinateType;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
@@ -727,11 +748,10 @@ namespace Dune::GFE
     const bool doDerivative_;
 
     // Values of all scalar shape functions at the point we are bound to
-    std::vector<FieldVector<double,1> > shapeFunctionValues_;
+    std::vector<FERangeType> shapeFunctionValues_;
 
     // Gradients of all scalar shape functions at the point we are bound to
-    // TODO: The second dimension must be WorldDim
-    std::vector<FieldMatrix<double,1,gridDim> > shapeFunctionGradients_;
+    std::vector<FEJacobianType> shapeFunctionGradients_;
 
     // TODO: Don't hardcode FieldMatrix
     std::vector<FieldMatrix<double,blocksize,embeddedBlocksize> > orthonormalFrames_;
@@ -763,10 +783,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -833,7 +852,7 @@ namespace Dune::GFE
       // First derivatives of the function gradient wrt to the FE coefficients
       for (size_t i=0; i<nDofs; ++i)
         for (int j=0; j<embeddedBlocksize; ++j)
-          for (int k=0; k<gridDim; ++k)
+          for (std::size_t k=0; k<gridDim; ++k)
             partialDerivative[embeddedBlocksize + j*gridDim + k][i*embeddedBlocksize+j] = shapeFunctionGradients_[i][0][k];
 
       // Euclidean derivative: Derivatives in the direction of the projected canonical vectors
@@ -1033,16 +1052,25 @@ namespace Dune::GFE
    * by hand with reasonable effort.  The corresponding implementation is faster than the one
    * based on ADOL-C.
    */
-  template <int gridDim, typename field_type, typename LocalFiniteElement,int dim>
-  class InterpolationDerivatives<LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, UnitVector<field_type,dim> > >
+  template <typename Basis, typename field_type,int dim>
+  class InterpolationDerivatives<LocalProjectedFEFunction<Basis, UnitVector<field_type,dim> > >
   {
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using FERangeType = typename LocalBasis::Traits::RangeType;
+    using FEJacobianType = typename LocalBasis::Traits::JacobianType;
+
     using TargetSpace = UnitVector<field_type,dim>;
     using CoordinateType = typename TargetSpace::CoordinateType;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
     constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
 
-    using LocalInterpolationRule = LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, TargetSpace>;
+    using LocalInterpolationRule = LocalProjectedFEFunction<Basis, TargetSpace>;
 
     /** \brief A vector that can be contracted with a derivative of the projection onto the unit sphere
      *
@@ -1226,11 +1254,10 @@ namespace Dune::GFE
     const bool doDerivative_;
 
     // Values of all scalar shape functions at the point we are bound to
-    std::vector<FieldVector<double,1> > shapeFunctionValues_;
+    std::vector<FERangeType> shapeFunctionValues_;
 
     // Gradients of all scalar shape functions at the point we are bound to
-    // TODO: The second dimension must be WorldDim
-    std::vector<FieldMatrix<double,1,gridDim> > shapeFunctionGradients_;
+    std::vector<FEJacobianType> shapeFunctionGradients_;
 
     // TODO: Is this needed at all?
     typename LocalInterpolationRule::DerivativeType euclideanDerivative_;
@@ -1283,10 +1310,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -1478,7 +1504,7 @@ namespace Dune::GFE
 
               if (doDerivative_)
               {
-                for (int beta=0; beta<gridDim; beta++)
+                for (std::size_t beta=0; beta<gridDim; beta++)
                 {
                   VectorToContractWith derivativeWeightsBeta(derivativeWeightsTransposed[beta]);
                   derivativeWeightsBeta.bind(projection_.x_);

dune/gfe/functions/localgeodesicfefunction.hh

+#ifndef DUNE_GFE_FUNCTIONS_LOCALGEODESICFEFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_LOCALGEODESICFEFUNCTION_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/type.hh>
+
+#include <dune/gfe/averagedistanceassembler.hh>
+#include <dune/gfe/targetspacertrsolver.hh>
+#include <dune/gfe/functions/localquickanddirtyfefunction.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+#include <dune/gfe/symmetricmatrix.hh>
+#include <dune/gfe/tensor3.hh>
+#include <dune/gfe/tensorssd.hh>
+#include <dune/gfe/linearalgebra.hh>
+
+
+namespace Dune::GFE
+{
+
+  /** \brief A function defined by geodesic interpolation
+             from the reference element to a Riemannian manifold.
+
+     \tparam Basis The scalar basis used as weight functions in the Riemannian mean
+     \tparam TS TargetSpace: The manifold that the function takes its values in
+   */
+  template <typename Basis, typename TS>
+  class LocalGeodesicFEFunction
+  {
+  public:
+    using TargetSpace = TS;
+
+  private:
+
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using WeightType = typename LocalBasis::Traits::RangeType;
+    using WeightJacobianType = typename LocalBasis::Traits::JacobianType;
+
+    typedef typename TargetSpace::ctype RT;
+
+    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+    static const int embeddedDim = EmbeddedTangentVector::dimension;
+
+    static const int spaceDim = TargetSpace::TangentVector::dimension;
+
+    /** \brief Short-cut to the local basis
+     */
+    const auto& localBasis() const
+    {
+      return localBasisView_.tree().finiteElement().localBasis();
+    }
+
+  public:
+
+    /** \brief The type used for derivatives */
+    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
+
+    /** \brief The type used for derivatives of the gradient with respect to coefficients */
+    typedef Tensor3<RT,embeddedDim,embeddedDim,dim> DerivativeOfGradientWRTCoefficientType;
+
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalGeodesicFEFunction(const Basis basis)
+      : basis_(basis)
+      , localBasisView_(basis_.localView())
+    {}
+
+    /** \brief Copy constructor
+     */
+    LocalGeodesicFEFunction(const LocalGeodesicFEFunction& other)
+      : basis_(other.basis_)
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(other.coefficients_)
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move constructor
+     */
+    LocalGeodesicFEFunction(LocalGeodesicFEFunction&& other)
+      : basis_(std::move(other.basis_))
+      // Do NOT move the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(std::move(other.coefficients_))
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Copy assignment
+     */
+    LocalGeodesicFEFunction& operator=(const LocalGeodesicFEFunction& other)
+    {
+      basis_ = other.basis_;
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = other.coefficients_;
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move assignment
+     */
+    LocalGeodesicFEFunction& operator=(LocalGeodesicFEFunction&& other)
+    {
+      basis_ = std::move(other.basis_);
+      // Do NOT move the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = std::move(other.coefficients_);
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Bind the local function to a grid element and a set of coefficients
+     *
+     * \param element The grid element
+     * \param coefficients Values to be interpolated
+     */
+    void bind(const Element& element,
+              const std::vector<TargetSpace>& coefficients)
+    {
+      localBasisView_.bind(element);
+      coefficients_ = coefficients;
+
+      assert(localBasisView_.tree().finiteElement().size() == coefficients.size());
+    }
+
+    /** \brief Rebind the FEFunction to another TargetSpace */
+    template<class U>
+    struct rebind
+    {
+      using other = LocalGeodesicFEFunction<Basis,U>;
+    };
+
+    /** \brief The number of coefficients */
+    unsigned int size() const
+    {
+      return localBasisView_.tree().finiteElement().size();
+    }
+
+    /** \brief The type of the reference element */
+    Dune::GeometryType type() const
+    {
+      return localBasisView_.element().type();
+    }
+
+    /** \brief The scalar finite element basis used as interpolation weights
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
+    /** \brief The scalar finite element used as the interpolation weights
+     *
+     * \note This method was added for InterpolationDerivatives, which needs it
+     * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
+     * number type.  This is not optimal, because the localFiniteElement
+     * really is an implementation detail of LocalGeodesicFEFunction and
+     * should not be needed just to copy an entire object.  Other non-Euclidean
+     * interpolation rules may not have such a finite element at all.
+     * Therefore, this method may disappear again eventually.
+     */
+    const LocalFiniteElement& localFiniteElement() const
+    {
+      return localBasisView_.tree().finiteElement();
+    }
+
+    /** \brief Evaluate the function */
+    TargetSpace evaluate(const LocalCoordinate& local) const;
+
+    /** \brief Evaluate the derivative of the function */
+    DerivativeType evaluateDerivative(const LocalCoordinate& 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 LocalCoordinate& local,
+                                      const TargetSpace& q) const;
+
+    /** \brief Evaluate the value and the derivative of the interpolation function
+     */
+    std::pair<TargetSpace,DerivativeType> evaluateValueAndDerivative(const LocalCoordinate& local) const;
+
+    /** \brief Evaluate the derivative of the function value with respect to a coefficient */
+    void evaluateDerivativeOfValueWRTCoefficient(const LocalCoordinate& local,
+                                                 int coefficient,
+                                                 Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& derivative) const;
+
+    /** \brief Evaluate the derivative of the function value with respect to a coefficient */
+    void evaluateFDDerivativeOfValueWRTCoefficient(const LocalCoordinate& local,
+                                                   int coefficient,
+                                                   Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& derivative) const;
+
+    /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
+    void evaluateDerivativeOfGradientWRTCoefficient(const LocalCoordinate& local,
+                                                    int coefficient,
+                                                    DerivativeOfGradientWRTCoefficientType& result) const;
+
+    /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
+    void evaluateFDDerivativeOfGradientWRTCoefficient(const LocalCoordinate& local,
+                                                      int coefficient,
+                                                      DerivativeOfGradientWRTCoefficientType& result) const;
+
+    /** \brief Get the i'th base coefficient. */
+    const TargetSpace& coefficient(int i) const
+    {
+      return coefficients_[i];
+    }
+  private:
+
+    static SymmetricMatrix<RT,embeddedDim> pseudoInverse(const SymmetricMatrix<RT,embeddedDim>& dFdq,
+                                                         const TargetSpace& q)
+    {
+      const int shortDim = TargetSpace::TangentVector::dimension;
+
+      // the orthonormal frame
+      Dune::FieldMatrix<RT,shortDim,embeddedDim> O = q.orthonormalFrame();
+
+      // compute A = O dFDq O^T
+      // TODO A really is symmetric, too
+      Dune::FieldMatrix<RT,shortDim,shortDim> A;
+      for (int i=0; i<shortDim; i++)
+        for (int j=0; j<shortDim; j++) {
+          A[i][j] = 0;
+          for (int k=0; k<embeddedDim; k++)
+            for (int l=0; l<embeddedDim; l++)
+              A[i][j] += O[i][k] * ((k>=l) ? dFdq(k,l) : dFdq(l,k)) *O[j][l];
+        }
+
+      A.invert();
+
+      SymmetricMatrix<RT,embeddedDim> result;
+      result = 0.0;
+      for (int i=0; i<embeddedDim; i++)
+        for (int j=0; j<=i; j++)
+          for (int k=0; k<shortDim; k++)
+            for (int l=0; l<shortDim; l++)
+              result(i,j) += O[k][i]*A[k][l]*O[l][j];
+
+      return result;
+    }
+
+    /** \brief Compute derivate of F(w,q) (the derivative of the weighted distance fctl) wrt to w */
+    Dune::Matrix<RT> computeDFdw(const TargetSpace& q) const
+    {
+      Dune::Matrix<RT> dFdw(embeddedDim,localBasisView_.tree().finiteElement().size());
+      for (size_t i=0; i<localBasisView_.tree().finiteElement().size(); i++) {
+        Dune::FieldVector<RT,embeddedDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
+        for (int j=0; j<embeddedDim; j++)
+          dFdw[j][i] = tmp[j];
+      }
+      return dFdw;
+    }
+
+    Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> computeDqDqF(const std::vector<WeightType>& w, const TargetSpace& q) const
+    {
+      Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> result;
+      result = Tensor3<RT,embeddedDim,embeddedDim,embeddedDim>(RT(0));
+      for (size_t i=0; i<w.size(); i++)
+        result.axpy(w[i][0], TargetSpace::thirdDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i],q));
+      return result;
+    }
+
+    /** \brief The scalar basis that implements the weights for the
+     * Riemannian center of mass
+     */
+    const Basis basis_;
+
+    /** \brief The scalar local finite element, which provides the weighting factors
+     */
+    typename Basis::LocalView localBasisView_;
+
+    /** \brief The coefficient vector */
+    std::vector<TargetSpace> coefficients_;
+
+  };
+
+  template <typename Basis, typename TargetSpace>
+  TargetSpace LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluate(const LocalCoordinate& local) const
+  {
+    // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+    std::vector<WeightType> w;
+    localBasis().evaluateFunction(local,w);
+
+    // The energy functional whose mimimizer is the value of the geodesic interpolation
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+
+    // Create a reasonable initial iterate for the iterative solver
+    GFE::LocalQuickAndDirtyFEFunction<dim,typename Basis::GridView::ctype,LocalFiniteElement,TargetSpace> localProjectedFEFunction(localBasisView_.tree().finiteElement(), coefficients_);
+    TargetSpace initialIterate = localProjectedFEFunction.evaluate(local);
+
+    // Iteratively solve the GFE minimization problem
+    TargetSpaceRiemannianTRSolver<TargetSpace> solver;
+
+    solver.setup(&assembler,
+                 initialIterate,
+                 1e-14,    // tolerance
+                 100       // maxNewtonSteps
+                 );
+
+    solver.solve();
+
+    return solver.getSol();
+  }
+
+  template <typename Basis, typename TargetSpace>
+  typename LocalGeodesicFEFunction<Basis,TargetSpace>::DerivativeType
+  LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateDerivative(const LocalCoordinate& 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 <typename Basis, typename TargetSpace>
+  typename LocalGeodesicFEFunction<Basis,TargetSpace>::DerivativeType
+  LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateDerivative(const LocalCoordinate& local, const TargetSpace& q) const
+  {
+    DerivativeType result;
+
+    // ////////////////////////////////////////////////////////////////////////
+    //  The derivative is evaluated using the implicit function theorem.
+    //  Hence we need to solve a small system of linear equations.
+    // ////////////////////////////////////////////////////////////////////////
+
+    // the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
+    std::vector<WeightJacobianType> B(coefficients_.size());
+    localBasis().evaluateJacobian(local, B);
+
+    // compute negative derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
+    Dune::Matrix<RT> dFdw = computeDFdw(q);
+    dFdw *= -1;
+
+    // multiply the two previous matrices: the result is the right hand side
+    // RHS = dFdw * B;
+    // We need to write this multiplication by hand, because B is actually an (#coefficients) times 1 times dim matrix,
+    // and we need to ignore the middle index.
+    Dune::Matrix<RT> RHS(dFdw.N(), dim);
+
+    for (size_t i=0; i<RHS.N(); i++)
+      for (size_t j=0; j<RHS.M(); j++) {
+        RHS[i][j] = 0;
+        for (size_t k=0; k<dFdw.M(); k++)
+          RHS[i][j] += dFdw[i][k]*B[k][0][j];
+      }
+
+    // the actual system matrix
+    std::vector<WeightType> w;
+    localBasis().evaluateFunction(local, w);
+
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+
+    SymmetricMatrix<RT,embeddedDim> dFdq;
+    assembler.assembleEmbeddedHessian(q,dFdq);
+
+    // We want to solve
+    //     dFdq * x = rhs
+    // We use the Gram-Schmidt solver because we know that dFdq is rank-deficient.
+
+    const int shortDim = TargetSpace::TangentVector::dimension;
+
+    // the orthonormal frame
+    Dune::FieldMatrix<RT,shortDim,embeddedDim> basis = q.orthonormalFrame();
+    GramSchmidtSolver<RT,shortDim,embeddedDim> gramSchmidtSolver(dFdq, basis);
+
+    for (std::size_t i=0; i<dim; i++) {
+
+      Dune::FieldVector<RT,embeddedDim> rhs;
+      for (int j=0; j<embeddedDim; j++)
+        rhs[j] = RHS[j][i];
+
+      Dune::FieldVector<RT,embeddedDim> x;
+      gramSchmidtSolver.solve(x, rhs);
+
+      for (int j=0; j<embeddedDim; j++)
+        result[j][i] = x[j];
+
+    }
+
+    return result;
+  }
+
+  template <typename Basis, typename TargetSpace>
+  std::pair<TargetSpace,typename LocalGeodesicFEFunction<Basis,TargetSpace>::DerivativeType>
+  LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateValueAndDerivative(const LocalCoordinate& local) const
+  {
+    std::pair<TargetSpace,DerivativeType> result;
+    result.first = evaluate(local);
+    result.second = evaluateDerivative(local,result.first);
+    return result;
+  }
+
+  template <typename Basis, typename TargetSpace>
+  void LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateDerivativeOfValueWRTCoefficient(const LocalCoordinate& local,
+                                          int coefficient,
+                                          Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& result) const
+  {
+    // the function value at the point where we are evaluating the derivative
+    TargetSpace q = evaluate(local);
+
+    // dFdq
+    std::vector<WeightType> w;
+    localBasis().evaluateFunction(local,w);
+
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+
+    SymmetricMatrix<RT,embeddedDim> dFdq;
+    assembler.assembleEmbeddedHessian(q,dFdq);
+
+    const int shortDim = TargetSpace::TangentVector::dimension;
+
+    // the orthonormal frame
+    Dune::FieldMatrix<RT,shortDim,embeddedDim> O = q.orthonormalFrame();
+
+    // compute A = O dFDq O^T
+    Dune::FieldMatrix<RT,shortDim,shortDim> A;
+    for (int i=0; i<shortDim; i++)
+      for (int j=0; j<shortDim; j++) {
+        A[i][j] = 0;
+        for (int k=0; k<embeddedDim; k++)
+          for (int l=0; l<embeddedDim; l++)
+            A[i][j] += O[i][k] * ((k>=l) ? dFdq(k,l) : dFdq(l,k)) * O[j][l];
+      }
+
+    //
+    Dune::FieldMatrix<RT,embeddedDim,embeddedDim> rhs = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
+    rhs *= -w[coefficient];
+
+    for (int i=0; i<embeddedDim; i++) {
+
+      Dune::FieldVector<RT,shortDim> shortRhs;
+      O.mv(rhs[i],shortRhs);
+
+      Dune::FieldVector<RT,shortDim> shortX;
+      A.solve(shortX,shortRhs);
+
+      O.mtv(shortX,result[i]);
+
+    }
+
+  }
+
+  template <typename Basis, typename TargetSpace>
+  void LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateFDDerivativeOfValueWRTCoefficient(const LocalCoordinate& local,
+                                            int coefficient,
+                                            Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& result) const
+  {
+    double eps = 1e-6;
+
+    // the function value at the point where we are evaluating the derivative
+    const Dune::FieldMatrix<RT,spaceDim,embeddedDim> B = coefficients_[coefficient].orthonormalFrame();
+
+    Dune::FieldMatrix<RT,spaceDim,embeddedDim> interimResult;
+
+    std::vector<TargetSpace> cornersPlus  = coefficients_;
+    std::vector<TargetSpace> cornersMinus = coefficients_;
+
+    for (int j=0; j<spaceDim; j++) {
+
+      typename TargetSpace::EmbeddedTangentVector forwardVariation = B[j];
+      forwardVariation *= eps;
+      typename TargetSpace::EmbeddedTangentVector backwardVariation = B[j];
+      backwardVariation *= -eps;
+
+      cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
+      cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
+
+      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus(basis_), fMinus(basis_);
+      fPlus.bind(localBasisView_.element(),cornersPlus);
+      fMinus.bind(localBasisView_.element(),cornersMinus);
+
+      TargetSpace hPlus  = fPlus.evaluate(local);
+      TargetSpace hMinus = fMinus.evaluate(local);
+
+      interimResult[j]  = hPlus.globalCoordinates();
+      interimResult[j] -= hMinus.globalCoordinates();
+      interimResult[j] /= 2*eps;
+
+    }
+
+    for (int i=0; i<embeddedDim; i++)
+      for (int j=0; j<embeddedDim; j++) {
+        result[i][j] = 0;
+        for (int k=0; k<spaceDim; k++)
+          result[i][j] += B[k][i]*interimResult[k][j];
+      }
+
+  }
+
+
+  template <typename Basis, typename TargetSpace>
+  void LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateDerivativeOfGradientWRTCoefficient(const LocalCoordinate& local,
+                                             int coefficient,
+                                             DerivativeOfGradientWRTCoefficientType& result) const
+  {
+    // the function value at the point where we are evaluating the derivative
+    TargetSpace q = evaluate(local);
+
+    // the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
+    std::vector<WeightJacobianType> BNested(coefficients_.size());
+    localBasis().evaluateJacobian(local, BNested);
+    Dune::Matrix<RT> B(coefficients_.size(), dim);
+    for (size_t i=0; i<coefficients_.size(); i++)
+      for (size_t j=0; j<dim; j++)
+        B[i][j] = BNested[i][0][j];
+
+    // the actual system matrix
+    std::vector<WeightType> w;
+    localBasis().evaluateFunction(local,w);
+
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+
+    /** \todo Use a symmetric matrix here */
+    SymmetricMatrix<RT,embeddedDim> dFdq;
+    assembler.assembleEmbeddedHessian(q,dFdq);
+
+
+    Dune::FieldMatrix<RT,embeddedDim,embeddedDim> mixedDerivative = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
+    TensorSSD<RT,embeddedDim,embeddedDim> dvDwF(coefficients_.size());
+    dvDwF = 0;
+    for (int i=0; i<embeddedDim; i++)
+      for (int j=0; j<embeddedDim; j++)
+        dvDwF(i, j, coefficient) = mixedDerivative[i][j];
+
+
+    // dFDq is not invertible, if the target space is embedded into a higher-dimensional
+    // Euclidean space.  Therefore we use its pseudo inverse.  I don't think that is the
+    // best way, though.
+    SymmetricMatrix<RT,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq,q);
+
+    //
+    Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> dvDqF
+      =  TargetSpace::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(coefficients_[coefficient], q);
+
+    dvDqF = RT(w[coefficient]) * dvDqF;
+
+    // Put it all together
+    // dvq[i][j] = \partial q_j / \partial v_i
+    Dune::FieldMatrix<RT,embeddedDim,embeddedDim> dvq;
+    evaluateDerivativeOfValueWRTCoefficient(local,coefficient,dvq);
+
+    Dune::FieldMatrix<RT, embeddedDim, dim> derivative = evaluateDerivative(local);
+
+    Tensor3<RT, embeddedDim,embeddedDim,embeddedDim> dqdqF;
+    dqdqF = computeDqDqF(w,q);
+
+    TensorSSD<RT, embeddedDim,embeddedDim> dqdwF(coefficients_.size());
+
+    for (size_t k=0; k<coefficients_.size(); k++) {
+      SymmetricMatrix<RT,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
+      for (int i=0; i<embeddedDim; i++)
+        for (int j=0; j<=i; j++)
+          dqdwF(i, j, k) = dqdwF(j, i, k) = hesse(i,j);
+
+    }
+
+    TensorSSD<RT, embeddedDim,embeddedDim> dqdwF_times_dvq(coefficients_.size());
+    for (int i=0; i<embeddedDim; i++)
+      for (int j=0; j<embeddedDim; j++)
+        for (size_t k=0; k<coefficients_.size(); k++) {
+          dqdwF_times_dvq(i, j, k) = 0;
+          for (int l=0; l<embeddedDim; l++)
+            dqdwF_times_dvq(i, j, k) += dqdwF(l, j, k) * dvq[i][l];
+        }
+
+    Tensor3<RT,embeddedDim,embeddedDim,dim> foo;
+    foo =  -1 * dvDqF*derivative - (dvq*dqdqF)*derivative;
+    TensorSSD<RT,embeddedDim,embeddedDim> bar(dim);
+    bar = dvDwF * B + dqdwF_times_dvq*B;
+
+    for (int i=0; i<embeddedDim; i++)
+      for (int j=0; j<embeddedDim; j++)
+        for (std::size_t k=0; k<dim; k++)
+          foo[i][j][k] -= bar(i, j, k);
+
+    result = RT(0);
+    for (int i=0; i<embeddedDim; i++)
+      for (int j=0; j<embeddedDim; j++)
+        for (std::size_t k=0; k<dim; k++)
+          for (int l=0; l<embeddedDim; l++)
+            // TODO Smarter implementation of the product with a symmetric matrix
+            result[i][j][k] += ((j>=l) ? dFdqPseudoInv(j,l) : dFdqPseudoInv(l,j)) * foo[i][l][k];
+
+  }
+
+
+  template <typename Basis, typename TargetSpace>
+  void LocalGeodesicFEFunction<Basis,TargetSpace>::
+  evaluateFDDerivativeOfGradientWRTCoefficient(const LocalCoordinate& local,
+                                               int coefficient,
+                                               DerivativeOfGradientWRTCoefficientType& result) const
+  {
+    double eps = 1e-6;
+
+    // loop over the different partial derivatives
+    for (int j=0; j<embeddedDim; j++) {
+
+      std::vector<TargetSpace> cornersPlus  = coefficients_;
+      std::vector<TargetSpace> cornersMinus = coefficients_;
+      typename TargetSpace::CoordinateType aPlus  = coefficients_[coefficient].globalCoordinates();
+      typename TargetSpace::CoordinateType aMinus = coefficients_[coefficient].globalCoordinates();
+      aPlus[j]  += eps;
+      aMinus[j] -= eps;
+      cornersPlus[coefficient]  = TargetSpace(aPlus);
+      cornersMinus[coefficient] = TargetSpace(aMinus);
+
+      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus(basis_), fMinus(basis_);
+      fPlus.bind(localBasisView_.element(),cornersPlus);
+      fMinus.bind(localBasisView_.element(),cornersMinus);
+
+      const auto hPlus  = fPlus.evaluateDerivative(local);
+      const auto hMinus = fMinus.evaluateDerivative(local);
+
+      result[j]  = hPlus;
+      result[j] -= hMinus;
+      result[j] /= 2*eps;
+
+    }
+
+    for (int j=0; j<embeddedDim; j++) {
+
+      TargetSpace q = evaluate(local);
+      Dune::FieldVector<RT,embeddedDim> foo;
+      for (std::size_t l=0; l<dim; l++) {
+
+        for (int k=0; k<embeddedDim; k++)
+          foo[k] = result[j][k][l];
+
+        foo = q.projectOntoTangentSpace(foo);
+
+        for (int k=0; k<embeddedDim; k++)
+          result[j][k][l] = foo[k];
+
+      }
+
+    }
+
+  }
+
+
+  /** \brief A function defined by geodesic interpolation
+             from the reference element to a ProductManifold<RealTuple,Rotation>.
+
+     This is a specialization for speeding up the code.
+
+     \tparam Basis The scalar basis used as weight functions in the Riemannian mean
+     \tparam field_type Type used for point coefficients in the target space
+   */
+  template <typename Basis, typename field_type>
+  class LocalGeodesicFEFunction<Basis,GFE::ProductManifold<GFE::RealTuple<field_type,3>, GFE::Rotation<field_type,3> > >
+  {
+    using TargetSpace = GFE::ProductManifold<GFE::RealTuple<field_type,3>, GFE::Rotation<field_type,3> >;
+
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using WeightType = typename LocalBasis::Traits::RangeType;
+    using WeightJacobianType = typename LocalBasis::Traits::JacobianType;
+
+    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+    static const int embeddedDim = EmbeddedTangentVector::dimension;
+
+    static const int spaceDim = TargetSpace::TangentVector::dimension;
+
+    /** \brief Short-cut to the local basis
+     */
+    const auto& localBasis() const
+    {
+      return localBasisView_.tree().finiteElement().localBasis();
+    }
+
+  public:
+
+    /** \brief The type used for derivatives */
+    typedef Dune::FieldMatrix<field_type, embeddedDim, dim> DerivativeType;
+
+    /** \brief The type used for derivatives of the gradient with respect to coefficients */
+    typedef Tensor3<field_type,embeddedDim,embeddedDim,dim> DerivativeOfGradientWRTCoefficientType;
+
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalGeodesicFEFunction(const Basis basis)
+      : basis_(basis)
+      , localBasisView_(basis_)
+      , orientationFEFunction_(basis)
+    {}
+
+    /** \brief Copy constructor
+     */
+    LocalGeodesicFEFunction(const LocalGeodesicFEFunction& other)
+      : basis_(other.basis_)
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(other.coefficients_)
+      , translationCoefficients_(other.translationCoefficients_)
+      , orientationFEFunction_(other.orientationFEFunction_)
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move constructor
+     */
+    LocalGeodesicFEFunction(LocalGeodesicFEFunction&& other)
+      : basis_(std::move(other.basis_))
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(std::move(other.coefficients_))
+      , translationCoefficients_(std::move(other.translationCoefficients_))
+      , orientationFEFunction_(std::move(other.orientationFEFunction_))
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Copy assignment
+     */
+    LocalGeodesicFEFunction& operator=(const LocalGeodesicFEFunction& other)
+    {
+      basis_ = other.basis_;
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = other.coefficients_;
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+      translationCoefficients_ = other.translationCoefficients_;
+      orientationFEFunction_ = other.orientationFEFunction_;
+    }
+
+    /** \brief Move assignment
+     */
+    LocalGeodesicFEFunction& operator=(LocalGeodesicFEFunction&& other)
+    {
+      basis_ = std::move(other.basis_);
+      // Do NOT move the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = std::move(other.coefficients_);
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+      translationCoefficients_ = std::move(other.translationCoefficients_);
+      orientationFEFunction_ = std::move(other.orientationFEFunction_);
+    }
+
+    /** \brief Bind the function to a particular weight function set and coefficients
+     */
+    void bind(const Element& element,
+              const std::vector<TargetSpace>& coefficients)
+    {
+      using namespace Dune::Indices;
+
+      localBasisView_.bind(element);
+      coefficients_ = coefficients;
+      assert(localBasis().size() == coefficients.size());
+
+      translationCoefficients_.resize(coefficients.size());
+      for (size_t i=0; i<coefficients.size(); i++)
+        translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
+
+      std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
+      for (size_t i=0; i<coefficients.size(); i++)
+        orientationCoefficients[i] = coefficients[i][_1];
+
+      orientationFEFunction_.bind(element,orientationCoefficients);
+    }
+
+    /** \brief Rebind the FEFunction to another TargetSpace */
+    template<class U>
+    struct rebind
+    {
+      using other = LocalGeodesicFEFunction<Basis,U>;
+    };
+
+    /** \brief The number of Lagrange points */
+    unsigned int size() const
+    {
+      return localBasisView_.size();
+    }
+
+    /** \brief The type of the reference element */
+    Dune::GeometryType type() const
+    {
+      return localBasisView_.element().type();
+    }
+
+    /** \brief The scalar finite element used as the interpolation weights
+     *
+     * \note This method was added for InterpolationDerivatives, which needs it
+     * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
+     * number type.  This is not optimal, because the localFiniteElement
+     * really is an implementation detail of LocalGeodesicFEFunction and
+     * should not be needed just to copy an entire object.  Other non-Euclidean
+     * interpolation rules may not have such a finite element at all.
+     * Therefore, this method may disappear again eventually.
+     */
+    const LocalFiniteElement& localFiniteElement() const
+    {
+      return localBasisView_.tree().finiteElement();
+    }
+
+    /** \brief Evaluate the function */
+    TargetSpace evaluate(const LocalCoordinate& local) const
+    {
+      using namespace Dune::Indices;
+
+      TargetSpace result;
+
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<WeightType> w;
+      localBasis().evaluateFunction(local,w);
+
+      result[_0] = Dune::FieldVector<field_type,3>(0.0);
+      for (size_t i=0; i<w.size(); i++)
+        result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);
+
+      result[_1] = orientationFEFunction_.evaluate(local);
+      return result;
+    }
+
+    /** \brief Evaluate the derivative of the function */
+    DerivativeType evaluateDerivative(const LocalCoordinate& local) const
+    {
+      DerivativeType result(0);
+
+      // get translation part
+      std::vector<WeightJacobianType> sfDer(translationCoefficients_.size());
+      localBasis().evaluateJacobian(local, sfDer);
+
+      for (size_t i=0; i<translationCoefficients_.size(); i++)
+        for (int j=0; j<3; j++)
+          result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
+
+      // get orientation part
+      Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_.evaluateDerivative(local);
+      for (int i=0; i<4; i++)
+        for (std::size_t j=0; j<dim; j++)
+          result[3+i][j] = qResult[i][j];
+
+      return result;
+    }
+
+    /** \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 LocalCoordinate& local,
+                                      const TargetSpace& q) const
+    {
+      using namespace Dune::Indices;
+
+      DerivativeType result(0);
+
+      // get translation part
+      std::vector<WeightJacobianType> sfDer(translationCoefficients_.size());
+      localBasis().evaluateJacobian(local, sfDer);
+
+      for (size_t i=0; i<translationCoefficients_.size(); i++)
+        for (int j=0; j<3; j++)
+          result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
+
+      // get orientation part
+      Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_.evaluateDerivative(local,q[_1]);
+      for (int i=0; i<4; i++)
+        for (std::size_t j=0; j<dim; j++)
+          result[3+i][j] = qResult[i][j];
+
+      return result;
+    }
+
+    /** \brief Evaluate the value and the derivative of the interpolation function
+     */
+    std::pair<TargetSpace,DerivativeType> evaluateValueAndDerivative(const LocalCoordinate& local) const
+    {
+      std::pair<TargetSpace,DerivativeType> result;
+      result.first = evaluate(local);
+      result.second = evaluateDerivative(local,result.first);
+      return result;
+    }
+
+    /** \brief Evaluate the derivative of the function value with respect to a coefficient */
+    void evaluateDerivativeOfValueWRTCoefficient(const LocalCoordinate& local,
+                                                 int coefficient,
+                                                 Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& derivative) const
+    {
+      derivative = 0;
+
+      // Translation part
+      std::vector<WeightType> w;
+      localBasis().evaluateFunction(local,w);
+      for (int i=0; i<3; i++)
+        derivative[i][i] = w[coefficient];
+
+      // Rotation part
+      Dune::FieldMatrix<field_type,4,4> qDerivative;
+      orientationFEFunction_->evaluateDerivativeOfValueWRTCoefficient(local,coefficient,qDerivative);
+      for (int i=0; i<4; i++)
+        for (int j=0; j<4; j++)
+          derivative[3+i][3+j] = qDerivative[i][j];
+    }
+
+    /** \brief Evaluate the derivative of the function value with respect to a coefficient */
+    void evaluateFDDerivativeOfValueWRTCoefficient(const LocalCoordinate& local,
+                                                   int coefficient,
+                                                   Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& derivative) const
+    {
+      derivative = 0;
+
+      // Translation part
+      std::vector<WeightType> w;
+      localBasis().evaluateFunction(local,w);
+      for (int i=0; i<3; i++)
+        derivative[i][i] = w[coefficient];
+
+      // Rotation part
+      Dune::FieldMatrix<field_type,4,4> qDerivative;
+      orientationFEFunction_->evaluateFDDerivativeOfValueWRTCoefficient(local,coefficient,qDerivative);
+      for (int i=0; i<4; i++)
+        for (int j=0; j<4; j++)
+          derivative[3+i][3+j] = qDerivative[i][j];
+    }
+
+    /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
+    void evaluateDerivativeOfGradientWRTCoefficient(const LocalCoordinate& local,
+                                                    int coefficient,
+                                                    DerivativeOfGradientWRTCoefficientType& derivative) const
+    {
+      derivative = field_type(0);
+
+      // Translation part
+      std::vector<WeightType> w;
+      localBasis().evaluateJacobian(local,w);
+      for (int i=0; i<3; i++)
+        derivative[i][i] = w[coefficient][0];
+
+      // Rotation part
+      Tensor3<field_type,4,4,dim> qDerivative;
+      orientationFEFunction_->evaluateDerivativeOfGradientWRTCoefficient(local,coefficient,qDerivative);
+      for (int i=0; i<4; i++)
+        for (int j=0; j<4; j++)
+          for (int k=0; k<dim; k++)
+            derivative[3+i][3+j][k] = qDerivative[i][j][k];
+    }
+
+    /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
+    void evaluateFDDerivativeOfGradientWRTCoefficient(const LocalCoordinate& local,
+                                                      int coefficient,
+                                                      DerivativeOfGradientWRTCoefficientType& derivative) const
+    {
+      derivative = 0;
+
+      // Translation part
+      std::vector<WeightType> w;
+      localBasis().evaluateJacobian(local,w);
+      for (int i=0; i<3; i++)
+        derivative[i][i] = w[coefficient][0];
+
+      // Rotation part
+      Tensor3<field_type,4,4,dim> qDerivative;
+      orientationFEFunction_->evaluateFDDerivativeOfGradientWRTCoefficient(local,coefficient,qDerivative);
+      for (int i=0; i<4; i++)
+        for (int j=0; j<4; j++)
+          for (int k=0; k<dim; k++)
+            derivative[3+i][3+j][k] = qDerivative[i][j][k];
+    }
+
+    const TargetSpace& coefficient(int i) const {
+      return coefficients_[i];
+    }
+
+  private:
+
+    /** \brief The scalar basis that implements the weights for the
+     * Riemannian center of mass
+     */
+    const Basis basis_;
+
+    /** \brief The scalar local finite element, which provides the weighting factors
+     */
+    typename Basis::LocalView localBasisView_;
+
+    // The coefficients of this interpolation rule
+    std::vector<TargetSpace> coefficients_;
+
+    // The coefficients again.  Yes, we store it twice:  To evaluate the interpolation rule efficiently
+    // we need access to the coefficients for the various factor spaces separately.
+    std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;
+
+    LocalGeodesicFEFunction<Basis,Rotation<field_type,3> > orientationFEFunction_;
+  };
+
+}  // namespace Dune::GFE
+
+#endif

dune/gfe/functions/localisometrycomponentfunction.hh

+#ifndef DUNE_GFE_FUNCTIONS_LOCALISOMETRYCOMPONENTFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_LOCALISOMETRYCOMPONENTFUNCTION_HH
+
+#include <dune/gfe/bendingisometryhelper.hh>
+
+namespace Dune::GFE::Impl
+{
+  /**
+   * \brief Local isometry component function that can be used as a Dune::Functions::DifferentiableFunction
+   *
+   * \tparam ctype Type used for coordinates on the reference element
+   * \tparam LocalInterpolation Local function mapping into a product manifold of type RealTuple x Rotations
+   */
+  template<class ctype, class LocalInterpolation>
+  class LocalIsometryComponentFunction
+  {
+  public:
+    LocalIsometryComponentFunction(LocalInterpolation& localInterpolation,
+                                   const int component)
+      : localInterpolation_(localInterpolation),
+      component_(component)
+    {}
+
+    auto operator() (const Dune::FieldVector<ctype,2>& x) const
+    {
+      return localInterpolation_.evaluate(x)[Dune::Indices::_0].globalCoordinates()[component_];
+    }
+
+    /**
+     * \brief Obtain derivative function
+     */
+    friend auto derivative(const LocalIsometryComponentFunction& p)
+    {
+      return [&p](const Dune::FieldVector<ctype,2>& x)
+             {
+               auto derivativeQuaternion = p.localInterpolation_.evaluate(x)[Dune::Indices::_1];
+               auto derivativeMatrix = Rotation<ctype,3>::quaternionToMatrix(derivativeQuaternion);
+               auto reducedDerivativeMatrix = cancelLastMatrixColumn(derivativeMatrix);
+               return reducedDerivativeMatrix[p.component_];
+             };
+    }
+    LocalInterpolation& localInterpolation_;
+    const int component_;
+  };
+
+} // end namespace Dune::GFE::Impl
+
+#endif   //#ifndef DUNE_GFE_FUNCTIONS_LOCALISOMETRYCOMPONENTFUNCTION_HH

dune/gfe/functions/localprojectedfefunction.hh

+#ifndef DUNE_GFE_FUNCTIONS_LOCALPROJECTEDFEFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_LOCALPROJECTEDFEFUNCTION_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/type.hh>
+
+#include <dune/gfe/linearalgebra.hh>
+#include <dune/gfe/polardecomposition.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+namespace Dune::GFE
+{
+
+  /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold
+   *
+   * \tparam Basis The basis used for the interpolation in the embedding space.
+   * \tparam TS The manifold that the function takes its values in
+   * \tparam conforming Geometrical conformity of the functions (omits projections when set to false)
+   *
+   * \note Even though the embedding space is typically more than one-dimensional,
+   * we allow this basis to be scalar-valued.
+   */
+  template <class Basis, class TS, bool conforming=true>
+  class LocalProjectedFEFunction
+  {
+  public:
+    using TargetSpace=TS;
+  private:
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using ScalarFERange = typename LocalBasis::Traits::RangeType;
+    using ScalarFEJacobian = typename LocalBasis::Traits::JacobianType;
+
+    typedef typename TargetSpace::ctype RT;
+
+    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+    static const int embeddedDim = EmbeddedTangentVector::dimension;
+
+    /** \brief Short-cut to the local basis
+     */
+    const auto& localBasis() const
+    {
+      return localBasisView_.tree().finiteElement().localBasis();
+    }
+
+  public:
+
+    /** \brief The type used for derivatives */
+    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
+
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalProjectedFEFunction(const Basis& basis)
+      : basis_(basis)
+      , localBasisView_(basis_)
+    {}
+
+    /** \brief Copy constructor
+     */
+    LocalProjectedFEFunction(const LocalProjectedFEFunction& other)
+      : basis_(other.basis_)
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(other.coefficients_)
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move constructor
+     */
+    LocalProjectedFEFunction(LocalProjectedFEFunction&& other)
+      : basis_(std::move(other.basis_))
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(std::move(other.coefficients_))
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Copy assignment
+     */
+    LocalProjectedFEFunction& operator=(const LocalProjectedFEFunction& other)
+    {
+      basis_ = other.basis_;
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = other.coefficients_;
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move assignment
+     */
+    LocalProjectedFEFunction& operator=(LocalProjectedFEFunction&& other)
+    {
+      basis_ = std::move(other.basis_);
+      // Do NOT move the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = std::move(other.coefficients_);
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Bind the local function to a grid element and a set of coefficients
+     *
+     * \param element The grid element
+     * \param coefficients Values to be interpolated
+     */
+    void bind(const Element& element,
+              const std::vector<TargetSpace>& coefficients)
+    {
+      localBasisView_.bind(element);
+      coefficients_ = coefficients;
+
+      assert(localBasisView_.tree().finiteElement().size() == coefficients.size());
+    }
+
+    /** \brief Rebind the FEFunction to another TargetSpace */
+    template<class U>
+    struct rebind
+    {
+      using other = LocalProjectedFEFunction<Basis,U,conforming>;
+    };
+
+    /** \brief The number of Lagrange points */
+    unsigned int size() const
+    {
+      return localBasisView_.size();
+    }
+
+    /** \brief The type of the reference element */
+    Dune::GeometryType type() const
+    {
+      return localBasisView_.element().type();
+    }
+
+    /** \brief The finite element basis used to interpolate in the embedding space
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
+    /** \brief The scalar finite element used as the interpolation weights
+     *
+     * \note This method was added for InterpolationDerivatives, which needs it
+     * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
+     * number type.  This is not optimal, because the localFiniteElement
+     * really is an implementation detail of LocalGeodesicFEFunction and
+     * should not be needed just to copy an entire object.  Other non-Euclidean
+     * interpolation rules may not have such a finite element at all.
+     * Therefore, this method may disappear again eventually.
+     */
+    const LocalFiniteElement& localFiniteElement() const
+    {
+      return localBasisView_.tree().finiteElement();
+    }
+
+    /** \brief Evaluate the function */
+    auto evaluate(const LocalCoordinate& local) const;
+
+    /** \brief Evaluate the derivative of the function */
+    DerivativeType evaluateDerivative(const LocalCoordinate& 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!
+     *
+     * \note This method is only usable in the conforming setting, because it requires the caller
+     * to hand over the interpolation value as a TargetSpace object.
+     */
+    DerivativeType evaluateDerivative(const LocalCoordinate& local,
+                                      const TargetSpace& q) const;
+
+    /** \brief Evaluate the value and the derivative of the interpolation function
+     *
+     * \return A std::pair containing the value and the first derivative of the interpolation function.
+     * If the interpolation is conforming then the first member of the pair will be a TargetSpace.
+     * Otherwise it will be a RealTuple.
+     */
+    auto evaluateValueAndDerivative(const LocalCoordinate& local) const;
+
+    /** \brief Get the i'th base coefficient. */
+    const TargetSpace& coefficient(int i) const
+    {
+      return coefficients_[i];
+    }
+  private:
+
+    /** \brief The finite element basis for the interpolation in the embedding space
+     */
+    const Basis basis_;
+
+    /** \brief The local finite element used for the interpolation in the embedding space
+     */
+    typename Basis::LocalView localBasisView_;
+
+    /** \brief The coefficient vector */
+    std::vector<TargetSpace> coefficients_;
+
+  };
+
+  template <class Basis, class TargetSpace, bool conforming>
+  auto LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
+  evaluate(const LocalCoordinate& local) const
+  {
+    // Compute value of FE function in embedding space
+    std::vector<ScalarFERange> w;
+    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());
+
+    // Possibly project onto target space
+    if constexpr (conforming)
+      return TargetSpace::projectOnto(c);
+    else
+      return (RealTuple<RT, TargetSpace::CoordinateType::dimension>)c;
+
+  }
+
+  template <class Basis, class TargetSpace,bool conforming>
+  typename LocalProjectedFEFunction<Basis,TargetSpace,conforming>::DerivativeType
+  LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
+  evaluateDerivative(const LocalCoordinate& local) const
+  {
+    if constexpr(conforming)
+    {
+      // the function value at the point where we are evaluating the derivative
+      TargetSpace q = evaluate(local);
+
+      // Actually compute the derivative
+      return evaluateDerivative(local, q);
+    }
+    else {
+      std::vector<ScalarFEJacobian> wDer;
+      localBasis().evaluateJacobian(local, wDer);
+
+      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];
+
+      return derivative;
+    }
+  }
+
+  template <class Basis, class TargetSpace,bool conforming>
+  typename LocalProjectedFEFunction<Basis,TargetSpace,conforming>::DerivativeType
+  LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
+  evaluateDerivative(const LocalCoordinate& local, const TargetSpace& q) const
+  {
+    // Compute value and derivative of FE approximation
+    std::vector<ScalarFERange> w;
+    localBasis().evaluateFunction(local,w);
+
+    std::vector<ScalarFEJacobian> wDer;
+    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);
+
+    return derivativeOfProjection*derivative;
+  }
+
+  template <class Basis, class TargetSpace,bool conforming>
+  auto
+  LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
+  evaluateValueAndDerivative(const LocalCoordinate& local) const
+  {
+    // Construct the type of the result -- it depends on whether the interpolation
+    // is conforming or not.
+    using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;
+
+    std::pair<Value,DerivativeType> result;
+
+    ///////////////////////////////////////////////////////////
+    //  Compute the value of the interpolation function
+    ///////////////////////////////////////////////////////////
+
+    std::vector<ScalarFERange> w;
+    localBasis().evaluateFunction(local,w);
+
+    typename TargetSpace::CoordinateType embeddedInterpolation(0);
+    for (size_t i=0; i<coefficients_.size(); i++)
+      embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());
+
+    if constexpr (conforming)
+      result.first = TargetSpace::projectOnto(embeddedInterpolation);
+    else
+      result.first = (RealTuple<RT, TargetSpace::CoordinateType::dimension>)embeddedInterpolation;
+
+    ///////////////////////////////////////////////////////////
+    //  Compute the derivative of the interpolation function
+    ///////////////////////////////////////////////////////////
+
+    // Compute the interpolation in the surrounding space
+    std::vector<ScalarFEJacobian> wDer;
+    localBasis().evaluateJacobian(local,wDer);
+
+    result.second = 0.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++)
+          result.second[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];
+
+    // The derivative of the projection onto the manifold
+    if constexpr(conforming)
+      result.second = TargetSpace::derivativeOfProjection(embeddedInterpolation) * result.second;
+
+    return result;
+  }
+
+
+  /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for SO(3)
+   *
+   * \tparam Basis The basis used for the interpolation in the embedding space.
+   * \tparam field_type The number type used for function values
+   * \tparam conforming Geometrical conformity of the functions (omits projections when set to false)
+   *
+   * \note Even though the embedding space is 9-dimensional, we allow this basis
+   * to be scalar-valued.
+   */
+  template <class Basis, class field_type, bool conforming>
+  class LocalProjectedFEFunction<Basis,Rotation<field_type,3>,conforming>
+  {
+  public:
+    typedef Rotation<field_type,3> TargetSpace;
+  private:
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using ScalarFERange = typename LocalBasis::Traits::RangeType;
+    using ScalarFEJacobian = typename LocalBasis::Traits::JacobianType;
+
+    typedef typename TargetSpace::ctype RT;
+
+    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+    static const int embeddedDim = EmbeddedTangentVector::dimension;
+
+    /**
+     * \param A The argument of the projection
+     * \param polar The image of the projection, i.e., the polar factor of A
+     */
+    static std::array<std::array<FieldMatrix<field_type,3,3>, 3>, 3> derivativeOfProjection(const FieldMatrix<field_type,3,3>& A,
+                                                                                            FieldMatrix<field_type,3,3>& polar)
+    {
+      std::array<std::array<FieldMatrix<field_type,3,3>, 3>, 3> result;
+
+      for (int i=0; i<3; i++)
+        for (int j=0; j<3; j++)
+          for (int k=0; k<3; k++)
+            for (int l=0; l<3; l++)
+              result[i][j][k][l] = (i==k) and (j==l);
+
+      polar = A;
+
+      // Use Heron's method
+      const size_t maxIterations = 100;
+      double tol = 0.001;
+      for (size_t iteration=0; iteration<maxIterations; iteration++)
+      {
+        auto oldPolar = polar;
+        auto polarInvert = polar;
+        polarInvert.invert();
+        for (size_t i=0; i<polar.N(); i++)
+          for (size_t j=0; j<polar.M(); j++)
+            polar[i][j] = 0.5 * (polar[i][j] + polarInvert[j][i]);
+
+        // Alternative name to align the code better with a description in a math text
+        const auto& dQT = result;
+
+        // Multiply from the right with Q^{-T}
+        decltype(result) tmp2;
+        for (int i=0; i<3; i++)
+          for (int j=0; j<3; j++)
+            for (int k=0; k<3; k++)
+              for (int l=0; l<3; l++)
+                tmp2[i][j][k][l] = 0.0;
+
+        for (int i=0; i<3; i++)
+          for (int j=0; j<3; j++)
+            for (int k=0; k<3; k++)
+              for (int l=0; l<3; l++)
+                for (int m=0; m<3; m++)
+                  for (int o=0; o<3; o++)
+                    tmp2[i][j][k][l] += polarInvert[m][i] * dQT[o][m][k][l] * polarInvert[j][o];
+
+        for (int i=0; i<3; i++)
+          for (int j=0; j<3; j++)
+            for (int k=0; k<3; k++)
+              for (int l=0; l<3; l++)
+                result[i][j][k][l] = 0.5 * (result[i][j][k][l] - tmp2[i][j][k][l]);
+
+        oldPolar -= polar;
+        if (oldPolar.frobenius_norm() < tol) {
+          break;
+        }
+      }
+
+      return result;
+    }
+
+    /** \brief Short-cut to the local basis
+     */
+    const auto& localBasis() const
+    {
+      return localBasisView_.tree().finiteElement().localBasis();
+    }
+
+  public:
+
+    /** \brief The type used for derivatives */
+    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
+
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalProjectedFEFunction(const Basis basis)
+      : basis_(basis)
+      , localBasisView_(basis_)
+    {}
+
+    /** \brief Copy constructor
+     */
+    LocalProjectedFEFunction(const LocalProjectedFEFunction& other)
+      : basis_(other.basis_)
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(other.coefficients_)
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move constructor
+     */
+    LocalProjectedFEFunction(LocalProjectedFEFunction&& other)
+      : basis_(std::move(other.basis_))
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(std::move(other.coefficients_))
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Copy assignment
+     */
+    LocalProjectedFEFunction& operator=(const LocalProjectedFEFunction& other)
+    {
+      basis_ = other.basis_;
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = other.coefficients_;
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move assignment
+     */
+    LocalProjectedFEFunction& operator=(LocalProjectedFEFunction&& other)
+    {
+      basis_ = std::move(other.basis_);
+      // Do NOT move the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = std::move(other.coefficients_);
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Bind the function to a specific finite element interpolation space and set of coefficients
+     *
+     * \param element The grid element to bind to
+     * \param coefficients Coefficients of the finite element function on that element
+     */
+    void bind(const Element& element,
+              const std::vector<TargetSpace>& coefficients)
+    {
+      localBasisView_.bind(element);
+      coefficients_ = coefficients;
+
+      assert(localBasisView_.size() == coefficients.size());
+    }
+
+    /** \brief Rebind the FEFunction to another TargetSpace */
+    template<class U>
+    struct rebind
+    {
+      using other = LocalProjectedFEFunction<Basis,U>;
+    };
+
+    /** \brief The number of Lagrange points */
+    unsigned int size() const
+    {
+      return localBasisView_.tree().finiteElement().size();
+    }
+
+    /** \brief The type of the reference element */
+    Dune::GeometryType type() const
+    {
+      return localBasisView_.element().type();
+    }
+
+    /** \brief The finite element basis used to interpolate in the embedding space
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
+    /** \brief The scalar finite element used as the interpolation weights
+     *
+     * \note This method was added for InterpolationDerivatives, which needs it
+     * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
+     * number type.  This is not optimal, because the localFiniteElement
+     * really is an implementation detail of LocalGeodesicFEFunction and
+     * should not be needed just to copy an entire object.  Other non-Euclidean
+     * interpolation rules may not have such a finite element at all.
+     * Therefore, this method may disappear again eventually.
+     */
+    const LocalFiniteElement& localFiniteElement() const
+    {
+      return localBasisView_.tree().finiteElement();
+    }
+
+    /** \brief Evaluate the function */
+    TargetSpace evaluate(const LocalCoordinate& local) const
+    {
+      Rotation<field_type,3> result;
+
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<ScalarFERange> w;
+      localBasis().evaluateFunction(local,w);
+
+      // Interpolate in R^{3x3}
+      FieldMatrix<field_type,3,3> interpolatedMatrix(0);
+      for (size_t i=0; i<coefficients_.size(); i++)
+      {
+        FieldMatrix<field_type,3,3> coefficientAsMatrix;
+        coefficients_[i].matrix(coefficientAsMatrix);
+        interpolatedMatrix.axpy(w[i][0], coefficientAsMatrix);
+      }
+
+      // Project back onto SO(3)
+      result.set(Dune::GFE::PolarDecomposition()(interpolatedMatrix));
+
+      return result;
+    }
+
+    /** \brief Evaluate the derivative of the function */
+    DerivativeType evaluateDerivative(const LocalCoordinate& 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);
+    }
+
+    /** \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 function at 'local'.  If you provide something wrong here the result will be wrong, too!
+     */
+    DerivativeType evaluateDerivative(const LocalCoordinate& local,
+                                      const TargetSpace& q) const
+    {
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<ScalarFERange> w;
+      localBasis().evaluateFunction(local,w);
+
+      std::vector<ScalarFEJacobian> wDer;
+      localBasis().evaluateJacobian(local,wDer);
+
+      // Compute matrix representations for all coefficients (we only have them in quaternion representation)
+      std::vector<Dune::FieldMatrix<field_type,3,3> > coefficientsAsMatrix(coefficients_.size());
+      for (size_t i=0; i<coefficients_.size(); i++)
+        coefficients_[i].matrix(coefficientsAsMatrix[i]);
+
+      // Interpolate in R^{3x3}
+      FieldMatrix<field_type,3,3> interpolatedMatrix(0);
+      for (size_t i=0; i<coefficients_.size(); i++)
+        interpolatedMatrix.axpy(w[i][0], coefficientsAsMatrix[i]);
+
+      Tensor3<RT,dim,3,3> derivative(0);
+
+      for (size_t dir=0; dir<dim; dir++)
+        for (size_t i=0; i<3; i++)
+          for (size_t j=0; j<3; j++)
+            for (size_t k=0; k<coefficients_.size(); k++)
+              derivative[dir][i][j] += wDer[k][0][dir] * coefficientsAsMatrix[k][i][j];
+
+      FieldMatrix<field_type,3,3> polarFactor;
+      auto derivativeOfProjection = this->derivativeOfProjection(interpolatedMatrix,polarFactor);
+
+      Tensor3<field_type,dim,3,3> intermediateResult(0);
+
+      for (size_t dir=0; dir<dim; dir++)
+        for (size_t i=0; i<3; i++)
+          for (size_t j=0; j<3; j++)
+            for (size_t k=0; k<3; k++)
+              for (size_t l=0; l<3; l++)
+                intermediateResult[dir][i][j] += derivativeOfProjection[i][j][k][l]*derivative[dir][k][l];
+
+      // One more application of the chain rule: we need to go from orthogonal matrices to quaternions
+      Tensor3<field_type,4,3,3> derivativeOfMatrixToQuaternion = Rotation<field_type,3>::derivativeOfMatrixToQuaternion(polarFactor);
+
+      DerivativeType result(0);
+
+      for (size_t dir0=0; dir0<4; dir0++)
+        for (size_t dir1=0; dir1<dim; dir1++)
+          for (size_t i=0; i<3; i++)
+            for (size_t j=0; j<3; j++)
+              result[dir0][dir1] += derivativeOfMatrixToQuaternion[dir0][i][j] * intermediateResult[dir1][i][j];
+
+      return result;
+    }
+
+    /** \brief Evaluate the value and the derivative of the interpolation function
+     *
+     * \return A std::pair containing the value and the first derivative of the interpolation function.
+     * If the interpolation is conforming then the first member of the pair will be a TargetSpace.
+     * Otherwise it will be a RealTuple.
+     */
+    auto evaluateValueAndDerivative(const LocalCoordinate& local) const
+    {
+      // Construct the type of the result -- it depends on whether the interpolation
+      // is conforming or not.
+      using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;
+
+      std::pair<Value,DerivativeType> result;
+
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<ScalarFERange> w;
+      localBasis().evaluateFunction(local,w);
+
+      // Interpolate in R^{3x3}
+      FieldMatrix<field_type,3,3> interpolatedMatrix(0);
+      for (size_t i=0; i<coefficients_.size(); i++)
+      {
+        FieldMatrix<field_type,3,3> coefficientAsMatrix;
+        coefficients_[i].matrix(coefficientAsMatrix);
+        interpolatedMatrix.axpy(w[i][0], coefficientAsMatrix);
+      }
+
+      // Project back onto SO(3)
+      static_assert(conforming, "Nonconforming interpolation into SO(3) is not implemented!");
+      result.first.set(Dune::GFE::PolarDecomposition()(interpolatedMatrix));
+
+      result.second = evaluateDerivative(local, result.first);
+
+      return result;
+    }
+
+    /** \brief Get the i'th base coefficient. */
+    const TargetSpace& coefficient(int i) const
+    {
+      return coefficients_[i];
+    }
+  private:
+
+    /** \brief The finite element basis for the interpolation in the embedding space
+     */
+    const Basis basis_;
+
+    /** \brief The scalar local finite element, for interpolating in the embedding space
+     */
+    typename Basis::LocalView localBasisView_;
+
+    /** \brief The coefficient vector */
+    std::vector<TargetSpace> coefficients_;
+
+  };
+
+
+  /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for R^3 x SO(3)
+   *
+   * \tparam Basis The basis used for the interpolation in the embedding space.
+   * \tparam field_type The number type used for function values
+   *
+   * \note Even though the embedding space is 12-dimensional, we allow this basis
+   * to be scalar-valued.
+   */
+  template <class Basis, class field_type>
+  class LocalProjectedFEFunction<Basis,ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
+  {
+  public:
+    using TargetSpace = ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >;
+  private:
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using ScalarFERange = typename LocalBasis::Traits::RangeType;
+    using ScalarFEJacobian = typename LocalBasis::Traits::JacobianType;
+
+    typedef typename TargetSpace::ctype RT;
+
+    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+    static const int embeddedDim = EmbeddedTangentVector::dimension;
+
+    /** \brief Short-cut to the local basis
+     */
+    const auto& localBasis() const
+    {
+      return localBasisView_.tree().finiteElement().localBasis();
+    }
+
+  public:
+
+    /** \brief The type used for derivatives */
+    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
+
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalProjectedFEFunction(const Basis basis)
+      : basis_(basis)
+      , localBasisView_(basis_)
+      , orientationFunction_(basis)
+    {}
+
+    /** \brief Copy constructor
+     */
+    LocalProjectedFEFunction(const LocalProjectedFEFunction& other)
+      : basis_(other.basis_)
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(other.coefficients_)
+      , translationCoefficients_(other.translationCoefficients_)
+      , orientationFunction_(other.orientationFunction_)
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Move constructor
+     */
+    LocalProjectedFEFunction(LocalProjectedFEFunction&& other)
+      : basis_(std::move(other.basis_))
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      , localBasisView_(basis_.localView())
+      , coefficients_(std::move(other.coefficients_))
+      , translationCoefficients_(std::move(other.translationCoefficients_))
+      , orientationFunction_(std::move(other.orientationFunction_))
+    {
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+    }
+
+    /** \brief Copy assignment
+     */
+    LocalProjectedFEFunction& operator=(const LocalProjectedFEFunction& other)
+    {
+      basis_ = other.basis_;
+      // Do NOT copy the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = other.coefficients_;
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+      translationCoefficients_ = other.translationCoefficients_;
+      orientationFunction_ = other.orientationFunction_;
+    }
+
+    /** \brief Move assignment
+     */
+    LocalProjectedFEFunction& operator=(LocalProjectedFEFunction&& other)
+    {
+      basis_ = std::move(other.basis_);
+      // Do NOT move the localView: It contains a reference to other.basis_.
+      localBasisView_ = basis_.localView();
+      coefficients_ = std::move(other.coefficients_);
+      if (other.localBasisView_.bound())
+        localBasisView_.bind(other.localBasisView_.element());
+      translationCoefficients_ = std::move(other.translationCoefficients_);
+      orientationFunction_ = std::move(other.orientationFunction_);
+    }
+
+    /** \brief Bind the function to a specific finite element interpolation space and set of coefficients
+     *
+     * \param element The grid element to bind to
+     * \param coefficients Coefficients of the finite element function
+     */
+    void bind(const Element& element,
+              const std::vector<TargetSpace>& coefficients)
+    {
+      using namespace Dune::Indices;
+
+      localBasisView_.bind(element);
+      coefficients_ = coefficients;
+
+      assert(localBasis().size() == coefficients.size());
+
+      translationCoefficients_.resize(coefficients.size());
+      for (size_t i=0; i<coefficients.size(); i++)
+        translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
+
+      std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
+      for (size_t i=0; i<coefficients.size(); i++)
+        orientationCoefficients[i] = coefficients[i][_1];
+
+      orientationFunction_.bind(element,orientationCoefficients);
+    }
+
+    /** \brief Rebind the FEFunction to another TargetSpace */
+    template<class U>
+    struct rebind
+    {
+      using other = LocalProjectedFEFunction<Basis,U>;
+    };
+
+    /** \brief The number of Lagrange points */
+    unsigned int size() const
+    {
+      return localBasisView_.size();
+    }
+
+    /** \brief The type of the reference element */
+    Dune::GeometryType type() const
+    {
+      return localBasisView_.element().type();
+    }
+
+    /** \brief The finite element basis used to interpolate in the embedding space
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
+    /** \brief Evaluate the function */
+    TargetSpace evaluate(const LocalCoordinate& local) const
+    {
+      using namespace Dune::Indices;
+
+      TargetSpace result;
+
+      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+      std::vector<ScalarFERange> w;
+      localBasis().evaluateFunction(local,w);
+
+      result[_0] = FieldVector<field_type,3>(0.0);
+      for (size_t i=0; i<w.size(); i++)
+        result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);
+
+      result[_1] = orientationFunction_.evaluate(local);
+
+      return result;
+    }
+
+    /** \brief Evaluate the derivative of the function */
+    DerivativeType evaluateDerivative(const LocalCoordinate& 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);
+    }
+
+    /** \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 function at 'local'.  If you provide something wrong here the result will be wrong, too!
+     */
+    DerivativeType evaluateDerivative(const LocalCoordinate& local,
+                                      const TargetSpace& q) const
+    {
+      using namespace Dune::Indices;
+
+      DerivativeType result(0);
+
+      // get translation part
+      std::vector<ScalarFEJacobian> sfDer(translationCoefficients_.size());
+      localBasis().evaluateJacobian(local, sfDer);
+
+      for (size_t i=0; i<translationCoefficients_.size(); i++)
+        for (int j=0; j<3; j++)
+          result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
+
+      // get orientation part
+      Dune::FieldMatrix<field_type,4,dim> qResult = orientationFunction_.evaluateDerivative(local,q[_1]);
+
+      for (int i=0; i<4; i++)
+        for (std::size_t j=0; j<dim; j++)
+          result[3+i][j] = qResult[i][j];
+
+      return result;
+    }
+
+    /** \brief Get the i'th base coefficient. */
+    const TargetSpace& coefficient(int i) const
+    {
+      return coefficients_[i];
+    }
+  private:
+
+    /** \brief The finite element basis for the interpolation in the embedding space
+     */
+    const Basis basis_;
+
+    /** \brief The local finite element, for the interpolation in the embedding space
+     */
+    typename Basis::LocalView localBasisView_;
+
+    // The coefficients of this interpolation rule
+    std::vector<TargetSpace> coefficients_;
+
+    // The coefficients again.  Yes, we store it twice:  To evaluate the interpolation rule efficiently
+    // we need access to the coefficients for the various factor spaces separately.
+    std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;
+
+    LocalProjectedFEFunction<Basis,Rotation<field_type, 3> > orientationFunction_;
+
+
+  };
+
+}  // namespace Dune::GFE
+
+#endif

dune/gfe/functions/localquickanddirtyfefunction.hh

+#ifndef DUNE_GFE_FUNCTIONS_LOCALQUICKANDDIRTYFEFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_LOCALQUICKANDDIRTYFEFUNCTION_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/type.hh>
+
+#include <dune/gfe/linearalgebra.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+namespace Dune::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 consistently 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;
+
+    /** \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);
+
+    // Special-casing for the Rotations.  This is quite a challenge:
+    //
+    // Projection-based interpolation in the
+    // unit quaternions is not a good idea, here, because you may end up on the
+    // wrong sheet of the covering of SO(3) by the unit quaternions.  The minimization
+    // solver for the weighted distance functional may then actually converge
+    // to the wrong local minimizer -- the interpolating functions wraps "the wrong way"
+    // around SO(3).
+    //
+    // Projection-based interpolation in SO(3) is not a good idea, because it is about as
+    // expensive as the geodesic interpolation itself.  That's not good if all we want
+    // is an initial iterate.
+    //
+    // Averaging in R^{3x3} followed by QR decomposition is not a good idea either:
+    // The averaging can give you singular matrices quite quickly (try two matrices that
+    // differ by a rotation of pi/2).  Then, the QR factorization of the identity may
+    // not be the identity (but differ in sign)!  I tried that with the implementation
+    // from Numerical Recipes.
+    //
+    // What do we do instead?  We start from the coefficient that produces the lowest
+    // value of the objective functional.
+    if constexpr (std::is_same<TargetSpace,Rotation<double,3> >::value)
+    {
+      AverageDistanceAssembler averageDistance(coefficients_,w);
+
+      auto minDistance = averageDistance.value(coefficients_[0]);
+      std::size_t minCoefficient = 0;
+
+      for (std::size_t i=1; i<coefficients_.size(); i++)
+      {
+        auto distance = averageDistance.value(coefficients_[i]);
+        if (distance < minDistance)
+        {
+          minDistance = distance;
+          minCoefficient = i;
+        }
+      }
+
+      return coefficients_[minCoefficient];
+    }
+    else
+    {
+      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);
+    }
+  }
+
+}   // namespace Dune::GFE
+#endif

dune/gfe/geodesicdifference.hh

@@ -5,23 +5,29 @@
 
 #include <dune/istl/bvector.hh>
 
-template <class TargetSpace>
-Dune::BlockVector<typename TargetSpace::TangentVector> computeGeodesicDifference(const std::vector<TargetSpace>& a,
-                                                                                 const std::vector<TargetSpace>& b)
+
+namespace Dune::GFE
 {
-  if (a.size() != b.size())
-    DUNE_THROW(Dune::Exception, "a and b have to have the same length!");
 
-  Dune::BlockVector<typename TargetSpace::TangentVector> result(a.size());
+  template <class TargetSpace>
+  Dune::BlockVector<typename TargetSpace::TangentVector> computeGeodesicDifference(const std::vector<TargetSpace>& a,
+                                                                                   const std::vector<TargetSpace>& b)
+  {
+    if (a.size() != b.size())
+      DUNE_THROW(Dune::Exception, "a and b have to have the same length!");
+
+    Dune::BlockVector<typename TargetSpace::TangentVector> result(a.size());
+
+    for (size_t i=0; i<result.size(); i++) {
 
-  for (size_t i=0; i<result.size(); i++) {
+      // Subtract orientations on the tangent space of 'a'
+      result[i] = TargetSpace::difference(a[i], b[i]);
 
-    // Subtract orientations on the tangent space of 'a'
-    result[i] = TargetSpace::difference(a[i], b[i]);
+    }
 
+    return result;
   }
 
-  return result;
-}
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/geodesicfefunctionadaptor.hh

@@ -6,206 +6,212 @@
 
 #include "localgeodesicfefunction.hh"
 
-template <class Basis, class TargetSpace>
-class GeodesicFEFunctionAdaptor
+
+namespace Dune::GFE
 {
-public:
-  /** \brief Refine a grid globally and prolong a given geodesic finite element function
-   */
-  template <class GridType>
-  static void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
+
+  template <class Basis, class TargetSpace>
+  class GeodesicFEFunctionAdaptor
   {
-    const int dim = GridType::dimension;
+  public:
+    /** \brief Refine a grid globally and prolong a given geodesic finite element function
+     */
+    template <class GridType>
+    static void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
+    {
+      const int dim = GridType::dimension;
 
-    assert(x.size() == grid.size(dim));
+      assert(x.size() == grid.size(dim));
 
-    // /////////////////////////////////////////////////////
-    //   Save leaf p1 data in a map
-    // /////////////////////////////////////////////////////
+      // /////////////////////////////////////////////////////
+      //   Save leaf p1 data in a map
+      // /////////////////////////////////////////////////////
 
-    const typename GridType::Traits::LocalIdSet&   idSet    = grid.localIdSet();
-    const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
+      const typename GridType::Traits::LocalIdSet&   idSet    = grid.localIdSet();
+      const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
 
-    std::map<typename GridType::Traits::LocalIdSet::IdType, TargetSpace> dofMap;
+      std::map<typename GridType::Traits::LocalIdSet::IdType, TargetSpace> dofMap;
 
-    for (const auto& vertex : vertices(grid.leafGridView()))
-      dofMap.insert(std::make_pair(idSet.id(vertex), x[indexSet.index(vertex)]));
+      for (const auto& vertex : vertices(grid.leafGridView()))
+        dofMap.insert(std::make_pair(idSet.id(vertex), x[indexSet.index(vertex)]));
 
 
 
-    // /////////////////////////////////////////////////////
-    //   Globally refine the grid
-    // /////////////////////////////////////////////////////
+      // /////////////////////////////////////////////////////
+      //   Globally refine the grid
+      // /////////////////////////////////////////////////////
 
-    grid.globalRefine(1);
+      grid.globalRefine(1);
 
 
-    // /////////////////////////////////////////////////////
-    //   Restore and interpolate the data
-    // /////////////////////////////////////////////////////
+      // /////////////////////////////////////////////////////
+      //   Restore and interpolate the data
+      // /////////////////////////////////////////////////////
 
-    P1NodalBasis<typename GridType::LeafGridView> p1Basis(grid.leafGridView());
-    x.resize(grid.size(dim));
+      P1NodalBasis<typename GridType::LeafGridView> p1Basis(grid.leafGridView());
+      x.resize(grid.size(dim));
 
-    for (const auto& element : elements(grid.leafGridView())) {
+      for (const auto& element : elements(grid.leafGridView())) {
 
-      // Set up a local gfe function on the father element
-      size_t nFatherDofs = element.father().subEntities(dim);
-      std::vector<TargetSpace> coefficients(nFatherDofs);
+        // Set up a local gfe function on the father element
+        size_t nFatherDofs = element.father().subEntities(dim);
+        std::vector<TargetSpace> coefficients(nFatherDofs);
 
-      for (int i=0; i<nFatherDofs; i++)
-        coefficients[i] = dofMap.find(idSet.subId(element.father(),i,dim))->second;
+        for (int i=0; i<nFatherDofs; i++)
+          coefficients[i] = dofMap.find(idSet.subId(element.father(),i,dim))->second;
 
-      typedef typename P1NodalBasis<typename GridType::LeafGridView>::LocalFiniteElement LocalFiniteElement;
-      LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fatherFunction(p1Basis.getLocalFiniteElement(element),
-                                                                                        coefficients);
+        typedef typename P1NodalBasis<typename GridType::LeafGridView>::LocalFiniteElement LocalFiniteElement;
+        LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fatherFunction(p1Basis.getLocalFiniteElement(element),
+                                                                                          coefficients);
 
-      // The embedding of this element into the father geometry
-      const auto& geometryInFather = element.geometryInFather();
+        // The embedding of this element into the father geometry
+        const auto& geometryInFather = element.geometryInFather();
 
-      size_t nDofs = element.subEntities(dim);
-      for (int i=0; i<nDofs; i++) {
+        size_t nDofs = element.subEntities(dim);
+        for (int i=0; i<nDofs; i++) {
 
-        if (dofMap.find(idSet.subId(element,i,dim)) != dofMap.end()) {
+          if (dofMap.find(idSet.subId(element,i,dim)) != dofMap.end()) {
 
-          // If the vertex exists on the coarser level we take the value from there.
-          // That should be faster and more accurate than interpolating
-          x[indexSet.subIndex(element,i,dim)] = dofMap[idSet.subId(element,i,dim)];
+            // If the vertex exists on the coarser level we take the value from there.
+            // That should be faster and more accurate than interpolating
+            x[indexSet.subIndex(element,i,dim)] = dofMap[idSet.subId(element,i,dim)];
 
-        } else {
+          } else {
 
-          // Interpolate
-          x[indexSet.subIndex(element,i,dim)] = fatherFunction.evaluate(geometryInFather.corner(i));
+            // Interpolate
+            x[indexSet.subIndex(element,i,dim)] = fatherFunction.evaluate(geometryInFather.corner(i));
+
+          }
 
         }
 
       }
 
-    }
 
+    }
 
-  }
 
+    /** \brief Coordinate function in one variable, constant in the others
 
-  /** \brief Coordinate function in one variable, constant in the others
+        This is used to extract the positions of the Lagrange nodes.
+     */
+    template <int dim>
+    struct CoordinateFunction
+      : public Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,1> >
+    {
+      CoordinateFunction(int d)
+        : d_(d)
+      {}
 
-      This is used to extract the positions of the Lagrange nodes.
-   */
-  template <int dim>
-  struct CoordinateFunction
-    : public Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,1> >
-  {
-    CoordinateFunction(int d)
-      : d_(d)
-    {}
+      void evaluate(const Dune::FieldVector<double, dim>& x, Dune::FieldVector<double,1>& out) const {
+        out[0] = x[d_];
+      }
 
-    void evaluate(const Dune::FieldVector<double, dim>& x, Dune::FieldVector<double,1>& out) const {
-      out[0] = x[d_];
-    }
+      //
+      int d_;
+    };
 
-    //
-    int d_;
-  };
 
+    /** \brief Refine a grid globally and prolong a given geodesic finite element function
+     *
+     * \tparam order Interpolation order of the function space.  Kinda stupid that I
+     *  have to provide this by hand.  Will change...
+     */
+    template <int order>
+    static void higherOrderGFEFunctionAdaptor(Basis& basis,
+                                              typename Basis::GridView::Grid& grid, std::vector<TargetSpace>& x)
+    {
+      typedef typename Basis::GridView::Grid GridType;
+      const int dim = GridType::dimension;
 
-  /** \brief Refine a grid globally and prolong a given geodesic finite element function
-   *
-   * \tparam order Interpolation order of the function space.  Kinda stupid that I
-   *  have to provide this by hand.  Will change...
-   */
-  template <int order>
-  static void higherOrderGFEFunctionAdaptor(Basis& basis,
-                                            typename Basis::GridView::Grid& grid, std::vector<TargetSpace>& x)
-  {
-    typedef typename Basis::GridView::Grid GridType;
-    const int dim = GridType::dimension;
+      typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
 
-    typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
+      // /////////////////////////////////////////////////////
+      //   Save leaf p1 data in a map
+      // /////////////////////////////////////////////////////
 
-    // /////////////////////////////////////////////////////
-    //   Save leaf p1 data in a map
-    // /////////////////////////////////////////////////////
+      const typename GridType::Traits::LocalIdSet&   idSet    = grid.localIdSet();
 
-    const typename GridType::Traits::LocalIdSet&   idSet    = grid.localIdSet();
+      typedef typename GridType::Traits::LocalIdSet::IdType IdType;
+      std::map<IdType, std::vector<TargetSpace> > dofMap;
 
-    typedef typename GridType::Traits::LocalIdSet::IdType IdType;
-    std::map<IdType, std::vector<TargetSpace> > dofMap;
+      assert(x.size() == basis.size());
 
-    assert(x.size() == basis.size());
+      ElementIterator eIt    = grid.template leafbegin<0>();
+      ElementIterator eEndIt = grid.template leafend<0>();
 
-    ElementIterator eIt    = grid.template leafbegin<0>();
-    ElementIterator eEndIt = grid.template leafend<0>();
+      for (; eIt!=eEndIt; ++eIt) {
 
-    for (; eIt!=eEndIt; ++eIt) {
+        const typename Basis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(*eIt);
+        std::vector<TargetSpace> coefficients(lfe.localCoefficients().size());
 
-      const typename Basis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(*eIt);
-      std::vector<TargetSpace> coefficients(lfe.localCoefficients().size());
+        for (size_t i=0; i<lfe.localCoefficients().size(); i++)
+          coefficients[i] = x[basis.index(*eIt, i)];
 
-      for (size_t i=0; i<lfe.localCoefficients().size(); i++)
-        coefficients[i] = x[basis.index(*eIt, i)];
+        IdType id = idSet.id(*eIt);
+        dofMap.insert(std::make_pair(id, coefficients));
 
-      IdType id = idSet.id(*eIt);
-      dofMap.insert(std::make_pair(id, coefficients));
+      }
 
-    }
 
+      // /////////////////////////////////////////////////////
+      //   Globally refine the grid
+      // /////////////////////////////////////////////////////
 
-    // /////////////////////////////////////////////////////
-    //   Globally refine the grid
-    // /////////////////////////////////////////////////////
+      grid.globalRefine(1);
 
-    grid.globalRefine(1);
 
+      // /////////////////////////////////////////////////////
+      //   Restore and interpolate the data
+      // /////////////////////////////////////////////////////
 
-    // /////////////////////////////////////////////////////
-    //   Restore and interpolate the data
-    // /////////////////////////////////////////////////////
+      basis.update(grid.leafGridView());
 
-    basis.update(grid.leafGridView());
+      x.resize(basis.size());
 
-    x.resize(basis.size());
+      for (eIt=grid.template leafbegin<0>(); eIt!=eEndIt; ++eIt) {
 
-    for (eIt=grid.template leafbegin<0>(); eIt!=eEndIt; ++eIt) {
+        const typename Basis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(*eIt);
 
-      const typename Basis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(*eIt);
+        typedef typename Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::FiniteElementType FatherFiniteElementType;
 
-      typedef typename Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::FiniteElementType FatherFiniteElementType;
+        auto fatherLFE = std::unique_ptr<FatherFiniteElementType>(Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::create(eIt->father()->type()));
 
-      auto fatherLFE = std::unique_ptr<FatherFiniteElementType>(Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::create(eIt->father()->type()));
+        // Set up a local gfe function on the father element
+        std::vector<TargetSpace> coefficients = dofMap[idSet.id(*eIt->father())];
 
-      // Set up a local gfe function on the father element
-      std::vector<TargetSpace> coefficients = dofMap[idSet.id(*eIt->father())];
+        LocalGeodesicFEFunction<dim,double,typename Basis::LocalFiniteElement,TargetSpace> fatherFunction(*fatherLFE, coefficients);
 
-      LocalGeodesicFEFunction<dim,double,typename Basis::LocalFiniteElement,TargetSpace> fatherFunction(*fatherLFE, coefficients);
+        // The embedding of this element into the father geometry
+        const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
 
-      // The embedding of this element into the father geometry
-      const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
+        // Generate position of the Lagrange nodes
+        std::vector<Dune::FieldVector<double,dim> > lagrangeNodes(lfe.localBasis().size());
 
-      // Generate position of the Lagrange nodes
-      std::vector<Dune::FieldVector<double,dim> > lagrangeNodes(lfe.localBasis().size());
+        for (int i=0; i<dim; i++) {
+          CoordinateFunction<dim> lFunction(i);
+          std::vector<Dune::FieldVector<double,1> > coordinates;
+          lfe.localInterpolation().interpolate(lFunction, coordinates);
 
-      for (int i=0; i<dim; i++) {
-        CoordinateFunction<dim> lFunction(i);
-        std::vector<Dune::FieldVector<double,1> > coordinates;
-        lfe.localInterpolation().interpolate(lFunction, coordinates);
+          for (size_t j=0; j<coordinates.size(); j++)
+            lagrangeNodes[j][i] = coordinates[j];
 
-        for (size_t j=0; j<coordinates.size(); j++)
-          lagrangeNodes[j][i] = coordinates[j];
+        }
 
-      }
+        for (int i=0; i<lfe.localCoefficients().size(); i++) {
 
-      for (int i=0; i<lfe.localCoefficients().size(); i++) {
+          unsigned int idx = basis.index(*eIt, i);
 
-        unsigned int idx = basis.index(*eIt, i);
+          x[idx] = fatherFunction.evaluate(geometryInFather.global(lagrangeNodes[i]));
 
-        x[idx] = fatherFunction.evaluate(geometryInFather.global(lagrangeNodes[i]));
+        }
 
       }
 
     }
 
-  }
+  };
 
-};
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/geometries/moebiusstrip.hh

@@ -7,7 +7,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   /// \brief Functor representing a Möbius strip in 3D, where the base circle is the circle in the x-y-plane with radius_ around 0,0,0
   template <class T = double>
@@ -101,6 +102,6 @@ namespace Dune {
     return analyticGridFunction<Grid>(MoebiusStripProjection<T>{radius});
   }
 
-} // end namespace Dune
+}  // namespace Dune::GFE
 
 #endif // DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH

dune/gfe/geometries/twistedstrip.hh

@@ -7,7 +7,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   /// \brief Functor representing a twisted strip in 3D, with length length_ and nTwists twists
   template <class T = double>
@@ -99,6 +100,6 @@ namespace Dune {
     return analyticGridFunction<Grid>(TwistedStripProjection<T>{length, nTwists});
   }
 
-} // end namespace Dune
+}  // namespace Dune::GFE
 
 #endif // DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH

dune/gfe/gfedifferencenormsquared.hh

@@ -15,274 +15,280 @@
 
 #include <dune/gfe/globalgeodesicfefunction.hh>
 
-template <class Basis, class TargetSpace>
-class GFEDifferenceNormSquared {
-
-  typedef typename Basis::GridView GridView;
-  typedef typename Basis::GridView::Grid GridType;
-
-  // Grid dimension
-  const static int dim = GridType::dimension;
-
-  /** \brief Check whether given local coordinates are contained in the reference element
-      \todo This method exists in the Dune grid interface!  But we need the eps.
-   */
-  static bool checkInside(const Dune::GeometryType& type,
-                          const Dune::FieldVector<double, dim> &loc,
-                          double eps)
-  {
-    switch (type.dim()) {
-
-    case 0 :    // vertex
-      return false;
-    case 1 :    // line
-      return -eps <= loc[0] && loc[0] <= 1+eps;
-    case 2 :
-
-      if (type.isSimplex()) {
-        return -eps <= loc[0] && -eps <= loc[1] && (loc[0]+loc[1])<=1+eps;
-      } else if (type.isCube()) {
-        return -eps <= loc[0] && loc[0] <= 1+eps
-               && -eps <= loc[1] && loc[1] <= 1+eps;
-      } else {
-        DUNE_THROW(Dune::GridError, "checkInside():  ERROR:  Unknown type " << type << " found!");
-      }
 
-    case 3 :
-      if (type.isSimplex()) {
-        return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
-               && (loc[0]+loc[1]+loc[2]) <= 1+eps;
-      } else if (type.isPyramid()) {
-        return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
-               && (loc[0]+loc[2]) <= 1+eps
-               && (loc[1]+loc[2]) <= 1+eps;
-      } else if (type.isPrism()) {
-        return -eps <= loc[0] && -eps <= loc[1]
-               && (loc[0]+loc[1])<= 1+eps
-               && -eps <= loc[2] && loc[2] <= 1+eps;
-      } else if (type.isCube()) {
-        return -eps <= loc[0] && loc[0] <= 1+eps
-               && -eps <= loc[1] && loc[1] <= 1+eps
-               && -eps <= loc[2] && loc[2] <= 1+eps;
-      }else {
+namespace Dune::GFE
+{
+
+  template <class Basis, class TargetSpace>
+  class GFEDifferenceNormSquared {
+
+    typedef typename Basis::GridView GridView;
+    typedef typename Basis::GridView::Grid GridType;
+
+    // Grid dimension
+    const static int dim = GridType::dimension;
+
+    /** \brief Check whether given local coordinates are contained in the reference element
+        \todo This method exists in the Dune grid interface!  But we need the eps.
+     */
+    static bool checkInside(const Dune::GeometryType& type,
+                            const Dune::FieldVector<double, dim> &loc,
+                            double eps)
+    {
+      switch (type.dim()) {
+
+      case 0 :  // vertex
+        return false;
+      case 1 :  // line
+        return -eps <= loc[0] && loc[0] <= 1+eps;
+      case 2 :
+
+        if (type.isSimplex()) {
+          return -eps <= loc[0] && -eps <= loc[1] && (loc[0]+loc[1])<=1+eps;
+        } else if (type.isCube()) {
+          return -eps <= loc[0] && loc[0] <= 1+eps
+                 && -eps <= loc[1] && loc[1] <= 1+eps;
+        } else {
+          DUNE_THROW(Dune::GridError, "checkInside():  ERROR:  Unknown type " << type << " found!");
+        }
+
+      case 3 :
+        if (type.isSimplex()) {
+          return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
+                 && (loc[0]+loc[1]+loc[2]) <= 1+eps;
+        } else if (type.isPyramid()) {
+          return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
+                 && (loc[0]+loc[2]) <= 1+eps
+                 && (loc[1]+loc[2]) <= 1+eps;
+        } else if (type.isPrism()) {
+          return -eps <= loc[0] && -eps <= loc[1]
+                 && (loc[0]+loc[1])<= 1+eps
+                 && -eps <= loc[2] && loc[2] <= 1+eps;
+        } else if (type.isCube()) {
+          return -eps <= loc[0] && loc[0] <= 1+eps
+                 && -eps <= loc[1] && loc[1] <= 1+eps
+                 && -eps <= loc[2] && loc[2] <= 1+eps;
+        }else {
+          DUNE_THROW(Dune::GridError, "checkInside():  ERROR:  Unknown type " << type << " found!");
+        }
+      default :
         DUNE_THROW(Dune::GridError, "checkInside():  ERROR:  Unknown type " << type << " found!");
       }
-    default :
-      DUNE_THROW(Dune::GridError, "checkInside():  ERROR:  Unknown type " << type << " found!");
+
     }
 
-  }
 
+    /** \brief Coordinate function in one variable, constant in the others
 
-  /** \brief Coordinate function in one variable, constant in the others
+       This is used to extract the positions of the Lagrange nodes.
+     */
+    struct CoordinateFunction
+      : public Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,1> >
+    {
+      CoordinateFunction(int d)
+        : d_(d)
+      {}
 
-     This is used to extract the positions of the Lagrange nodes.
-   */
-  struct CoordinateFunction
-    : public Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,1> >
-  {
-    CoordinateFunction(int d)
-      : d_(d)
-    {}
+      void evaluate(const Dune::FieldVector<double, dim>& x, Dune::FieldVector<double,1>& out) const {
+        out[0] = x[d_];
+      }
 
-    void evaluate(const Dune::FieldVector<double, dim>& x, Dune::FieldVector<double,1>& out) const {
-      out[0] = x[d_];
-    }
+      //
+      int d_;
+    };
 
-    //
-    int d_;
-  };
+  public:
 
-public:
+    static void interpolate(const GridType& sourceGrid,
+                            const std::vector<TargetSpace>& source,
+                            const GridType& targetGrid,
+                            std::vector<TargetSpace>& target)
+    {
+      // Create a leaf function, which we need to be able to call 'evalall()'
+      Basis sourceBasis(sourceGrid.leafGridView());
+      GlobalGeodesicFEFunction<Basis,TargetSpace> sourceFunction(sourceBasis, source);
 
-  static void interpolate(const GridType& sourceGrid,
-                          const std::vector<TargetSpace>& source,
-                          const GridType& targetGrid,
-                          std::vector<TargetSpace>& target)
-  {
-    // Create a leaf function, which we need to be able to call 'evalall()'
-    Basis sourceBasis(sourceGrid.leafGridView());
-    GlobalGeodesicFEFunction<Basis,TargetSpace> sourceFunction(sourceBasis, source);
+      Basis targetBasis(targetGrid.leafGridView());
 
-    Basis targetBasis(targetGrid.leafGridView());
+      // ///////////////////////////////////////////////////////////////////////////////////////////
+      //   Prolong the adaptive solution onto the uniform grid in order to make it comparable
+      // ///////////////////////////////////////////////////////////////////////////////////////////
 
-    // ///////////////////////////////////////////////////////////////////////////////////////////
-    //   Prolong the adaptive solution onto the uniform grid in order to make it comparable
-    // ///////////////////////////////////////////////////////////////////////////////////////////
+      target.resize(targetBasis.size());
 
-    target.resize(targetBasis.size());
+      // handle each dof only once
+      Dune::BitSetVector<1> handled(targetBasis.size(), false);
 
-    // handle each dof only once
-    Dune::BitSetVector<1> handled(targetBasis.size(), false);
+      typename GridView::template Codim<0>::Iterator eIt    = targetGrid.template leafbegin<0>();
+      typename GridView::template Codim<0>::Iterator eEndIt = targetGrid.template leafend<0>();
 
-    typename GridView::template Codim<0>::Iterator eIt    = targetGrid.template leafbegin<0>();
-    typename GridView::template Codim<0>::Iterator eEndIt = targetGrid.template leafend<0>();
+      // Loop over all dofs by looping over all elements
+      for (; eIt!=eEndIt; ++eIt) {
 
-    // Loop over all dofs by looping over all elements
-    for (; eIt!=eEndIt; ++eIt) {
+        const typename Basis::LocalFiniteElement& targetLFE = targetBasis.getLocalFiniteElement(*eIt);
 
-      const typename Basis::LocalFiniteElement& targetLFE = targetBasis.getLocalFiniteElement(*eIt);
+        // Generate position of the Lagrange nodes
+        std::vector<Dune::FieldVector<double,dim> > lagrangeNodes(targetLFE.localBasis().size());
 
-      // Generate position of the Lagrange nodes
-      std::vector<Dune::FieldVector<double,dim> > lagrangeNodes(targetLFE.localBasis().size());
+        for (int i=0; i<dim; i++) {
+          CoordinateFunction lFunction(i);
+          std::vector<Dune::FieldVector<double,1> > coordinates;
+          targetLFE.localInterpolation().interpolate(lFunction, coordinates);
 
-      for (int i=0; i<dim; i++) {
-        CoordinateFunction lFunction(i);
-        std::vector<Dune::FieldVector<double,1> > coordinates;
-        targetLFE.localInterpolation().interpolate(lFunction, coordinates);
+          for (size_t j=0; j<coordinates.size(); j++)
+            lagrangeNodes[j][i] = coordinates[j];
 
-        for (size_t j=0; j<coordinates.size(); j++)
-          lagrangeNodes[j][i] = coordinates[j];
+        }
 
-      }
+        for (size_t i=0; i<targetLFE.localCoefficients().size(); i++) {
 
-      for (size_t i=0; i<targetLFE.localCoefficients().size(); i++) {
+          typename GridType::template Codim<0>::EntityPointer element(eIt);
 
-        typename GridType::template Codim<0>::EntityPointer element(eIt);
+          Dune::FieldVector<double,dim> pos = lagrangeNodes[i];
 
-        Dune::FieldVector<double,dim> pos = lagrangeNodes[i];
+          if (handled[targetBasis.index(*eIt,i)][0])
+            continue;
 
-        if (handled[targetBasis.index(*eIt,i)][0])
-          continue;
+          handled[targetBasis.index(*eIt,i)] = true;
 
-        handled[targetBasis.index(*eIt,i)] = true;
+          assert(checkInside(element->type(), pos, 1e-7));
 
-        assert(checkInside(element->type(), pos, 1e-7));
+          // ////////////////////////////////////////////////////////////////////
+          // Get an element on the coarsest grid which contains the vertex and
+          // its local coordinates there
+          // ////////////////////////////////////////////////////////////////////
+          while (element->level() != 0) {
 
-        // ////////////////////////////////////////////////////////////////////
-        // Get an element on the coarsest grid which contains the vertex and
-        // its local coordinates there
-        // ////////////////////////////////////////////////////////////////////
-        while (element->level() != 0) {
+            pos = element->geometryInFather().global(pos);
+            element = element->father();
 
-          pos = element->geometryInFather().global(pos);
-          element = element->father();
+            assert(checkInside(element->type(), pos, 1e-7));
+          }
 
-          assert(checkInside(element->type(), pos, 1e-7));
-        }
+          // ////////////////////////////////////////////////////////////////////
+          //   Find the corresponding coarse grid element on the adaptive grid.
+          //   This is a linear algorithm, but we expect the coarse grid to be small.
+          // ////////////////////////////////////////////////////////////////////
+          Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> uniformP0Mapper (targetGrid,  0);
+          Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> adaptiveP0Mapper(sourceGrid, 0);
+          int coarseIndex = uniformP0Mapper.map(*element);
 
-        // ////////////////////////////////////////////////////////////////////
-        //   Find the corresponding coarse grid element on the adaptive grid.
-        //   This is a linear algorithm, but we expect the coarse grid to be small.
-        // ////////////////////////////////////////////////////////////////////
-        Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> uniformP0Mapper (targetGrid,  0);
-        Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> adaptiveP0Mapper(sourceGrid, 0);
-        int coarseIndex = uniformP0Mapper.map(*element);
+          typename GridType::template Codim<0>::LevelIterator adaptEIt    = sourceGrid.template lbegin<0>(0);
+          typename GridType::template Codim<0>::LevelIterator adaptEEndIt = sourceGrid.template lend<0>(0);
 
-        typename GridType::template Codim<0>::LevelIterator adaptEIt    = sourceGrid.template lbegin<0>(0);
-        typename GridType::template Codim<0>::LevelIterator adaptEEndIt = sourceGrid.template lend<0>(0);
+          for (; adaptEIt!=adaptEEndIt; ++adaptEIt)
+            if (adaptiveP0Mapper.map(*adaptEIt) == coarseIndex)
+              break;
 
-        for (; adaptEIt!=adaptEEndIt; ++adaptEIt)
-          if (adaptiveP0Mapper.map(*adaptEIt) == coarseIndex)
-            break;
+          element = adaptEIt;
 
-        element = adaptEIt;
+          // ////////////////////////////////////////////////////////////////////////
+          //   Find a corresponding point (not necessarily vertex) on the leaf level
+          //   of the adaptive grid.
+          // ////////////////////////////////////////////////////////////////////////
 
-        // ////////////////////////////////////////////////////////////////////////
-        //   Find a corresponding point (not necessarily vertex) on the leaf level
-        //   of the adaptive grid.
-        // ////////////////////////////////////////////////////////////////////////
+          while (!element->isLeaf()) {
 
-        while (!element->isLeaf()) {
 
+            typename GridType::template Codim<0>::Entity::HierarchicIterator hIt    = element->hbegin(element->level()+1);
+            typename GridType::template Codim<0>::Entity::HierarchicIterator hEndIt = element->hend(element->level()+1);
 
-          typename GridType::template Codim<0>::Entity::HierarchicIterator hIt    = element->hbegin(element->level()+1);
-          typename GridType::template Codim<0>::Entity::HierarchicIterator hEndIt = element->hend(element->level()+1);
+            Dune::FieldVector<double,dim> childPos;
+            assert(checkInside(element->type(), pos, 1e-7));
 
-          Dune::FieldVector<double,dim> childPos;
-          assert(checkInside(element->type(), pos, 1e-7));
+            for (; hIt!=hEndIt; ++hIt) {
 
-          for (; hIt!=hEndIt; ++hIt) {
+              childPos = hIt->geometryInFather().local(pos);
+              if (checkInside(hIt->type(), childPos, 1e-7))
+                break;
 
-            childPos = hIt->geometryInFather().local(pos);
-            if (checkInside(hIt->type(), childPos, 1e-7))
-              break;
+            }
+
+            assert(hIt!=hEndIt);
+            element = hIt;
+            pos     = childPos;
 
           }
 
-          assert(hIt!=hEndIt);
-          element = hIt;
-          pos     = childPos;
+          // ////////////////////////////////////////////////////////////////////////
+          //   Sample adaptive function
+          // ////////////////////////////////////////////////////////////////////////
+          sourceFunction.evaluateLocal(*element, pos,
+                                       target[targetBasis.index(*eIt,i)]);
 
         }
 
-        // ////////////////////////////////////////////////////////////////////////
-        //   Sample adaptive function
-        // ////////////////////////////////////////////////////////////////////////
-        sourceFunction.evaluateLocal(*element, pos,
-                                     target[targetBasis.index(*eIt,i)]);
-
       }
 
     }
 
-  }
 
+    template <int blocksize>
+    static double computeNormSquared(const GridType& grid,
+                                     const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& x,
+                                     const LocalOperatorAssembler<GridType,
+                                         typename Basis::LocalFiniteElement,
+                                         typename Basis::LocalFiniteElement,
+                                         Dune::FieldMatrix<double,1,1> >* localStiffness)
+    {
+      // ///////////////////////////////////////////////////////////////////////////////////////////
+      //   Compute the energy norm
+      // ///////////////////////////////////////////////////////////////////////////////////////////
 
-  template <int blocksize>
-  static double computeNormSquared(const GridType& grid,
-                                   const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& x,
-                                   const LocalOperatorAssembler<GridType,
-                                       typename Basis::LocalFiniteElement,
-                                       typename Basis::LocalFiniteElement,
-                                       Dune::FieldMatrix<double,1,1> >* localStiffness)
-  {
-    // ///////////////////////////////////////////////////////////////////////////////////////////
-    //   Compute the energy norm
-    // ///////////////////////////////////////////////////////////////////////////////////////////
+      double energyNormSquared = 0;
+      Basis basis(grid.leafGridView());
+      Dune::Matrix<Dune::FieldMatrix<double,1,1> > localMatrix;
 
-    double energyNormSquared = 0;
-    Basis basis(grid.leafGridView());
-    Dune::Matrix<Dune::FieldMatrix<double,1,1> > localMatrix;
+      typename GridType::template Codim<0>::LeafIterator eIt    = grid.template leafbegin<0>();
+      typename GridType::template Codim<0>::LeafIterator eEndIt = grid.template leafend<0>();
 
-    typename GridType::template Codim<0>::LeafIterator eIt    = grid.template leafbegin<0>();
-    typename GridType::template Codim<0>::LeafIterator eEndIt = grid.template leafend<0>();
+      for (; eIt!=eEndIt; ++eIt) {
 
-    for (; eIt!=eEndIt; ++eIt) {
+        size_t nLocalDofs = basis.getLocalFiniteElement(*eIt).localCoefficients().size();
 
-      size_t nLocalDofs = basis.getLocalFiniteElement(*eIt).localCoefficients().size();
+        localMatrix.setSize(nLocalDofs,nLocalDofs);
 
-      localMatrix.setSize(nLocalDofs,nLocalDofs);
+        localStiffness->assemble(*eIt, localMatrix,
+                                 basis.getLocalFiniteElement(*eIt),
+                                 basis.getLocalFiniteElement(*eIt));
 
-      localStiffness->assemble(*eIt, localMatrix,
-                               basis.getLocalFiniteElement(*eIt),
-                               basis.getLocalFiniteElement(*eIt));
+        for (size_t i=0; i<nLocalDofs; i++)
+          for (size_t j=0; j<nLocalDofs; j++)
+            energyNormSquared += localMatrix[i][j][0][0] * (x[basis.index(*eIt,i)] * x[basis.index(*eIt,j)]);
 
-      for (size_t i=0; i<nLocalDofs; i++)
-        for (size_t j=0; j<nLocalDofs; j++)
-          energyNormSquared += localMatrix[i][j][0][0] * (x[basis.index(*eIt,i)] * x[basis.index(*eIt,j)]);
+      }
 
+      return energyNormSquared;
     }
 
-    return energyNormSquared;
-  }
+    static double compute(const GridType& uniformGrid,
+                          const std::vector<TargetSpace>& uniformSolution,
+                          const GridType& adaptiveGrid,
+                          const std::vector<TargetSpace>& adaptiveSolution,
+                          const LocalOperatorAssembler<GridType,
+                              typename Basis::LocalFiniteElement,
+                              typename Basis::LocalFiniteElement,
+                              Dune::FieldMatrix<double,1,1> >* localStiffness)
+    {
+      std::vector<TargetSpace> uniformAdaptiveSolution;
 
-  static double compute(const GridType& uniformGrid,
-                        const std::vector<TargetSpace>& uniformSolution,
-                        const GridType& adaptiveGrid,
-                        const std::vector<TargetSpace>& adaptiveSolution,
-                        const LocalOperatorAssembler<GridType,
-                            typename Basis::LocalFiniteElement,
-                            typename Basis::LocalFiniteElement,
-                            Dune::FieldMatrix<double,1,1> >* localStiffness)
-  {
-    std::vector<TargetSpace> uniformAdaptiveSolution;
+      interpolate(adaptiveGrid, adaptiveSolution, uniformGrid, uniformAdaptiveSolution);
 
-    interpolate(adaptiveGrid, adaptiveSolution, uniformGrid, uniformAdaptiveSolution);
+      // Compute difference in the embedding space
+      Dune::BlockVector<typename TargetSpace::CoordinateType> difference(uniformSolution.size());
 
-    // Compute difference in the embedding space
-    Dune::BlockVector<typename TargetSpace::CoordinateType> difference(uniformSolution.size());
+      for (size_t i=0; i<difference.size(); i++)
+        difference[i] = uniformAdaptiveSolution[i].globalCoordinates() - uniformSolution[i].globalCoordinates();
 
-    for (size_t i=0; i<difference.size(); i++)
-      difference[i] = uniformAdaptiveSolution[i].globalCoordinates() - uniformSolution[i].globalCoordinates();
+      //std::cout << "new difference:\n" << difference << std::endl;
+      //         std::cout << "new projected:\n" << std::endl;
+      //         for (size_t i=0; i<difference.size(); i++)
+      //             std::cout << uniformAdaptiveSolution[i].globalCoordinates() << std::endl;
 
-    //std::cout << "new difference:\n" << difference << std::endl;
-    //         std::cout << "new projected:\n" << std::endl;
-    //         for (size_t i=0; i<difference.size(); i++)
-    //             std::cout << uniformAdaptiveSolution[i].globalCoordinates() << std::endl;
+      return computeNormSquared(uniformGrid, difference, localStiffness);
+    }
+  };
 
-    return computeNormSquared(uniformGrid, difference, localStiffness);
-  }
-};
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/gramschmidtsolver.hh

@@ -5,130 +5,136 @@
 
 #include <dune/gfe/symmetricmatrix.hh>
 
-/** \brief Direct solver for a dense symmetric linear system, using an orthonormal basis
- *
- * This solver computes an A-orthonormal basis, and uses that to compute the solution
- * of a linear system.  The advantage of this is that it works even if the matrix is
- * known to have a non-trivial kernel.  This happens in GFE applications when the Hessian
- * of a functional on a manifold is given in coordinates of the surrounding space.
- * The method is efficient for the small systems that we consider in GFE applications.
- *
- * \tparam field_type Type used to store scalars
- * \tparam embeddedDim Number of rows and columns of the linear system
- * \tparam rank The rank of the matrix
- */
-template <class field_type, int rank, int embeddedDim>
-class GramSchmidtSolver
-{
-  /** \brief Normalize a vector to unit length, measured in a matrix norm
-   * \param matrix The matrix inducing the matrix norm
-   * \param[in,out] v The vector to normalize
-   */
-  static void normalize(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
-                        Dune::FieldVector<field_type,embeddedDim>& v)
-  {
-    using std::sqrt;
-    v /= sqrt(matrix.energyScalarProduct(v,v));
-  }
 
+namespace Dune::GFE
+{
 
-  /** \brief Project vj on vi, and subtract the result from vj
+  /** \brief Direct solver for a dense symmetric linear system, using an orthonormal basis
+   *
+   * This solver computes an A-orthonormal basis, and uses that to compute the solution
+   * of a linear system.  The advantage of this is that it works even if the matrix is
+   * known to have a non-trivial kernel.  This happens in GFE applications when the Hessian
+   * of a functional on a manifold is given in coordinates of the surrounding space.
+   * The method is efficient for the small systems that we consider in GFE applications.
    *
-   * \param matrix The matrix the defines the scalar product
+   * \tparam field_type Type used to store scalars
+   * \tparam embeddedDim Number of rows and columns of the linear system
+   * \tparam rank The rank of the matrix
    */
-  static void project(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
-                      const Dune::FieldVector<field_type,embeddedDim>& vi,
-                      Dune::FieldVector<field_type,embeddedDim>& vj)
+  template <class field_type, int rank, int embeddedDim>
+  class GramSchmidtSolver
   {
+    /** \brief Normalize a vector to unit length, measured in a matrix norm
+     * \param matrix The matrix inducing the matrix norm
+     * \param[in,out] v The vector to normalize
+     */
+    static void normalize(const SymmetricMatrix<field_type,embeddedDim>& matrix,
+                          Dune::FieldVector<field_type,embeddedDim>& v)
+    {
+      using std::sqrt;
+      v /= sqrt(matrix.energyScalarProduct(v,v));
+    }
 
-    field_type energyScalarProduct = matrix.energyScalarProduct(vi,vj);
 
-    for (size_t i=0; i<vj.size(); i++)
-      vj[i] -= energyScalarProduct * vi[i];
+    /** \brief Project vj on vi, and subtract the result from vj
+     *
+     * \param matrix The matrix the defines the scalar product
+     */
+    static void project(const SymmetricMatrix<field_type,embeddedDim>& matrix,
+                        const Dune::FieldVector<field_type,embeddedDim>& vi,
+                        Dune::FieldVector<field_type,embeddedDim>& vj)
+    {
 
-  }
+      field_type energyScalarProduct = matrix.energyScalarProduct(vi,vj);
 
-public:
+      for (size_t i=0; i<vj.size(); i++)
+        vj[i] -= energyScalarProduct * vi[i];
 
-  /** \brief Constructor computing the A-orthogonal basis
-   *
-   * This constructor uses Gram-Schmidt orthogonalization to compute an A-orthogonal basis.
-   * All (non-static) calls to 'solve' will use that basis.  Since computing the basis
-   * is the expensive part, the calls to 'solve' will be comparatively cheap.
-   *
-   * \param basis Any basis of the orthogonal complement of the kernel,
-   *              used as the input for the Gram-Schmidt orthogonalization process
-   */
-  GramSchmidtSolver(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
-                    const Dune::FieldMatrix<field_type,rank,embeddedDim>& basis)
-    : orthonormalBasis_(basis)
-  {
-    // Use the Gram-Schmidt algorithm to compute a basis that is orthonormal
-    // with respect to the given matrix.
-    normalize(matrix, orthonormalBasis_[0]);
+    }
 
-    for (int i=1; i<rank; i++) {
+  public:
+
+    /** \brief Constructor computing the A-orthogonal basis
+     *
+     * This constructor uses Gram-Schmidt orthogonalization to compute an A-orthogonal basis.
+     * All (non-static) calls to 'solve' will use that basis.  Since computing the basis
+     * is the expensive part, the calls to 'solve' will be comparatively cheap.
+     *
+     * \param basis Any basis of the orthogonal complement of the kernel,
+     *              used as the input for the Gram-Schmidt orthogonalization process
+     */
+    GramSchmidtSolver(const SymmetricMatrix<field_type,embeddedDim>& matrix,
+                      const Dune::FieldMatrix<field_type,rank,embeddedDim>& basis)
+      : orthonormalBasis_(basis)
+    {
+      // Use the Gram-Schmidt algorithm to compute a basis that is orthonormal
+      // with respect to the given matrix.
+      normalize(matrix, orthonormalBasis_[0]);
+
+      for (int i=1; i<rank; i++) {
+
+        for (int j=0; j<i; j++)
+          project(matrix, orthonormalBasis_[j], orthonormalBasis_[i]);
+
+        normalize(matrix, orthonormalBasis_[i]);
+      }
 
-      for (int j=0; j<i; j++)
-        project(matrix, orthonormalBasis_[j], orthonormalBasis_[i]);
+    }
 
-      normalize(matrix, orthonormalBasis_[i]);
+    void solve(Dune::FieldVector<field_type,embeddedDim>& x,
+               const Dune::FieldVector<field_type,embeddedDim>& rhs) const
+    {
+      // Solve the system in the orthonormal basis
+      Dune::FieldVector<field_type,rank> orthoCoefficient;
+      for (int i=0; i<rank; i++)
+        orthoCoefficient[i] = rhs*orthonormalBasis_[i];
+
+      // Solution in canonical basis
+      x = 0;
+      for (int i=0; i<rank; i++)
+        x.axpy(orthoCoefficient[i], orthonormalBasis_[i]);
     }
 
-  }
+    /** Solve linear system by constructing an energy-orthonormal basis
 
-  void solve(Dune::FieldVector<field_type,embeddedDim>& x,
-             const Dune::FieldVector<field_type,embeddedDim>& rhs) const
-  {
-    // Solve the system in the orthonormal basis
-    Dune::FieldVector<field_type,rank> orthoCoefficient;
-    for (int i=0; i<rank; i++)
-      orthoCoefficient[i] = rhs*orthonormalBasis_[i];
+     * \param basis Any basis of the space, used as the input for the Gram-Schmidt orthogonalization process
+     */
+    static void solve(const SymmetricMatrix<field_type,embeddedDim>& matrix,
+                      Dune::FieldVector<field_type,embeddedDim>& x,
+                      const Dune::FieldVector<field_type,embeddedDim>& rhs,
+                      const Dune::FieldMatrix<field_type,rank,embeddedDim>& basis)
+    {
+      // Use the Gram-Schmidt algorithm to compute a basis that is orthonormal
+      // with respect to the given matrix.
+      Dune::FieldMatrix<field_type,rank,embeddedDim> orthonormalBasis(basis);
+      normalize(matrix, orthonormalBasis[0]);
 
-    // Solution in canonical basis
-    x = 0;
-    for (int i=0; i<rank; i++)
-      x.axpy(orthoCoefficient[i], orthonormalBasis_[i]);
-  }
+      for (int i=1; i<rank; i++) {
 
-  /** Solve linear system by constructing an energy-orthonormal basis
+        for (int j=0; j<i; j++)
+          project(matrix, orthonormalBasis[j], orthonormalBasis[i]);
 
-   * \param basis Any basis of the space, used as the input for the Gram-Schmidt orthogonalization process
-   */
-  static void solve(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
-                    Dune::FieldVector<field_type,embeddedDim>& x,
-                    const Dune::FieldVector<field_type,embeddedDim>& rhs,
-                    const Dune::FieldMatrix<field_type,rank,embeddedDim>& basis)
-  {
-    // Use the Gram-Schmidt algorithm to compute a basis that is orthonormal
-    // with respect to the given matrix.
-    Dune::FieldMatrix<field_type,rank,embeddedDim> orthonormalBasis(basis);
-    normalize(matrix, orthonormalBasis[0]);
+        normalize(matrix, orthonormalBasis[i]);
+      }
 
-    for (int i=1; i<rank; i++) {
+      // Solve the system in the orthonormal basis
+      Dune::FieldVector<field_type,rank> orthoCoefficient;
+      for (int i=0; i<rank; i++)
+        orthoCoefficient[i] = rhs*orthonormalBasis[i];
 
-      for (int j=0; j<i; j++)
-        project(matrix, orthonormalBasis[j], orthonormalBasis[i]);
+      // Solution in canonical basis
+      x = 0;
+      for (int i=0; i<rank; i++)
+        x.axpy(orthoCoefficient[i], orthonormalBasis[i]);
 
-      normalize(matrix, orthonormalBasis[i]);
     }
 
-    // Solve the system in the orthonormal basis
-    Dune::FieldVector<field_type,rank> orthoCoefficient;
-    for (int i=0; i<rank; i++)
-      orthoCoefficient[i] = rhs*orthonormalBasis[i];
-
-    // Solution in canonical basis
-    x = 0;
-    for (int i=0; i<rank; i++)
-      x.axpy(orthoCoefficient[i], orthonormalBasis[i]);
-
-  }
+  private:
 
-private:
+    Dune::FieldMatrix<field_type,rank,embeddedDim> orthonormalBasis_;
 
-  Dune::FieldMatrix<field_type,rank,embeddedDim> orthonormalBasis_;
+  };
 
-};
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/linearalgebra.hh

@@ -2,234 +2,235 @@
 #define DUNE_GFE_LINEARALGEBRA_HH
 
 #include <random>
+#include <type_traits>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/scaledidmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>  // For PromotedType
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////
 //  Various matrix methods
 ///////////////////////////////////////////////////////////////////////////////////////////
 
-namespace Dune {
-
-  namespace GFE {
+namespace Dune::GFE
+{
 
 #if ADOLC_ADOUBLE_H
-    /** \brief Calculates ret = s*A, where A has as field_type of adouble.
-     *
-     * The function template is disabled if s isn't a scalar or adolc type.
-     */
-    template<typename T1,int m, int n,class = typename std::enable_if< std::is_scalar_v<T1> || std::is_base_of_v<badouble,T1> >::type >
-    auto operator* ( const T1& s, const Dune::FieldMatrix<adouble, m, n> &A)
-    {
-      typedef typename Dune::FieldMatrix<adouble,m,n> :: size_type size_type;
-      Dune::FieldMatrix<adouble,m,n> ret;
-
-      for( size_type i = 0; i < m; ++i )
-        for( size_type j = 0; j < n; ++j )
-          ret[i][j] = s * A[i][j];
-
-      return ret;
-    }
-
-    /** \brief Calculates ret = A*v, where A has as field_type of adouble.
-     *
-     * The function template is disabled if the field_type of v is no an adolc type
-     */
-    template<typename T1,int m, int n,class = typename std::enable_if<  std::is_base_of_v<badouble,T1> >::type >
-    auto operator* (const Dune::FieldMatrix<double, m, n> &A,  const Dune::FieldVector<T1,n>& v )
-    {
-      typedef typename Dune::FieldMatrix<adouble,m,n> :: size_type size_type;
-      Dune::FieldVector<adouble,m> ret(0.0);
-
-      for( size_type i = 0; i < m; ++i )
-        for( size_type j = 0; j < n; ++j )
-          ret[i] +=  A[i][j]*v[j];
-
-      return ret;
-    }
-
-
-    /** \brief Calculates ret = A*s, where A has as field_type of adouble.
-     *
-     * The function template is disabled if s isn't a scalar or adolc type.
-     */
-    template<typename T1,int m, int n,class = typename std::enable_if< std::is_scalar_v<T1> || std::is_base_of_v<badouble,T1> >::type >
-    auto operator* (const Dune::FieldMatrix<adouble, m, n> &A, const T1& s )
-    {
-      return s*A;
-    }
+  /** \brief Calculates ret = s*A, where A has as field_type of adouble.
+   *
+   * The function template is disabled if s isn't a scalar or adolc type.
+   */
+  template<typename T1,int m, int n,class = typename std::enable_if< std::is_scalar_v<T1> || std::is_base_of_v<badouble,T1> >::type >
+  auto operator* ( const T1& s, const Dune::FieldMatrix<adouble, m, n> &A)
+  {
+    typedef typename Dune::FieldMatrix<adouble,m,n> :: size_type size_type;
+    Dune::FieldMatrix<adouble,m,n> ret;
+
+    for( size_type i = 0; i < m; ++i )
+      for( size_type j = 0; j < n; ++j )
+        ret[i][j] = s * A[i][j];
+
+    return ret;
+  }
+
+  /** \brief Calculates ret = A*v, where A has as field_type of adouble.
+   *
+   * The function template is disabled if the field_type of v is no an adolc type
+   */
+  template<typename T1,int m, int n,class = typename std::enable_if<  std::is_base_of_v<badouble,T1> >::type >
+  auto operator* (const Dune::FieldMatrix<double, m, n> &A,  const Dune::FieldVector<T1,n>& v )
+  {
+    typedef typename Dune::FieldMatrix<adouble,m,n> :: size_type size_type;
+    Dune::FieldVector<adouble,m> ret(0.0);
+
+    for( size_type i = 0; i < m; ++i )
+      for( size_type j = 0; j < n; ++j )
+        ret[i] +=  A[i][j]*v[j];
+
+    return ret;
+  }
+
+
+  /** \brief Calculates ret = A*s, where A has as field_type of adouble.
+   *
+   * The function template is disabled if s isn't a scalar or adolc type.
+   */
+  template<typename T1,int m, int n,class = typename std::enable_if< std::is_scalar_v<T1> || std::is_base_of_v<badouble,T1> >::type >
+  auto operator* (const Dune::FieldMatrix<adouble, m, n> &A, const T1& s )
+  {
+    return s*A;
+  }
 #endif
 
-    /** \brief Return the trace of a matrix */
-    template <class T, int n>
-    static T trace(const FieldMatrix<T,n,n>& A)
-    {
-      T trace = 0;
-      for (int i=0; i<n; i++)
-        trace += A[i][i];
-      return trace;
-    }
-
-    /** \brief Return the square of the trace of a matrix */
-    template <class T, int n>
-    static T traceSquared(const FieldMatrix<T,n,n>& A)
-    {
-      T trace = 0;
-      for (int i=0; i<n; i++)
-        trace += A[i][i];
-      return trace*trace;
-    }
-
-    /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
-    template <class T, int n>
-    static FieldMatrix<T,n,n> sym(const FieldMatrix<T,n,n>& A)
-    {
-      FieldMatrix<T,n,n> result;
-      for (int i=0; i<n; i++)
-        for (int j=0; j<n; j++)
-          result[i][j] = 0.5 * (A[i][j] + A[j][i]);
-      return result;
-    }
-
-    /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
-    template <class T, int n>
-    static FieldMatrix<T,n,n> skew(const FieldMatrix<T,n,n>& A)
-    {
-      FieldMatrix<T,n,n> result;
-      for (int i=0; i<n; i++)
-        for (int j=0; j<n; j++)
-          result[i][j] = 0.5 * (A[i][j] - A[j][i]);
-      return result;
-    }
-
-    /** \brief Compute the deviator of a matrix A */
-    template <class T, int n>
-    static FieldMatrix<T,n,n> dev(const FieldMatrix<T,n,n>& A)
-    {
-      FieldMatrix<T,n,n> result = A;
-      auto t = trace(A);
-      for (int i=0; i<n; i++)
-        result[i][i] -= t / n;
-      return result;
-    }
-
-    /** \brief Return the transposed matrix */
-    template <class T, int n, int m>
-    static FieldMatrix<T,m,n> transpose(const FieldMatrix<T,n,m>& A)
-    {
-      FieldMatrix<T,m,n> result;
-
-      for (int i=0; i<m; i++)
-        for (int j=0; j<n; j++)
-          result[i][j] = A[j][i];
-
-      return result;
-    }
-
-    /** \brief The Frobenius (i.e., componentwise) product of two matrices */
-    template <class T, int n>
-    static T frobeniusProduct(const FieldMatrix<T,n,n>& A, const FieldMatrix<T,n,n>& B)
-    {
-      T result(0.0);
-
-      for (int i=0; i<n; i++)
-        for (int j=0; j<n; j++)
-          result += A[i][j] * B[i][j];
-
-      return result;
-    }
-
-
-    /** \brief Return a*b^T */
-    template <class T1,class T2, int n, int m>
-    static auto dyadicProduct(const FieldVector<T1,n>& a, const FieldVector<T2,m>& b)
-    {
-      using ScalarResultType = typename Dune::PromotionTraits<T1,T2>::PromotedType;
-      FieldMatrix<ScalarResultType,n,m> result;
-      for (int i=0; i<n; i++)
-        for (int j=0; j<m; j++)
-          result[i][j] = a[i]*b[j];
-
-      return result;
-    }
-
-
-    /** \brief Get the requested column of fieldmatrix */
-    template<typename field_type, int cols, int rows>
-    auto col(const Dune::FieldMatrix<field_type, rows, cols> &mat, const int requestedCol)
-    {
-      Dune::FieldVector<field_type, rows> col;
-
-      for (int i = 0; i < rows; ++i)
-        col[i] = mat[i][requestedCol];
-
-      return col;
-    }
-
-    /** \brief Return a segment of a FieldVector from lower up to lower+size-1 */
-    template< int lower, int size,typename field_type,int n>
-    static FieldVector<field_type,size> segment(const FieldVector<field_type,n>& v)
-    {
-      FieldVector<field_type,size> res;
-      std::copy(v.begin()+lower,v.begin()+lower+size,res.begin());
-      return res;
-    }
-
-    /** \brief Return a segment of a FieldVector from lower up to lower+size-1
-     * lower is unknown at compile time*/
-    template< int size,typename field_type,int n>
-    static FieldVector<field_type,size> segmentAt(const FieldVector<field_type,n>& v,const size_t lower)
-    {
-      FieldVector<field_type,size> res;
-      std::copy(v.begin()+lower,v.begin()+lower+size,res.begin());
-      return res;
-    }
-
-    /** \brief Return a block of a FieldMatrix  (lower1...lower1+size1-1,lower2...lower2+size2-1 */
-    template< int lower1, int size1, int lower2,int size2,typename field_type,int n,int m>
-    static auto block(const FieldMatrix<field_type,n,m>& v)
-    {
-      static_assert(lower1+size1<=n && lower2+size2<=m, "Size mismatch for Block!");
-      FieldMatrix<field_type,size1,size2> res;
-
-      for(int i=lower1; i<lower1+size1; ++i)
-        for(int j=lower2; j<lower2+size2; ++j)
-          res[i-lower1][j-lower2] = v[i][j];
-      return res;
-    }
-
-    /** \brief Return a block of a FieldMatrix  (lower1...lower1+size1-1,lower2...lower2+size2-1
-     * * lower1 and lower2 are unknown at compile time*/
-    template< int size1,int size2,typename field_type,int n,int m>
-    static auto blockAt(const FieldMatrix<field_type,n,m>& v, const size_t& lower1, const size_t& lower2)
-    {
-      assert(lower1+size1<=n && lower2+size2<=m);
-      FieldMatrix<field_type,size1,size2> res;
-
-      for(size_t i=lower1; i<lower1+size1; ++i)
-        for(size_t j=lower2; j<lower2+size2; ++j)
-          res[i-lower1][j-lower2] = v[i][j];
-      return res;
-    }
-
-    /** \brief Generates FieldVector with random entries in the range -1..1 */
-    template<typename field_type,int n>
-    auto randomFieldVector(field_type lower=-1, field_type upper=1)
-    {
-      std::random_device rd;
-      std::mt19937 mt(rd());
-      std::uniform_real_distribution<field_type> dist(lower, upper);
-      auto rand = [&dist,&mt](){
-                    return dist(mt);
-                  };
-      FieldVector<field_type,n> vec;
-      std::generate(vec.begin(), vec.end(), rand);
-      return vec;
-    }
-  }
-}
+  /** \brief Return the trace of a matrix */
+  template <class T, int n>
+  static T trace(const FieldMatrix<T,n,n>& A)
+  {
+    T trace = 0;
+    for (int i=0; i<n; i++)
+      trace += A[i][i];
+    return trace;
+  }
+
+  /** \brief Return the square of the trace of a matrix */
+  template <class T, int n>
+  static T traceSquared(const FieldMatrix<T,n,n>& A)
+  {
+    T trace = 0;
+    for (int i=0; i<n; i++)
+      trace += A[i][i];
+    return trace*trace;
+  }
+
+  /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
+  template <class T, int n>
+  static FieldMatrix<T,n,n> sym(const FieldMatrix<T,n,n>& A)
+  {
+    FieldMatrix<T,n,n> result;
+    for (int i=0; i<n; i++)
+      for (int j=0; j<n; j++)
+        result[i][j] = 0.5 * (A[i][j] + A[j][i]);
+    return result;
+  }
+
+  /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
+  template <class T, int n>
+  static FieldMatrix<T,n,n> skew(const FieldMatrix<T,n,n>& A)
+  {
+    FieldMatrix<T,n,n> result;
+    for (int i=0; i<n; i++)
+      for (int j=0; j<n; j++)
+        result[i][j] = 0.5 * (A[i][j] - A[j][i]);
+    return result;
+  }
+
+  /** \brief Compute the deviator of a matrix A */
+  template <class T, int n>
+  static FieldMatrix<T,n,n> dev(const FieldMatrix<T,n,n>& A)
+  {
+    FieldMatrix<T,n,n> result = A;
+    auto t = trace(A);
+    for (int i=0; i<n; i++)
+      result[i][i] -= t / n;
+    return result;
+  }
+
+  /** \brief Return the transposed matrix */
+  template <class T, int n, int m>
+  static FieldMatrix<T,m,n> transpose(const FieldMatrix<T,n,m>& A)
+  {
+    FieldMatrix<T,m,n> result;
+
+    for (int i=0; i<m; i++)
+      for (int j=0; j<n; j++)
+        result[i][j] = A[j][i];
+
+    return result;
+  }
+
+  /** \brief The Frobenius (i.e., componentwise) product of two matrices */
+  template <class T, int n>
+  static T frobeniusProduct(const FieldMatrix<T,n,n>& A, const FieldMatrix<T,n,n>& B)
+  {
+    T result(0.0);
+
+    for (int i=0; i<n; i++)
+      for (int j=0; j<n; j++)
+        result += A[i][j] * B[i][j];
+
+    return result;
+  }
+
+
+  /** \brief Return a*b^T */
+  template <class T1,class T2, int n, int m>
+  static auto dyadicProduct(const FieldVector<T1,n>& a, const FieldVector<T2,m>& b)
+  {
+    using ScalarResultType = typename Dune::PromotionTraits<T1,T2>::PromotedType;
+    FieldMatrix<ScalarResultType,n,m> result;
+    for (int i=0; i<n; i++)
+      for (int j=0; j<m; j++)
+        result[i][j] = a[i]*b[j];
+
+    return result;
+  }
+
+
+  /** \brief Get the requested column of fieldmatrix */
+  template<typename field_type, int cols, int rows>
+  auto col(const Dune::FieldMatrix<field_type, rows, cols> &mat, const int requestedCol)
+  {
+    Dune::FieldVector<field_type, rows> col;
+
+    for (int i = 0; i < rows; ++i)
+      col[i] = mat[i][requestedCol];
+
+    return col;
+  }
+
+  /** \brief Return a segment of a FieldVector from lower up to lower+size-1 */
+  template< int lower, int size,typename field_type,int n>
+  static FieldVector<field_type,size> segment(const FieldVector<field_type,n>& v)
+  {
+    FieldVector<field_type,size> res;
+    std::copy(v.begin()+lower,v.begin()+lower+size,res.begin());
+    return res;
+  }
+
+  /** \brief Return a segment of a FieldVector from lower up to lower+size-1
+   * lower is unknown at compile time*/
+  template< int size,typename field_type,int n>
+  static FieldVector<field_type,size> segmentAt(const FieldVector<field_type,n>& v,const size_t lower)
+  {
+    FieldVector<field_type,size> res;
+    std::copy(v.begin()+lower,v.begin()+lower+size,res.begin());
+    return res;
+  }
+
+  /** \brief Return a block of a FieldMatrix  (lower1...lower1+size1-1,lower2...lower2+size2-1 */
+  template< int lower1, int size1, int lower2,int size2,typename field_type,int n,int m>
+  static auto block(const FieldMatrix<field_type,n,m>& v)
+  {
+    static_assert(lower1+size1<=n && lower2+size2<=m, "Size mismatch for Block!");
+    FieldMatrix<field_type,size1,size2> res;
+
+    for(int i=lower1; i<lower1+size1; ++i)
+      for(int j=lower2; j<lower2+size2; ++j)
+        res[i-lower1][j-lower2] = v[i][j];
+    return res;
+  }
+
+  /** \brief Return a block of a FieldMatrix  (lower1...lower1+size1-1,lower2...lower2+size2-1
+   * * lower1 and lower2 are unknown at compile time*/
+  template< int size1,int size2,typename field_type,int n,int m>
+  static auto blockAt(const FieldMatrix<field_type,n,m>& v, const size_t& lower1, const size_t& lower2)
+  {
+    assert(lower1+size1<=n && lower2+size2<=m);
+    FieldMatrix<field_type,size1,size2> res;
+
+    for(size_t i=lower1; i<lower1+size1; ++i)
+      for(size_t j=lower2; j<lower2+size2; ++j)
+        res[i-lower1][j-lower2] = v[i][j];
+    return res;
+  }
+
+  /** \brief Generates FieldVector with random entries in the range -1..1 */
+  template<typename field_type,int n>
+  auto randomFieldVector(field_type lower=-1, field_type upper=1)
+  {
+    std::random_device rd;
+    std::mt19937 mt(rd());
+    std::uniform_real_distribution<field_type> dist(lower, upper);
+    auto rand = [&dist,&mt](){
+                  return dist(mt);
+                };
+    FieldVector<field_type,n> vec;
+    std::generate(vec.begin(), vec.end(), rand);
+    return vec;
+  }
 
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/localgeodesicfefunction.hh

-#ifndef LOCAL_GEODESIC_FE_FUNCTION_HH
-#define LOCAL_GEODESIC_FE_FUNCTION_HH
-
-#include <vector>
-
-#include <dune/common/fvector.hh>
-
-#include <dune/geometry/type.hh>
-
-#include <dune/gfe/averagedistanceassembler.hh>
-#include <dune/gfe/targetspacertrsolver.hh>
-#include <dune/gfe/localquickanddirtyfefunction.hh>
-#include <dune/gfe/spaces/productmanifold.hh>
-#include <dune/gfe/spaces/realtuple.hh>
-#include <dune/gfe/spaces/rotation.hh>
-
-#include <dune/gfe/tensor3.hh>
-#include <dune/gfe/tensorssd.hh>
-#include <dune/gfe/linearalgebra.hh>
-
-// forward declaration
-template <class LocalFiniteElement, class TargetSpace>
-class LocalGfeTestFunctionBasis;
-
-/** \brief A function defined by simplicial geodesic interpolation
-           from the reference element to a Riemannian manifold.
-
-   \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 TS TargetSpace: The manifold that the function takes its values in
- */
-template <int dim, class ctype, class LocalFiniteElement, class TS>
-class LocalGeodesicFEFunction
-{
-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;
-
-  friend class LocalGfeTestFunctionBasis<LocalFiniteElement,TargetSpace>;
-public:
-
-  /** \brief The type used for derivatives */
-  typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
-
-  /** \brief The type used for derivatives of the gradient with respect to coefficients */
-  typedef Tensor3<RT,embeddedDim,embeddedDim,dim> DerivativeOfGradientWRTCoefficientType;
-
-  /** \brief Constructor
-   * \param localFiniteElement A Lagrangian finite element that provides the interpolation points
-   * \param coefficients Values of the function at the Lagrange points
-   */
-  LocalGeodesicFEFunction(const LocalFiniteElement& localFiniteElement,
-                          const std::vector<TargetSpace>& coefficients)
-    : localFiniteElement_(localFiniteElement),
-    coefficients_(coefficients)
-  {
-    assert(localFiniteElement_.localBasis().size() == coefficients_.size());
-  }
-
-  /** \brief Rebind the FEFunction to another TargetSpace */
-  template<class U>
-  struct rebind
-  {
-    using other = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,U>;
-  };
-
-  /** \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 The scalar finite element used as the interpolation weights
-   *
-   * \note This method was added for InterpolationDerivatives, which needs it
-   * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
-   * number type.  This is not optimal, because the localFiniteElement
-   * really is an implementation detail of LocalGeodesicFEFunction and
-   * should not be needed just to copy an entire object.  Other non-Euclidean
-   * interpolation rules may not have such a finite element at all.
-   * Therefore, this method may disappear again eventually.
-   */
-  const LocalFiniteElement& localFiniteElement() const
-  {
-    return localFiniteElement_;
-  }
-
-  /** \brief Evaluate the function */
-  TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const;
-
-  /** \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;
-
-  /** \brief Evaluate the value and the derivative of the interpolation function
-   */
-  std::pair<TargetSpace,DerivativeType> evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const;
-
-  /** \brief Evaluate the derivative of the function value with respect to a coefficient */
-  void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                               int coefficient,
-                                               Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& derivative) const;
-
-  /** \brief Evaluate the derivative of the function value with respect to a coefficient */
-  void evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                                 int coefficient,
-                                                 Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& derivative) const;
-
-  /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
-  void evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                                  int coefficient,
-                                                  DerivativeOfGradientWRTCoefficientType& result) const;
-
-  /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
-  void evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                                    int coefficient,
-                                                    DerivativeOfGradientWRTCoefficientType& result) const;
-
-  /** \brief Get the i'th base coefficient. */
-  const TargetSpace& coefficient(int i) const
-  {
-    return coefficients_[i];
-  }
-private:
-
-  static Dune::SymmetricMatrix<RT,embeddedDim> pseudoInverse(const Dune::SymmetricMatrix<RT,embeddedDim>& dFdq,
-                                                             const TargetSpace& q)
-  {
-    const int shortDim = TargetSpace::TangentVector::dimension;
-
-    // the orthonormal frame
-    Dune::FieldMatrix<RT,shortDim,embeddedDim> O = q.orthonormalFrame();
-
-    // compute A = O dFDq O^T
-    // TODO A really is symmetric, too
-    Dune::FieldMatrix<RT,shortDim,shortDim> A;
-    for (int i=0; i<shortDim; i++)
-      for (int j=0; j<shortDim; j++) {
-        A[i][j] = 0;
-        for (int k=0; k<embeddedDim; k++)
-          for (int l=0; l<embeddedDim; l++)
-            A[i][j] += O[i][k] * ((k>=l) ? dFdq(k,l) : dFdq(l,k)) *O[j][l];
-      }
-
-    A.invert();
-
-    Dune::SymmetricMatrix<RT,embeddedDim> result;
-    result = 0.0;
-    for (int i=0; i<embeddedDim; i++)
-      for (int j=0; j<=i; j++)
-        for (int k=0; k<shortDim; k++)
-          for (int l=0; l<shortDim; l++)
-            result(i,j) += O[k][i]*A[k][l]*O[l][j];
-
-    return result;
-  }
-
-  /** \brief Compute derivate of F(w,q) (the derivative of the weighted distance fctl) wrt to w */
-  Dune::Matrix<RT> computeDFdw(const TargetSpace& q) const
-  {
-    Dune::Matrix<RT> dFdw(embeddedDim,localFiniteElement_.localBasis().size());
-    for (size_t i=0; i<localFiniteElement_.localBasis().size(); i++) {
-      Dune::FieldVector<RT,embeddedDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
-      for (int j=0; j<embeddedDim; j++)
-        dFdw[j][i] = tmp[j];
-    }
-    return dFdw;
-  }
-
-  Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> computeDqDqF(const std::vector<Dune::FieldVector<ctype,1> >& w, const TargetSpace& q) const
-  {
-    Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> result;
-    result = Tensor3<RT,embeddedDim,embeddedDim,embeddedDim>(RT(0));
-    for (size_t i=0; i<w.size(); i++)
-      result.axpy(w[i][0], TargetSpace::thirdDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i],q));
-    return result;
-  }
-
-  /** \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 LocalGeodesicFEFunction<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);
-
-  // The energy functional whose mimimizer is the value of the geodesic interpolation
-  AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
-
-  // Create a reasonable initial iterate for the iterative solver
-  Dune::GFE::LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> localProjectedFEFunction(localFiniteElement_, coefficients_);
-  TargetSpace initialIterate = localProjectedFEFunction.evaluate(local);
-
-  // Iteratively solve the GFE minimization problem
-  TargetSpaceRiemannianTRSolver<TargetSpace> solver;
-
-  solver.setup(&assembler,
-               initialIterate,
-               1e-14,      // tolerance
-               100         // maxNewtonSteps
-               );
-
-  solver.solve();
-
-  return solver.getSol();
-}
-
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-typename LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
-LocalGeodesicFEFunction<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 LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
-LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace& q) const
-{
-  Dune::FieldMatrix<RT, embeddedDim, dim> result;
-
-#if 0  // this is probably faster than the general implementation, but we leave it out for testing purposes
-  if (dim==1) {
-
-    EmbeddedTangentVector tmp = TargetSpace::interpolateDerivative(coefficients_[0], coefficients_[1], local[0]);
-
-    for (int i=0; i<embeddedDim; i++)
-      result[i][0] = tmp[i];
-
-  }
-#endif
-
-  // ////////////////////////////////////////////////////////////////////////
-  //  The derivative is evaluated using the implicit function theorem.
-  //  Hence we need to solve a small system of linear equations.
-  // ////////////////////////////////////////////////////////////////////////
-
-  // the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
-  std::vector<Dune::FieldMatrix<ctype,1,dim> > B(coefficients_.size());
-  localFiniteElement_.localBasis().evaluateJacobian(local, B);
-
-  // compute negative derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
-  Dune::Matrix<RT> dFdw = computeDFdw(q);
-  dFdw *= -1;
-
-  // multiply the two previous matrices: the result is the right hand side
-  // RHS = dFdw * B;
-  // We need to write this multiplication by hand, because B is actually an (#coefficients) times 1 times dim matrix,
-  // and we need to ignore the middle index.
-  Dune::Matrix<RT> RHS(dFdw.N(), dim);
-
-  for (size_t i=0; i<RHS.N(); i++)
-    for (size_t j=0; j<RHS.M(); j++) {
-      RHS[i][j] = 0;
-      for (size_t k=0; k<dFdw.M(); k++)
-        RHS[i][j] += dFdw[i][k]*B[k][0][j];
-    }
-
-  // the actual system matrix
-  std::vector<Dune::FieldVector<ctype,1> > w;
-  localFiniteElement_.localBasis().evaluateFunction(local, w);
-
-  AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
-
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
-  assembler.assembleEmbeddedHessian(q,dFdq);
-
-  // We want to solve
-  //     dFdq * x = rhs
-  // We use the Gram-Schmidt solver because we know that dFdq is rank-deficient.
-
-  const int shortDim = TargetSpace::TangentVector::dimension;
-
-  // the orthonormal frame
-  Dune::FieldMatrix<RT,shortDim,embeddedDim> basis = q.orthonormalFrame();
-  GramSchmidtSolver<RT,shortDim,embeddedDim> gramSchmidtSolver(dFdq, basis);
-
-  for (int i=0; i<dim; i++) {
-
-    Dune::FieldVector<RT,embeddedDim> rhs;
-    for (int j=0; j<embeddedDim; j++)
-      rhs[j] = RHS[j][i];
-
-    Dune::FieldVector<RT,embeddedDim> x;
-    gramSchmidtSolver.solve(x, rhs);
-
-    for (int j=0; j<embeddedDim; j++)
-      result[j][i] = x[j];
-
-  }
-
-  return result;
-}
-
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-std::pair<TargetSpace,typename LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType>
-LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
-{
-  std::pair<TargetSpace,DerivativeType> result;
-  result.first = evaluate(local);
-  result.second = evaluateDerivative(local,result.first);
-  return result;
-}
-
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                        int coefficient,
-                                        Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& result) const
-{
-  // the function value at the point where we are evaluating the derivative
-  TargetSpace q = evaluate(local);
-
-  // dFdq
-  std::vector<Dune::FieldVector<ctype,1> > w;
-  localFiniteElement_.localBasis().evaluateFunction(local,w);
-
-  AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
-
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
-  assembler.assembleEmbeddedHessian(q,dFdq);
-
-  const int shortDim = TargetSpace::TangentVector::dimension;
-
-  // the orthonormal frame
-  Dune::FieldMatrix<RT,shortDim,embeddedDim> O = q.orthonormalFrame();
-
-  // compute A = O dFDq O^T
-  Dune::FieldMatrix<RT,shortDim,shortDim> A;
-  for (int i=0; i<shortDim; i++)
-    for (int j=0; j<shortDim; j++) {
-      A[i][j] = 0;
-      for (int k=0; k<embeddedDim; k++)
-        for (int l=0; l<embeddedDim; l++)
-          A[i][j] += O[i][k] * ((k>=l) ? dFdq(k,l) : dFdq(l,k)) * O[j][l];
-    }
-
-  //
-  Dune::FieldMatrix<RT,embeddedDim,embeddedDim> rhs = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
-  rhs *= -w[coefficient];
-
-  for (int i=0; i<embeddedDim; i++) {
-
-    Dune::FieldVector<RT,shortDim> shortRhs;
-    O.mv(rhs[i],shortRhs);
-
-    Dune::FieldVector<RT,shortDim> shortX;
-    A.solve(shortX,shortRhs);
-
-    O.mtv(shortX,result[i]);
-
-  }
-
-}
-
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                          int coefficient,
-                                          Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& result) const
-{
-  double eps = 1e-6;
-
-  // the function value at the point where we are evaluating the derivative
-  const Dune::FieldMatrix<RT,spaceDim,embeddedDim> B = coefficients_[coefficient].orthonormalFrame();
-
-  Dune::FieldMatrix<RT,spaceDim,embeddedDim> interimResult;
-
-  std::vector<TargetSpace> cornersPlus  = coefficients_;
-  std::vector<TargetSpace> cornersMinus = coefficients_;
-
-  for (int j=0; j<spaceDim; j++) {
-
-    typename TargetSpace::EmbeddedTangentVector forwardVariation = B[j];
-    forwardVariation *= eps;
-    typename TargetSpace::EmbeddedTangentVector backwardVariation = B[j];
-    backwardVariation *= -eps;
-
-    cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
-    cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
-
-    LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> fPlus(localFiniteElement_,cornersPlus);
-    LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> fMinus(localFiniteElement_,cornersMinus);
-
-    TargetSpace hPlus  = fPlus.evaluate(local);
-    TargetSpace hMinus = fMinus.evaluate(local);
-
-    interimResult[j]  = hPlus.globalCoordinates();
-    interimResult[j] -= hMinus.globalCoordinates();
-    interimResult[j] /= 2*eps;
-
-  }
-
-  for (int i=0; i<embeddedDim; i++)
-    for (int j=0; j<embeddedDim; j++) {
-      result[i][j] = 0;
-      for (int k=0; k<spaceDim; k++)
-        result[i][j] += B[k][i]*interimResult[k][j];
-    }
-
-}
-
-
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                           int coefficient,
-                                           DerivativeOfGradientWRTCoefficientType& result) const
-{
-  // the function value at the point where we are evaluating the derivative
-  TargetSpace q = evaluate(local);
-
-  // the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
-  std::vector<Dune::FieldMatrix<ctype,1,dim> > BNested(coefficients_.size());
-  localFiniteElement_.localBasis().evaluateJacobian(local, BNested);
-  Dune::Matrix<RT> B(coefficients_.size(), dim);
-  for (size_t i=0; i<coefficients_.size(); i++)
-    for (size_t j=0; j<dim; j++)
-      B[i][j] = BNested[i][0][j];
-
-  // the actual system matrix
-  std::vector<Dune::FieldVector<ctype,1> > w;
-  localFiniteElement_.localBasis().evaluateFunction(local,w);
-
-  AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
-
-  /** \todo Use a symmetric matrix here */
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
-  assembler.assembleEmbeddedHessian(q,dFdq);
-
-
-  Dune::FieldMatrix<RT,embeddedDim,embeddedDim> mixedDerivative = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
-  TensorSSD<RT,embeddedDim,embeddedDim> dvDwF(coefficients_.size());
-  dvDwF = 0;
-  for (int i=0; i<embeddedDim; i++)
-    for (int j=0; j<embeddedDim; j++)
-      dvDwF(i, j, coefficient) = mixedDerivative[i][j];
-
-
-  // dFDq is not invertible, if the target space is embedded into a higher-dimensional
-  // Euclidean space.  Therefore we use its pseudo inverse.  I don't think that is the
-  // best way, though.
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq,q);
-
-  //
-  Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> dvDqF
-    =  TargetSpace::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(coefficients_[coefficient], q);
-
-  dvDqF = RT(w[coefficient]) * dvDqF;
-
-  // Put it all together
-  // dvq[i][j] = \partial q_j / \partial v_i
-  Dune::FieldMatrix<RT,embeddedDim,embeddedDim> dvq;
-  evaluateDerivativeOfValueWRTCoefficient(local,coefficient,dvq);
-
-  Dune::FieldMatrix<RT, embeddedDim, dim> derivative = evaluateDerivative(local);
-
-  Tensor3<RT, embeddedDim,embeddedDim,embeddedDim> dqdqF;
-  dqdqF = computeDqDqF(w,q);
-
-  TensorSSD<RT, embeddedDim,embeddedDim> dqdwF(coefficients_.size());
-
-  for (size_t k=0; k<coefficients_.size(); k++) {
-    Dune::SymmetricMatrix<RT,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
-    for (int i=0; i<embeddedDim; i++)
-      for (int j=0; j<=i; j++)
-        dqdwF(i, j, k) = dqdwF(j, i, k) = hesse(i,j);
-
-  }
-
-  TensorSSD<RT, embeddedDim,embeddedDim> dqdwF_times_dvq(coefficients_.size());
-  for (int i=0; i<embeddedDim; i++)
-    for (int j=0; j<embeddedDim; j++)
-      for (size_t k=0; k<coefficients_.size(); k++) {
-        dqdwF_times_dvq(i, j, k) = 0;
-        for (int l=0; l<embeddedDim; l++)
-          dqdwF_times_dvq(i, j, k) += dqdwF(l, j, k) * dvq[i][l];
-      }
-
-  Tensor3<RT,embeddedDim,embeddedDim,dim> foo;
-  foo =  -1 * dvDqF*derivative - (dvq*dqdqF)*derivative;
-  TensorSSD<RT,embeddedDim,embeddedDim> bar(dim);
-  bar = dvDwF * B + dqdwF_times_dvq*B;
-
-  for (int i=0; i<embeddedDim; i++)
-    for (int j=0; j<embeddedDim; j++)
-      for (int k=0; k<dim; k++)
-        foo[i][j][k] -= bar(i, j, k);
-
-  result = RT(0);
-  for (int i=0; i<embeddedDim; i++)
-    for (int j=0; j<embeddedDim; j++)
-      for (int k=0; k<dim; k++)
-        for (int l=0; l<embeddedDim; l++)
-          // TODO Smarter implementation of the product with a symmetric matrix
-          result[i][j][k] += ((j>=l) ? dFdqPseudoInv(j,l) : dFdqPseudoInv(l,j)) * foo[i][l][k];
-
-}
-
-
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                             int coefficient,
-                                             DerivativeOfGradientWRTCoefficientType& result) const
-{
-  double eps = 1e-6;
-
-  // loop over the different partial derivatives
-  for (int j=0; j<embeddedDim; j++) {
-
-    std::vector<TargetSpace> cornersPlus  = coefficients_;
-    std::vector<TargetSpace> cornersMinus = coefficients_;
-    typename TargetSpace::CoordinateType aPlus  = coefficients_[coefficient].globalCoordinates();
-    typename TargetSpace::CoordinateType aMinus = coefficients_[coefficient].globalCoordinates();
-    aPlus[j]  += eps;
-    aMinus[j] -= eps;
-    cornersPlus[coefficient]  = TargetSpace(aPlus);
-    cornersMinus[coefficient] = TargetSpace(aMinus);
-    LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> fPlus(localFiniteElement_,cornersPlus);
-    LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> fMinus(localFiniteElement_,cornersMinus);
-
-    Dune::FieldMatrix<RT,embeddedDim,dim> hPlus  = fPlus.evaluateDerivative(local);
-    Dune::FieldMatrix<RT,embeddedDim,dim> hMinus = fMinus.evaluateDerivative(local);
-
-    result[j]  = hPlus;
-    result[j] -= hMinus;
-    result[j] /= 2*eps;
-
-  }
-
-  for (int j=0; j<embeddedDim; j++) {
-
-    TargetSpace q = evaluate(local);
-    Dune::FieldVector<RT,embeddedDim> foo;
-    for (int l=0; l<dim; l++) {
-
-      for (int k=0; k<embeddedDim; k++)
-        foo[k] = result[j][k][l];
-
-      foo = q.projectOntoTangentSpace(foo);
-
-      for (int k=0; k<embeddedDim; k++)
-        result[j][k][l] = foo[k];
-
-    }
-
-  }
-
-}
-
-
-/** \brief A function defined by geodesic interpolation
-           from the reference element to a ProductManifold<RealTuple,Rotation>.
-
-   This is a specialization for speeding up the code.
-
-   \tparam dim Dimension of the reference element
-   \tparam ctype Type used for coordinates on the reference element
- */
-template <int dim, class ctype, class LocalFiniteElement, class field_type>
-class LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement, Dune::GFE::ProductManifold<RealTuple<field_type,3>, Rotation<field_type,3> > >
-{
-  using TargetSpace = Dune::GFE::ProductManifold<RealTuple<field_type,3>, Rotation<field_type,3> >;
-
-  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<field_type, embeddedDim, dim> DerivativeType;
-
-  /** \brief The type used for derivatives of the gradient with respect to coefficients */
-  typedef Tensor3<field_type,embeddedDim,embeddedDim,dim> DerivativeOfGradientWRTCoefficientType;
-
-  /** \brief Constructor */
-  LocalGeodesicFEFunction(const LocalFiniteElement& localFiniteElement,
-                          const std::vector<TargetSpace>& coefficients)
-    : localFiniteElement_(localFiniteElement),
-    coefficients_(coefficients),
-    translationCoefficients_(coefficients.size())
-  {
-    using namespace Dune::Indices;
-    assert(localFiniteElement.localBasis().size() == coefficients.size());
-
-    for (size_t i=0; i<coefficients.size(); i++)
-      translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
-
-    std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
-    for (size_t i=0; i<coefficients.size(); i++)
-      orientationCoefficients[i] = coefficients[i][_1];
-
-    orientationFEFunction_ = std::unique_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (new LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> >(localFiniteElement,orientationCoefficients));
-
-  }
-
-  /** \brief Rebind the FEFunction to another TargetSpace */
-  template<class U>
-  struct rebind
-  {
-    using other = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,U>;
-  };
-
-  /** \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
-  {
-    using namespace Dune::Indices;
-
-    TargetSpace result;
-
-    // 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);
-
-    result[_0] = Dune::FieldVector<field_type,3>(0.0);
-    for (size_t i=0; i<w.size(); i++)
-      result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);
-
-    result[_1] = orientationFEFunction_->evaluate(local);
-    return result;
-  }
-
-  /** \brief Evaluate the derivative of the function */
-  DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
-  {
-    DerivativeType result(0);
-
-    // get translation part
-    std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
-    localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);
-
-    for (size_t i=0; i<translationCoefficients_.size(); i++)
-      for (int j=0; j<3; j++)
-        result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
-
-    // get orientation part
-    Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_->evaluateDerivative(local);
-    for (int i=0; i<4; i++)
-      for (int j=0; j<dim; j++)
-        result[3+i][j] = qResult[i][j];
-
-    return result;
-  }
-
-  /** \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
-  {
-    using namespace Dune::Indices;
-
-    DerivativeType result(0);
-
-    // get translation part
-    std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
-    localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);
-
-    for (size_t i=0; i<translationCoefficients_.size(); i++)
-      for (int j=0; j<3; j++)
-        result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
-
-    // get orientation part
-    Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_->evaluateDerivative(local,q[_1]);
-    for (int i=0; i<4; i++)
-      for (int j=0; j<dim; j++)
-        result[3+i][j] = qResult[i][j];
-
-    return result;
-  }
-
-  /** \brief Evaluate the value and the derivative of the interpolation function
-   */
-  std::pair<TargetSpace,DerivativeType> evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
-  {
-    std::pair<TargetSpace,DerivativeType> result;
-    result.first = evaluate(local);
-    result.second = evaluateDerivative(local,result.first);
-    return result;
-  }
-
-  /** \brief Evaluate the derivative of the function value with respect to a coefficient */
-  void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                               int coefficient,
-                                               Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& derivative) const
-  {
-    derivative = 0;
-
-    // Translation part
-    std::vector<Dune::FieldVector<ctype,1> > w;
-    localFiniteElement_.localBasis().evaluateFunction(local,w);
-    for (int i=0; i<3; i++)
-      derivative[i][i] = w[coefficient];
-
-    // Rotation part
-    Dune::FieldMatrix<field_type,4,4> qDerivative;
-    orientationFEFunction_->evaluateDerivativeOfValueWRTCoefficient(local,coefficient,qDerivative);
-    for (int i=0; i<4; i++)
-      for (int j=0; j<4; j++)
-        derivative[3+i][3+j] = qDerivative[i][j];
-  }
-
-  /** \brief Evaluate the derivative of the function value with respect to a coefficient */
-  void evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                                 int coefficient,
-                                                 Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& derivative) const
-  {
-    derivative = 0;
-
-    // Translation part
-    std::vector<Dune::FieldVector<ctype,1> > w;
-    localFiniteElement_.localBasis().evaluateFunction(local,w);
-    for (int i=0; i<3; i++)
-      derivative[i][i] = w[coefficient];
-
-    // Rotation part
-    Dune::FieldMatrix<ctype,4,4> qDerivative;
-    orientationFEFunction_->evaluateFDDerivativeOfValueWRTCoefficient(local,coefficient,qDerivative);
-    for (int i=0; i<4; i++)
-      for (int j=0; j<4; j++)
-        derivative[3+i][3+j] = qDerivative[i][j];
-  }
-
-  /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
-  void evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                                  int coefficient,
-                                                  DerivativeOfGradientWRTCoefficientType& derivative) const
-  {
-    derivative = field_type(0);
-
-    // Translation part
-    std::vector<Dune::FieldMatrix<ctype,1,dim> > w;
-    localFiniteElement_.localBasis().evaluateJacobian(local,w);
-    for (int i=0; i<3; i++)
-      derivative[i][i] = w[coefficient][0];
-
-    // Rotation part
-    Tensor3<field_type,4,4,dim> qDerivative;
-    orientationFEFunction_->evaluateDerivativeOfGradientWRTCoefficient(local,coefficient,qDerivative);
-    for (int i=0; i<4; i++)
-      for (int j=0; j<4; j++)
-        for (int k=0; k<dim; k++)
-          derivative[3+i][3+j][k] = qDerivative[i][j][k];
-  }
-
-  /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
-  void evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
-                                                    int coefficient,
-                                                    DerivativeOfGradientWRTCoefficientType& derivative) const
-  {
-    derivative = 0;
-
-    // Translation part
-    std::vector<Dune::FieldMatrix<ctype,1,dim> > w;
-    localFiniteElement_.localBasis().evaluateJacobian(local,w);
-    for (int i=0; i<3; i++)
-      derivative[i][i] = w[coefficient][0];
-
-    // Rotation part
-    Tensor3<ctype,4,4,dim> qDerivative;
-    orientationFEFunction_->evaluateFDDerivativeOfGradientWRTCoefficient(local,coefficient,qDerivative);
-    for (int i=0; i<4; i++)
-      for (int j=0; j<4; j++)
-        for (int k=0; k<dim; k++)
-          derivative[3+i][3+j][k] = qDerivative[i][j][k];
-  }
-
-  const 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_;
-
-  // The coefficients of this interpolation rule
-  std::vector<TargetSpace> coefficients_;
-
-  // The coefficients again.  Yes, we store it twice:  To evaluate the interpolation rule efficiently
-  // we need access to the coefficients for the various factor spaces separately.
-  std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;
-
-  std::unique_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > orientationFEFunction_;
-};
-
-
-#endif