# Different results of area of a deformed mesh in parallel

**URL:** <https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105>\
**Category:** dolfinx\
**Tags:** mesh\
**Created:** [December 20, 2023, 7:10am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105 "2023-12-20T07:10:56Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 20, 2023, 7:10am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/1 "2023-12-20T07:10:56Z")

</div>

I deformed a mesh by appending a deformation array to the mesh with `mesh.geometry.x[:, d] += deformation_array.x.array[:]`, where `d` is the dimension of mesh to be moved. I got the correct result with a serial run. **The result was wrong with parallel runs but the visualization was correct.**

To brief the problem: I have two non-matching meshes in a FSI problem. Every time I need to interpolate the fields from the background (bigger) domain to the foreground (smaller) domain and then form the deformation array. In the following code, serial run returns the expected result 1, but it gives 1.30419921875 with  
`mpirun -n 4 python parallel_area.py`  
The doflinx version is 0.6.0.

MWE:

```auto
# parallel_area.py
# Test to measure area of mesh after deformation in parallel 

from mpi4py import MPI
from dolfinx import fem, io, mesh
from ufl import dx

# Dimension of the problem
ndims = 2

# Foreground mesh and function space
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
mesh_square.geometry.x[:, :2] -= 0.5 # move center to (0, 0) with size (-0.5, -0.5) -> (0.5, 0.5)
V_square = fem.VectorFunctionSpace(mesh_square, ("Lagrange", 1))
Q_square = fem.FunctionSpace(mesh_square, ("Lagrange", 1))

# Background mesh and function space
mesh_background = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_background.geometry.x[:, :2] -= 0.5 # move center to (0, 0)
mesh_background.geometry.x[:, :2] *= 2 # double the size: (-1, -1) -> (1, 1)
V_background = fem.VectorFunctionSpace(mesh_background, ("Lagrange", 1))
Q_background = fem.FunctionSpace(mesh_background, ("Lagrange", 1))

# Original area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area0 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(one*dx)), op = MPI.SUM)
print("Original area:", area0) # always returns 1 as expected
with io.XDMFFile(MPI.COMM_WORLD, "original_area.xdmf", "w") as file:
    file.write_mesh(mesh_square)

## Deformation occurs
# Define background fields
u_background = fem.Function(V_background)
u_background.sub(0).interpolate(lambda x: x[1]) # like a shear flow

# Interpolate to foreground fields
u_square = fem.Function(V_square)
for d in range(ndims):
    u_square.sub(d).interpolate(u_background.sub(d))

# Deform the foreground mesh
dt = 1 # an coefficient which mimics the one in the real problem
for d in range(ndims):
    mesh_square.geometry.x[:, d] += u_square.sub(d).collapse().x.array[:] * dt

# Compute new area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area1 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(one*dx)), op = MPI.SUM)
with io.XDMFFile(MPI.COMM_WORLD, "changed_area.xdmf", "w") as file:
    file.write_mesh(mesh_square)

print("Changed area:", area1) # mpirun -n 1: 1; -n 2: 1; -n 3: 1.45068359375; -n 4: 1.30419921875

```

 ![area](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/2cb3719ca617b56f52216f2feea6683ae4bc5742.jpeg)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 20, 2023, 10:57am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/2 "2023-12-20T10:57:23Z")

</div>

Running:

```python
with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)

```

you can see that the mapping actually isnt correct in parallel.

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/41679a4e093cfe4c7cd5194b981b1a7b81214426.png)

To fix this consider:

