I’m following the demo for the Stokes equations with an iterative solver and I can get that to run. However, if I try to reformulate the problem as below, it appears that my velocity boundary condition gets applied correctly on the elements that share a face with the boundary, but the pressure and the remainder of the velocity solution remain close to 0. My question is, is there a way to modify the existing code below to solve the problem correctly?

```
from fenics import *
mesh = UnitCubeMesh(25, 25, 25)
CG2 = VectorElement("CG", mesh.ufl_cell(), 3)
CG1 = FiniteElement("CG", mesh.ufl_cell(), 2)
TH = MixedElement([CG2, CG1])
W = FunctionSpace(mesh, TH)
noslip = DirichletBC(W.sub(0), (0.0, 0.0, 0.0), "on_boundary && ( x[1] < DOLFIN_EPS | x[1] > 1.0 - DOLFIN_EPS | x[2] < DOLFIN_EPS | x[2] > 1 - DOLFIN_EPS)")
inflow = DirichletBC(W.sub(0), (1.0, 0.0, 0.0), "x[0] < DOLFIN_EPS")
ambient_pressure = DirichletBC(W.sub(1), 0.0, "x[0] > 1 - DOLFIN_EPS")
bcs = [noslip,
inflow,
ambient_pressure]
mu = Constant(0.89e-3)
rho = Constant(997.0)
f = Constant((0.0, 0.0, 0.0))
nu = Constant(mu.values()[0]/rho.values()[0])
w = Function(W)
(u_vel, u_pre) = split(w)
(v_vel, v_pre) = TestFunctions(W)
F = nu*inner(grad(u_vel), grad(v_vel))*dx + div(v_vel)*u_pre*dx + v_pre*div(u_vel)*dx - inner(f, v_vel)*dx
solve(F == 0, w, bcs, solver_parameters={"newton_solver":
{"relative_tolerance": 1e-4,
'linear_solver': 'gmres',
'preconditioner': 'amg'}})
u_vel, u_pre = w.split(deepcopy=True)
# Store velocity field data
ufile = File(f"test/velocity_stokes.pvd")
ufile << u_vel
# Store pressure field data
pfile = File(f"test/pressure_stokes.pvd")
pfile << u_pre
```