From e2d97ab8f95e1bfa63d88ab9b5ae886a581781ab Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sat, 15 Jun 2024 13:40:53 +0200
Subject: [PATCH] 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  | 15 ++++++++++-----
 dune/gfe/localprojectedfefunction.hh | 14 ++++++++++----
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index cb36d6aa5..073c15f35 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 5c394feec..ba16da308 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_;
-- 
GitLab