diff --git a/dune/gfe/densities/CMakeLists.txt b/dune/gfe/densities/CMakeLists.txt
index a43ea970523a3bd520be7a63aac3542b52276c8c..d8b787647800ca37e99681957072679044c85485 100644
--- a/dune/gfe/densities/CMakeLists.txt
+++ b/dune/gfe/densities/CMakeLists.txt
@@ -3,6 +3,7 @@ install(FILES
         bulkcosseratdensity.hh
         chiralskyrmiondensity.hh
         cosseratshelldensity.hh
+        cosseratvolumeloaddensity.hh
         duneelasticitydensity.hh
         harmonicdensity.hh
         localdensity.hh
diff --git a/dune/gfe/densities/cosseratvolumeloaddensity.hh b/dune/gfe/densities/cosseratvolumeloaddensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f90a2f2f775c778fb8c95d7e9e533c41e5dfb24d
--- /dev/null
+++ b/dune/gfe/densities/cosseratvolumeloaddensity.hh
@@ -0,0 +1,74 @@
+#ifndef DUNE_GFE_DENSITIES_COSSERATVOLUMELOADDENSITY_HH
+#define DUNE_GFE_DENSITIES_COSSERATVOLUMELOADDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/gfe/densities/localdensity.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+namespace Dune::GFE
+{
+  template<class Position, class field_type>
+  class CosseratVolumeLoadDensity final
+    : public GFE::LocalDensity<Position, GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
+  {
+    static constexpr int gridDim = Position::size();
+
+    // The target space with 'adouble' as the number type
+    using ATargetSpace = GFE::ProductManifold<RealTuple<adouble,3>,Rotation<adouble,3> >;
+
+    /** \brief The function implementing a volume load */
+    const std::function<FieldVector<double,3>(FieldVector<double,gridDim>)> volumeLoad_;
+
+  public:
+
+    /** \brief Constructor with a given load density
+     */
+    CosseratVolumeLoadDensity(const std::function<FieldVector<double,3>(FieldVector<double,gridDim>)>& volumeLoad)
+      : volumeLoad_(volumeLoad)
+    {}
+
+    /** \brief Evaluate the density
+     */
+    virtual field_type operator() (const Position& x,
+                                   const typename GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >::CoordinateType& value,
+                                   const FieldMatrix<field_type,7,gridDim>& derivative) const override
+    {
+      field_type density = 0;
+
+      // Value of the volume load density at the current position
+      auto loadVector = volumeLoad_(x);
+
+      // In this implementation, only translational dofs are affected by the volume load
+      for (size_t i=0; i<loadVector.size(); i++)
+        density += loadVector[i] * value[i];
+
+      return density;
+    }
+
+    // Construct a copy of this density but using 'adouble' as the number type
+    virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
+    {
+      return std::make_unique<CosseratVolumeLoadDensity<Position,adouble> >(volumeLoad_);
+    }
+
+    /** \brief The density depends on the deformation value, but not on the microrotation value
+     */
+    virtual bool dependsOnValue(int factor=-1) const override
+    {
+      return factor==-1 || factor==0;
+    }
+
+    /** \brief The density neither depends on the deformation gradient nor on the microrotation gradient
+     */
+    virtual bool dependsOnDerivative([[maybe_unused]] int factor=-1) const override
+    {
+      return false;
+    }
+  };
+
+} // namespace GFE
+
+#endif   //#ifndef DUNE_GFE_BULKCOSSERATDENSITY_HH
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 1d75dd1cd717ca6cb5d87f102ff81bec0c7f1b0c..5c1c451d3428a1a315d8cc6b5dd66d50e66b0bd1 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -69,6 +69,7 @@
 #include <dune/gfe/assemblers/mixedgfeassembler.hh>
 #include <dune/gfe/assemblers/sumenergy.hh>
 #include <dune/gfe/densities/bulkcosseratdensity.hh>
+#include <dune/gfe/densities/cosseratvolumeloaddensity.hh>
 #include <dune/gfe/densities/planarcosseratshelldensity.hh>
 
 #if MIXED_SPACE
@@ -498,12 +499,6 @@ int main (int argc, char *argv[]) try
                              return nV;
                            };
 
-    if (parameterSet.get<std::string>("volumeLoadPythonFunction", "zero-volume-load") != "zero-volume-load")
-    {
-      std::cerr << "cosserat-continuum.cc: Volume loads are not fully implemented yet." << std::endl;
-      std::abort();
-    }
-
     Python::Reference volumeLoadClass = Python::import(parameterSet.get<std::string>("volumeLoadPythonFunction", "zero-volume-load"));
 
     Python::Callable volumeLoadCallable = volumeLoadClass.get("VolumeLoad");
@@ -574,6 +569,16 @@ int main (int argc, char *argv[]) try
       auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,3>, Rotation<adouble,3> > >(neumannBoundary,neumannFunction);
       sumEnergy->addLocalEnergy(neumannEnergy);
 
+      // The volume load term
+      // TODO: There is a bug here: The volumeLoad function is currently defined in global coordinates,
+      // but the density expects it to be in local coordinates with respect to the element being
+      // integrated over.  I have to think about where the binding should happen.
+      // Practically, the bug does not really show, because all our volume loads
+      // are constant in space anyway.
+      auto volumeLoadDensity = std::make_shared<GFE::CosseratVolumeLoadDensity<LocalCoordinate,adouble> >(volumeLoad);
+      auto volumeLoadEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, AInterpolationRule, ATargetSpace> >(volumeLoadDensity);
+      sumEnergy->addLocalEnergy(volumeLoadEnergy);
+
       // The local assembler
       LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy,
                                                                                        adolcScalarMode);