Let phi be a function given as a matrix N*N and representing the values of the function at the vertices of a unit square mesh. I want to interpolate this function on the finite element basis. The interpolated function is called psi. I tried this:
mesh = RectangleMesh(0, 0, 1, 1, N, N)
V = FunctionSpace(mesh, "Lagrange", 1)
psi = Function(V)
phi_vector = phi.reshape((N+1)*(N+1))
psi.vector()[:] = phi_vector
but the resulting function psi is not correct. I guess this is because the basis functions are not indexed like the mesh vertices. Any idea how to do the interpolation then? Is there a way to find the indexing of the basis functions with respect to the mesh vertices?