diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 56ad0096c0996a114dfed091d749f12b2f7d441f..7b312f83fff82a5fb4efb641e0ccd928a41b3b82 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -4,6 +4,7 @@ The methods to run the parts of the form compiler
 Should also contain the entrypoint methods.
 """
 from __future__ import absolute_import
+from os.path import splitext
 
 # Disable loopy caching before we do anything else!
 import loopy
@@ -32,69 +33,64 @@ def read_ufl(uflfile):
 
     # Extract the form from the given data
     data = interpret_ufl_namespace(namespace)
-    form = data.forms[0]
-    formdata = compute_form_data(form)
-    form = formdata.preprocessed_form
+    formdatas = []
+    forms = data.forms
+    for index, form in enumerate(forms):
+        formdatas.append(compute_form_data(form))
+        forms[index] = formdatas[index].preprocessed_form
 
-    # We do not expect more than one form
-    assert len(data.forms) == 1
+    # We expect at least one form
+    assert len(data.forms) >= 1
 
     # apply some transformations unconditionally!
     from dune.perftool.ufl.transformations import transform_form
     from dune.perftool.ufl.transformations.indexpushdown import pushdown_indexed
     from dune.perftool.ufl.transformations.reindexing import reindexing
+    for i in range(len(forms)):
+        forms[i] = transform_form(forms[i], pushdown_indexed)
+        forms[i] = transform_form(forms[i], reindexing)
+        formdatas[i].preprocessed_form = forms[i]
 
-    form = transform_form(form, pushdown_indexed)
-    form = transform_form(form, reindexing)
-
-    formdata.preprocessed_form = form
-
-    return formdata, data
-
-
-def generate_driver(form, filename):
-    # The driver module uses a global variable for ease of use
-    from dune.perftool.pdelab.driver import set_form
-    set_form(form)
-
-    # The vtkoutput is the generating method that triggers all others.
-    # Alternatively, one could use the `solve` method.
-    from dune.perftool.pdelab.driver import vtkoutput
-    vtkoutput()
-
-    from dune.perftool.generation import retrieve_cache_items
-    from cgen import FunctionDeclaration, FunctionBody, Block, Value
-    driver_signature = FunctionDeclaration(Value('void', 'driver'), [Value('int', 'argc'), Value('char**', 'argv')])
-    driver_body = Block(contents=[i for i in retrieve_cache_items("driver and preamble", make_generable=True)])
-    driver = FunctionBody(driver_signature, driver_body)
-
-    from dune.perftool.file import generate_file
-    generate_file(filename, "driver", [driver])
-
-    # Reset the caching data structure
-    from dune.perftool.generation import delete_cache_items
-    delete_cache_items()
+    return formdatas, data
 
 
 # This function is the entrypoint of the ufl2pdelab executable
 def compile_form():
     from dune.perftool.options import get_option
-    formdata, data = read_ufl(get_option("uflfile"))
+    formdatas, data = read_ufl(get_option("uflfile"))
 
     from dune.perftool.generation import cache_context, global_context
-    with global_context(formdata=formdata, data=data):
+    with global_context(data=data, formdatas=formdatas):
+        # Generate driver file
         if get_option("driver_file"):
             with cache_context('driver', delete=True):
-                generate_driver(formdata.preprocessed_form, get_option("driver_file"))
-
-        if get_option("operator_file"):
-            from dune.perftool.pdelab.localoperator import generate_localoperator_kernels
-            kernels = generate_localoperator_kernels(formdata, data)
-
-        # TODO insert sophisticated analysis/feedback loops here
-        if get_option("interactive"):
-            from dune.perftool.interactive import start_interactive_session
-            start_interactive_session(kernels)
-
-        from dune.perftool.pdelab.localoperator import generate_localoperator_file
-        generate_localoperator_file(kernels)
+                from dune.perftool.pdelab.driver import generate_driver
+                generate_driver(formdatas, data)
+
+    # In case of multiple forms: Genarate one file that includes all localoperator files
+    if len(formdatas) > 1:
+        from dune.perftool.pdelab.localoperator import generate_localoperator_basefile
+        generate_localoperator_basefile(formdatas, data)
+
+    # Generate local operator files
+    for formdata in formdatas:
+        with global_context(data=data, formdata=formdata):
+            # Make sure cache is empty
+            from dune.perftool.generation import delete_cache_items
+            delete_cache_items()
+
+            if get_option("operator_file"):
+                from dune.perftool.pdelab.localoperator import generate_localoperator_kernels
+                kernels = generate_localoperator_kernels(formdata, data)
+
+            # TODO insert sophisticated analysis/feedback loops here
+            if get_option("interactive"):
+                from dune.perftool.interactive import start_interactive_session
+                start_interactive_session(kernels)
+
+            if get_option("operator_file"):
+                from dune.perftool.pdelab.localoperator import name_localoperator_file
+                filename = name_localoperator_file(formdata, data)
+
+                from dune.perftool.pdelab.localoperator import generate_localoperator_file
+                generate_localoperator_file(kernels, filename)
diff --git a/python/dune/perftool/generation/cache.py b/python/dune/perftool/generation/cache.py
index ba5e357d91f358cda2705527610013573174d1e3..c4255edf0873c67956be787c2219446c987deb94 100644
--- a/python/dune/perftool/generation/cache.py
+++ b/python/dune/perftool/generation/cache.py
@@ -196,6 +196,18 @@ class _ConditionDict(dict):
             self.tags = tags
 
         def __getitem__(self, i):
+            # If we do not add these special cases the dictionary will return False
+            # when we execute the following code:
+            #
+            # eval ("True", _ConditionDict(v.tags)
+            #
+            # But in this case we want to return True! A normal dictionary would not attempt
+            # to replace "True" if "True" is not a key. The _ConditionDict obviously has no
+            # such concerns ;).
+            if i == "True":
+                return True
+            if i == "False":
+                return False
             return i in self.tags
 
 
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index bbea79cfe4e129bc19fc3e96563dd2c34cc2ff43..d6761605aaaaf6527dc2a08be44267a317201fbd 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -14,6 +14,7 @@ def get_form_compiler_arguments():
     parser.add_argument('--version', action='version', version='%(prog)s 0.1')
     parser.add_argument("--numerical-jacobian", action="store_true", help="use numerical jacobians (only makes sense, if uflpdelab for some reason fails to generate analytic jacobians)")
     parser.add_argument("--matrix-free", action="store_true", help="use iterative solver with matrix free jacobian application")
+    parser.add_argument("--explicit-time-stepping", action="store_true", help="use explicit time stepping")
     parser.add_argument(
         "--exact-solution-expression",
         type=str,
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index ff87bd114c7ee54107f25ece4e5d0d83531e3c87..11c66909900f0189e700d0a24148864585dd1e37 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -13,6 +13,40 @@ from dune.perftool.generation import (generator_factory,
                                       )
 from dune.perftool.options import get_option
 
+# Have a global variable with the entire form data. This allows functions that depend
+# deterministically on the entire data set to directly access it instead of passing it
+# through the entire generator chain.
+_driver_data = {}
+
+
+# Have a function access this global data structure
+def set_driver_data(formdatas, data):
+    assert (len(formdatas) <= 2)
+    if len(formdatas) == 1:
+        _driver_data['form'] = formdatas[0].preprocessed_form
+        _driver_data['formdata'] = formdatas[0]
+    else:
+        mass_index = mass_form_index(formdatas, data)
+        if mass_index is None:
+            raise NotImplementedError("Form for mass matrix needs to have name 'mass' in ufl file.")
+        _driver_data['mass_form'] = formdatas[mass_index].preprocessed_form
+        _driver_data['mass_formdata'] = formdatas[mass_index]
+        _driver_data['form'] = formdatas[1 - mass_index].preprocessed_form
+        _driver_data['formdata'] = formdatas[1 - mass_index]
+
+    _driver_data['data'] = data
+
+
+def is_stationary():
+    return 'mass_form' not in _driver_data
+
+
+def form_name_suffix(name, formdata):
+    from dune.perftool.pdelab.localoperator import name_form
+    data = _driver_data['data']
+    form_name = name_form(formdata, data)
+    return name + '_' + form_name
+
 
 def has_constraints(element):
     from ufl import MixedElement, VectorElement
@@ -39,16 +73,13 @@ fem_metadata_dependent_symbol = generator_factory(item_tags=("symbol",),
                                                   )
 
 
-# Have a global variable with the entire form data. This allows functions that depend
-# deterministically on the entire data set to directly access it instead of passing it
-# through the entire generator chain.
-_form = None
-
-
-# Have a function access this global data structure
-def set_form(f):
-    global _form
-    _form = f
+def mass_form_index(formdatas, data):
+    for index, formdata in enumerate(formdatas):
+        try:
+            if data.object_names[id(formdata.original_form)] == 'mass':
+                return index
+        except KeyError:
+            continue
 
 
 def is_linear(form):
@@ -134,7 +165,7 @@ def name_initree():
 
 @preamble
 def define_dimension(name):
-    return "static const int {} = {};".format(name, _form.ufl_cell().geometric_dimension())
+    return "static const int {} = {};".format(name, _driver_data['form'].ufl_cell().geometric_dimension())
 
 
 @symbol
@@ -146,11 +177,11 @@ def name_dimension():
 @preamble
 def typedef_grid(name):
     dim = name_dimension()
-    if any(_form.ufl_cell().cellname() in x for x in ["vertex", "interval", "quadrilateral", "hexahedron"]):
+    if any(_driver_data['form'].ufl_cell().cellname() in x for x in ["vertex", "interval", "quadrilateral", "hexahedron"]):
         gridt = "Dune::YaspGrid<{}>".format(dim)
         include_file("dune/grid/yaspgrid.hh", filetag="driver")
     else:
-        if any(_form.ufl_cell().cellname() in x for x in ["triangle", "tetrahedron"]):
+        if any(_driver_data['form'].ufl_cell().cellname() in x for x in ["triangle", "tetrahedron"]):
             # gridt = "Dune::UGGrid<{}>".format(dim)
             # include_file("dune/grid/uggrid.hh", filetag="driver")
             gridt = "Dune::ALUGrid<{}, {}, Dune::simplex, Dune::conforming>".format(dim, dim)
@@ -537,11 +568,11 @@ def name_gfs(expr):
 @preamble
 def define_dofestimate(name):
     # Provide a worstcase estimate for the number of entries per row based on the given gridfunction space and cell geometry
-    if isQuadrilateral(_form.coefficients()[0].ufl_element()):
+    if isQuadrilateral(_driver_data['form'].coefficients()[0].ufl_element()):
         geo_factor = "4"
     else:
         geo_factor = "6"
-    gfs = name_gfs(_form.coefficients()[0].ufl_element())
+    gfs = name_gfs(_driver_data['form'].coefficients()[0].ufl_element())
     ini = name_initree()
     return ["int generic_dof_estimate =  {} * {}.maxLocalSize();".format(geo_factor, gfs),
             "int {} = {}.get<int>(\"istl.number_of_nnz\", generic_dof_estimate);".format(name, ini)]
@@ -579,61 +610,72 @@ def name_matrixbackend():
 
 
 @symbol
-def type_parameters():
-    return "LocalOperatorParams"
+def type_parameters(formdata):
+    from dune.perftool.generation import get_global_context_value
+    data = get_global_context_value("data")
+    from dune.perftool.pdelab.parameter import parameterclass_basename
+    name = parameterclass_basename(formdata, data)
+    return name
 
 
 @preamble
-def define_parameters(name):
-    partype = type_parameters()
+def define_parameters(name, formdata):
+    partype = type_parameters(formdata)
     return "{} {};".format(partype, name)
 
 
 @symbol
-def name_parameters():
-    define_parameters("params")
-    return "params"
+def name_parameters(formdata):
+    name = form_name_suffix("params", formdata)
+    define_parameters(name, formdata)
+    return name
 
 
 @preamble
-def typedef_localoperator(name):
+def typedef_localoperator(name, formdata):
     # No Parameter class here, yet
     # params = type_parameters()
     # return "typedef LocalOperator<{}> {};".format(params, name)
-    ugfs = type_gfs(_form.coefficients()[0].ufl_element())
-    vgfs = type_gfs(_form.arguments()[0].ufl_element())
-    from dune.perftool.options import get_option
-    include_file(get_option('operator_file'), filetag="driver")
-    return "using {} = LocalOperator<{}, {}>;".format(name, ugfs, vgfs)
+    ugfs = type_gfs(_driver_data['form'].coefficients()[0].ufl_element())
+    vgfs = type_gfs(_driver_data['form'].arguments()[0].ufl_element())
+    from dune.perftool.generation import get_global_context_value
+    data = get_global_context_value("data")
+    filename = get_option("operator_file")
+    include_file(filename, filetag="driver")
+    from dune.perftool.pdelab.localoperator import localoperator_basename
+    lopname = localoperator_basename(formdata, data)
+    return "using {} = {}<{}, {}>;".format(name, lopname, ugfs, vgfs)
 
 
 @symbol
-def type_localoperator():
-    typedef_localoperator("LOP")
-    return "LOP"
+def type_localoperator(formdata):
+    name = form_name_suffix("LOP", formdata).upper()
+    typedef_localoperator(name, formdata)
+    return name
 
 
 @preamble
-def define_localoperator(name):
-    loptype = type_localoperator()
+def define_localoperator(name, formdata):
+    loptype = type_localoperator(formdata)
     ini = name_initree()
-    params = name_parameters()
+    params = name_parameters(formdata)
     return "{} {}({}, {});".format(loptype, name, ini, params)
 
 
 @symbol
-def name_localoperator():
-    define_localoperator("lop")
-    return "lop"
+def name_localoperator(formdata):
+    name = form_name_suffix("lop", formdata)
+    define_localoperator(name, formdata)
+    return name
 
 
 @preamble
-def typedef_gridoperator(name):
-    ugfs = type_gfs(_form.coefficients()[0].ufl_element())
-    vgfs = type_gfs(_form.arguments()[0].ufl_element())
-    lop = type_localoperator()
-    ucc = type_constraintscontainer(_form.coefficients()[0].ufl_element())
-    vcc = type_constraintscontainer(_form.arguments()[0].ufl_element())
+def typedef_gridoperator(name, formdata):
+    ugfs = type_gfs(_driver_data['form'].coefficients()[0].ufl_element())
+    vgfs = type_gfs(_driver_data['form'].arguments()[0].ufl_element())
+    lop = type_localoperator(formdata)
+    ucc = type_constraintscontainer(_driver_data['form'].coefficients()[0].ufl_element())
+    vcc = type_constraintscontainer(_driver_data['form'].arguments()[0].ufl_element())
     mb = type_matrixbackend()
     df = type_domainfield()
     r = type_range()
@@ -642,19 +684,21 @@ def typedef_gridoperator(name):
 
 
 @symbol
-def type_gridoperator():
-    typedef_gridoperator("GO")
-    return "GO"
+def type_gridoperator(formdata):
+    name = form_name_suffix("GO", formdata).upper()
+    typedef_gridoperator(name, formdata)
+    return name
 
 
 @preamble
-def define_gridoperator(name):
-    gotype = type_gridoperator()
-    ugfs = name_gfs(_form.coefficients()[0].ufl_element())
-    ucc = name_assembled_constraints(_form.coefficients()[0].ufl_element())
-    vgfs = name_gfs(_form.arguments()[0].ufl_element())
-    vcc = name_assembled_constraints(_form.arguments()[0].ufl_element())
-    lop = name_localoperator()
+def define_gridoperator(name, formdata):
+    gotype = type_gridoperator(formdata)
+    # TODO use formdata insteat of _driver_data object
+    ugfs = name_gfs(_driver_data['form'].coefficients()[0].ufl_element())
+    ucc = name_assembled_constraints(_driver_data['form'].coefficients()[0].ufl_element())
+    vgfs = name_gfs(_driver_data['form'].arguments()[0].ufl_element())
+    vcc = name_assembled_constraints(_driver_data['form'].arguments()[0].ufl_element())
+    lop = name_localoperator(formdata)
     mb = name_matrixbackend()
     return ["{} {}({}, {}, {}, {}, {}, {});".format(gotype, name, ugfs, ucc, vgfs, vcc, lop, mb),
             "std::cout << \"gfs with \" << {}.size() << \" dofs generated  \"<< std::endl;".format(ugfs),
@@ -662,27 +706,29 @@ def define_gridoperator(name):
 
 
 @symbol
-def name_gridoperator():
-    define_gridoperator("go")
-    return "go"
+def name_gridoperator(formdata):
+    name = form_name_suffix("go", formdata)
+    define_gridoperator(name, formdata)
+    return name
 
 
 @preamble
-def typedef_vector(name):
-    gotype = type_gridoperator()
-    return "typedef {}::Traits::Domain {};".format(gotype, name)
+def typedef_vector(name, formdata):
+    go_type = type_gridoperator(formdata)
+    return "typedef {}::Traits::Domain {};".format(go_type, name)
 
 
 @symbol
-def type_vector():
-    typedef_vector("V")
-    return "V"
+def type_vector(formdata):
+    name = form_name_suffix("V", formdata).upper()
+    typedef_vector(name, formdata)
+    return name
 
 
 @preamble
-def define_vector(name):
-    vtype = type_vector()
-    gfs = name_gfs(_form.coefficients()[0].ufl_element())
+def define_vector(name, formdata):
+    vtype = type_vector(formdata)
+    gfs = name_gfs(_driver_data['form'].coefficients()[0].ufl_element())
     return ["{} {}({});".format(vtype, name, gfs), "{} = 0.0;".format(name)]
 
 
@@ -707,10 +753,17 @@ def define_boundary_function(boundary, name):
     gv = name_leafview()
     lambdaname = name_boundary_lambda(boundary, name)
     include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
-    return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
-                                                                                  gv,
-                                                                                  lambdaname,
-                                                                                  )
+    if is_stationary():
+        return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
+                                                                                      gv,
+                                                                                      lambdaname,
+                                                                                      )
+    else:
+        params = name_parameters(_driver_data['formdata'])
+        return "auto {} = Dune::PDELab::makeInstationaryGridFunctionFromCallable({}, {}, {});".format(name,
+                                                                                                      gv,
+                                                                                                      lambdaname,
+                                                                                                      params)
 
 
 @preamble
@@ -760,9 +813,9 @@ def name_solution_function(expr):
 
 
 @preamble
-def interpolate_vector(name):
-    define_vector(name)
-    element = _form.coefficients()[0].ufl_element()
+def interpolate_vector(name, formdata):
+    define_vector(name, formdata)
+    element = _driver_data['form'].coefficients()[0].ufl_element()
     bf = name_boundary_function(element)
     gfs = name_gfs(element)
     return "Dune::PDELab::interpolate({}, {}, {});".format(bf,
@@ -773,8 +826,9 @@ def interpolate_vector(name):
 
 @preamble
 def interpolate_solution_expression(name):
-    define_vector(name)
-    element = _form.coefficients()[0].ufl_element()
+    formdata = _driver_data['formdata']
+    define_vector(name, formdata)
+    element = _driver_data['form'].coefficients()[0].ufl_element()
     sol = name_solution_function(element)
     gfs = name_gfs(element)
     return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
@@ -783,18 +837,19 @@ def interpolate_solution_expression(name):
                                                            )
 
 
-def maybe_interpolate_vector(name):
-    element = _form.coefficients()[0].ufl_element()
+def maybe_interpolate_vector(name, formdata):
+    element = _driver_data['form'].coefficients()[0].ufl_element()
     if has_constraints(element):
-        interpolate_vector(name)
+        interpolate_vector(name, formdata)
     else:
-        define_vector(name)
+        define_vector(name, formdata)
 
 
 @symbol
-def name_vector():
-    maybe_interpolate_vector("x")
-    return "x"
+def name_vector(formdata):
+    name = form_name_suffix("x", formdata)
+    maybe_interpolate_vector(name, formdata)
+    return name
 
 
 @symbol
@@ -842,82 +897,214 @@ def name_reduction():
 @preamble
 def typedef_stationarylinearproblemsolver(name):
     include_file("dune/pdelab/stationary/linearproblem.hh", filetag="driver")
-    gotype = type_gridoperator()
+    gotype = type_gridoperator(_driver_data['formdata'])
     lstype = type_linearsolver()
-    xtype = type_vector()
+    xtype = type_vector(_driver_data['formdata'])
     return "typedef Dune::PDELab::StationaryLinearProblemSolver<{}, {}, {}> {};".format(gotype, lstype, xtype, name)
 
 
-@preamble
-def typedef_stationarynonlinearproblemsolver(name):
-    include_file("dune/pdelab/newton/newton.hh", filetag="driver")
-    gotype = type_gridoperator()
-    lstype = type_linearsolver()
-    xtype = type_vector()
-    return "typedef Dune::PDELab::Newton<{}, {}, {}> {};".format(gotype, lstype, xtype, name)
-
-
 @symbol
 def type_stationarylinearproblemsolver():
     typedef_stationarylinearproblemsolver("SLP")
     return "SLP"
 
 
-@symbol
-def type_stationarynonlinearproblemssolver():
-    typedef_stationarynonlinearproblemsolver("SNP")
-    return "SNP"
-
-
 @preamble
 def define_stationarylinearproblemsolver(name):
     slptype = type_stationarylinearproblemsolver()
-    go = name_gridoperator()
+    formdata = _driver_data['formdata']
+    go = name_gridoperator(formdata)
     ls = name_linearsolver()
-    x = name_vector()
+    x = name_vector(formdata)
     red = name_reduction()
     return "{} {}({}, {}, {}, {});".format(slptype, name, go, ls, x, red)
 
 
-@preamble
-def define_stationarynonlinearproblemsolver(name):
-    snptype = type_stationarynonlinearproblemssolver()
-    go = name_gridoperator()
-    x = name_vector()
-    ls = name_linearsolver()
-    return "{}  {}({}, {}, {});".format(snptype, name, go, x, ls)
-
-
 @symbol
 def name_stationarylinearproblemsolver():
     define_stationarylinearproblemsolver("slp")
     return "slp"
 
 
+@preamble
+def typedef_stationarynonlinearproblemsolver(name, go_type):
+    include_file("dune/pdelab/newton/newton.hh", filetag="driver")
+    ls_type = type_linearsolver()
+    x_type = type_vector(_driver_data['formdata'])
+    return "typedef Dune::PDELab::Newton<{}, {}, {}> {};".format(go_type, ls_type, x_type, name)
+
+
 @symbol
-def name_stationarynonlinearproblemsolver():
-    define_stationarynonlinearproblemsolver("snp")
+def type_stationarynonlinearproblemssolver(go_type):
+    typedef_stationarynonlinearproblemsolver("SNP", go_type)
+    return "SNP"
+
+
+@preamble
+def define_stationarynonlinearproblemsolver(name, go_type, go):
+    snptype = type_stationarynonlinearproblemssolver(go_type)
+    x = name_vector(_driver_data['formdata'])
+    ls = name_linearsolver()
+    return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
+
+
+@symbol
+def name_stationarynonlinearproblemsolver(go_type, go):
+    define_stationarynonlinearproblemsolver("snp", go_type, go)
     return "snp"
 
 
+@preamble
+def typedef_timesteppingmethod(name):
+    r_type = type_range()
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        return "typedef Dune::PDELab::ExplicitEulerParameter<{}> {};".format(r_type, name)
+    else:
+        return "typedef Dune::PDELab::OneStepThetaParameter<{}> {};".format(r_type, name)
+
+
+@symbol
+def type_timesteppingmethod():
+    typedef_timesteppingmethod("TSM")
+    return "TSM"
+
+
+@preamble
+def define_timesteppingmethod(name):
+    tsm_type = type_timesteppingmethod()
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        return "{} {};".format(tsm_type, name)
+    else:
+        return "{} {}(1.0);".format(tsm_type, name)
+
+
+@symbol
+def name_timesteppingmethod():
+    define_timesteppingmethod("tsm")
+    return "tsm"
+
+
+@preamble
+def typedef_instationarygridoperator(name):
+    include_file("dune/pdelab/gridoperator/onestep.hh", filetag="driver")
+    go_type = type_gridoperator(_driver_data['formdata'])
+    mass_go_type = type_gridoperator(_driver_data['mass_formdata'])
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        return "typedef Dune::PDELab::OneStepGridOperator<{},{},false> {};".format(go_type, mass_go_type, name)
+    else:
+        return "typedef Dune::PDELab::OneStepGridOperator<{},{}> {};".format(go_type, mass_go_type, name)
+
+
+@symbol
+def type_instationarygridoperator():
+    typedef_instationarygridoperator("IGO")
+    return "IGO"
+
+
+@preamble
+def define_instationarygridoperator(name):
+    igo_type = type_instationarygridoperator()
+    go = name_gridoperator(_driver_data['formdata'])
+    mass_go = name_gridoperator(_driver_data['mass_formdata'])
+    return "{} {}({}, {});".format(igo_type, name, go, mass_go)
+
+
+@symbol
+def name_instationarygridoperator():
+    define_instationarygridoperator("igo")
+    return "igo"
+
+
+@preamble
+def typedef_onestepmethod(name):
+    r_type = type_range()
+    igo_type = type_instationarygridoperator()
+    snp_type = type_stationarynonlinearproblemssolver(igo_type)
+    vector_type = type_vector(_driver_data['formdata'])
+    return "typedef Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}> {};".format(r_type, igo_type, snp_type, vector_type, vector_type, name)
+
+
+@symbol
+def type_onestepmethod():
+    typedef_onestepmethod("OSM")
+    return "OSM"
+
+
+@preamble
+def define_onestepmethod(name):
+    ilptype = type_onestepmethod()
+    tsm = name_timesteppingmethod()
+    igo_type = type_instationarygridoperator()
+    igo = name_instationarygridoperator()
+    snp = name_stationarynonlinearproblemsolver(igo_type, igo)
+    return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)
+
+
+@symbol
+def name_onestepmethod():
+    define_onestepmethod("osm")
+    return "osm"
+
+
+@preamble
+def typedef_explicitonestepmethod(name):
+    r_type = type_range()
+    igo_type = type_instationarygridoperator()
+    ls_type = type_linearsolver()
+    vector_type = type_vector(_driver_data['formdata'])
+    return "typedef Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}> {};".format(r_type, igo_type, ls_type, vector_type, name)
+
+
+@symbol
+def type_explicitonestepmethod():
+    typedef_explicitonestepmethod("EOSM")
+    return "EOSM"
+
+
+@preamble
+def define_explicitonestepmethod(name):
+    eosm_type = type_explicitonestepmethod()
+    tsm = name_timesteppingmethod()
+    igo = name_instationarygridoperator()
+    ls = name_linearsolver()
+    return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)
+
+
+@symbol
+def name_explicitonestepmethod():
+    define_explicitonestepmethod("eosm")
+    return "eosm"
+
+
 @preamble
 def dune_solve():
     # Test if form is linear in ansatzfunction
-    if is_linear(_form):
-        if get_option("matrix_free"):
-            go = name_gridoperator()
-            x = name_vector()
-            include_file("dune/perftool/matrixfree.hh", filetag="driver")
-            return "solveMatrixFree({},{});".format(go, x)
-        else:
-            slp = name_stationarylinearproblemsolver()
-            return "{}.apply();".format(slp)
+    linear = is_linear(_driver_data['form'])
+
+    # Test wether we want to do matrix free operator evaluation
+    matrix_free = get_option('matrix_free')
+
+    if linear and matrix_free:
+        formdata = _driver_data['formdata']
+        go = name_gridoperator(formdata)
+        x = name_vector(formdata)
+        include_file("dune/perftool/matrixfree.hh", filetag="driver")
+        return "solveMatrixFree({},{});".format(go, x)
+    elif linear and not matrix_free:
+        slp = name_stationarylinearproblemsolver()
+        return "{}.apply();".format(slp)
+    elif not linear and matrix_free:
+        raise NotImplementedError("Nonlinear and matrix free is not yet implemented")
+    elif not linear and not matrix_free:
+        go_type = type_gridoperator(_driver_data['formdata'])
+        go = name_gridoperator(_driver_data['formdata'])
+        snp = name_stationarynonlinearproblemsolver(go_type, go)
+        return "{}.apply();".format(snp)
     else:
-        if get_option("matrix_free"):
-            raise NotImplementedError("Nonlinear and matrix free is not yet implemented")
-        else:
-            snp = name_stationarynonlinearproblemsolver()
-            return "{}.apply();".format(snp)
+        assert False
 
 
 @preamble
@@ -935,7 +1122,7 @@ def name_vtkfile():
 
 @preamble
 def compare_dofs():
-    v = name_vector()
+    v = name_vector(_driver_data['formdata'])
     solution = name_vector_from_solution_expression()
     return ["",
             "// Maximum norm of difference between dof vectors of the",
@@ -965,8 +1152,9 @@ def name_discrete_grid_function(gfs, vector_name):
 
 @preamble
 def compare_L2_squared():
-    element = _form.coefficients()[0].ufl_element()
-    v = name_vector()
+    element = _driver_data['form'].coefficients()[0].ufl_element()
+    formdata = _driver_data['formdata']
+    v = name_vector(formdata)
     gfs = name_gfs(element)
     vdgf = name_discrete_grid_function(gfs, v)
     solution_function = name_solution_function(element)
@@ -989,9 +1177,10 @@ def compare_L2_squared():
 @preamble
 def print_residual():
     ini = name_initree()
-    n_go = name_gridoperator()
-    v = name_vector()
-    t_v = type_vector()
+    formdata = _driver_data['formdata']
+    n_go = name_gridoperator(formdata)
+    v = name_vector(formdata)
+    t_v = type_vector(formdata)
 
     return ["if ({}.get<bool>(\"printresidual\", false)) {{".format(ini),
             "  using Dune::PDELab::Backend::native;",
@@ -1004,11 +1193,12 @@ def print_residual():
 
 @preamble
 def print_matrix():
+    formdata = _driver_data['formdata']
     ini = name_initree()
-    t_go = type_gridoperator()
-    n_go = name_gridoperator()
-    v = name_vector()
-    t_v = type_vector()
+    t_go = type_gridoperator(formdata)
+    n_go = name_gridoperator(formdata)
+    v = name_vector(formdata)
+    t_v = type_vector(formdata)
 
     return ["if ({}.get<bool>(\"printmatrix\", false)) {{".format(ini),
             "  typedef typename {}::Traits::Jacobian M;".format(t_go),
@@ -1054,12 +1244,12 @@ def name_predicate():
 
 @preamble
 def vtkoutput():
-    element = _form.coefficients()[0].ufl_element()
+    element = _driver_data['form'].coefficients()[0].ufl_element()
     define_gfs_name(element)
     include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
     vtkwriter = name_vtkwriter()
     gfs = name_gfs(element)
-    vec = name_vector()
+    vec = name_vector(_driver_data['formdata'])
     vtkfile = name_vtkfile()
     predicate = name_predicate()
     dune_solve()
@@ -1068,7 +1258,7 @@ def vtkoutput():
 
     if get_option("exact_solution_expression"):
         from ufl import MixedElement, VectorElement, TensorElement
-        if isinstance(_form.coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
+        if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
             raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
 
         if get_option("compare_dofs"):
@@ -1078,3 +1268,159 @@ def vtkoutput():
 
     return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate),
             "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]
+
+
+@preamble
+def define_time(name):
+    return "double {} = 0.0;".format(name)
+
+
+@symbol
+def name_time():
+    define_time("time")
+    return "time"
+
+
+@preamble
+def typedef_vtk_sequence_writer(name):
+    include_file("dune/grid/io/file/vtk/vtksequencewriter.hh", filetag="driver")
+    gv_type = type_leafview()
+    return "typedef Dune::VTKSequenceWriter<{}> {};".format(gv_type, name)
+
+
+@symbol
+def type_vtk_sequence_writer():
+    typedef_vtk_sequence_writer("VTKSW")
+    return "VTKSW"
+
+
+@preamble
+def define_vtk_sequence_writer(name):
+    vtksw_type = type_vtk_sequence_writer()
+    vtkw_type = type_vtkwriter()
+    vtkw = name_vtkwriter()
+    vtkfile = name_vtkfile()
+    return "{} {}(std::make_shared<{}>({}), {});".format(vtksw_type, name, vtkw_type, vtkw, vtkfile)
+
+
+@symbol
+def name_vtk_sequence_writer():
+    define_vtk_sequence_writer("vtkSequenceWriter")
+    return "vtkSequenceWriter"
+
+
+@preamble
+def visualize_initial_condition():
+    include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
+    vtkwriter = name_vtk_sequence_writer()
+    element = _driver_data['form'].coefficients()[0].ufl_element()
+    define_gfs_name(element)
+    gfs = name_gfs(element)
+    vector = name_vector(_driver_data['formdata'])
+    predicate = name_predicate()
+    time = name_time()
+    return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vector, predicate),
+            "{}.write({}, Dune::VTK::appendedraw);".format(vtkwriter, time)]
+
+
+@preamble
+def implicit_time_loop():
+    ini = name_initree()
+    formdata = _driver_data['formdata']
+    params = name_parameters(formdata)
+    time = name_time()
+    expr = _driver_data['form'].coefficients()[0].ufl_element()
+    bctype = name_bctype_function(expr)
+    gfs = name_gfs(expr)
+    cc = name_constraintscontainer(expr)
+    vector_type = type_vector(formdata)
+    vector = name_vector(formdata)
+
+    # TODO -> formcompiler argument
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        osm = name_explicitonestepmethod()
+        apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
+    else:
+        boundary = name_boundary_function(expr)
+        osm = name_onestepmethod()
+        apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
+    visualize_initial_condition()
+    vtk_sequence_writer = name_vtk_sequence_writer()
+
+    return ["",
+            "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
+            "double dt = {}.get<double>(\"instat.dt\", 0.1);".format(ini),
+            "while (time<T-1e-8){",
+            "  // Assemble constraints for new time step",
+            "  {}.setTime({}+dt);".format(params, time),
+            "  Dune::PDELab::constraints({}, {}, {});".format(bctype, gfs, cc),
+            "",
+            "  // Do time step",
+            "  {} {}new({});".format(vector_type, vector, vector),
+            "  {}".format(apply_call),
+            "",
+            "  // Accept new time step",
+            "  {} = {}new;".format(vector, vector),
+            "  time += dt;",
+            "",
+            "  // Output to VTK File",
+            "  {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
+            "}",
+            ""]
+
+
+def time_loop():
+    # Test if form is linear in ansatzfunction
+    linear = is_linear(_driver_data['form'])
+
+    # Test wether we want to do matrix free operator evaluation
+    matrix_free = get_option('matrix_free')
+
+    # Create time loop
+    if linear and matrix_free:
+        assert False
+    elif linear and not matrix_free:
+        implicit_time_loop()
+    if not linear and matrix_free:
+        assert False
+    elif not linear and not matrix_free:
+        assert False
+
+    print_residual()
+    print_matrix()
+    if get_option("exact_solution_expression"):
+        from ufl import MixedElement, VectorElement, TensorElement
+        if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
+            raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
+
+        if get_option("compare_dofs"):
+            compare_dofs()
+        if get_option("compare_l2errorsquared"):
+            compare_L2_squared()
+
+def generate_driver(formdatas, data):
+    # The driver module uses a global dictionary for storing necessary data
+    set_driver_data(formdatas, data)
+
+    # The vtkoutput is the generating method that triggers all others.
+    # Alternatively, one could use the `solve` method.
+    if is_stationary():
+        vtkoutput()
+    else:
+        time_loop()
+
+    from dune.perftool.generation import retrieve_cache_items
+    from cgen import FunctionDeclaration, FunctionBody, Block, Value
+    driver_signature = FunctionDeclaration(Value('void', 'driver'), [Value('int', 'argc'), Value('char**', 'argv')])
+    driver_body = Block(contents=[i for i in retrieve_cache_items("driver and preamble", make_generable=True)])
+    driver = FunctionBody(driver_signature, driver_body)
+
+    filename = get_option("driver_file")
+
+    from dune.perftool.file import generate_file
+    generate_file(filename, "driver", [driver])
+
+    # Reset the caching data structure
+    from dune.perftool.generation import delete_cache_items
+    delete_cache_items()
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index b539e14fdce6942cf288968e47cc2effc7928a65..965efb395b598756ac2be3178da515476eb59bf3 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -1,4 +1,5 @@
 from __future__ import absolute_import
+from os.path import splitext
 
 from dune.perftool.options import get_option
 from dune.perftool.generation import (base_class,
@@ -18,6 +19,31 @@ from dune.perftool.cgen.clazz import (AccessModifier,
 from dune.perftool import Restriction
 
 
+def name_form(formdata, data):
+    # Check wether the formdata has a name in UFL
+    try:
+        name = data.object_names[id(formdata.original_form)]
+        return name
+    except:
+        for index, form in enumerate(data.forms):
+            if formdata.preprocessed_form.equals(form):
+                name = str(index)
+                return name
+    # If the form has no name and can not be found in data.forms something went wrong
+    assert False
+
+
+def name_localoperator_file(formdata, data):
+    from dune.perftool.options import get_option
+    if len(data.forms) == 1:
+        filename = get_option("operator_file")
+    else:
+        suffix = '_' + name_form(formdata, data)
+        basename, extension = splitext(get_option("operator_file"))
+        filename = basename + suffix + extension
+    return filename
+
+
 @template_parameter("operator")
 def lop_template_ansatz_gfs():
     return "GFSU"
@@ -103,8 +129,9 @@ def name_initree_member():
 
 
 @class_basename("operator")
-def localoperator_basename():
-    return "LocalOperator"
+def localoperator_basename(formdata, data):
+    form_name = name_form(formdata, data)
+    return "LocalOperator" + form_name.capitalize()
 
 
 def class_type_from_cache(classtag):
@@ -299,11 +326,11 @@ def generate_localoperator_kernels(formdata, data):
     include_file('dune/pdelab/localoperator/pattern.hh', filetag="operatorfile")
 
     # Trigger this one once early on to avoid wrong stuff happening
-    localoperator_basename()
+    localoperator_basename(formdata, data)
     lop_template_ansatz_gfs()
     lop_template_test_gfs()
     from dune.perftool.pdelab.parameter import parameterclass_basename
-    parameterclass_basename()
+    parameterclass_basename(formdata, data)
 
     # Make sure there is always the same constructor arguments (even if parameter class is empty)
     from dune.perftool.pdelab.localoperator import name_initree_member
@@ -312,6 +339,14 @@ def generate_localoperator_kernels(formdata, data):
     name_paramclass()
 
     base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
+    from dune.perftool.pdelab.driver import is_stationary
+    if not is_stationary():
+        # TODO replace double with clever typename stuff
+        base_class('Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>', classtag="operator")
+
+        # Create set time method in parameter class
+        from dune.perftool.pdelab.parameter import define_set_time_method
+        define_set_time_method()
 
     # Have a data structure collect the generated kernels
     operator_kernels = {}
@@ -406,7 +441,7 @@ def generate_localoperator_kernels(formdata, data):
     return operator_kernels
 
 
-def generate_localoperator_file(kernels):
+def generate_localoperator_file(kernels, filename):
     operator_methods = []
 
     # Make generables from the given kernels
@@ -422,4 +457,14 @@ def generate_localoperator_file(kernels):
     param = cgen_class_from_cache("parameterclass")
     # TODO take the name of this thing from the UFL file
     lop = cgen_class_from_cache("operator", members=operator_methods)
-    generate_file(get_option("operator_file"), "operatorfile", [param, lop])
+    generate_file(filename, "operatorfile", [param, lop])
+
+
+def generate_localoperator_basefile(formdatas, data):
+    filename = get_option("operator_file")
+    for formdata in formdatas:
+        lop_filename = name_localoperator_file(formdata, data)
+        include_file(lop_filename, filetag="operatorbasefile")
+
+    from dune.perftool.file import generate_file
+    generate_file(filename, "operatorbasefile", [])
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index c4e3c5488015445dc954258525370cece2c82d65..05c21e823571df0e5973dc34357df2f17fc476d8 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -24,8 +24,8 @@ from dune.perftool.pdelab.localoperator import (class_type_from_cache,
 
 
 @class_basename("parameterclass")
-def parameterclass_basename():
-    lopbase = localoperator_basename()
+def parameterclass_basename(formdata, data):
+    lopbase = localoperator_basename(formdata, data)
     return "{}Params".format(lopbase)
 
 
@@ -43,6 +43,32 @@ def name_paramclass():
     return "param"
 
 
+@class_member(classtag="parameterclass", access=AccessModifier.PRIVATE)
+def define_time(name):
+    initializer_list(name, ["0.0"], classtag="parameterclass")
+    return "double {};".format(name)
+
+
+@symbol
+def name_time():
+    define_time("t")
+    return "t"
+
+
+@class_member("parameterclass", access=AccessModifier.PUBLIC)
+def define_set_time_method():
+    time_name = name_time()
+    # TODO double?
+    result = ["// Set time in instationary case",
+              "void setTime (double t_)",
+              "{",
+              "  {} = t_;".format(time_name),
+              "}"
+              ]
+
+    return result
+
+
 @class_member("parameterclass", access=AccessModifier.PUBLIC)
 def define_parameter_function_class_member(name, expr, t, cell):
     geot = "E" if cell else "I"
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 76611b318c120c0173096585ccacdeb1b3a6a058..d41b93944ed9bdad71cf18b0d3fbd337589c6012 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,4 +1,5 @@
+add_subdirectory(heatequation)
 add_subdirectory(laplace)
+add_subdirectory(nonlinear)
 add_subdirectory(poisson)
 add_subdirectory(stokes)
-add_subdirectory(nonlinear)
diff --git a/test/heatequation/CMakeLists.txt b/test/heatequation/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..543dd2d08e6b2c862c2f990b1cc992bb559a0c26
--- /dev/null
+++ b/test/heatequation/CMakeLists.txt
@@ -0,0 +1,27 @@
+add_generated_executable(UFLFILE heatequation.ufl
+                         TARGET heatequation_implicit
+                         FORM_COMPILER_ARGS --exact-solution-expression g --compare-l2errorsquared 1e-7
+                         )
+
+
+dune_add_system_test(TARGET heatequation_implicit
+                     INIFILE heatequation_implicit.mini
+                     )
+
+add_generated_executable(UFLFILE heatequation.ufl
+                         TARGET heatequation_explicit
+                         FORM_COMPILER_ARGS --explicit-time-stepping --exact-solution-expression g --compare-l2errorsquared 1e-7
+                         )
+
+
+dune_add_system_test(TARGET heatequation_explicit
+                     INIFILE heatequation_explicit.mini
+                     )
+
+# Hand written reference solutions
+add_executable(heatequation_implicit_ref heatequation_implicit_ref.cc)
+dune_symlink_to_source_files(FILES heatequation_implicit_ref.ini)
+set_target_properties(heatequation_implicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
+
+add_executable(heatequation_explicit_ref heatequation_explicit_ref.cc)
+set_target_properties(heatequation_explicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/heatequation/driver.hh b/test/heatequation/driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1aaf9f171747325885cfebf6a74103e30575d138
--- /dev/null
+++ b/test/heatequation/driver.hh
@@ -0,0 +1,164 @@
+/********************************************************/
+// Beware of line number changes, they may corrupt docu!
+//! \brief Driver function to set up and solve the problem
+/********************************************************/
+
+#include <dune/perftool/vtkpredicate.hh>
+template<typename GV, typename FEM>
+void driver (const GV& gv, const FEM& fem, Dune::ParameterTree& ptree)
+{
+  // dimension and important types
+  const int dim = GV::dimension;
+  typedef double RF;                   // type for computations
+
+  // make user functions and set initial time
+  RF time = 0.0;
+  RF eta = ptree.get("problem.eta",(RF)1.0);
+  Problem<RF> problem(eta);
+  problem.setTime(time);
+  auto glambda = [&](const auto& x){ Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2()); };
+  // auto glambda = [&](const auto& e, const auto& x)
+  //   {return problem.g(e,x);};
+  auto g = Dune::PDELab::
+    makeInstationaryGridFunctionFromCallable(gv,glambda,problem);
+  auto blambda = [&](const auto& i, const auto& x)
+    {return true;};
+  // auto blambda = [&](const auto& i, const auto& x)
+  //   {return problem.b(i,x);};
+  auto b = Dune::PDELab::
+    makeBoundaryConditionFromCallable(gv,blambda);
+
+  // Make grid function space
+  typedef Dune::PDELab::ConformingDirichletConstraints CON;
+  typedef Dune::PDELab::istl::VectorBackend<> VBE;
+  typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE> GFS;
+  GFS gfs(gv,fem);
+  gfs.name("Vh");
+
+  // Assemble constraints
+  typedef typename GFS::template
+    ConstraintsContainer<RF>::Type CC;
+  CC cc;
+  Dune::PDELab::constraints(b,gfs,cc); // assemble constraints
+  std::cout << "constrained dofs=" << cc.size() << " of "
+            << gfs.globalSize() << std::endl;
+
+  // A coefficient vector
+  using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
+  Z z(gfs); // initial value
+
+  // Make a grid function out of it
+  typedef Dune::PDELab::DiscreteGridFunction<GFS,Z> ZDGF;
+  ZDGF zdgf(gfs,z);
+
+  // initialize simulation time,  the coefficient vector
+  Dune::PDELab::interpolate(g,gfs,z);
+
+
+  // Make instationary grid operator
+  typedef NonlinearHeatFEM<Problem<RF>,FEM> LOP;
+  LOP lop(problem);
+  // using LOP = LocalOperatorPoisson<GFS, GFS>;
+  // LocalOperatorPoissonParams params_poisson;
+  // LOP lop(ptree, params_poisson);
+
+
+  typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
+  int degree = ptree.get("fem.degree",(int)1);
+  MBE mbe((int)pow(1+2*degree,dim));
+  typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,
+                                     RF,RF,RF,CC,CC> GO0;
+  GO0 go0(gfs,cc,gfs,cc,lop,mbe);
+
+
+  typedef L2<FEM> TLOP;
+  TLOP tlop;
+  // using TLOP = LocalOperatorMass<GFS, GFS>;
+  // LocalOperatorMassParams params_mass;
+  // TLOP tlop(ptree, params_mass);
+
+
+  typedef Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,
+                                     RF,RF,RF,CC,CC> GO1;
+  GO1 go1(gfs,cc,gfs,cc,tlop,mbe);
+  typedef Dune::PDELab::OneStepGridOperator<GO0,GO1> IGO;
+  IGO igo(go0,go1);
+
+  // Select a linear solver backend
+  typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<IGO> LS;
+  LS ls(100,0);
+
+  // solve nonlinear problem
+  typedef Dune::PDELab::Newton<IGO,LS,Z> PDESOLVER;
+  PDESOLVER newton(igo,z,ls);
+  newton.setReassembleThreshold(0.0);
+  newton.setVerbosityLevel(2);
+  newton.setReduction(1e-8);
+  newton.setMinLinearReduction(1e-4);
+  newton.setMaxIterations(25);
+  newton.setLineSearchMaxIterations(10);
+
+  // select and prepare time-stepping scheme
+  Dune::PDELab::OneStepThetaParameter<RF> method1(1.0);
+  Dune::PDELab::Alexander2Parameter<RF> method2;
+  Dune::PDELab::Alexander3Parameter<RF> method3;
+  int torder = ptree.get("fem.torder",(int)1);
+  Dune::PDELab::TimeSteppingParameterInterface<RF>*
+    pmethod=&method1;
+  if (torder==1) pmethod = &method1;
+  if (torder==2) pmethod = &method2;
+  if (torder==3) pmethod = &method3;
+  if (torder<1||torder>3)
+    std::cout<<"torder not in [1,3]"<<std::endl;
+  typedef Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,Z,Z> OSM;
+  OSM  osm(*pmethod,igo,newton);
+  osm.setVerbosityLevel(2);
+
+  // prepare VTK writer and write first file
+  std::string filename=ptree.get("output.filename","output");
+  struct stat st;
+  if( stat( filename.c_str(), &st ) != 0 )
+    {
+      int stat = 0;
+      stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+      if( stat != 0 && stat != -1)
+        std::cout << "Error: Cannot create directory "
+                  << filename << std::endl;
+    }
+  int subsampling=ptree.get("output.subsampling",(int)0);
+  typedef Dune::SubsamplingVTKWriter<GV> VTKWRITER;
+  VTKWRITER vtkwriter(gv,subsampling);
+  typedef Dune::VTKSequenceWriter<GV> VTKSEQUENCEWRITER;
+  VTKSEQUENCEWRITER vtkSequenceWriter(
+    std::make_shared<VTKWRITER>(vtkwriter),filename,filename,"");
+
+
+  // typedef Dune::PDELab::VTKGridFunctionAdapter<ZDGF> VTKF;
+  // vtkSequenceWriter.addVertexData(std::shared_ptr<VTKF>(
+  //                            new VTKF(zdgf,"solution")));
+  CuttingPredicate predicate;
+  Dune::PDELab::addSolutionToVTKWriter(vtkSequenceWriter, gfs, z, Dune::PDELab::vtk::defaultNameScheme(), predicate);
+  vtkSequenceWriter.write(time,Dune::VTK::appendedraw);
+
+
+  // time loop
+  RF T = ptree.get("problem.T",(RF)1.0);
+  RF dt = ptree.get("fem.dt",(RF)0.1);
+  while (time<T-1e-8)
+    {
+      // assemble constraints for new time step
+      problem.setTime(time+dt);
+      Dune::PDELab::constraints(b,gfs,cc);
+
+      // do time step
+      Z znew(z);
+      osm.apply(time,dt,z,g,znew);
+
+      // accept time step
+      z = znew;
+      time+=dt;
+
+      // output to VTK file
+      vtkSequenceWriter.write(time,Dune::VTK::appendedraw);
+    }
+}
diff --git a/test/heatequation/heatequation.ufl b/test/heatequation/heatequation.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..b6d3f65c7e107483029acd67226bd8e8ec6ab455
--- /dev/null
+++ b/test/heatequation/heatequation.ufl
@@ -0,0 +1,11 @@
+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
+
+V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+mass = (u*v)*dx
+poisson = (inner(grad(u), grad(v)) - f*v)*dx
+
+forms = [mass, poisson]
\ No newline at end of file
diff --git a/test/heatequation/heatequation_explicit.mini b/test/heatequation/heatequation_explicit.mini
new file mode 100644
index 0000000000000000000000000000000000000000..2360fb6bc8436deab3b3cfb68f7146689cbe4d25
--- /dev/null
+++ b/test/heatequation/heatequation_explicit.mini
@@ -0,0 +1,15 @@
+__name = heatequation_explicit
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 8 8
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = heatequation_ref
+extension = vtu
+
+[instat]
+T = 1.0
+dt = 0.01
\ No newline at end of file
diff --git a/test/heatequation/heatequation_explicit_ref.cc b/test/heatequation/heatequation_explicit_ref.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9a07f04958156b7a93d1f17106695c6c81d76312
--- /dev/null
+++ b/test/heatequation/heatequation_explicit_ref.cc
@@ -0,0 +1,263 @@
+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/pdelab/boilerplate/pdelab.hh>
+#include <dune/pdelab/localoperator/convectiondiffusionfem.hh>
+#include <dune/pdelab/localoperator/l2.hh>
+
+//***********************************************************************
+//***********************************************************************
+// diffusion problem with time dependent coefficients
+//***********************************************************************
+//***********************************************************************
+
+const double kx = 2.0, ky = 2.0;
+
+template<typename GV, typename RF>
+class GenericProblem
+{
+  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
+
+public:
+  typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
+
+  GenericProblem () : time(0.0) {}
+
+  //! tensor diffusion coefficient
+  typename Traits::PermTensorType
+  A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::PermTensorType I;
+    for (std::size_t i=0; i<Traits::dimDomain; i++)
+      for (std::size_t j=0; j<Traits::dimDomain; j++)
+        I[i][j] = (i==j) ? 1.0 : 0.0;
+    return I;
+  }
+
+  //! velocity field
+  typename Traits::RangeType
+  b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::RangeType v(0.0);
+    return v;
+  }
+
+  //! sink term
+  typename Traits::RangeFieldType
+  c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! source term
+  typename Traits::RangeFieldType
+  f (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
+  {
+    typename Traits::DomainType x = e.geometry().global(xlocal);
+    Dune::FieldVector<double,2> c(0.5);
+    c-= x;
+    return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+    // return 0.0;
+  }
+
+  //! boundary condition type function
+  BCType
+  bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
+  }
+
+  //! Dirichlet boundary condition value
+  typename Traits::RangeFieldType
+  g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
+  {
+    typename Traits::DomainType x = e.geometry().global(xlocal);
+    Dune::FieldVector<double,2> c(0.5);
+    c-= x;
+    return std::exp(-1.*c.two_norm2());
+    // return std::exp(-(kx*kx+ky*ky)*M_PI*M_PI*time) * sin(kx*M_PI*x[0]) * sin(ky*M_PI*x[1]);
+  }
+
+  //! Neumann boundary condition
+  typename Traits::RangeFieldType
+  j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! Neumann boundary condition
+  typename Traits::RangeFieldType
+  o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! set time for subsequent evaluation
+  void setTime (RF t)
+  {
+    time = t;
+    //std::cout << "setting time to " << time << std::endl;
+  }
+
+private:
+  RF time;
+};
+
+//***********************************************************************
+//***********************************************************************
+// a function that does the simulation on a given grid
+//***********************************************************************
+//***********************************************************************
+
+template<typename GM, unsigned int degree, Dune::GeometryType::BasicType elemtype,
+         Dune::PDELab::MeshType meshtype, Dune::SolverCategory::Category solvertype>
+void do_simulation (double T, double dt, GM& grid, std::string basename)
+{
+  // define parameters
+  typedef double NumberType;
+
+  // make problem parameters
+  typedef GenericProblem<typename GM::LeafGridView,NumberType> Problem;
+  Problem problem;
+  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem> BCType;
+  BCType bctype(grid.leafGridView(),problem);
+
+  // make a finite element space
+  typedef Dune::PDELab::CGSpace<GM,NumberType,degree,BCType,elemtype,meshtype,solvertype> FS;
+  FS fs(grid,bctype);
+
+  // assemblers for finite element problem
+  typedef Dune::PDELab::ConvectionDiffusionFEM<Problem,typename FS::FEM> LOP;
+  LOP lop(problem,4);
+  typedef Dune::PDELab::GalerkinGlobalAssembler<FS,LOP,solvertype> SASS;
+  SASS sass(fs,lop,6);
+  typedef Dune::PDELab::L2 MLOP;
+  MLOP mlop(2*degree);
+  typedef Dune::PDELab::GalerkinGlobalAssembler<FS,MLOP,solvertype> TASS;
+  TASS tass(fs,mlop,6);
+  typedef Dune::PDELab::OneStepGlobalAssembler<SASS,TASS,false> ASSEMBLER;
+  ASSEMBLER assembler(sass,tass);
+
+  // make a degree of freedom vector and set initial value
+  typedef typename FS::DOF V;
+  V x(fs.getGFS(),0.0);
+  typedef Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem> G;
+  G g(grid.leafGridView(),problem);
+  problem.setTime(0.0);
+  Dune::PDELab::interpolate(g,fs.getGFS(),x);
+
+  // linear solver backend
+  typedef Dune::PDELab::ISTLSolverBackend_CG_AMG_SSOR<FS,ASSEMBLER,solvertype> SBE;
+  SBE sbe(fs,assembler,5000,1);
+
+  // linear problem solver
+  typedef Dune::PDELab::StationaryLinearProblemSolver<typename ASSEMBLER::GO,typename SBE::LS,V> PDESOLVER;
+  PDESOLVER pdesolver(*assembler,*sbe,1e-6);
+
+  // // time-stepper
+  // Dune::PDELab::OneStepThetaParameter<NumberType> method(1.0);
+  // typedef Dune::PDELab::OneStepMethod<NumberType,typename ASSEMBLER::GO,PDESOLVER,V> OSM;
+  // OSM osm(method,*assembler,pdesolver);
+  // osm.setVerbosityLevel(2);
+
+  Dune::PDELab::ExplicitEulerParameter<NumberType> method;
+  typedef Dune::PDELab::SimpleTimeController<NumberType> TC;
+  TC tc;
+  typedef Dune::PDELab::ExplicitOneStepMethod<NumberType,typename ASSEMBLER::GO,typename SBE::LS,V,V,TC> OSM;
+  OSM osm(method,*assembler,*sbe,tc);
+  osm.setVerbosityLevel(2);
+
+  // graphics for initial guess
+  Dune::PDELab::FilenameHelper fn(basename);
+  { // start a new block to automatically delete the VTKWriter object
+    Dune::SubsamplingVTKWriter<typename GM::LeafGridView> vtkwriter(grid.leafGridView(),degree-1);
+    typename FS::DGF xdgf(fs.getGFS(),x);
+    vtkwriter.addVertexData(std::make_shared<typename FS::VTKF>(xdgf,"x_h"));
+    vtkwriter.write(fn.getName(),Dune::VTK::appendedraw);
+    fn.increment();
+  }
+
+  // time loop
+  NumberType time = 0.0;
+  while (time<T-1e-10)
+    {
+      // assemble constraints for new time step (assumed to be constant for all substeps)
+      problem.setTime(time+dt);
+      fs.assembleConstraints(bctype);
+
+      // do time step
+      V xnew(fs.getGFS(),0.0);
+      osm.apply(time,dt,x,xnew);
+
+      // output to VTK file
+      {
+        Dune::SubsamplingVTKWriter<typename GM::LeafGridView> vtkwriter(grid.leafGridView(),degree-1);
+        typename FS::DGF xdgf(fs.getGFS(),xnew);
+        vtkwriter.addVertexData(std::make_shared<typename FS::VTKF>(xdgf,"x_h"));
+        vtkwriter.write(fn.getName(),Dune::VTK::appendedraw);
+        fn.increment();
+      }
+
+      // accept time step
+      x = xnew;
+      time += dt;
+    }
+}
+
+//***********************************************************************
+//***********************************************************************
+// the main function
+//***********************************************************************
+//***********************************************************************
+
+int main(int argc, char **argv)
+{
+  // initialize MPI, finalize is done automatically on exit
+  Dune::MPIHelper::instance(argc,argv);
+
+  // read command line arguments
+  if (argc!=4)
+    {
+      std::cout << "usage: " << argv[0] << " <T> <dt> <cells>" << std::endl;
+      return 0;
+    }
+  double T; sscanf(argv[1],"%lg",&T);
+  double dt; sscanf(argv[2],"%lg",&dt);
+  int cells; sscanf(argv[3],"%d",&cells);
+
+  // start try/catch block to get error messages from dune
+  try {
+
+    const int dim=2;
+    const int degree=1;
+    const Dune::SolverCategory::Category solvertype = Dune::SolverCategory::sequential;
+    const Dune::GeometryType::BasicType elemtype = Dune::GeometryType::cube;
+    const Dune::PDELab::MeshType meshtype = Dune::PDELab::MeshType::conforming;
+
+    typedef Dune::YaspGrid<dim> GM;
+    typedef Dune::PDELab::StructuredGrid<GM> Grid;
+    Grid grid(elemtype,cells);
+    grid->loadBalance();
+
+    std::stringstream basename;
+    basename << "heatequation_explicit_ref";
+    do_simulation<GM,degree,elemtype,meshtype,solvertype>(T,dt,*grid,basename.str());
+  }
+  catch (std::exception & e) {
+    std::cout << "STL ERROR: " << e.what() << std::endl;
+    return 1;
+  }
+  catch (Dune::Exception & e) {
+    std::cout << "DUNE ERROR: " << e.what() << std::endl;
+    return 1;
+  }
+  catch (...) {
+    std::cout << "Unknown ERROR" << std::endl;
+    return 1;
+  }
+
+  // done
+  return 0;
+}
diff --git a/test/heatequation/heatequation_implicit.mini b/test/heatequation/heatequation_implicit.mini
new file mode 100644
index 0000000000000000000000000000000000000000..c49499f7833c2af8acbeca3ac86d3f65097476ec
--- /dev/null
+++ b/test/heatequation/heatequation_implicit.mini
@@ -0,0 +1,11 @@
+__name = heatequation_implicit
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = heatequation_ref
+extension = vtu
diff --git a/test/heatequation/heatequation_implicit_ref.cc b/test/heatequation/heatequation_implicit_ref.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a486b85b28941399008026e482e7a22b7bfb76f
--- /dev/null
+++ b/test/heatequation/heatequation_implicit_ref.cc
@@ -0,0 +1,192 @@
+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+// Beware of line number changes, they may corrupt docu!
+/** \file
+    \brief Solve Poisson equation with P1 conforming finite elements
+*/
+// always include the config file
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<math.h>
+#include<iostream>
+// dune-common includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.hh>
+// dune-geometry includes
+#include<dune/geometry/referenceelements.hh>
+#include<dune/geometry/quadraturerules.hh>
+// dune-grid includes
+#include<dune/grid/onedgrid.hh>
+#include<dune/grid/yaspgrid.hh>
+#include<dune/grid/utility/structuredgridfactory.hh>
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/grid/io/file/gmshreader.hh>
+#if HAVE_UG
+#include<dune/grid/uggrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
+#include<dune/alugrid/grid.hh>
+#include<dune/grid/io/file/dgfparser/dgfalu.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+// dune-istl included by pdelab
+// dune-pdelab includes
+#include<dune/pdelab/common/function.hh>
+#include<dune/pdelab/common/vtkexport.hh>
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+#include<dune/pdelab/finiteelementmap/qkfem.hh>
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/common/constraintsparameters.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+#include<dune/pdelab/function/callableadapter.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include<dune/pdelab/gridfunctionspace/interpolate.hh>
+#include<dune/pdelab/gridfunctionspace/vtk.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/gridoperator/onestep.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/backend/istl.hh>
+#include<dune/pdelab/stationary/linearproblem.hh>
+#include<dune/pdelab/instationary/onestep.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+// include all components making up this tutorial
+#include"nonlinearheatfem.hh"
+#include"problem.hh"
+#include"driver.hh"
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+int main(int argc, char** argv)
+{
+  try{
+    // Maybe initialize Mpi
+    Dune::MPIHelper&
+      helper = Dune::MPIHelper::instance(argc, argv);
+    if(Dune::MPIHelper::isFake)
+      std::cout<< "This is a sequential program." << std::endl;
+    else
+      std::cout << "Parallel code run on "
+                << helper.size() << " process(es)" << std::endl;
+
+    // open ini file
+    Dune::ParameterTree ptree;
+    Dune::ParameterTreeParser ptreeparser;
+    ptreeparser.readINITree("heatequation_implicit_ref.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    // read ini file
+    const int dim = ptree.get<int>("grid.dim");
+    const int refinement = ptree.get<int>("grid.refinement");
+    const int degree = ptree.get<int>("fem.degree");
+
+    // in 1d use OneDGrid
+    if (dim==1)
+      {
+        // read grid parameters from input file
+        typedef Dune::OneDGrid::ctype DF;
+        const DF a = ptree.get<DF>("grid.oned.a");
+        const DF b = ptree.get<DF>("grid.oned.b");
+        const unsigned int N = ptree.get<int>("grid.oned.elements");
+
+        // create equidistant intervals
+        typedef std::vector<DF> Intervals;
+        Intervals intervals(N+1);
+        for(unsigned int i=0; i<N+1; ++i)
+          intervals[i] = a + DF(i)*(b-a)/DF(N);
+
+        // Construct grid
+        typedef Dune::OneDGrid Grid;
+        Grid grid(intervals);
+        grid.globalRefine(refinement);
+
+        // call generic function
+        typedef Dune::OneDGrid::LeafGridView GV;
+        GV gv=grid.leafGridView();
+        if (degree==1) {
+          typedef Dune::PDELab::PkLocalFiniteElementMap<GV,DF,double,1> FEM;
+          FEM fem(gv);
+          driver(gv,fem,ptree);
+        }
+        if (degree==2) {
+          typedef Dune::PDELab::PkLocalFiniteElementMap<GV,DF,double,2> FEM;
+          FEM fem(gv);
+          driver(gv,fem,ptree);
+        }
+      }
+
+    // YaspGrid section
+    if (dim==2)
+      {
+        const int dim=2;
+        typedef Dune::YaspGrid<dim> Grid;
+        typedef Grid::ctype DF;
+        Dune::FieldVector<DF,dim> L;
+        L[0] = ptree.get("grid.structured.LX",(double)1.0);
+        L[1] = ptree.get("grid.structured.LY",(double)1.0);
+        std::array<int,dim> N;
+        N[0] = ptree.get("grid.structured.NX",(int)10);
+        N[1] = ptree.get("grid.structured.NY",(int)10);
+        std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>(new Grid(L,N));
+        gridp->globalRefine(refinement);
+        typedef Grid::LeafGridView GV;
+        GV gv=gridp->leafGridView();
+        if (degree==1) {
+          typedef Dune::PDELab::QkLocalFiniteElementMap<GV,DF,double,1> FEM;
+          FEM fem(gv);
+          driver(gv,fem,ptree);
+        }
+        if (degree==2) {
+          typedef Dune::PDELab::QkLocalFiniteElementMap<GV,DF,double,2> FEM;
+          FEM fem(gv);
+          driver(gv,fem,ptree);
+        }
+      }
+    if (dim==3)
+      {
+        const int dim=3;
+        typedef Dune::YaspGrid<dim> Grid;
+        typedef Grid::ctype DF;
+        Dune::FieldVector<DF,dim> L;
+        L[0] = ptree.get("grid.structured.LX",(double)1.0);
+        L[1] = ptree.get("grid.structured.LY",(double)1.0);
+        L[2] = ptree.get("grid.structured.LZ",(double)1.0);
+        std::array<int,dim> N;
+        N[0] = ptree.get("grid.structured.NX",(int)1);
+        N[1] = ptree.get("grid.structured.NY",(int)1);
+        N[2] = ptree.get("grid.structured.NZ",(int)1);
+        std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>(new Grid(L,N));
+        gridp->globalRefine(refinement);
+        typedef Grid::LeafGridView GV;
+        GV gv=gridp->leafGridView();
+        if (degree==1) {
+          typedef Dune::PDELab::QkLocalFiniteElementMap<GV,DF,double,1> FEM;
+          FEM fem(gv);
+          driver(gv,fem,ptree);
+        }
+        if (degree==2) {
+          typedef Dune::PDELab::QkLocalFiniteElementMap<GV,DF,double,2> FEM;
+          FEM fem(gv);
+          driver(gv,fem,ptree);
+        }
+      }
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+}
diff --git a/test/heatequation/heatequation_implicit_ref.ini b/test/heatequation/heatequation_implicit_ref.ini
new file mode 100644
index 0000000000000000000000000000000000000000..7783f5b05e8eb7a431f428329eada972aebb3474
--- /dev/null
+++ b/test/heatequation/heatequation_implicit_ref.ini
@@ -0,0 +1,29 @@
+[grid]
+dim=2
+refinement=3
+
+[grid.structured]
+LX=1.0
+LY=1.0
+LZ=1.0
+NX=4
+NY=4
+NZ=4
+
+[grid.oned]
+a=0.0
+b=3.0
+elements=6
+
+[fem]
+degree=1
+torder=1
+dt=0.1
+
+[problem]
+eta=5.0
+T=1.0
+
+[output]
+filename=tuttut
+subsampling=2
diff --git a/test/heatequation/nonlinearheatfem.hh b/test/heatequation/nonlinearheatfem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bf16ba55bfc1318ffd543667b2f3b9b738eba3b1
--- /dev/null
+++ b/test/heatequation/nonlinearheatfem.hh
@@ -0,0 +1,151 @@
+#include<dune/geometry/referenceelements.hh>
+#include<dune/geometry/quadraturerules.hh>
+#include<dune/pdelab/common/geometrywrapper.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/idefault.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+
+// include the operator from tutorial 01
+#include"nonlinearpoissonfem.hh"
+
+/** a local operator for solving the equation
+ *
+ *  \partial_t u - \Delta u + q(u) = f   in \Omega
+ *                               u = g   on \Gamma_D\subseteq\partial\Omega
+ *              - \nabla u \cdot n = j   on \Gamma_N = \partial\Omega\setminus\Gamma_D
+ *                               u = u_0 at t=t_0
+ *
+ * (spatial part!) with conforming finite elements on all types of grids in any dimension
+ *
+ * \tparam Param    A parameter class providing all coefficient functions in the PDE
+ * \tparam FEM      Type of a finite element map
+ */
+template<typename Param, typename FEM>
+class NonlinearHeatFEM :
+  public NonlinearPoissonFEM<Param,FEM>,
+  public Dune::PDELab::
+    InstationaryLocalOperatorDefaultMethods<double>
+{
+  Param& param;
+public:
+  //! Pass on constructor
+  NonlinearHeatFEM (Param& param_, int incrementorder_=0)
+    : NonlinearPoissonFEM<Param,FEM>(param_,incrementorder_),
+    param(param_)
+  {}
+  //! set time for subsequent evaluation
+  void setTime (double t) {
+    param.setTime(t);
+  }
+};
+
+/** a local operator for the mass operator (L_2 integral)
+ *
+ * \f{align*}{
+ \int_\Omega uv dx
+ * \f}
+ * \tparam FEM      Type of a finite element map
+ */
+template<typename FEM>
+class L2
+  : public Dune::PDELab::FullVolumePattern,
+    public Dune::PDELab::LocalOperatorDefaultFlags,
+    public Dune::PDELab::
+  InstationaryLocalOperatorDefaultMethods<double>
+{
+  typedef typename FEM::Traits::FiniteElementType::
+     Traits::LocalBasisType LocalBasis;
+  Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+
+public:
+  // pattern assembly flags
+  enum { doPatternVolume = true };
+
+  // residual assembly flags
+  enum { doAlphaVolume = true };
+
+  // volume integral depending on test and ansatz functions
+  template<typename EG, typename LFSU, typename X,
+           typename LFSV, typename R>
+  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
+                     const LFSV& lfsv, R& r) const
+  {
+    // types & dimension
+    typedef decltype(makeZeroBasisFieldValue(lfsu)) RF;
+
+    // select quadrature rule
+    auto geo = eg.geometry();
+    const int order = 2*lfsu.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(geo,order);
+
+    // loop over quadrature points
+    for (const auto& ip : rule)
+      {
+        // evaluate basis functions
+        auto& phihat =
+          cache.evaluateFunction(ip.position(),
+                                 lfsu.finiteElement().localBasis());
+
+        // evaluate u
+        RF u=0.0;
+        for (size_t i=0; i<lfsu.size(); i++) u += x(lfsu,i)*phihat[i];
+
+        // integrate u*phi_i
+        RF factor = ip.weight() * geo.integrationElement(ip.position());
+        for (size_t i=0; i<lfsv.size(); i++)
+          r.accumulate(lfsv,i,u*phihat[i]*factor);
+      }
+  }
+
+  //! jacobian contribution of volume term
+  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
+  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
+                        M& mat) const
+  {
+    // types & dimension
+    typedef decltype(makeZeroBasisFieldValue(lfsu)) RF;
+
+    // select quadrature rule
+    auto geo = eg.geometry();
+    const int order = 2*lfsu.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(geo,order);
+
+    // loop over quadrature points
+    for (const auto& ip : rule)
+      {
+        // evaluate basis functions
+        auto& phihat =
+          cache.evaluateFunction(ip.position(),
+                                 lfsu.finiteElement().localBasis());
+
+        // integrate phi_j*phi_i
+        RF factor = ip.weight() * geo.integrationElement(ip.position());
+        for (size_t j=0; j<lfsu.size(); j++)
+          for (size_t i=0; i<lfsv.size(); i++)
+            mat.accumulate(lfsv,i,lfsu,j,phihat[j]*phihat[i]*factor);
+      }
+  }
+
+  //! apply local jacobian of the volume term -> nonlinear variant
+  template<typename EG, typename LFSU, typename X,
+           typename LFSV, typename R>
+  void jacobian_apply_volume (const EG& eg, const LFSU& lfsu,
+                              const X& x, const X& z, const LFSV& lfsv,
+                              R& r) const
+  {
+    alpha_volume(eg,lfsu,z,lfsv,r);
+  }
+
+  //! apply local jacobian of the volume term -> linear variant
+  template<typename EG, typename LFSU, typename X,
+           typename LFSV, typename R>
+  void jacobian_apply_volume (const EG& eg, const LFSU& lfsu,
+                              const X& x, const LFSV& lfsv,
+                              R& r) const
+  {
+    alpha_volume(eg,lfsu,x,lfsv,r);
+  }
+
+};
diff --git a/test/heatequation/nonlinearpoissonfem.hh b/test/heatequation/nonlinearpoissonfem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..032f737b9028532fd388adc203e35e15f803c715
--- /dev/null
+++ b/test/heatequation/nonlinearpoissonfem.hh
@@ -0,0 +1,171 @@
+#include<vector>
+
+#include<dune/common/exceptions.hh>
+#include<dune/common/fvector.hh>
+#include<dune/geometry/type.hh>
+
+#include<dune/pdelab/common/referenceelements.hh>
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+
+/** a local operator for solving the nonlinear Poisson equation with conforming FEM
+ *
+ * \f{align*}{
+ *   -\Delta u(x) + q(u(x)) &=& f(x) x\in\Omega,  \\
+ *                     u(x) &=& g(x) x\in\partial\Omega_D \\
+ *  -\nabla u(x) \cdot n(x) &=& j(x) x\in\partial\Omega_N \\
+ * \f}
+ *
+ */
+template<typename Param, typename FEM>
+class NonlinearPoissonFEM :
+  public Dune::PDELab::
+    NumericalJacobianVolume<NonlinearPoissonFEM<Param,FEM> >,
+  public Dune::PDELab::
+    NumericalJacobianApplyVolume<NonlinearPoissonFEM<Param,FEM> >,
+  public Dune::PDELab::FullVolumePattern,
+  public Dune::PDELab::LocalOperatorDefaultFlags
+{
+  typedef typename FEM::Traits::FiniteElementType::
+     Traits::LocalBasisType LocalBasis;
+  Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+  Param& param; // parameter functions
+  int incrementorder; // increase of integration order
+
+public:
+  // pattern assembly flags
+  enum { doPatternVolume = true };
+
+  // residual assembly flags
+  enum { doLambdaVolume = true };
+  enum { doLambdaBoundary = true };
+  enum { doAlphaVolume = true };
+
+  //! constructor stores a copy of the parameter object
+  NonlinearPoissonFEM (Param& param_, int incrementorder_=0)
+    : param(param_), incrementorder(incrementorder_)
+  {}
+
+  //! right hand side integral
+  template<typename EG, typename LFSV, typename R>
+  void lambda_volume (const EG& eg, const LFSV& lfsv,
+                      R& r) const
+  {
+    // select quadrature rule
+    auto geo = eg.geometry();
+    const int order = incrementorder+
+      2*lfsv.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(geo,order);
+
+    // loop over quadrature points
+    for (const auto& ip : rule)
+      {
+        // evaluate basis functions
+        auto& phihat = cache.evaluateFunction(ip.position(),
+                             lfsv.finiteElement().localBasis());
+
+        // integrate -f*phi_i
+        decltype(ip.weight()) factor = ip.weight()*
+          geo.integrationElement(ip.position());
+        auto f=param.f(eg.entity(),ip.position());
+        for (size_t i=0; i<lfsv.size(); i++)
+          r.accumulate(lfsv,i,-f*phihat[i]*factor);
+      }
+  }
+
+  // Neumann boundary integral
+  template<typename IG, typename LFSV, typename R>
+  void lambda_boundary (const IG& ig, const LFSV& lfsv,
+                        R& r) const
+  {
+    // evaluate boundary condition type
+    auto localgeo = ig.geometryInInside();
+    auto facecenterlocal = Dune::PDELab::
+      referenceElement(localgeo).position(0,0);
+    bool isdirichlet=param.b(ig.intersection(),facecenterlocal);
+
+    // skip rest if we are on Dirichlet boundary
+    if (isdirichlet) return;
+
+    // select quadrature rule
+    auto globalgeo = ig.geometry();
+    const int order = incrementorder+
+      2*lfsv.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(globalgeo,order);
+
+    // loop over quadrature points and integrate normal flux
+    for (const auto& ip : rule)
+      {
+        // quadrature point in local coordinates of element
+        auto local = localgeo.global(ip.position());
+
+        // evaluate shape functions (assume Galerkin method)
+        auto& phihat = cache.evaluateFunction(local,
+                             lfsv.finiteElement().localBasis());
+
+        // integrate j
+        decltype(ip.weight()) factor = ip.weight()*
+          globalgeo.integrationElement(ip.position());
+        auto j = param.j(ig.intersection(),ip.position());
+        for (size_t i=0; i<lfsv.size(); i++)
+          r.accumulate(lfsv,i,j*phihat[i]*factor);
+      }
+  }
+
+  //! volume integral depending on test and ansatz functions
+  template<typename EG, typename LFSU, typename X,
+           typename LFSV, typename R>
+  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
+                     const LFSV& lfsv, R& r) const
+  {
+    // types & dimension
+    const int dim = EG::Entity::dimension;
+    typedef decltype(Dune::PDELab::
+                     makeZeroBasisFieldValue(lfsu)) RF;
+
+    // select quadrature rule
+    auto geo = eg.geometry();
+    const int order = incrementorder+
+      2*lfsu.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(geo,order);
+
+    // loop over quadrature points
+    for (const auto& ip : rule)
+      {
+        // evaluate basis functions
+        auto& phihat = cache.evaluateFunction(ip.position(),
+                             lfsu.finiteElement().localBasis());
+
+        // evaluate u
+        RF u=0.0;
+        for (size_t i=0; i<lfsu.size(); i++)
+          u += x(lfsu,i)*phihat[i];
+
+        // evaluate gradient of shape functions
+        auto& gradphihat = cache.evaluateJacobian(ip.position(),
+                             lfsu.finiteElement().localBasis());
+
+        // transform gradients of shape functions to real element
+        const auto S = geo.jacobianInverseTransposed(ip.position());
+        auto gradphi = makeJacobianContainer(lfsu);
+        for (size_t i=0; i<lfsu.size(); i++)
+          S.mv(gradphihat[i][0],gradphi[i][0]);
+
+        // compute gradient of u
+        Dune::FieldVector<RF,dim> gradu(0.0);
+        for (size_t i=0; i<lfsu.size(); i++)
+          gradu.axpy(x(lfsu,i),gradphi[i][0]);
+
+        // integrate (grad u)*grad phi_i + q(u)*phi_i
+        auto factor = ip.weight()*
+          geo.integrationElement(ip.position());
+        auto q = param.q(u);
+        for (size_t i=0; i<lfsu.size(); i++)
+          r.accumulate(lfsu,i,(gradu*gradphi[i][0]+
+                               q*phihat[i])*factor);
+      }
+  }
+};
diff --git a/test/heatequation/problem.hh b/test/heatequation/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8c1a13c6d53b2c822b1a4c276f23c4f47a80f0d3
--- /dev/null
+++ b/test/heatequation/problem.hh
@@ -0,0 +1,70 @@
+template<typename Number>
+class Problem
+{
+  Number eta;
+  Number t;
+
+public:
+  typedef Number value_type;
+
+  //! Constructor without arg sets nonlinear term to zero
+  Problem () : eta(0.0), t(0.0) {}
+
+  //! Constructor takes eta parameter
+  Problem (const Number& eta_) : eta(eta_), t(0.0) {}
+
+  //! nonlinearity
+  Number q (Number u) const
+  {
+    return 0.0;
+  }
+
+  //! derivative of nonlinearity
+  Number qprime (Number u) const
+  {
+    return 0.0;
+  }
+
+  //! right hand side
+  template<typename E, typename X>
+  Number f (const E& e, const X& local) const
+  {
+    // return 0.0;
+    auto x = e.geometry().global(local);
+    Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+  }
+
+  //! boundary condition type function (true = Dirichlet)
+  template<typename I, typename X>
+  bool b (const I& i, const X& x) const
+  {
+    auto global = i.geometry().global(x);
+    return (global[0]<=1e-7) ? true : false;
+  }
+
+  //! Dirichlet extension
+  template<typename E, typename X>
+  Number g (const E& e, const X& x) const
+  {
+    auto global = e.geometry().global(x);
+    Number s=sin(2.0*M_PI*t);
+    for (std::size_t i=1; i<global.size(); i++)
+      s*=sin(global[i]*M_PI)*sin(global[i]*M_PI);
+    for (std::size_t i=1; i<global.size(); i++)
+      s*=sin(10*global[i]*M_PI)*sin(10*global[i]*M_PI);
+    return s;
+  }
+
+  //! Neumann boundary condition
+  template<typename I, typename X>
+  Number j (const I& i, const X& x) const
+  {
+    return 0.0;
+  }
+
+  //! Set time in instationary case
+  void setTime (Number t_)
+  {
+    t = t_;
+  }
+};