Hi dokken, sorry for the late reply and here is a minimal example.

```
import dolfinx.io.gmshio as gmshio
from mpi4py import MPI
from dolfinx import mesh, fem, io
import pygmsh
import numpy as np
# Create 2d mesh of a unit square
domain_2d = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
# Use pygmsh to generate 1/8 of a cylinder with tetrahedron mesh
with pygmsh.geo.Geometry() as geom:
p1 = geom.add_point([0.0, 0.0, 0.0], 0.1)
p2 = geom.add_point([1.0, 0.0, 0.0], 0.1)
p3 = geom.add_point([1.0, 1.0, 0.0], 0.1)
p4 = geom.add_point([0.0, 1.0, 0.0], 0.1)
l1 = geom.add_line(p1, p2)
l2 = geom.add_line(p2, p3)
l3 = geom.add_line(p3, p4)
l4 = geom.add_line(p4, p1)
cl = geom.add_curve_loop([l1, l2, l3, l4])
plane = geom.add_plane_surface(cl)
left, cylinder, lateral = geom.revolve(plane, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0], -0.5 * np.pi)
geom.generate_mesh()
pygmsh.write("cylinder.msh")
# Load mesh and meshtags
domain_3d, cell_markers, facet_markers = gmshio.read_from_msh("cylinder.msh", MPI.COMM_WORLD, gdim=3)
# Define function spaces
V_2d = fem.functionspace(domain_2d, ("CG", 1))
V_3d = fem.functionspace(domain_3d, ("CG", 1))
# Define a function in 2d
u_2d = fem.Function(V_2d)
u_expr = lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2
u_2d.interpolate(u_expr)
# Define the interpolation data between non-matching meshes
interpolation_data = fem.create_nonmatching_meshes_interpolation_data(domain_3d._cpp_object, V_3d.element, domain_2d._cpp_object)
# Create the function in 3d function space and interpolate
u_3d = fem.Function(V_3d)
u_3d.interpolate(u_2d, nmm_interpolation_data = interpolation_data)
# Write to .xdmf file and view in paraview:
if not os.path.exists("mve"):
os.makedirs("mve")
with io.XDMFFile(MPI.COMM_WORLD, "mve/u_2d.xdmf", "w") as file:
file.write_mesh(domain_2d)
file.write_function(u_2d)
with io.XDMFFile(MPI.COMM_WORLD, "mve/u_3d.xdmf", "w") as file:
file.write_mesh(domain_3d)
file.write_function(u_3d)
```

The result is as follows

As one can see that this interpolation works perfectly for the plane

z=0 (front face). What I would like to do is to interpolate this to all cross-sections of this cylinder, i.e. the function

`u_3d`

is axisymmetric with respect to

y-axis.