```python
# parallel_area.py
# Test to measure area of mesh after deformation in parallel

import numpy as np
from mpi4py import MPI
from dolfinx import fem, io, mesh
from ufl import dx

# Dimension of the problem
ndims = 2

# Foreground mesh and function space
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
# move center to (0, 0) with size (-0.5, -0.5) -> (0.5, 0.5)
mesh_square.geometry.x[:, :2] -= 0.5
V_square = fem.VectorFunctionSpace(mesh_square, ("Lagrange", 1))
Q_square = fem.FunctionSpace(mesh_square, ("Lagrange", 1))

# Background mesh and function space
mesh_background = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_background.geometry.x[:, :2] -= 0.5 # move center to (0, 0)
# double the size: (-1, -1) -> (1, 1)
mesh_background.geometry.x[:, :2] *= 2
V_background = fem.VectorFunctionSpace(mesh_background, ("Lagrange", 1))
Q_background = fem.FunctionSpace(mesh_background, ("Lagrange", 1))

# Original area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area0 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(one*dx)), op=MPI.SUM)
# always returns 1 as expected
print("Original area:", area0)
with io.XDMFFile(MPI.COMM_WORLD, "original_area.xdmf", "w") as file:
    file.write_mesh(mesh_square)

# Deformation occurs
# Define background fields
u_background = fem.Function(V_background)
u_background.sub(0).interpolate(lambda x: x[1]) # like a shear flow

# Interpolate to foreground fields
u_square = fem.Function(V_square)
for d in range(ndims):
    nms = fem.create_nonmatching_meshes_interpolation_data(
        u_square.sub(d).function_space.mesh._cpp_object,
        u_square.sub(d).function_space.element,
        u_background.sub(d).function_space.mesh._cpp_object, padding=1e-6)
    u_square.sub(d).interpolate(
        u_background.sub(d), nmm_interpolation_data=nms)
u_square.x.scatter_forward()
# Deform the foreground mesh
# an coefficient which mimics the one in the real problem
dt = 1
for d in range(ndims):
    V_sub, sub_to_parent = u_square.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_square.sub(
        d).collapse().x.array * dt
# Compute new area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area1 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(one*dx)), op=MPI.SUM)
with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)
print("Changed area:", area1)

```

which I’ve executed with DOLFINx v0.7.2  
yielding:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d11277de13e051421ca19303b1ab0188092cc7c4.png)  
for the 0th process.

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 20, 2023, 3:36pm UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/3 "2023-12-20T15:36:15Z")

</div>

I see. I thought the `>1` result is because I added the contribution of the ghost cells for serveral times in parallel. The real problem in my code is the wrong mapping between mesh and function. So I replace with the following part and I can get 1 now with DOLFINx 0.6.0:

