diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index b1fde82b59a6fbd2d471e63b547c6f4b03649d16..5e60f0276eabbfc1d1113710c2fab2cb3643d7dc 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -14,3 +14,8 @@ add_generated_executable(UFLFILE stokes_dg.ufl
                          TARGET stokes_dg_numdiff
                          FORM_COMPILER_ARGS --numerical-jacobian
                          )
+
+dune_add_system_test(TARGET stokes_dg_numdiff
+                     INIFILE stokes_dg_numdiff.mini
+                     SCRIPT dune_vtkcompare.py
+                     )
diff --git a/test/stokes/hagenpoiseuille_dg_ref.vtu b/test/stokes/hagenpoiseuille_dg_ref.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..bb365f28eab2b6b5a792a503df57b7d335c2dcf6
--- /dev/null
+++ b/test/stokes/hagenpoiseuille_dg_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a8fe239de245d906b2c6cd6be54d0f57c93cc93928a94e416f210794e2b81d80
+size 91807
diff --git a/test/stokes/stokes_dg.ufl b/test/stokes/stokes_dg.ufl
index 71ff6674ec3f1c09b115c11ecc748a39d0c311c5..7b15e24ce48c5aac90f274581d3e13f44ae80fb5 100644
--- a/test/stokes/stokes_dg.ufl
+++ b/test/stokes/stokes_dg.ufl
@@ -1,27 +1,33 @@
+g = VectorExpression("Dune::FieldVector<double, 2> y(0.0); y[0]=4*x[1]*(1.-x[1]); return y;", on_intersection=True)
+bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
+
 cell = triangle
-P2 = VectorElement("Lagrange", cell, 2)
-P1 = FiniteElement("Lagrange", cell, 1)
+P2 = VectorElement("DG", cell, 2)
+P1 = FiniteElement("DG", cell, 1)
 TH = P2 * P1
 
 v, q = TestFunctions(TH)
 u, p = TrialFunctions(TH)
 
+ds = ds(subdomain_id=1, subdomain_data=bctype)
+
 n = FacetNormal(cell)('+')
 eps = -1.0
 sigma = 1.0
 
 r = inner(grad(u), grad(v))*dx \
-  - inner(avg(grad(u))*n, jump(v))*dS \
-  + eps * inner(avg(grad(v))*n, jump(u))*dS \
+  + inner(avg(grad(u))*n, jump(v))*dS \
+  - eps * inner(avg(grad(v))*n, jump(u))*dS \
   - inner(grad(u)*n, v)*ds \
-  + eps * inner(grad(v)*n, u)*ds \
+  + eps * inner(grad(v)*n, u-g)*ds \
   + sigma * inner(jump(u), jump(v))*dS \
-  + sigma * inner(u, v)*ds \
+  + sigma * inner(u-g, v)*ds \
   - p*div(v)*dx \
-  - q*div(u)*dx \
-  + avg(p)*inner(jump(v), n)*dS \
-  + avg(q)*inner(jump(u), n)*dS \
+  - avg(p)*inner(jump(v), n)*dS \
   + p*inner(v, n)*ds \
-  + q*inner(u, n)*ds
+  - q*div(u)*dx \
+  - avg(q)*inner(jump(u), n)*dS \
+  + q*inner(u, n)*ds \
+  - q*inner(g, n)*ds
 
 forms = [r]
diff --git a/test/stokes/stokes_dg_numdiff.mini b/test/stokes/stokes_dg_numdiff.mini
new file mode 100644
index 0000000000000000000000000000000000000000..86749cea0190742b1621835922f6469acc463f2a
--- /dev/null
+++ b/test/stokes/stokes_dg_numdiff.mini
@@ -0,0 +1,12 @@
+__name = stokes_dg_numdiff
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 16 16
+elementType = simplical
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = hagenpoiseuille_dg_ref
+extension = vtu