P2 = VectorElement("CG", tetrahedron, 2)
P1 = FiniteElement("CG", tetrahedron, 1)

u = Coefficient(P2)
v = TestFunction(P2)
du = TrialFunction(P2)
p = Coefficient(P1)

mass = div(u) * v[1] * dx
conv_x = inner( dot(u, grad(u[0])), v[0]) * dx
conv_z = inner( dot(u, grad(u[2])), v[2]) * dx
press_x = p.dx(0)*v[0] * dx
press_z = p.dx(2)*v[2] * dx
visc_x = u[0].dx(1) * v[0].dx(1) * dx
visc_z = u[2].dx(1) * v[2].dx(1) * dx

F = mass + ( conv_x + press_x + visc_x) + ( conv_z + press_z + visc_z)

J = derivative(F,u,du)