diff --git a/CHANGELOG.md b/CHANGELOG.md
index f452d9055c0fd7990f1adf7c8824f81c6d5ba3b1..fb0526297f23a1807dc9b5e1b8297ef1bd73e083 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,7 +15,10 @@
 
 - The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
   now take a `dune-functions` basis as their first arguments, instead
-  of a `LocalFiniteElement`.
+  of a `LocalFiniteElement`. In addition, they need to be given an object
+  of this type in the constructor.  This is the basis to be used for
+  interpolation weights (for `LocalGeodesicFEFunction`) and for
+  interpolating in the embedding space (for `LocalProjectedFEFunction`).
 
 - The entire code is now in the `Dune::GFE` namespace.
 
diff --git a/dune/gfe/assemblers/l2distancesquaredenergy.hh b/dune/gfe/assemblers/l2distancesquaredenergy.hh
index 4d625a47446ac4a2766cc3fe77a1b95b45219470..bb8b4db193cd0ccdd46a1f5fd5d71ffd777470a5 100644
--- a/dune/gfe/assemblers/l2distancesquaredenergy.hh
+++ b/dune/gfe/assemblers/l2distancesquaredenergy.hh
@@ -41,7 +41,7 @@ namespace Dune::GFE
 
       const auto& localFiniteElement = localView.tree().finiteElement();
       typedef LocalGeodesicFEFunction<Basis, TargetSpace> LocalGFEFunctionType;
-      LocalGFEFunctionType localGeodesicFEFunction;
+      LocalGFEFunctionType localGeodesicFEFunction(localView.globalBasis());
       localGeodesicFEFunction.bind(localFiniteElement,localSolution);
 
       const auto element = localView.element();
diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index f4e30ac65fb31db6efe4fa46d293e7ca7664f309..6127c16eea0782292b3344705991edcb215539e3 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -175,7 +175,7 @@ namespace Dune::GFE
     //  Set up the local nonlinear finite element function
     ////////////////////////////////////////////////////////////////////////////////////
     typedef LocalGeodesicFEFunction<decltype(scalarBasis), TargetSpace> LocalGFEFunctionType;
-    LocalGFEFunctionType localGeodesicFEFunction;
+    LocalGFEFunctionType localGeodesicFEFunction(scalarBasis);
     localGeodesicFEFunction.bind(localFiniteElement,localSolution);
 
     // Bind the density to the current element
@@ -286,8 +286,8 @@ namespace Dune::GFE
     ////////////////////////////////////////////////////////////////////////////////////
     typedef LocalGeodesicFEFunction<decltype(deformationScalarBasis), RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
     typedef LocalGeodesicFEFunction<decltype(orientationScalarBasis), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
-    LocalDeformationGFEFunctionType localDeformationGFEFunction;
-    LocalOrientationGFEFunctionType localOrientationGFEFunction;
+    LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationScalarBasis);
+    LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationScalarBasis);
     localDeformationGFEFunction.bind(deformationLocalFiniteElement,localConfiguration[_0]);
     localOrientationGFEFunction.bind(orientationLocalFiniteElement,localConfiguration[_1]);
 
diff --git a/dune/gfe/assemblers/simofoxenergy.hh b/dune/gfe/assemblers/simofoxenergy.hh
index 75ff2888ac5cc01cdd31dfe83a3833892f988bbc..a00cae84495649fad5cdb45d258e540a2c261e3d 100644
--- a/dune/gfe/assemblers/simofoxenergy.hh
+++ b/dune/gfe/assemblers/simofoxenergy.hh
@@ -286,11 +286,11 @@ namespace Dune::GFE
     using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<decltype(midSurfaceBasis), RealTuple<double, 3> >;
     using LocalDirectorReferenceFunctionType   = LocalFEFunction<decltype(directorBasis), UnitVector<double, 3> >;
 
