Why do non-interacting cylinders bend in unintuitive direction under identical forces?

My full sigma is this:

def sigma(u):
    return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)

Is this correct and if not what is the mistake?

When changing the mesh by setting these numbers I get better results. Is this the way to go or should I try something else?

gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1) 
gmsh.option.setNumber("Mesh.MeshSizeMin", 0.0001)  


Thanks for helping me with this!