From e794a07f46fa74cb86c4c9bd39f27ab558b634b7 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 7 Feb 2025 10:18:06 +0100
Subject: [PATCH] Let local GFE functions store a basis

This is the basis used as interpolation weights (for geodesic
FE functions) or for interpolating in the embedding space
(for projection-based FE functions).  The basis is not actually
used yet, but it needs to be given to the class in the
constructor.
---
 CHANGELOG.md                                  |  5 +-
 .../gfe/assemblers/l2distancesquaredenergy.hh |  2 +-
 .../nonplanarcosseratshellenergy.hh           |  6 +-
 dune/gfe/assemblers/simofoxenergy.hh          | 10 +--
 dune/gfe/assemblers/surfacecosseratenergy.hh  |  4 +-
 .../discretekirchhoffbendingisometry.hh       |  2 +-
 dune/gfe/functions/globalgfefunction.hh       |  2 +
 .../gfe/functions/interpolationderivatives.hh |  2 +-
 dune/gfe/functions/localgeodesicfefunction.hh | 40 +++++++++++-
 .../gfe/functions/localprojectedfefunction.hh | 61 +++++++++++++++++++
 src/cosserat-continuum.cc                     |  4 +-
 src/cosserat-rod.cc                           |  8 +--
 src/film-on-substrate-stress-plot.cc          |  2 +-
 src/film-on-substrate.cc                      |  3 +-
 src/gradient-flow.cc                          |  2 +-
 src/harmonicmaps.cc                           |  4 +-
 test/cosseratcontinuumtest.cc                 |  3 +-
 test/cosseratrodenergytest.cc                 |  4 +-
 test/cosseratrodtest.cc                       |  8 +--
 test/filmonsubstratetest.cc                   |  7 ++-
 test/harmonicmaptest.cc                       |  2 +-
 test/interpolationderivativestest.cc          | 14 +++--
 test/localadolcstiffnesstest.cc               |  4 +-
 test/localgeodesicfefunctiontest.cc           |  2 +-
 test/localintegralstiffnesstest.cc            | 20 +++---
 test/localprojectedfefunctiontest.cc          |  4 +-
 test/mixedriemannianpnsolvertest.cc           |  4 +-
 test/planarcosseratshelltest.cc               |  3 +-
 test/surfacecosseratstressassemblertest.cc    |  2 +-
 29 files changed, 178 insertions(+), 56 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f452d905..fb052629 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 4d625a47..bb8b4db1 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 f4e30ac6..6127c16e 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 75ff2888..a00cae84 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 3d49755a..87914ede 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 754f50d6..0e511e1c 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 73f3451e..81fed5c2 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 18caba50..74dfff52 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 526cb294..567394fd 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 a4e6a49f..4fa37d6a 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 3c2a8e59..43ad8325 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 646dc942..a5de9ead 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 78b3c21e..0e4d5a83 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 f1526a5a..ce4acc8f 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 c98c6a4c..1403b364 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 18b387fa..d17f1fc5 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 d9d4096d..7901b09c 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 0ef37b2a..a288fc35 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 630b8e53..16715291 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 a7e8d902..6dc814af 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 60b3025e..162c0d14 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 480d26cc..ee448ca5 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 cb63d4be..75f53544 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 30f1ce39..f2d0bfb5 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 7a4e8c7c..2cad70ba 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 56143cf5..de421777 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 0ddaa68d..05c4a611 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 0db97578..5b0543f7 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 508146eb..9f619649 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),
-- 
GitLab