FTreeCoordinate.hpp 15 KB
Newer Older
1
// ===================================================================================
2 3 4 5 6 7 8 9 10 11 12 13 14
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
15
// ===================================================================================
16 17
#ifndef FTREECOORDINATE_HPP
#define FTREECOORDINATE_HPP
18

COULAUD Olivier's avatar
COULAUD Olivier committed
19
#include <string>
20

21
#include "../Utils/FGlobal.hpp"
22
#include "../Utils/FMath.hpp"
23

24 25
#include "../Components/FAbstractSerializable.hpp"

26 27 28 29 30 31 32
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FTreeCoordinate
* Please read the license
*
* This class represents tree coordinate. It is used to save
* the position in "box unit" (not system/space unit!).
33 34
* It is directly related to morton index, as interleaves
* bits from this coordinate make the morton index
35
*/
36
class FTreeCoordinate : public FAbstractSerializable {
37
private:
38
    int data[3];	//< all box-th position
39 40 41

public:	
    /** Default constructor (position = {0,0,0})*/
42 43
    FTreeCoordinate() {
        data[0] = data[1] = data[2] = 0;
44 45
    }

46 47 48 49 50
    /** Default constructor (position = {0,0,0})*/
    explicit FTreeCoordinate(const MortonIndex mindex, const int indexLevel) {
        setPositionFromMorton(mindex, indexLevel);
    }

51 52 53 54 55 56
    /**
        * Default constructor
        * @param inX the x
        * @param inY the y
        * @param inZ the z
        */
57
    explicit FTreeCoordinate(const int inX,const int inY,const int inZ) {
58 59 60
        data[0] = inX;
        data[1] = inY;
        data[2] = inZ;
61 62
    }

63 64 65 66 67 68
    explicit FTreeCoordinate(const int inPosition[3]) {
        data[0] = inPosition[0];
        data[1] = inPosition[1];
        data[2] = inPosition[2];
    }

69
    /**
70 71 72
    * Copy constructor
    * @param other the source class to copy
    */
73 74 75 76
    FTreeCoordinate(const FTreeCoordinate& other) {
        data[0] = other.data[0];
        data[1] = other.data[1];
        data[2] = other.data[2];
77 78
    }

berenger-bramas's avatar
berenger-bramas committed
79 80 81 82
    /**
        * Copy constructor
        * @param other the source class to copy
        */
83 84 85 86
    FTreeCoordinate(const FTreeCoordinate& other, const int inOffset) {
        data[0] = other.data[0] + inOffset;
        data[1] = other.data[1] + inOffset;
        data[2] = other.data[2] + inOffset;
berenger-bramas's avatar
berenger-bramas committed
87 88
    }

89
    /**
90 91 92 93
    * Copy constructor
    * @param other the source class to copy
    * @return this a reference to the current object
    */
94
    FTreeCoordinate& operator=(const FTreeCoordinate& other){
95 96 97
        data[0] = other.data[0];
        data[1] = other.data[1];
        data[2] = other.data[2];
98 99 100 101
        return *this;
    }

    /**
102
    * Position setter
103 104 105
        * @param inX the new x
        * @param inY the new y
        * @param inZ the new z
106
    */
107
    void setPosition(const int inX,const int inY,const int inZ){
108 109 110
        data[0] = inX;
        data[1] = inY;
        data[2] = inZ;
111 112 113
    }

    /**
114
    * X Getter
115
        * @return data[0]
116
    */
117
    int getX() const{
118
        return data[0];
119 120 121
    }

    /**
122
    * Y Getter
123
        * @return data[1]
124
    */
125
    int getY() const{
126
        return data[1];
127 128 129
    }

    /**
130
    * Z Getter
131
        * @return data[2]
132
    */
133
    int getZ() const{
134
        return data[2];
135 136 137
    }

    /**
138 139 140
    * X Setter, simply change x position
    * @param the new x
    */
141
    void setX(const int inX){
142
        data[0] = inX;
143 144 145
    }

