diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index 41fc5e9885128430faa47584d603376768dca910..0db5a2c0106cc822dfc646e5cf9842d0f328e25d 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -23,7 +23,6 @@
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/prolongboundarypatch.hh>
 #include <dune/fufem/sampleonbitfield.hh>
-#include <dune/fufem/neumannassembler.hh>
 #include <dune/fufem/computestress.hh>
 
 #include <dune/fufem/functionspacebases/q1nodalbasis.hh>
@@ -392,10 +391,17 @@ int main (int argc, char *argv[]) try
                                               rodX[0].r,
                                               neumannValues);
 
-            rhs3d = 0;
+/*            rhs3d = 0;
             assembleAndAddNeumannTerm<GridType::LevelGridView, VectorType>(interfaceBoundary[complex.continuumGrids_["continuum"]->maxLevel()],
                                                         neumannValues,
-                                                        rhs3d);
+                                                        rhs3d);*/
+            /** \todo The LevelBasis is a hack.  The interfaceBoundary should really by a LeafBoundaryPatch anyways */
+            typedef P1NodalBasis<GridType::LevelGridView,double> LevelBasis;
+            LevelBasis levelBasis(complex.continuumGrids_["continuum"]->levelView(toplevel));
+            BoundaryFunctionalAssembler<LevelBasis> boundaryFunctionalAssembler(levelBasis, interfaceBoundary.back());
+            BasisGridFunction<LevelBasis, VectorType> neumannValuesFunction(levelBasis, neumannValues);
+            NeumannBoundaryAssembler<GridType, FieldVector<double,dim> > localNeumannAssembler(neumannValuesFunction);
+            boundaryFunctionalAssembler.assemble(localNeumannAssembler, rhs3d, true);
 
             // ///////////////////////////////////////////////////////////
             //   Solve the Neumann problem for the continuum
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index 7d0d95f797ace6f40293012ede407ad7b378b431..faffd1756c95c5b58ab21cd423016ab41430095c 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -11,6 +11,9 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+
 #include <dune/gfe/coupling/rodcontinuumcomplex.hh>
 
 template <class GridView, class MatrixType, class VectorType>
@@ -78,9 +81,11 @@ public:
 
         // The weak form of the Neumann data
         VectorType rhs = weakVolumeAndNeumannTerm_;
-        assembleAndAddNeumannTerm<GridView, VectorType>(interface_,
-                                                        neumannValues,
-                                                        rhs);
+
+        BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(basis, interface_);
+        BasisGridFunction<P1Basis,VectorType> neumannValuesFunction(basis,neumannValues);
+        NeumannBoundaryAssembler<typename GridView::Grid, Dune::FieldVector<double,3> > localNeumannAssembler(neumannValuesFunction);
+        boundaryFunctionalAssembler.assemble(localNeumannAssembler, rhs, false);
 
         //   Solve the Neumann problem for the continuum
         VectorType x = dirichletValues_;