Based on the above your post, I still didn’t get it how to deal with the DOFs of a mixed element. Will be waiting for your code of map method here 
I think I can work around the mixed element DOF problem with the following code, which is based on the facts:
(1) The DOFs of VectorFunctionSpace are ordered;
(2) The DOFs of a mixed element are not always ordered;
(3) The vertex-to-DOF mapping is only valid with P1 element.
So firstly I interpolate a higher order mixed element function (P2) to a vector function (P1) and then move mesh. Any comments and/or suggestions?
from mpi4py import MPI
from dolfinx import fem, mesh
import ufl
import numpy as np
# Dimension of the problem and mesh
ndims = 2
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
# Original area
area0 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(1*ufl.dx(mesh_square))),
op = MPI.SUM)
print("Original area:", area0) # always returns 1 as expected
## Deformation occurs: replace (Vector)Function(Space) to a MixedElement with P1-P1 element
deg = 2
FE = ufl.FiniteElement("Lagrange", mesh_square.ufl_cell(), deg)
ME = ufl.MixedElement([FE, FE])
V_ME = fem.FunctionSpace(mesh_square, ME)
u_square = fem.Function(V_ME)
u_square.sub(0).interpolate(lambda x: x[1]) # like a shear flow
# To interpolate to a Vector Function with P1 element as a deformation array
V = fem.VectorFunctionSpace(mesh_square, ("Lagrange", 1))
u = fem.Function(V)
for d in range(ndims):
u.sub(d).interpolate(u_square.sub(d))
# Deform the mesh with a deformation array
for d in range(ndims):
V_sub, sub_to_parent = u.sub(d).function_space.collapse()
sub_to_parent = np.asarray(sub_to_parent)
sub_to_parent = sub_to_parent // ndims
mesh_square.geometry.x[sub_to_parent, d] += u.sub(d).collapse().x.array
# Compute new area
area1 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(1*ufl.dx(mesh_square))),
op = MPI.SUM)
print("Changed area:", area1)