Please find below the code I have written for this purpose. The weak variational formulation is provided above, in the original post. Since I do not know how to specify the boundary conditions for vorticity, I am obtaining a zero value for the solution everywhere. Any help is greatly appreciated.

```
## Solving Stokes flow in a rectangular pipe
## using streamfunction-vorticity formulation
from dolfin import *
import numpy as np
# Define domain and discretization
delta = 0.5
length = 10.0
diameter = 4.0
# Create mesh
p0 = Point(np.array([0.0, 0.0]))
p1 = Point(np.array([length, diameter]))
nx = int(length/delta)
ny = int(diameter/delta)
mesh = RectangleMesh(p0, p1, nx, ny)
# Making a mark dictionary
# Note: the values should be UNIQUE identifiers.
mark = {"generic": 0,
"left": 1,
"right": 2,
"top" : 3,
"bottom" : 4 }
subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(mark["generic"])
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], length)
class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1],diameter)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1],0)
top = Top()
top.mark(subdomains, mark["top"])
bottom = Bottom()
bottom.mark(subdomains, mark["bottom"])
geom_name="geometry.pvd"
File(geom_name) << subdomains
# Define function spaces
PSI = FiniteElement("Lagrange", triangle, 1)
OMEGA = FiniteElement("Lagrange", triangle, 1)
W = FunctionSpace(mesh, PSI*OMEGA)
# Define variational problem
(psi, omega) = TrialFunctions(W)
(v, q) = TestFunctions(W)
# Surface normal
n = FacetNormal(mesh)
# Pressures. First define the numbers (for later use):
omega_top = (diameter/4.)
omega_bottom = -(diameter/4.)
# ...and then the DOLFIN constants:
om_top = Constant(omega_top)
om_bottom = Constant(omega_bottom)
# Body force:
force = Constant((0.0))
a = inner(grad(psi), grad(v))*dx - (omega*v)*dx + inner(grad(omega), grad(q))*dx
L = inner(force, v) * dx
##Applying boundary conditions
noslip = Constant((0.0))
bc_psi_top = DirichletBC(W.sub(0), noslip, subdomains, mark["top"])
bc_psi_bottom = DirichletBC(W.sub(0), noslip, subdomains, mark["bottom"])
#bc_omega_top = DirichletBC(W.sub(1), om_top, subdomains, mark["top"])
#bc_omega_bottom = DirichletBC(W.sub(1), om_bottom, subdomains, mark["bottom"])
bcs = [bc_psi_top, bc_psi_bottom]
#
## Solve variational problem
w = Function(W)
solve(a == L, w, bcs)
# Split using deep copy
(psi, omega) = w.split(True)
## Save solution in VTK format
psi_file = File("psi.pvd")
psi_file << psi
omg_file = File("omega.pvd")
omg_file << omega
vel_x=project(psi.dx(1),W.sub(0).collapse())
vel_y=-project(psi.dx(0),W.sub(0).collapse())
vel_profile=as_vector((vel_x,vel_y))
vel_x_file = File("velocity_x.pvd")
vel_x_file << vel_x
#vel_y_file = File("velocity_y.pvd")
#vel_y_file << vel_y
```

Solution for \psi as obtained from code: