diff --git a/python/dune/codegen/pdelab/driver/error.py b/python/dune/codegen/pdelab/driver/error.py
index 21cc1fba0f20edd30ce2894b2137a386c37b87e8..cf288249a4d4f784d96e32c22f3353dfbd4bb93f 100644
--- a/python/dune/codegen/pdelab/driver/error.py
+++ b/python/dune/codegen/pdelab/driver/error.py
@@ -12,10 +12,11 @@ from dune.codegen.pdelab.driver import (get_form_ident,
 from dune.codegen.pdelab.driver.gridfunctionspace import (main_type_trial_gfs,
                                                           name_leafview,
                                                           name_trial_subgfs,
-                                                          type_range,
+                                                          main_type_range,
                                                           )
 from dune.codegen.pdelab.driver.interpolate import (interpolate_vector,
-                                                    name_boundary_function,
+                                                    main_name_boundary_grid_function,
+                                                    main_type_boundary_grid_function,
                                                     )
 from dune.codegen.pdelab.driver.solve import (define_vector,
                                               dune_solve,
@@ -43,7 +44,7 @@ def name_exact_solution_gridfunction(treepath):
         index = treepath_to_index(element, treepath)
         func = (func[index],)
         element = element.extract_component(index)[1]
-    return name_boundary_function(element, func)
+    return main_name_boundary_grid_function(element, func)
 
 
 def type_discrete_grid_function(gfs):
@@ -68,12 +69,18 @@ def name_discrete_grid_function(gfs, vector_name):
 
 @preamble(section="error", kernel="main")
 def typedef_difference_squared_adapter(name, treepath):
-    sol = name_exact_solution_gridfunction(treepath)
+    element = get_trial_element()
+    func = preprocess_leaf_data(element, "exact_solution")
+    if isinstance(element, MixedElement):
+        index = treepath_to_index(element, treepath)
+        func = (func[index],)
+        element = element.extract_component(index)[1]
+    bgf_type = main_type_boundary_grid_function(func)
     vector = main_name_vector(get_form_ident())
     gfs = name_trial_subgfs(treepath)
     dgf = name_discrete_grid_function(gfs, vector)
     dgf_type = type_discrete_grid_function(gfs)
-    return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), {}>;'.format(name, sol, dgf_type)
+    return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<{}, {}>;'.format(name, bgf_type, dgf_type)
 
 
 def type_difference_squared_adapter(treepath):
@@ -90,7 +97,7 @@ def define_difference_squared_adapter(name, treepath):
     gfs = name_trial_subgfs(treepath)
     dgf = name_discrete_grid_function(gfs, vector)
 
-    return '{} {}({}, {});'.format(t, name, sol, dgf)
+    return '{} {}(*{}, {});'.format(t, name, sol, dgf)
 
 
 def name_difference_squared_adapter(treepath):
@@ -164,7 +171,7 @@ def accumulate_L2_squared():
 
 @preamble(section="error", kernel="main")
 def define_accumulated_L2_error(name):
-    t = type_range()
+    t = main_type_range()
     return "Dune::FieldVector<{}, 1> {}(0.0);".format(t, name)
 
 
diff --git a/python/dune/codegen/pdelab/driver/gridfunctionspace.py b/python/dune/codegen/pdelab/driver/gridfunctionspace.py
index cef1832b80c8adc96169807bcc9dab0498e9f47f..e86077a6415da80acb2d7ee2faf132bdfb5b6eb7 100644
--- a/python/dune/codegen/pdelab/driver/gridfunctionspace.py
+++ b/python/dune/codegen/pdelab/driver/gridfunctionspace.py
@@ -46,6 +46,17 @@ def type_range():
     return name
 
 
+@preamble(section="init", kernel="main")
+def main_typedef_range(name):
+    return "using {} = {};".format(name, type_floatingpoint())
+
+
+def main_type_range():
+    name = "RangeType"
+    main_typedef_range(name)
+    return name
+
+
 @preamble(section="grid", kernel="main")
 def typedef_grid(name):
     dim = get_dimension()
diff --git a/python/dune/codegen/pdelab/driver/gridoperator.py b/python/dune/codegen/pdelab/driver/gridoperator.py
index 7ddbb50f0f7f4d40f2c02a24e35721c84e8641d5..2780344d3ddf6f6027bee455c09f0f524ac4d3cf 100644
--- a/python/dune/codegen/pdelab/driver/gridoperator.py
+++ b/python/dune/codegen/pdelab/driver/gridoperator.py
@@ -81,6 +81,20 @@ def name_gridoperator(form_ident):
     return name
 
 
+@preamble(section="postprocessing", kernel="main")
+def main_typedef_gridoperator(name, form_ident):
+    driver_block_type = type_driver_block()
+    db_gridoperator_type = type_gridoperator(form_ident)
+    gridoperator_type = "using {} = {}::{};".format(name, driver_block_type, db_gridoperator_type)
+    return gridoperator_type
+
+
+def main_type_gridoperator(form_ident):
+    name = "GO_{}".format(form_ident)
+    main_typedef_gridoperator(name, form_ident)
+    return name
+
+
 @class_member(classtag="driver_block")
 def driver_block_get_gridoperator(form_ident):
     gridoperator_type = type_gridoperator(form_ident)
@@ -97,20 +111,6 @@ def main_define_gridoperator(name, form_ident):
     return "auto {} = {}.get_gridoperator();".format(name, driver_block_name)
 
 
-@preamble(section="postprocessing", kernel="main")
-def main_typedef_gridoperator(name, form_ident):
-    driver_block_type = type_driver_block()
-    db_gridoperator_type = type_gridoperator(form_ident)
-    gridoperator_type = "using {} = {}::{};".format(name, driver_block_type, db_gridoperator_type)
-    return gridoperator_type
-
-
-def main_type_gridoperator(form_ident):
-    name = "GO_{}".format(form_ident)
-    main_typedef_gridoperator(name, form_ident)
-    return name
-
-
 def main_name_gridoperator(form_ident):
     name = "go_{}".format(form_ident)
     main_define_gridoperator(name, form_ident)
diff --git a/python/dune/codegen/pdelab/driver/instationary.py b/python/dune/codegen/pdelab/driver/instationary.py
index eb01fe89f519d25cbfaa9c9196a705bf4b9b87c0..6b7d19f432c5525248233bbef31dc9f8e98d77c5 100644
--- a/python/dune/codegen/pdelab/driver/instationary.py
+++ b/python/dune/codegen/pdelab/driver/instationary.py
@@ -19,7 +19,7 @@ from dune.codegen.pdelab.driver.constraints import (has_dirichlet_constraints,
                                                     name_constraintscontainer,
                                                     )
 from dune.codegen.pdelab.driver.interpolate import (interpolate_dirichlet_data,
-                                                    name_boundary_function,
+                                                    name_boundary_grid_function,
                                                     )
 from dune.codegen.pdelab.driver.solve import (print_matrix,
                                               print_residual,
@@ -80,7 +80,7 @@ def time_loop():
         osm = name_onestepmethod()
     if has_dirichlet_constraints(is_dirichlet):
         dirichlet = preprocess_leaf_data(element, "interpolate_expression")
-        boundary = name_boundary_function(element, dirichlet)
+        boundary = name_boundary_grid_function(element, dirichlet)
         apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
     else:
         apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
diff --git a/python/dune/codegen/pdelab/driver/interpolate.py b/python/dune/codegen/pdelab/driver/interpolate.py
index fdc8cd2afea3369cbbdeb809544902c719e0922c..9814d59e58f6f6f7b31164238485e2e4a3f052cd 100644
--- a/python/dune/codegen/pdelab/driver/interpolate.py
+++ b/python/dune/codegen/pdelab/driver/interpolate.py
@@ -1,4 +1,5 @@
 from dune.codegen.generation import (cached,
+                                     class_member,
                                      get_counted_variable,
                                      global_context,
                                      include_file,
@@ -10,8 +11,12 @@ from dune.codegen.pdelab.driver import (FEM_name_mangling,
                                         is_stationary,
                                         preprocess_leaf_data,
                                         )
+from dune.codegen.pdelab.driver.driverblock import (name_driver_block,
+                                                    type_driver_block)
 from dune.codegen.pdelab.driver.gridfunctionspace import (name_trial_gfs,
                                                           name_leafview,
+                                                          type_leafview,
+                                                          type_range,
                                                           )
 from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
 
@@ -20,37 +25,17 @@ def interpolate_dirichlet_data(name):
     element = get_trial_element()
     func = preprocess_leaf_data(element, "interpolate_expression", applyZeroDefault=False)
     if func is not None:
-        bf = name_boundary_function(element, func)
+        bf = name_boundary_grid_function(element, func)
         gfs = name_trial_gfs()
         interpolate_vector(bf, gfs, name)
 
 
 @preamble(section="vector", kernel="driver_block")
 def interpolate_vector(func, gfs, name):
-    return "Dune::PDELab::interpolate({}, *{}, *{});".format(func,
-                                                             gfs,
-                                                             name,
-                                                             )
-
-
-@cached
-def name_boundary_function(element, func):
-    if isinstance(element, MixedElement):
-        k = 0
-        childs = []
-        for subel in element.sub_elements():
-            childs.append(name_boundary_function(subel, func[k:k + subel.value_size()]))
-            k = k + subel.value_size()
-        name = "_".join(childs)
-        if len(childs) == 1:
-            name = "{}_dummy".format(name)
-        define_compositegfs_parameterfunction(name, tuple(childs))
-        return name
-    else:
-        assert isinstance(element, (FiniteElement, TensorProductElement))
-        name = get_counted_variable("func")
-        define_boundary_function(name, func[0])
-        return name
+    return "Dune::PDELab::interpolate(*{}, *{}, *{});".format(func,
+                                                              gfs,
+                                                              name,
+                                                              )
 
 
 @preamble(section="vector", kernel="driver_block")
@@ -60,17 +45,43 @@ def define_compositegfs_parameterfunction(name, children):
                                                                     ', '.join(children))
 
 
+@class_member(classtag="driver_block")
+def typedef_boundary_grid_function(name, func):
+    leafview_type = type_leafview()
+    range_type = type_range()
+    boundary_function_type = type_boundary_function(func)
+    # palpo TODO: 1 in the format below!
+    return "using {} = Dune::PDELab::LocalCallableToGridFunctionAdapter<{}, {}, {}, {}>;".format(name,
+                                                                                                 leafview_type,
+                                                                                                 range_type,
+                                                                                                 1,
+                                                                                                 boundary_function_type)
+
+
+def type_boundary_grid_function(func):
+    name = "BoundaryGridFunction"
+    typedef_boundary_grid_function(name, func)
+    return name
+
+
+@class_member(classtag="driver_block")
+def declare_boundary_grid_function(name, func):
+    bgf_type = type_boundary_grid_function(func)
+    return "std::shared_ptr<{}> {};".format(bgf_type, name)
+
+
 @preamble(section="vector", kernel="driver_block")
-def define_boundary_function(name, dirichlet):
+def define_boundary_grid_function(name, func):
+    declare_boundary_grid_function(name, func)
     gv = name_leafview()
-    lambdaname = name_boundary_lambda(dirichlet)
+    boundary_function = name_boundary_function(func)
+    bgf_type = type_boundary_grid_function(func)
     include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
     if is_stationary():
-        return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
-                                                                                      gv,
-                                                                                      lambdaname,
-                                                                                      )
+        return "{} = std::make_shared<{}>({}, *{});".format(name, bgf_type, gv, boundary_function)
     else:
+        # palpo TODO
+        assert False
         from dune.codegen.pdelab.driver.gridoperator import name_localoperator
         lop = name_localoperator(get_form_ident())
         return "auto {} = Dune::PDELab::makeInstationaryGridFunctionFromCallable({}, {}, {});".format(name,
@@ -81,21 +92,128 @@ def define_boundary_function(name, dirichlet):
 
 
 @cached
-def name_boundary_lambda(boundary):
-    name = get_counted_variable("lambda")
-    define_boundary_lambda(name, boundary)
-    return name
+def name_boundary_grid_function(element, func):
+    assert isinstance(func, tuple)
+
+    if isinstance(element, MixedElement):
+        # palpo TODO
+        assert False
+        k = 0
+        childs = []
+        for subel in element.sub_elements():
+            childs.append(name_boundary_grid_function(subel, func[k:k + subel.value_size()]))
+            k = k + subel.value_size()
+        name = "_".join(childs)
+        if len(childs) == 1:
+            name = "{}_dummy".format(name)
+        define_compositegfs_parameterfunction(name, tuple(childs))
+        return name
+    else:
+        assert isinstance(element, (FiniteElement, TensorProductElement))
+        name = "boundary_grid_function"
+        define_boundary_grid_function(name, func)
+        return name
 
 
-@preamble(section="vector", kernel="driver_block")
-def define_boundary_lambda(name, boundary):
-    if boundary is None:
-        boundary = 0.0
+def boundary_lambda(func):
+    assert isinstance(func, tuple)
+    assert len(func) == 1
+    func = func[0]
+
+    if func is None:
+        func = 0.0
     from ufl.classes import Expr
-    if isinstance(boundary, (int, float)):
-        return "auto {} = [&](const auto& x){{ return {}; }};".format(name, float(boundary))
-    elif isinstance(boundary, Expr):
+    if isinstance(func, (int, float)):
+        return "[&](const auto& x){{ return {}; }}".format(float(func))
+    elif isinstance(func, Expr):
         from dune.codegen.pdelab.driver.visitor import ufl_to_code
-        return "auto {} = [&](const auto& is, const auto& xl){{ {}; }};".format(name, ufl_to_code(boundary))
+        return "[&](const auto& is, const auto& xl){{ {}; }}".format(ufl_to_code(func))
+    else:
+        raise NotImplementedError("What is this?")
+
+
+@class_member(classtag="driver_block")
+def typedef_boundary_function(name, func):
+    assert isinstance(func, tuple)
+    assert len(func) == 1
+    func = func[0]
+
+    range_type = type_range()
+    leafview_type = type_leafview()
+    if func is None:
+        func = 0.0
+    from ufl.classes import Expr
+    if isinstance(func, (int, float)):
+        coordinate = "typename {}::template Codim<0>::Entity::Geometry::GlobalCoordinate".format(leafview_type)
+        return "using {} = std::function<{}({})>;".format(name, range_type, coordinate)
+    elif isinstance(func, Expr):
+        entity = "typename {}::template Codim<0>::Entity".format(leafview_type)
+        coordinate = "typename {}::template Codim<0>::Entity::Geometry::LocalCoordinate".format(leafview_type)
+        return "using {} = std::function<{}({}, {})>;".format(name, range_type, entity, coordinate)
     else:
         raise NotImplementedError("What is this?")
+
+
+def type_boundary_function(func):
+    name = "BoundaryFunction"
+    typedef_boundary_function(name, func)
+    return name
+
+
+@class_member(classtag="driver_block")
+def declare_boundary_function(name, func):
+    bf_type = type_boundary_function(func)
+    return "std::shared_ptr<{}> {};".format(bf_type, name)
+
+
+@preamble(section="vector", kernel="driver_block")
+def define_boundary_function(name, func):
+    declare_boundary_function(name, func)
+    bf_type = type_boundary_function(func)
+    bf_lambda = boundary_lambda(func)
+    return "{} = std::make_shared<{}> ({});".format(name, bf_type, bf_lambda)
+
+
+@cached
+def name_boundary_function(func):
+    name = get_counted_variable("function")
+    define_boundary_function(name, func)
+    return name
+
+
+@class_member(classtag="driver_block")
+def driver_block_get_boundarygridfunction(element, func):
+    name = name_boundary_grid_function(element, func)
+    bgf_type = type_boundary_grid_function(func)
+    return ["std::shared_ptr<{}> get_boundarygridfunction(){{".format(bgf_type),
+            "  return {};".format(name),
+            "}"]
+
+
+@preamble(section="postprocessing", kernel="main")
+def main_typedef_boundary_grid_function(name, func):
+    driver_block_type = type_driver_block()
+    bgf_type = type_boundary_grid_function(func)
+    return "using {} = {}::{};".format(name, driver_block_type, bgf_type)
+
+
+def main_type_boundary_grid_function(func):
+    name = "BoundaryGridFunction"
+    main_typedef_boundary_grid_function(name, func)
+    return name
+
+
+@preamble(section="postprocessing", kernel="main")
+def main_define_boundar_grid_function(name, element, func):
+    driver_block_name = name_driver_block()
+    driver_block_get_boundarygridfunction(element, func)
+    return "auto {} = {}.get_boundarygridfunction();".format(name, driver_block_name)
+
+
+@cached
+def main_name_boundary_grid_function(element, func):
+    assert isinstance(func, tuple)
+
+    name = "boundary_grid_function"
+    main_define_boundar_grid_function(name, element, func)
+    return name