diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 0f33ad4b3ca81141d48a09b86dbf357bdaed2a3c..989e5e50b7b0047dec566d3a7950df250c37fc7d 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -189,7 +189,8 @@ public:
      */
     CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
                                  const BoundaryPatch<GridView>* neumannBoundary,
-                                 const Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> >* neumannFunction)
+                                 const Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> >* neumannFunction,
+                                 const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
     : neumannBoundary_(neumannBoundary),
       neumannFunction_(neumannFunction)
     {
@@ -363,6 +364,9 @@ public:
 
     /** \brief The function implementing the Neumann data */
     const Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> >* neumannFunction_;
+
+    /** \brief The function implementing a volume load */
+    const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
 };
 
 template <class Basis, int dim, class field_type>
@@ -448,6 +452,24 @@ energy(const typename Basis::LocalView& localView,
         } else
             DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
 
+        ///////////////////////////////////////////////////////////
+        // Volume load contribution
+        ///////////////////////////////////////////////////////////
+
+        if (not volumeLoad_)
+            continue;
+
+        // Value of the volume load density at the current position
+        Dune::FieldVector<double,3> volumeLoadDensity;
+
+        if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
+            std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
+        else
+            volumeLoad_->evaluate(element.geometry().global(quad[pt].position()), volumeLoadDensity);
+
+        // Only translational dofs are affected by the volume load
+        for (size_t i=0; i<volumeLoadDensity.size(); i++)
+            energy -= thickness_ * (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
     }
 
 
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index c2126be16a9e51ff588a9f58dd77f2ccc98e384f..75eea0156923c331b665229540474f24317d4647 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -76,6 +76,25 @@ struct NeumannFunction
     double homotopyParameter_;
 };
 
+/** \brief A constant vector-valued function, for simple volume loads */
+struct VolumeLoad
+    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
+{
+    VolumeLoad(const FieldVector<double,3> values,
+               double homotopyParameter)
+    : values_(values),
+      homotopyParameter_(homotopyParameter)
+    {}
+
+    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
+        out = 0;
+        out.axpy(homotopyParameter_, values_);
+    }
+
+    FieldVector<double,3> values_;
+    double homotopyParameter_;
+};
+
 
 int main (int argc, char *argv[]) try
 {
@@ -260,6 +279,10 @@ int main (int argc, char *argv[]) try
         neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
                                                        homotopyParameter);
 
+        shared_ptr<VolumeLoad> volumeLoad;
+        if (parameterSet.hasKey("volumeLoad"))
+            volumeLoad = make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
+                                                                                          homotopyParameter);
 
         if (mpiHelper.rank() == 0) {
             std::cout << "Material parameters:" << std::endl;
@@ -270,7 +293,8 @@ int main (int argc, char *argv[]) try
     CosseratEnergyLocalStiffness<FEBasis,
                                  3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
                                                                               &neumannBoundary,
-                                                                              neumannFunction.get());
+                                                                              neumannFunction.get(),
+                                                                              volumeLoad);
     LocalGeodesicFEADOLCStiffness<FEBasis,
                                   TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);