diff --git a/dune/gfe/assemblers/cosseratrodenergy.hh b/dune/gfe/assemblers/cosseratrodenergy.hh
index 7d6c6e67cb497b9a50df3a348e88f39a15bd5d81..d8dde8cd048f7246a1a0f9900f3228ece1d0d08d 100644
--- a/dune/gfe/assemblers/cosseratrodenergy.hh
+++ b/dune/gfe/assemblers/cosseratrodenergy.hh
@@ -179,7 +179,8 @@ namespace Dune::GFE
          const std::vector<TargetSpace>& localCoefficients) const
   {
     const auto& localFiniteElement = localView.tree().finiteElement();
-    LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
+    LocalInterpolationRule localConfiguration;
+    localConfiguration.bind(localFiniteElement, localCoefficients);
 
     const auto& element = localView.element();
 
@@ -187,7 +188,8 @@ namespace Dune::GFE
 
     std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
     using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::other;
-    InactiveLocalInterpolationRule localReferenceConfiguration(localFiniteElement, localReferenceCoefficients);
+    InactiveLocalInterpolationRule localReferenceConfiguration;
+    localReferenceConfiguration.bind(localFiniteElement, localReferenceCoefficients);
 
     // ///////////////////////////////////////////////////////////////////////////////
     //   The following two loops are a reduced integration scheme.  We integrate
diff --git a/dune/gfe/assemblers/l2distancesquaredenergy.hh b/dune/gfe/assemblers/l2distancesquaredenergy.hh
index dcc8ba194d33a62cb1fd86eb8f21d658d95a96fa..4d625a47446ac4a2766cc3fe77a1b95b45219470 100644
--- a/dune/gfe/assemblers/l2distancesquaredenergy.hh
+++ b/dune/gfe/assemblers/l2distancesquaredenergy.hh
@@ -41,7 +41,8 @@ namespace Dune::GFE
 
       const auto& localFiniteElement = localView.tree().finiteElement();
       typedef LocalGeodesicFEFunction<Basis, TargetSpace> LocalGFEFunctionType;
-      LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+      LocalGFEFunctionType localGeodesicFEFunction;
+      localGeodesicFEFunction.bind(localFiniteElement,localSolution);
 
       const auto element = localView.element();
       auto localOrigin = localFunction(*origin_);
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index 8458ad0e173dd809d9592ddbbc283d0b324b1e3f..f0656a40439eb6b049798b696482053b48c6d9c5 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -96,7 +96,8 @@ namespace Dune::GFE
 #else
         const auto& localFiniteElement = Impl::LocalFiniteElementFactory<Basis>::get(localView);
 #endif
-        LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
+        LocalInterpolationRule localInterpolationRule;
+        localInterpolationRule.bind(localFiniteElement,localConfiguration);
 
         // Bind density to the element
         const auto& element = localView.element();
@@ -179,8 +180,10 @@ namespace Dune::GFE
         using LocalGFEFunctionType0 = typename std::tuple_element<0, LocalInterpolationRule>::type;
         using LocalGFEFunctionType1 = typename std::tuple_element<1, LocalInterpolationRule>::type;
 
-        LocalGFEFunctionType0 localGFEFunction0(localFiniteElement0, coefficients[_0]);
-        LocalGFEFunctionType1 localGFEFunction1(localFiniteElement1, coefficients[_1]);
+        LocalGFEFunctionType0 localGFEFunction0;
+        LocalGFEFunctionType1 localGFEFunction1;
+        localGFEFunction0.bind(localFiniteElement0, coefficients[_0]);
+        localGFEFunction1.bind(localFiniteElement1, coefficients[_1]);
 
         // Bind density to the element
         const auto& element = localView.element();
diff --git a/dune/gfe/assemblers/localintegralstiffness.hh b/dune/gfe/assemblers/localintegralstiffness.hh
index 67e20c6a72ccbb73b5367286349bbe009e8fabdc..dab17a1e17756cdd0c8bb025fc36d8bbff523b82 100644
--- a/dune/gfe/assemblers/localintegralstiffness.hh
+++ b/dune/gfe/assemblers/localintegralstiffness.hh
@@ -68,7 +68,8 @@ namespace Dune::GFE
       if constexpr (not Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
       {
         const auto& localFiniteElement = localView.tree().finiteElement();
-        LocalInterpolationRule localInterpolationRule(localFiniteElement,localCoefficients);
+        LocalInterpolationRule localInterpolationRule;
+        localInterpolationRule.bind(localFiniteElement,localCoefficients);
 
         // Bind density to the element
         const auto& element = localView.element();
@@ -152,7 +153,8 @@ namespace Dune::GFE
                                             std::vector<double>& localGradient,
                                             typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override
     {
-      LocalInterpolationRule localGFEFunction(localView.tree().finiteElement(),localCoefficients);
+      LocalInterpolationRule localGFEFunction;
+      localGFEFunction.bind(localView.tree().finiteElement(),localCoefficients);
 
       InterpolationDerivatives<LocalInterpolationRule> interpolationDerivatives(localGFEFunction,
                                                                                 localDensity_->dependsOnValue(),
@@ -532,8 +534,10 @@ namespace Dune::GFE
     using DeformationLocalInterpolationRule = typename std::tuple_element<0,LocalInterpolationRule>::type;
     using OrientationLocalInterpolationRule = typename std::tuple_element<1,LocalInterpolationRule>::type;
 
-    DeformationLocalInterpolationRule localDeformationGFEFunction(localView.tree().child(_0,0).finiteElement(),localCoefficients[_0]);
-    OrientationLocalInterpolationRule localOrientationGFEFunction(localView.tree().child(_1,0).finiteElement(),localCoefficients[_1]);
+    DeformationLocalInterpolationRule localDeformationGFEFunction;
+    OrientationLocalInterpolationRule localOrientationGFEFunction;
+    localDeformationGFEFunction.bind(localView.tree().child(_0,0).finiteElement(),localCoefficients[_0]);
+    localOrientationGFEFunction.bind(localView.tree().child(_1,0).finiteElement(),localCoefficients[_1]);
 
     InterpolationDerivatives<DeformationLocalInterpolationRule> deformationInterpolationDerivatives(localDeformationGFEFunction,
                                                                                                     localDensity_->dependsOnValue(0),
diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index e3f43afff597ba0ebf60dc7c1eab7ef98550ce9e..f4e30ac65fb31db6efe4fa46d293e7ca7664f309 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -175,7 +175,8 @@ namespace Dune::GFE
     //  Set up the local nonlinear finite element function
     ////////////////////////////////////////////////////////////////////////////////////
     typedef LocalGeodesicFEFunction<decltype(scalarBasis), TargetSpace> LocalGFEFunctionType;
-    LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+    LocalGFEFunctionType localGeodesicFEFunction;
+    localGeodesicFEFunction.bind(localFiniteElement,localSolution);
 
     // Bind the density to the current element
     density_->bind(element);
@@ -285,8 +286,10 @@ namespace Dune::GFE
     ////////////////////////////////////////////////////////////////////////////////////
     typedef LocalGeodesicFEFunction<decltype(deformationScalarBasis), RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
     typedef LocalGeodesicFEFunction<decltype(orientationScalarBasis), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
-    LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localConfiguration[_0]);
-    LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localConfiguration[_1]);
+    LocalDeformationGFEFunctionType localDeformationGFEFunction;
+    LocalOrientationGFEFunctionType localOrientationGFEFunction;
+    localDeformationGFEFunction.bind(deformationLocalFiniteElement,localConfiguration[_0]);
+    localOrientationGFEFunction.bind(orientationLocalFiniteElement,localConfiguration[_1]);
 
     // Bind the density to the current element
     density_->bind(element);
diff --git a/dune/gfe/assemblers/simofoxenergy.hh b/dune/gfe/assemblers/simofoxenergy.hh
index 821dd1d198fefa6a4413ec6c89947754914cfc96..75ff2888ac5cc01cdd31dfe83a3833892f988bbc 100644
--- a/dune/gfe/assemblers/simofoxenergy.hh
+++ b/dune/gfe/assemblers/simofoxenergy.hh
@@ -286,11 +286,17 @@ namespace Dune::GFE
     using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<decltype(midSurfaceBasis), RealTuple<double, 3> >;
     using LocalDirectorReferenceFunctionType   = LocalFEFunction<decltype(directorBasis), UnitVector<double, 3> >;
 
-    const LocalMidSurfaceFunctionType localMidSurfaceFunction(midSurfaceElement, localMidSurfaceConfiguration);
-    const LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(midSurfaceElement, displacements);
-    const LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction(midSurfaceElement, localRefMidSurfaceConfiguration);
-    const LocalDirectorFunctionType localDirectorFunction(directorElement, localDirectorConfiguration);
-    const LocalDirectorReferenceFunctionType localDirectorReferenceFunction(directorElement, localRefDirectorConfiguration);
+    LocalMidSurfaceFunctionType localMidSurfaceFunction;
+    LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction;
+    LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction;
+    LocalDirectorFunctionType localDirectorFunction;
+    LocalDirectorReferenceFunctionType localDirectorReferenceFunction;
+
+    localMidSurfaceFunction.bind(midSurfaceElement, localMidSurfaceConfiguration);
+    localMidSurfaceDisplacementFunction.bind(midSurfaceElement, displacements);
+    localMidSurfaceReferenceFunction.bind(midSurfaceElement, localRefMidSurfaceConfiguration);
+    localDirectorFunction.bind(directorElement, localDirectorConfiguration);
+    localDirectorReferenceFunction.bind(directorElement, localRefDirectorConfiguration);
 
     const int quadOrder = (element.type().isSimplex()) ? std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order())
                                                        : std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order()) + 1;
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index cf8ec4cf1eea7f0065629e71ea1e33f5ddf1fab8..3d49755a77904427a4f0c4bb11d3d31bcd32c2cf 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -87,8 +87,10 @@ namespace Dune::GFE
 
       typedef LocalGeodesicFEFunction<decltype(deformationScalarBasis), RBM0> LocalGFEFunctionType0;
       typedef LocalGeodesicFEFunction<decltype(rotationScalarBasis), RBM1> LocalGFEFunctionType1;
-      LocalGFEFunctionType0 localGeodesicFEFunction0(deformationLocalFiniteElement,localCoefficients[_0]);
-      LocalGFEFunctionType1 localGeodesicFEFunction1(orientationLocalFiniteElement,localCoefficients[_1]);
+      LocalGFEFunctionType0 localGeodesicFEFunction0;
+      LocalGFEFunctionType1 localGeodesicFEFunction1;
+      localGeodesicFEFunction0.bind(deformationLocalFiniteElement,localCoefficients[_0]);
+      localGeodesicFEFunction1.bind(orientationLocalFiniteElement,localCoefficients[_1]);
 
       RT energy = 0;
 
diff --git a/dune/gfe/assemblers/surfacecosseratstressassembler.hh b/dune/gfe/assemblers/surfacecosseratstressassembler.hh
index 56ea6d0b36eae127f9a12b0ed0315403f514320e..3f0b16cff90c1f1ea94c040d1d1bc1a5944c88ec 100644
--- a/dune/gfe/assemblers/surfacecosseratstressassembler.hh
+++ b/dune/gfe/assemblers/surfacecosseratstressassembler.hh
@@ -1,8 +1,6 @@
 #ifndef DUNE_GFE_SURFACECOSSERATSTRESSASSEMBLER_HH
 #define DUNE_GFE_SURFACECOSSERATSTRESSASSEMBLER_HH
 
-#include <dune/functions/functionspacebases/subspacebasis.hh>
-
 #include <dune/fufem/boundarypatch.hh>
 
 #include <dune/gfe/linearalgebra.hh>
@@ -238,7 +236,8 @@ namespace Dune::GFE
           for (std::size_t i=0; i<localConfigurationRot.size(); i++)
             localConfigurationRot[i] = rot[localViewOrderR.index(i)[0]];  //localViewOrderR.index(i) is a multiindex, its first entry is the actual index
 
-          LocalGFEFunctionR localGeodesicFEFunction(lFEOrderR,localConfigurationRot);
+          LocalGFEFunctionR localGeodesicFEFunction;
+          localGeodesicFEFunction.bind(lFEOrderR,localConfigurationRot);
 
           auto evaluateAtPoint = [&](FieldVector<double,3> pointGlobal, FieldVector<double,3> pointLocal3d) -> FieldMatrix<double,dim,dim> {
                                    Dune::FieldMatrix<double,dim,dim> nablaTheta;
diff --git a/dune/gfe/functions/discretekirchhoffbendingisometry.hh b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
index 25fcceb71e90935edf17157d329782d3f3638ef0..754f50d6ecc64cae430f23c72bbf4746edab790e 100644
--- a/dune/gfe/functions/discretekirchhoffbendingisometry.hh
+++ b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
@@ -233,7 +233,8 @@ namespace Dune::GFE
 
       //Create a LocalProjected-Finite element from the local coefficients used for interpolation.
       auto P1LagrangeLFE = localViewCoefficient_.tree().finiteElement();
-      LocalInterpolationRule localPBfunction(P1LagrangeLFE,localIsometryCoefficients);
+      LocalInterpolationRule localPBfunction;
+      localPBfunction.bind(P1LagrangeLFE,localIsometryCoefficients);
 
       /**
           Interpolate into the local hermite space for each component.
diff --git a/dune/gfe/functions/embeddedglobalgfefunction.hh b/dune/gfe/functions/embeddedglobalgfefunction.hh
index 44cb890589eda3c4dc8bb05a786d411057b95a3d..38ccb5c62d5c0a2824c4fa084e7a5c55a5433d6d 100644
--- a/dune/gfe/functions/embeddedglobalgfefunction.hh
+++ b/dune/gfe/functions/embeddedglobalgfefunction.hh
@@ -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
@@ -319,7 +319,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);
diff --git a/dune/gfe/functions/globalgfefunction.hh b/dune/gfe/functions/globalgfefunction.hh
index 49ceff7f4fe23e5de6e82e324ae403b30dc49626..73f3451ea72c7df82e0ee2def92784d5d16a09ca 100644
--- a/dune/gfe/functions/globalgfefunction.hh
+++ b/dune/gfe/functions/globalgfefunction.hh
@@ -160,8 +160,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(this->localView_.tree().finiteElement(),localDoFs_);
         }
 
         //! Unbind the local-function.
@@ -191,7 +190,7 @@ namespace Dune::GFE
 #endif
         LocalView localView_;
         std::vector<Coefficient> localDoFs_;
-        std::unique_ptr<LocalInterpolationRule> localInterpolationRule_;
+        LocalInterpolationRule localInterpolationRule_;
       };
 
     protected:
@@ -308,7 +307,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
@@ -477,7 +476,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);
diff --git a/dune/gfe/functions/interpolationderivatives.hh b/dune/gfe/functions/interpolationderivatives.hh
index 4a0cd5d5defa87d97be17fb9ac1e3cabc9c3ad05..18caba50b67d99d69c195e0d65edde08d59883f0 100644
--- a/dune/gfe/functions/interpolationderivatives.hh
+++ b/dune/gfe/functions/interpolationderivatives.hh
@@ -224,7 +224,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;
+      localGFEFunction.bind(scalarFiniteElement,localAConfiguration);
 
       if (doValue_)
       {
diff --git a/dune/gfe/functions/localgeodesicfefunction.hh b/dune/gfe/functions/localgeodesicfefunction.hh
index 8424b3c4d55ab2e9968e8fcf55fe66d4bda1b66f..d405286d93e3e1710d34094fc93921638bb3bc89 100644
--- a/dune/gfe/functions/localgeodesicfefunction.hh
+++ b/dune/gfe/functions/localgeodesicfefunction.hh
@@ -60,16 +60,19 @@ namespace Dune::GFE
     /** \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
+    /** \brief Bind the local function to a particular scalar finite element
+     * and a set of coefficients
+     *
+     * \param localFiniteElement A scalar finite element that provides the weight functions
+     * \param coefficients Values to be interpolated
      */
-    LocalGeodesicFEFunction(const LocalFiniteElement& localFiniteElement,
-                            const std::vector<TargetSpace>& coefficients)
-      : localFiniteElement_(localFiniteElement),
-      coefficients_(coefficients)
+    void bind(const LocalFiniteElement& localFiniteElement,
+              const std::vector<TargetSpace>& coefficients)
     {
-      assert(localFiniteElement_.localBasis().size() == coefficients_.size());
+      assert(localFiniteElement.localBasis().size() == coefficients.size());
+
+      localFiniteElement_ = localFiniteElement;
+      coefficients_ = coefficients;
     }
 
     /** \brief Rebind the FEFunction to another TargetSpace */
@@ -206,7 +209,7 @@ namespace Dune::GFE
     /** \brief The scalar local finite element, which provides the weighting factors
         \todo We really only need the local basis
      */
-    const LocalFiniteElement& localFiniteElement_;
+    LocalFiniteElement localFiniteElement_;
 
     /** \brief The coefficient vector */
     std::vector<TargetSpace> coefficients_;
@@ -411,8 +414,9 @@ namespace Dune::GFE
       cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
       cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
 
-      const LocalGeodesicFEFunction<Basis,TargetSpace> fPlus(localFiniteElement_,cornersPlus);
-      const LocalGeodesicFEFunction<Basis,TargetSpace> fMinus(localFiniteElement_,cornersMinus);
+      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus, fMinus;
+      fPlus.bind(localFiniteElement_,cornersPlus);
+      fMinus.bind(localFiniteElement_,cornersMinus);
 
       TargetSpace hPlus  = fPlus.evaluate(local);
       TargetSpace hMinus = fMinus.evaluate(local);
@@ -550,8 +554,9 @@ namespace Dune::GFE
       cornersPlus[coefficient]  = TargetSpace(aPlus);
       cornersMinus[coefficient] = TargetSpace(aMinus);
 
-      const LocalGeodesicFEFunction<Basis,TargetSpace> fPlus(localFiniteElement_,cornersPlus);
-      const LocalGeodesicFEFunction<Basis,TargetSpace> fMinus(localFiniteElement_,cornersMinus);
+      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus,fMinus;
+      fPlus.bind(localFiniteElement_,cornersPlus);
+      fMinus.bind(localFiniteElement_,cornersMinus);
 
       const auto hPlus  = fPlus.evaluateDerivative(local);
       const auto hMinus = fMinus.evaluateDerivative(local);
@@ -618,16 +623,18 @@ namespace Dune::GFE
     /** \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())
+    /** \brief Bind the function to a particular weight function set and coefficients
+     */
+    void bind(const LocalFiniteElement& localFiniteElement,
+              const std::vector<TargetSpace>& coefficients)
     {
       using namespace Dune::Indices;
       assert(localFiniteElement.localBasis().size() == coefficients.size());
 
+      localFiniteElement_ = localFiniteElement;
+      coefficients_ = coefficients;
+
+      translationCoefficients_.resize(coefficients.size());
       for (size_t i=0; i<coefficients.size(); i++)
         translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
 
@@ -635,7 +642,8 @@ namespace Dune::GFE
       for (size_t i=0; i<coefficients.size(); i++)
         orientationCoefficients[i] = coefficients[i][_1];
 
-      orientationFEFunction_ = std::make_unique<LocalGeodesicFEFunction<Basis,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
+      orientationFEFunction_ = std::make_unique<LocalGeodesicFEFunction<Basis,Rotation<field_type,3> > >();
+      orientationFEFunction_->bind(localFiniteElement,orientationCoefficients);
     }
 
     /** \brief Rebind the FEFunction to another TargetSpace */
@@ -831,7 +839,7 @@ namespace Dune::GFE
     /** \brief The scalar local finite element, which provides the weighting factors
         \todo We really only need the local basis
      */
-    const LocalFiniteElement& localFiniteElement_;
+    LocalFiniteElement localFiniteElement_;
 
     // The coefficients of this interpolation rule
     std::vector<TargetSpace> coefficients_;
diff --git a/dune/gfe/functions/localprojectedfefunction.hh b/dune/gfe/functions/localprojectedfefunction.hh
index f17e3229923a8f126e14eb6feb4c1c7362cf9c23..fa6c6bd1e60ec8d11115e6edf00fdc4d8f059ee6 100644
--- a/dune/gfe/functions/localprojectedfefunction.hh
+++ b/dune/gfe/functions/localprojectedfefunction.hh
@@ -50,16 +50,18 @@ namespace Dune::GFE
     /** \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
+    /** \brief Bind the function to a specific finite element interpolation space and set of coefficients
+     *
+     * \param localFiniteElement A finite element that defines the interpolation in the embedding space
+     * \param coefficients Coefficients of the finite element function
      */
-    LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
-                             const std::vector<TargetSpace>& coefficients)
-      : localFiniteElement_(localFiniteElement),
-      coefficients_(coefficients)
+    void bind(const LocalFiniteElement& localFiniteElement,
+              const std::vector<TargetSpace>& coefficients)
     {
-      assert(localFiniteElement_.localBasis().size() == coefficients_.size());
+      assert(localFiniteElement.localBasis().size() == coefficients.size());
+
+      localFiniteElement_ = localFiniteElement;
+      coefficients_ = coefficients;
     }
 
     /** \brief Rebind the FEFunction to another TargetSpace */
@@ -130,7 +132,7 @@ namespace Dune::GFE
     /** \brief The scalar local finite element, which provides the weighting factors
      *        \todo We really only need the local basis
      */
-    const LocalFiniteElement& localFiniteElement_;
+    LocalFiniteElement localFiniteElement_;
 
     /** \brief The coefficient vector */
     std::vector<TargetSpace> coefficients_;
@@ -358,16 +360,18 @@ namespace Dune::GFE
     /** \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
+    /** \brief Bind the function to a specific finite element interpolation space and set of coefficients
+     *
+     * \param localFiniteElement A finite element that defines the interpolation in the embedding space
+     * \param coefficients Coefficients of the finite element function
      */
-    LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
-                             const std::vector<TargetSpace>& coefficients)
-      : localFiniteElement_(localFiniteElement),
-      coefficients_(coefficients)
+    void bind(const LocalFiniteElement& localFiniteElement,
+              const std::vector<TargetSpace>& coefficients)
     {
-      assert(localFiniteElement_.localBasis().size() == coefficients_.size());
+      assert(localFiniteElement.localBasis().size() == coefficients.size());
+
+      localFiniteElement_ = localFiniteElement;
+      coefficients_ = coefficients;
     }
 
     /** \brief Rebind the FEFunction to another TargetSpace */
@@ -542,7 +546,7 @@ namespace Dune::GFE
     /** \brief The scalar local finite element, which provides the weighting factors
      *        \todo We really only need the local basis
      */
-    const LocalFiniteElement& localFiniteElement_;
+    LocalFiniteElement localFiniteElement_;
 
     /** \brief The coefficient vector */
     std::vector<TargetSpace> coefficients_;
@@ -583,19 +587,21 @@ namespace Dune::GFE
     /** \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
+    /** \brief Bind the function to a specific finite element interpolation space and set of coefficients
+     *
+     * \param localFiniteElement A finite element that defines the interpolation in the embedding space
+     * \param coefficients Coefficients of the finite element function
      */
-    LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
-                             const std::vector<TargetSpace>& coefficients)
-      : localFiniteElement_(localFiniteElement),
-      coefficients_(coefficients),
-      translationCoefficients_(coefficients.size())
+    void bind(const LocalFiniteElement& localFiniteElement,
+              const std::vector<TargetSpace>& coefficients)
     {
       assert(localFiniteElement.localBasis().size() == coefficients.size());
       using namespace Dune::Indices;
 
+      localFiniteElement_ = localFiniteElement;
+      coefficients_ = coefficients;
+
+      translationCoefficients_.resize(coefficients.size());
       for (size_t i=0; i<coefficients.size(); i++)
         translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
 
@@ -603,7 +609,8 @@ namespace Dune::GFE
       for (size_t i=0; i<coefficients.size(); i++)
         orientationCoefficients[i] = coefficients[i][_1];
 
-      orientationFunction_ = std::make_unique<LocalProjectedFEFunction<Basis,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
+      orientationFunction_ = std::make_unique<LocalProjectedFEFunction<Basis,Rotation<field_type,3> > > ();
+      orientationFunction_->bind(localFiniteElement,orientationCoefficients);
     }
 
     /** \brief Rebind the FEFunction to another TargetSpace */
