| /* --------------------------------------------------------------------- |
| Author: H. Carter Edwards |
| hcedwar@sandia.gov |
| |
| Copyright: Copyright (C) 1997 H. Carter Edwards |
| Graduate Student |
| University of Texas |
| |
| Re-release: Copyright (C) 2011-2012 H. Carter Edwards |
| |
| Purpose: Domain paritioning based upon Hilbert Space-Filling Curve |
| ordering. |
| |
| License: Re-release under the less-restrictive CLAMR software terms. |
| Permitted by email with H. Carter Edwards on 9/13/2011 |
| |
| Disclaimer: |
| |
| These routines comes with ABSOLUTELY NO WARRANTY; |
| This is free software, and you are welcome to redistribute it |
| under certain conditions. See License terms in file 'LICENSE'. |
| --------------------------------------------------------------------- */ |
| |
| /*---------------------------------------------------------------------- |
| Description: |
| Inverse of the Hilbert Space-Filling Curve Map from a 2D or 3D |
| domain to the 1D domain. Two different 2D and 3D domains are |
| supported. |
| |
| For the routines 'hsfc2d' and 'hsfc3d' the 2D and 3D domains are |
| defined as follows. |
| Note that |
| * 0 is the minimum value of an unsigned integer |
| * ~(0u) is the maximum value of an unsigned integer - all bits set |
| thus the 2D and 3D domains are |
| * [0,~(0u)] x [0,~(0u)] |
| * [0,~(0u)] x [0,~(0u)] x [0,~(0u)] |
| respectively. |
| |
| For the routines 'fhsfc2d' and 'fhsfc3d' the 2D and 3D domains are |
| defines as: |
| * [0.0,1.0] x [0.0,1.0] |
| * [0.0,1.0] x [0.0,1.0] x [0.0,1.0] |
| respectively. |
| |
| The 1D domain is a multiword (array of unsigned integers) key. |
| This key is essentially an unsigned integer of an arbitrary |
| number of bits. The most significant bit is the leading bit |
| of the first (0th) word of the key. The least significant |
| bit is the trailing bit of the last word. |
| |
| ----------------------------------------------------------------------*/ |
| |
| #include <stdlib.h> |
| #include <limits.h> |
| |
| /* Bits per unsigned word */ |
| |
| #define MaxBits ( sizeof(unsigned) * CHAR_BIT ) |
| |
| /*--------------------------------------------------------------------*/ |
| /* 2D Hilbert Space-filling curve */ |
| |
| void hsfc2d( |
| unsigned coord[] , /* IN: Normalized integer coordinates */ |
| unsigned nkey , /* IN: Word length of key */ |
| unsigned key[] ) /* OUT: space-filling curve key */ |
| { |
| static int init = 0 ; |
| static unsigned char gray_inv[ 2 * 2 ] ; |
| |
| const unsigned NKey = ( 2 < nkey ) ? 2 : (nkey) ; |
| const unsigned NBits = ( MaxBits * NKey ) / 2 ; |
| |
| unsigned i ; |
| unsigned char order[2+2] ; |
| unsigned char reflect ; |
| |
| /* GRAY coding */ |
| |
| if ( ! init ) { |
| unsigned char gray[ 2 * 2 ] ; |
| register unsigned k ; |
| register unsigned j ; |
| |
| gray[0] = 0 ; |
| for ( k = 1 ; k < sizeof(gray) ; k <<= 1 ) { |
| for ( j = 0 ; j < k ; j++ ) gray[k+j] = k | gray[k-(j+1)] ; |
| } |
| for ( k = 0 ; k < sizeof(gray) ; k++ ) gray_inv[ gray[k] ] = k ; |
| init = 1 ; |
| } |
| |
| /* Zero out the key */ |
| |
| for ( i = 0 ; i < NKey ; ++i ) key[i] = 0 ; |
| |
| order[0] = 0 ; |
| order[1] = 1 ; |
| reflect = ( 0 << 0 ) | ( 0 ); |
| |
| for ( i = 1 ; i <= NBits ; i++ ) { |
| const unsigned s = MaxBits - i ; |
| const unsigned c = gray_inv[ reflect ^ ( |
| ( ( ( coord[0] >> s ) & 01 ) << order[0] ) | |
| ( ( ( coord[1] >> s ) & 01 ) << order[1] ) ) ]; |
| |
| const unsigned off = 2 * i ; /* Bit offset */ |
| const unsigned which = off / MaxBits ; /* Which word to update */ |
| const unsigned shift = MaxBits - off % MaxBits ; /* Which bits to update */ |
| |
| /* Set the two bits */ |
| |
| if ( shift == MaxBits ) { /* Word boundary */ |
| key[ which - 1 ] |= c ; |
| } |
| else { |
| key[ which ] |= c << shift ; |
| } |
| |
| /* Determine the recursive quadrant */ |
| |
| switch( c ) { |
| case 3: |
| reflect ^= 03 ; |
| case 0: |
| order[2+0] = order[0] ; |
| order[2+1] = order[1] ; |
| order[0] = order[2+1] ; |
| order[1] = order[2+0] ; |
| break ; |
| } |
| } |
| } |
| |
| /*--------------------------------------------------------------------*/ |
| /* 3D Hilbert Space-filling curve */ |
| |
| void hsfc3d( |
| unsigned coord[] , /* IN: Normalized integer coordinates */ |
| unsigned nkey , /* IN: Word length of 'key' */ |
| unsigned key[] ) /* OUT: space-filling curve key */ |
| { |
| static int init = 0 ; |
| static unsigned char gray_inv[ 2*2*2 ] ; |
| |
| const unsigned NKey = ( 3 < nkey ) ? 3 : (nkey) ; |
| const unsigned NBits = ( MaxBits * NKey ) / 3 ; |
| |
| unsigned i ; |
| unsigned char axis[3+3] ; |
| |
| /* GRAY coding */ |
| |
| if ( ! init ) { |
| unsigned char gray[ 2*2*2 ] ; |
| register unsigned k ; |
| register unsigned j ; |
| |
| gray[0] = 0 ; |
| for ( k = 1 ; k < sizeof(gray) ; k <<= 1 ) { |
| for ( j = 0 ; j < k ; j++ ) gray[k+j] = k | gray[k-(j+1)] ; |
| } |
| for ( k = 0 ; k < sizeof(gray) ; k++ ) gray_inv[ gray[k] ] = k ; |
| init = 1 ; |
| } |
| |
| /* Zero out the key */ |
| |
| for ( i = 0 ; i < NKey ; ++i ) key[i] = 0 ; |
| |
| axis[0] = 0 << 1 ; |
| axis[1] = 1 << 1 ; |
| axis[2] = 2 << 1 ; |
| |
| for ( i = 1 ; i <= NBits ; i++ ) { |
| const unsigned s = MaxBits - i ; |
| const unsigned c = gray_inv[ |
| (((( coord[ axis[0] >> 1 ] >> s ) ^ axis[0] ) & 01 ) << 0 ) | |
| (((( coord[ axis[1] >> 1 ] >> s ) ^ axis[1] ) & 01 ) << 1 ) | |
| (((( coord[ axis[2] >> 1 ] >> s ) ^ axis[2] ) & 01 ) << 2 ) ]; |
| unsigned n ; |
| |
| /* Set the 3bits */ |
| |
| for ( n = 0 ; n < 3 ; ++n ) { |
| const unsigned bit = 01 & ( c >> ( 2 - n ) ); /* Bit value */ |
| const unsigned off = 3 * i + n ; /* Bit offset */ |
| const unsigned which = off / MaxBits ; /* Which word */ |
| const unsigned shift = MaxBits - off % MaxBits ; /* Which bits */ |
| |
| if ( MaxBits == shift ) { /* Word boundary */ |
| key[ which - 1 ] |= bit ; |
| } |
| else { |
| key[ which ] |= bit << shift ; |
| } |
| } |
| |
| /* Determine the recursive quadrant */ |
| |
| axis[3+0] = axis[0] ; |
| axis[3+1] = axis[1] ; |
| axis[3+2] = axis[2] ; |
| |
| switch( c ) { |
| case 0: |
| axis[0] = axis[3+2]; |
| axis[1] = axis[3+1]; |
| axis[2] = axis[3+0]; |
| break ; |
| case 1: |
| axis[0] = axis[3+0]; |
| axis[1] = axis[3+2]; |
| axis[2] = axis[3+1]; |
| break ; |
| case 2: |
| axis[0] = axis[3+0]; |
| axis[1] = axis[3+1]; |
| axis[2] = axis[3+2]; |
| break ; |
| case 3: |
| axis[0] = axis[3+2] ^ 01 ; |
| axis[1] = axis[3+0] ^ 01 ; |
| axis[2] = axis[3+1]; |
| break ; |
| case 4: |
| axis[0] = axis[3+2]; |
| axis[1] = axis[3+0] ^ 01 ; |
| axis[2] = axis[3+1] ^ 01 ; |
| break ; |
| case 5: |
| axis[0] = axis[3+0]; |
| axis[1] = axis[3+1]; |
| axis[2] = axis[3+2]; |
| break ; |
| case 6: |
| axis[0] = axis[3+0]; |
| axis[1] = axis[3+2] ^ 01 ; |
| axis[2] = axis[3+1] ^ 01 ; |
| break ; |
| case 7: |
| axis[0] = axis[3+2] ^ 01 ; |
| axis[1] = axis[3+1]; |
| axis[2] = axis[3+0] ^ 01 ; |
| break ; |
| default: |
| exit(-1); |
| } |
| } |
| } |
| |
| /*--------------------------------------------------------------------*/ |
| |
| void fhsfc2d( |
| double coord[] , /* IN: Normalized floating point coordinates */ |
| unsigned nkey , /* IN: Word length of key */ |
| unsigned key[] ) /* OUT: space-filling curve key */ |
| { |
| const double imax = ~(0u); |
| unsigned c[2] ; |
| c[0] = coord[0] * imax ; |
| c[1] = coord[1] * imax ; |
| hsfc2d( c , nkey , key ); |
| } |
| |
| void fhsfc3d( |
| double coord[] , /* IN: Normalized floating point coordinates */ |
| unsigned nkey , /* IN: Word length of key */ |
| unsigned key[] ) /* OUT: space-filling curve key */ |
| { |
| const double imax = ~(0u); |
| unsigned c[3] ; |
| c[0] = coord[0] * imax ; |
| c[1] = coord[1] * imax ; |
| c[2] = coord[2] * imax ; |
| hsfc3d( c , nkey , key ); |
| } |
| |