diff --git a/demo/finitevolume.py b/demo/finitevolume.py
new file mode 100644
index 0000000000000000000000000000000000000000..d70950527bb5cae971e898284c0eb8d1a4c4d09d
--- /dev/null
+++ b/demo/finitevolume.py
@@ -0,0 +1,129 @@
+import numpy as np
+from mpi4py import MPI
+
+from dune.common import FieldVector
+import dune.grid as grid
+from dune.grid import yaspGrid, gridFunction
+
+# return a vector with the initial concentration by evaluating
+# the grid function c0 at the barycenter of each element using
+# global coordinates
+def initialize(gridView, c0):
+    # mapper with one value on each element
+    mapper = gridView.mapper(lambda gt: gt.dim == gridView.dimension)
+    # numpy vector with data on each element
+    c = np.zeros(len(mapper))
+    for e in gridView.allPartition.elements:
+        center = e.geometry.center      # global coordinate of cell center
+        c[mapper.index(e)] = c0(center) # initialize cell concentration
+    return c, mapper
+
+# main evolution method - the concentration is update in place
+# using 'b' at the boundary and 'u' as transport direction
+# c_T^{n+1} -= dt sum_j |e^T_j| g^T_j(c_T,c_{N_j})
+# where g^T_j is an upwind flux over the edges (e^T_j)_j of the element T
+# using the neighboring concentration value c_{N_j}. The numerical flux is
+# consistent with the linear analytical flux f(c) = uc with a given
+# constant transport direction u. The time step size dt is computed from
+# the maximum wave speed.
+def evolve(gridView, mapper, concentration, t, b, u):
+    # update vector and time step
+    update = np.zeros(len(concentration))
+    dt = 1e100
+
+    # iterate over the elements
+    for entity in gridView.allPartition.elements:
+        volume    = entity.geometry.volume
+        indexi    = mapper.index(entity)
+        sumfactor = 0.0
+
+        # iterate over the intersections
+        for intersection in gridView.intersections(entity):
+            igeo = intersection.geometry
+
+            faceglobal = igeo.center
+
+            normalvelocity = (u(faceglobal, t) * intersection.centerUnitOuterNormal) * igeo.volume
+
+            factor = normalvelocity / volume
+
+            if factor > 0:
+                sumfactor += factor
+
+            outside = intersection.outside
+            if outside is not None:              # inner intersection so use data on neighboring element
+                indexj = mapper.index(outside)
+
+                if indexi < indexj:
+                    nbfactor = normalvelocity / outside.geometry.volume
+
+                    # upwind concentration
+                    c = concentration[indexj if factor < 0 else indexi]
+
+                    update[indexi] -= c * factor
+                    update[indexj] += c * nbfactor
+            elif intersection.boundary:     # boundary intersection
+                c = b(faceglobal, t) if factor < 0 else concentration[indexi]
+                update[indexi] -= c * factor
+        dt = min(dt, 1.0 / sumfactor)
+
+    # reduce time step a bit
+    dt = 0.9 * dt
+
+    concentration += update*dt # concentration vector updated by side effect
+
+    return dt
+
+
+if __name__ == "__main__":
+    domain = grid.cartesianDomain([0, 0], [1, 1], [40, 40])
+    gridView = grid.yaspGrid(domain)
+
+    @gridFunction(gridView)
+    def c0(x):
+        if x.two_norm > 0.125 and x.two_norm < 0.5:
+            return 1.0
+        else:
+            return 0.0
+
+    def u(x, t):
+        return FieldVector([1.0, 1.0])
+
+    def b(x, t):
+        return c0(x-t*u(x,t))
+
+    gridView.hierarchicalGrid.globalRefine(1)
+
+    concentration, mapper = initialize(gridView, c0)
+
+    vtkwriter = gridView.vtkWriter()
+
+    # make a grid function from the piecewise constant data for output
+    @gridFunction(gridView)
+    def solution(element,point):
+        return [concentration[mapper.index(element)]]
+    solution.addToVTKWriter("concentration", vtkwriter, grid.DataType.CellData)
+    # solution.plot()     # uncomment this for matplot output
+
+    counter = 0
+    t, dt = 0.0, 0.0
+    k = 0
+    saveInterval = 0.1
+    saveStep = 0.1
+
+    # output the initial data to a vtk file
+    vtkwriter.write("concentration", counter, grid.OutputType.appendedraw)
+    counter += 1
+
+    # time loop
+    while t < 1.:
+        k += 1
+
+        dt = evolve(gridView, mapper, concentration, t, b, u)
+        print("k=" + str(k) + " ; dt=" + str(dt) + " ; time=" + str(t))
+
+        t += dt
+        if t >= saveStep:
+            vtkwriter.write("concentration", counter, grid.OutputType.appendedraw)
+            saveStep += saveInterval
+            counter += 1