How to get the first order derivative of the solution and store it in a numpy array

I would appreciate if @dokken would comment. This seems to be the usual way I had done this (see for example Nedelec elements or Lagrange? What about the order? - #2 by raboynics) however this lead me to strange results (result seems to converge with mesh refinement but not with higher polynomial degrees).

This way of computing derivatives seems easier to implement than the one you just linked to, dokken (Numerical values from ufl.SpatialCoordinate - #2 by dokken), where it is question to define a quadrature scheme with basix.