@@ -694,7 +701,7 @@ namespace Dune::GFE
     /** \brief The scalar local finite element, which provides the weighting factors
      *        \todo We really only need the local basis
      */
-    const LocalFiniteElement& localFiniteElement_;
+    LocalFiniteElement localFiniteElement_;
 
     // The coefficients of this interpolation rule
     std::vector<TargetSpace> coefficients_;
diff --git a/test/interpolationderivativestest.cc b/test/interpolationderivativestest.cc
index 3860337a9116fdac0390663ea3183e3d419149f9..480d26cc92038a4623e67b7460cdb5990b767ce6 100644
--- a/test/interpolationderivativestest.cc
+++ b/test/interpolationderivativestest.cc
@@ -143,8 +143,9 @@ public:
         cornersPlus [coefficient] = TargetSpace(coefficients_[coefficient].globalCoordinates() + variation);
         cornersMinus[coefficient] = TargetSpace(coefficients_[coefficient].globalCoordinates() - variation);
 
-        LocalInterpolationRule fPlus(localInterpolationRule_.localFiniteElement(),cornersPlus);
-        LocalInterpolationRule fMinus(localInterpolationRule_.localFiniteElement(),cornersMinus);
+        LocalInterpolationRule fPlus, fMinus;
+        fPlus.bind(localInterpolationRule_.localFiniteElement(),cornersPlus);
+        fMinus.bind(localInterpolationRule_.localFiniteElement(),cornersMinus);
 
         /////////////////////////////////////////////////////////////
         //  Compute first derivative of the interpolation value