-    LocalMidSurfaceFunctionType localMidSurfaceFunction;
-    LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction;
-    LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction;
-    LocalDirectorFunctionType localDirectorFunction;
-    LocalDirectorReferenceFunctionType localDirectorReferenceFunction;
+    LocalMidSurfaceFunctionType localMidSurfaceFunction(midSurfaceBasis);
+    LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(midSurfaceBasis);
+    LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction(midSurfaceBasis);
+    LocalDirectorFunctionType localDirectorFunction(directorBasis);
+    LocalDirectorReferenceFunctionType localDirectorReferenceFunction(directorBasis);
 
     localMidSurfaceFunction.bind(midSurfaceElement, localMidSurfaceConfiguration);
     localMidSurfaceDisplacementFunction.bind(midSurfaceElement, displacements);
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index 3d49755a77904427a4f0c4bb11d3d31bcd32c2cf..87914ede134c9f08811b4b01c512362f7e8bab92 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -87,8 +87,8 @@ namespace Dune::GFE
 
       typedef LocalGeodesicFEFunction<decltype(deformationScalarBasis), RBM0> LocalGFEFunctionType0;
       typedef LocalGeodesicFEFunction<decltype(rotationScalarBasis), RBM1> LocalGFEFunctionType1;
-      LocalGFEFunctionType0 localGeodesicFEFunction0;
-      LocalGFEFunctionType1 localGeodesicFEFunction1;
+      LocalGFEFunctionType0 localGeodesicFEFunction0(deformationScalarBasis);
+      LocalGFEFunctionType1 localGeodesicFEFunction1(rotationScalarBasis);
       localGeodesicFEFunction0.bind(deformationLocalFiniteElement,localCoefficients[_0]);
       localGeodesicFEFunction1.bind(orientationLocalFiniteElement,localCoefficients[_1]);
 
diff --git a/dune/gfe/functions/discretekirchhoffbendingisometry.hh b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
index 754f50d6ecc64cae430f23c72bbf4746edab790e..0e511e1c4ee88d8d3c608497bbd0a7f806d09372 100644
--- a/dune/gfe/functions/discretekirchhoffbendingisometry.hh
+++ b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
@@ -233,7 +233,7 @@ namespace Dune::GFE
 
       //Create a LocalProjected-Finite element from the local coefficients used for interpolation.
       auto P1LagrangeLFE = localViewCoefficient_.tree().finiteElement();
-      LocalInterpolationRule localPBfunction;
+      LocalInterpolationRule localPBfunction{coefficientBasis_};
       localPBfunction.bind(P1LagrangeLFE,localIsometryCoefficients);
 
       /**
diff --git a/dune/gfe/functions/globalgfefunction.hh b/dune/gfe/functions/globalgfefunction.hh
index 73f3451ea72c7df82e0ee2def92784d5d16a09ca..81fed5c29c3d11bce981b8085fcb8d30f0e008f8 100644
--- a/dune/gfe/functions/globalgfefunction.hh
+++ b/dune/gfe/functions/globalgfefunction.hh
@@ -102,6 +102,7 @@ namespace Dune::GFE
 #endif
           : data_(data)
           , localView_(data_->basis->localView())
+          , localInterpolationRule_(*data_->basis)
         {
           localDoFs_.reserve(localView_.maxSize());
         }
@@ -115,6 +116,7 @@ namespace Dune::GFE
         LocalFunctionBase(const LocalFunctionBase& other)
           : data_(other.data_)
           , localView_(other.localView_)
+          , localInterpolationRule_(other.localInterpolationRule_)
         {
           localDoFs_.reserve(localView_.maxSize());
           if (bound())
diff --git a/dune/gfe/functions/interpolationderivatives.hh b/dune/gfe/functions/interpolationderivatives.hh
index 18caba50b67d99d69c195e0d65edde08d59883f0..74dfff5208b9665adf9fa2139b646f699ceed015 100644
--- a/dune/gfe/functions/interpolationderivatives.hh
+++ b/dune/gfe/functions/interpolationderivatives.hh
@@ -224,7 +224,7 @@ 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;
+      ALocalInterpolationRule localGFEFunction(localInterpolationRule_.globalBasis());
       localGFEFunction.bind(scalarFiniteElement,localAConfiguration);
 
       if (doValue_)
diff --git a/dune/gfe/functions/localgeodesicfefunction.hh b/dune/gfe/functions/localgeodesicfefunction.hh
index 526cb2942a57efc81ca6a6df63544e71d7c62b0b..567394fd8382f0aa0bd62dc986cdaf41243920d6 100644
--- a/dune/gfe/functions/localgeodesicfefunction.hh
+++ b/dune/gfe/functions/localgeodesicfefunction.hh
@@ -60,6 +60,15 @@ 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 with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalGeodesicFEFunction(const Basis basis)
+      : basis_(basis)
+    {}
+
     /** \brief Bind the local function to a particular scalar finite element
      * and a set of coefficients
      *
@@ -94,6 +103,13 @@ namespace Dune::GFE
       return localFiniteElement_.type();
     }
 
+    /** \brief The scalar finite element basis used as interpolation weights
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
     /** \brief The scalar finite element used as the interpolation weights
      *
      * \note This method was added for InterpolationDerivatives, which needs it
@@ -206,6 +222,11 @@ namespace Dune::GFE
       return result;
     }
 
+    /** \brief The scalar basis that implements the weights for the
+     * Riemannian center of mass
+     */
+    const Basis basis_;
+
     /** \brief The scalar local finite element, which provides the weighting factors
         \todo We really only need the local basis
      */
