Access GridFunction with Coordinates

For the reference, I’ve posted a code example to evaluate the solution function on a 2D mesh grid, using flattened numpy arrays: How to display equipotential lines for a 2D problem - #4 by pierre-haessig