Assigning values to function: order issues

General
4

Consider something like the following code:

from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import fem, mesh, io

L = 2
W= 3
nx  = 3
ny = 5
nz = 7
msh = mesh.create_box(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0, 0.), (L, L, W)), n=(nx-1, ny-1, nz-1),
                            cell_type=mesh.CellType.hexahedron,)
dx, dy, dz = L/(nx-1), L/(ny-1), W/(nz-1)
values = np.arange(nx*ny*nz).reshape(nx, ny, nz)

P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = fem.FunctionSpace(msh, P1)
dof_coordinates = V.tabulate_dof_coordinates()
# dof index map to i,j,k
map_vec = np.zeros((len(dof_coordinates), 3), dtype=np.int32)
eps = 2e-14
for dof_index, coordinate in enumerate(dof_coordinates):
    map_vec[dof_index] = [(coordinate[0]+eps)//dx, (coordinate[1]+eps)//dy, (coordinate[2]+eps)//dz]
u = fem.Function(V)
for i, map_i in enumerate(map_vec):
    u.x.array[i] = values[map_i[0], map_i[1], map_i[2]]   

with io.XDMFFile(msh.comm, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    xdmf.write_function(u)