Extract solution values on surface nodes

Hello,

I am trying to extract the displacement field and plot the result on specific nodes.
I successfully extracted the coordinates of the volume and the values of the vectors on each node:

xyz = V.tabulate_dof_coordinates()
x_coor = xyz[:,0]
y_coor = xyz[:,1]
z_coor = xyz[:,2]

x_value, y_value, z_value = [], [], []
for i in range(len(x_coor)):
    x_value.append(u.x.array[3*i:3*(i+1)][0])
    y_value.append(u.x.array[3*i:3*(i+1)][1])
    z_value.append(u.x.array[3*i:3*(i+1)][2])

However, is it possible to do the same on a specific area of the volume ? I wanted to extract the nodes coordinates corresponding to an application surface of a force, but I couldn’t succeed in getting the solution value on these nodes:

fdim = mesh.topology.dim-1
num_facets_owned_by_proc = mesh.topology.index_map(fdim).size_local
app_surface_bis = ft.find(6)
geometry_entitites = dolfinx.cpp.mesh.entities_to_geometry(mesh, fdim, np.arange(num_facets_owned_by_proc, dtype=np.int32), False)
surf_entitites = dolfinx.cpp.mesh.entities_to_geometry(mesh, fdim, app_surface_bis, False)

points = mesh.geometry.x
e_list = []

for e, entity in enumerate(geometry_entitites):
    # print(e, points[entity])
    surf_list = [points[entity] for e, entity in enumerate(surf_entitites)]
    for i in range(len(surf_list)):
        if np.any(np.all(points[entity] == surf_list[i])):
            e_list.append(e)

surf_x_coor, surf_y_coor, surf_z_coor = [], [], []
for x in surf_list:
    for i in range(3):
        surf_x_coor.append(x[i][0])
        surf_y_coor.append(x[i][1])
        surf_z_coor.append(x[i][2])

See for instance: Fenicsx: use an array of data for boundary conditions - #5 by dokken