diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 48e91fdd51ebde345d1dcefdefa8e99de2c41f77..97f3ce47af0d543141672b08f1e837cbc2ff0cc6 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -115,6 +115,7 @@ class CodegenFormOptionsArray(ImmutableRecord):
     control_variable = CodegenOption(default=None, helpstr="Name of control variable in UFL file")
     block_preconditioner_diagonal = CodegenOption(default=False, helpstr="Whether this operator should implement the diagonal part of a block preconditioner")
     block_preconditioner_offdiagonal = CodegenOption(default=False, helpstr="Whether this operator should implement the off-diagonal part of a block preconditioner")
+    block_preconditioner_pointdiagonal = CodegenOption(default=False, helpstr="Whether this operator should implement the point diagonal part of a block preconditioner")
     geometry_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for geometries. Currently implemented mixins: generic, axiparallel, equidistant, sumfact_multilinear, sumfact_axiparallel, sumfact_equidistant")
     quadrature_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for quadrature. Currently implemented: generic, sumfact")
     basis_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for basis function evaluation. Currently implemented: generic, sumfact")
@@ -227,10 +228,16 @@ def process_form_options(opt, form):
     if opt.filename is None:
         opt = opt.copy(filename="{}_{}_file.hh".format(get_option("target_name"), opt.classname))
 
+    if opt.block_preconditioner_pointdiagonal:
+        opt = opt.copy(generate_jacobians=False,
+                       basis_mixins="sumfact_pointdiagonal",
+                       accumulation_mixins="sumfact_pointdiagonal",
+                       )
+
     if opt.block_preconditioner_diagonal or opt.block_preconditioner_offdiagonal:
         assert opt.numerical_jacobian is False
         opt = opt.copy(generate_residuals=False,
-                       generate_jacobians=False,
+                       generate_jacobians=True,
                        matrix_free=True,
                        )
 
diff --git a/python/dune/codegen/pdelab/argument.py b/python/dune/codegen/pdelab/argument.py
index c3cc48298c1b2efb70c6c7716b29b221486044e1..5124c77608fa78fe1ea5f72a75cc7c097ef178d9 100644
--- a/python/dune/codegen/pdelab/argument.py
+++ b/python/dune/codegen/pdelab/argument.py
@@ -19,6 +19,7 @@ from dune.codegen.pdelab.spaces import (lfs_iname,
                                         )
 from dune.codegen.pdelab.restriction import restricted_name
 from dune.codegen.ufl.modified_terminals import Restriction
+from dune.codegen.options import get_form_option
 
 from pymbolic.primitives import Call, Subscript, Variable
 
@@ -144,6 +145,8 @@ def name_accumulation_variable(restrictions=None):
                 restrictions = (Restriction.NONE,)
             else:
                 restrictions = (Restriction.POSITIVE,)
+        if get_form_option("block_preconditioner_pointdiagonal"):
+            restrictions = restrictions[:1]
         return name_residual(*restrictions)
     if ft == 'jacobian':
         if restrictions is None:
diff --git a/python/dune/codegen/pdelab/driver/__init__.py b/python/dune/codegen/pdelab/driver/__init__.py
index b60544c1c78242f1490c76c46d4be6a2c4448501..50effdf1f6730c8a6241292b3dd8375e7d472f4c 100644
--- a/python/dune/codegen/pdelab/driver/__init__.py
+++ b/python/dune/codegen/pdelab/driver/__init__.py
@@ -39,7 +39,10 @@ def get_form_ident():
 
 def get_form():
     data = get_global_context_value("data")
-    return data.object_by_name[get_form_option("form", get_form_ident())]
+    form = get_form_option("form")
+    if form is None:
+        form = get_form_ident()
+    return data.object_by_name[form]
 
 
 def get_dimension():
diff --git a/python/dune/codegen/pdelab/localoperator.py b/python/dune/codegen/pdelab/localoperator.py
index dfc553fbefee7fa8c2e6aa1d92ae7ea047acdd38..94d62d1cd4df704c06f1f668756be018697fe02f 100644
--- a/python/dune/codegen/pdelab/localoperator.py
+++ b/python/dune/codegen/pdelab/localoperator.py
@@ -153,6 +153,11 @@ def enum_alpha():
     return _enum_alpha(ufl_measure_to_pdelab_measure(integral_type))
 
 
+@class_member(classtag="operator")
+def enum_skeleton_twosided():
+    return "enum { doSkeletonTwoSided = true };"
+
+
 def name_initree_member():
     define_initree("_iniParams")
     return "_iniParams"
@@ -841,6 +846,14 @@ def local_operator_default_settings(operator, form):
         base_class('Dune::PDELab::InstationaryLocalOperatorDefaultMethods<{}>'
                    .format(rf), classtag="operator")
 
