diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index 0dc07d7d7cbe2e0e828ed46812a32329dbfb42de..aba5964acafc955f3cc34b61cc949dbb8c2fdcd9 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -79,21 +79,19 @@ class NonplanarCosseratShellEnergy
   typedef typename GridView::ctype DT;
   typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
   typedef typename TargetSpace::ctype RT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
+  typedef typename GridView::template Codim<0>::Entity Element;
 
   // some other sizes
   constexpr static int gridDim = GridView::dimension;
   constexpr static int dimworld = GridView::dimensionworld;
 
-  using Position = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
-
 public:
 
   /** \brief Constructor with a set of material parameters
    * \param parameters                  The material parameters
    * \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
    */
-  NonplanarCosseratShellEnergy(const std::shared_ptr<Dune::GFE::CosseratShellDensity<Position,field_type> >& density,
+  NonplanarCosseratShellEnergy(const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> >& density,
                                const StressFreeStateGridFunction* stressFreeStateGridFunction)
     : stressFreeStateGridFunction_(stressFreeStateGridFunction),
     density_(density)
@@ -123,7 +121,7 @@ public:
   const StressFreeStateGridFunction* stressFreeStateGridFunction_;
 
   /** \brief The energy density of a Cosserat shell with nonplanar reference configuration */
-  const std::shared_ptr<Dune::GFE::CosseratShellDensity<Position,field_type> > density_;
+  const std::shared_ptr<Dune::GFE::CosseratShellDensity<Element,field_type> > density_;
 };
 
 template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index b9c57d73b27dc63a99c96515a37f1c1c583f4be2..bbd18f12a4ba6c43d403181944af2bca612a4a01 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -37,6 +37,8 @@ namespace Dune::GFE
     constexpr static int gridDim = GridView::dimension;
     static constexpr int boundaryDim = gridDim - 1;
 
+    using Intersection = typename GridView::Intersection;
+
   public:
 
     /** \brief Constructor with a set of material parameters
@@ -45,7 +47,7 @@ namespace Dune::GFE
      * \param curvedGeometryGridFunction The curvedGeometryGridFunction gives the geometry of the shell in stress-free state.
               When assembling, we deform the intersections using the curvedGeometryGridFunction and then use the deformed geometries.
      */
-    SurfaceCosseratEnergy(const std::shared_ptr<CosseratShellDensity<FieldVector<DT,3>,RT> >& density,
+    SurfaceCosseratEnergy(const std::shared_ptr<CosseratShellDensity<Intersection,RT> >& density,
                           const BoundaryPatch<GridView>* shellBoundary,
                           const CurvedGeometryGridFunction& curvedGeometryGridFunction)
       : shellBoundary_(shellBoundary),
@@ -199,7 +201,7 @@ namespace Dune::GFE
     const CurvedGeometryGridFunction curvedGeometryGridFunction_;
 
     /** \brief The density that is being integrated */
-    const std::shared_ptr<CosseratShellDensity<FieldVector<DT,3>,RT> > density_;
+    const std::shared_ptr<CosseratShellDensity<Intersection,RT> > density_;
 
   };
 }  // namespace Dune::GFE
diff --git a/dune/gfe/densities/cosseratshelldensity.hh b/dune/gfe/densities/cosseratshelldensity.hh
index 447b5090fd98e93fea2d2c69053869457ea25879..a8027211aeeb44b7f81079a3dba9db278cb28dc9 100644
--- a/dune/gfe/densities/cosseratshelldensity.hh
+++ b/dune/gfe/densities/cosseratshelldensity.hh
@@ -24,16 +24,21 @@ namespace Dune::GFE
    *
    * (and earlier papers by Neff and Bîrsan).
    *
-   * \tparam Position Type for the points where the density will be evaluated
+   * \tparam ElementOrIntersection The domain of the density.
+   *   Can be either a grid element (i.e., a codimension-0 Entity)
+   *   or an Intersection.
    * \tparam field_type Type used for numbers
    */
-  template<class Position, class field_type>
+  template<class ElementOrIntersection, class field_type>
   class CosseratShellDensity
   {
     constexpr static int dimWorld = 3;
     constexpr static int domainDim = 2;
 
+    using Geometry = typename ElementOrIntersection::Geometry;
+
     // TODO: This is likely not a good long-term solution
+    using Position = Geometry::GlobalCoordinate;
     static_assert(Position::size()==dimWorld, "Position must be in world-coordinates.");
 
     using RigidBodyMotion = ProductManifold<RealTuple<field_type,dimWorld>, Rotation<field_type,dimWorld> >;
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index f656ae22e382ada3b830110d66e850be3571be60..64c11c35f118ec60c673e6a588747031b28de1bb 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -132,8 +132,8 @@ auto createCosseratEnergy(const ParameterTree& materialParameters,
   }
   else if constexpr (LocalCoordinate::size()==2 && dimworld==3)
   {
-    using GlobalCoordinate = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
-    auto density = std::make_shared<GFE::CosseratShellDensity<GlobalCoordinate, adouble> >(materialParameters);
+    using Element = typename GridView::template Codim<0>::Entity;
+    auto density = std::make_shared<GFE::CosseratShellDensity<Element, adouble> >(materialParameters);
 
     return std::make_shared<NonplanarCosseratShellEnergy<Basis, 3, adouble, decltype(creator)> >(density, &creator);
   }
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index e9684c2d9dd490ced157fc9d8caf41770c9f2c6f..e822c73f6e7df17bf9f570fea6775282526ef400 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -531,8 +531,8 @@ int main (int argc, char *argv[]) try
     auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunctionPtr);
 
     // The energy of the surface shell
-    auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<
-        FieldVector<double,3>, adouble> >(
+    using Intersection = typename GridView::Intersection;
+    auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<Intersection, adouble> >(
       materialParameters,
       fThickness,
       fLame);
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index d3f8b83c91137767c8873095841d438d0cc8e895..80e885799a6efda21c78f8669026793df67bbd4e 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -356,8 +356,8 @@ int main (int argc, char *argv[])
   auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
 
-  auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<
-      FieldVector<double,3>, adouble> >(
+  using Intersection = typename GridView::Intersection;
+  auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<Intersection, adouble> >(
     materialParameters,
     fThickness,
     fLame);
diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc
index 8f3ab0872b76b052260844fbb85f4d648b73e23b..6566c170a14b7c625bd9334cafc6663fccd91740 100644
--- a/test/nonplanarcosseratenergytest.cc
+++ b/test/nonplanarcosseratenergytest.cc
@@ -116,8 +116,8 @@ double calculateEnergy(const FlatGridView& flatGridView,
   //  Construct the energy functional
   ///////////////////////////////////////////////////
 
-  using GlobalCoordinate = typename FlatGridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
-  auto density = std::make_shared<GFE::CosseratShellDensity<GlobalCoordinate, double> >(materialParameters);
+  using Element = typename FlatGridView::template Codim<0>::Entity;
+  auto density = std::make_shared<GFE::CosseratShellDensity<Element, double> >(materialParameters);
 
   using ShellEnergy = NonplanarCosseratShellEnergy<FlatFEBasis,
       3,                                               // Dimension of the target space