Error: Nonlinear solver for elasticity problem

Hi,
I am solving elasticity problem directly from energy form. I have got the solution for method 1. but when i have split energy into deviatoric and dilation part I am getting error: ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator to expression depending on form argument v_1.

any idea to fix this issue.

Thanks in advance

# test elasticity 
# purpose:
# to solve directly from energy form
#=======================================
from dolfin import *
mesh = UnitSquareMesh(2,2)
E, nu = Constant(10.0), Constant(0.3)
plot(mesh, title = 'mesh', interactive = True)
ndim = mesh.geometry().dim()
#------------------------
# Define sigma & strain
#------------------------
def eps(v):
     return sym(grad(v))
def eps_dev(du):
    return eps(du) - tr(eps(du))*Identity(ndim)

zero_v=Constant((0.,)*ndim)   
f=zero_v

V = VectorFunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
du = Function(V)

mu    = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu )*(1.0-2.0*nu))
#EDef = 0.5*inner(sigma(du),eps(du))*dx- dot(f,du)*dx
#------------------------------------------
# energy (method: 1)
#------------------------------------------
#EDef = ( mu*tr(eps(du)*eps(du)) + 0.5*lmbda*tr(eps(du))**2 )*dx - dot(f,du)*dx
#------------------------------------------------
# split energy into deviatoric and dialation part (method: 2)
#------------------------------------------------
EDef = ( 0.5*(lmbda + mu)*(  0.5*( tr(eps(du)) + abs(tr(eps(du))))   )**2 +      mu*(tr(eps_dev(du)*eps_dev(du))) + 0.5*(lmbda + mu)+(  0.5*( tr(eps(du)) - abs(tr(eps(du))))   )**2 )*dx - dot(f,du)*dx
EDer = derivative(EDef,du,v)
EDereP = replace(EDer,{du:u})
ED = derivative(EDereP,du,u)
class top(SubDomain):
     def inside(self,x,on_boundary):
             tol = 1e-10
             return abs(x[1]-1.0) < tol and on_boundary

class bottom(SubDomain):
     def inside(self,x,on_boundary):
           tol = 1e-10
           return abs(x[1]) < tol and on_boundary
Top = top()
Bottom = bottom()
bclx= DirichletBC(V.sub(0), Constant(0.0), Bottom)
bcly = DirichletBC(V.sub(1), Constant(0.0), Bottom)
bcty = DirichletBC(V.sub(1), Constant(1.0), Top)
bc_u = [bclx, bcly, bcty]
 # solver for method: 1
#problem_u = LinearVariationalProblem(lhs(EDereP), rhs(EDereP), du, bc_u)
#solver_u = LinearVariationalSolver(problem_u)
#solver_u.solve()
# solver for method: 2
solve(EDereP==0,du, bcs=bc_u)

point = (1.0,1.0)
print 'numerical u at the corner point:', du(point)
plot(du, interactive  = True)
u1, u2 = split(du)
plot(u1, interactive = True)
plot(u2, interactive = True)