@@ -414,7 +435,7 @@ namespace Dune::GFE
       cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
       cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
 
-      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus, fMinus;
+      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus(basis_), fMinus(basis_);
       fPlus.bind(localFiniteElement_,cornersPlus);
       fMinus.bind(localFiniteElement_,cornersMinus);
 
@@ -554,7 +575,7 @@ namespace Dune::GFE
       cornersPlus[coefficient]  = TargetSpace(aPlus);
       cornersMinus[coefficient] = TargetSpace(aMinus);
 
-      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus,fMinus;
+      LocalGeodesicFEFunction<Basis,TargetSpace> fPlus(basis_), fMinus(basis_);
       fPlus.bind(localFiniteElement_,cornersPlus);
       fMinus.bind(localFiniteElement_,cornersMinus);
 
@@ -623,6 +644,16 @@ 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 with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalGeodesicFEFunction(const Basis basis)
+      : basis_(basis)
+      , orientationFEFunction_(basis)
+    {}
+
     /** \brief Bind the function to a particular weight function set and coefficients
      */
     void bind(const LocalFiniteElement& localFiniteElement,
@@ -835,6 +866,11 @@ namespace Dune::GFE
 
   private:
 
+    /** \brief The scalar basis that implements the weights for the
+     * Riemannian center of mass
+     */
+    const Basis basis_;
+
     /** \brief The scalar local finite element, which provides the weighting factors
         \todo We really only need the local basis
      */
diff --git a/dune/gfe/functions/localprojectedfefunction.hh b/dune/gfe/functions/localprojectedfefunction.hh
index a4e6a49f70988689adc15c43481e73a732570139..4fa37d6a9266e3b7a6ec58334223896839ac6c4d 100644
--- a/dune/gfe/functions/localprojectedfefunction.hh
+++ b/dune/gfe/functions/localprojectedfefunction.hh
@@ -50,6 +50,15 @@ namespace Dune::GFE
     /** \brief The type used for derivatives */
     typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
 
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalProjectedFEFunction(const Basis& basis)
+      : basis_(basis)
+    {}
+
     /** \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
@@ -83,6 +92,13 @@ namespace Dune::GFE
       return localFiniteElement_.type();
     }
 
+    /** \brief The finite element basis used to interpolate in the embedding space
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
     /** \brief The scalar finite element used as the interpolation weights
      *
      * \note This method was added for InterpolationDerivatives, which needs it
@@ -129,6 +145,10 @@ namespace Dune::GFE
     }
   private:
 
+    /** \brief The finite element basis for the interpolation in the embedding space
+     */
+    const Basis basis_;
+
     /** \brief The scalar local finite element, which provides the weighting factors
      *        \todo We really only need the local basis
      */
@@ -360,6 +380,15 @@ namespace Dune::GFE
     /** \brief The type used for derivatives */
     typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
 
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalProjectedFEFunction(const Basis basis)
+      : basis_(basis)
+    {}
+
     /** \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
@@ -393,6 +422,13 @@ namespace Dune::GFE
       return localFiniteElement_.type();
     }
 
+    /** \brief The finite element basis used to interpolate in the embedding space
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
     /** \brief The scalar finite element used as the interpolation weights
      *
      * \note This method was added for InterpolationDerivatives, which needs it
@@ -543,6 +579,10 @@ namespace Dune::GFE
     }
   private:
 
+    /** \brief The finite element basis for the interpolation in the embedding space
+     */
+    const Basis basis_;
+
     /** \brief The scalar local finite element, which provides the weighting factors
      *        \todo We really only need the local basis
      */
@@ -587,6 +627,16 @@ namespace Dune::GFE
     /** \brief The type used for derivatives */
     typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
 
+    /** \brief Constructor with a given scalar basis
+     *
+     * \param basis The basis that implements the weights for the Riemannian
+     * center of mass
+     */
+    LocalProjectedFEFunction(const Basis basis)
+      : basis_(basis)
+      , orientationFunction_(basis)
+    {}
+
     /** \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
@@ -631,6 +681,13 @@ namespace Dune::GFE
       return localFiniteElement_.type();
     }
 
+    /** \brief The finite element basis used to interpolate in the embedding space
+     */
+    const Basis& globalBasis() const
+    {
+      return basis_;
+    }
+
     /** \brief Evaluate the function */
     TargetSpace evaluate(const LocalCoordinate& local) const
     {
@@ -697,6 +754,10 @@ namespace Dune::GFE
     }
   private:
 
+    /** \brief The finite element basis for the interpolation in the embedding space
+     */
+    const Basis basis_;
+
     /** \brief The scalar local finite element, which provides the weighting factors
      *        \todo We really only need the local basis
      */
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 3c2a8e59a584955516e02fb56f274eca4ada0e21..43ad83258078457bb9b50c7460e0e55195c65d29 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -576,7 +576,9 @@ 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>();
+    auto localGFEFunction
+      = std::make_shared<AInterpolationRule>(GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >(deformationFEBasis),
+                                             GFE::LocalGeodesicFEFunction<OrientationFEBasis, GFE::Rotation<adouble,3> >(orientationFEBasis));
 
     using ATargetSpace = typename TargetSpace::rebind<adouble>::other;
 
diff --git a/src/cosserat-rod.cc b/src/cosserat-rod.cc
index 646dc9427ed9506aa082705283a6b754922e5100..a5de9eade86c7e0d709b7e113b4496ed43adfef6 100644
--- a/src/cosserat-rod.cc
+++ b/src/cosserat-rod.cc
@@ -222,14 +222,14 @@ int main (int argc, char *argv[]) try
   // GFE function describing the reference configuration
   using ReferenceGeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<ScalarBasis, TargetSpace>;
   using ReferenceProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, TargetSpace>;
-  ReferenceGeodesicInterpolationRule referenceGeodesicFEFunction;
-  ReferenceProjectedInterpolationRule referenceProjectedFEFunction;
+  ReferenceGeodesicInterpolationRule referenceGeodesicFEFunction(scalarBasis);
+  ReferenceProjectedInterpolationRule referenceProjectedFEFunction(scalarBasis);
 
   using ATargetSpace = TargetSpace::rebind<adouble>::other;
   using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<ScalarBasis, ATargetSpace>;
   using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, ATargetSpace>;
-  GeodesicInterpolationRule localGeodesicFEFunction;
-  ProjectedInterpolationRule localProjectedFEFunction;
+  GeodesicInterpolationRule localGeodesicFEFunction(scalarBasis);
+  ProjectedInterpolationRule localProjectedFEFunction(scalarBasis);
 
   // Assembler using ADOL-C
   std::shared_ptr<GFE::LocalEnergy<ScalarBasis,ATargetSpace> > localRodEnergy;
diff --git a/src/film-on-substrate-stress-plot.cc b/src/film-on-substrate-stress-plot.cc
index 78b3c21ed07c9f20ca252499d58b2c8757c26e8b..0e4d5a8360bb3723626c9b811cff03e76c76f1a2 100644
--- a/src/film-on-substrate-stress-plot.cc
+++ b/src/film-on-substrate-stress-plot.cc
@@ -258,7 +258,7 @@ int main (int argc, char *argv[]) try
 
   using LocalGFEFunctionR = GFE::LocalGeodesicFEFunction<decltype(scalarBasisR),GFE::Rotation<double,dim> >;
 
-  LocalGFEFunctionR localGFEFunction;
+  LocalGFEFunctionR localGFEFunction(scalarBasisR);
 
   auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), LocalGFEFunctionR, FieldVector<double,dim>, GFE::Rotation<double,dim> >
                            (basisOrderD, basisOrderR);
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index f1526a5a85a619c19eb1c4254e4a8a22bf934ed9..ce4acc8f3a5017e66a23b20c984799e5ee059661 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -527,7 +527,8 @@ 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;
+    LocalInterpolationRule localGFEFunction{LocalDeformationInterpolationRule(deformationFEBasis),
+                                            LocalOrientationInterpolationRule(orientationFEBasis)};
 
     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);
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index c98c6a4c6ac195d86d7e02014cc44e4e45b5ca0b..1403b3646e390b721286f415133333e1509cf10e 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -210,7 +210,7 @@ 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;
+  GeodesicInterpolationRule localGFEFunctionA(feBasis);
 
   auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity, ATargetSpace> >();
   addends[0] = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(std::move(localGFEFunctionA), harmonicDensity);
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 18b387fad6dfff3228c5c80ac02ba955f901780c..d17f1fc514e022ec06379032f47fe349bfd036c2 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -290,8 +290,8 @@ 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;
+  GeodesicInterpolationRule localGeodesicFEFunctionA(feBasis);
+  ProjectedInterpolationRule localProjectedFEFunctionA(feBasis);
 
   std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
 
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index d9d4096d930dbc44e97fb8d7071158e79fc96ac5..7901b09c4179713e134a060af2ca9c8d829623e6 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -209,7 +209,8 @@ int main (int argc, char *argv[])
   using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
   using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
