I have a specific problem where I have to create a data set from two solution vectors(both live on the same mesh), one vector is a true solution which can be broadcasted as let us say ist column of the data set and from solution2 I have to create features of the dataset as follows:
for each nodal value of true solution, I have to extract surrounding 9 nodal points from solution2 which will be features corresponding to that nodal value.
ofcourse, there will be issue with the boundary points because it will not have enough surrounding points. We can either do it only for interior nodal values of true solution
or we can attach zero values for missing points(incase we do it for boundary points as well).
I am unable to wrap my head around this problem, so far what I have reached is posted below:

# define mesh and function space
N = 4
degree = 1
msh = create_unit_square(MPI.COMM_WORLD,N,N)
V = functionspace(msh, ("Lagrange", degree))
x = SpatialCoordinate(msh)
u_true = fem.Function(V)
u_true.interpolate(lambda x : x[0] + x[1])
arr_true = u_true.x.array
print(arr_true)
u_true.function_space.tabulate_dof_coordinates()
u2 = fem.Function(V)
u2.interpolate(lambda x : 2*x[0] + 2*x[1])
print(u2.x.array)
X_features = np.zeros((len(array),9))

I am also attaching picture for more clarity. I would appreciate your help!
True nodal value is in red and surrounding values in yellow, I have to do it for each interior true nodal value.

I guess you are assuming a structured grid (as depicted), as it wouldn’t make sense to do this on an unstructured mesh that is not aligned with the coordinate axis.

However, as your grid is structured, I guess you know the location of the neighbouring points. If so, I would simply search for them in tabulate_dof_coordinates and extract the index of the matching degree of freedom (accessed in u.x.array).

dr @dokken I tried something on this but it will give all nodal values from solution2 but I want only neighbouring nodal values whose location I dont know. How can I find location of neighbouring values corresponding to each nodal value of u_true. Also I don’t want to use boundary nodal values from u_true I want only interior values.

X_features = np.zeros((len(array),9))
for i in range(dof.shape[0]):
if dof_true[i].all() == dof[i].all():
X_features[i,:] = u2.x.array

dr @dokken Apologies for the interruption, can you please outline steps to do this, I am still unable to proceed with this, I need to grab 9 points(including the center point) from solution2 corresponding to each internal nodal value of u_true. Do I need to specify the neighbouring points manually? I think that is not feasible for large no. of dofs. How can I locate neighbouring points for each internal nodal value? I would appreciate if you can provide MWE.

dr @dokken I am just working on the built in unit square. All I want is to create a data set with u_true(internal nodal values as first column) and for each internal nodal value I want 9 neighbouring points as features from solution2, essentially getting a matrix of size (no. of internal nodal values by no. features(9)) . I have to approximate(predict) a function actually from finite element method not finite difference / I appreciate all the work you are doing. It is really inspiring for we people and we will contribute to others in future, once we gain some experience in fenicsx.