### How do I integrate including the facet normal with quadrilateral elements

Let $\Omega$Ω be the unit square. I wish to compute the flux of a finite element function, $u\in V^{h,p}$`u`∈`V`^{h,p} , on the top boundary,

$\int_{\partial\Omega_{\text{top}}}\nabla u\cdot n\mathrm{\mathrm{\quad d}}s$∫_{∂Ωtop}∇`u`·`n` `d``s` where $\partial\Omega_{\text{top}}:=\left[0,1\right]\times\left\{1\right\}$∂Ω_{top}:=[0,1]×{1}. In my understanding, this is equivalent to computing the functional $\int_{\partial\Omega_{\text{top}}}\nabla u\cdot\left(0,1\right)^{\top}\quad\mathrm{d}s$∫_{∂Ωtop}∇`u`·(0,1)^{⊤} `d``s` . I am experiencing unexpected results, am I doing something wrong?

Using a triangle mesh, the results are as expected.

Using a quad mesh yields results which do not match that expected when using the facet normal.

Is this issue related to https://bitbucket.org/fenics-project/ffc/pull-requests/91/fail-on-compilation-of-dq-element/diff ?

Example code follows:

```
from dolfin import *
for ctype in [CellType.Type.triangle, CellType.Type.quadrilateral]:
mesh = RectangleMesh.create(MPI.comm_world,
[Point(0, 0), Point(1, 1)],
[16, 16], ctype)
V = FunctionSpace(mesh, "Lagrange", 1)
u = interpolate(Expression("pow(x[0], 2) + pow(x[1], 2)", element=V.ufl_element()), V)
# Mark the top edge
ff = MeshFunction("size_t", mesh, 1, 0)
CompiledSubDomain("near(x[1], 1.0)").mark(ff, 1)
ds = Measure("ds", subdomain_data=ff)
# Compute functionals
n = FacetNormal(mesh)
j = assemble(dot(grad(u), n)*ds(1))
k = assemble(dot(grad(u), Constant((0.0, 1.0)))*ds(1))
print("CellType.%s, grad(u).n: %.3e, grad(u).(0, 1): %.3e" % (ctype, j, k))
```