-  LocalInterpolationRule localGFEFunctionA;
+  LocalInterpolationRule localGFEFunctionA{LocalDeformationInterpolationRule(deformationFEBasis),
+                                           LocalOrientationInterpolationRule(orientationFEBasis)};
 
   using Element = GridView::Codim<0>::Entity;
   auto bulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<Element,adouble> >(parameters);
diff --git a/test/cosseratrodenergytest.cc b/test/cosseratrodenergytest.cc
index 0ef37b2aed8875a28246f6010a7d117b793ea2ad..a288fc3538179913790919fbb974b2786ed9663e 100644
--- a/test/cosseratrodenergytest.cc
+++ b/test/cosseratrodenergytest.cc
@@ -73,8 +73,8 @@ int main (int argc, char *argv[]) try
 
   using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis,
       TargetSpace>;
-  GeodesicInterpolationRule localGFEFunction;
-  GeodesicInterpolationRule referenceGFEFunction;
+  GeodesicInterpolationRule localGFEFunction(feBasis);
+  GeodesicInterpolationRule referenceGFEFunction(feBasis);
 
   GFE::CosseratRodEnergy<FEBasis,
       GeodesicInterpolationRule,
diff --git a/test/cosseratrodtest.cc b/test/cosseratrodtest.cc
index 630b8e536fb0fc1400eef24dd478c72c86fc1a36..167152916b7b55e53de0473bfe16d3e1c9b938a6 100644
--- a/test/cosseratrodtest.cc
+++ b/test/cosseratrodtest.cc
@@ -134,13 +134,13 @@ int main (int argc, char *argv[]) try
   using ATargetSpace = TargetSpace::rebind<adouble>::other;
   using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<ScalarBasis, ATargetSpace>;
   using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, ATargetSpace>;