> [@dokken](#):
>
> ```auto
> for d in range(ndims):
> V_sub, sub_to_parent = u_square.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_square.sub(
> d).collapse().x.array * dt
> 
> ```

Some follow-up questions:  
(1) `VTXWriter` vs `XDMFFile`. Why `XDMFFile` still gives correct visualization even though the mapping is actually not correct in parallel but `VTXWriter` can dectect that?  
(2) The result of area 1 with `mpirun -n 1` or `mpirun -n 2` in my original code is evaluated coincidentally？The wrong mapping should be independent of number of processes.  
(3) I see that you use `sub_to_parent = sub_to_parent // ndims` is because I used a (Vector) Function whose DOFs are ordered. What if I use a MixedElement and FunctionSpace instead of VectorFunctionSpace here? [The dofs are most usually ordered for a mixed element.](https://fenicsproject.discourse.group/t/dofmap-of-fem-vectorfunctionspace-vs-fem-functionspace-mesh-mixedelement/9411/4)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 20, 2023, 4:17pm UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/4 "2023-12-20T16:17:28Z")

</div>

> [@chenyongxin](#):
>
> `VTXWriter` vs `XDMFFile`. Why `XDMFFile` still gives correct visualization even though the mapping is actually not correct in parallel but `VTXWriter` can dectect that?

There are subtle differences in the implementation of the two. Also note that your XDMFFile is using

> [@chenyongxin](#):
>
> `with io.XDMFFile(MPI.COMM_WORLD, "changed_area.xdmf", "w") as file:`

while VTXWriter is in this case using

> [@dokken](#):
>
> `with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:`

writing data per process, which is possible due to how VTXWriter is created.

> [@chenyongxin](#):
>
> (2) The result of area 1 with `mpirun -n 1` or `mpirun -n 2` in my original code is evaluated coincidentally？The wrong mapping should be independent of number of processes.

The issue is with Ghost degrees of freedom. These depend on the number of processes.  
Thus the mapping depend on the number of processes.

> [@chenyongxin](#):
>
> 1. I see that you use `sub_to_parent = sub_to_parent // ndims` is because I used a (Vector) Function whose DOFs are ordered. What if I use a MixedElement and FunctionSpace instead of VectorFunctionSpace here?

Then things become more complicated. Please make a minimal reproducible example

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 21, 2023, 6:35am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/5 "2023-12-21T06:35:35Z")

</div>

> [@dokken](#):
>
> > [@chenyongxin](#):
> >
> > 1. I see that you use `sub_to_parent = sub_to_parent // ndims` is because I used a (Vector) Function whose DOFs are ordered. What if I use a MixedElement and FunctionSpace instead of VectorFunctionSpace here?
> 
> Then things become more complicated. Please make a minimal reproducible example

In the real applications, I use MixedElement instead of (Vector)Function. As the DOFs of a MixedElement are not always ordered, using the same way with a MixedElement gives wrong result as follows:

```auto
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
FE = ufl.FiniteElement("Lagrange", mesh_square.ufl_cell(), 1)
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

# Deform the mesh with a deformation array
dt = 1 # mimics the coefficient in the real applications
for d in range(ndims):

    V_sub, sub_to_parent = u_square.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_square.sub(
        d).collapse().x.array * dt

# 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)       
# mpirun -n 1: 0.99951171875
# mpirun -n 2: 1.02880859375

```

My understanding: for `mesh`, the vertices are ordered with vertex\_id: `mesh.geometry.x[:, 0]` → `[vertex_0_x, vertex_1_x, vertex_2_x, ...]`. For `Function`, `function.x.array[:]` → `[dof_0, dof_1, dof_2, ...]`. For a general usage, I guess here I need an equivalent `vertex_to_dof` functionality instead of `sub_to_parent = sub_to_parent // ndims`.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 21, 2023, 6:45am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/6 "2023-12-21T06:45:52Z")

</div>

> [@chenyongxin](#):
>
> In the real applications, I use MixedElement instead of (Vector)Function. As the DOFs of a MixedElement are not always ordered, using the same way with a MixedElement gives wrong result as follows:

Is there a specific motivation for using a Mixed element? I get that if you use it, you end up with this issue, but why not use a vector-element?

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 21, 2023, 6:57am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/7 "2023-12-21T06:57:49Z")

</div>

It is a fully-coupled FSI problem. All variables are solved simultaneously but not solved in a multi-stage decouple model, so I need to use MixedElement to include all variables (ux-uy-uz-p).

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 21, 2023, 7:00am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/8 "2023-12-21T07:00:06Z")

</div>

You could use a Mixed element of a vector element + pressure element. I think that would simplify the code alot.  
This is similar to what is done in the stokes demo in dolfinx

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 21, 2023, 7:20am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/9 "2023-12-21T07:20:01Z")

</div>

Irregardless, ill see if I have time to code up a dof to vertex map for you tonight.

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 28, 2023, 7:44am UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/10 "2023-12-28T07:44:27Z")

</div>

> [@Application of point forces, mapping vertex indices to corresponding DOFs](https://fenicsproject.discourse.group/t/application-of-point-forces-mapping-vertex-indices-to-corresponding-dofs/9646/2):
>
> Consider the following code, that can map from any mesh node to any Lagrange space: import ufl import dolfinx from mpi4py import MPI import numpy as np import scipy def f(x): return x[0] + 3 \* x[1] N = 4 P = 1 Q = 2 mesh = dolfinx.mesh.create\_unit\_cube(MPI.COMM\_WORLD, N, N, N) el0 = ufl.VectorElement("Lagrange", mesh.ufl\_cell(), P) el1 = ufl.FiniteElement("Lagrange", mesh.ufl\_cell(), Q) el = ufl.MixedElement([el0, el1]) V = dolfinx.fem.FunctionSpace(mesh, el) u = dolfinx.fem.Function(V…

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?

```auto
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)       

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 29, 2023, 6:52pm UTC](https://fenicsproject.discourse.group/t/different-results-of-area-of-a-deformed-mesh-in-parallel/13105/11 "2023-12-29T18:52:46Z")

</div>

> [@chenyongxin](#):
>
> 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?

As long as you work with function spaces that has a point evaluation dof at the vertices, you can skip the interpolation and do the following:

```python
# parallel_area.py
# Test to measure area of mesh after deformation in parallel

from IPython import embed
import numpy as np
from mpi4py import MPI
from dolfinx import fem, io, mesh, cpp
import ufl

# Dimension of the problem
ndims = 2

# Foreground mesh and function space
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
# move center to (0, 0) with size (-0.5, -0.5) -> (0.5, 0.5)
mesh_square.geometry.x[:, :2] -= 0.5
el = ufl.FiniteElement("Lagrange", mesh_square.ufl_cell(), 2)
mixed_el = ufl.MixedElement([el, el])
V_square = fem.FunctionSpace(mesh_square, mixed_el)

# Background mesh and function space
mesh_background = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_background.geometry.x[:, :2] -= 0.5 # move center to (0, 0)
# double the size: (-1, -1) -> (1, 1)
mesh_background.geometry.x[:, :2] *= 2
V_background = fem.FunctionSpace(mesh_background, mixed_el)

# Original area
f = 1*ufl.dx(domain=mesh_square)
area0 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(f)), op=MPI.SUM)

# always returns 1 as expected
print("Original area:", area0)

# Deformation occurs
# Define background fields
u_background = fem.Function(V_background)
u_background.sub(0).interpolate(lambda x: x[1]) # like a shear flow

# Interpolate to foreground fields
u_square = fem.Function(V_square)
for d in range(ndims):
    nms = fem.create_nonmatching_meshes_interpolation_data(
        u_square.sub(d).function_space.mesh._cpp_object,
        u_square.sub(d).function_space.element,
        u_background.sub(d).function_space.mesh._cpp_object, padding=1e-6)
    u_square.sub(d).interpolate(
        u_background.sub(d), nmm_interpolation_data=nms)
u_square.x.scatter_forward()

with io.VTXWriter(u_square.function_space.mesh.comm, f"u_square.bp", [u_square.sub(0).collapse()], engine="BP4") as file:
    file.write(0.0)
with io.VTXWriter(u_background.function_space.mesh.comm, f"u_background.bp", [u_background.sub(0).collapse()], engine="BP4") as file:
    file.write(0.0)

with io.VTXWriter(MPI.COMM_SELF, f"orginal_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)

# Create a vertex to dof map for the mixed element
num_vertices = mesh_square.topology.index_map(
    0).size_local + mesh_square.topology.index_map(0).num_ghosts
tdim = mesh_square.topology.dim
num_cells = mesh_square.topology.index_map(
    tdim).size_local + mesh_square.topology.index_map(tdim).num_ghosts
vertex_to_dofmaps = np.zeros(
    (V_square.num_sub_spaces, num_vertices), dtype=np.int32)
num_vertices_per_cell = cpp.mesh.cell_num_vertices(
    mesh_square.topology.cell_types[0])
c_to_v = mesh_square.topology.connectivity(tdim, 0)
vertices_to_mesh_nodes = cpp.mesh.entities_to_geometry(
    mesh_square._cpp_object, 0, np.arange(num_vertices, dtype=np.int32), False).reshape(-1)
for d in range(vertex_to_dofmaps.shape[0]):
    dofmap = V_square.sub(d).dofmap
    local_vertex_map = np.zeros(num_vertices_per_cell, dtype=np.int32)
    # For given sub space, find the local vertex to dof map for a cell
    for v in range(num_vertices_per_cell):
        v_to_d = dofmap.dof_layout.entity_dofs(0, v)
        assert len(v_to_d) == 1
        local_vertex_map[v] = v_to_d[0]

    for cell in range(num_cells):
        vertices = c_to_v.links(cell)
        for v in range(num_vertices_per_cell):
            vertex_to_dofmaps[d, vertices[v]] = dofmap.cell_dofs(
                cell)[local_vertex_map[v]]

# Deform the foreground mesh
# an coefficient which mimics the one in the real problem

dt = 1
for d in range(ndims):
    mesh_square.geometry.x[vertices_to_mesh_nodes,
                           d] += u_square.x.array[vertex_to_dofmaps[d]]

# Compute new area
f = 1*ufl.dx(domain=mesh_square)
area1 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(f)), op=MPI.SUM)
with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)
print("Changed area:", area1)

```

Some neat things in this code snippet:

1. For a time dependent problem, you can create `vertices_to_mesh_nodes` and `vertex_to_dofmaps` once, and use them in a loop.
2. No interpolation or collapsing is required inside the mesh update loop,i.e.

```python
for d in range(ndims):
    mesh_square.geometry.x[vertices_to_mesh_nodes,
                           d] += u_square.x.array[vertex_to_dofmaps[d]]

```

as we used the un-collapsed spaces to make the vertex to dof-map.