Numerical relativity - coordinates, discretization and shift vector

I would like to understand how, in the numerical implementation of the relativistic 3+1 decomposition, the spatial coordinates are evolved from slice to slice and what is their relation to the shift vector.

Let me elaborate on the problem. In 3+1 decomposition, the evolution of spacetime $mathcal{M}$ can be understood as evolving spatial fields on a hypersurface $Sigma$ that changes its topology and geometry. In other words, fields and coordinates are evolved simultaneously.
In the 3+1 framework, there exists a concept of a Eulerian observer.The tangent vector field to Eulerian observers’ worldlines is the vector field $boldsymbol{n}$ normal to $Sigma$.

Then, there is a concept of a coordinate time observer congruence, whose tangent vector field is $$boldsymbol{partial_{t}} = alpha boldsymbol{n} + boldsymbol{beta} equiv boldsymbol{m} + boldsymbol{beta},$$ where $alpha$ is the lapse and $boldsymbol{beta}$ is the shift vector.
The worldlines of coordinate time observers are curves along which the spatial coordinate values are constant. If the curve $gamma$ corresponds to $boldsymbol{partial_{t}}$, then solving the
$dot{gamma}=boldsymbol{partial_{t}}$ yields a family of curves $gamma(tau)=(t(tau),x^{i}_{0})$, where the spatial coordinates remain constant.

If I am correct, the coordinate values of the Eulerian observers spatial position can be recovered (by virtue of $boldsymbol{partial_{t}} =boldsymbol{m} + boldsymbol{beta}$) via $x^{i}_{boldsymbol{m}}(s)= x^{i}_{0} – x^{i}_{beta}(s)$.

Here is when I have a conceptual problem and I would like to know if my thinking is correct.

Let us assume we want to discretize the hypersurface, i.e. represent our hypersurface as a grid.
I will use a three dimensional box and for each point use indices $p=(i,j,k)$.

1. Now, what is the relation between these indices and values of some coordinates?
   Do we assign to grid indices the constant coordinate values? In other words, is it the case that $x^{i}(p)=(i,j,k) = x^{i}_{0}$ and remain constant through the computations?
2. I suspect that what could be done is assigning the initial spatial coordinate values to grid points (i,j,k), setting up the initial data (where the spatial initial coordinates of Eulerian and coordinate time observers coincide) and then, at each step in time ($trightarrow t + delta t)$, updating the spatial coordinate values of the grid through the shift vector as prescribed by $x^{i}_{boldsymbol{m}}(s)= x^{i}_{0} – x^{i}_{beta}(s)$. Then we would (of course?) have to assume that ANY other quantities are represented in the code in such a way, that they depend on the changing "background" coordinate grid.
3. Otherwise, we can assume that the numerical grid positions (i,j,k) always correspond to points that are occupied by coordinate time observers, therefore they never change, as by definition such observers have constant spatial coordinates.
   What is the role played by the shift vector in the simulation then, except for influencing the evolution of other quantities (spatial metric $gamma_{ij}$, extrinsic curvature $K_{ij}$)?
   Surely not changing the coordinate values of grid points anymore?
   The numerical grid stays constant, but it represents the evolving spacetime intrinsically in a sense, implicitly, through the constant spatial coordinate values of the moving (coordinate time) observers.
4. If question 2 has it correct, are the values of spatial coordinates on a grid position (i,j,k) monitored in some way by following the coordinate time observers?

It is difficult to express my question in a consistent way, so I hope it is understandable. I would be grateful for a detailed explanation of how it is done in practice and if my ideas are consistent.

Here is an illustration:
https://imgur.com/a/rAfH0f2

EDIT: I have decided to cross out points 2 and 4, as they seem not to be the case, based on answers in the post, the provided literature and after some consideration.