diff --git a/CHANGELOG.md b/CHANGELOG.md
index 556901d5fb71f9468ef36b402fb3c8ced40c84da..ccc354c82fb15ad0b4b2a2f6f74c07a5138882bc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,10 @@
   does not accept a surface load density anymore. Use the
   `NeumannEnergy` class for assembling surface loads.
 
+- The class `LocalIntegralEnergy` now does not construct the local interpolation
+  rule itself anymore.  Rather, it expects to be given one in its constructor.
+  This will allow interpolation rules with state.
+
 - The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
   now take a `dune-functions` basis as their first arguments, instead
   of a `LocalFiniteElement`.
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index f0656a40439eb6b049798b696482053b48c6d9c5..03d49c76a129d8a4037150d75ca9d7ce6466ba52 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -61,10 +61,28 @@ namespace Dune::GFE
 
   public:
 
-    /** \brief Constructor with a Dune::GFE::LocalDensity
+    /** \brief Constructor from a GFE function as a shared pointer
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param ld The density that will be integrated
      */
-    LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
-      : localDensity_(ld)
+    LocalIntegralEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+                        const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
+      : localGFEFunction_(localGFEFunction),
+      localDensity_(ld)
+    {}
+
+    /** \brief Constructor from a GFE function r-value reference
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param ld The density that will be integrated
+     */
+    LocalIntegralEnergy(LocalInterpolationRule&& localGFEFunction,
+                        const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
+      : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction))),
+      localDensity_(ld)
     {}
 
   private:
@@ -96,8 +114,7 @@ namespace Dune::GFE
 #else
         const auto& localFiniteElement = Impl::LocalFiniteElementFactory<Basis>::get(localView);
 #endif
-        LocalInterpolationRule localInterpolationRule;
-        localInterpolationRule.bind(localFiniteElement,localConfiguration);
+        localGFEFunction_->bind(localFiniteElement,localConfiguration);
 
         // Bind density to the element
         const auto& element = localView.element();
@@ -123,13 +140,13 @@ namespace Dune::GFE
           {
             if (localDensity_->dependsOnDerivative())
             {
-              auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos);
+              auto [value, derivative] = localGFEFunction_->evaluateValueAndDerivative(quadPos);
               derivative = derivative * geometryJacobianInverse;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
             }
             else
             {
-              auto value = localInterpolationRule.evaluate(quadPos);
+              const auto value = localGFEFunction_->evaluate(quadPos);
               typename LocalInterpolationRule::DerivativeType dummyDerivative;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
             }
@@ -139,7 +156,7 @@ namespace Dune::GFE
             if (localDensity_->dependsOnDerivative())
             {
               typename TargetSpace::CoordinateType dummyValue;
-              auto derivative = localInterpolationRule.evaluateDerivative(quadPos);
+              auto derivative = localGFEFunction_->evaluateDerivative(quadPos);
               derivative = derivative * geometryJacobianInverse;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
             }
@@ -177,13 +194,8 @@ namespace Dune::GFE
         const auto& localFiniteElement0 = localView.tree().child(_0,0).finiteElement();
         const auto& localFiniteElement1 = localView.tree().child(_1,0).finiteElement();
 
-        using LocalGFEFunctionType0 = typename std::tuple_element<0, LocalInterpolationRule>::type;
-        using LocalGFEFunctionType1 = typename std::tuple_element<1, LocalInterpolationRule>::type;
-
-        LocalGFEFunctionType0 localGFEFunction0;
-        LocalGFEFunctionType1 localGFEFunction1;
-        localGFEFunction0.bind(localFiniteElement0, coefficients[_0]);
-        localGFEFunction1.bind(localFiniteElement1, coefficients[_1]);
+        std::get<0>(*localGFEFunction_).bind(localFiniteElement0, coefficients[_0]);
+        std::get<1>(*localGFEFunction_).bind(localFiniteElement1, coefficients[_1]);
 
         // Bind density to the element
         const auto& element = localView.element();
@@ -212,19 +224,19 @@ namespace Dune::GFE
           // A value is needed either directly, or for computing the derivative.
           TargetSpace value;
           if (localDensity_->dependsOnValue(0) || localDensity_->dependsOnDerivative(0))
-            value[_0] = localGFEFunction0.evaluate(quadPos);
+            value[_0] = std::get<0>(*localGFEFunction_).evaluate(quadPos);
           if (localDensity_->dependsOnValue(1) || localDensity_->dependsOnDerivative(1))
