[sbaffini]

Dear all,

can you suggest a dummy proof book (or alternative material) on the space filling curves, especially their use for parallel grid partitioning?

Thanks

[mprinkey]

Berger at NASA has done work on this.

http://people.nas.nasa.gov/aftosmis/...a2004_1232.pdf

The approach is pretty simple:

- Assign a space-filling curve index to each entity based on its geometric location
- Sort the entities based on the space-filling curve index
- Walk through the in the sorted list, assigning an ordered block of entities to each PE

[sbaffini]

Dear Michael thank you for your help.

Thus, long story short, all the matter is using a black box, say hilbert_map(x,y,z), which for every different location in a grid will return a distinguished unique index?

What i don't get (probably because i need to delve more in it) is how for a general grid i can do this, considering that:

- i have to choose a finite scale for the sfc to stop at (do i?)
- i imagine a finite scale sfc as a specific geometric object which might or not unambiguously mark every generic cell in my domain

[mprinkey]

Normalization of the domain is pretty easy. Generally you want your smallest length scale choice to be isotropic, unless you have a good reason to do otherwise. I always take the largest coordinate extent of the domain as the maximum length scale and pick the smallest scale based on the number of bits per coordinate direction I have in the SFC index. So, then it becomes a question of dynamic range of our index variable type.

If you use a 32-bit integer (signed or unsigned) for your SFC index in 3D, you may well run into dynamic range issues. That only provides 10 bits per coordinate direction (assuming an isotropic SFC), so you have dx of 1mm for 1m domain, or a billion cell 3D Cartesian mesh to drop entities into. That may be a problem, but probably not as often as you might think.

However, if you use 64-bit integers in 3D instead, you have 21 bits per coordinate direction. That is a dynamic range ~1:2,000,000. So, it represents scales from 1 millimeter to 2 kilometers. I've never used anything but 64-bit integers, even though 32-bit would generally be sufficient.

To your last point, I think you go to far in asserting that the SFC index must be unique to be useful. Not at all. You can treat them as "bins" into which to throw objects. In fact, this technique is sometimes called "spatial hashing" because it does produce something like a hash. But there is no guarantee that the hash is unique and therefore "collision" treatments may be necessary. With 64-bit integers, the chance of collision is very low unless you are doing something pathological. Moreover, the sorted list with non-unique entities is STILL spatially ordered leading to good cache coherence, reasonable parallel partitions (and sane/easy RE-partitions if you are doing adaption or particle tracking), and all of the other good stuff. It just so happens that a few entity indices are not unique...and it won't matter.

The SFC index just helps you put things in better order and better bunches. It doesn't replace proper unique indexing and the rest of the book keeping that comes with it.

One final point. If you use 64-bit signed integers, make sure that you use 0 offset (so C/C++ style counting). If you start the first bin as 1, the last bin (2^63) will overflow. 64-bit signed ints are only good to 2^63-1. So, either count C-style or use unsigned 64-bit ints.

[sbaffini]

Dear Michael, thank you again, you really enlightened me the SFC on several aspects.

Can then we say that a SFC is a (binning) tool that helps you ordering in 1D a given space partition?

Then, for a given (bin) partition, according to the resolution coming from the spatial (index) coordinates and the case dependent object clustering, one might or not end up with one or more objects in the same partition (bin)? Then these few cases are trivially resolved by heuristics?

I guess i now have to learn how to move from the Z order to the more appealing Hilbert one

[mprinkey]

Originally Posted by sbaffini

Can then we say that a SFC is a (binning) tool that helps you ordering in 1D a given space partition?

Yes, that is a good way to look at it, and that is the way it is often described. It is a mechanism to map multidimensional data into one dimension while carrying over some sort of spatial coherence.

Originally Posted by sbaffini

Then, for a given (bin) partition, according to the resolution coming from the spatial (index) coordinates and the case dependent object clustering, one might or not end up with one or more objects in the same partition (bin)? Then these few cases are trivially resolved by heuristics?

