.. _program_listing_file_include_hilbert.h: Program Listing for File hilbert.h ================================== |exhale_lsh| :ref:`Return to documentation for file ` (``include/hilbert.h``) .. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS .. code-block:: cpp #ifndef hilbert_h #define hilbert_h #include #include #include "exafmm_t.h" #define EXAFMM_HILBERT 0 namespace exafmm_t { inline uint64_t levelOffset(int level) { return (((uint64_t)1 << 3 * level) - 1) / 7; } int getLevel(uint64_t i) { int level = -1; uint64_t offset = 0; while (i >= offset) { level++; offset += (uint64_t)1 << 3 * level; } return level; } uint64_t getParent(uint64_t i) { int level = getLevel(i); return (i - levelOffset(level)) / 8 + levelOffset(level-1); } uint64_t getChild(uint64_t i) { int level = getLevel(i); return (i - levelOffset(level)) * 8 + levelOffset(level+1); } int getOctant(uint64_t key, bool offset=true) { int level = getLevel(key); if (offset) key -= levelOffset(level); return key & 7; } uint64_t getKey(ivec3 iX, int level, bool offset=true) { #if EXAFMM_HILBERT int M = 1 << (level - 1); for (int Q=M; Q>1; Q>>=1) { int R = Q - 1; for (int d=0; d<3; d++) { if (iX[d] & Q) iX[0] ^= R; else { int t = (iX[0] ^ iX[d]) & R; iX[0] ^= t; iX[d] ^= t; } } } for (int d=1; d<3; d++) iX[d] ^= iX[d-1]; int t = 0; for (int Q=M; Q>1; Q>>=1) if (iX[2] & Q) t ^= Q - 1; for (int d=0; d<3; d++) iX[d] ^= t; #endif uint64_t i = 0; for (int l=0; l> 2*l; iX[1] |= (i & (uint64_t)1 << (3*l + 1)) >> (2*l + 1); iX[0] |= (i & (uint64_t)1 << (3*l + 2)) >> (2*l + 2); } #if EXAFMM_HILBERT int N = 2 << (level - 1); int t = iX[2] >> 1; for (int d=2; d>0; d--) iX[d] ^= iX[d-1]; iX[0] ^= t; for (int Q=2; Q!=N; Q<<=1) { int R = Q - 1; for (int d=2; d>=0; d--) { if (iX[d] & Q) iX[0] ^= R; else { t = (iX[0] ^ iX[d]) & R; iX[0] ^= t; iX[d] ^= t; } } } #endif return iX; } ivec3 get3DIndex(uint64_t i, int level) { ivec3 iX = 0; for (int l=0; l> 2*l; iX[1] |= (i & (uint64_t)1 << (3*l + 1)) >> (2*l + 1); iX[0] |= (i & (uint64_t)1 << (3*l + 2)) >> (2*l + 2); } #if EXAFMM_HILBERT int N = 2 << (level - 1); int t = iX[2] >> 1; for (int d=2; d>0; d--) iX[d] ^= iX[d-1]; iX[0] ^= t; for (int Q=2; Q!=N; Q<<=1) { int R = Q - 1; for (int d=2; d>=0; d--) { if (iX[d] & Q) iX[0] ^= R; else { t = (iX[0] ^ iX[d]) & R; iX[0] ^= t; iX[d] ^= t; } } } #endif return iX; } ivec3 get3DIndex(vec3 X, int level, vec3 x0, real_t r0) { vec3 Xmin = x0 - r0; real_t dx = 2 * r0 / (1 << level); ivec3 iX; for (int d=0; d<3; d++) { iX[d] = floor((X[d] - Xmin[d]) / dx); } return iX; } vec3 getCoordinates(ivec3 iX, int level, vec3 x0, real_t r0) { vec3 Xmin = x0 - r0; real_t dx = 2 * r0 / (1 << level); vec3 X; for (int d=0; d<3; d++) { X[d] = (iX[d] + 0.5) * dx + Xmin[d]; } return X; } } #endif