-            value[_1] = localGFEFunction1.evaluate(quadPos);
+            value[_1] = std::get<1>(*localGFEFunction_).evaluate(quadPos);
 
           // Compute the derivatives of the interpolation function factors
-          typename LocalGFEFunctionType0::DerivativeType derivative0;
-          typename LocalGFEFunctionType1::DerivativeType derivative1;
+          typename std::tuple_element_t<0, LocalInterpolationRule>::DerivativeType derivative0;
+          typename std::tuple_element_t<1, LocalInterpolationRule>::DerivativeType derivative1;
 
           if (localDensity_->dependsOnDerivative(0))
-            derivative0 = localGFEFunction0.evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
+            derivative0 = std::get<0>(*localGFEFunction_).evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
 
           if (localDensity_->dependsOnDerivative(1))
-            derivative1 = localGFEFunction1.evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
+            derivative1 = std::get<1>(*localGFEFunction_).evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
 
           // Copy the two derivatives into a joint matrix object
           // TODO: I am not sure about this.  May the densities should get the
@@ -249,6 +261,11 @@ namespace Dune::GFE
     }
 
   protected:
+
+    // The value and derivative of this function are evaluated at the quadrature points,
+    // and given to the density.
+    const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
+
     const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> > localDensity_;
   };
 
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 858e2e138c6bfac97bfaab2ee989ec05bc1eaae7..3c2a8e59a584955516e02fb56f274eca4ada0e21 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -115,7 +115,8 @@ using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<
 //
 // TODO: Do we really need the 'creator' argument?
 template <typename Basis, typename InterpolationRule, typename TargetSpace, typename Creator>
-auto createCosseratEnergy(const ParameterTree& materialParameters,
+auto createCosseratEnergy(std::shared_ptr<InterpolationRule>& localGFEFunction,
+                          const ParameterTree& materialParameters,
                           const Creator creator)
 {
   using GridView = typename Basis::GridView;
@@ -128,7 +129,7 @@ auto createCosseratEnergy(const ParameterTree& materialParameters,
   {
     auto density = std::make_shared<GFE::PlanarCosseratShellDensity<Element, adouble> >(materialParameters);
 
-    return std::make_shared<GFE::LocalIntegralEnergy<Basis,InterpolationRule,TargetSpace> >(density);
+    return std::make_shared<GFE::LocalIntegralEnergy<Basis,InterpolationRule,TargetSpace> >(localGFEFunction, density);
   }
   else if constexpr (dim==2 && dimworld==3)
   {
@@ -141,7 +142,7 @@ auto createCosseratEnergy(const ParameterTree& materialParameters,
   {
     auto density = std::make_shared<GFE::BulkCosseratDensity<Element, adouble> >(materialParameters);
 
-    return std::make_shared<GFE::LocalIntegralEnergy<Basis,InterpolationRule,TargetSpace> >(density);
+    return std::make_shared<GFE::LocalIntegralEnergy<Basis,InterpolationRule,TargetSpace> >(localGFEFunction, density);
   }
 }
 
@@ -575,6 +576,7 @@ int main (int argc, char *argv[]) try
     // Construct the interpolation rule, i.e., the geometric finite element
     using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >,
         GFE::LocalGeodesicFEFunction<OrientationFEBasis, GFE::Rotation<adouble,3> > >;
+    auto localGFEFunction = std::make_shared<AInterpolationRule>();
 
     using ATargetSpace = typename TargetSpace::rebind<adouble>::other;
 
@@ -582,8 +584,9 @@ int main (int argc, char *argv[]) try
     auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<adouble,3>,GFE::Rotation<adouble,3> > >();
 
     // The actual Cosserat energy
-    auto localCosseratEnergy = createCosseratEnergy<CompositeBasis,AInterpolationRule,ATargetSpace,decltype(creator)>(materialParameters,
-                                                                                                                      creator);
+    auto localCosseratEnergy = createCosseratEnergy<CompositeBasis,AInterpolationRule,ATargetSpace>(localGFEFunction,
+                                                                                                    materialParameters,
+                                                                                                    creator);
 
     sumEnergy->addLocalEnergy(localCosseratEnergy);
 
@@ -594,7 +597,7 @@ int main (int argc, char *argv[]) try
     // The volume load term
     using Element = GridView::Codim<0>::Entity;
     auto volumeLoadDensity = std::make_shared<GFE::CosseratVolumeLoadDensity<Element,adouble> >(volumeLoad);
-    auto volumeLoadEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, AInterpolationRule, ATargetSpace> >(volumeLoadDensity);
+    auto volumeLoadEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, AInterpolationRule, ATargetSpace> >(localGFEFunction, volumeLoadDensity);
     sumEnergy->addLocalEnergy(volumeLoadEnergy);
 
     // The local assembler
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index b8e3eed544ac6b0624c111f0115430dced78bf4b..f1526a5a85a619c19eb1c4254e4a8a22bf934ed9 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -527,8 +527,9 @@ int main (int argc, char *argv[]) try
     using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
+    LocalInterpolationRule localGFEFunction;
 
-    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
+    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(std::move(localGFEFunction), elasticDensityWrapped);
     auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<ValueType,targetDim>, GFE::Rotation<ValueType,dim> > >(neumannBoundary,neumannFunctionPtr);
 
     // The energy of the surface shell
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 1fec4da9dc5a97c3546312fa8880d7ba8c4f059f..c98c6a4c6ac195d86d7e02014cc44e4e45b5ca0b 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -210,9 +210,10 @@ int main (int argc, char *argv[]) try
 
   std::vector<std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > > addends(2);
   using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis, ATargetSpace>;
+  GeodesicInterpolationRule localGFEFunctionA;
 
   auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity, ATargetSpace> >();
-  addends[0] = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(harmonicDensity);
+  addends[0] = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(std::move(localGFEFunctionA), harmonicDensity);
   addends[1] = l2DistanceSquaredEnergy;
 
   std::vector<double> weights = {1.0, 1.0/(2*timeStepSize)};
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 667b184b833d0e3d7a8748a67f9afe0abcdfd22b..18b387fad6dfff3228c5c80ac02ba955f901780c 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -290,13 +290,15 @@ int main (int argc, char *argv[])
   // Next: The local energy, i.e., the integral of the density over one element
   using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<FEBasis, ATargetSpace>;
   using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis, ATargetSpace>;
+  GeodesicInterpolationRule localGeodesicFEFunctionA;
+  ProjectedInterpolationRule localProjectedFEFunctionA;
 
   std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
 
   if (parameterSet["interpolationMethod"] == "geodesic")
-    localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(density);
+    localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(std::move(localGeodesicFEFunctionA), density);
   else if (parameterSet["interpolationMethod"] == "projected")
-    localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >(density);
+    localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >(std::move(localProjectedFEFunctionA), density);
   else
     DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
 
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index 3af4851ba2415bfec4b2a812348235c7defff047..d9d4096d930dbc44e97fb8d7071158e79fc96ac5 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -209,10 +209,11 @@ int main (int argc, char *argv[])
   using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
   using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
+  LocalInterpolationRule localGFEFunctionA;
 
   using Element = GridView::Codim<0>::Entity;
   auto bulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<Element,adouble> >(parameters);
-  auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ARigidBodyMotion> >(bulkCosseratDensity);
+  auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ARigidBodyMotion> >(std::move(localGFEFunctionA), bulkCosseratDensity);
   sumEnergy->addLocalEnergy(bulkCosseratEnergy);
   sumEnergy->addLocalEnergy(neumannEnergy);
 
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index f6d3117e5484318a92e593dd5df2cfdc5cabacdd..a7e8d902007f1964e96a2e6b451ac0f34fd909ae 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -353,8 +353,9 @@ int main (int argc, char *argv[])
   using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
   using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
+  LocalInterpolationRule localGFEFunction;
 
-  auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
+  auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(std::move(localGFEFunction), elasticDensityWrapped);
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<ValueType,targetDim>, GFE::Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
 
   using Intersection = typename GridView::Intersection;
diff --git a/test/localadolcstiffnesstest.cc b/test/localadolcstiffnesstest.cc
index ad44807ecdcbc8166673688f21c34f26dbaab28d..cb63d4bea3008f3a2c40190e71f7c14e293524a1 100644
--- a/test/localadolcstiffnesstest.cc
+++ b/test/localadolcstiffnesstest.cc
@@ -188,10 +188,11 @@ int main (int argc, char *argv[]) try
 
   // Select geometric finite element interpolation method
   using AInterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis, ATargetSpace>;
+  AInterpolationRule localGFEFunctionA;
 
   auto activeDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(materialParameters);
 
-  auto activeCosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,AInterpolationRule,ATargetSpace> >(activeDensity);
+  auto activeCosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,AInterpolationRule,ATargetSpace> >(std::move(localGFEFunctionA), activeDensity);
 
   // The actual assembler
   GFE::LocalGeodesicFEADOLCStiffness<TangentBasis,
@@ -203,10 +204,11 @@ int main (int argc, char *argv[]) try
 
   // Select geometric finite element interpolation method
   using InterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis, TargetSpace>;
+  InterpolationRule localGFEFunction;
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, double> >(materialParameters);
 
-  auto cosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,InterpolationRule,TargetSpace> >(cosseratDensity);
+  auto cosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,InterpolationRule,TargetSpace> >(std::move(localGFEFunction), cosseratDensity);
 
   // The actual assembler
   GFE::LocalGeodesicFEFDStiffness<TangentBasis,
diff --git a/test/localintegralstiffnesstest.cc b/test/localintegralstiffnesstest.cc
index f44ef299ff98f896ff67aceddd5def6b4f7e3c89..10783d17ce85c595468d50bc968b70539d46b78a 100644
--- a/test/localintegralstiffnesstest.cc
+++ b/test/localintegralstiffnesstest.cc
@@ -112,9 +112,10 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
 
     using ALocalInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(feBasis),
         ATargetSpace>;
+    ALocalInterpolationRule localGFEFunctionA;
 
     // Assemble using the old assembler
-    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensityA);
+    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(std::move(localGFEFunctionA), harmonicDensityA);
 
     auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
     assemblerADOLC = std::make_shared<GFE::GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
@@ -134,9 +135,10 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
 
     using ALocalInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis,
         ATargetSpace, interpolationType!=Nonconforming>;
+    ALocalInterpolationRule localGFEFunctionA;
 
     // Assemble using the old assembler
-    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensityA);
+    auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(std::move(localGFEFunctionA), harmonicDensityA);
 
     auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
     assemblerADOLC = std::make_shared<GFE::GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
@@ -322,9 +324,10 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     using ALocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
+    ALocalInterpolationRule localGFEFunctionA;
 
     // Assemble using the ADOL-C assembler
-    auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(aBulkCosseratDensity);
+    auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(std::move(localGFEFunctionA), aBulkCosseratDensity);
 
     auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
 
@@ -346,9 +349,10 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
+    ALocalInterpolationRule localGFEFunctionA;
 
     // Assemble using the ADOL-C assembler
-    auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(aBulkCosseratDensity);
+    auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(std::move(localGFEFunctionA), aBulkCosseratDensity);
 
     auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
 
diff --git a/test/mixedriemannianpnsolvertest.cc b/test/mixedriemannianpnsolvertest.cc
index 32acd8780be07de0b9b72ca62b05abce9194389e..0ddaa68dec41e2a9db7a23ca5c6523dee11d4ebc 100644
--- a/test/mixedriemannianpnsolvertest.cc
+++ b/test/mixedriemannianpnsolvertest.cc
@@ -174,10 +174,11 @@ int main (int argc, char *argv[])
   // The Cosserat shell energy
   using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >,
       GFE::LocalGeodesicFEFunction<RotationFEBasis, GFE::Rotation<adouble,3> > >;
+  AInterpolationRule localGFEFunctionA;
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
 
-  auto planarCosseratShellEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ARBM> >(cosseratDensity);
+  auto planarCosseratShellEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ARBM> >(std::move(localGFEFunctionA), cosseratDensity);
 
   sumEnergy->addLocalEnergy(planarCosseratShellEnergy);
 
diff --git a/test/planarcosseratshelltest.cc b/test/planarcosseratshelltest.cc
index a79d6765aec85549a8d1d5a7cb35ff21af90ef08..0db975787afe8656d89f6ba98561ddc85358c503 100644
--- a/test/planarcosseratshelltest.cc
+++ b/test/planarcosseratshelltest.cc
@@ -201,10 +201,11 @@ int main (int argc, char *argv[])
   // The Cosserat shell energy
   using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >,
       GFE::LocalGeodesicFEFunction<OrientationFEBasis, GFE::Rotation<adouble,3> > >;
+  AInterpolationRule localGFEFunctionA;
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
 
-  auto planarCosseratShellEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ARigidBodyMotion> >(cosseratDensity);
+  auto planarCosseratShellEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ARigidBodyMotion> >(std::move(localGFEFunctionA), cosseratDensity);
 
   sumEnergy->addLocalEnergy(planarCosseratShellEnergy);