diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 0324629b6f6603a7d0c6d56932b0f645e80dcd55..4937720c4adabd0abb1ce46a66eed4ab1a2605be 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -22,10 +22,10 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 #endif
 
-template <class Basis, class TargetSpace>
-void RiemannianProximalNewtonSolver<Basis,TargetSpace>::
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
 setup(const GridType& grid,
-      const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+      const Assembler* assembler,
       const SolutionType& x,
       const Dune::BitSetVector<blocksize>& dirichletNodes,
       double tolerance,
@@ -150,8 +150,8 @@ setup(const GridType& grid,
 }
 
 
-template <class Basis, class TargetSpace>
-void RiemannianProximalNewtonSolver<Basis,TargetSpace>::solve()
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
 {
     int rank = grid_->comm().rank();
 
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index a1c1aca5605b2e1c6fd8a579f562295ecfc39f9a..1024a72a7aaa1fe36268c2a84a02f98af279c8b4 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -17,7 +17,6 @@
 
 #include <dune/gfe/periodic1dpq1nodalbasis.hh>
 
-#include "geodesicfeassembler.hh"
 #include "riemanniantrsolver.hh"
 #include <dune/grid/utility/globalindexset.hh>
 #include <dune/gfe/parallel/globalmapper.hh>
@@ -26,7 +25,7 @@
 #include <dune/gfe/parallel/p2mapper.hh>
 
 /** \brief Riemannian proximal-newton solver for geodesic finite-element problems */
-template <class Basis, class TargetSpace>
+template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
 class RiemannianProximalNewtonSolver
     : public IterativeSolver<std::vector<TargetSpace>,
                              Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
@@ -70,7 +69,7 @@ public:
 
     /** \brief Set up the solver using a choldmod solver as the inner solver */
     void setup(const GridType& grid,
-               const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+               const Assembler* assembler,
                const SolutionType& x,
                const Dune::BitSetVector<blocksize>& dirichletNodes,
                double tolerance,
@@ -121,7 +120,7 @@ protected:
     std::unique_ptr<MatrixType> hessianMatrix_;
 
     /** \brief The assembler for the material law */
-    const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
+    const Assembler* assembler_;
 
     /** \brief The solver for the quadratic inner problems */
     std::shared_ptr<typename Dune::Solvers::CholmodSolver<MatrixType,CorrectionType>> innerSolver_;
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index b8604dfc6c94f033342cbc61addc33454fd577fa..b45d502d01d5cf19495d1430d9e31c572112b456 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -28,10 +28,10 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 #endif
 
-template <class Basis, class TargetSpace>
-void RiemannianTrustRegionSolver<Basis,TargetSpace>::
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
 setup(const GridType& grid,
-      const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+      const Assembler* assembler,
          const SolutionType& x,
          const Dune::BitSetVector<blocksize>& dirichletNodes,
          double tolerance,
@@ -306,8 +306,8 @@ setup(const GridType& grid,
 }
 
 
-template <class Basis, class TargetSpace>
-void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
 {
     int rank = grid_->comm().rank();
 
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 70c69ffa7a518e134910d8a9a23c439205e7bec1..3184c86f91dfbcbbd56177908ae5366f53a4fb03 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -69,7 +69,7 @@ struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
 };
 
 /** \brief Riemannian trust-region solver for geodesic finite-element problems */
-template <class Basis, class TargetSpace>
+template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
 class RiemannianTrustRegionSolver
     : public IterativeSolver<std::vector<TargetSpace>,
                              Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
@@ -115,7 +115,7 @@ public:
 
     /** \brief Set up the solver using a monotone multigrid method as the inner solver */
     void setup(const GridType& grid,
-               const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
+               const Assembler* assembler,
                const SolutionType& x,
                const Dune::BitSetVector<blocksize>& dirichletNodes,
                double tolerance,
@@ -188,7 +188,7 @@ protected:
     std::unique_ptr<MatrixType> hessianMatrix_;
 
     /** \brief The assembler for the material law */
-    const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
+    const Assembler* assembler_;
 
     /** \brief The solver for the quadratic inner problems */
     std::shared_ptr<Solver> innerSolver_;