I cannot understand why such a simple problem is giving an error. My code is given below:

```
from fenics import *
from dolfin import *
import ufl
import dolfin.cpp as cpp
import numpy
# Scaled variables
L = 1; W = 1
mu = 7.6923e10
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 11.53846e10
lambda_ = beta
g = gamma
#Create mesh and define function space
mesh = RectangleMesh(Point(0,0), Point(100, 100), 100, 100, "left")
V = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
class bottom(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-14
return (near(x[1],0.00)) and on_boundary
class top(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-10
return near(x[1],100.00) and on_boundary
class Marked(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-10
return near(x[0]+x[1],100,0.001) and x[0]>25 and x[1]>25
boundaries = MeshFunction("size_t",mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
bottom= bottom()
top = top()
bottom.mark(boundaries,1)
top.mark(boundaries,2)
Marked=Marked()
Marked.mark(boundaries,3)
meshfile=File('mesh.pvd')
meshfile<<boundaries
bc_bottom = DirichletBC(V.sub(1),Constant((0)),bottom)
bc = [bc_bottom]
ex = Constant((1.0, 0.0))
ey = Constant((0.0, 1.0))
dummy1 = as_tensor([\
[1.0, 0.0],\
[0.0, 0.0]\
])
dummy2 = as_tensor([\
[0.0, 1.0],\
[0.0, 0.0]\
])
dummy3 = as_tensor([\
[0.0, 0.0],\
[1.0, 0.0]\
])
dummy1 = as_tensor([\
[0.0, 0.0],\
[0.0, 1.0]\
])
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)
# Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
def sigel1(u):
return inner(sigma(u),dummy1)
def sigel2(u):
return inner(sigma(u),dummy2)
def sigel3(u):
return inner(sigma(u),dummy3)
def sigel4(u):
return inner(sigma(u),dummy4)
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0,0))
T = Constant((0, 1000))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds(2)
# Compute solution
u = Function(V)
solve(a == L, u, bc)
Tx = assemble(0.7071*(sigel1(u)+sigel2(u))*dS(3))
print Tx
Ty = assemble(0.7071*(sigel3(u)+sigel4(u))*dS(3))
print Ty
```

Please, please look into it. Thank you!`Preformatted text`