Putting unknowns in user defined expressions

I am solving a registration problem where the displacemente field appears in an image argument: T( x + u(x) ), which makes a nonlinear term appear in the formulation. I am solving it with a Newton method (at the algebraic level, as the documentation once said) and the forms definition are VERY slow, and also don't allow for a parallel implementation. I attach the relevant piece of code, and would greatly thank any kind of help.

First load the images into the problem instance:

    image_R = self.R
    image_T = self.T
    image_dT= self.dT
    image_HT = self.HT
    element = V.ufl_element()

Here define the new expressions which incorporate the displacement

    class scalarExpression(Expression):
        def __init__(self, im, element, u = Expression(("0.","0."), degree = 0)):
            self.im = im
            self.element = element
            self.u = u
        def eval(self, value,x):
            value[0] = self.im(x + self.u(x))
    class vectorExpression(Expression):
        def __init__(self, im, element, u = Expression(("0.","0."), degree = 0)):
            self.im = im
            self.element = element
            self.u = u
        def eval(self, value, x):
            value[0] = self.im[0](x+self.u(x))
            value[1] = self.im[1](x+self.u(x))
    class tensorExpression(Expression):
        def __init__(self, im, element, u = Expression(("0.","0."), degree = 0)):
            self.im = im
            self.element = element
            self.u = u
        def eval(self, value, x):
            value[0] = self.im[0][0](x+self.u(x))
            value[1] = self.im[0][1](x+self.u(x))
            value[2] = self.im[1][0](x+self.u(x))
            value[3] = self.im[1][1](x+self.u(x))

Now define functionals which will give the expressions in every iteration

    gen_T = lambda u: scalarExpression(im = image_T, 
                                       element = FiniteElement("CG", mesh.ufl_cell(), N), u = u)
    gen_dT = lambda u: vectorExpression(im = image_dT,
                                        element = VectorElement("CG", mesh.ufl_cell(), N), u = u)
    gen_HT = lambda u: tensorExpression(im = image_HT,
                                        element = TensorElement("CG", mesh.ufl_cell(), N), u = u)

The image R won't need to be updated

    R = scalarExpression(im = image_R, element = FiniteElement("CG", mesh.ufl_cell(), N))
    T = gen_T(u0)
    dT = gen_dT(u0)
    HT = gen_HT(u0)

Solve a linearized version for initial guess

    a = (inner( self.C(eps(du_n)), eps(v) ) + 
            self.alpha*( (R-T)*dot(du_n, dot(HT,v)) - dot(dT,v)*dot(dT,du_n) ) )*dx
    f = dot(Constant((0.,0.)),v)*dx
    solve(a == f, u_n, solver_parameters = self.solver_params)
    it, error = 0, 1.
    du = u_n.copy()
    print "Begin iterations"
    while error > tol and it < maxiter: 
        print it, error
        T, dT, HT = gen_T(u_n), gen_dT(u_n), gen_HT(u_n)
        a = ( inner(C(eps(du_n)),eps(v)) + 
                self.alpha*(-dot(dT, v)*dot(dT, du_n) + (R - T)*dot(du_n, dot(HT, v))) )*dx
        f = -( inner(C(eps(u_n)), eps(v)) - self.alpha*(R - T)*dot(dT, v) )*dx
        solve( a == f, du, solver_parameters = self.solver_params)
        it += 1
        error = norm(du, mesh = mesh)
        u_n.assign(du + u_n)
    print "End"

Note: the images I defined are python callables or lists of callables depending on the dimensions. Also, if a natural way of incorporation function composition were possible, I would just use the nonlinear solver and save all the trouble. This is only possible with functions such as log, exp, etc (it is shown in the neo-hookean example) but not with expressions.