diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index 50ff86fcc7385f13e29cb49e0bb1dd093fb2c1b6..37a29a84889f74c3ad0c96450649133b5144c694 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -21,17 +21,15 @@
 
 namespace Dune::GFE {
 
-  /**
-     \brief Assembles the elastic energy for a single element integrating the localdensity over one element.
-
-           This class works similarly to the class Dune::Elasticity::LocalIntegralEnergy, where Dune::Elasticity::LocalIntegralEnergy extends
-           Dune::Elasticity::LocalEnergy and Dune::GFE::LocalIntegralEnergy extends Dune::GFE::LocalEnergy.
+  /** \brief An energy given as an integral over a density
+   *
+   * \tparam Basis The scalar finite element basis used to construct the interpolation rule
+   * \tparam TargetSpace The space that the geometric finite element function maps into
    */
-  template<class Basis, class ... TargetSpaces>
+  template<class Basis, class TargetSpace>
   class LocalIntegralEnergy
-    : public Dune::GFE::LocalEnergy<Basis,ProductManifold<TargetSpaces...> >
+    : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
   {
-    using TargetSpace = ProductManifold<TargetSpaces...>;
     using LocalView = typename Basis::LocalView;
     using GridView = typename LocalView::GridView;
     using DT = typename GridView::Grid::ctype;
@@ -39,10 +37,6 @@ namespace Dune::GFE {
 
     constexpr static int gridDim = GridView::dimension;
 
-    static_assert(sizeof...(TargetSpaces) == 2, "LocalGeodesicIntegralEnergy needs two TargetSpaces!");
-    using TargetSpaceDeformation = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
-    using TargetSpaceRotation = typename std::tuple_element<1, std::tuple<TargetSpaces...> >::type;
-
   public:
 
     /** \brief Constructor with a Dune::Elasticity::LocalDensity
@@ -53,7 +47,7 @@ namespace Dune::GFE {
 
     /** \brief Constructor with a Dune::GFE::LocalDensity
      */
-    LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,ProductManifold<TargetSpaces...> > >& ld)
+    LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> >& ld)
       : localDensityGFE_(ld)
     {}
 
@@ -69,6 +63,11 @@ namespace Dune::GFE {
     RT energy (const typename Basis::LocalView& localView,
                const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
     {
+      // TODO: Cosserat materials are hard-wired here for historical reasons.
+      static_assert(TargetSpace::size() == 2, "LocalGeodesicIntegralEnergy needs two TargetSpaces!");
+      using TargetSpaceDeformation = typename std::tuple_element<0, TargetSpace>::type;
+      using TargetSpaceRotation = typename std::tuple_element<1, TargetSpace>::type;
+
       const auto& element = localView.element();
 
       using namespace Indices;
@@ -148,7 +147,7 @@ namespace Dune::GFE {
 
   protected:
     const std::shared_ptr<Elasticity::LocalDensity<gridDim,RT,DT> > localDensityElasticity_ = nullptr;
-    const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,ProductManifold<TargetSpaces...> > > localDensityGFE_ = nullptr;
+    const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> > localDensityGFE_ = nullptr;
   };
 
 }  // namespace Dune::GFE
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index da1337cfe0c26d0105a9fc5253f01dc953ef8ea4..51795b97a500a0a70b58fbfa19c8df679292b14e 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -508,7 +508,9 @@ int main (int argc, char *argv[]) try
     if(!elasticDensity)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
-    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(elasticDensity);
+    using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+
+    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ActiveRigidBodyMotion> >(elasticDensity);
     auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,*neumannFunctionPtr);
     auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<
         decltype(stressFreeShellFunction), CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> > >(
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index 90ab7653c90c593ef6a9768e854b0db6ca628c75..2b5ed2e175f7cff37031267fd5d5199d57d9231c 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -198,15 +198,17 @@ int main (int argc, char *argv[])
   //  Create an assembler
   // ////////////////////////////
 
+  using RigidBodyMotion = GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> >;
+  using ARigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>,Rotation<adouble,dim> >;
+
   GFE::SumEnergy<CompositeBasis, RealTuple<adouble,dim>,Rotation<adouble,dim> > sumEnergy;
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
   auto bulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<FieldVector<double,dim>,adouble> >(parameters);
-  auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(bulkCosseratDensity);
+  auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ARigidBodyMotion> >(bulkCosseratDensity);
   sumEnergy.addLocalEnergy(bulkCosseratEnergy);
   sumEnergy.addLocalEnergy(neumannEnergy);
 
-  LocalGeodesicFEADOLCStiffness<CompositeBasis,
-      GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> > > localGFEADOLCStiffness(&sumEnergy);
+  LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(&sumEnergy);
   MixedGFEAssembler<CompositeBasis,
       RealTuple<double,dim>,
       Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index 81116e7ceb9904d7bfa7d06feab2713d20219b14..0c1988588ee8a94a2c41d2ddeda5ca14402e42f1 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -333,7 +333,9 @@ int main (int argc, char *argv[])
   std::shared_ptr<Elasticity::LocalDensity<dim,ValueType> > elasticDensity;
   elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType> >(materialParameters);
 
-  auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(elasticDensity);
+  using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+
+  auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ActiveRigidBodyMotion> >(elasticDensity);
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
 
   auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<