Vector function spaces

Assembling into vector function spaces is currently unsupported. There are a few reasons why:

1. The `ufl_element` on a `VectorFunctionSpace` needs to be a `VectorElement`.

2. If this is done, the extruded code breaks, I think because ufl doesn't know how to build a `MixedElement` (of which a `VectorElement` is a special case) on top of an `OuterProductCell`. This is a problem for a combination of Andrew and Doru I think.

3. Ignoring the extrusion failures, our ffc fork currently generates the wrong code for assembling into vector spaces. I think this is a known problem (Florian?)

Corollaries:

1. Vector `Expression`s don't work, that is: u.interpolate(Expression(("0", "1"))) doesn't work.  Given all the other bits and pieces being in place, I think that is fixed by the following:

````patch
@@ -1008,6 +1012,9 @@ the :class:`FunctionSpace`.
 
         to_pts = []
 
+        if expression.rank() != self.function_space().rank:
+            raise RuntimeError('Rank mismatch between Expression and FunctionSpace')
+
         for dual in to_element.dual_basis():
             if not isinstance(dual, FIAT.functional.PointEvaluation):
                 raise NotImplementedError("Can only interpolate onto point evaluation operators. Try projecting instead")
@@ -1019,8 +1026,9 @@ the :class:`FunctionSpace`.
         # Produce C array notation of X.
         X_str = "{{"+"},\n{".join([ ",".join(map(str,x)) for x in X.T])+"}}"
 
+        assign_expression = ";\n".join(["A[%(i)d] = %(code)s" % { 'i': i, 'code': code } for i, code in enumerate(expression.code)])
         _expression_template = """
-void expression_kernel(double A[1], double **x_, int k)
+void expression_kernel(double A[%(rank)d], double **x_, int k)
 {
   const double X[%(ndof)d][%(xndof)d] = %(x_array)s;
 
@@ -1034,15 +1042,15 @@ void expression_kernel(double A[1], double **x_, int k)
     };
   };
 
-  // Note this currently only does the scalar case.
-  A[0] = %(expression)s;
+  %(assign_expression)s;
 }
 """
         kernel = op2.Kernel(_expression_template % { "x_array" : X_str,
                                                      "dim" : coords_space.dim,
                                                      "xndof" : coords_element.space_dimension(),
                                                      "ndof" : to_element.space_dimension(),
-                                                     "expression" : expression.code },
+                                                     "assign_expression" : assign_expression,
+                                                     "rank" : expression.rank() },
                             "expression_kernel")
 
         op2.par_loop(kernel, self.cell_set,
diff --git a/python/firedrake/expression.py b/python/firedrake/expression.py
index ea8e2e6..56b558a 100644
--- a/python/firedrake/expression.py
+++ b/python/firedrake/expression.py
@@ -1,11 +1,20 @@
 import ufl
-
+import numpy as np
 
 class Expression(ufl.Coefficient):
 
     def __init__(self, code=None, element=None, cell=None, degree=None, **kwargs):
 
-        self.code = code
+        self.code = np.array(code)
+        shape = self.code.shape
+        self._rank = len(shape)
+        self._shape = shape
         self.element = element
         self.cell = cell
         self.degree = degree
+
+    def rank(self):
+        return self._rank
+
+    def shape(self):
+        return self._shape
````