    /**
146 147 148
    * Y Setter, simply change y position
    * @param the new y
    */
149
    void setY(const int inY){
150
        data[1] = inY;
151 152 153
    }

    /**
154 155 156
    * Z Setter, simply change z position
    * @param the new z
    */
157
    void setZ(const int inZ){
158
        data[2] = inZ;
159 160 161
    }

    /**
162 163 164 165 166
    * To get the morton index of the current position
    * @complexity inLevel
    * @param inLevel the level of the component
    * @return morton index
    */
167 168 169 170
    MortonIndex getMortonIndex(const int inLevel) const{
        MortonIndex index = 0x0LL;
        MortonIndex mask = 0x1LL;
        // the ordre is xyz.xyz...
171 172 173
        MortonIndex mx = data[0] << 2;
        MortonIndex my = data[1] << 1;
        MortonIndex mz = data[2];
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197

        for(int indexLevel = 0; indexLevel < inLevel ; ++indexLevel){
            index |= (mz & mask);
            mask <<= 1;
            index |= (my & mask);
            mask <<= 1;
            index |= (mx & mask);
            mask <<= 1;

            mz <<= 2;
            my <<= 2;
            mx <<= 2;
        }

        return index;
    }

    /** This function set the position of the current object using a morton index
          * @param inIndex the morton index to compute position
          * @param the level of the morton index
          */
    void setPositionFromMorton(MortonIndex inIndex, const int inLevel){
        MortonIndex mask = 0x1LL;

198 199 200
        data[0] = 0;
        data[1] = 0;
        data[2] = 0;
201 202

        for(int indexLevel = 0; indexLevel < inLevel ; ++indexLevel){
203
            data[2] |= int(inIndex & mask);
204
            inIndex >>= 1;
205
            data[1] |= int(inIndex & mask);
206
            inIndex >>= 1;
207
            data[0] |= int(inIndex & mask);
208 209 210 211 212 213 214 215 216 217

            mask <<= 1;
        }

    }

    /** Test equal operator
          * @param other the coordinate to compare
          * @return true if other & current object have same position
          */
218
    bool operator==(const FTreeCoordinate& other) const {
219
        return data[0] == other.data[0] && data[1] == other.data[1] && data[2] == other.data[2];
220
    }
221

222 223 224 225 226 227 228 229
    /** Test equal operator
          * @param other the coordinate to compare
          * @return true if other & current object have same position
          */
    bool equals(const int inX, const int inY, const int inZ) const {
        return data[0] == inX && data[1] == inY && data[2] == inZ;
    }

230 231 232 233
    /** To test difference
      *
      */
    bool operator!=(const FTreeCoordinate& other) const{
234
        return data[0] != other.data[0] || data[1] != other.data[1] || data[2] != other.data[2];
235
    }
236 237 238 239 240 241 242

    /**
     * Operator stream FTreeCoordinate to std::ostream
     * This can be used to simply write out a tree coordinate
     * @param[in,out] output where to write the coordinate
     * @param[in] inCoordinate the coordinate to write out
     * @return the output for multiple << operators
243
     */
244 245
    template <class StreamClass>
    friend StreamClass& operator<<(StreamClass& output, const FTreeCoordinate& inCoordinate){
246 247 248 249
        output << "(" <<  inCoordinate.getX() << ", " << inCoordinate.getY() << ", " << inCoordinate.getZ() <<")";
        return output;  // for multiple << operators.
    }

250
    /** Save current object */
BRAMAS Berenger's avatar
BRAMAS Berenger committed
251 252
    template <class BufferWriterClass>
    void save(BufferWriterClass& buffer) const {
253
        buffer << data[0] << data[1] << data[2];
254 255
    }
    /** Retrieve current object */
BRAMAS Berenger's avatar
BRAMAS Berenger committed
256 257
    template <class BufferReaderClass>
    void restore(BufferReaderClass& buffer) {
258
        buffer >> data[0] >> data[1] >> data[2];
259
    }
260

261
    /** To know the size when we save it */
262 263
    FSize getSavedSize() const {
        return FSize(sizeof(data[0]) + sizeof(data[1]) + sizeof(data[2]));
264 265
    }

266

267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
    static std::string MortonToBinary(MortonIndex index, int level){
        std::string str;
        int bits = 1 << ((level * 3) - 1);
        int dim = 0;
        while(bits){
            if(index & bits) str.append("1");
            else str.append("0");
            bits >>= 1;
            // we put a dot each 3 values
            if(++dim == 3){
                str.append(".");
                dim = 0;
            }
        }
        return str;
    }
283 284

