diff --git a/CHANGELOG.md b/CHANGELOG.md
index b7a71f58d5c52cc4112fddb1bc74a5f80b2965f9..5f2f6dd8e3609ae8df829bd6d8633349c52af710 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,10 @@
 # Master
 
+- Do not scale the density functions  with the thickness to avoid confusion
+  since some densities need to be scaled and some do not need to be scaled
+  with the thickness depending on the dimension of the grid, their direction
+  and their kind (Neumann or volume load).
+
 - Fix bug in the `RealTuple::log` method: Calling `log(a,b)`returned `a-b`
   instead of `b-a`.
 
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 6f308de6207222da4808764ee1089b9b0489d9a8..3ddc696106637e0629a5d5f1efcba15aeaed56d1 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -371,7 +371,7 @@ energy(const typename Basis::LocalView& localView,
 
         // 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;
+            energy += (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
     }
 
 
@@ -405,7 +405,7 @@ energy(const typename Basis::LocalView& localView,
 
             // Only translational dofs are affected by the Neumann force
             for (size_t i=0; i<neumannValue.size(); i++)
-                energy += thickness_ * (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+                energy += (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
 
         }
 
@@ -544,7 +544,7 @@ energy(const typename Basis::LocalView& localView,
 
             // Only translational dofs are affected by the Neumann force
             for (size_t i=0; i<neumannValue.size(); i++)
-                energy += thickness_ * (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
+                energy += (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
 
         }
 
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index e8a3678eed7a6c7880bea11b2a400e9ab451bb53..01575c5e28d6ba0b10eadb1a14eed83164771247 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -332,7 +332,7 @@ energy(const typename Basis::LocalView& localView,
 
     // 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;
+      energy += (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
   }
 
 
@@ -365,7 +365,7 @@ energy(const typename Basis::LocalView& localView,
 
       // Only translational dofs are affected by the Neumann force
       for (size_t i=0; i<neumannValue.size(); i++)
-        energy += thickness_ * (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+        energy += (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
     }
 
   }