diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index 23c8dd8d2b34656fdb1dd3e07e3694d63f05a46e..f2de00369bfa3f26906baf33a37c06a4b92755ac 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -38,6 +38,10 @@ 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
@@ -48,7 +52,7 @@ namespace Dune::GFE {
 
     /** \brief Constructor with a Dune::GFE::LocalDensity
      */
-    LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<gridDim,RT,DT> >& ld)
+    LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,ProductManifold<TargetSpaces...> > >& ld)
       : localDensityGFE_(ld)
     {}
 
@@ -60,10 +64,6 @@ namespace Dune::GFE {
     {
       const auto& element = localView.element();
 
-      static_assert(sizeof...(TargetSpaces) > 1, "LocalGeodesicIntegralEnergy needs at least two TargetSpace!");
-
-      using TargetSpaceDeformation = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
-      using TargetSpaceRotation = typename std::tuple_element<1, std::tuple<TargetSpaces...> >::type;
 
       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 ...));
@@ -142,7 +142,7 @@ namespace Dune::GFE {
 
   protected:
     const std::shared_ptr<Elasticity::LocalDensity<gridDim,RT,DT> > localDensityElasticity_ = nullptr;
-    const std::shared_ptr<GFE::LocalDensity<gridDim,RT,DT> > localDensityGFE_ = nullptr;
+    const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,ProductManifold<TargetSpaces...> > > localDensityGFE_ = nullptr;
   };
 
 }  // namespace Dune::GFE
diff --git a/dune/gfe/densities/bulkcosseratdensity.hh b/dune/gfe/densities/bulkcosseratdensity.hh
index 9f4cddbda67da91eeef45c656f0b63ab8bd1cd84..5a1ce6f8a08a9dd239d1711d4c11513cefa3c2ba 100644
--- a/dune/gfe/densities/bulkcosseratdensity.hh
+++ b/dune/gfe/densities/bulkcosseratdensity.hh
@@ -9,16 +9,18 @@
 #include <dune/gfe/cosseratstrain.hh>
 #include <dune/gfe/linearalgebra.hh>
 #include <dune/gfe/densities/localdensity.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
 
 namespace Dune::GFE {
 
-  template<class field_type = double, class ctype = double>
+  template<class Position, class field_type>
   class BulkCosseratDensity final
-    : public GFE::LocalDensity<3,field_type,ctype>
+    : public GFE::LocalDensity<Position, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
   {
   private:
-    // gridDim = 3 - this can be hardwired here, because BulkCosseratDensity only works for 3d->3d problems
-    static const int gridDim = 3;
+    // BulkCosseratDensity only works for 3d->3d problems
+    static_assert(Position::size()==3);
+    static const int gridDim = Position::size();
     static const int embeddedDim = Rotation<field_type,gridDim>::embeddedDim;
     using OrientationDerivativeType = FieldMatrix<field_type, embeddedDim, gridDim>;
 
@@ -116,7 +118,7 @@ namespace Dune::GFE {
      * \param orientationValue The orientation at the current position
      * \param orientationDerivative The derivative of the orientation at the current position
      */
-    field_type operator() (const FieldVector<ctype,gridDim>& x,
+    field_type operator() (const Position& x,
                            const RealTuple<field_type,gridDim>& deformationValue,
                            const FieldMatrix<field_type,gridDim,gridDim>& deformationDerivative,
                            const Rotation<field_type,gridDim>& orientationValue,
diff --git a/dune/gfe/densities/localdensity.hh b/dune/gfe/densities/localdensity.hh
index 8f29d58df9f483e5d9fd2876cbaed5ee2fe405f8..80df27fca2ec160fae496d1a15d7b7ac3d13d782 100644
--- a/dune/gfe/densities/localdensity.hh
+++ b/dune/gfe/densities/localdensity.hh
@@ -9,14 +9,13 @@ namespace Dune::GFE {
 
   /** \brief A base class for energy densities to be evaluated in an integral energy
    *
-   * \tparam field_type type of the gradient entries
-   * \tparam ctype type of the coordinates
+   * \tparam Position The evaluation point in the integration domain
+   * \tparam TargetSpace Type for the function value
    */
-  template<int dim, class field_type = double, class ctype = double>
+  template<class Position, class TargetSpace>
   class LocalDensity
   {
-  private:
-    static const int embeddedDim = Rotation<field_type,dim>::embeddedDim;
+    using field_type = typename TargetSpace::field_type;
   public:
 
     /** \brief Evaluation with the current position, the deformation function, the deformation gradient, the rotation and the rotation gradient
@@ -27,11 +26,11 @@ namespace Dune::GFE {
      * \param orientationValue The orientation at the current position
      * \param orientationDerivative The derivative of the orientation at the current position
      */
-    virtual field_type operator() (const FieldVector<ctype,dim>& x,
-                                   const RealTuple<field_type,dim>& deformation,
-                                   const FieldMatrix<field_type,dim,dim>& gradient,
-                                   const Rotation<field_type,dim>& rotation,
-                                   const FieldMatrix<field_type, embeddedDim, dim>& rotationGradient) const = 0;
+    virtual field_type operator() (const Position& x,
+                                   const RealTuple<field_type,3>& deformation,
+                                   const FieldMatrix<field_type,3,3>& gradient,
+                                   const Rotation<field_type,3>& rotation,
+                                   const FieldMatrix<field_type, 4, 3>& rotationGradient) const = 0;
 
   };
 
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index 54f8dd910af829852572460eb1ca0448fe2341e7..2680242dde6e28d3ffd1fd2bb17aca510ca7980a 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -202,7 +202,7 @@ int main (int argc, char *argv[])
 
   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<adouble,double> >(parameters);
+  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);
   sumEnergy.addLocalEnergy(bulkCosseratEnergy);
   sumEnergy.addLocalEnergy(neumannEnergy);