PABLO  0.1
PArallel Balanced Linear Octree
 All Classes Functions Variables Pages
Class_Octant_2D.tpp
1 
30 // =================================================================================== //
31 // CLASS SPECIALIZATION //
32 // =================================================================================== //
33 template<>
34 class Class_Octant<2>{
35  // ------------------------------------------------------------------------------- //
36  // FRIENDSHIPS ------------------------------------------------------------------- //
37  template<int dim> friend class Class_Local_Tree;
38  template<int dim> friend class Class_Para_Tree;
39 
40  // ------------------------------------------------------------------------------- //
41  // TYPEDEFS ----------------------------------------------------------------------- //
42 
43  typedef vector<Class_Octant<2> > OctantsType;
44  typedef vector<double> dvector;
45  typedef vector<uint32_t> u32vector;
46  typedef vector<vector<uint32_t> > u32vector2D;
47  typedef vector<vector<uint64_t> > u64vector2D;
48 
49  // ------------------------------------------------------------------------------- //
50  // MEMBERS ----------------------------------------------------------------------- //
51 
52 private:
53  uint32_t x;
54  uint32_t y;
55  uint8_t level;
56  int8_t marker;
57  bitset<12> info;
63  // ------------------------------------------------------------------------------- //
64  // CONSTRUCTORS AND OPERATORS----------------------------------------------------- //
65 
66 public:
67  Class_Octant(){
68  x = y = 0;
69  level = 0;
70  marker = 0;
71  for (int i=0; i<global2D.nfaces; i++){
72  info[i] = true;
73  }
74  };
75  Class_Octant(uint8_t level, uint32_t x, uint32_t y){
76  this->x = x;
77  this->y = y;
78  this->level = level;
79  marker = 0;
80  if (level==0){
81  for (int i=0; i<global2D.nfaces; i++){
82  info[i] = true;
83  }
84  }
85 
86  };
87 
88  Class_Octant(uint8_t level, uint32_t x, uint32_t y, bool bound){
89  this->x = x;
90  this->y = y;
91  this->level = level;
92  marker = 0;
93  if (level==0){
94  for (int i=0; i<global2D.nfaces; i++){
95  info[i] = bound;
96  }
97  }
98 
99  };
100 
101  Class_Octant(const Class_Octant<2> &octant){
102  x = octant.x;
103  y = octant.y;
104  level = octant.level;
105  marker = octant.marker;
106  info = octant.info;
107  };
108 
111  bool operator ==(const Class_Octant<2> & oct2){
112  bool check = true;
113  check = check && (x == oct2.x);
114  check = check && (y == oct2.y);
115  check = check && (level == oct2.level);
116  return check;
117  }
118 
119  // ------------------------------------------------------------------------------- //
120  // METHODS ----------------------------------------------------------------------- //
121 
122  // Basic Get/Set methods --------------------------------------------------------- //
123 
124 public:
128  uint32_t getX() const{return x;};
129 
133  uint32_t getY() const{return y;};
134 
138  uint32_t getZ() const{return 0;};
139 
143  uint8_t getLevel() const{return level;};
144 
148  int8_t getMarker() const{return marker;};
149 
154  bool getBound(uint8_t face) const{
155  return info[face];
156  };
157 
158 private:
159  void setBound(uint8_t face) { // Set if face is boundary
160  info[face] = true;
161  };
162 
163 public:
168  bool getPbound(uint8_t face) const{
169  return info[global2D.nfaces+face];
170  };
171 
175  bool getIsNewR() const{return info[8];};
176 
180  bool getIsNewC() const{return info[9];};
181 
185  bool getNotBalance() const{return info[10];};
186 
190  bool getBalance() const{return (!info[10]);};
191 
195  void setMarker(int8_t marker){
196  this->marker = marker;
197  };
198 
202  void setBalance(bool balance){
203  info[10] = balance;
204  };
205 
206 private:
207  void setLevel(uint8_t level){
208  this->level = level;
209  };
210  void setPbound(uint8_t face, bool flag){
211  info[global2D.nfaces+face] = flag;
212  };
213 
214  //-------------------------------------------------------------------------------- //
215  // Other Get/Set methods --------------------------------------------------------- //
216 
217 public:
218 
222  uint32_t getSize() const{
223  uint32_t size = uint32_t(pow(double(2),double(MAX_LEVEL_2D-level)));
224  return size;
225  };//TODO change pow to myPow
226 
230  uint32_t getArea() const{
231  uint32_t area =getSize();
232  return area;
233  };
234 
238  uint64_t getVolume() const{
239  uint64_t volume = uint64_t(pow(double(getSize()),2.0));
240  return volume;
241  };
242 
243  // ------------------------------------------------------------------------------- //
244 
248  dvector getCenter(){
249  double dh;
250 
251  dh = double(getSize())/2.0;
252  vector<double> center(3);
253 
254  center[0] = (double)x + dh;
255  center[1] = (double)y + dh;
256  center[2] = 0.0;
257  return center;
258  };
259 
260  // ------------------------------------------------------------------------------- //
261 
265  dvector getFaceCenter(uint8_t iface){
266  double dh_2;
267 
268  int A[4][2] = { {0,1} , {2,1} , {1,0} , {1,2} };
269 
270  dh_2 = double(getSize())/2.0;
271  vector<double> center(3);
272 
273  if (iface < global2D.nfaces){
274  center[0] = (double)x + (double)A[iface][0] * dh_2;
275  center[1] = (double)y + (double)A[iface][1] * dh_2;
276  center[2] = 0.0;
277  }
278  return center;
279  };
280 
281  // ------------------------------------------------------------------------------- //
282 
286  void getNodes(u32vector2D & nodes){
287  uint8_t i, cx, cy;
288  uint32_t dh;
289 
290  dh = getSize();
291  nodes.clear();
292  nodes.resize(global2D.nnodes);
293 
294  for (i = 0; i < global2D.nnodes; i++){
295  nodes[i].resize(3);
296  cx = uint8_t(i%2);
297  cy = uint8_t(i/2);
298  nodes[i][0] = x + cx*dh;
299  nodes[i][1] = y + cy*dh;
300  nodes[i][2] = 0;
301 #if defined(__INTEL_COMPILER) || defined(__ICC)
302 #else
303  nodes[i].shrink_to_fit();
304 #endif
305  }
306 #if defined(__INTEL_COMPILER) || defined(__ICC)
307 #else
308  nodes.shrink_to_fit();
309 #endif
310  };
311 
316  void getNode(u32vector & node, uint8_t inode){
317  uint8_t cx, cy;
318  uint32_t dh;
319 
320  dh = getSize();
321  node.clear();
322  node.resize(3);
323  cx = uint8_t(inode%2);
324  cy = uint8_t(inode/2);
325  node[0] = x + cx*dh;
326  node[1] = y + cy*dh;
327  };
328 
333  u32vector getNode(uint8_t inode){
334  u32vector node(3);
335  uint8_t cx, cy;
336  uint32_t dh;
337 
338  dh = getSize();
339  cx = uint8_t(inode%2);
340  cy = uint8_t(inode/2);
341  node[0] = x + cx*dh;
342  node[1] = y + cy*dh;
343  return node;
344  };
345 
346 
347  // ------------------------------------------------------------------------------- //
348 
353  void getNormal(uint8_t & iface,
354  vector<int8_t> & normal){
355  uint8_t i;
356 
357  normal.clear();
358  normal.resize(3);
359  for (i = 0; i < 3; i++){
360  normal[i] = global2D.normals[iface][i];
361  }
362 #if defined(__INTEL_COMPILER) || defined(__ICC)
363 #else
364  normal.shrink_to_fit();
365 #endif
366  };
367 
368  // ------------------------------------------------------------------------------- //
369 
373  uint64_t computeMorton() const{
374  uint64_t morton = 0;
375  morton = mortonEncode_magicbits(this->x,this->y);
376  return morton;
377  };
378 
379  // ------------------------------------------------------------------------------- //
380 
384  uint64_t computeMorton(){
385  uint64_t morton = 0;
386  morton = mortonEncode_magicbits(this->x,this->y);
387  return morton;
388  };
389 
390  //-------------------------------------------------------------------------------- //
391  // Other methods ----------------------------------------------------------------- //
392 
393 private:
394  Class_Octant<2> buildLastDesc(){ // Build last descendant of octant and return the last descendant octant (no info update)
395  uint32_t delta = (uint32_t)pow(2.0,(double)((uint8_t)MAX_LEVEL_2D - level)) - 1;
396  Class_Octant<2> last_desc(MAX_LEVEL_2D,x+delta,y+delta);
397  return last_desc;
398  };
399 
400  // ------------------------------------------------------------------------------- //
401 
402  Class_Octant<2> buildFather(){ // Build father of octant and return the father octant (no info update)
403  uint32_t deltax = x%(uint32_t(pow(2.0,(double)((uint8_t)MAX_LEVEL_2D - (level-1)))));
404  uint32_t deltay = y%(uint32_t(pow(2.0,(double)((uint8_t)MAX_LEVEL_2D - (level-1)))));
405  Class_Octant<2> father(level-1, (x-deltax), (y-deltay));
406  return father;
407  };
408 
409  // ------------------------------------------------------------------------------- //
410 
414  vector<Class_Octant<2> > buildChildren(){
415  uint8_t xf,yf;
416 
417  if (this->level < MAX_LEVEL_2D){
418  vector<Class_Octant<2> > children(global2D.nchildren);
419  for (int i=0; i<global2D.nchildren; i++){
420  switch (i) {
421  case 0 :
422  {
423  Class_Octant<2> oct(*this);
424  oct.setMarker(max(0,oct.marker-1));
425  oct.setLevel(oct.level+1);
426  oct.info[8]=true;
427  // Update interior face bound and pbound
428  xf=1; yf=3;
429  oct.info[xf] = oct.info[xf+global2D.nfaces] = false;
430  oct.info[yf] = oct.info[yf+global2D.nfaces] = false;
431  children[0] = oct;
432  }
433  break;
434  case 1 :
435  {
436  Class_Octant<2> oct(*this);
437  oct.setMarker(max(0,oct.marker-1));
438  oct.setLevel(oct.level+1);
439  oct.info[8]=true;
440  uint32_t dh = oct.getSize();
441  oct.x += dh;
442  // Update interior face bound and pbound
443  xf=0; yf=3;
444  oct.info[xf] = oct.info[xf+global2D.nfaces] = false;
445  oct.info[yf] = oct.info[yf+global2D.nfaces] = false;
446  children[1] = oct;
447  }
448  break;
449  case 2 :
450  {
451  Class_Octant<2> oct(*this);
452  oct.setMarker(max(0,oct.marker-1));
453  oct.setLevel(oct.level+1);
454  oct.info[8]=true;
455  uint32_t dh = oct.getSize();
456  oct.y += dh;
457  // Update interior face bound and pbound
458  xf=1; yf=2;
459  oct.info[xf] = oct.info[xf+global2D.nfaces] = false;
460  oct.info[yf] = oct.info[yf+global2D.nfaces] = false;
461  children[2] = oct;
462  }
463  break;
464  case 3 :
465  {
466  Class_Octant<2> oct(*this);
467  oct.setMarker(max(0,oct.marker-1));
468  oct.setLevel(oct.level+1);
469  oct.info[8]=true;
470  uint32_t dh = oct.getSize();
471  oct.x += dh;
472  oct.y += dh;
473  // Update interior face bound and pbound
474  xf=0; yf=2;
475  oct.info[xf] = oct.info[xf+global2D.nfaces] = false;
476  oct.info[yf] = oct.info[yf+global2D.nfaces] = false;
477  children[3] = oct;
478  }
479  break;
480  }
481  }
482  return children;
483  }
484  else{
485  vector<Class_Octant<2> > children(0);
486  //writeLog("Max level reached ---> No Children Built");
487  return children;
488  }
489  };
490 
491  // ------------------------------------------------------------------------------- //
492 
493  vector<uint64_t > computeHalfSizeMorton(uint8_t iface, // Computes Morton index (without level) of "n=sizehf" half-size (or same size if level=maxlevel)
494  uint32_t & sizehf){ // possible neighbours of octant throught face iface (sizehf=0 if boundary octant)
495  uint32_t dh,dh2;
496  uint32_t nneigh;
497  uint32_t i,cx,cy;
498 
499  nneigh = (level < MAX_LEVEL_2D) ? global2D.nchildren/2 : 1;
500  dh = (level < MAX_LEVEL_2D) ? getSize()/2 : getSize();
501  dh2 = getSize();
502 
503  if (info[iface]){
504  sizehf = 0;
505  vector<uint64_t> Morton(0);
506  return Morton;
507  }
508  else{
509  vector<uint64_t> Morton(nneigh);
510  switch (iface) {
511  case 0 :
512  {
513  for (i=0; i<nneigh; i++){
514  cy = (i==1);
515  Morton[i] = mortonEncode_magicbits(this->x-dh,this->y+dh*cy);
516  }
517  }
518  break;
519  case 1 :
520  {
521  for (i=0; i<nneigh; i++){
522  cy = (i==1);
523  Morton[i] = mortonEncode_magicbits(this->x+dh2,this->y+dh*cy);
524  }
525  }
526  break;
527  case 2 :
528  {
529  for (i=0; i<nneigh; i++){
530  cx = (i==1);
531  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y-dh);
532  }
533  }
534  break;
535  case 3 :
536  {
537  for (i=0; i<nneigh; i++){
538  cx = (i==1);
539  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2);
540  }
541  }
542  break;
543  }
544  sizehf = nneigh;
545  return Morton;
546  }
547 
548 
549  };
550 
551  // ------------------------------------------------------------------------------- //
552 
553  vector<uint64_t> computeMinSizeMorton(uint8_t iface, // Computes Morton index (without level) of "n=sizem" min-size (or same size if level=maxlevel)
554  const uint8_t & maxdepth, // possible neighbours of octant throught face iface (sizem=0 if boundary octant)
555  uint32_t & sizem){
556  uint32_t dh,dh2;
557  uint32_t nneigh, nline;
558  uint32_t i,cx,cy;
559 
560  nneigh = (level < MAX_LEVEL_2D) ? uint32_t(pow(2.0,double(maxdepth-level))) : 1;
561  dh = (level < MAX_LEVEL_2D) ? uint32_t(pow(2.0,double(MAX_LEVEL_2D - maxdepth))) : getSize();
562  dh2 = getSize();
563  nline = uint32_t(pow(2.0,double((maxdepth-level))));
564 
565  if (info[iface]){
566  sizem = 0;
567  vector<uint64_t> Morton(0);
568  return Morton;
569  }
570  else{
571  vector<uint64_t> Morton(nneigh);
572  switch (iface) {
573  case 0 :
574  {
575  for (i=0; i<nneigh; i++){
576  cy = (i%nline);
577  Morton[i] = mortonEncode_magicbits(this->x-dh,this->y+dh*cy);
578  }
579  }
580  break;
581  case 1 :
582  {
583  for (i=0; i<nneigh; i++){
584  cy = (i%nline);
585  Morton[i] = mortonEncode_magicbits(this->x+dh2,this->y+dh*cy);
586  }
587  }
588  break;
589  case 2 :
590  {
591  for (i=0; i<nneigh; i++){
592  cx = (i%nline);
593  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y-dh);
594  }
595  }
596  break;
597  case 3 :
598  {
599  for (i=0; i<nneigh; i++){
600  cx = (i%nline);
601  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2);
602  }
603  }
604  break;
605  }
606  sizem = nneigh;
607  sort(Morton.begin(), Morton.end());
608  return Morton;
609  }
610 
611  };
612 
613  // ------------------------------------------------------------------------------- //
614 
615  vector<uint64_t> computeVirtualMorton(uint8_t iface, // Computes Morton index (without level) of possible (virtual) neighbours of octant throught iface
616  const uint8_t & maxdepth, // Checks if balanced or not and uses half-size or min-size method (sizeneigh=0 if boundary octant)
617  uint32_t & sizeneigh){
618  if (getNotBalance()){
619  return computeMinSizeMorton(iface,
620  maxdepth,
621  sizeneigh);
622  }
623  else{
624  return computeHalfSizeMorton(iface,
625  sizeneigh);
626  }
627  };
628 
629  // ------------------------------------------------------------------------------- //
630 
631  uint64_t computeNodeHalfSizeMorton(uint8_t inode, // Computes Morton index (without level) of "n=sizehf" half-size (or same size if level=maxlevel)
632  uint32_t & sizehf){ // possible neighbours of octant throught face iface (sizehf=0 if boundary octant)
633  uint32_t dh,dh2;
634  uint32_t nneigh;
635  int8_t cx,cy;
636  uint8_t iface1, iface2;
637  nneigh = 1;
638  dh = (level < MAX_LEVEL_2D) ? getSize()/2 : getSize();
639  dh2 = getSize();
640  iface1 = global3D.nodeface[inode][0];
641  iface2 = global3D.nodeface[inode][1];
642 
643  if (info[iface1] || info[iface2]){
644  sizehf = 0;
645  return this->computeMorton();
646  }
647  else{
648  uint64_t Morton;
649  switch (inode) {
650  case 0 :
651  {
652  cx = -1;
653  cy = -1;
654  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy);
655  }
656  break;
657  case 1 :
658  {
659  cx = 1;
660  cy = -1;
661  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy);
662  }
663  break;
664  case 2 :
665  {
666  cx = -1;
667  cy = 1;
668  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy);
669  }
670  break;
671  case 3 :
672  {
673  cx = 1;
674  cy = 1;
675  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy);
676  }
677  break;
678  }
679  sizehf = nneigh;
680  return Morton;
681  }
682  };
683 
684  // ------------------------------------------------------------------------------- //
685 
686  uint64_t computeNodeMinSizeMorton(uint8_t inode, // Computes Morton index (without level) of "n=sizem" min-size (or same size if level=maxlevel)
687  const uint8_t & maxdepth, // possible neighbours of octant throught face iface (sizem=0 if boundary octant)
688  uint32_t & sizehf){
689  uint32_t dh,dh2;
690  uint32_t nneigh;
691  int8_t cx,cy;
692  uint8_t iface1, iface2;
693 
694  nneigh = 1;
695  dh = (level < MAX_LEVEL_2D) ? uint32_t(pow(2.0,double(MAX_LEVEL_2D - maxdepth))) : getSize();
696  dh2 = getSize();
697  iface1 = global3D.nodeface[inode][0];
698  iface2 = global3D.nodeface[inode][1];
699 
700  if (info[iface1] || info[iface2]){
701  sizehf = 0;
702  return this->computeMorton();
703  }
704  else{
705  uint64_t Morton;
706  switch (inode) {
707  case 0 :
708  {
709  cx = -1;
710  cy = -1;
711  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy);
712  }
713  break;
714  case 1 :
715  {
716  cx = 1;
717  cy = -1;
718  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy);
719  }
720  break;
721  case 2 :
722  {
723  cx = -1;
724  cy = 1;
725  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy);
726  }
727  break;
728  case 3 :
729  {
730  cx = 1;
731  cy = 1;
732  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy);
733  }
734  break;
735  }
736  sizehf = nneigh;
737  return Morton;
738  }
739 
740  };
741 
742  // ------------------------------------------------------------------------------- //
743 
744  uint64_t computeNodeVirtualMorton(uint8_t inode, // Computes Morton index (without level) of possible (virtual) neighbours of octant throught iface
745  const uint8_t & maxdepth, // Checks if balanced or not and uses half-size or min-size method (sizeneigh=0 if boundary octant)
746  uint32_t & sizeneigh){
747 
748  return computeNodeMinSizeMorton(inode,
749  maxdepth,
750  sizeneigh);
751 
752  };
753 
754  // ------------------------------------------------------------------------------- //
755 
756 };
void setBalance(bool balance)
uint32_t getX() const
int8_t getMarker() const
Parallel Octree Manager Class.
u32vector getNode(uint8_t inode)
int8_t normals[4][3]
Octant class definition - 2D specialization.
uint64_t computeMorton() const
uint32_t getSize() const
uint32_t getArea() const
dvector getFaceCenter(uint8_t iface)
void setMarker(int8_t marker)
void getNodes(u32vector2D &nodes)
uint64_t getVolume() const
bool getBound(uint8_t face) const
bool getIsNewC() const
void getNormal(uint8_t &iface, vector< int8_t > &normal)
bool getPbound(uint8_t face) const
bool getBalance() const
uint32_t getY() const
uint8_t getLevel() const
uint8_t nodeface[8][3]
uint32_t getZ() const
Local octree portion for each process.
bool getIsNewR() const
void getNode(u32vector &node, uint8_t inode)
bool getNotBalance() const
Octant class definition.