Since it is possible to interpolate one function defined on a function space into another one defined on another function space (from a non matching mesh), like so:

I wonder if I can use the interpolate function to build an interpolation matrix that links the two function vectors u_1.vector and u_2.vector such that

\bf{u_2} = \bf{I} \bf{u_1}

where \bf{I} Is a rectangular matrix, such that \bf{u_2} is the same that I would obtain from the interpolation function

Oh, sorry, i didn’t see that your are doing interpolation across different meshes.

There is currently not a way or creating that matrix without doing it by hand, but you can cache alot of the data/communication required for the cross mesh interpolation, see:

create_nonmatching_meshes_interpolation_data() seems to associate the point where the function has to be interpolated with the cell that contains it. It is indeed not varying when the two meshes don’t change and saves a lot of computations for multiple interpolations.

Is there a way that allows me, given a cell and a point contained in it, to retrieve the shape functions and the corresponding cell nodes?

You can check out how dolfinx.fem.Function.eval works, as it does all the steps needed to pull back points from physical space to the reference, and use the underlying dolfinx.fem.FiniteElement to tabulate the basis functions as these points.

The working slice of code for multiple points is the following:

for i in range(0, len(cells)):
geom_dofs = msh_a.geometry.dofmap.links(cells[i])
x_ref[i,:] = msh_a.geometry.cmap.pull_back([points[i,:]], x[geom_dofs])

Sorry for the lack of context, but after I clean the complete code for the interpolation matrix, I am going to put it here.

If we wanted to interpolate the function, we could use the command:

p_1.interpolate(p_0)

But now we have an alternative. We can indeed compute the interpolation matrix and compute the p_1 array from a matrix multiplication between I_matrix and the p_0 array: