diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index d36e95bf33323f370570ed7ff91d8f68264a0341..9d3dd53bf6822563b84d158278e33687222494d9 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -7,6 +7,11 @@ dune_add_formcompiler_system_test(UFLFILE stokes.ufl
                                   INIFILE stokes.mini
                                   )
 
+dune_add_formcompiler_system_test(UFLFILE stokes_quadrilateral.ufl
+                                  BASENAME stokes_quadrilateral
+                                  INIFILE stokes_quadrilateral.mini
+                                  )
+
 dune_add_formcompiler_system_test(UFLFILE stokes_sym.ufl
                                   BASENAME stokes_sym
                                   INIFILE stokes_sym.mini
diff --git a/test/stokes/stokes_quadrilateral.mini b/test/stokes/stokes_quadrilateral.mini
new file mode 100644
index 0000000000000000000000000000000000000000..0c2a459d2ebdec9279635e47c679d13aff548c89
--- /dev/null
+++ b/test/stokes/stokes_quadrilateral.mini
@@ -0,0 +1,16 @@
+__name = stokes_quadrilateral_{__exec_suffix}
+
+__exec_suffix = {diff_suffix}
+diff_suffix = numdiff, symdiff | expand num
+
+cells = 8 8
+extension = 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-11
diff --git a/test/stokes/stokes_quadrilateral.ufl b/test/stokes/stokes_quadrilateral.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..0d7c58aa180a9d4586184d1eda6c53353bbbb67a
--- /dev/null
+++ b/test/stokes/stokes_quadrilateral.ufl
@@ -0,0 +1,18 @@
+cell = quadrilateral
+
+x = SpatialCoordinate(cell)
+v_bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0))
+g_p = 8.*(1.-x[0])
+g = (g_v, g_p)
+
+P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
+P1 = FiniteElement("Lagrange", cell, 1)
+TH = P2 * P1
+
+v, q = TestFunctions(TH)
+u, p = TrialFunctions(TH)
+
+r = (inner(grad(v), grad(u)) - div(v)*p - q*div(u))*dx
+
+forms = [r]