diff --git a/python/dune/perftool/__init__.py b/python/dune/perftool/__init__.py
index 02ef874ff476ba7fc9bca5ecf261f0fcf7073edb..b907a8008731bcb55225ce08fc250d821da4a26e 100644
--- a/python/dune/perftool/__init__.py
+++ b/python/dune/perftool/__init__.py
@@ -1,4 +1,6 @@
-class Restriction:
-    NONE = 0
-    NEGATIVE = 1
-    POSITIVE = 2
+from dune.perftool.options import get_option
+
+# Trigger some imports that are needed to have all backend implementations visible
+# to the selection mechanisms
+import dune.perftool.pdelab
+import dune.perftool.sumfact
diff --git a/python/dune/perftool/generation/__init__.py b/python/dune/perftool/generation/__init__.py
index 32463b4b1ccd0af0ececa244c5578b5ccc4a6d67..ea94e8258aed4595a900ebe64bebf27bdd721780 100644
--- a/python/dune/perftool/generation/__init__.py
+++ b/python/dune/perftool/generation/__init__.py
@@ -37,7 +37,6 @@ from dune.perftool.generation.loopy import (constantarg,
                                             iname,
                                             instruction,
                                             noop_instruction,
-                                            pymbolic_expr,
                                             silenced_warning,
                                             temporary_variable,
                                             valuearg,
diff --git a/python/dune/perftool/generation/backend.py b/python/dune/perftool/generation/backend.py
index 27356c937f1c8617297736a34e74d540dd0996d2..0e7b7ec2c705bf7f71fa642544f10861257e2352 100644
--- a/python/dune/perftool/generation/backend.py
+++ b/python/dune/perftool/generation/backend.py
@@ -1,4 +1,5 @@
 from dune.perftool.generation.cache import _RegisteredFunction
+from dune.perftool.options import option_switch
 from pytools import Record
 
 
@@ -9,10 +10,12 @@ class FuncProxy(Record):
     def __init__(self, interface, name, func):
         Record.__init__(self, interface=interface, name=name, func=func)
 
+    def __call__(self, *args, **kwargs):
+        return self.func(*args, **kwargs)
+
 
 def register_backend(interface, name, func):
     _backend_mapping.setdefault(interface, {})
-    assert name not in _backend_mapping[interface]
     _backend_mapping[interface][name] = func
 
 
@@ -24,15 +27,15 @@ def backend(interface=None, name='default'):
             # Allow order independence of the generator decorators
             # and the backend decorator by delaying the registration
             func = FuncProxy(interface, name, func)
-        else:
-            register_backend(interface, name, func)
+
+        register_backend(interface, name, func)
 
         return func
 
     return _dec
 
 
-def get_backend(interface=None, selector=None, **kwargs):
+def get_backend(interface=None, selector=option_switch("sumfact"), **kwargs):
     assert interface and selector
 
     select = selector(**kwargs)
diff --git a/python/dune/perftool/generation/loopy.py b/python/dune/perftool/generation/loopy.py
index f62f309c069b3cbec2ca3a9ebce6f13cd7653646..7acc9a0eb4a4ea3ab7203e04608e70a6fd5bb99a 100644
--- a/python/dune/perftool/generation/loopy.py
+++ b/python/dune/perftool/generation/loopy.py
@@ -11,7 +11,6 @@ import loopy
 import numpy
 
 iname = generator_factory(item_tags=("iname",))
-pymbolic_expr = generator_factory(item_tags=("kernel", "pymbolic"))
 function_mangler = generator_factory(item_tags=("mangler",))
 silenced_warning = generator_factory(item_tags=("silenced_warning",), no_deco=True)
 
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 60f66248d2dee8d521f04cfa55acd3337f70454a..a4c2e21070b45cbe6d3d6273744fb940181399e1 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -111,3 +111,12 @@ def get_option(key, default=None):
         init_option_dict()
 
         return _option_dict.get(key, default)
+
+
+def option_switch(opt):
+    def _switch():
+        if get_option(opt):
+            return opt
+        else:
+            return "default"
+    return _switch
\ No newline at end of file
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index abbc10ffa31ea40dc38fa46bdcfe32a7fdbf26b7..76aaff919dd3d63045f811ebe9896a5910b529d5 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -1,30 +1,106 @@
 """ The pdelab specific parts of the code generation process """
 
-from dune.perftool.generation import (preamble,
-                                      cached,
-                                      )
-
-
-# Now define some commonly used generators that do not fall into a specific category
-@cached
-def name_index(index):
-    from ufl.classes import MultiIndex, Index
-    if isinstance(index, Index):
-        # This failed for index > 9 because ufl placed curly brackets around
-        # return str(index)
-        return "i_{}".format(index.count())
-    if isinstance(index, MultiIndex):
-        assert len(index) == 1
-        # return str(index._indices[0])
-        return "i_{}".format(index._indices[0].count())
-    raise NotImplementedError
-
-
-def restricted_name(name, restriction):
-    from dune.perftool import Restriction
-    if restriction == Restriction.NONE:
-        return name
-    if restriction == Restriction.POSITIVE:
-        return name + '_n'
-    if restriction == Restriction.NEGATIVE:
-        return name + '_s'
+# Trigger some imports that are needed to have all backend implementations visible
+# to the selection mechanisms
+from dune.perftool.pdelab.argument import (pymbolic_apply_function,
+                                           pymbolic_apply_function_gradient,
+                                           pymbolic_trialfunction,
+                                           pymbolic_trialfunction_gradient,
+                                           )
+from dune.perftool.pdelab.basis import (name_basis,
+                                        name_reference_gradient,
+                                        )
+from dune.perftool.pdelab.geometry import (dimension_iname,
+                                           name_facet_jacobian_determinant,
+                                           name_jacobian_determinant,
+                                           name_jacobian_inverse_transposed,
+                                           name_unit_inner_normal,
+                                           name_unit_outer_normal,
+                                           )
+from dune.perftool.pdelab.index import (name_index,
+                                        )
+from dune.perftool.pdelab.parameter import (cell_parameter_function,
+                                            intersection_parameter_function,
+                                            )
+from dune.perftool.pdelab.quadrature import (name_quadrature_weight,
+                                             )
+from dune.perftool.pdelab.spaces import (lfs_iname,
+                                         )
+
+
+class PDELabInterface(object):
+    #
+    # TODO: The following ones are actually entirely PDELab independent!
+    # They should be placed elsewhere and be used directly in the visitor.
+    #
+
+    def dimension_iname(self, context=None, count=None):
+        return dimension_iname(context=context, count=count)
+
+    def name_index(self, ind):
+        return name_index(ind)
+
+    #
+    # Local function space related generator functions
+    #
+
+    def lfs_iname(self, element, restriction, number):
+        return lfs_iname(element, restriction, number)
+
+    #
+    # Test and trial function related generator functions
+    #
+
+    def name_basis(self, element, restriction):
+        return name_basis(element, restriction)
+
+    def name_reference_gradient(self, element, restriction):
+        return name_reference_gradient(element, restriction)
+
+    def pymbolic_trialfunction_gradient(self, element, restriction, component):
+        return pymbolic_trialfunction_gradient(element, restriction, component)
+
+    def pymbolic_apply_function_gradient(self, element, restriction, component):
+        return pymbolic_apply_function_gradient(element, restriction, component)
+
+    def pymbolic_trialfunction(self, element, restriction, component):
+        return pymbolic_trialfunction(element, restriction, component)
+
+    def pymbolic_apply_function(self, element, restriction, component):
+        return pymbolic_apply_function(element, restriction, component)
+
+    #
+    # Parameter function related generator functions
+    #
+
+    def intersection_parameter_function(self, name, expr, cellwise_constant):
+        return intersection_parameter_function(name, expr, cellwise_constant)
+
+    def cell_parameter_function(self, name, expr, restriction, cellwise_constant):
+        return cell_parameter_function(name, expr, restriction, cellwise_constant)
+
+    #
+    # Geometry related generator functions
+    #
+
+    def name_facet_jacobian_determinant(self):
+        return name_facet_jacobian_determinant()
+
+    def name_jacobian_determinant(self):
+        return name_jacobian_determinant()
+
+    def name_jacobian_inverse_transposed(self, restriction):
+        return name_jacobian_inverse_transposed(restriction)
+
+    def name_unit_inner_normal(self):
+        return name_unit_inner_normal()
+
+    def name_unit_outer_normal(self):
+        return name_unit_outer_normal()
+
+    #
+    # Quadrature related generator functions
+    #
+
+    def name_quadrature_weight(self):
+        return name_quadrature_weight()
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 06d4a75dc02bfb21051eb16311ed2c7d03c0c55a..afcb0e956c1597d526251dcf12efdc86def82d8f 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -6,17 +6,15 @@ Namely:
 """
 
 from dune.perftool.options import get_option
-from dune.perftool.generation import (domain,
+from dune.perftool.generation import (cached,
+                                      domain,
                                       function_mangler,
                                       iname,
-                                      pymbolic_expr,
                                       globalarg,
                                       valuearg,
                                       get_global_context_value
                                       )
-from dune.perftool.pdelab import (name_index,
-                                  restricted_name,
-                                  )
+from dune.perftool.pdelab.index import name_index
 from dune.perftool.pdelab.basis import (evaluate_coefficient,
                                         evaluate_coefficient_gradient,
                                         name_basis,
@@ -24,7 +22,8 @@ from dune.perftool.pdelab.basis import (evaluate_coefficient,
 from dune.perftool.pdelab.spaces import (lfs_iname,
                                          name_lfs_bound,
                                          )
-from dune.perftool import Restriction
+from dune.perftool.pdelab.restriction import restricted_name
+from dune.perftool.ufl.modified_terminals import Restriction
 
 from pymbolic.primitives import Call, Subscript, Variable
 
@@ -91,46 +90,36 @@ def accumulation_mangler(target, func, dtypes):
                                   )
 
 
-def name_trialfunction_gradient(element, restriction, component):
+def pymbolic_trialfunction_gradient(element, restriction, component):
     rawname = "gradu" + "_".join(str(c) for c in component)
     name = restricted_name(rawname, restriction)
     container = name_coefficientcontainer(restriction)
-
-    # TODO
-    #
-    # This is just a temporary test used to create an A-matrix as
-    # local operator class member. Right now it doesn't evaluate
-    # anything.
-    if get_option("sumfact") and restriction == Restriction.NONE:
-        from dune.perftool.sumfact import start_sumfactorization
-        start_sumfactorization(element, container, restriction, component)
-
     evaluate_coefficient_gradient(element, name, container, restriction, component)
-    return name
+    return Variable(name)
 
 
-def name_trialfunction(element, restriction, component):
+def pymbolic_trialfunction(element, restriction, component):
     rawname = "u" + "_".join(str(c) for c in component)
     name = restricted_name(rawname, restriction)
     container = name_coefficientcontainer(restriction)
     evaluate_coefficient(element, name, container, restriction, component)
-    return name
+    return Variable(name)
 
 
-def name_apply_function_gradient(element, restriction, component):
+def pymbolic_apply_function_gradient(element, restriction, component):
     rawname = "gradz_func" + "_".join(str(c) for c in component)
     name = restricted_name(rawname, restriction)
     container = name_applycontainer(restriction)
     evaluate_coefficient_gradient(element, name, container, restriction, component)
-    return name
+    return Variable(name)
 
 
-def name_apply_function(element, restriction, component):
+def pymbolic_apply_function(element, restriction, component):
     rawname = "z_func" + "_".join(str(c) for c in component)
     name = restricted_name(rawname, restriction)
     container = name_applycontainer(restriction)
     evaluate_coefficient(element, name, container, restriction, component)
-    return name
+    return Variable(name)
 
 
 def name_coefficientcontainer(restriction):
@@ -143,7 +132,7 @@ def name_applycontainer(restriction):
     return name
 
 
-@pymbolic_expr
+@cached
 def pymbolic_coefficient(container, lfs, index):
     # TODO introduce a proper type for local function spaces!
     if isinstance(lfs, str):
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index ced7d2da58613b7ac8e6d583861ff0220fe3ace8..ad587091096fff8a51debf2002332da2412890bf 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -1,6 +1,7 @@
 """ Generators for basis evaluations """
 
-from dune.perftool.generation import (cached,
+from dune.perftool.generation import (backend,
+                                      cached,
                                       class_member,
                                       generator_factory,
                                       include_file,
@@ -27,8 +28,7 @@ from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs,
                                                 lop_template_test_gfs,
                                                 )
 from dune.perftool.pdelab.driver import FEM_name_mangling
-
-from dune.perftool.pdelab import restricted_name
+from dune.perftool.pdelab.restriction import restricted_name
 from pymbolic.primitives import Product, Subscript, Variable
 from loopy import Reduction
 
@@ -124,6 +124,7 @@ def shape_as_pymbolic(shape):
     return tuple(_shape_as_pymbolic(s) for s in shape)
 
 
+@backend(interface="eval_coefficient")
 @cached
 def evaluate_coefficient(element, name, container, restriction, component):
     from ufl.functionview import select_subelement
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index f4d476bc6d2c1f07291e12a6b8d88393f029e9ae..8dd41c1237c43e76d11a10d4b94365d98fbdd472 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -1,5 +1,5 @@
-from dune.perftool import Restriction
-from dune.perftool.pdelab import restricted_name
+from dune.perftool.ufl.modified_terminals import Restriction
+from dune.perftool.pdelab.restriction import restricted_name
 from dune.perftool.generation import (cached,
                                       domain,
                                       get_global_context_value,
diff --git a/python/dune/perftool/pdelab/index.py b/python/dune/perftool/pdelab/index.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc87edfdcb3ff2be39ec3fe4fa136c49e389b5cf
--- /dev/null
+++ b/python/dune/perftool/pdelab/index.py
@@ -0,0 +1,17 @@
+from dune.perftool.generation import cached
+
+from ufl.classes import MultiIndex, Index
+
+
+# Now define some commonly used generators that do not fall into a specific category
+@cached
+def name_index(index):
+    if isinstance(index, Index):
+        # This failed for index > 9 because ufl placed curly brackets around
+        # return str(index)
+        return "i_{}".format(index.count())
+    if isinstance(index, MultiIndex):
+        assert len(index) == 1
+        # return str(index._indices[0])
+        return "i_{}".format(index._indices[0].count())
+    raise NotImplementedError
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 3c008a173b98bcfcd06fb77f2948d4b6a76e2462..5f9ac0d9bd952d9e21cbba86a4db53477be97944 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -20,7 +20,7 @@ from dune.perftool.cgen.clazz import (AccessModifier,
                                       BaseClass,
                                       ClassMember,
                                       )
-from dune.perftool import Restriction
+from dune.perftool.ufl.modified_terminals import Restriction
 from pymbolic.primitives import Variable
 from pytools import Record
 
@@ -355,7 +355,7 @@ def boundary_predicates(expr, measure, subdomain_id):
 
 @iname
 def grad_iname(index, dim):
-    from dune.perftool.pdelab import name_index
+    from dune.perftool.pdelab.index import name_index
     name = name_index(index)
     domain(name, dim)
     return name
@@ -403,11 +403,13 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
                 (ansatz_lfs.get_args() + test_lfs.get_args() + (pymbolic_expr,))
                 )
 
-    from dune.perftool.generation import instruction
-    from dune.perftool.pdelab.quadrature import quadrature_iname
+    from dune.perftool.generation import get_backend, instruction
+    from dune.perftool.options import option_switch
+    quad_inames = get_backend(interface="quadinames", selector=option_switch("sumfac"))()
+
     instruction(assignees=(),
                 expression=expr,
-                forced_iname_deps=additional_inames.union(frozenset(visitor.inames).union(frozenset({quadrature_iname()}))),
+                forced_iname_deps=additional_inames.union(frozenset(visitor.inames).union(quad_inames)),
                 forced_iname_deps_is_final=True,
                 predicates=predicates
                 )
@@ -457,8 +459,10 @@ def generate_kernel(integrals):
                     indexmap[j] = indexmap[i]
 
             # Get a transformer instance for this kernel
+            from dune.perftool.pdelab import PDELabInterface
+            interface = PDELabInterface()
             from dune.perftool.ufl.visitor import UFL2LoopyVisitor
-            visitor = UFL2LoopyVisitor(measure, indexmap)
+            visitor = UFL2LoopyVisitor(interface, measure, indexmap)
             generate_accumulation_instruction(visitor, term, measure, subdomain_id)
 
     # Extract the information, which is needed to create a loopy kernel.
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index bf9431779060dbb33bd1c49a9cd84c2a4aad9441..b50c8248512ce363cb35b735c199698102ddf110 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -1,5 +1,5 @@
-from dune.perftool import Restriction
-from dune.perftool.generation import (cached,
+from dune.perftool.generation import (backend,
+                                      cached,
                                       domain,
                                       get_global_context_value,
                                       iname,
@@ -8,6 +8,7 @@ from dune.perftool.generation import (cached,
                                       temporary_variable,
                                       )
 from dune.perftool.options import get_option
+from dune.perftool.ufl.modified_terminals import Restriction
 
 
 @iname
@@ -16,6 +17,11 @@ def quadrature_iname():
     return "q"
 
 
+@backend(interface="quadinames")
+def quadrature_inames():
+    return frozenset({quadrature_iname()})
+
+
 def quadrature_preamble(code, **kw):
     return instruction(inames=quadrature_iname(), code=code, **kw)
 
diff --git a/python/dune/perftool/pdelab/restriction.py b/python/dune/perftool/pdelab/restriction.py
new file mode 100644
index 0000000000000000000000000000000000000000..e889fcd64c743e6cc301aa46fa1d71162c2662d7
--- /dev/null
+++ b/python/dune/perftool/pdelab/restriction.py
@@ -0,0 +1,10 @@
+from dune.perftool.ufl.modified_terminals import Restriction
+
+
+def restricted_name(name, restriction):
+    if restriction == Restriction.NONE:
+        return name
+    if restriction == Restriction.POSITIVE:
+        return name + '_n'
+    if restriction == Restriction.NEGATIVE:
+        return name + '_s'
\ No newline at end of file
diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py
index 3a40e636de6ba3224fdd1aaab2102419d3b4f1e6..b329245d2d6d2c5e04f60cc4ebcdca1439df9666 100644
--- a/python/dune/perftool/pdelab/signatures.py
+++ b/python/dune/perftool/pdelab/signatures.py
@@ -1,6 +1,6 @@
 """ Signatures for PDELab local opreator assembly functions """
 
-from dune.perftool import Restriction
+from dune.perftool.ufl.modified_terminals import Restriction
 from dune.perftool.pdelab.geometry import (name_geometry_wrapper,
                                            type_geometry_wrapper,
                                            )
diff --git a/python/dune/perftool/pdelab/spaces.py b/python/dune/perftool/pdelab/spaces.py
index e439425f36e8964d407814a2042337578b12dc34..07403a8641ba8701ccd3a89a8c9c619ef8167df9 100644
--- a/python/dune/perftool/pdelab/spaces.py
+++ b/python/dune/perftool/pdelab/spaces.py
@@ -6,7 +6,7 @@ from dune.perftool.generation import (domain,
                                       include_file,
                                       preamble,
                                       )
-from dune.perftool.pdelab import restricted_name
+from dune.perftool.pdelab.restriction import restricted_name
 
 from loopy import CallMangleInfo
 from loopy.symbolic import FunctionIdentifier
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index 793cbde027e0e7c6a9b8278adbfa66cd60d38156..4f1ce7c9d8dbe7946c664dc9ebe2f1b51fd340f7 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -1 +1,7 @@
+# Trigger some imports that are needed to have all backend implementations visible
+# to the selection mechanisms
+import dune.perftool.sumfact.amatrix
+import dune.perftool.sumfact.quadrature
+import dune.perftool.sumfact.sumfact
+
 from dune.perftool.sumfact.sumfact import start_sumfactorization
diff --git a/python/dune/perftool/sumfact/amatrix.py b/python/dune/perftool/sumfact/amatrix.py
index 9d2e249886c0075fe318449cfede12183ec24d2d..968f4f16de04ce7daace2e19f74185ff24380222 100644
--- a/python/dune/perftool/sumfact/amatrix.py
+++ b/python/dune/perftool/sumfact/amatrix.py
@@ -1,4 +1,4 @@
-from dune.perftool import Restriction
+from dune.perftool.ufl.modified_terminals import Restriction
 
 from dune.perftool.options import get_option
 
diff --git a/python/dune/perftool/sumfact/sumfact.py b/python/dune/perftool/sumfact/sumfact.py
index a7ee0a042c1fab3ac5b3662877c5a4b3a841f8a4..5a3e364d2705281f8184ceb1ab4b013956bba50a 100644
--- a/python/dune/perftool/sumfact/sumfact.py
+++ b/python/dune/perftool/sumfact/sumfact.py
@@ -1,5 +1,6 @@
 from dune.perftool.pdelab.argument import pymbolic_coefficient
-from dune.perftool.generation import (domain,
+from dune.perftool.generation import (backend,
+                                      domain,
                                       get_counter,
                                       iname,
                                       instruction,
@@ -60,7 +61,8 @@ def setup_theta(element, container, restriction, component, a_matrices):
 
 
 # TODO this code is WIP and mainly used for experiments.
-def start_sumfactorization(element, container, restriction, component):
+@backend(interface="eval_coefficient", name="sumfact")
+def start_sumfactorization(element, name, container, restriction, component):
     theta = name_theta()
     rows = quadrature_points_per_direction()
     cols = basis_functions_per_direction()
diff --git a/python/dune/perftool/ufl/modified_terminals.py b/python/dune/perftool/ufl/modified_terminals.py
index 33319b273143455b599bcc1665ea4ff81d35182f..4ba1d0a3d34673b19cf04a3819f12ebaee03e353 100644
--- a/python/dune/perftool/ufl/modified_terminals.py
+++ b/python/dune/perftool/ufl/modified_terminals.py
@@ -1,11 +1,16 @@
 """ A module mimicking some functionality of uflacs' modified terminals """
 
 from ufl.algorithms import MultiFunction
-from dune.perftool import Restriction
 from ufl.classes import MultiIndex
 from pytools import Record
 
 
+class Restriction:
+    NONE = 0
+    NEGATIVE = 1
+    POSITIVE = 2
+
+
 class ModifiedArgument(Record):
     def __init__(self,
                  expr=None,
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index 70d77358bdb105eeb5d09b7635214190d2517510..d87c5dcf57f395c6f6ef47359b44825557c776d8 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -3,32 +3,40 @@ This module defines the main visitor algorithm transforming ufl expressions
 to pymbolic and loopy.
 """
 
-from dune.perftool import Restriction
-from dune.perftool.ufl.modified_terminals import ModifiedTerminalTracker
 from dune.perftool.generation import (domain,
-                                      get_temporary_name,
-                                      global_context,
-                                      globalarg,
-                                      iname,
-                                      instruction,
-                                      temporary_variable,
-                                      valuearg,
+                                      get_global_context_value,
                                       )
+from dune.perftool.ufl.flatoperators import get_operands
+from dune.perftool.ufl.modified_terminals import (ModifiedTerminalTracker,
+                                                  Restriction,
+                                                  )
+from dune.perftool.ufl.execution import Expression
+
+from loopy import Reduction
+
+from pymbolic.primitives import (Call,
+                                 Product,
+                                 Quotient,
+                                 Subscript,
+                                 Sum,
+                                 Variable,
+                                 )
 
-from dune.perftool.pdelab.spaces import (lfs_iname,
-                                         name_leaf_lfs,
-                                         name_lfs,
-                                         name_lfs_bound,
-                                         traverse_lfs_tree,
-                                         )
-from dune.perftool.pdelab.quadrature import quadrature_iname
-from pymbolic.primitives import Subscript, Variable
 from ufl.algorithms import MultiFunction
+from ufl.checks import is_cellwise_constant
+from ufl.functionview import select_subelement
+from ufl import (VectorElement,
+                 TensorElement,
+                 )
+from ufl.classes import (FixedIndex,
+                         IndexSum,
+                         JacobianDeterminant,
+                         )
 
 
 class UFL2LoopyVisitor(ModifiedTerminalTracker):
-    def __init__(self, measure, dimension_indices):
-        # Some variables describing the integral measure of this integral
+    def __init__(self, interface, measure, dimension_indices):
+        self.interface = interface
         self.measure = measure
         self.dimension_indices = dimension_indices
 
@@ -58,14 +66,12 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             restriction = Restriction.NEGATIVE
 
         # Select the correct subtree of the finite element
-        from ufl.functionview import select_subelement
         element = select_subelement(o.ufl_element(), self.component)
         leaf_element = element
 
         # Now treat the case of this being a vector finite element
         if element.num_sub_elements() > 0:
             # I cannot handle general mixed elements here...
-            from ufl import VectorElement, TensorElement
             assert isinstance(element, (VectorElement, TensorElement))
 
             # Determine whether this is a non-scalar subargument. This information is later
@@ -73,9 +79,8 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             self.argshape = len(element.value_shape())
 
             # If this is a vector element, we need add an additional accumulation loop iname
-            from dune.perftool.pdelab.geometry import dimension_iname
             for i in range(self.argshape):
-                self.inames.append(dimension_iname(context='arg', count=i))
+                self.inames.append(self.interface.dimension_iname(context='arg', count=i))
 
             # For the purpose of basis evaluation, we need to take the leaf element
             leaf_element = element.sub_elements()[0]
@@ -84,16 +89,13 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             raise ValueError("Gradients should have been transformed to reference gradients!!!")
 
         # Have the issued instruction depend on the iname for this localfunction space
-        from dune.perftool.pdelab.spaces import lfs_iname
-        iname = lfs_iname(leaf_element, restriction, o.number())
+        iname = self.interface.lfs_iname(leaf_element, restriction, o.number())
         self.inames.append(iname)
 
         if self.reference_grad:
-            from dune.perftool.pdelab.basis import name_reference_gradient
-            return Subscript(Variable(name_reference_gradient(leaf_element, restriction)), (Variable(iname), 0))
+            return Subscript(Variable(self.interface.name_reference_gradient(leaf_element, restriction)), (Variable(iname), 0))
         else:
-            from dune.perftool.pdelab.basis import name_basis
-            return Subscript(Variable(name_basis(leaf_element, restriction)), (Variable(iname),))
+            return Subscript(Variable(self.interface.name_basis(leaf_element, restriction)), (Variable(iname),))
 
     def coefficient(self, o):
         # Do something different for trial function and coefficients from jacobian apply
@@ -108,40 +110,30 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
 
             if self.reference_grad:
                 if o.count() == 0:
-                    from dune.perftool.pdelab.argument import name_trialfunction_gradient
-                    return Variable(name_trialfunction_gradient(o.ufl_element(), restriction, self.component))
+                    return self.interface.pymbolic_trialfunction_gradient(o.ufl_element(), restriction, self.component)
                 else:
-                    from dune.perftool.pdelab.argument import name_apply_function_gradient
-                    return Variable(name_apply_function_gradient(o.ufl_element(), restriction, self.component))
+                    return self.interface.pymbolic_apply_function_gradient(o.ufl_element(), restriction, self.component)
             else:
                 if o.count() == 0:
-                    from dune.perftool.pdelab.argument import name_trialfunction
-                    return Variable(name_trialfunction(o.ufl_element(), restriction, self.component))
+                    return self.interface.pymbolic_trialfunction(o.ufl_element(), restriction, self.component)
                 else:
-                    from dune.perftool.pdelab.argument import name_apply_function
-                    return Variable(name_apply_function(o.ufl_element(), restriction, self.component))
+                    return self.interface.pymbolic_apply_function(o.ufl_element(), restriction, self.component)
 
         # Check if this is a parameter function
         else:
             # We expect all coefficients to be of type Expression!
-            from dune.perftool.ufl.execution import Expression
             assert isinstance(o, Expression)
 
             # Determine the name of the parameter function
-            from dune.perftool.generation import get_global_context_value
             name = get_global_context_value("data").object_names[id(o)]
 
-            from ufl.checks import is_cellwise_constant
             cellwise_constant = is_cellwise_constant(o)
 
             # Trigger the generation of code for this thing in the parameter class
-            from dune.perftool.pdelab.parameter import (cell_parameter_function,
-                                                        intersection_parameter_function,
-                                                        )
             if o.on_intersection:
-                intersection_parameter_function(name, o, cellwise_constant)
+                self.interface.intersection_parameter_function(name, o, cellwise_constant)
             else:
-                cell_parameter_function(name, o, self.restriction, cellwise_constant)
+                self.interface.cell_parameter_function(name, o, self.restriction, cellwise_constant)
 
             # And return a symbol
             return Variable(name)
@@ -185,39 +177,31 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         ind = o.ufl_operands[1][0]
         redinames = additional_inames + (ind,)
         shape = o.ufl_operands[0].ufl_index_dimensions[0]
-        from dune.perftool.pdelab import name_index
-        domain(name_index(ind), shape)
+        domain(self.interface.name_index(ind), shape)
 
         # If the left operand is an index sum to, we do it in one reduction
-        from ufl.classes import IndexSum
         if isinstance(o.ufl_operands[0], IndexSum):
             return self.index_sum(o.ufl_operands[0], additional_inames=redinames)
         else:
-            from loopy import Reduction
-
             # Recurse to get the summation expression
             term = self.call(o.ufl_operands[0])
             redinames = tuple(i for i in redinames if i not in self.dimension_indices)
             if len(redinames) > 0:
-                ret = Reduction("sum", tuple(name_index(ind) for ind in redinames), term)
+                ret = Reduction("sum", tuple(self.interface.name_index(ind) for ind in redinames), term)
             else:
                 ret = term
 
             return ret
 
     def _index_or_fixed_index(self, index):
-        from ufl.classes import FixedIndex
         if isinstance(index, FixedIndex):
             return index._value
         else:
-            from pymbolic.primitives import Variable
-            from dune.perftool.pdelab import name_index
             if index in self.dimension_indices:
-                from dune.perftool.pdelab.geometry import dimension_iname
                 self.inames.append(self.dimension_indices[index])
                 return Variable(self.dimension_indices[index])
             else:
-                return Variable(name_index(index))
+                return Variable(self.interface.name_index(index))
 
     def multi_index(self, o):
         return tuple(self._index_or_fixed_index(i) for i in o)
@@ -231,8 +215,6 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
     #
 
     def product(self, o):
-        from dune.perftool.ufl.flatoperators import get_operands
-        from pymbolic.primitives import Product
         return Product(tuple(self.call(op) for op in get_operands(o)))
 
     def float_value(self, o):
@@ -242,25 +224,18 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         return o.value()
 
     def division(self, o):
-        assert len(o.ufl_operands) == 2
-
-        from pymbolic.primitives import Quotient
         return Quotient(self.call(o.ufl_operands[0]), self.call(o.ufl_operands[1]))
 
     def sum(self, o):
-        from dune.perftool.ufl.flatoperators import get_operands
-        from pymbolic.primitives import Sum
         return Sum(tuple(self.call(op) for op in get_operands(o)))
 
     def zero(self, o):
         return 0
 
     def abs(self, o):
-        from ufl.classes import JacobianDeterminant
         if isinstance(o.ufl_operands[0], JacobianDeterminant):
             return self.call(o.ufl_operands[0])
         else:
-            from pymbolic.primitives import Call
             return Call('abs', self.call(o.ufl_operands[0]))
 
     #
@@ -271,29 +246,20 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         # The normal must be restricted to be well-defined
         assert self.restriction is not Restriction.NONE
 
-        from pymbolic.primitives import Variable
         if self.restriction == Restriction.POSITIVE:
-            from dune.perftool.pdelab.geometry import name_unit_outer_normal
-            return Variable(name_unit_outer_normal())
+            return Variable(self.interface.name_unit_outer_normal())
         if self.restriction == Restriction.NEGATIVE:
             # It is highly unnatural to have this generator function,
             # but I do run into subtle trouble with return -1*outer
             # as the indexing into the normal happens only later.
             # Not investing more time into this cornercase right now.
-            from dune.perftool.pdelab.geometry import name_unit_inner_normal
-            return Variable(name_unit_inner_normal())
-
-    def facet_area(self, o):
-        from dune.perftool.pdelab.geometry import name_facetarea
-        return Variable(name_facetarea())
+            return Variable(self.interface.name_unit_inner_normal())
 
     def quadrature_weight(self, o):
-        from dune.perftool.pdelab.quadrature import name_quadrature_weight
-        return Variable(name_quadrature_weight())
+        return Variable(self.interface.name_quadrature_weight())
 
     def jacobian_determinant(self, o):
-        from dune.perftool.pdelab.geometry import name_jacobian_determinant
-        return Variable(name_jacobian_determinant())
+        return Variable(self.interface.name_jacobian_determinant())
 
     def jacobian_inverse(self, o):
         restriction = self.restriction
@@ -301,12 +267,10 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             restriction = Restriction.NEGATIVE
 
         self.transpose_necessary = True
-        from dune.perftool.pdelab.geometry import name_jacobian_inverse_transposed
-        return Variable(name_jacobian_inverse_transposed(restriction))
+        return Variable(self.interface.name_jacobian_inverse_transposed(restriction))
 
     def jacobian(self, o):
         raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
 
     def facet_jacobian_determinant(self, o):
-        from dune.perftool.pdelab.geometry import name_facet_jacobian_determinant
-        return Variable(name_facet_jacobian_determinant())
+        return Variable(self.interface.name_facet_jacobian_determinant())
diff --git a/test/sumfact/CMakeLists.txt b/test/sumfact/CMakeLists.txt
index 6ef2765c02a5e22303da1bfe440916c3fce78d26..c3435819caf68bfd22d03a728cb97d1c75cc8290 100644
--- a/test/sumfact/CMakeLists.txt
+++ b/test/sumfact/CMakeLists.txt
@@ -1 +1,2 @@
+add_subdirectory(mass)
 add_subdirectory(poisson)
diff --git a/test/sumfact/mass/CMakeLists.txt b/test/sumfact/mass/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a3357befce6cbf08249c3c56af6f8ceb4476d0fd
--- /dev/null
+++ b/test/sumfact/mass/CMakeLists.txt
@@ -0,0 +1,5 @@
+# A mass matrix test: Simplest possible example
+dune_add_formcompiler_system_test(UFLFILE mass.ufl
+                                  BASENAME sumfact_mass
+                                  INIFILE mass.mini
+                                  )
diff --git a/test/sumfact/mass/mass.mini b/test/sumfact/mass/mass.mini
new file mode 100644
index 0000000000000000000000000000000000000000..acab6ddcf5f3dc76546726807d155488c6b78aba
--- /dev/null
+++ b/test/sumfact/mass/mass.mini
@@ -0,0 +1,13 @@
+__name = sumfact_mass_{__exec_suffix}
+__exec_suffix = numdiff, symdiff | expand num
+
+cells = 1 1
+extension = 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+sumfact = 1
diff --git a/test/sumfact/mass/mass.ufl b/test/sumfact/mass/mass.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..6434e09f62b28ca9ea003c936997233d31449431
--- /dev/null
+++ b/test/sumfact/mass/mass.ufl
@@ -0,0 +1,10 @@
+cell = "quadrilateral"
+
+V = FiniteElement("DG", cell, 1)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+r = u * v * dx
+
+forms = [r]
diff --git a/test/sumfact/poisson/CMakeLists.txt b/test/sumfact/poisson/CMakeLists.txt
index 6a7006ce4c84db13bd9376b56be9eb69770d466d..07d4ab70639052991c9b9a1d72454a73822f9403 100644
--- a/test/sumfact/poisson/CMakeLists.txt
+++ b/test/sumfact/poisson/CMakeLists.txt
@@ -3,9 +3,9 @@
 #                                   BASENAME sumfact_poisson
 #                                   INIFILE poisson.mini
 #                                   )
-
-# 2. Poisson Test Case: DG, f + pure dirichlet
-dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
-                                  BASENAME sumfact_poisson_dg
-                                  INIFILE poisson_dg.mini
-                                  )
+#
+# # 2. Poisson Test Case: DG, f + pure dirichlet
+#dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
+#                                  BASENAME sumfact_poisson_dg
+#                                  INIFILE poisson_dg.mini
+#                                  )