Yes. But I'd argue that the "binning" happens when you cast the spatial locations (double x, double y, double z) into integer space bins (int i, int j, int k) via the normalization I discussed above. The act of converting (i,j,k) into the SFC index is always one-to-one. If a collision happens in SFC index, it will occur in (i,j,k) too. SFC simply orders the discrete bins.

Originally Posted by sbaffini

I guess i now have to learn how to move from the Z order to the more appealing Hilbert one

I would not be in much of a hurry to use other curves. I slaved through implementing other ordering schemes and didn't see significant improvements over Z. The multi-resolution nature seems to be the powerful component more so than the particular stencil. Z ordering is really very easy to implement using bit twiddling by the way:

http://www.forceflow.be/2013/10/07/m...plementations/

[sbaffini]

Dear Michael,

despite your suggetsion, i really wanted to delve into the Hilbert order and i came up with the following Fortran implementation, which i adapted from a template found in:

G. Zumbusch: Parallel Multilevel Methods - Adaptive Mesh Refinement and Load Balancing.

I made some tests and it seems to work for the most obvious cases. However, the original reference had a 12 in fortran position (9,7) of the matrix state, which appeared to be wrong to me. Still, i had no difference in the results of the obvious cases when using 21 as below.

I was wandering if you could point to me some obvious mistakes i could have done in the implementation below.

Thank you very much.

Note that "ishift" in Fortran performs a left shift for positive arguments and a right shift for negative ones. 0 padding is done in all the cases.

program hilbert_sfc
implicit none
real*4 :: x, y, z
integer*4 :: level, box, lh, max_level
integer*4 :: i, j, k
integer*8 :: c
integer*4, dimension(24,8) :: state, curve

state = transpose(reshape(&

[ 16, 8, 8, 6, 6, 11, 11, 19 ,& !1
18, 10, 10, 7, 7, 9, 9, 17 ,& !2
12, 17, 17, 4, 4, 18, 18, 15 ,& !3
14, 19, 19, 5, 5, 16, 16, 13 ,& !4
9, 20, 20, 2, 2, 23, 23, 10 ,& !5
11, 22, 22, 3, 3, 21, 21, 8 ,& !6
21, 13, 13, 0, 0, 14, 14, 22 ,& !7
23, 15, 15, 1, 1, 12, 12, 20 ,& !8
0, 16, 16, 14, 14, 21, 21, 5 ,& !9
4, 20, 20, 15, 15, 17, 17, 1 ,& !10
18, 1, 1, 12, 12, 4, 4, 23 ,& !11
22, 5, 5, 13, 13, 0, 0, 19 ,& !12
17, 2, 2, 10, 10, 7, 7, 20 ,& !13
21, 6, 6, 11, 11, 3, 3, 16 ,& !14
3, 19, 19, 8, 8, 22, 22, 6 ,& !15
7, 23, 23, 9, 9, 18, 18, 2 ,& !16
0, 8, 8, 22, 22, 13, 13, 3 ,& !17
2, 12, 12, 23, 23, 9, 9, 1 ,& !18
10, 1, 1, 20, 20, 2, 2, 15 ,& !19
14, 3, 3, 21, 21, 0, 0, 11 ,& !20
9, 4, 4, 18, 18, 7, 7, 12 ,& !21
13, 6, 6, 19, 19, 5, 5, 8 ,& !22
5, 11, 11, 16, 16, 14, 14, 6 ,& !23
7, 15, 15, 17, 17, 10, 10, 4 ],& !24

(/8, 24/)))