-  GeodesicInterpolationRule localGeodesicFEFunction;
-  ProjectedInterpolationRule localProjectedFEFunction;
+  GeodesicInterpolationRule localGeodesicFEFunction(scalarBasis);
+  ProjectedInterpolationRule localProjectedFEFunction(scalarBasis);
 
   using ReferenceGeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<ScalarBasis, TargetSpace>;
   using ReferenceProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, TargetSpace>;
-  ReferenceGeodesicInterpolationRule localReferenceGeodesicFEFunction;
-  ReferenceProjectedInterpolationRule localReferenceProjectedFEFunction;
+  ReferenceGeodesicInterpolationRule localReferenceGeodesicFEFunction(scalarBasis);
+  ReferenceProjectedInterpolationRule localReferenceProjectedFEFunction(scalarBasis);
 
   // Assembler using ADOL-C
   std::shared_ptr<GFE::LocalEnergy<ScalarBasis,ATargetSpace> > localRodEnergy;
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index a7e8d902007f1964e96a2e6b451ac0f34fd909ae..6dc814afce42643aeaa6edba3c8a6be1291d9279 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -353,7 +353,12 @@ int main (int argc, char *argv[])
   using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
   using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
-  LocalInterpolationRule localGFEFunction;
+
+  LocalDeformationInterpolationRule localDeformationGFEFunction(deformationFEBasis);
+  LocalOrientationInterpolationRule localOrientationGFEFunction(orientationFEBasis);
+
+  auto localGFEFunction = std::make_tuple(localDeformationGFEFunction,
+                                          localOrientationGFEFunction);
 
   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);
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
index 60b3025e64d34d626a74f6b98ca0283d9031275e..162c0d1471f4eb0e7f5e3a0153d746d16d68ee07 100644
--- a/test/harmonicmaptest.cc
+++ b/test/harmonicmaptest.cc
@@ -158,7 +158,7 @@ int main (int argc, char *argv[])
   using InterpolationRule = GFE::LocalProjectedFEFunction<FEBasis, TargetSpace,false>;
 #endif
 #endif
