diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index cb36d6aa583e3a05e00eab11d88c9ad730938ac5..073c15f35bcb45b30af1ddfa59e5f3012eacd1e0 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -118,7 +118,7 @@ public:
                                                     DerivativeOfGradientWRTCoefficientType& result) const;
 
   /** \brief Get the i'th base coefficient. */
-  TargetSpace coefficient(int i) const
+  const TargetSpace& coefficient(int i) const
   {
     return coefficients_[i];
   }
@@ -597,6 +597,7 @@ public:
   LocalGeodesicFEFunction(const LocalFiniteElement& localFiniteElement,
                           const std::vector<TargetSpace>& coefficients)
     : localFiniteElement_(localFiniteElement),
+    coefficients_(coefficients),
     translationCoefficients_(coefficients.size())
   {
     using namespace Dune::Indices;
@@ -787,10 +788,10 @@ public:
           derivative[3+i][3+j][k] = qDerivative[i][j][k];
   }
 
-  TargetSpace coefficient(int i) const {
-    return TargetSpace(translationCoefficients_[i],orientationFEFunction_->coefficient(i));
-
+  const TargetSpace& coefficient(int i) const {
+    return coefficients_[i];
   }
+
 private:
 
   /** \brief The scalar local finite element, which provides the weighting factors
@@ -798,7 +799,11 @@ private:
    */
   const LocalFiniteElement& localFiniteElement_;
 
-  // Coefficients for the two factors of the product manifold
+  // 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_;
diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index 5c394feec193d6670e8cf46bcf62bb39256389c2..ba16da3087b4607235c2d452b543e930e1700967 100644
--- a/dune/gfe/localprojectedfefunction.hh
+++ b/dune/gfe/localprojectedfefunction.hh
@@ -97,7 +97,7 @@ namespace Dune {
       auto evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const;
 
       /** \brief Get the i'th base coefficient. */
-      TargetSpace coefficient(int i) const
+      const TargetSpace& coefficient(int i) const
       {
         return coefficients_[i];
       }
@@ -438,7 +438,7 @@ namespace Dune {
       }
 
       /** \brief Get the i'th base coefficient. */
-      TargetSpace coefficient(int i) const
+      const TargetSpace& coefficient(int i) const
       {
         return coefficients_[i];
       }
@@ -484,6 +484,7 @@ namespace Dune {
       LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
                                const std::vector<TargetSpace>& coefficients)
         : localFiniteElement_(localFiniteElement),
+        coefficients_(coefficients),
         translationCoefficients_(coefficients.size())
       {
         assert(localFiniteElement.localBasis().size() == coefficients.size());
@@ -578,9 +579,9 @@ namespace Dune {
       }
 
       /** \brief Get the i'th base coefficient. */
-      TargetSpace coefficient(int i) const
+      const TargetSpace& coefficient(int i) const
       {
-        return TargetSpace(translationCoefficients_[i],orientationFunction_->coefficient(i));
+        return coefficients_[i];
       }
     private:
 
@@ -589,6 +590,11 @@ namespace Dune {
        */
       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<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type, 3> > > orientationFunction_;