curve = transpose(reshape(&
[ 0, 7, 1, 6, 3, 4, 2, 5 ,& !1
7, 0, 6, 1, 4, 3, 5, 2 ,& !2
3, 4, 0, 7, 2, 5, 1, 6 ,& !3
4, 3, 7, 0, 5, 2, 6, 1 ,& !4
1, 6, 2, 5, 0, 7, 3, 4 ,& !5
6, 1, 5, 2, 7, 0, 4, 3 ,& !6
2, 5, 3, 4, 1, 6, 0, 7 ,& !7
5, 2, 4, 3, 6, 1, 7, 0 ,& !8
0, 1, 3, 2, 7, 6, 4, 5 ,& !9
7, 6, 4, 5, 0, 1, 3, 2 ,& !10
3, 0, 2, 1, 4, 7, 5, 6 ,& !11
4, 7, 5, 6, 3, 0, 2, 1 ,& !12
1, 2, 0, 3, 6, 5, 7, 4 ,& !13
6, 5, 7, 4, 1, 2, 0, 3 ,& !14
2, 3, 1, 0, 5, 4, 6, 7 ,& !15
5, 4, 6, 7, 2, 3, 1, 0 ,& !16
0, 1, 7, 6, 3, 2, 4, 5 ,& !17
7, 6, 0, 1, 4, 5, 3, 2 ,& !18
3, 0, 4, 7, 2, 1, 5, 6 ,& !19
4, 7, 3, 0, 5, 6, 2, 1 ,& !20
1, 2, 6, 5, 0, 3, 7, 4 ,& !21
6, 5, 1, 2, 7, 4, 0, 3 ,& !22
2, 3, 5, 4, 1, 0, 6, 7 ,& !23
5, 4, 2, 3, 6, 7, 1, 0 ],& !24
(/8, 24/)))

max_level = 21

x=1.0
y=0.0
z=0.0

lh=0
c=0

i=MIN(FLOOR(x*2**max_level),2**max_level-1)
j=MIN(FLOOR(y*2**max_level),2**max_level-1)
k=MIN(FLOOR(z*2**max_level),2**max_level-1)

DO level=max_level,0,-1

box=ior(ior((iand(ishft(i,-level),1)),(iand(ishft(j,-level+1),2))),(iand(ishft(k,-level+2),4)))

c = curve(lh+1,box+1)+c*8

lh = state(lh+1,box+1)

ENDDO

WRITE(*,*) 'index =', c

end program hilbert_sfc

[mprinkey]

Sorry, I cooked up my own algorithm for U-curves (aka C-curves). I don't remember all of the details, but I know that it was a pain to sort out. I am not familiar with the approach being used in your code and I honestly can't spend the time to crawl through all of those state/curve indices.

I don't have my 3D U-curve code in Fortran--I did port the 2D U-curve, but that likely won't help you much. I wrote them all originally in C. Here is my 3D U-curve C function--maybe it can help you. It starts with an initial template in T, identifies the octant at each level, and then rotates/mirrors (with XFORM) the template as it descends levels.

typedef unsigned long uint64;

#define XFORM(A,B,C,D,E,F,G,H)                      \
{TMP[0]=T[A];TMP[1]=T[B];TMP[2]=T[C];TMP[3]=T[D];   \
 TMP[4]=T[E];TMP[5]=T[F];TMP[6]=T[G];TMP[7]=T[H];   \
{char *pT=T;T=TMP;TMP=pT;}}