@@ -195,8 +196,9 @@ public:
         cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
         cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
 
-        LocalInterpolationRule fPlus(localInterpolationRule_.localFiniteElement(),cornersPlus);
-        LocalInterpolationRule fMinus(localInterpolationRule_.localFiniteElement(),cornersMinus);
+        LocalInterpolationRule fPlus, fMinus;
+        fPlus.bind(localInterpolationRule_.localFiniteElement(),cornersPlus);
+        fMinus.bind(localInterpolationRule_.localFiniteElement(),cornersMinus);
 
         /////////////////////////////////////////////////////////////
         //  Compute first derivative of the interpolation value
@@ -254,8 +256,9 @@ public:
         forwardSolution[i]  = TargetSpace::exp(coefficients_[i], eps * xi);
         backwardSolution[i] = TargetSpace::exp(coefficients_[i], -1 * eps * xi);
 
-        LocalInterpolationRule fPlus(localInterpolationRule_.localFiniteElement(),forwardSolution);
-        LocalInterpolationRule fMinus(localInterpolationRule_.localFiniteElement(),backwardSolution);
+        LocalInterpolationRule fPlus, fMinus;
+        fPlus.bind(localInterpolationRule_.localFiniteElement(),forwardSolution);
+        fMinus.bind(localInterpolationRule_.localFiniteElement(),backwardSolution);
 
         forwardValue[i][i2] = fPlus.evaluate(localPos_);
         backwardValue[i][i2] = fMinus.evaluate(localPos_);
