Something like this could perhaps do the trick?
from dolfin import *
mesh = UnitSquareMesh(4, 4)
Q = FunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh, "CG", 1)
T = TensorFunctionSpace(mesh, "CG", 1)
u = interpolate(Constant((1,2)), V)
uiuj = Function(T)
dummy = Function(Q)
u0 = Function(Q)
u1 = Function(Q)
assign(u0, u.sub(0))
assign(u1, u.sub(1))
dummy.vector()[:] = u0.vector()*u0.vector()
assign(uiuj.sub(0), dummy)
dummy.vector()[:] = u0.vector()*u1.vector()
assign(uiuj.sub(1), dummy)
assign(uiuj.sub(2), dummy) # symmetric
dummy.vector()[:] = u1.vector()*u1.vector()
assign(uiuj.sub(3), dummy)