While transformation of 3D coordinates to a z-order curve was relatively straightforward (Efficient z-order transformation in Fortran) I am having difficulties to wrap my head around around the math for using different space-filling curves, for example Peano or Hilbert. Any hints on how an actual code to do the transformation might look like would be appreciated. The goal is to have a subroutine that takes xyz coordinates as an input with whatever normalization is necessary and returns the index of the space-filling curve.
subroutine(x, y, z, space_filling_index)
And related to this: I read that there are many ways to define Hilbert curves in 3D space, which would be best in terms of locality? If there is a definitive answer to this...
The application would be a reordering of cells in a Cartesian computational grid with the goal of increasing cache hits when a cell accesses its neighbors cells.
The Hilbert curve works by dividing recursively a cube (for 3D) using the same basic shape at every step, by rotating the curve so that the exit point of a sub-cube matches the entry point of the next cube.
A fantastic resource is the technical report Compact Hilbert Indices by C. Hamilton. The report also introduces the compact Hilbert index for non-cubic systems.
As I wrapped my head around it I wrote a blog post in 2015: Understanding the Hilbert curve with sample Python code for the Hilbert index and an illustration of the rotation of Hilbert "sub-cubes". As part of a particle-based Molecular Dynamics simulation code I wrote, I implemented the compact Hilbert index in Fortran, see here.
Re-discussing at length the details is "beyond the scope" of a SO answer I believe but the resources above should help you very much.