Interpolation is done on a cell level. This means that to interpolate a function into a function space, one computes it cell by cell.

As you have 4 cells, with 4 dofs in each cell, this means that you have 16 interpolation points to consider.

The reason for this is that not every basis function is defined with a functional that is a point evaluation, there are basis functions defined by integral moments, which needs co-variant or contra-variant Piola mappings, which are defined on a cell level.

Locate the degree of freedom based on its dof coordinate (use V.tabulate_dof_coordinates() to find the relevant coordinate index to access in u.x.array