diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index c6958e50a890e5a995849d64dbe119a447750e6f..36640199e34b47338c444911ddca64527fb807cc 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -8,6 +8,7 @@
 
 #include "localgeodesicfestiffness.hh"
 
+#include <dune/solvers/common/wrapownshare.hh>
 
 /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
  */
@@ -15,6 +16,7 @@ template <class Basis, class TargetSpace>
 class GeodesicFEAssembler {
 
     typedef typename Basis::GridView GridView;
+    using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
 
     //! Dimension of the grid.
     enum { gridDim = GridView::dimension };
@@ -25,22 +27,51 @@ class GeodesicFEAssembler {
     //!
     typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
 
-public:
-    const Basis basis_;
 
 protected:
 
-    LocalGeodesicFEStiffness<Basis,TargetSpace>* localStiffness_;
+    //! The global basis
+    const Basis basis_;
+
+    //! The local stiffness operator
+    std::shared_ptr<LocalStiffness> localStiffness_;
 
 public:
 
     /** \brief Constructor for a given grid */
+    GeodesicFEAssembler(const Basis& basis)
+        : basis_(basis)
+    {}
+
+    /** \brief Constructor for a given grid */
+    template <class LocalStiffnessT>
     GeodesicFEAssembler(const Basis& basis,
-                        LocalGeodesicFEStiffness<Basis, TargetSpace>* localStiffness)
+                        LocalStiffnessT&& localStiffness)
         : basis_(basis),
-          localStiffness_(localStiffness)
+          localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
     {}
 
+    /** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
+    template <class LocalStiffnessT>
+    void setLocalStiffness(LocalStiffnessT&& localStiffness) {
+        localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
+    }
+
+    /** \brief Get the local stiffness operator. */
+    const LocalStiffness& getLocalStiffness() const {
+        return *localStiffness_;
+    }
+
+    /** \brief Get the local stiffness operator. */
+    LocalStiffness& getLocalStiffness() {
+        return *localStiffness_;
+    }
+
+    /** \brief Get the basis. */
+    const Basis& getBasis() const {
+        return basis_;
+    }
+
     /** \brief Assemble the tangent stiffness matrix and the functional gradient together
      *
      * This is more efficient than computing them separately, because you need the gradient
diff --git a/dune/gfe/geodesicfeassemblerwrapper.hh b/dune/gfe/geodesicfeassemblerwrapper.hh
index 8dee1f1ffc276fadf8180a460c13b10eff6f3cc0..4570b4077798ce8fbe3e14a55444e91adf4ec1bc 100644
--- a/dune/gfe/geodesicfeassemblerwrapper.hh
+++ b/dune/gfe/geodesicfeassemblerwrapper.hh
@@ -62,6 +62,11 @@ public:
     /** \brief Get the occupation structure of the Hessian */
     virtual void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
 
+    /** \brief Get the basis. */
+    const ScalarBasis& getBasis() const {
+        return basis_;
+    }
+
 private:
     Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> splitVector(const std::vector<TargetSpace>& sol) const;
     std::unique_ptr<MatrixType> hessianMixed_;
@@ -186,4 +191,4 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
     auto solutionSplit = splitVector(sol);
     return mixedAssembler_->computeEnergy(solutionSplit[_0], solutionSplit[_1]);
 }
-#endif //GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
\ No newline at end of file
+#endif //GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index b45d502d01d5cf19495d1430d9e31c572112b456..88d4e7713aa48b760478a93cbd3ddeb323979acc 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -117,7 +117,7 @@ setup(const GridType& grid,
     // //////////////////////////////////////////////////////////////////////////////////////
 
     typedef DuneFunctionsBasis<Basis> FufemBasis;
-    FufemBasis basis(assembler_->basis_);
+    FufemBasis basis(assembler_->getBasis());
     OperatorAssembler<FufemBasis,FufemBasis> operatorAssembler(basis, basis);
 
     LaplaceAssembler<GridType, typename FufemBasis::LocalFiniteElement, typename FufemBasis::LocalFiniteElement> laplaceStiffness;
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index 797c006e8e167509ee3b61731ecfbbdd6e9b8da8..32300bd3fe334971b8bd850a69255938bb5afbfa 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -98,7 +98,7 @@ getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
 
             double weight = quad[pt].weight() * element.geometry().integrationElement(quadPos);
 
-            FieldVector<double,blocksize> localStrain = dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->getStrain(localSolution, element, quad[pt].position());
+            FieldVector<double,blocksize> localStrain = std::dynamic_pointer_cast<RodLocalStiffness<GridView, double> >(this->localStiffness_)->getStrain(localSolution, element, quad[pt].position());
 
             // Sum it all up
             strain[elementIdx].axpy(weight, localStrain);
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
index a800cc8ba1bc1626a4ca33ce16483b228fea3507..0b113a1d4229cdcf23595c326be9649e35746729 100644
--- a/dune/gfe/rodassembler.hh
+++ b/dune/gfe/rodassembler.hh
@@ -62,7 +62,7 @@ public:
     auto rodEnergy()
     {
       // TODO: Does not work for other stiffness implementations
-      auto localFDStiffness = dynamic_cast<LocalGeodesicFEFDStiffness<Basis, RigidBodyMotion<double,3> >*>(this->localStiffness_);
+      auto localFDStiffness = std::dynamic_pointer_cast<LocalGeodesicFEFDStiffness<Basis, RigidBodyMotion<double,3> > >(this->localStiffness_);
       return const_cast<RodLocalStiffness<GridView,double>*>(dynamic_cast<const RodLocalStiffness<GridView,double>*>(localFDStiffness->localEnergy_));
     }