diff --git a/dune/gfe/localfirstordermodel.hh b/dune/gfe/localfirstordermodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a098c5bacc3c549d4caea909e6815fb3111c0fab
--- /dev/null
+++ b/dune/gfe/localfirstordermodel.hh
@@ -0,0 +1,30 @@
+#ifndef DUNE_GFE_LOCALFIRSTORDERMODEL_HH
+#define DUNE_GFE_LOCALFIRSTORDERMODEL_HH
+
+#include <dune/gfe/localenergy.hh>
+
+namespace Dune {
+
+namespace GFE {
+
+template<class Basis, class TargetSpace>
+class LocalFirstOrderModel
+: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+{
+public:
+
+    /** \brief Assemble the element gradient of the energy functional
+
+    The default implementation in this class uses a finite difference approximation */
+    virtual void assembleGradient(const typename Basis::LocalView& localView,
+                                  const std::vector<TargetSpace>& solution,
+                                  std::vector<typename TargetSpace::TangentVector>& gradient) const = 0;
+
+};
+
+}  // namespace GFE
+
+}  // namespace Dune
+
+#endif   // DUNE_GFE_LOCALFIRSTORDERMODEL_HH
+
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
index ffc5d222a0f6dd737532e5cde0fecd1570bb4b1b..cf93563fc6770751eee0957d85411165065d1c61 100644
--- a/dune/gfe/rodassembler.hh
+++ b/dune/gfe/rodassembler.hh
@@ -8,6 +8,7 @@
 
 #include <dune/fufem/boundarypatch.hh>
 
+#include <dune/gfe/localgeodesicfefdstiffness.hh>
 #include "rigidbodymotion.hh"
 #include "rodlocalstiffness.hh"
 #include "geodesicfeassembler.hh"
@@ -41,8 +42,11 @@ public:
         //! ???
     RodAssembler(const Basis& basis,
                  RodLocalStiffness<GridView,double>* localStiffness)
-        : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,localStiffness)
+        : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,nullptr)
+        , rodEnergy_(localStiffness)
         {
+            this->localStiffness_ = new LocalGeodesicFEFDStiffness<Basis,RigidBodyMotion<double,3>, double>(localStiffness);
+
             std::vector<RigidBodyMotion<double,3> > referenceConfiguration(basis.size());
 
     for (const auto vertex : Dune::vertices(basis.gridView()))
@@ -55,11 +59,12 @@ public:
                 referenceConfiguration[idx].q = Rotation<double,3>::identity();
             }
 
-            dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->setReferenceConfiguration(referenceConfiguration);
+    rodEnergy_->setReferenceConfiguration(referenceConfiguration);
         }
 
         std::vector<RigidBodyMotion<double,3> > getRefConfig()
-        {   return  dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->referenceConfiguration_;
+    {
+      return rodEnergy_->referenceConfiguration_;
         }
 
   virtual void assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
@@ -78,6 +83,7 @@ public:
         Dune::FieldVector<double,6> getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
                                                       const std::vector<RigidBodyMotion<double,3> >& sol) const;
 
+    RodLocalStiffness<GridView,double>* rodEnergy_;
     }; // end class
 
 
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index 508e18f032a542bd3efd97af9e08c1a807588e08..5157a30e3f7a4e4d67f292c938a13ccc505b8f2d 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -9,12 +9,12 @@
 
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 
-#include "localgeodesicfestiffness.hh"
+#include <dune/gfe/localfirstordermodel.hh>
 #include "rigidbodymotion.hh"
 
 template<class GridView, class RT>
 class RodLocalStiffness
-    : public LocalGeodesicFEStiffness<Dune::Functions::LagrangeBasis<GridView,1>, RigidBodyMotion<RT,3> >
+: public Dune::GFE::LocalFirstOrderModel<Dune::Functions::LagrangeBasis<GridView,1>, RigidBodyMotion<RT,3> >
 {
     typedef RigidBodyMotion<RT,3> TargetSpace;
     typedef Dune::Functions::LagrangeBasis<GridView,1> Basis;