diff --git a/linelast.cc b/linelast.cc
index 572abaf3e396808a520ed581802c036a5d4efaef..b09d547c2fa63113d697ce3a30e31d4353b02c98 100644
--- a/linelast.cc
+++ b/linelast.cc
@@ -16,6 +16,8 @@
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/sampleonbitfield.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
 #include <dune/fufem/boundarypatchprolongator.hh>
 #include <dune/fufem/readbitfield.hh>
 #include <dune/fufem/estimators/fractionalmarking.hh>
@@ -74,6 +76,8 @@ int main (int argc, char *argv[]) try
     std::string parFile             = parameterSet.get("parFile", "");
     std::string dirichletFile       = parameterSet.get<std::string>("dirichletFile");
     std::string dirichletValuesFile = parameterSet.get<std::string>("dirichletValuesFile");
+    std::string neumannFile       = parameterSet.get<std::string>("neumannFile");
+    std::string neumannValuesFile = parameterSet.get<std::string>("neumannValuesFile");
 
     // /////////////////////////////
     //   Generate the grid
@@ -97,7 +101,13 @@ int main (int argc, char *argv[]) try
     VectorType coarseDirichletValues(grid->size(0, dim));
     AmiraMeshReader<GridType>::readFunction(coarseDirichletValues, path + dirichletValuesFile);
 
+    LevelBoundaryPatch coarseNeumannBoundary;
+    coarseNeumannBoundary.setup(grid->levelView(0));
+    readBoundaryPatch<GridType>(coarseNeumannBoundary, path + neumannFile);
 
+    VectorType coarseNeumannValues(grid->size(0, dim));
+    AmiraMeshReader<GridType>::readFunction(coarseNeumannValues, path + neumannFile);
+    P1NodalBasis<GridType::LevelGridView> coarseBasis(grid->levelView(0));
 
     for (int i=0; i<minLevel; i++)
         grid->globalRefine(1);
@@ -119,6 +129,8 @@ int main (int argc, char *argv[]) try
         LeafBoundaryPatch leafDirichletBoundary(grid->leafGridView());
         BoundaryPatchProlongator<GridType>::prolong(coarseDirichletBoundary, leafDirichletBoundary);
         
+        LeafBoundaryPatch leafNeumannBoundary(grid->leafGridView());
+        BoundaryPatchProlongator<GridType>::prolong(coarseNeumannBoundary, leafNeumannBoundary);
 
         BitSetVector<dim> dirichletNodes(grid->size(dim),false);
         for (int i=0; i<grid->size(dim);i++)
@@ -126,6 +138,10 @@ int main (int argc, char *argv[]) try
                 for (int j=0; j<dim; j++)
                     dirichletNodes[i][j]=true;
 
+        VectorType dirichletValues, neumannValues;
+        Functions::interpolate(p1NodalBasis, dirichletValues, Functions::makeFunction(coarseBasis,coarseDirichletValues));
+
+        Functions::interpolate(p1NodalBasis, neumannValues, Functions::makeFunction(coarseBasis,coarseNeumannValues));
 
         VectorType rhs(grid->size(dim));
         rhs = 0;
@@ -143,9 +159,14 @@ int main (int argc, char *argv[]) try
         OperatorType stiffnessMatrix;
                 
         p1Assembler.assemble(p1LocalAssembler, stiffnessMatrix);
-        
-        
-        
+
+        // Assemble Neumann forces
+        BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(p1NodalBasis,neumannBoundary);
+        BasisGridFunction<P1Basis,VectorType> neumannFunction(p1NodalBasis, neumannValues);
+
+        NeumannBoundaryAssembler<GridType, Dune::FieldVector<double,3> > localNeumannAssembler(neumannFunction);
+        boundaryFunctionalAssembler.assemble(localNeumannAssembler, rhs);
+
         // ///////////////////////////
         //   Create a solver
         // ///////////////////////////