diff --git a/test/sumfact/preconditioner/CMakeLists.txt b/test/sumfact/preconditioner/CMakeLists.txt
index 71a4460d63f1e1074505de1f964a10806859acac..ea77ba5dca71ce5afdff567be1f872dc0dc604e0 100644
--- a/test/sumfact/preconditioner/CMakeLists.txt
+++ b/test/sumfact/preconditioner/CMakeLists.txt
@@ -3,3 +3,9 @@ dune_add_formcompiler_system_test(UFLFILE poisson_dg_2d.ufl
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_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_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_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;
+ }
+}
+