From 3b39a00cea57e76f17555f02868dc0b55c00bf8f Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Tue, 17 Jan 2023 18:55:20 +0100
Subject: [PATCH] Repair the call to the volume load in
 cosseratenergystiffness.hh

---
 dune/gfe/cosseratenergystiffness.hh | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 3ddc6961..2f7c2523 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/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;
+
     }
 
     //////////////////////////////////////////////////////////////////////////////
-- 
GitLab