Hello,

I was attempting to implement a time stepping scheme, but kept getting weird solutions. After some trial and error I reduced the problem to the following:

```
from dolfin import *
import numpy as np
file = File('data/data.pvd')
mesh = UnitSquareMesh(100, 100)
element = {
"f": FiniteElement("CG", mesh.ufl_cell(), 2),
"sigma": VectorElement("CG", mesh.ufl_cell(), 2),
}
X = FunctionSpace(mesh, MixedElement(*element.values()))
initial_f = Expression('exp(-(pow(x[0]-0.5, 2) + pow(x[1]-0.5, 2))/0.1)', degree = 2)
U0 = Function(X)
f0 = project(initial_f, X.sub(0).collapse())
assign(U0.sub(0), f0)
f0, sigma = split(U0)
file << U0.split(deepcopy=True)[0]
bc = DirichletBC(X.sub(0), Constant(0), "on_boundary")
U = Function(X)
f, sigma = split(U)
phi = TestFunctions(X)
v, tau = phi
# Define the variational form
a = f * v *dx - f0 * v * dx + inner(sigma, tau) * dx
# This should be solved by f = f0 and sigma = 0
solve(a == 0, U, bcs = bc)
f0, sigma0 = split(U)
solutions_f = U.sub(0, deepcopy=True).vector()
solutions_s = U.sub(1, deepcopy=True).vector()
# Update U0
U0.vector().set_local(np.concatenate((solutions_f, solutions_s)))
f0, sigma0 = split(U0)
file << U0.split(deepcopy=True)[0]
```

As expected a solution is found in a single iteration. I would expect the solution to the variational problem to be f = f0 and sigma = 0. However, plotting the solution in Paraview results in this:

I am using version 2019.2.0.dev0 of dolfin, installed using the Ubuntu ppa. I should also note that this does not happen, when I use a regular Function Space. This only happens when I introduce the mixed function space.