Solving 3D elasticity

Question

I am solving 3D elasticity problem. left side fixed and uniform pressure on top. code works but result is not matching with abaqus
sample code:
class Left(SubDomain):
def inside(self,x,on_boundary):
return (near (x[0],0.0) and on_boundary)

left = Left()

class Top(SubDomain ):
def inside(self , x, on_boundary ):
return (near (x[1],1.0) and on_boundary)

top = Top()
bc_fixed = DirichletBC(V, Constant ((0,0,0)), left)
bcs = [bc_fixed]
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all (0)
left.mark(boundaries ,1)
top.mark(boundaries ,2)
ds = ds(domain=mesh , subdomain_data=boundaries)

load = -1

g = Constant ((0,load, 0 ))
def epsilon(v):
return 0.5 * (grad(v) + transpose(grad(v)))
def sigma(v):
return 2 * mu * epsilon(v) + lmbda * tr(epsilon(v)) * Identity (3)
a = inner(grad(v), sigma(u))dx
L = inner(v,g)ds(2)
u = Function(V)
solve(a==L,u,bcs)
(u1, u2,u3) = u.split(True)

v_nodal_values = v.vector().array()

extract solution

u_nodal_values = u.vector()

u_array = u_nodal_values.array()

coor = mesh.coordinates()

for i in range(len(u_array)):

print"u(%8g,%8g,%8g) =%g" % \

(coor[i][0], coor[i][1], coor[i][2],u_array[i])

print solution

plot(u)

Comments

You really need to clean this up and write properly before you expect an answer.

---

Yes, this will not even compile since you have not defined the constants, the mesh or even included done basic debugging, e.g.

L = inner(v,g)ds(2)

---

I have solved this problem, Thank you for your reply