uint64 uCurve3DIndex (double x, double y, double z, int L, double xl, double xh,
                      double yl, double yh, double zl, double zh)
{
int l;
uint64 index = 0;
double Xm, Ym, Zm, Xl = xl, Xh = xh, Yl = yl, Yh = yh, Zl = zl, Zh = zh;

/* Using notation derived from:
        http://www.math.uwaterloo.ca/~wgilbert/Research/HilbertCurve/HilbertCurve.html
*/

char aT[8]={5,2,4,3,6,1,7,0};
char aTMP[8];
char delta;
char *T=aT,*TMP=aTMP;

/* Descend levels and add increments at each level */
for (l=0;l<L;l++)
    {
    index *= 8;
    Xm = 0.5*(Xl+Xh);
    Ym = 0.5*(Yl+Yh);
    Zm = 0.5*(Zl+Zh);

    /* Identify which octant contains the point and load delta from template. */
    if (x < Xm)
        {
        Xh = Xm;
        if (y < Ym)
            {
            Yh = Ym;
            if (z < Zm)
                { Zh = Zm; delta = T[0]; }
            else
                { Zl = Zm; delta = T[1]; }
            }
        else
            {
            Yl = Ym;
            if (z < Zm)
                { Zh = Zm; delta = T[2]; }
            else
                { Zl = Zm; delta = T[3]; }
            }
        }
    else
        {
        Xl = Xm;
        if (y < Ym)
            {
            Yh = Ym;
            if (z < Zm)
                { Zh = Zm; delta = T[4]; }
            else
                { Zl = Zm; delta = T[5]; }
            }
        else
            {
            Yl = Ym;
            if (z < Zm)
                { Zh = Zm; delta = T[6]; }
            else
                { Zl = Zm; delta = T[7]; }
            }
        }

    /* Rotate/mirror template based on octant */
    switch (delta)
        {
        case 0: XFORM(0,3,4,7,6,5,2,1) break;
        case 1:
        case 2: XFORM(0,7,6,1,2,5,4,3) break;
        case 3:
        case 4: XFORM(2,3,0,1,6,7,4,5) break;
        case 5:
        case 6: XFORM(4,3,2,5,6,1,0,7) break;
        case 7: XFORM(6,5,2,1,0,3,4,7) break;
        }

    /* Add increment for current level */
    index += delta;
    }
return index;
}
#undef XFORM

[sbaffini]

Dear Michael,

thank you again. Could you maybe confirm that in your switch on delta, cases 1, 3 and 5 have the same output than, respectively, 2, 4 and 6? Thanks

[sbaffini]

I finally solved the issue by using the tables of code_saturne:

http://code-saturne.org/viewvc/satur...27&view=markup

and starting the DO loop above from max_level-1.

I'm still unsure about the overall correctness but it seems to work now

[mprinkey]

I'm glad that you solved your problem. To answer your question, yes...case statements "fall through" and have to be ended with an explicit break statement.

[sbaffini]

Just to update on the matter. I was clearly unhappy with my previous black box solution, so i started building my own state tables by hand (for the 12 terminals case). It is by this process that i finally found the second mistake in the code above. The line:

lh = state(lh+1,box+1)

should be instead substituted by:

lh = state(lh+1,curve(lh+1,box+1)+1)

due to the way the states are listed in the matrix.

As a final question (i promise ), do you know Michael how the 3D Hilbert curve should behave for points distributed on a plane? Should it preserve all its features and, maybe, reduce to the 2D curve?

Thanks

[mprinkey]

I can't answer your question precisely. I don't think that you can easily say uCurve3D(i,j,fixedlK) = someFunction(uCurve2D(i,j),fixedK). Such a form may well exist. I just don't know it. Those sorts of things are not terribly useful anyway, because it couldn't say anything about planes that aren't aligned with a coordinate axis.

Now, with regard to properties in a more general sense, I think you can assume the same index distribution tendencies will hold due solely to the multiresolution nature of the curve--points that are close on any plane (or any 2D surface) embedded in a 3D space should tend to have close Hilbert index values. That is making a statement analogous to "points close together tend to fall into the same branches of an octtree."

I know that is wishy-washy, but it is the best I can offer. Honestly, since you have a working function, maybe just feed points from any arbitrary plane into it, and plot the results.

[sbaffini]

Thank you again. I was indeed making some tests and i simply found out that these curves, obviously, fo not work well for arbitrary point distributions, holes being the prime cause of self intersections.

To be honest, i really forced some very unusual cases, so i think i will be happy with my present implementation which, by the way, will be used on cartesian grids with local anisotropic refinement