diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index aba5964acafc955f3cc34b61cc949dbb8c2fdcd9..ec740a9d4c69e0ad55123749aecba19356b639bb 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -160,6 +160,9 @@ energy(const typename Basis::LocalView& localView,
   typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
   LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
 
+  // Bind the density to the current element
+  density_->bind(element);
+
   RT energy = 0;
 
   auto quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
@@ -172,9 +175,6 @@ energy(const typename Basis::LocalView& localView,
     // Local position of the quadrature point
     const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
 
-    // Global position of the quadrature point
-    auto quadPosGlobal = element.geometry().global(quadPos);
-
     const DT integrationElement = geometry.integrationElement(quadPos);
 
     // The value of the local function
@@ -216,7 +216,7 @@ energy(const typename Basis::LocalView& localView,
     // Add the local energy density
     //////////////////////////////////////////////////////////
 
-    const auto energyDensity = (*density_)(quadPosGlobal,
+    const auto energyDensity = (*density_)(quadPos,
                                            aCovariant,
                                            normalGradient,
                                            value,
@@ -267,6 +267,9 @@ energy(const typename Basis::LocalView& localView,
   LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localConfiguration[_0]);
   LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localConfiguration[_1]);
 
+  // Bind the density to the current element
+  density_->bind(element);
+
   RT energy = 0;
 
   auto quadOrder = (deformationLocalFiniteElement.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
@@ -279,9 +282,6 @@ energy(const typename Basis::LocalView& localView,
     // Local position of the quadrature point
     const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
 
-    // Global position of the quadrature point
-    auto quadPosGlobal = element.geometry().global(quadPos);
-
     const DT integrationElement = geometry.integrationElement(quadPos);
 
     // The value of the local function
@@ -335,7 +335,7 @@ energy(const typename Basis::LocalView& localView,
     // Add the local energy density
     //////////////////////////////////////////////////////////
 
-    const auto energyDensity = (*density_)(quadPosGlobal,
+    const auto energyDensity = (*density_)(quadPos,
                                            aCovariant,
                                            normalGradient,
                                            value,
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index bbd18f12a4ba6c43d403181944af2bca612a4a01..be8dc783aad93be3dac5ad4719dd8078627213f2 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -94,6 +94,9 @@ namespace Dune::GFE
         if (not shellBoundary_->contains(it))
           continue;
 
+        // Bind the density to the current intersection
+        density_->bind(it);
+
 #if HAVE_DUNE_CURVEDGEOMETRY
         auto localGridFunction = localFunction(curvedGeometryGridFunction_);
         auto curvedGeometryGridFunctionOrder = deformationLocalFiniteElement.localBasis().order();
@@ -123,9 +126,6 @@ namespace Dune::GFE
           // Local position of the quadrature point
           const auto& quadPos = it.geometryInInside().global(qp.position());
 
-          // Global position of the quadrature point
-          auto quadPosGlobal = it.geometry().global(qp.position());
-
           const DT integrationElement = boundaryGeometry.integrationElement(qp.position());
 
           FieldMatrix<RT, TargetSpace::embeddedDim, boundaryDim> derivative2D;
@@ -181,7 +181,7 @@ namespace Dune::GFE
           // Add the local energy density
           //////////////////////////////////////////////////////////
 
-          const auto energyDensity = (*density_)(quadPosGlobal,
+          const auto energyDensity = (*density_)(qp.position(),
                                                  aCovariant,
                                                  normalGradient,
                                                  value,
diff --git a/dune/gfe/densities/cosseratshelldensity.hh b/dune/gfe/densities/cosseratshelldensity.hh
index a8027211aeeb44b7f81079a3dba9db278cb28dc9..c7c069c6db1c49ea8cb69e14557b0422dc432424 100644
--- a/dune/gfe/densities/cosseratshelldensity.hh
+++ b/dune/gfe/densities/cosseratshelldensity.hh
@@ -36,10 +36,7 @@ namespace Dune::GFE
     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 Position = typename Geometry::LocalCoordinate;
 
     using RigidBodyMotion = ProductManifold<RealTuple<field_type,dimWorld>, Rotation<field_type,dimWorld> >;
     using Derivative = FieldMatrix<field_type,RigidBodyMotion::embeddedDim,domainDim>;
@@ -190,6 +187,15 @@ namespace Dune::GFE
       useAlternativeEnergyWCoss_ = parameters.template get<bool>("useAlternativeEnergyWCoss", false);
     }
 
+    /** \brief Bind the density to a `ElementOrIntersection` object
+     *
+     * Stores a copy of the `ElementOrIntersection`'s geometry.
+     **/
+    void bind(const ElementOrIntersection& elementOrIntersection)
+    {
+      geometry_.emplace(elementOrIntersection.geometry());
+    }
+
     /** \brief Evaluation with the current position, the deformation function, the deformation gradient, the rotation and the rotation gradient
      *
      * \param x The current position
@@ -204,8 +210,9 @@ namespace Dune::GFE
                            const RigidBodyMotion& value,
                            const Derivative& derivative) const
     {
-      double thickness = thicknessF_(x);
-      auto lameConstants = lameF_(x);
+      auto xGlobal = geometry_->global(x);
+      double thickness = thicknessF_(xGlobal);
+      auto lameConstants = lameF_(xGlobal);
       auto mu = lameConstants[0];
       auto lambda = lameConstants[1];
       // TODO: Use structured binding here
@@ -345,6 +352,9 @@ namespace Dune::GFE
     /** \brief Whether to use the alternative energy W_Coss from Birsan 2021:
                "Alternative derivation of the higher-order constitutive model for six-parameter elastic shells", equations (119) and (126). */
     bool useAlternativeEnergyWCoss_;
+
+    /** \brief The geometry that this density is defined on */
+    std::optional<Geometry> geometry_ = std::nullopt;
   };
 
 }  // namespace Dune::GFE