Error: Interpolate coordinates to high order GridFunction

I want to obtain the spatial coordinates (x, y) corresponding to the DOFs of a higher-order finite element function. The code is as follows:

from ngsolve import *
from ngsolve.webgui import Draw

import matplotlib.pyplot as plt

# create mesh
mesh = Mesh(unit_square.GenerateMesh(maxh=1))

init_refine =  1
for i in range(init_refine):
    mesh.Refine()

L = 2**init_refine
dh = 1
for p in mesh.ngmesh.Points():
    p[0] = p[0]*L - 0.5*L
    p[1] = p[1]*L - 0.5*L

mesh.ngmesh.Update() 

# Declare a FE space 
k = 2
Q_phi = H1(mesh, order=k)

# Declare grid function to represent coordinates x and y
dof_x = GridFunction(Q_phi)
dof_y = GridFunction(Q_phi)

# interpolate/Set x,y coordinates to dofs
# dof_x.Set(x)
# dof_y.Set(y)
dof_x.Interpolate(x)
dof_y.Interpolate(y)

print("ndof = ", Q_phi.ndof)
print("x = ", dof_x.vec.data)
print("y = ", dof_y.vec.data)

plt.scatter(dof_x.vec.FV().NumPy(), dof_y.vec.FV().NumPy(), alpha=0.2)
plt.axis('equal')
plt.grid()
plt.show()

Draw(dof_x, deformation=True)
Draw(dof_y, deformation=True)

Both Set() and Interpolate() functions do not work.

From the resulting plot, it can be seen that the spatial coordinates of the DOFs defined at the vertices are correct. However, the spatial coordinates of the DOFs defined at the edges are incorrect; they are all set to be near in (0,0). So, the color appears darker at the (0,0) location in the image.

The correct result should be:

However, the results displayed by the Draw() function in ngsolve are correct. Why is that?

Attachment:
high_order_coordinate_interpolate.ipynb (22.1 KB)

H1 space doesn’t have a nodal high order basis, but a hierarchic

best, Christopher

My question is how to obtain the spatial coordinates of DOFs, especially the spatial coordinates of DOFs with COUPLING_TYPE belonging to INTERFACE_DOF and LOCAL_DOF on higher-order elements.