Result of Steady state dynamic linear elastic isnt correct

Hello,
I am new to FEnics, I want to solve the steady state dynamic linear elastic model in solid. my equation is function of frequency and the strong form is:
Divergence(stress(vec(x), w)) + w^2 * density * u( vec(x),w)=0
BC: Stress(vec(x),w) n(x)=T(x,w)
u(x,w)=U0
the physical problem is a plate with dimension of 1*1*0.1 with a harmonic load on the all the top and also the bottom of the plate is fixed.
I solved this as a test with the fixed frequency (w=200)

The weak form of the this equation is :
inner(sigma(u), sym(grad(v)))*dx -po*w*w *inner(v,u)*dx-inner(f,v)*ds(2) =0

I compared the results with the model with converged mesh from Abaqus, I almost checked everything in the code the order of the displacement is mostly correct but the signs of the displacements are wrong and also the error is large between the solution of the Abaqus and this code.

I checked the variation form so I cant find any error in it. So the only thing that I suspected to be wrong is the neumann boundary condition but I cant find any bug in it.

I would be very thankful if someone can help me to figure out what is wrong in my code that I get different answers from the code.

Here is the code::

from dolfin import *
import pickle
import numpy
import csv
po=2700
w=200
Magnitude=-100

mesh = BoxMesh(0.0,0.0,0.0 , 1.0,1.0,0.1 , 10,10,10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Defining Domain
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)
# Initialize sub-domain instances
top = Top()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
top.mark(boundaries, 2)
bottom.mark(boundaries, 4)
bc = DirichletBC(V, (0.0, 0.0, 0.0), boundaries,4)
# Define new measures associated with the interior domains and exterior boundaries
ds = Measure("ds")[boundaries]
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
f = Expression(("0.0", "scale","0.0"), w0 = w, scale = Magnitude )
# Elasticity parameters
E, nu = 69000000000, 0.3
mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

def sigma(u):
    return 2.0*mu*sym(grad(u)) + lmbda*tr(sym(grad(u)))*Identity(v.cell().d)

a = inner(sigma(u), sym(grad(v)))*dx -po*w*w *inner(v,u)*dx
L = inner(f,v)*ds(2)
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Save solution in VTK format
file = File("ElasticSolution.pvd")
file << u

plot(u, interactive=True)

I already wrote the weak form but since I couldn't find any example or demo similar to what I want to do . does any one can help me how should i deal with this problem and solve the problem in a range of frequency or give me some tips about this problem.

Thanks