Visualising shape functions

In your code you are interpolating an arbitrary function into a CG 1 function space. In most cases you will then loose a lot of information about the basis function.

Note that using pvd results in a CG-1 interpolation of your field. For arbitrary order CG/DG see: Loading xdmf data back in - #4 by dokken

Also, you can find plots of shape functions and explicit formulas for shape functions at: https://defelement.com/