+    # *always* add the volume pattern, PDELab cannot handle matrices without diagonal blocks
+    with global_context(integral_type="cell"):
+        enum_pattern()
+        pattern_baseclass()
+
+    if get_form_option("block_preconditioner_diagonal") or get_form_option("block_preconditioner_pointdiagonal"):
+        enum_skeleton_twosided()
+
 
 def measure_is_enabled(measure):
     option_dict = {"cell": "enable_volume",
@@ -855,6 +868,16 @@ def generate_residual_kernels(form, original_form):
     if not get_form_option("generate_residuals"):
         return {}
 
+    if get_form_option("block_preconditioner_pointdiagonal"):
+        from ufl import derivative
+        jacform = derivative(original_form, original_form.coefficients()[0])
+
+        from dune.codegen.ufl.preprocess import preprocess_form
+        jacform = preprocess_form(jacform).preprocessed_form
+
+        from dune.codegen.ufl.transformations.blockpreconditioner import diagonal_block_jacobian
+        form = diagonal_block_jacobian(jacform)
+
     logger = logging.getLogger(__name__)
     with global_context(form_type='residual'):
         operator_kernels = {}
@@ -957,6 +980,12 @@ def generate_jacobian_kernels(form, original_form):
                     from dune.codegen.pdelab.signatures import assembler_routine_name
                     with global_context(kernel=assembler_routine_name()):
                         kernel = [k for k in generate_kernels_per_integral(jac_apply_form.integrals_by_type(measure))]
+
+                        if kernel:
+                            enum_pattern()
+                            pattern_baseclass()
+                            enum_alpha()
+
                 operator_kernels[(measure, 'jacobian_apply')] = kernel
 
                 # Generate dummy functions for those kernels, that vanished in the differentiation process
@@ -986,6 +1015,12 @@ def generate_jacobian_kernels(form, original_form):
                         from dune.codegen.pdelab.signatures import assembler_routine_name
                         with global_context(kernel=assembler_routine_name()):
                             kernel = [k for k in generate_kernels_per_integral(jacform.integrals_by_type(measure))]
+
+                            if kernel:
+                                enum_pattern()
+                                pattern_baseclass()
+                                enum_alpha()
+
                     operator_kernels[(measure, 'jacobian')] = kernel
 
                 # Generate dummy functions for those kernels, that vanished in the differentiation process
diff --git a/python/dune/codegen/pdelab/tensors.py b/python/dune/codegen/pdelab/tensors.py
index a924a39ae1bba19a94c47e5542567d8717cedf57..59bb314da1cca79a90f4284cb869a059186f78eb 100644
--- a/python/dune/codegen/pdelab/tensors.py
+++ b/python/dune/codegen/pdelab/tensors.py
@@ -22,7 +22,7 @@ def define_assembled_tensor(name, expr, visitor):
         visitor.indices = indices
         instruction(assignee=prim.Subscript(prim.Variable(name), indices),
                     expression=visitor.call(expr),
-                    forced_iname_deps=frozenset(visitor.interface.quadrature_inames()),
+                    forced_iname_deps=frozenset(visitor.quadrature_inames()),
                     depends_on=frozenset({lp.match.Tagged("sumfact_stage1")}),
                     tags=frozenset({"quad"}),
                     )
diff --git a/python/dune/codegen/sumfact/accumulation.py b/python/dune/codegen/sumfact/accumulation.py
index c1ccd13a5552980db0be3788651065dc8a39654c..90fd01fd5b90b722c0a3ee194e1da82ed44f6eab 100644
--- a/python/dune/codegen/sumfact/accumulation.py
+++ b/python/dune/codegen/sumfact/accumulation.py
@@ -55,6 +55,7 @@ import pymbolic.primitives as prim
 from loopy.symbolic import WalkMapper
 import ufl.classes as uc
 from ufl import FiniteElement, MixedElement, TensorProductElement
+import itertools
 
 
 basis_sf_kernels = generator_factory(item_tags=("basis_sf_kernels",), context_tags='kernel', no_deco=True)
@@ -161,8 +162,10 @@ class AccumulationOutput(SumfactKernelInterfaceBase, ImmutableCuttingRecord):
         if self.trial_element is None:
             return ()
         else:
-            from dune.codegen.sumfact.basis import SumfactBasisMixin
-            return SumfactBasisMixin.lfs_inames(SumfactBasisMixin(), get_leaf(self.trial_element, self.trial_element_index), self.restriction)
+            mixin = get_form_option("basis_mixins")
+            from dune.codegen.generation import construct_from_mixins
+            MixinType = construct_from_mixins(mixins=[mixin], mixintype="basis", name="MixinType")
+            return MixinType.lfs_inames(MixinType(), get_leaf(self.trial_element, self.trial_element_index), self.restriction)
 
     def realize_input(self, sf, inames, shape, vec_iname, vec_shape, buf, ftags):
         # The result of stage 2 has the correct quadrature permutation but no
@@ -351,11 +354,68 @@ class SumfactAccumulationMixin(AccumulationMixinBase):
         return get_accumulation_info(expr, self)
 
     def list_accumulation_infos(self, expr):
-        return list_accumulation_infos(expr, self)
+        return itertools.product(_gradsplitting_generator(expr, self),
+                                 _trial_generator(expr, self),
+                                 )
 
     def generate_accumulation_instruction(self, expr):
         return generate_accumulation_instruction(expr, self)
 
+    def additional_matrix_sequence(self):
+        return None
+
+    @property
+    def prohibit_jacobian(self):
+        return False
+
+
+@accumulation_mixin("sumfact_pointdiagonal")
+class SumfactPointDiagonalAccumulationMixin(SumfactAccumulationMixin):
+    def additional_matrix_sequence(self):
+        info = self.current_info[1]
+        return construct_basis_matrix_sequence(transpose=True,
+                                               derivative=info.grad_index,
+                                               facedir=get_facedir(info.restriction),
+                                               facemod=get_facemod(info.restriction),
+                                               basis_size=get_basis_size(info),
+                                               )
+
+    def get_accumulation_info(self, expr):
+        element = expr.ufl_element()
+        leaf_element = element
+        element_index = 0
+        from ufl import MixedElement
+        if isinstance(expr.ufl_element(), MixedElement):
+            element_index = self.indices[0]
+            leaf_element = element.extract_component(element_index)[1]
+
+        restriction = self.restriction
+        if self.measure == 'exterior_facet':
+            from dune.codegen.pdelab.restriction import Restriction
+            restriction = Restriction.POSITIVE
+
+        grad_index = None
+        if self.reference_grad:
+            if isinstance(expr.ufl_element(), MixedElement):
+                grad_index = self.indices[1]
+            else:
+                grad_index = self.indices[0]
+
+        return SumfactAccumulationInfo(element=expr.ufl_element(),
+                                       element_index=element_index,
+                                       restriction=restriction,
+                                       grad_index=grad_index,
+                                       )
+
+    def list_accumulation_infos(self, expr):
+        return itertools.product(_gradsplitting_generator(expr, self, number=0),
+                                 _gradsplitting_generator(expr, self, number=1),
+                                 )
+
+    @property
+    def prohibit_jacobian(self):
+        return True
+
 
 class SumfactAccumulationInfo(ImmutableRecord):
     def __init__(self,
@@ -422,9 +482,9 @@ def _get_childs(element):
             yield (i, element.extract_component(i)[1])
 
 
-def _test_generator(expr, visitor):
+def _gradsplitting_generator(expr, visitor, number=0):
     from dune.codegen.ufl.modified_terminals import extract_modified_arguments
-    ma = extract_modified_arguments(expr, argnumber=0)
+    ma = extract_modified_arguments(expr, argnumber=number)
     if len(ma) == 0:
         return
     element = ma[0].argexpr.ufl_element()
@@ -440,7 +500,8 @@ def _test_generator(expr, visitor):
     for res in restrictions:
         for ei, e in _get_childs(element):
             for grad in (None,) + tuple(range(dim)):
-                yield SumfactAccumulationInfo(element_index=ei,
+                yield SumfactAccumulationInfo(element=element,
+                                              element_index=ei,
                                               restriction=res,
                                               grad_index=grad)
 
@@ -465,13 +526,20 @@ def _trial_generator(expr, visitor):
             yield SumfactAccumulationInfo(element_index=ei, restriction=res, element=e)
 
 
-def list_accumulation_infos(expr, visitor):
-    import itertools
-    return itertools.product(_test_generator(expr, visitor), _trial_generator(expr, visitor))
+def get_basis_size(info):
+    leaf_element = info.element
+    element_index = info.element_index
+    dim = world_dimension()
+    from ufl import MixedElement
+    if isinstance(leaf_element, MixedElement):
+        leaf_element = leaf_element.extract_component(element_index)[1]
+    degree = leaf_element._degree
+    if isinstance(degree, int):
+        degree = (degree,) * dim
+    return tuple(deg + 1 for deg in degree)
 
 
 def generate_accumulation_instruction(expr, visitor):
-    dim = world_dimension()
     test_info = visitor.test_info
     trial_info = visitor.trial_info
 
@@ -480,14 +548,7 @@ def generate_accumulation_instruction(expr, visitor):
     count_quadrature_point_operations(expr)
 
     # Number of basis functions per direction
-    leaf_element = test_info.element
-    from ufl import MixedElement
-    if isinstance(leaf_element, MixedElement):
-        leaf_element = leaf_element.extract_component(test_info.element_index)[1]
-    degree = leaf_element._degree
-    if isinstance(degree, int):
-        degree = (degree,) * dim
-    basis_size = tuple(deg + 1 for deg in degree)
+    basis_size = get_basis_size(test_info)
 
     # Anisotropic finite elements are not (yet) supported by Dune
     assert(size == basis_size[0] for size in basis_size)
@@ -523,20 +584,27 @@ def generate_accumulation_instruction(expr, visitor):
         derivative=test_info.grad_index,
         facedir=get_facedir(test_info.restriction),
         facemod=get_facemod(test_info.restriction),
-        basis_size=basis_size)
+        basis_size=basis_size,
+        additional_sequence=visitor.additional_matrix_sequence())
 
     jacobian_inames = trial_info.inames
     priority = test_info.grad_index
     if priority is None:
         priority = 3
 
+    trial_element = trial_info.element
+    trial_element_index = trial_info.element_index
+    if visitor.prohibit_jacobian:
+        trial_element = None
+        trial_element_index = 0
+
     output = AccumulationOutput(matrix_sequence=matrix_sequence,
                                 accumvar=accumvar,
                                 restriction=(test_info.restriction, trial_info.restriction),
                                 test_element=test_info.element,
                                 test_element_index=test_info.element_index,
-                                trial_element=trial_info.element,
-                                trial_element_index=trial_info.element_index,
+                                trial_element=trial_element,
+                                trial_element_index=trial_element_index,
                                 )
 
     sf = SumfactKernel(matrix_sequence=matrix_sequence,
diff --git a/python/dune/codegen/sumfact/basis.py b/python/dune/codegen/sumfact/basis.py
index 757d51870aa7ccf690b9362b9b0522a4892f792b..c49e052a90faac34ad6f64b539f93b3bb1d560cb 100644
--- a/python/dune/codegen/sumfact/basis.py
+++ b/python/dune/codegen/sumfact/basis.py
@@ -235,6 +235,28 @@ class SumfactBasisMixin(GenericBasisMixin):
         return prim.Subscript(var, vsf.quadrature_index(sf, self))
 
 
+@basis_mixin("sumfact_pointdiagonal")
+class SumfactPointDiagonalBasisMixin(SumfactBasisMixin):
+    def lfs_inames(self, element, restriction, number=1):
+        return ()
+
+    def implement_basis(self, element, restriction, number):
+        info = self.current_info[number]
+        if element == info.element and restriction == info.restriction:
+            return 1
+        else:
+            return 0
+
+    def implement_reference_gradient(self, element, restriction, number):
+        index, = self.indices
+        self.indices = None
+        info = self.current_info[number]
+        if element == info.element and restriction == info.restriction and index == info.grad_index:
+            return 1
+        else:
+            return 0
+
+
 class LFSSumfactKernelInput(SumfactKernelInterfaceBase, ImmutableCuttingRecord):
     def __init__(self,
                  matrix_sequence=None,
diff --git a/python/dune/codegen/sumfact/tabulation.py b/python/dune/codegen/sumfact/tabulation.py
index 9def97eb3ba4fdae280cc65e7f12ca73164d1146..e2e4a771ec9918f960e5ea88821aa7acb885a0c2 100644
--- a/python/dune/codegen/sumfact/tabulation.py
+++ b/python/dune/codegen/sumfact/tabulation.py
@@ -50,6 +50,7 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
                  direction=None,
                  slice_size=None,
                  slice_index=None,
+                 additional_tabulation=None,
                  ):
         """
         Arguments:
@@ -61,6 +62,7 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
         direction: Direction corresponding to this matrix
         slice_size: Number of slices for this direction
         slice_index: To which slice does this belong
+        additional_tabulation: A factor to be multiplied with this basis tabulation matrix.
         """
         assert(isinstance(basis_size, int))
         ImmutableRecord.__init__(self,
@@ -71,6 +73,7 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
                                  direction=direction,
                                  slice_size=slice_size,
                                  slice_index=slice_index,
+                                 additional_tabulation=additional_tabulation,
                                  )
 
     @property
@@ -90,6 +93,9 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
         if self.slice_size is not None:
             infos.append("s{}".format(self.slice_index))
 
+        if self.additional_tabulation is not None:
+            infos.append("prod{}".format(self.additional_tabulation._shortname))
+
         return "".join(infos)
 
     def __str__(self):
@@ -122,7 +128,9 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
     def pymbolic(self, indices):
         name = str(self)
         define_theta(name, self)
-        return prim.Subscript(prim.Variable(name), indices)
+        ret = prim.Subscript(prim.Variable(name), indices)
+
+        return ret
 
     @property
     def vectorized(self):
@@ -461,22 +469,31 @@ def define_theta(name, tabmat, additional_indices=(), width=None):
     else:
         args.append(tabmat.face)
 
+    # Get right hand side of basis evaluation matrix assignment
+    expr = prim.Call(PolynomialLookup(polynomials, tabmat.derivative), tuple(args))
+
+    # Maybe multiply another matrix (needed for the very special case of assembling point diagonals)
+    if tabmat.additional_tabulation is not None:
+        expr = prim.Product((expr, prim.Call(PolynomialLookup(polynomials, tabmat.additional_tabulation.derivative), tuple(args))))
+
     instruction(assignee=prim.Subscript(prim.Variable(name), (i, j) + additional_indices),
-                expression=prim.Call(PolynomialLookup(polynomials, tabmat.derivative), tuple(args)),
+                expression=expr,
                 kernel="operator",
                 depends_on=dep,
                 )
 
 
-def construct_basis_matrix_sequence(transpose=False, derivative=None, facedir=None, facemod=None, basis_size=None):
+def construct_basis_matrix_sequence(transpose=False, derivative=None, facedir=None, facemod=None, basis_size=None, additional_sequence=None):
     dim = world_dimension()
     result = [None] * dim
+    if additional_sequence is None:
+        additional_sequence = [None] * dim
     quadrature_size = quadrature_points_per_direction()
     assert (basis_size is not None)
     if facedir is not None:
         quadrature_size = quadrature_size[:facedir] + (1,) + quadrature_size[facedir:]
 
-    for i in range(dim):
+    for i, add_seq in zip(range(dim), additional_sequence):
         onface = None
         if facedir == i:
             onface = facemod
@@ -485,6 +502,7 @@ def construct_basis_matrix_sequence(transpose=False, derivative=None, facedir=No
                                           basis_size=basis_size[i],
                                           transpose=transpose,
                                           derivative=derivative == i,
-                                          face=onface)
+                                          face=onface,
+                                          additional_tabulation=add_seq)
 
     return tuple(result)
diff --git a/python/dune/codegen/ufl/transformations/blockpreconditioner.py b/python/dune/codegen/ufl/transformations/blockpreconditioner.py
index c992af12e0cae3499f6f6694027610dfd53c535e..b3fd11249e766bdc8abec4c0d01defbf5de9e59c 100644
--- a/python/dune/codegen/ufl/transformations/blockpreconditioner.py
+++ b/python/dune/codegen/ufl/transformations/blockpreconditioner.py
@@ -21,7 +21,7 @@ class OffDiagonalBlockSwitcher(MultiFunction):
     def positive_restricted(self, o):
         self.res = Restriction.POSITIVE
         ret = self(o.ufl_operands[0])
-        self.rest = Restriction.NONE
+        self.res = Restriction.NONE
         if isinstance(ret, uc.Zero):
             return ret
         else:
@@ -55,6 +55,8 @@ class OffDiagonalBlockSwitcher(MultiFunction):
 def list_restriction_tuples(diagonal):
     if diagonal:
         yield (Restriction.NONE, Restriction.NONE)
+        yield (Restriction.POSITIVE, Restriction.POSITIVE)
+        return
 
     res = (Restriction.POSITIVE, Restriction.NEGATIVE)
     amount = 1 if diagonal else 2
diff --git a/test/sumfact/CMakeLists.txt b/test/sumfact/CMakeLists.txt
index 187aee212ed125a99bc82670fe34e680379542f7..61b81b76082cfa3c83b04d2632194d5a4a6e6447 100644
--- a/test/sumfact/CMakeLists.txt
+++ b/test/sumfact/CMakeLists.txt
@@ -2,3 +2,4 @@ add_subdirectory(hyperbolic)
 add_subdirectory(mass)
 add_subdirectory(poisson)
 add_subdirectory(stokes)
+add_subdirectory(preconditioner)
diff --git a/test/sumfact/preconditioner/CMakeLists.txt b/test/sumfact/preconditioner/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ea77ba5dca71ce5afdff567be1f872dc0dc604e0
--- /dev/null
+++ b/test/sumfact/preconditioner/CMakeLists.txt
@@ -0,0 +1,11 @@
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_2d.ufl
+                                  BASENAME sumfact_preconditioner_2d
+                                  INIFILE preconditioner_2d.mini
+                                  SOURCE test_preconditioner_2d.cc
+                                  )
+
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_3d.ufl
+                                  BASENAME sumfact_preconditioner_3d
+                                  INIFILE preconditioner_3d.mini
+                                  SOURCE test_preconditioner_3d.cc
+                                  )
diff --git a/test/sumfact/preconditioner/poisson_dg_2d.ufl b/test/sumfact/preconditioner/poisson_dg_2d.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..ff41f6164df46aa1909db03768aef4041011afae
--- /dev/null
+++ b/test/sumfact/preconditioner/poisson_dg_2d.ufl
@@ -0,0 +1,38 @@
+cell = "quadrilateral"
+dim = 2
+degree = 1
+
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 2*(2.-2*c)*g
+
+V = FiniteElement("DG", cell, degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+# penalty factor
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+  - f*v*dx \
+  - inner(n, avg(grad(u)))*jump(v)*dS \
+  + gamma_int*jump(u)*jump(v)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - inner(n, grad(u))*v*ds \
+  + gamma_ext*u*v*ds \
+  + theta*u*inner(grad(v), n)*ds \
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(grad(v), n)*ds
+
+exact_solution = g
diff --git a/test/sumfact/preconditioner/poisson_dg_3d.ufl b/test/sumfact/preconditioner/poisson_dg_3d.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..cc6594573762ee12e1b11206eeba07a4bd038726
--- /dev/null
+++ b/test/sumfact/preconditioner/poisson_dg_3d.ufl
@@ -0,0 +1,38 @@
+cell = hexahedron
+dim = 3
+degree = 1
+
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2 + (0.5-x[2])**2
+g = exp(-1.*c)
+f = 2*(3.-2*c)*g
+
+V = FiniteElement("DG", cell, degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+# penalty factor
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+  - f*v*dx \
+  - inner(n, avg(grad(u)))*jump(v)*dS \
+  + gamma_int*jump(u)*jump(v)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - inner(n, grad(u))*v*ds \
+  + gamma_ext*u*v*ds \
+  + theta*u*inner(grad(v), n)*ds \
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(grad(v), n)*ds
+
+exact_solution = g
diff --git a/test/sumfact/preconditioner/preconditioner_2d.mini b/test/sumfact/preconditioner/preconditioner_2d.mini
new file mode 100644
index 0000000000000000000000000000000000000000..adcf214dde7018a9c5aa093cff89c51cf67d6a3a
--- /dev/null
+++ b/test/sumfact/preconditioner/preconditioner_2d.mini
@@ -0,0 +1,46 @@
+__name = preconditioner_2d
+__exec_suffix = exec
+
+cells = 2 2
+extension = 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+operators = r, blockdiag, blockoffdiag, pointdiag
+
+[formcompiler.r]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+classname = FullOperator
+filename = full_2d_operator.hh
+
+[formcompiler.blockdiag]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+block_preconditioner_diagonal = 1
+form = r
+classname = BlockDiagonalOperator
+filename = block_diagonal_2d_operator.hh
+
+[formcompiler.blockoffdiag]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+block_preconditioner_offdiagonal = 1
+form = r
+classname = BlockOffDiagonalOperator
+filename = block_offdiagonal_2d_operator.hh
+
+[formcompiler.pointdiag]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+block_preconditioner_pointdiagonal = 1
+form = r
+classname = PointDiagonalOperator
+filename = point_diagonal_2d_operator.hh
diff --git a/test/sumfact/preconditioner/preconditioner_3d.mini b/test/sumfact/preconditioner/preconditioner_3d.mini
new file mode 100644
index 0000000000000000000000000000000000000000..304d6e28f0e6f152d6327929984c0121dd55e978
--- /dev/null
+++ b/test/sumfact/preconditioner/preconditioner_3d.mini
@@ -0,0 +1,46 @@
+__name = preconditioner_3d
+__exec_suffix = exec
+
+cells = 2 2 2
+extension = 1. 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+operators = r, blockdiag, blockoffdiag, pointdiag
+
+[formcompiler.r]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+classname = FullOperator
+filename = full_3d_operator.hh
+
+[formcompiler.blockdiag]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+block_preconditioner_diagonal = 1
+form = r
+classname = BlockDiagonalOperator
+filename = block_diagonal_3d_operator.hh
+
+[formcompiler.blockoffdiag]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+block_preconditioner_offdiagonal = 1
+form = r
+classname = BlockOffDiagonalOperator
+filename = block_offdiagonal_3d_operator.hh
+
+[formcompiler.pointdiag]
+sumfact = 1
+fastdg = 1
+geometry_mixins = sumfact_equidistant
+block_preconditioner_pointdiagonal = 1
+form = r
+classname = PointDiagonalOperator
+filename = point_diagonal_3d_operator.hh
diff --git a/test/sumfact/preconditioner/test_preconditioner_2d.cc b/test/sumfact/preconditioner/test_preconditioner_2d.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2781b32399f800d1d6334565aefbede360fdbfa9
--- /dev/null
+++ b/test/sumfact/preconditioner/test_preconditioner_2d.cc
@@ -0,0 +1,129 @@
+#include "config.h"
+
+#include "dune/common/parallel/mpihelper.hh"
+#include "dune/pdelab/stationary/linearproblem.hh"
+#include "dune/pdelab/backend/istl.hh"
+#include "dune/grid/yaspgrid.hh"
+#include "dune/pdelab/finiteelementmap/qkdg.hh"
+#include "dune/pdelab/gridoperator/fastdg.hh"
+#include "dune/testtools/gridconstruction.hh"
+#include "dune/common/parametertree.hh"
+#include "dune/common/parametertreeparser.hh"
+#include <random>
+#include "dune/pdelab/gridfunctionspace/vtk.hh"
+#include "dune/grid/io/file/vtk/subsamplingvtkwriter.hh"
+#include "string"
+#include "dune/codegen/vtkpredicate.hh"
+
+// Include all the generated operators
+#include "full_2d_operator.hh"
+#include "block_diagonal_2d_operator.hh"
+#include "block_offdiagonal_2d_operator.hh"
+#include "point_diagonal_2d_operator.hh"
+
+int main(int argc, char** argv){  
+  try
+  {    
+    // Initialize basic stuff...    
+    Dune::MPIHelper& mpihelper = Dune::MPIHelper::instance(argc, argv);
+    using RangeType = double;
+    Dune::ParameterTree initree;
+    Dune::ParameterTreeParser::readINITree(argv[1], initree);
+    
+    // Setup grid (view)...    
+    using Grid = Dune::YaspGrid<2, Dune::EquidistantCoordinates<RangeType, 2>>;
+    using GV = Grid::LeafGridView;
+    using DF = Grid::ctype;
+    IniGridFactory<Grid> factory(initree);
+    std::shared_ptr<Grid> grid = factory.getGrid();
+    GV gv = grid->leafGridView();
+    
+    // Set up finite element maps...    
+    using DG2_FEM = Dune::PDELab::QkDGLocalFiniteElementMap<DF, RangeType, 1, 2>;
+    DG2_FEM dg2_fem;
+    
+    // Set up grid function spaces...    
+    using VectorBackendDG2 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>;
+    using NoConstraintsAssembler = Dune::PDELab::NoConstraints;
+    using DG2_GFS = Dune::PDELab::GridFunctionSpace<GV, DG2_FEM, NoConstraintsAssembler, VectorBackendDG2>;
+    DG2_GFS dg2_gfs_(gv, dg2_fem);
+    dg2_gfs_.name("dg2_gfs_");
+    
+    // Set up constraints container...    
+    using DG2_GFS_CC = DG2_GFS::ConstraintsContainer<RangeType>::Type;
+    DG2_GFS_CC dg2_gfs__cc;
+    dg2_gfs__cc.clear();
+    Dune::PDELab::constraints(dg2_gfs_, dg2_gfs__cc);
+    
+    // Set up grid grid operators...    
+    using FullLOP = FullOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
+    using FullGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, FullLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    FullLOP fulllop(dg2_gfs_, dg2_gfs_, initree);
+    dg2_gfs_.update();
+    MatrixBackend mb(5);
+    FullGO fullgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, fulllop, mb);
+
+    // Additional grid operators for preconditioner
+    using BDLOP = BlockDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    BDLOP bdlop(dg2_gfs_, dg2_gfs_, initree);
+    BDGO bdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bdlop, mb);
+
+    using BODLOP = BlockOffDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BODGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BODLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    BODLOP bodlop(dg2_gfs_, dg2_gfs_, initree);
+    BODGO bodgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bodlop, mb);
+
+    using PDLOP = PointDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using PDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, PDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    PDLOP pdlop(dg2_gfs_, dg2_gfs_, initree);
+    PDGO pdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, pdlop, mb);
+
+    // Set up solution vectors...    
+    using V_R = Dune::PDELab::Backend::Vector<DG2_GFS,DF>;
+    V_R x(dg2_gfs_, 0.0);
+    
+    // Testing!
+    
+    // Assemble all those matrices
+    using Dune::PDELab::Backend::native;
+    using M = typename FullGO::Traits::Jacobian;
+    M m(fullgo);
+    fullgo.jacobian(x, m);
+    Dune::printmatrix(std::cout, native(m),"full matrix","row",9,1);
+
+    using BDM = typename BDGO::Traits::Jacobian;
+    BDM bdm(bdgo);
+    bdgo.jacobian(x, bdm);
+    Dune::printmatrix(std::cout, native(bdm),"blockdiagonal matrix","row",9,1);
+
+    using BODM = typename BODGO::Traits::Jacobian;
+    BODM bodm(bodgo);
+    bodgo.jacobian(x, bodm);
+    Dune::printmatrix(std::cout, native(bodm),"blockoffdiagonal matrix","row",9,1);
+
+    V_R pd(dg2_gfs_, 0.0);
+    pdgo.residual(x, pd);
+    Dune::printvector(std::cout, native(pd), "point diagonal vector", "row");
+
+    // test failure boolean
+    bool testfail(false);
+
+    // TODO: Properly test this stuff given the above matrices.
+    //       Right now, visuals need to suffice.
+
+    // Return statement...    
+    return testfail;
+    
+  }  
+  catch (Dune::Exception& e)
+  {    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }  
+  catch (std::exception& e)
+  {    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }  
+}
+
diff --git a/test/sumfact/preconditioner/test_preconditioner_3d.cc b/test/sumfact/preconditioner/test_preconditioner_3d.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ec84a1efbf7e19f46d8531622374f2902a3139eb
--- /dev/null
+++ b/test/sumfact/preconditioner/test_preconditioner_3d.cc
@@ -0,0 +1,129 @@
+#include "config.h"
+
+#include "dune/common/parallel/mpihelper.hh"
+#include "dune/pdelab/stationary/linearproblem.hh"
+#include "dune/pdelab/backend/istl.hh"
+#include "dune/grid/yaspgrid.hh"
+#include "dune/pdelab/finiteelementmap/qkdg.hh"
+#include "dune/pdelab/gridoperator/fastdg.hh"
+#include "dune/testtools/gridconstruction.hh"
+#include "dune/common/parametertree.hh"
+#include "dune/common/parametertreeparser.hh"
+#include <random>
+#include "dune/pdelab/gridfunctionspace/vtk.hh"
+#include "dune/grid/io/file/vtk/subsamplingvtkwriter.hh"
+#include "string"
+#include "dune/codegen/vtkpredicate.hh"
+
+// Include all the generated operators
+#include "full_3d_operator.hh"
+#include "block_diagonal_3d_operator.hh"
+#include "block_offdiagonal_3d_operator.hh"
+#include "point_diagonal_3d_operator.hh"
+
+int main(int argc, char** argv){  
+  try
+  {    
+    // Initialize basic stuff...    
+    Dune::MPIHelper& mpihelper = Dune::MPIHelper::instance(argc, argv);
+    using RangeType = double;
+    Dune::ParameterTree initree;
+    Dune::ParameterTreeParser::readINITree(argv[1], initree);
+    
+    // Setup grid (view)...    
+    using Grid = Dune::YaspGrid<3, Dune::EquidistantCoordinates<RangeType, 3>>;
+    using GV = Grid::LeafGridView;
+    using DF = Grid::ctype;
+    IniGridFactory<Grid> factory(initree);
+    std::shared_ptr<Grid> grid = factory.getGrid();
+    GV gv = grid->leafGridView();
+    
+    // Set up finite element maps...    
+    using DG2_FEM = Dune::PDELab::QkDGLocalFiniteElementMap<DF, RangeType, 1, 3>;
+    DG2_FEM dg2_fem;
+    
+    // Set up grid function spaces...    
+    using VectorBackendDG2 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>;
+    using NoConstraintsAssembler = Dune::PDELab::NoConstraints;
+    using DG2_GFS = Dune::PDELab::GridFunctionSpace<GV, DG2_FEM, NoConstraintsAssembler, VectorBackendDG2>;
+    DG2_GFS dg2_gfs_(gv, dg2_fem);
+    dg2_gfs_.name("dg2_gfs_");
+    
+    // Set up constraints container...    
+    using DG2_GFS_CC = DG2_GFS::ConstraintsContainer<RangeType>::Type;
+    DG2_GFS_CC dg2_gfs__cc;
+    dg2_gfs__cc.clear();
+    Dune::PDELab::constraints(dg2_gfs_, dg2_gfs__cc);
+    
+    // Set up grid grid operators...    
+    using FullLOP = FullOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
+    using FullGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, FullLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    FullLOP fulllop(dg2_gfs_, dg2_gfs_, initree);
+    dg2_gfs_.update();
+    MatrixBackend mb(5);
+    FullGO fullgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, fulllop, mb);
+
+    // Additional grid operators for preconditioner
+    using BDLOP = BlockDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    BDLOP bdlop(dg2_gfs_, dg2_gfs_, initree);
+    BDGO bdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bdlop, mb);
+
+    using BODLOP = BlockOffDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BODGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BODLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    BODLOP bodlop(dg2_gfs_, dg2_gfs_, initree);
+    BODGO bodgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bodlop, mb);
+
+    using PDLOP = PointDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using PDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, PDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
+    PDLOP pdlop(dg2_gfs_, dg2_gfs_, initree);
+    PDGO pdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, pdlop, mb);
+
+    // Set up solution vectors...    
+    using V_R = Dune::PDELab::Backend::Vector<DG2_GFS,DF>;
+    V_R x(dg2_gfs_, 0.0);
+    
+    // Testing!
+    
+    // Assemble all those matrices
+    using Dune::PDELab::Backend::native;
+    using M = typename FullGO::Traits::Jacobian;
+    M m(fullgo);
+    fullgo.jacobian(x, m);
+    Dune::printmatrix(std::cout, native(m),"full matrix","row",9,1);
+
+    using BDM = typename BDGO::Traits::Jacobian;
+    BDM bdm(bdgo);
+    bdgo.jacobian(x, bdm);
+    Dune::printmatrix(std::cout, native(bdm),"blockdiagonal matrix","row",9,1);
+
+    using BODM = typename BODGO::Traits::Jacobian;
+    BODM bodm(bodgo);
+    bodgo.jacobian(x, bodm);
+    Dune::printmatrix(std::cout, native(bodm),"blockoffdiagonal matrix","row",9,1);
+
+    V_R pd(dg2_gfs_, 0.0);
+    pdgo.residual(x, pd);
+    Dune::printvector(std::cout, native(pd), "point diagonal vector", "row");
+
+    // test failure boolean
+    bool testfail(false);
+
+    // TODO: Properly test this stuff given the above matrices.
+    //       Right now, visuals need to suffice.
+
+    // Return statement...    
+    return testfail;
+    
+  }  
+  catch (Dune::Exception& e)
+  {    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }  
+  catch (std::exception& e)
+  {    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }  
+}
+