InterpolationRules: Return coefficients by const reference

Previously they were returned by value, but I do not remember any
actual reason for that.  Return by const reference leads to a
measurable speed increase when computing the derivatives of
projection-based FE onto the sphere.  It also fixed some
undefined behavior in that code, which took const-references
of the returned temporaries.

Returning the coefficients by const reference is more difficult
for product manifolds, because previously the interpolation rules
stored the coefficients separately for each factor space.
With this patch, these rules now store the coefficients twice:
once in the old separate format, and once as needed for returning
them by reference.  Interpolation rules are not meant to exist
in large numbers, and therefore I do not think that this
extra space consumption matters.

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_;

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_;