Repair the call to the volume load in cosseratenergystiffness.hh

dune/gfe/cosseratenergystiffness.hh

@@ -99,7 +99,8 @@ public:
                                  const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
                                  const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
     : neumannBoundary_(neumannBoundary),
-      neumannFunction_(neumannFunction)
+      neumannFunction_(neumannFunction),
+      volumeLoad_(volumeLoad)
     {
         // The shell thickness
         thickness_ = parameters.template get<double>("thickness");
@@ -513,6 +514,20 @@ 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
+        auto volumeLoadDensity = volumeLoad_(element.geometry().global(quad[pt].position()));
+
+
+        // Only translational dofs are affected by the volume load
+        for (size_t i=0; i<volumeLoadDensity.size(); i++)
+            energy += (volumeLoadDensity[i] * deformationValue[i]) * quad[pt].weight() * integrationElement;
+
     }
 
     //////////////////////////////////////////////////////////////////////////////