diff --git a/dune/gfe/chiralskyrmionenergy.hh b/dune/gfe/chiralskyrmionenergy.hh
index 75b550b40da3c32bf9f961983450ff2536212e04..5756c3ec23c5d837006ec59ee84aff675b0de25c 100644
--- a/dune/gfe/chiralskyrmionenergy.hh
+++ b/dune/gfe/chiralskyrmionenergy.hh
@@ -4,7 +4,7 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localenergy.hh>
 
 namespace Dune {
 
@@ -17,7 +17,7 @@ namespace GFE {
  */
 template<class Basis, class LocalInterpolationRule, class field_type>
 class ChiralSkyrmionEnergy
-: public LocalGeodesicFEStiffness<Basis,UnitVector<field_type,3> >
+: public GFE::LocalEnergy<Basis,UnitVector<field_type,3> >
 {
   // various useful types
   typedef UnitVector<field_type,3> TargetSpace;
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index bcc3b95092dbabde8eb16a88d23ee2c6995eb0a0..ead0fa54cd5816b2e022bc61c48b6be73e4ea2aa 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -8,7 +8,7 @@
 #include <dune/fufem/functions/virtualgridfunction.hh>
 #include <dune/fufem/boundarypatch.hh>
 
-#include "localgeodesicfestiffness.hh"
+#include <dune/gfe/localenergy.hh>
 #include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
 #ifdef PROJECTED_INTERPOLATION
 #include <dune/gfe/localprojectedfefunction.hh>
@@ -58,7 +58,7 @@ public:
 
 template<class Basis, int dim, class field_type=double>
 class CosseratEnergyLocalStiffness
-    : public LocalGeodesicFEStiffness<Basis,RigidBodyMotion<field_type,dim> >,
+    : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >,
       public MixedLocalGeodesicFEStiffness<Basis,
                                            RealTuple<field_type,dim>,
                                            Rotation<field_type,dim> >
diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index 253bd102fda9c3fd1f60311e1e84379a306d24ad..b9f561c9009977a13c54483563bce55786437b39 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -4,11 +4,11 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/geometry/quadraturerules.hh>
 
-#include "localgeodesicfestiffness.hh"
+#include <dune/gfe/localenergy.hh>
 
 template<class Basis, class LocalInterpolationRule, class TargetSpace>
 class HarmonicEnergyLocalStiffness
-    : public LocalGeodesicFEStiffness<Basis,TargetSpace>
+    : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
 {
     // grid types
     typedef typename Basis::GridView GridView;
diff --git a/dune/gfe/l2distancesquaredenergy.hh b/dune/gfe/l2distancesquaredenergy.hh
index 5e1e5a988c52a65d2b47c601c585238ab62c24c3..911abfabf3823b96425a02a13b958ff1e6390c20 100644
--- a/dune/gfe/l2distancesquaredenergy.hh
+++ b/dune/gfe/l2distancesquaredenergy.hh
@@ -5,12 +5,12 @@
 
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 
 template<class Basis, class TargetSpace>
 class L2DistanceSquaredEnergy
-  : public LocalGeodesicFEStiffness<Basis,TargetSpace>
+  : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
 {
   // grid types
   typedef typename Basis::GridView GridView;
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index c8c48f7b076e95fd6093144b37cfb1792d74552d..c1d098809c9b3faee26622d99d1bb2e768f67899 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -10,7 +10,7 @@
 #include <dune/fufem/functions/virtualgridfunction.hh>
 #include <dune/fufem/boundarypatch.hh>
 
-#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 #include <dune/gfe/rigidbodymotion.hh>
 #include <dune/gfe/unitvector.hh>
@@ -20,7 +20,7 @@
 
 template<class Basis, int dim, class field_type=double>
 class NonplanarCosseratShellEnergy
-  : public LocalGeodesicFEStiffness<Basis,RigidBodyMotion<field_type,dim> >
+  : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >
 {
   // grid types
   typedef typename Basis::GridView GridView;
diff --git a/dune/gfe/sumcosseratenergy.hh b/dune/gfe/sumcosseratenergy.hh
index f3e1e090b940435f15208a68a9ece867b6306571..2011d909604194c646996056ea814ad55253d8ca 100644
--- a/dune/gfe/sumcosseratenergy.hh
+++ b/dune/gfe/sumcosseratenergy.hh
@@ -3,14 +3,14 @@
 
 #include <dune/elasticity/assemblers/localfestiffness.hh>
 
-#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localenergy.hh>
 
 namespace Dune {
 
 namespace GFE {
 template<class Basis, class TargetSpace, class field_type=double>
 class SumCosseratEnergy
-: public LocalGeodesicFEStiffness<Basis, TargetSpace>
+: public GFE::LocalEnergy<Basis, TargetSpace>
 {
  // grid types
   typedef typename Basis::GridView GridView;
@@ -27,8 +27,8 @@ public:
    * \param elasticEnergy The elastic energy
    * \param cosseratEnergy The cosserat energy
    */
-  SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
-            std::shared_ptr<LocalGeodesicFEStiffness<Basis, TargetSpace>> cosseratEnergy)
+  SumCosseratEnergy(std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
+            std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace>> cosseratEnergy)
 
   : elasticEnergy_(elasticEnergy),
     cosseratEnergy_(cosseratEnergy)
@@ -49,9 +49,9 @@ public:
 
 private:
 
-  std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
+  std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
 
-  std::shared_ptr<LocalGeodesicFEStiffness<Basis, TargetSpace> > cosseratEnergy_;
+  std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace> > cosseratEnergy_;
 };
 
 }  // namespace GFE
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index 905ab45ddfe66282d160f4e5e290078dd87db5f5..7d8e85c1d91248abd241cda9dc849282e04b652f 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -7,6 +7,7 @@
 #include <dune/fufem/boundarypatch.hh>
 
 #include <dune/gfe/cosseratstrain.hh>
+#include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 #include <dune/gfe/localprojectedfefunction.hh>
 #include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
@@ -42,7 +43,7 @@ public:
 
 template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
 class SurfaceCosseratEnergy
-: public LocalGeodesicFEStiffness<Basis,TargetSpace>
+: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
 {
  // grid types
   typedef typename Basis::GridView GridView;
diff --git a/dune/gfe/weightedsumenergy.hh b/dune/gfe/weightedsumenergy.hh
index 38d5cbef8ed3f96f9486d087cff060dbdda8f899..57795057f3bd89ceb50e7c786511d13dbec86def 100644
--- a/dune/gfe/weightedsumenergy.hh
+++ b/dune/gfe/weightedsumenergy.hh
@@ -6,12 +6,12 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 
 template<class Basis, class TargetSpace>
 class WeightedSumEnergy
-  : public LocalGeodesicFEStiffness<Basis,TargetSpace>
+  : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
 {
   // grid types
   typedef typename Basis::GridView GridView;
@@ -23,11 +23,11 @@ class WeightedSumEnergy
 
 public:
 
-  std::vector<std::shared_ptr<LocalGeodesicFEStiffness<Basis,TargetSpace> > > addends_;
+  std::vector<std::shared_ptr<Dune::GFE::LocalEnergy<Basis,TargetSpace> > > addends_;
 
   std::vector<double> weights_;
 
-  WeightedSumEnergy(std::vector<std::shared_ptr<LocalGeodesicFEStiffness<Basis,TargetSpace> > > addends,
+  WeightedSumEnergy(std::vector<std::shared_ptr<Dune::GFE::LocalEnergy<Basis,TargetSpace> > > addends,
                     std::vector<double> weights)
   : addends_(addends),
     weights_(weights)
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 32f7d64b274bef137f7f622e0f17c16d5762fab3..832d3955894819ccfb705c659bfe8a9cf5870082 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -452,29 +452,27 @@ int main (int argc, char *argv[]) try
                       RealTuple<double,3>,
                       Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
 #else
-    using LocalEnergyBase = LocalGeodesicFEStiffness<FEBasis,RigidBodyMotion<adouble,3> >;
-
-    std::shared_ptr<LocalEnergyBase> cosseratEnergyADOLCLocalStiffness;
+    std::shared_ptr<GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble,3> > > localCosseratEnergy;
 
     if (dim==dimworld)
     {
-      cosseratEnergyADOLCLocalStiffness = std::make_shared<CosseratEnergyLocalStiffness<FEBasis,3,adouble> >(materialParameters,
-                                                                                                             &neumannBoundary,
-                                                                                                             neumannFunction,
-                                                                                                             volumeLoad);
+      localCosseratEnergy = std::make_shared<CosseratEnergyLocalStiffness<FEBasis,3,adouble> >(materialParameters,
+                                                                                               &neumannBoundary,
+                                                                                               neumannFunction,
+                                                                                               volumeLoad);
     }
     else
     {
       std::vector<UnitVector<double,3> > vertexNormals = computeVertexNormals(gridView);
-      cosseratEnergyADOLCLocalStiffness = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
-                                                                                                             std::move(vertexNormals),
-                                                                                                             &neumannBoundary,
-                                                                                                             neumannFunction,
-                                                                                                             volumeLoad);
+      localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
+                                                                                               std::move(vertexNormals),
+                                                                                               &neumannBoundary,
+                                                                                               neumannFunction,
+                                                                                               volumeLoad);
     }
 
     LocalGeodesicFEADOLCStiffness<FEBasis,
-                                  TargetSpace> localGFEADOLCStiffness(cosseratEnergyADOLCLocalStiffness.get());
+                                  TargetSpace> localGFEADOLCStiffness(localCosseratEnergy.get());
 
     GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
 #endif
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 8ff8f6de893868bbed3d59a23e2a925a64dfbf6c..fd9e1fbe4fdaa05ccc06a5eac834b2d739a3b221 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -336,7 +336,7 @@ int main (int argc, char *argv[]) try
 
     // Assembler using ADOL-C
     std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-    std::shared_ptr<LocalFEStiffness<GridView,
+    std::shared_ptr<Elasticity::LocalEnergy<GridView,
                                      FEBasis::LocalView::Tree::FiniteElement,
                                      std::vector<Dune::FieldVector<ValueType, dim>> > > elasticEnergy;
 
@@ -376,7 +376,7 @@ int main (int argc, char *argv[]) try
           ValueType>>(elasticEnergy, neumannEnergy);
 
 
-    using LocalEnergyBase = LocalGeodesicFEStiffness<FEBasis,RigidBodyMotion<adouble, dim> >;
+    using LocalEnergyBase = GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble, dim> >;
 
     std::shared_ptr<LocalEnergyBase> surfaceCosseratEnergy;
     std::vector<UnitVector<double,3> > vertexNormals(gridView.size(3));
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index d74c21a2b6237013323989dfff01bf0b8e0fde5f..e23ec346899f4642cd8f8a85d0297903527b46de 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -201,7 +201,7 @@ int main (int argc, char *argv[]) try
 
   auto l2DistanceSquaredEnergy = std::make_shared<L2DistanceSquaredEnergy<FEBasis, ATargetSpace> >();
 
-  std::vector<std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > > addends(2);
+  std::vector<std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > > addends(2);
   using GeodesicInterpolationRule  = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
   addends[0] = std::make_shared<HarmonicEnergyLocalStiffness<FEBasis, GeodesicInterpolationRule, ATargetSpace> >();
   addends[1] = l2DistanceSquaredEnergy;
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 0700d120b9630af7629c905425094923d5ce34ed..d216f704bf9f0a6721b6a0660429ce61c2d2545b 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -272,7 +272,7 @@ int main (int argc, char *argv[])
     using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
 
     // Assembler using ADOL-C
-    std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > localEnergy;
+    std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
 
     std::string energy = parameterSet.get<std::string>("energy");
     if (energy == "harmonic")