    /* @brief Compute the index of the cells in neighborhood of a given cell
285 286
   * @param OtreeHeight Height of the Octree
   * @param indexes target array to store the MortonIndexes computed
287
   * @param indexInArray store
288
   */
289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
    int getNeighborsIndexes(const int OctreeHeight, MortonIndex indexes[26], int indexInArray[26]) const{
        int idxNeig = 0;
        int limite = 1 << (OctreeHeight - 1);
        // We test all cells around
        for(int idxX = -1 ; idxX <= 1 ; ++idxX){
            if(!FMath::Between(this->getX() + idxX,0, limite)) continue;

            for(int idxY = -1 ; idxY <= 1 ; ++idxY){
                if(!FMath::Between(this->getY() + idxY,0, limite)) continue;

                for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
                    if(!FMath::Between(this->getZ() + idxZ,0, limite)) continue;

                    // if we are not on the current cell
                    if( idxX || idxY || idxZ ){
                        const FTreeCoordinate other(this->getX() + idxX, this->getY() + idxY, this->getZ() + idxZ);
                        indexes[ idxNeig ] = other.getMortonIndex(OctreeHeight - 1);
                        indexInArray[ idxNeig ] = ((idxX+1)*3 + (idxY+1)) * 3 + (idxZ+1);
                        ++idxNeig;
                    }
                }
            }
        }
        return idxNeig;
313 314
    }

315 316

    /* @brief Compute the indexes of the neighborhood of the calling cell
317 318 319
   * @param OtreeHeight Height of the Octree
   * @param indexes target array to store the MortonIndexes computed
   */
320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
    int getNeighborsIndexes(const int OctreeHeight, MortonIndex indexes[26]) const{
        int idxNeig = 0;
        int limite = 1 << (OctreeHeight - 1);
        // We test all cells around
        for(int idxX = -1 ; idxX <= 1 ; ++idxX){
            if(!FMath::Between(this->getX() + idxX,0, limite)) continue;

            for(int idxY = -1 ; idxY <= 1 ; ++idxY){
                if(!FMath::Between(this->getY() + idxY,0, limite)) continue;

                for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
                    if(!FMath::Between(this->getZ() + idxZ,0, limite)) continue;

                    // if we are not on the current cell
                    if( idxX || idxY || idxZ ){
                        const FTreeCoordinate other(this->getX() + idxX, this->getY() + idxY, this->getZ() + idxZ);
                        indexes[ idxNeig ] = other.getMortonIndex(OctreeHeight - 1);
                        ++idxNeig;
                    }
                }
            }
        }
        return idxNeig;
343 344
    }

