I am new to NGSolve and I’d like to ask two basic questions.

How can I get the ordering of the basis functions on each element?
I know we can use fes.GetDofNrs(el) to get the DOF number for each element. But I wonder how can I know what basis function corresponds to a given DOF. For example, I am using L2 space for quadrilateral meshes. I know the L2 space uses a tensor product of Legendre Polynomials in x and y directions. But I notice that in some elements, the basis functions are ordering like 1, y, x, xy…, while in some other element, they are ordering like 1, x, y, xy. Is there an easy way to know the ordering of the basis functions in each element?

Can I do interpolation with a set of points?
I know we can use GridFunction.Set to interpolate a CoefficientFunction, but I wonder if I have a set of points and values (x_i, y_i, eval_i), can I interpolate (or do least squares) to an element in the finite element space?

the transformation from local coordinates to global coordinates depends on the element vertex numbers. This ensures consistent orientation on edges and faces.
It would not be necessary for L2, but we made L2 consistent with the other spaces.

The face f is ordered with the vertex numbers, f[0] becomes the smallest.

Here is a pure Python-solution for least squares approximation of arbitrary point data.
If you need performance, you can recode the same in C++.
If you want to interpolate data from a structured grid, you can use the VoxelCoefficientFunction.
If your interpolation points are given on the reference element, you can use user-defined integration-rules.