diff --git a/dune/gfe/assemblers/chiralskyrmionenergy.hh b/dune/gfe/assemblers/chiralskyrmionenergy.hh
index d3e02a67495370f8439e14985ea167f2ccb905a0..d88ef7e3fd1f3a388a893c626d5ffc9726934a10 100644
--- a/dune/gfe/assemblers/chiralskyrmionenergy.hh
+++ b/dune/gfe/assemblers/chiralskyrmionenergy.hh
@@ -44,6 +44,12 @@ namespace Dune {
       RT energy (const typename Basis::LocalView& localView,
                  const std::vector<TargetSpace>& localConfiguration) const override;
 
+      RT energy (const typename Basis::LocalView& localView,
+                 const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+      {
+        DUNE_THROW(NotImplemented, "!");
+      }
+
       field_type h_;
       field_type kappa_;
     };
diff --git a/dune/gfe/assemblers/cosseratenergystiffness.hh b/dune/gfe/assemblers/cosseratenergystiffness.hh
index 2bbac9fb0315842dab4755a470ce1fc87b993e4f..9153bef823830a37936e899a394a169918b488c7 100644
--- a/dune/gfe/assemblers/cosseratenergystiffness.hh
+++ b/dune/gfe/assemblers/cosseratenergystiffness.hh
@@ -143,6 +143,12 @@ public:
   RT energy (const typename Basis::LocalView& localView,
              const std::vector<TargetSpace>& localSolution) const override;
 
+  RT energy (const typename Basis::LocalView& localView,
+             const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
+
   /** \brief Assemble the energy for a single element */
   RT energy (const typename Basis::LocalView& localView,
              const std::vector<RealTuple<field_type,dim> >& localDisplacementConfiguration,
diff --git a/dune/gfe/assemblers/cosseratrodenergy.hh b/dune/gfe/assemblers/cosseratrodenergy.hh
index 5a3361354d1429d7bdfa5c23aea02443c07682a6..27bc7e6581e01720db259a61682a26b5037da4ec 100644
--- a/dune/gfe/assemblers/cosseratrodenergy.hh
+++ b/dune/gfe/assemblers/cosseratrodenergy.hh
@@ -109,6 +109,12 @@ namespace Dune::GFE {
     virtual RT energy (const typename Basis::LocalView& localView,
                        const std::vector<TargetSpace>& localSolution) const override;
 
+    virtual RT energy (const typename Basis::LocalView& localView,
+                       const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
     /** \brief Get the rod strain at one point in the rod
      *
      * \tparam Number This is a member template because the method has to work for double and adouble
diff --git a/dune/gfe/assemblers/harmonicenergy.hh b/dune/gfe/assemblers/harmonicenergy.hh
index 7e8a7139aa8150f78c13ba2fe737cd8ec384f3d7..90c47a93ca7627f286ace775b9639bae253ef40e 100644
--- a/dune/gfe/assemblers/harmonicenergy.hh
+++ b/dune/gfe/assemblers/harmonicenergy.hh
@@ -26,6 +26,12 @@ public:
   RT energy (const typename Basis::LocalView& localView,
              const std::vector<TargetSpace>& localSolution) const override;
 
+  virtual RT energy (const typename Basis::LocalView& localView,
+                     const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
+
 };
 
 template <class Basis, class LocalInterpolationRule, class TargetSpace>
diff --git a/dune/gfe/assemblers/l2distancesquaredenergy.hh b/dune/gfe/assemblers/l2distancesquaredenergy.hh
index 1fe1ff611133b7866b113a92eb4b0ed59d7f4294..cffbb01af51704bccce10e6b1a185307c4d46eea 100644
--- a/dune/gfe/assemblers/l2distancesquaredenergy.hh
+++ b/dune/gfe/assemblers/l2distancesquaredenergy.hh
@@ -80,6 +80,11 @@ public:
     return energy;
   }
 
+  virtual RT energy (const typename Basis::LocalView& localView,
+                     const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
 };
 
 #endif
diff --git a/dune/gfe/assemblers/localenergy.hh b/dune/gfe/assemblers/localenergy.hh
index 58c122c5f5e2043e458892e98082c95bff98b9d5..e188993fbead5ee7fc97c0a715f926d3a815550d 100644
--- a/dune/gfe/assemblers/localenergy.hh
+++ b/dune/gfe/assemblers/localenergy.hh
@@ -3,16 +3,43 @@
 
 #include <vector>
 
+#include <dune/common/tuplevector.hh>
+
+#include <dune/gfe/spaces/productmanifold.hh>
+
 namespace Dune {
 
+  namespace GFE:: Impl
+  {
+    /** \brief A class exporting container types for coefficient sets
+     *
+     * This generic template handles TargetSpaces that are not product manifolds.
+     */
+    template <class TargetSpace>
+    struct LocalEnergyTypes
+    {
+      using Coefficients = std::vector<TargetSpace>;
+      using CompositeCoefficients = TupleVector<std::vector<TargetSpace> >;
+    };
+
+    /** \brief A class exporting container types for coefficient sets -- specialization for product manifolds
+     */
+    template <class ... Factors>
+    struct LocalEnergyTypes<ProductManifold<Factors...> >
+    {
+      using Coefficients = std::vector<ProductManifold<Factors...> >;
+      using CompositeCoefficients = TupleVector<std::vector<Factors>... >;
+    };
+  }
+
   namespace GFE {
 
     /** \brief Base class for energies defined by integrating over one grid element */
-    template<class Basis, class ... TargetSpaces>
+    template<class Basis, class TargetSpace>
     class LocalEnergy
     {
     public:
-      using RT = typename std::common_type<typename TargetSpaces::ctype...>::type;
+      using RT = typename TargetSpace::ctype;
 
       /** \brief Compute the energy
        *
@@ -21,7 +48,14 @@ namespace Dune {
        */
       virtual RT
       energy (const typename Basis::LocalView& localView,
-              const std::vector<TargetSpaces>& ... localSolution) const = 0;
+              const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
+
+      /** \brief ProductManifolds: Compute the energy from coefficients in separate containers
+       * for each factor
+       */
+      virtual RT
+      energy (const typename Basis::LocalView& localView,
+              const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
 
       /** Empty virtual default destructor
        *
diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index 74565e23058cf5d090cfba195e2334118e05253a..fb59c5cd9437a4a758ea331daf3578031e08c929 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -49,6 +49,15 @@ public:
   virtual RT energy (const typename Basis::LocalView& localView,
                      const std::vector<TargetSpace>& localSolution) const override;
 
+  /** \brief ProductManifolds: Compute the energy from coefficients in separate containers
+   * for each factor
+   */
+  virtual RT energy (const typename Basis::LocalView& localView,
+                     const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
+
   /** \brief Assemble the element gradient of the energy functional
 
      This uses the automatic differentiation toolbox ADOL_C.
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index 442fc88bc4a7f475bcc592628909e788b5e12d5a..cf400921d0c5636a164ea64f85533cabf78ba8c9 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -42,6 +42,12 @@ public:
     return localEnergy_->energy(localView,localSolution);
   }
 
+  RT energy (const typename Basis::LocalView& localView,
+             const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
+
   /** \brief Assemble the element gradient of the energy functional
 
      The default implementation in this class uses a finite difference approximation */
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index f2de00369bfa3f26906baf33a37c06a4b92755ac..50ff86fcc7385f13e29cb49e0bb1dd093fb2c1b6 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -29,12 +29,13 @@ namespace Dune::GFE {
    */
   template<class Basis, class ... TargetSpaces>
   class LocalIntegralEnergy
-    : public Dune::GFE::LocalEnergy<Basis,TargetSpaces...>
+    : public Dune::GFE::LocalEnergy<Basis,ProductManifold<TargetSpaces...> >
   {
+    using TargetSpace = ProductManifold<TargetSpaces...>;
     using LocalView = typename Basis::LocalView;
     using GridView = typename LocalView::GridView;
     using DT = typename GridView::Grid::ctype;
-    using RT = typename GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+    using RT = typename GFE::LocalEnergy<Basis,TargetSpace>::RT;
 
     constexpr static int gridDim = GridView::dimension;
 
@@ -60,15 +61,20 @@ namespace Dune::GFE {
 
     /** \brief Assemble the energy for a single element */
     RT energy(const typename Basis::LocalView& localView,
-              const std::vector<TargetSpaces>& ... localSolutions) const
+              const std::vector<TargetSpace>& localSolutions) const
     {
-      const auto& element = localView.element();
-
+      DUNE_THROW(NotImplemented, "!");
+    }
 
-      const std::vector<TargetSpaceDeformation>& localDeformationConfiguration = std::get<0>(std::forward_as_tuple(localSolutions ...));
-      const std::vector<TargetSpaceRotation>& localOrientationConfiguration = std::get<1>(std::forward_as_tuple(localSolutions ...));
+    RT energy (const typename Basis::LocalView& localView,
+               const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      const auto& element = localView.element();
 
       using namespace Indices;
+      const std::vector<TargetSpaceDeformation>& localDeformationConfiguration = coefficients[_0];
+      const std::vector<TargetSpaceRotation>& localOrientationConfiguration = coefficients[_1];
+
       // composite Basis: grab the finite element of the first child
       const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
       const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index 54b8a1b5cea33a274cf731cc722eec1e89c3060e..0140348d8a90bef221d61820cb30bd81b8030c50 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -95,6 +95,12 @@ public:
   RT energy (const typename Basis::LocalView& localView,
              const std::vector<TargetSpace>& localSolution) const override;
 
+  RT energy (const typename Basis::LocalView& localView,
+             const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
+
   /** \brief Assemble the energy for a single element */
   RT energy (const typename Basis::LocalView& localView,
              const std::vector<RealTuple<field_type,dim> >& localDisplacementConfiguration,
diff --git a/dune/gfe/assemblers/simofoxenergy.hh b/dune/gfe/assemblers/simofoxenergy.hh
index dddf5eb897807911557b145dc16ea8982f60a196..4102c08c56f68137c16bad826bb0eac089085d39 100644
--- a/dune/gfe/assemblers/simofoxenergy.hh
+++ b/dune/gfe/assemblers/simofoxenergy.hh
@@ -52,11 +52,13 @@ namespace Dune::GFE {
    */
   template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type = double>
   class SimoFoxEnergyLocalStiffness
-    : public Dune::GFE::LocalEnergy<Basis, RealTuple<field_type, 3>,
-          UnitVector<field_type, 3> >,                             // inheritance to allow usage with LocalGeodesicFEADOLCStiffness
+    : public Dune::GFE::LocalEnergy<Basis, ProductManifold<RealTuple<field_type, 3>,
+          UnitVector<field_type, 3> > >,                             // inheritance to allow usage with LocalGeodesicFEADOLCStiffness
       public MixedLocalGeodesicFEStiffness<Basis, ProductManifold<RealTuple<field_type, 3>,
           UnitVector<field_type, 3> > >                                    // inheritance to allow usage with MixedGFEAssembler
   {
+    using TargetSpace = ProductManifold<RealTuple<field_type, 3>, UnitVector<field_type, 3> >;
+
     // grid types
     typedef typename Basis::GridView GridView;
     typedef typename GridView::ctype DT;
@@ -122,6 +124,19 @@ namespace Dune::GFE {
       CMat_[6][6] = CMat_[7][7] = fac3;
     }
 
+    /** \brief Assemble the energy for a single element */
+    RT energy(const typename Basis::LocalView &localView,
+              const std::vector<TargetSpace> &localConfiguration) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
+    RT energy (const typename Basis::LocalView& localView,
+               const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
     /** \brief Assemble the energy for a single element */
     RT energy(const typename Basis::LocalView &localView, const std::vector<RealTuple<field_type, 3> > &localMidSurfaceConfiguration,
               const std::vector<UnitVector<field_type, 3> > &localDirectorConfiguration) const override;
diff --git a/dune/gfe/assemblers/sumenergy.hh b/dune/gfe/assemblers/sumenergy.hh
index 46f84e7b88735f9deb58ada751b857de1ba66c92..312bd4bbfefd0d14eb7156be2d646983a4897c2b 100644
--- a/dune/gfe/assemblers/sumenergy.hh
+++ b/dune/gfe/assemblers/sumenergy.hh
@@ -17,17 +17,18 @@ namespace Dune::GFE {
 
   template<class Basis, class ... TargetSpaces>
   class SumEnergy
-    : public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>,
+    : public Dune::GFE::LocalEnergy<Basis, ProductManifold<TargetSpaces...> >,
       public MixedLocalGeodesicFEStiffness<Basis, ProductManifold<TargetSpaces...> >
       //Inheriting from MixedLocalGeodesicFEStiffness is hack, and will be replaced eventually; once MixedLocalGFEADOLCStiffness
       //will be removed and its functionality will be included in LocalGeodesicFEADOLCStiffness this is not needed anymore!
 
   {
+    using TargetSpace = ProductManifold<TargetSpaces...>;
 
     using LocalView = typename Basis::LocalView;
     using GridView = typename LocalView::GridView;
     using DT = typename GridView::Grid::ctype;
-    using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+    using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpace>::RT;
 
   public:
 
@@ -35,24 +36,44 @@ namespace Dune::GFE {
     SumEnergy()
     {}
 
-    void addLocalEnergy(std::shared_ptr<GFE::LocalEnergy<Basis,TargetSpaces...> > localEnergy)
+    void addLocalEnergy(std::shared_ptr<GFE::LocalEnergy<Basis,TargetSpace> > localEnergy)
     {
       localEnergies_.push_back(localEnergy);
     }
 
     RT energy(const typename Basis::LocalView& localView,
-              const std::vector<TargetSpaces>& ... localSolution) const
+              const std::vector<TargetSpace>& localSolution) const override
     {
       RT sum = 0.;
       for ( const auto& localEnergy : localEnergies_ )
-        sum += localEnergy->energy( localView, localSolution ... );
+        sum += localEnergy->energy(localView, localSolution);
+
+      return sum;
+    }
+
+    RT energy(const typename Basis::LocalView& localView,
+              const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      RT sum = 0.;
+      for ( const auto& localEnergy : localEnergies_ )
+        sum += localEnergy->energy(localView, coefficients);
+
+      return sum;
+    }
+
+    RT energy(const typename Basis::LocalView& localView,
+              const std::vector<TargetSpaces>& ... localSolution) const override
+    {
+      RT sum = 0.;
+      for ( const auto& localEnergy : localEnergies_ )
+        sum += localEnergy->energy(localView, TupleVector<std::vector<TargetSpaces>...>(localSolution ...));
 
       return sum;
     }
 
   private:
     // vector containing pointers to the local energies
-    std::vector<std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpaces...> > > localEnergies_;
+    std::vector<std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace> > > localEnergies_;
   };
 
 }  // namespace Dune::GFE
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index ca634df28e12e4b99105383b904854912a9bae63..1deda7d4b8d0dad93e551fe5c1b3a5dae126e792 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -28,11 +28,12 @@ namespace Dune::GFE {
    */
   template<class CurvedGeometryGridFunction, class Basis, class ... TargetSpaces>
   class SurfaceCosseratEnergy
-    : public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>
+    : public Dune::GFE::LocalEnergy<Basis, ProductManifold<TargetSpaces...> >
   {
+    using TargetSpace = ProductManifold<TargetSpaces...>;
     using GridView = typename Basis::GridView;
     using DT = typename GridView::ctype ;
-    using RT = typename Dune::GFE::LocalEnergy<Basis, TargetSpaces...>::RT ;
+    using RT = typename Dune::GFE::LocalEnergy<Basis, TargetSpace>::RT ;
     using Entity = typename GridView::template Codim<0>::Entity ;
     using RBM0 = RealTuple<RT,GridView::dimensionworld> ;
     using RBM1 = Rotation<RT,GridView::dimensionworld> ;
@@ -105,6 +106,19 @@ namespace Dune::GFE {
       b3_ = parameters.template get<double>("b3");
     }
 
+    /** \brief Assemble the energy for a single element */
+    RT energy(const typename Basis::LocalView& localView,
+              const std::vector<TargetSpace>& localSolutions) const
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
+    RT energy (const typename Basis::LocalView& localView,
+               const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
     RT energy(const typename Basis::LocalView& localView,
               const std::vector<TargetSpaces>& ... localSolutions) const
     {
diff --git a/dune/gfe/assemblers/weightedsumenergy.hh b/dune/gfe/assemblers/weightedsumenergy.hh
index 9f7d7179e8a82349edf4d43281a94b7201e31778..7dfab8366b2aa6857cf68bd4014e417668911c4f 100644
--- a/dune/gfe/assemblers/weightedsumenergy.hh
+++ b/dune/gfe/assemblers/weightedsumenergy.hh
@@ -48,6 +48,11 @@ public:
     return energy;
   }
 
+  RT energy (const typename Basis::LocalView& localView,
+             const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+  {
+    DUNE_THROW(Dune::NotImplemented, "!");
+  }
 };
 
 #endif
diff --git a/dune/gfe/neumannenergy.hh b/dune/gfe/neumannenergy.hh
index d9e0551748defb5cb394434298b0c5f2ec6b34dd..f91c44e4d59072f9249a73883e4077c37a17c878 100644
--- a/dune/gfe/neumannenergy.hh
+++ b/dune/gfe/neumannenergy.hh
@@ -8,20 +8,20 @@
 #include <dune/elasticity/assemblers/localenergy.hh>
 
 namespace Dune::GFE {
-  /**
-     \brief Assembles the Neumann energy for a single element on the Neumann Boundary using the Neumann Function.
-
-           This class works similarly to the class Dune::Elasticity::NeumannEnergy, where Dune::Elasticity::NeumannEnergy extends
-           Dune::Elasticity::LocalEnergy and Dune::GFE::NeumannEnergy extends Dune::GFE::LocalEnergy.
+  /** \brief Integrate a density over a part of the domain boundary
+   *
+   * This is typically used to implement Neumann boundary conditions.  Hence the name.
    */
   template<class Basis, class ... TargetSpaces>
   class NeumannEnergy
-    : public Dune::GFE::LocalEnergy<Basis,TargetSpaces...>
+    : public Dune::GFE::LocalEnergy<Basis,ProductManifold<TargetSpaces...> >
   {
+    using TargetSpace = ProductManifold<TargetSpaces...>;
+
     using LocalView = typename Basis::LocalView;
     using GridView = typename LocalView::GridView;
     using DT = typename GridView::Grid::ctype;
-    using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+    using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpace>::RT;
 
     constexpr static int dim = GridView::dimension;
 
@@ -38,13 +38,18 @@ namespace Dune::GFE {
 
     /** \brief Assemble the energy for a single element */
     RT energy(const typename Basis::LocalView& localView,
-              const std::vector<TargetSpaces>& ... localSolutions) const
+              const std::vector<TargetSpace>& localSolutions) const
     {
-      static_assert(sizeof...(TargetSpaces) > 0, "NeumannEnergy needs at least one TargetSpace!");
+      DUNE_THROW(NotImplemented, "!");
+    }
 
+    RT energy (const typename Basis::LocalView& localView,
+               const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
       using namespace Dune::Indices;
+      // TODO: Remove the hard-coded first factor space!
       using TargetSpace = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
-      const std::vector<TargetSpace>& localSolution = std::get<0>(std::forward_as_tuple(localSolutions ...));
+      const std::vector<TargetSpace>& localSolution = coefficients[_0];
 
       const auto& localFiniteElement = localView.tree().child(_0,0).finiteElement();
       const auto& element = localView.element();