diff --git a/src/06-interpolation.cc b/src/06-interpolation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2e4b83edff7b27836b74a4d0a2eeaafff2b072a6
--- /dev/null
+++ b/src/06-interpolation.cc
@@ -0,0 +1,200 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+
+
+// included standard library headers
+#include <iostream>
+#include <array>
+#include <math.h>
+
+// included dune-common headers
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/stringutility.hh>
+
+// included dune-geometry headers
+#include <dune/geometry/quadraturerules.hh>
+
+// included dune-grid headers
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/uggrid.hh>
+
+// included dune-istl headers
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+
+
+// included dune-localfunctions headers
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+// included dune-functions headers
+#include <dune/functions/common/functionfromcallable.hh>
+
+// included dune-fu-tutorial headers
+#include <dune/fu-tutorial/referenceelementutility.hh>
+
+
+// included dune-fu-tutorial headers
+#include "04-gridviews.hh"
+#include "05-poisson-problem.hh"
+
+
+
+template<class GridView, class Vector, class F>
+void interpolateFunction(const GridView& gridView, Vector& v, const F& f)
+{
+ static const int dim = GridView::dimension;
+ const auto& indexSet = gridView.indexSet();
+ std::size_t size = indexSet.size(dim);
+
+ v.resize(size);
+ v = 0;
+
+ Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;
+
+ for(const auto& e: Dune::elements(gridView))
+ {
+ const auto& geometry = e.geometry();
+ const auto& localFiniteElement = feCache.get(e.type());
+ auto localSize = localFiniteElement.localBasis().size();
+
+ using Element = std::decay_t<decltype(e)>;
+ using LocalDomain = typename Element::Geometry::LocalCoordinate;
+ using Range = double;
+
+ auto localF = [&](const LocalDomain& x) {
+ return f(geometry.global(x));
+ };
+
+ std::vector<Range> localV;
+
+// localFiniteElement.localInterpolation().interpolate(localF, localV);
+
+ using FunctionFromCallable = typename Dune::Functions::FunctionFromCallable<Range(LocalDomain), decltype(localF)>;
+ auto localFWithEvaluate = FunctionFromCallable(localF);
+ localFiniteElement.localInterpolation().interpolate(localFWithEvaluate, localV);
+
+ for(std::size_t i=0; i<localSize; ++i)
+ {
+ std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
+ v[globalI] = localV[i];
+ }
+ }
+}
+
+
+
+int main(int argc, char** argv)
+{
+ try{
+ // Maybe initialize MPI
+ Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+
+ // Print process rank
+ if(Dune::MPIHelper::isFake)
+ std::cout<< "This is a sequential program." << std::endl;
+ else
+ std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
+ <<" processes!"<<std::endl;
+
+// using Grid = Dune::YaspGrid<2>;
+// auto gridPointer = createCubeGrid<Grid>();
+
+// using Grid = Dune::YaspGrid<3>;
+// auto gridPointer = createCubeGrid<Grid>();
+
+// using Grid = Dune::UGGrid<2>;
+// auto gridPointer = createCubeGrid<Grid>();
+
+ using Grid = Dune::UGGrid<2>;
+ auto gridPointer = createSimplexGrid<Grid>();
+
+ Grid& grid = *gridPointer;
+
+ grid.globalRefine(5);
+ for(int i=0; i<2; ++i)
+ {
+ auto gridView = grid.levelGridView(grid.maxLevel());
+ for(const auto& e : Dune::elements(gridView))
+ if (e.geometry().center().two_norm() < .5)
+ grid.mark(1, e);
+ grid.adapt();
+ }
+ auto gridView = grid.leafGridView();
+ using GridView = decltype(gridView);
+
+
+ using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
+ using Vector = Dune::BlockVector<Dune::FieldVector<double,1>>;
+
+ Matrix A;
+ Vector rhs;
+ Vector x;
+
+
+ auto rhsFunction = [](auto x) {
+ return (x.two_norm() < .5)*100.0;
+ };
+ assemblePoissonProblem(gridView, A, rhs, rhsFunction);
+
+ x.resize(rhs.size(), 0);
+
+ auto dirichletFunction = [](auto x) {
+ return sin(x[0] * 2.0*M_PI)*.1;
+ };
+ interpolateFunction(gridView, x, dirichletFunction);
+
+ std::vector<bool> isBoundary;
+ isBoundary.resize(rhs.size(), false);
+ computeBoundaryVertices(gridView, isBoundary);
+
+ for(std::size_t i=0; i<rhs.size(); ++i)
+ {
+ if (isBoundary[i])
+ {
+ for(auto& entry: A[i])
+ entry = 0;
+ A[i][i] = 1;
+ rhs[i] = x[i];
+ }
+ }
+
+ Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);
+ Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,1.0);
+ Dune::SeqSSOR<Matrix,Vector,Vector> preconditioner(A,1,1.0);
+ Dune::CGSolver<Vector> cg(op,
+ preconditioner, // preconditione
+ 1e-4, // desired residual reduction factor
+ 200, // maximum number of iterations
+ 2); // verbosity of the solver
+
+ // Object storing some statistics about the solving process
+ Dune::InverseOperatorResult statistics;
+
+ // Solve!
+ cg.apply(x, rhs, statistics);
+
+ Dune::VTKWriter<GridView> vtkWriter(gridView);
+ vtkWriter.addVertexData(x, "solution");
+ vtkWriter.write("06-interpolation");
+
+ return 0;
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ }
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d50032e6f78b2324df4517d536d7ca4fa9d21b1a..067c003e6a4aec0120a46dd8f0145d27de950c21 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,3 +15,6 @@ target_link_dune_default_libraries("04-gridviews")
add_executable("05-poisson-problem" 05-poisson-problem.cc)
target_link_dune_default_libraries("05-poisson-problem")
+
+add_executable("06-interpolation" 06-interpolation.cc)
+target_link_dune_default_libraries("06-interpolation")