Odd strains with linear elasticity

I'm trying to solve a linear elasticity problem with varying dilation in the domain. (Random in this case). In order for the solver to converge, I need to fix translation and rotation which I do so by fixing the displacement of one corner, and the vertical displacement of the other. It solves fine, but gives exaggerated (wrong) strains at the corners where I have the conditions (as visible in a plot of the strain energy density).
I haven't been able to fix it with small tolerances, different solvers, higher resolutions etc. Anyone have an explaination or a better way to solve this problem?
Thanks!

from dolfin import *
from random import uniform

comm = mpi_comm_world()

class InitialCondition(Expression):
    def eval(self, values, x):
        values[0] = .5+1*uniform(-.1,.1) 

mesh = RectangleMesh(-25,-25,25,25, 50, 50, 'crossed')

V = FunctionSpace(mesh, "Lagrange", 1)
V_vec = VectorFunctionSpace(mesh, "Lagrange", 1)

c = Function(V)
u = TrialFunction(V_vec)
test_u = TestFunction(V_vec)
u_old = Function(V_vec)

def eps(u):
            return as_vector([u[i].dx(i) for i in range(2)] +
                     [u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])


boundary = CompiledSubDomain("on_boundary")
left = CompiledSubDomain("(std::abs(x[0]-d) < DOLFIN_EPS) && on_boundary", d = -25)
right = CompiledSubDomain("(std::abs(x[0]-d) < DOLFIN_EPS) && on_boundary", d = 25)
bottom = CompiledSubDomain("(std::abs(x[1]-d) < DOLFIN_EPS) && on_boundary", d = -25)
top = CompiledSubDomain("(std::abs(x[1]-d) < DOLFIN_EPS) && on_boundary", d = 25)
LeftBottom = CompiledSubDomain("(std::abs(x[0]-w) < DOLFIN_EPS) && (std::abs(x[1]-h) < DOLFIN_EPS)", w = -25, h = -25)   
LeftTop = CompiledSubDomain("(std::abs(x[0]-w) < DOLFIN_EPS) && (std::abs(x[1]-h) < DOLFIN_EPS)", w = -25, h = 25)

bcs = [
     DirichletBC(V_vec, [0,0], LeftBottom, 'pointwise'),
     DirichletBC(V_vec.sub(1), 0, LeftTop, 'pointwise')]

au = inner(eps(u), eps(test_u))*dx
Lu = 1e-0*(c-.5)*inner(as_vector((.04,.04, 0)), eps(test_u))*dx
u = Function(V_vec)

c = interpolate(InitialCondition(), V)

A, b = assemble_system(au, Lu, bcs)

file_u = XDMFFile(comm, 'output/disp.xdmf')

file_u << (u, 0.)

solve(A, u.vector(), b)
file_u << (u, 1.)

plot(c)
plot(inner(eps(u),eps(u)), interactive = True)