CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Space Filling Curves

Register Blogs Community New Posts Updated Threads Search

Like Tree2Likes
  • 1 Post By mprinkey
  • 1 Post By mprinkey

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 15, 2015, 09:02
Default Space Filling Curves
  #1
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Old   January 15, 2015, 11:20
Default
  #2
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25
mprinkey will become famous soon enough
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 N_{entities} in the sorted list, assigning an ordered block of N_{entities}/N_{processes} entities to each PE
sbaffini likes this.
mprinkey is offline   Reply With Quote

Old   January 15, 2015, 17:20
Default
  #3
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Old   January 15, 2015, 21:28
Default
  #4
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25
mprinkey will become famous soon enough
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 likes this.
mprinkey is offline   Reply With Quote

Old   January 19, 2015, 18:41
Default
  #5
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Old   January 19, 2015, 19:09
Default
  #6
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25
mprinkey will become famous soon enough
Quote:
Originally Posted by sbaffini View Post
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.

Quote:
Originally Posted by sbaffini View Post
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.

Quote:
Originally Posted by sbaffini View Post
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/
mprinkey is offline   Reply With Quote

Old   January 28, 2015, 06:40
Default
  #7
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Old   January 28, 2015, 08:22
Default
  #8
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25
mprinkey will become famous soon enough
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.

Code:
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

Last edited by mprinkey; January 28, 2015 at 20:45.
mprinkey is offline   Reply With Quote

Old   January 28, 2015, 09:11
Default
  #9
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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 is offline   Reply With Quote

Old   January 28, 2015, 12:48
Default
  #10
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Old   January 28, 2015, 20:43
Default
  #11
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25
mprinkey will become famous soon enough
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.
mprinkey is offline   Reply With Quote

Old   January 30, 2015, 04:31
Default
  #12
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Old   January 30, 2015, 08:14
Default
  #13
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25
mprinkey will become famous soon enough
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.
mprinkey is offline   Reply With Quote

Old   January 30, 2015, 15:37
Default
  #14
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to 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
sbaffini is offline   Reply With Quote

Reply


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[ICEM] Error 2D mesh - Separate curves seen as one when associating bbmorales ANSYS Meshing & Geometry 3 June 11, 2013 12:56
Version 7.04 : disappearing feature curves and other woes nomad STAR-CCM+ 2 September 4, 2012 22:07
how to track/view the filling of resin inside a mold using viscid/laminar modelFluent gtuncol FLUENT 0 August 7, 2010 00:49
How to model the NR eqns in a domain with empty space Vasilis Main CFD Forum 1 April 14, 2009 05:35
Fatal error error writing to tmp No space left on device maka OpenFOAM Installation 2 April 3, 2006 09:48


All times are GMT -4. The time now is 13:57.