Thank you for the prompt response. Given below is the complete code with the mesh definitions and the definitions of “r” and “theta”, and executing it would reproduce the error mentioned above.

```
from dolfin import *
from mshr import *
import math
import numpy as np
Re = 5.
Ri = 1.
rect = Rectangle(Point(0., 0.), Point(Re, Re))
num_segments=200
domain = Circle(Point(0., 0.), Re, num_segments) - Circle(Point(0., 0.), Ri, num_segments)
#creating semi-circular domain
domain = (domain - Rectangle(Point(0., 0.), Point(Re, -Re))
- Rectangle(Point(0., 0.), Point(-Re, -Re)))
mesh_resolution = 128
mesh = generate_mesh(domain, mesh_resolution)
x = SpatialCoordinate(mesh)
strm = FiniteElement("CG",triangle,2)
vort = FiniteElement("CG",triangle,2)
W = FunctionSpace(mesh,strm*vort)
r = Expression("sqrt(x[0]*x[0]+x[1]*x[1])", degree=2)
theta = Expression("atan2(x[1],x[0])", degree=2)
U = Constant(1.0)
## Constructing known solution for \psi and vorticity and applying on boundaries
psi_mfg = Expression("U*sin(theta)*sin(theta)*(0.5*r*r-0.75*r+(0.25/r))",r=r,
theta=theta,U=U,degree=2)
omeg_mfg = Expression("-1.5*U*sin(theta)/(r*r)",r=r,
theta=theta,U=U,degree=2)
### Defining boundary conditions
### and marking subdomains
class Inside(SubDomain):
def inside(self, x, on_boundary):
tol = 1.5*Ri * (1 - np.cos(np.pi / num_segments))
return near(sqrt(x[0]**2 + x[1]**2), Ri, tol) and on_boundary
class Outside(SubDomain):
def inside(self, x, on_boundary):
tol = 1.5*Re * (1 - np.cos(np.pi / num_segments))
return near(sqrt(x[0]**2 + x[1]**2), Re, tol) and on_boundary
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
Inside().mark(sub_domains, 1)
Outside().mark(sub_domains, 2)
## Assembling the BCs
inner_circ_psi = DirichletBC(W.sub(0), psi_mfg, sub_domains, 1)
outer_circ_psi = DirichletBC(W.sub(0), psi_mfg, sub_domains, 2)
inner_circ_omeg = DirichletBC(W.sub(1), omeg_mfg, sub_domains, 1)
outer_circ_omeg = DirichletBC(W.sub(1), omeg_mfg, sub_domains, 2)
bcs = [inner_circ_psi, outer_circ_psi, inner_circ_omeg, outer_circ_omeg]
### Defining test and trial functions
(psi, omeg) = TrialFunctions(W)
(u, v) = TestFunctions(W)
## Assembling bilinear equation
a1 = (inner(grad(omeg),grad(v))-(Constant(1.0)/(r*sin(theta)))*v*(omeg.dx(1))
+(omeg*v/(r*r*sin(theta)*sin(theta))))
a2 = (inner(grad(psi),grad(u))+(Constant(1.0)/(r*sin(theta)))*u*(psi.dx(1))
-(omeg*r*sin(theta)*u))
a = (a1+a2)
f = Constant(0.0)
g = Constant(0.0)
a = (a1+a2)*dx
L = f*v*dx+g*u*dx
## Compute solution
w = Function(W)
solve(a == L, w, bcs)
# split using deep copy
(psi, omeg) = w.split(True)
w_eval =-1.0/(r*sin(theta))*(psi.dx(0).dx(0)+psi.dx(1).dx(1))+((x[0]*x[0])/(r*x[1]*x[1]))*psi.dx(1)-(x[0]/(x[1]*r*r))*psi.dx(1)
w_chk= ((x[0]*w_eval.dx(0)) + (x[1]*w_eval.dx(1)) - w_eval)
M = assemble(w_chk*ds(1))
```