-  InterpolationRule localGFEFunction;
+  InterpolationRule localGFEFunction(feBasis);
 
   auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity, TargetSpace> >();
 
diff --git a/test/interpolationderivativestest.cc b/test/interpolationderivativestest.cc
index 480d26cc92038a4623e67b7460cdb5990b767ce6..ee448ca5fbc254a91589efc6649408d4213d0886 100644
--- a/test/interpolationderivativestest.cc
+++ b/test/interpolationderivativestest.cc
@@ -143,7 +143,8 @@ public:
         cornersPlus [coefficient] = TargetSpace(coefficients_[coefficient].globalCoordinates() + variation);
         cornersMinus[coefficient] = TargetSpace(coefficients_[coefficient].globalCoordinates() - variation);
 
-        LocalInterpolationRule fPlus, fMinus;
+        LocalInterpolationRule fPlus(localInterpolationRule_.globalBasis());
+        LocalInterpolationRule fMinus(localInterpolationRule_.globalBasis());
         fPlus.bind(localInterpolationRule_.localFiniteElement(),cornersPlus);
         fMinus.bind(localInterpolationRule_.localFiniteElement(),cornersMinus);
 
@@ -196,7 +197,8 @@ public:
         cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
         cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
 
-        LocalInterpolationRule fPlus, fMinus;
+        LocalInterpolationRule fPlus(localInterpolationRule_.globalBasis());
+        LocalInterpolationRule fMinus(localInterpolationRule_.globalBasis());
         fPlus.bind(localInterpolationRule_.localFiniteElement(),cornersPlus);
         fMinus.bind(localInterpolationRule_.localFiniteElement(),cornersMinus);
 
@@ -256,7 +258,8 @@ public:
         forwardSolution[i]  = TargetSpace::exp(coefficients_[i], eps * xi);
         backwardSolution[i] = TargetSpace::exp(coefficients_[i], -1 * eps * xi);
 
-        LocalInterpolationRule fPlus, fMinus;
+        LocalInterpolationRule fPlus(localInterpolationRule_.globalBasis());
+        LocalInterpolationRule fMinus(localInterpolationRule_.globalBasis());
         fPlus.bind(localInterpolationRule_.localFiniteElement(),forwardSolution);
         fMinus.bind(localInterpolationRule_.localFiniteElement(),backwardSolution);
 
@@ -315,7 +318,8 @@ public:
               backwardSolutionXiEta[j] = TargetSpace::exp(coefficients_[j], (-1)*epsEta);
             }
 
-            LocalInterpolationRule fPlus, fMinus;
+            LocalInterpolationRule fPlus(localInterpolationRule_.globalBasis());
+            LocalInterpolationRule fMinus(localInterpolationRule_.globalBasis());
             fPlus.bind(localInterpolationRule_.localFiniteElement(),forwardSolutionXiEta);
             fMinus.bind(localInterpolationRule_.localFiniteElement(),backwardSolutionXiEta);
 
@@ -422,7 +426,7 @@ TestSuite checkDerivatives()
 
   auto localView = scalarBasis.localView();
   localView.bind(*gridView.begin<0>());
