diff --git a/notebooks/dunefunctions.ipynb b/notebooks/dunefunctions.ipynb
index cee9febac2d32d43233ff48f28679cded1cf5796..4fd9b37918577595ac2a9482569a05145d5c618b 100644
--- a/notebooks/dunefunctions.ipynb
+++ b/notebooks/dunefunctions.ipynb
@@ -194,6 +194,23 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.5.2"
+  },
   "livereveal": {
    "center": false,
    "controls": false,
diff --git a/notebooks/dunefunctions.py b/notebooks/dunefunctions.py
index 688912e283a158c8f9069451769e16f7b159c503..bb5282f97830044717286255239ea9075952f749 100644
--- a/notebooks/dunefunctions.py
+++ b/notebooks/dunefunctions.py
@@ -3,7 +3,7 @@
 
 # # Integrating dune-functions
 
-# In[1]:
+# In[ ]:
 
 
 import time
@@ -95,6 +95,7 @@ rhs_dofs = numpy.zeros(taylorHoodBasis.dimension)
 
 # In[ ]:
 
+
 from dune.generator import algorithm
 def onBoundary(x):
     return any(x[i] < 1e-10 or x[i] > (1-1e-10) for i in range(len(x)))
diff --git a/notebooks/finiteelements.ipynb b/notebooks/finiteelements.ipynb
index 01b3816affd4752bfd4f69073e1a7bd8add7fb83..05564b206ec4fbec9aa767aa0b1c9742640df396 100644
--- a/notebooks/finiteelements.ipynb
+++ b/notebooks/finiteelements.ipynb
@@ -261,8 +261,7 @@
    },
    "outputs": [],
    "source": [
-    "from dune.plotting import plotPointData\n",
-    "plotPointData(lgf, figsize=(9,9), gridLines=None)"
+    "lgf.plot(figsize=(9,9), gridLines=None)"
    ]
   },
   {
@@ -274,6 +273,23 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.5.2"
+  },
   "livereveal": {
    "center": false,
    "controls": false,
diff --git a/notebooks/finiteelements.py b/notebooks/finiteelements.py
index 9893f618d558f91a4d24118e29dc8a7937e2e517..e1a8092ce71a9d721930e2164690759022b0ddd2 100644
--- a/notebooks/finiteelements.py
+++ b/notebooks/finiteelements.py
@@ -25,7 +25,7 @@ import math
 # In[ ]:
 
 
-from dune.grid import cartesianDomain
+from dune.grid import cartesianDomain, gridFunction
 from dune.alugrid import aluConformGrid
 vertices = numpy.array([(0,0), (1,0), (1,1), (0,1), (-1,1), (-1,0), (-1,-1), (0,-1)])
 triangles = numpy.array([(2,0,1), (0,2,3), (4,0,3), (0,4,5), (6,0,5), (0,6,7)])
@@ -41,7 +41,7 @@ class LinearShapeFunction:
         self.ofs = ofs
         self.grad = grad
     def evaluate(self, local):
-        return self.ofs + sum([x*y for x, y in zip(self.grad, local)])
+        return [self.ofs + sum([x*y for x, y in zip(self.grad, local)])]
     def gradient(self, local):
         return self.grad
 
@@ -62,7 +62,7 @@ for i in range(dim):
 from dune.geometry import quadratureRules
 from dune.istl import blockVector
 
-f = lambda v: sum(2.0 * x * (1 - x) for x in v)
+f = lambda v: [sum(2.0 * x * (1 - x) for x in v)]
 
 dim, indexSet = aluView.dimension, aluView.indexSet
 rhs = blockVector(indexSet.size(dim))
@@ -83,10 +83,10 @@ for e in aluView.elements:
 # In[ ]:
 
 
-from dune.istl import bcrsMatrix, BCRSMatrix11
+from dune.istl import BuildMode, bcrsMatrix, BCRSMatrix
 
 dim, indexSet = aluView.dimension, aluView.indexSet
-matrix = bcrsMatrix((indexSet.size(dim), indexSet.size(dim)), 8, 0.1, BCRSMatrix11.implicit)
+matrix = bcrsMatrix((indexSet.size(dim), indexSet.size(dim)), 8, 0.1, BuildMode.implicit)
 
 quadrature = quadratureRules(1)
 for e in aluView.elements:
@@ -145,13 +145,11 @@ _ = solver(u, rhs)
 
 @gridFunction(aluView)
 def lgf( element, local ):
-    return [ sum( phi.evaluate(local) * u[indexSet.subIndex(element, i, dim)]\
-            for i, phi in enumerate(p1ShapeFunctionSet)) ]
+    return [ sum( phi.evaluate(local) * u[indexSet.subIndex(element, i, dim)]                  for i, phi in enumerate(p1ShapeFunctionSet)) ]
 _ = aluView.writeVTK("fem2d", pointdata={"u": lgf})
 
 
 # In[ ]:
 
 
-from dune.plotting import plotPointData
-plotPointData(lgf, figsize=(9,9), gridLines=None)
+lgf.plot(figsize=(9,9), gridLines=None)
diff --git a/notebooks/finitevolume.ipynb b/notebooks/finitevolume.ipynb
index 2b461b59ae2d5757314ebc12ed6427e3bb146c7f..278c15a138d2574183bfd3edb6cd512a8b4651cb 100644
--- a/notebooks/finitevolume.ipynb
+++ b/notebooks/finitevolume.ipynb
@@ -189,9 +189,8 @@
    },
    "outputs": [],
    "source": [
-    "from dune.plotting import plotPointData\n",
     "cgrid = yaspView.function(lambda e,p: [c[mapper.index(e)]])\n",
-    "plotPointData(cgrid, figsize=(9,9))"
+    "cgrid.plot(figsize=(9,9),on=\"cells\")"
    ]
   },
   {
@@ -234,8 +233,7 @@
    },
    "outputs": [],
    "source": [
-    "from dune.plotting import plotPointData\n",
-    "plotPointData(cgrid, figsize=(9,9))"
+    "cgrid.plot( figsize=(9,9), on=\"cells\")"
    ]
   },
   {
@@ -259,7 +257,7 @@
     "while t < 0.5:\n",
     "    t += evolve(yaspView, mapper, c, b, t)\n",
     "print(\"time used:\", time.time()-start)\n",
-    "plotPointData(cgrid, figsize=(9,9),gridLines=\"\")"
+    "cgrid.plot(figsize=(9,9),gridLines=\"\",on=\"cells\")"
    ]
   },
   {
@@ -271,6 +269,23 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.5.2"
+  },
   "livereveal": {
    "center": false,
    "controls": false,
diff --git a/notebooks/finitevolume.py b/notebooks/finitevolume.py
index 7838c68cda5a6715ecfa1bf28f022fbd3f33dc69..e2c95bb0f369ae03fa64e0c6f232d4be7887a569 100644
--- a/notebooks/finitevolume.py
+++ b/notebooks/finitevolume.py
@@ -110,9 +110,8 @@ print("time used:", time.time()-start)
 # In[ ]:
 
 
-from dune.plotting import plotPointData
 cgrid = yaspView.function(lambda e,p: [c[mapper.index(e)]])
-plotPointData(cgrid, figsize=(9,9))
+cgrid.plot(figsize=(9,9),on="cells")
 
 
 # Let's do this on the C++ side instead:
@@ -136,8 +135,7 @@ print("time used:", time.time()-start)
 # In[ ]:
 
 
-from dune.plotting import plotPointData
-plotPointData(cgrid, figsize=(9,9))
+cgrid.plot( figsize=(9,9), on="cells")
 
 
 # Now that its a bit faster let us recompute the solution on a finer grid:
@@ -153,4 +151,4 @@ t = 0.0
 while t < 0.5:
     t += evolve(yaspView, mapper, c, b, t)
 print("time used:", time.time()-start)
-plotPointData(cgrid, figsize=(9,9),gridLines="")
+cgrid.plot(figsize=(9,9),gridLines="",on="cells")
diff --git a/python/dune/grid/grid_generator.py b/python/dune/grid/grid_generator.py
index d72a64a7cf17c4992cb893ac54b72d5e8a8f6596..82fe64a3d4e6d1087522431475914992e4e6624f 100644
--- a/python/dune/grid/grid_generator.py
+++ b/python/dune/grid/grid_generator.py
@@ -97,9 +97,9 @@ def plot(self, function=None, *args, **kwargs):
     else:
         try:
             grid = function.grid
-            dune.plotting.plotPointData(solution=function,*args,**kwargs)
+            dune.plotting.plot(solution=function,*args,**kwargs)
         except AttributeError:
-            dune.plotting.plotPointData(solution=self.function(function),*args,**kwargs)
+            dune.plotting.plot(solution=self.function(function),*args,**kwargs)
 
 @deprecated("use the `gridFunction` decorator")
 def globalGridFunction(gv, evaluator):
diff --git a/python/dune/plotting.py b/python/dune/plotting.py
index 19f19e3d9d28cefc98519f0bf13f211daeb540a9..75f6cdd5c4a4d4c6a8003ac3d99d4d72c0bd8c70 100644
--- a/python/dune/plotting.py
+++ b/python/dune/plotting.py
@@ -34,15 +34,17 @@ def plotGrid(grid, gridLines="black", figure=None,
     pyplot.show(block=block)
 
 
-def _plotPointData(fig, grid, solution, level=0, gridLines="black",
+def _plotData(fig, grid, solution, level=0, gridLines="black",
         component=None, vectors=None, nofVectors=None,
-        xlim=None, ylim=None, clim=None, cmap=None, colorbar=True):
+        xlim=None, ylim=None, clim=None, cmap=None, colorbar=True,
+        on="cell"):
 
     if (gridLines is not None) and (gridLines != ""):
         _plotGrid(fig, grid, gridLines=gridLines)
 
     if solution is not None:
-        if not any(g.isNone for g in grid.indexSet.types(0)):
+        if on == "points":
+            assert not any(gt.isNone for gt in grid.indexSet.types(0)), "Can't plot point data with polygonal grids, use `on=\"cells\" in plotting command"
             triangulation = grid.triangulation(level)
             data = solution.pointData(level)
             try:
@@ -152,7 +154,7 @@ def plotPointData(solution, level=0, gridLines="black",
         grid = solution
         solution = None
     if not grid.dimension == 2:
-        print("inline plotting so far only available for 2d grids")
+        raise ValueError("inline plotting so far only available for 2d grids")
         return
 
     if figure is None:
@@ -160,7 +162,33 @@ def plotPointData(solution, level=0, gridLines="black",
         show = True
     else:
         show = False
-    _plotPointData(figure,grid,solution,level,gridLines,None,vectors,nofVectors,xlim,ylim,clim,cmap,colorbar=colorbar)
+    _plotData(figure,grid,solution,level,gridLines,None,
+            vectors,nofVectors,xlim,ylim,clim,cmap,
+            colorbar=colorbar,on="points")
+
+    if show:
+        pyplot.show(block=block)
+
+def plotCellData(solution, level=0, gridLines="black",
+        vectors=None, nofVectors=None, figure=None,
+        xlim=None, ylim=None, clim=None, figsize=None, cmap=None,
+        colorbar=True):
+    try:
+        grid = solution.grid
+    except:
+        grid = solution
+        solution = None
+    if not grid.dimension == 2:
+        raise ValueError("inline plotting so far only available for 2d grids")
+        return
+
+    if figure is None:
+        figure = pyplot.figure(figsize=figsize)
+        show = True
+    else:
+        show = False
+    _plotData(figure,grid,solution,level,gridLines,None,vectors,nofVectors,xlim,ylim,clim,cmap,
+            colorbar=colorbar,on="cells")
 
     if show:
         pyplot.show(block=block)
@@ -173,7 +201,7 @@ def plotComponents(solution, level=0, show=None, gridLines="black", figure=None,
         grid = solution
         solution = None
     if not grid.dimension == 2:
-        print("inline plotting so far only available for 2d grids")
+        raise ValueError("inline plotting so far only available for 2d grids")
         return
 
     if not show:
@@ -187,15 +215,33 @@ def plotComponents(solution, level=0, show=None, gridLines="black", figure=None,
     # first the grid if required
     if (gridLines is not None) and (gridLines != ""):
         pyplot.subplot(subfig)
-        _plotPointData(figure,grid,None,level,gridLines,None,False,None,xlim,ylim,clim,cmap)
+        _plotData(figure,grid,None,level,gridLines,None,False,None,xlim,ylim,clim,cmap,
+                on="points")
 
     # add the data
     for p in show:
         pyplot.subplot(subfig+offset+p)
-        _plotPointData(figure,grid,solution,level,"",p,False,None,xlim,ylim,clim,cmap,False)
+        _plotData(figure,grid,solution,level,"",p,False,None,xlim,ylim,clim,cmap,False,
+                on="points")
 
     pyplot.show(block=block)
 
+def plot(solution,*args,**kwargs):
+    try:
+        grid = solution.grid
+    except:
+        grid = solution
+    defaultOn = "cells" if any(gt.isNone for gt in grid.indexSet.types(0)) else "points"
+    use = kwargs.pop("on",defaultOn)
+    if use == "points":
+        plotPointData(solution,*args,**kwargs)
+    elif use == "components-points":
+        plotComponents(solution,*args,**kwargs)
+    elif use == "cells":
+        plotCellData(solution,*args,**kwargs)
+    else:
+        raise ValueError("wrong value for 'on' parameter should be one of 'points','cells','components-points'")
+
 def mayaviPointData(solution, level=0, component=0):
     grid = solution.grid
     from mayavi import mlab