345 346
    int getInteractionNeighbors(const int inLevel, MortonIndex inNeighbors[/*189+26+1*/216], int inNeighborsPosition[/*189+26+1*/216],
                            const int neighSeparation = 1) const{
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
        // Then take each child of the parent's neighbors if not in directNeighbors
        // Father coordinate
        const FTreeCoordinate parentCell(this->getX()>>1,this->getY()>>1,this->getZ()>>1);

        // Limite at parent level number of box (split by 2 by level)
        const int limite = FMath::pow2(inLevel-1);

        int idxNeighbors = 0;
        // We test all cells around
        for(int idxX = -1 ; idxX <= 1 ; ++idxX){
            if(!FMath::Between(parentCell.getX() + idxX,0,limite)) continue;

            for(int idxY = -1 ; idxY <= 1 ; ++idxY){
                if(!FMath::Between(parentCell.getY() + idxY,0,limite)) continue;

                for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
                    if(!FMath::Between(parentCell.getZ() + idxZ,0,limite)) continue;

                    // if we are not on the current cell
366
                    if(neighSeparation<1 || idxX || idxY || idxZ ){
367 368 369 370 371 372 373 374 375 376
                        const FTreeCoordinate otherParent(parentCell.getX() + idxX,parentCell.getY() + idxY,parentCell.getZ() + idxZ);
                        const MortonIndex mortonOther = otherParent.getMortonIndex(inLevel-1);

                        // For each child
                        for(int idxCousin = 0 ; idxCousin < 8 ; ++idxCousin){
                            const int xdiff  = ((otherParent.getX()<<1) | ( (idxCousin>>2) & 1)) - this->getX();
                            const int ydiff  = ((otherParent.getY()<<1) | ( (idxCousin>>1) & 1)) - this->getY();
                            const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1)) - this->getZ();

                            // Test if it is a direct neighbor
377
                            if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
378 379 380 381 382 383 384 385 386 387 388
                                // add to neighbors
                                inNeighborsPosition[idxNeighbors] = ((( (xdiff+3) * 7) + (ydiff+3))) * 7 + zdiff + 3;
                                inNeighbors[idxNeighbors++] = (mortonOther << 3) | idxCousin;
                            }
                        }
                    }
                }
            }
        }

        return idxNeighbors;
389 390
    }

391
    int getInteractionNeighbors(const int inLevel, MortonIndex inNeighbors[/*189+26+1*/216], const int neighSeparation = 1) const{
392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410
        // Then take each child of the parent's neighbors if not in directNeighbors
        // Father coordinate
        const FTreeCoordinate parentCell(this->getX()>>1,this->getY()>>1,this->getZ()>>1);

        // Limite at parent level number of box (split by 2 by level)
        const int limite = FMath::pow2(inLevel-1);

        int idxNeighbors = 0;
        // We test all cells around
        for(int idxX = -1 ; idxX <= 1 ; ++idxX){
            if(!FMath::Between(parentCell.getX() + idxX,0,limite)) continue;

            for(int idxY = -1 ; idxY <= 1 ; ++idxY){
                if(!FMath::Between(parentCell.getY() + idxY,0,limite)) continue;

                for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
                    if(!FMath::Between(parentCell.getZ() + idxZ,0,limite)) continue;

                    // if we are not on the current cell
411
                    if(neighSeparation<1 || idxX || idxY || idxZ ){
412 413 414 415 416 417 418 419 420 421
                        const FTreeCoordinate otherParent(parentCell.getX() + idxX,parentCell.getY() + idxY,parentCell.getZ() + idxZ);
                        const MortonIndex mortonOther = otherParent.getMortonIndex(inLevel-1);

                        // For each child
                        for(int idxCousin = 0 ; idxCousin < 8 ; ++idxCousin){
                            const int xdiff  = ((otherParent.getX()<<1) | ( (idxCousin>>2) & 1)) - this->getX();
                            const int ydiff  = ((otherParent.getY()<<1) | ( (idxCousin>>1) & 1)) - this->getY();
                            const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1)) - this->getZ();

                            // Test if it is a direct neighbor
422
                            if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
423 424 425 426 427 428 429 430 431 432 433
                                // add to neighbors
                                inNeighbors[idxNeighbors++] = (mortonOther << 3) | idxCousin;
                            }
                        }
                    }
                }
            }
        }

        return idxNeighbors;
    }
434

435 436 437 438 439 440
};



#endif //FTREECOORDINATE_HPP

441