diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index 9d3dd53bf6822563b84d158278e33687222494d9..dc0e4d36d6d3fb94f6fbaf0643c428850c36e4ec 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -22,6 +22,11 @@ dune_add_formcompiler_system_test(UFLFILE stokes_dg.ufl
                                   INIFILE stokes_dg.mini
                                   )
 
+dune_add_formcompiler_system_test(UFLFILE stokes_dg_quadrilateral.ufl
+                                  BASENAME stokes_dg_quadrilateral
+                                  INIFILE stokes_dg_quadrilateral.mini
+                                  )
+
 dune_add_formcompiler_system_test(UFLFILE stokes_stress.ufl
                                   BASENAME stokes_stress
                                   INIFILE stokes_stress.mini
diff --git a/test/stokes/stokes_dg_quadrilateral.mini b/test/stokes/stokes_dg_quadrilateral.mini
new file mode 100644
index 0000000000000000000000000000000000000000..e8b76346031f7993fb5e95681c1949b29979b574
--- /dev/null
+++ b/test/stokes/stokes_dg_quadrilateral.mini
@@ -0,0 +1,15 @@
+__name = stokes_dg_quadrilateral_{__exec_suffix}
+__exec_suffix = symdiff, numdiff | expand num
+
+cells = 8 8
+extension = 1. 1.
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-8
diff --git a/test/stokes/stokes_dg_quadrilateral.ufl b/test/stokes/stokes_dg_quadrilateral.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..e3c63e85a3ace9ae58545f9839147e752429ce7d
--- /dev/null
+++ b/test/stokes/stokes_dg_quadrilateral.ufl
@@ -0,0 +1,37 @@
+cell = quadrilateral
+
+x = SpatialCoordinate(cell)
+g_v = as_vector((4*x[1]*(1.-x[1]), 0.0))
+g_p = 8*(1.-x[0])
+g = (g_v, g_p)
+bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+
+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(grad(u)*n, v)*ds \
+  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + sigma * inner(jump(u), jump(v))*dS \
+  + sigma * inner(u-g_v, v)*ds \
+  - p*div(v)*dx \
+  - avg(p)*inner(jump(v), n)*dS \
+  + p*inner(v, n)*ds \
+  - q*div(u)*dx \
+  - avg(q)*inner(jump(u), n)*dS \
+  + q*inner(u, n)*ds \
+  - q*inner(g_v, n)*ds
+
+forms = [r]