Efficient way to move an internal domain

Hello,

I have two domains including a square and a rectangle which is inside the square. I want to move the internal rectangle after calculation of displacements of the nodes corresponding to this domain by solving a linear elasticity problem separately on this domain. Here are the domains:

I created the second mesh which is exactly the same as the red domain and solved a linear elasticity problem. Then I moved the coordinates of the red domain in the first mesh according to the displacements obtained from the second mesh. The problem is that the displacements are different along the length of the red domain. After moving the coordinated the elements go inside each other (Mostly at the bottom of the red domain). Here is what happens:

And here is the code reproducing the above results:

from dolfin import *
parameters['allow_extrapolation'] = True

mesh = UnitSquareMesh(30,30)

# Defining the second mesh which is consistent to the internal rectangle domain with same number of elements
mesh2 = RectangleMesh(Point(0.4,0.3), Point(0.6, 1),6, 21)

# Marking the top edge
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1],1)

top = Top()
boundaries = MeshFunction('size_t', mesh2, mesh.topology().dim()-1)
boundaries.set_all(0)
top.mark(boundaries, 1)

File("boundaries.pvd") << boundaries

# Solving the displacement filed for the internal rectangle
E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))

def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

V = VectorFunctionSpace(mesh2, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(Constant((1000, 0)), u_)*dx

bc = DirichletBC(V, Constant((0.,0.)), boundaries, 1)
u = Function(V)
solve(a == l, u, bc)

# Moving the nodes of the internal rectangle in the original mesh according to the obtained displacements
x = mesh.coordinates()

A = []
for i in range(len(x)):
    if x[i][0] >= 0.4 and x[i][0] <= 0.6 and x[i][1] >= 0.3 and x[i][1] <= 1:
        A.append(list(x[i]))
c = []
for i in range(len(A)):
    for j in range(len(x)):
        if A[i][0] == x[j][0] and A[i][1] == x[j][1]:
            c.append(j)
        else:
            pass
for i in c:
    x[i][0] = x[i][0] + u(x[i][0], x[i][1])[0]
    x[i][1] = x[i][1] + u(x[i][0], x[i][1])[1]

File("moved-mesh.pvd") << mesh

What I am looking at is to preserve the aspect ratio of the elements as much as possible. To be more specific I need the elements in blue domain also move to avoid this issue.
I was wondering if there is any efficient way so the quality of the elements is preserved after moving the internal domain.

See for instance: How to use ALE for moving boundary problems? - #2 by dokken