FTreeCoordinate.hpp 14.9 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
366
367
368
369
370
371
372
373
374
375
376
        // 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
                    if( idxX || idxY || idxZ ){
                        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
411
412
413
414
415
416
417
418
419
420
421
        // 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
                    if( idxX || idxY || idxZ ){
                        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