@@ -312,8 +315,9 @@ public:
               backwardSolutionXiEta[j] = TargetSpace::exp(coefficients_[j], (-1)*epsEta);
             }
 
-            LocalInterpolationRule fPlus(localInterpolationRule_.localFiniteElement(),forwardSolutionXiEta);
-            LocalInterpolationRule fMinus(localInterpolationRule_.localFiniteElement(),backwardSolutionXiEta);
+            LocalInterpolationRule fPlus, fMinus;
+            fPlus.bind(localInterpolationRule_.localFiniteElement(),forwardSolutionXiEta);
+            fMinus.bind(localInterpolationRule_.localFiniteElement(),backwardSolutionXiEta);
 
             /////////////////////////////////////////////////////////////////////////////////////
             //  Compute second derivative of the adjoint vector times the interpolation value
@@ -418,7 +422,8 @@ TestSuite checkDerivatives()
 
   auto localView = scalarBasis.localView();
   localView.bind(*gridView.begin<0>());
-  LocalInterpolationRule localGFEFunction(localView.tree().finiteElement(),localCoefficients);
+  LocalInterpolationRule localGFEFunction;
+  localGFEFunction.bind(localView.tree().finiteElement(),localCoefficients);
 
   GFE::InterpolationDerivatives<LocalInterpolationRule> interpolationDerivatives(localGFEFunction,
                                                                                  true,   // doValue
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index bf0187b6abe5724e89fcf89ad743616580d75c9a..30f1ce3926d2e7a7e721160829e29845afe9509a 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -324,8 +324,8 @@ void test(const GeometryType& element)
     localBasisView.bind(*gridView.template begin<0>());
     const auto& localFiniteElement = localBasisView.tree().finiteElement();
 
-    GFE::LocalGeodesicFEFunction<InterpolationBasis,TargetSpace> f(localFiniteElement,
-                                                                   corners);
+    GFE::LocalGeodesicFEFunction<InterpolationBasis,TargetSpace> f;
+    f.bind(localFiniteElement,corners);
 
     //testPermutationInvariance(corners);
     testDerivative(f);
diff --git a/test/localprojectedfefunctiontest.cc b/test/localprojectedfefunctiontest.cc
index 8b5473bf3bcc0e8727c4a9b1227b1f4bb39646ee..56143cf537bf2100a76a830bee3d2df79cacd683 100644
--- a/test/localprojectedfefunctiontest.cc
+++ b/test/localprojectedfefunctiontest.cc
@@ -281,8 +281,10 @@ void test(const GeometryType& element)
     localBasisView.bind(*gridView.template begin<0>());
     const auto& localFiniteElement = localBasisView.tree().finiteElement();
 
-    GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f(localFiniteElement,corners);
-    GFE::LocalProjectedFEFunction<InterpolationBasis, TargetSpace,false> f_nonconforming(localFiniteElement,corners);
+    GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f;
+    f.bind(localFiniteElement,corners);
+    GFE::LocalProjectedFEFunction<InterpolationBasis, TargetSpace,false> f_nonconforming;
+    f_nonconforming.bind(localFiniteElement,corners);
 
     //testPermutationInvariance(corners);
     testDerivative(f);