Evaluate solution at grid points of structured mesh and reshape to tensor

Hello everyone,

I have a structured 3D grid and solve linear elasticity on that grid. After computing the solution vector uh, I would like to reshape it into an array of shape (Nx+1, Ny+1, Nz+1, 3) where Nx is the number of cells in x direction and so one, 3 is the number of component functions of the solution. The same should be done for a vector valued coefficient function C.

I am using the following domain and function spaces:

Nx = 20
Ny = 10
Nz = 10
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 1, 1])],
                         [Nx, Ny, Nz], cell_type=mesh.CellType.tetrahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
# solution space 
u = ufl.TrialFunction(V)

# space for coefficient function in PDE
space_Coeff = fem.functionspace(domain, ("CG", 1, (6,6)))
C = fem.Function(space_Coeff)

# after computing the solution vector uh from the variational problem we reshape
uh = uh.x.array.reshape(Nx+1, Ny+1, Nz+1, 3)
C_ = C.x.array.reshape(Nx+1, Ny+1, Nz+1, 6*6)

My question:
Is this the correct way of reshaping such that the tensors uh and C_ contain the values of the solution/coefficient function u and C at the grid points, i.e. uh[i_x, i_y, i_z, k] = u(i_x/Nx, i_y/Ny, i_z/Nz, k) where k=0,1,2 indices the component. Or does fenicsx use a different numbering when computing the solution vector uh over a structured grid such that the standard reshape does not fulfill this?

Thanks in advance and kind regards:)

See for instance:

DOLFINx in general does not consider a mesh structured. All meshes are unstructured, and data is reordered with respect to data-locality.