-  LocalInterpolationRule localGFEFunction;
+  LocalInterpolationRule localGFEFunction(scalarBasis);
   localGFEFunction.bind(localView.tree().finiteElement(),localCoefficients);
 
   GFE::InterpolationDerivatives<LocalInterpolationRule> interpolationDerivatives(localGFEFunction,
diff --git a/test/localadolcstiffnesstest.cc b/test/localadolcstiffnesstest.cc
index cb63d4bea3008f3a2c40190e71f7c14e293524a1..75f5354488b9bd8fcb947d1042fd3e6f0ac25758 100644
--- a/test/localadolcstiffnesstest.cc
+++ b/test/localadolcstiffnesstest.cc
@@ -188,7 +188,7 @@ int main (int argc, char *argv[]) try
 
   // Select geometric finite element interpolation method
   using AInterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis, ATargetSpace>;
-  AInterpolationRule localGFEFunctionA;
+  AInterpolationRule localGFEFunctionA(feBasis);
 
   auto activeDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(materialParameters);
 
@@ -204,7 +204,7 @@ int main (int argc, char *argv[]) try
 
   // Select geometric finite element interpolation method
   using InterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis, TargetSpace>;
-  InterpolationRule localGFEFunction;
+  InterpolationRule localGFEFunction(feBasis);
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, double> >(materialParameters);
 
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index 30f1ce3926d2e7a7e721160829e29845afe9509a..f2d0bfb5ebdde52167c649a0dd8329db5d9bede7 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -324,7 +324,7 @@ void test(const GeometryType& element)
     localBasisView.bind(*gridView.template begin<0>());
     const auto& localFiniteElement = localBasisView.tree().finiteElement();
 
-    GFE::LocalGeodesicFEFunction<InterpolationBasis,TargetSpace> f;
+    GFE::LocalGeodesicFEFunction<InterpolationBasis,TargetSpace> f(interpolationBasis);
     f.bind(localFiniteElement,corners);
 
     //testPermutationInvariance(corners);
diff --git a/test/localintegralstiffnesstest.cc b/test/localintegralstiffnesstest.cc
index 7a4e8c7cf0e57ef0ae9a9fd6c8372b7213f2a636..2cad70ba316005f663caea76c103508bfaa4887a 100644
--- a/test/localintegralstiffnesstest.cc
+++ b/test/localintegralstiffnesstest.cc
@@ -109,11 +109,11 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
     std::cout << "Using geodesic interpolation" << std::endl;
     using LocalInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(feBasis),
         TargetSpace>;
-    LocalInterpolationRule localGFEFunction;
+    LocalInterpolationRule localGFEFunction(feBasis);
 
     using ALocalInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(feBasis),
         ATargetSpace>;
-    ALocalInterpolationRule localGFEFunctionA;
+    ALocalInterpolationRule localGFEFunctionA(feBasis);
 
     // Assemble using the old assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(std::move(localGFEFunctionA), harmonicDensityA);
@@ -133,11 +133,11 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
 
     using LocalInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis,
         TargetSpace, interpolationType!=Nonconforming>;
-    LocalInterpolationRule localGFEFunction;
+    LocalInterpolationRule localGFEFunction(feBasis);
 
     using ALocalInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis,
         ATargetSpace, interpolationType!=Nonconforming>;
-    ALocalInterpolationRule localGFEFunctionA;
+    ALocalInterpolationRule localGFEFunctionA(feBasis);
 
     // Assemble using the old assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(std::move(localGFEFunctionA), harmonicDensityA);
@@ -321,13 +321,15 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<double,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
-    LocalInterpolationRule localGFEFunction;
+    LocalInterpolationRule localGFEFunction{LocalDeformationInterpolationRule(deformationFEBasis),
+                                            LocalOrientationInterpolationRule(orientationFEBasis)};
 
     using ALocalDeformationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(deformationFEBasis), GFE::RealTuple<adouble,dim> >;
     using ALocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
-    ALocalInterpolationRule localGFEFunctionA;
+    ALocalInterpolationRule localGFEFunctionA{ALocalDeformationInterpolationRule(deformationFEBasis),
+                                              ALocalOrientationInterpolationRule(orientationFEBasis)};
 
     // Assemble using the ADOL-C assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(std::move(localGFEFunctionA), aBulkCosseratDensity);
@@ -347,13 +349,15 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     using LocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(orientationFEBasis), GFE::Rotation<double,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
-    LocalInterpolationRule localGFEFunction;
+    LocalInterpolationRule localGFEFunction{LocalDeformationInterpolationRule(deformationFEBasis),
+                                            LocalOrientationInterpolationRule(orientationFEBasis)};
 
     using ALocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(deformationFEBasis), GFE::RealTuple<adouble,dim> >;
     using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
-    ALocalInterpolationRule localGFEFunctionA;
+    ALocalInterpolationRule localGFEFunctionA{ALocalDeformationInterpolationRule(deformationFEBasis),
+                                              ALocalOrientationInterpolationRule(orientationFEBasis)};
 
     // Assemble using the ADOL-C assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(std::move(localGFEFunctionA), aBulkCosseratDensity);
diff --git a/test/localprojectedfefunctiontest.cc b/test/localprojectedfefunctiontest.cc
index 56143cf537bf2100a76a830bee3d2df79cacd683..de421777a4a108a26b7315b53f367759773a5c9a 100644
--- a/test/localprojectedfefunctiontest.cc
+++ b/test/localprojectedfefunctiontest.cc
@@ -281,9 +281,9 @@ void test(const GeometryType& element)
     localBasisView.bind(*gridView.template begin<0>());
     const auto& localFiniteElement = localBasisView.tree().finiteElement();
 
-    GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f;
+    GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f(interpolationBasis);
     f.bind(localFiniteElement,corners);
-    GFE::LocalProjectedFEFunction<InterpolationBasis, TargetSpace,false> f_nonconforming;
+    GFE::LocalProjectedFEFunction<InterpolationBasis, TargetSpace,false> f_nonconforming(interpolationBasis);
     f_nonconforming.bind(localFiniteElement,corners);
 
     //testPermutationInvariance(corners);
diff --git a/test/mixedriemannianpnsolvertest.cc b/test/mixedriemannianpnsolvertest.cc
index 0ddaa68dec41e2a9db7a23ca5c6523dee11d4ebc..05c4a6114da19aeb8a0c45a514a8fe046fbc4a0c 100644
--- a/test/mixedriemannianpnsolvertest.cc
+++ b/test/mixedriemannianpnsolvertest.cc
@@ -100,6 +100,7 @@ int main (int argc, char *argv[])
   using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
   DeformationFEBasis deformationFEBasis(gridView);
   using RotationFEBasis = Functions::LagrangeBasis<GridView,rotationOrder>;
+  RotationFEBasis rotationFEBasis(gridView);
 
   /////////////////////////////////////////////////////////////////////////
   //  Create the Neumann and Dirichlet boundary
@@ -174,7 +175,8 @@ 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;
+  AInterpolationRule localGFEFunctionA{GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >(deformationFEBasis),
+                                       GFE::LocalGeodesicFEFunction<RotationFEBasis, GFE::Rotation<adouble,3> >(rotationFEBasis)};
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
 
diff --git a/test/planarcosseratshelltest.cc b/test/planarcosseratshelltest.cc
index 0db975787afe8656d89f6ba98561ddc85358c503..5b0543f7baf60368a22aa001d2be829f4f456465 100644
--- a/test/planarcosseratshelltest.cc
+++ b/test/planarcosseratshelltest.cc
@@ -201,7 +201,8 @@ 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;
+  AInterpolationRule localGFEFunctionA{GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >(deformationFEBasis),
+                                       GFE::LocalGeodesicFEFunction<OrientationFEBasis, GFE::Rotation<adouble,3> >(orientationFEBasis)};
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
 
diff --git a/test/surfacecosseratstressassemblertest.cc b/test/surfacecosseratstressassemblertest.cc
index 508146eb2c6f5b75dc6e51970b449736cf166389..9f61964925acf1aba4cbee05cbaf5f66daac3aff 100644
--- a/test/surfacecosseratstressassemblertest.cc
+++ b/test/surfacecosseratstressassemblertest.cc
@@ -243,7 +243,7 @@ int main (int argc, char *argv[])
 
   using LocalGFEFunctionR = GFE::LocalGeodesicFEFunction<decltype(scalarBasisR),GFE::Rotation<double,dim> >;
 
-  LocalGFEFunctionR localGFEFunction;
+  LocalGFEFunctionR localGFEFunction(scalarBasisR);
 
   auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),
       decltype(basisOrderR),