PABLO  0.1
PArallel Balanced Linear Octree
 All Classes Functions Variables Pages
Class_Octant_3D.tpp
1 
28 // =================================================================================== //
29 // CLASS SPECIALIZATION //
30 // =================================================================================== //
31 template<>
32 class Class_Octant<3>{
33  // ------------------------------------------------------------------------------- //
34  // FRIENDSHIPS ------------------------------------------------------------------- //
35  template<int dim> friend class Class_Local_Tree;
36  template<int dim> friend class Class_Para_Tree;
37 
38  // ------------------------------------------------------------------------------- //
39  // TYPEDEFS ----------------------------------------------------------------------- //
40 public:
41  typedef vector<Class_Octant<3> > OctantsType;
42  typedef vector<double> dvector;
43  typedef vector<uint32_t> u32vector;
44  typedef vector<vector<uint32_t> > u32vector2D;
45  typedef vector<vector<uint64_t> > u64vector2D;
46 
47  // ------------------------------------------------------------------------------- //
48  // MEMBERS ----------------------------------------------------------------------- //
49 
50 private:
51  uint32_t x;
52  uint32_t y;
53  uint32_t z;
54  uint8_t level;
55  int8_t marker;
56  bitset<16> info;
63  // ------------------------------------------------------------------------------- //
64  // CONSTRUCTORS AND OPERATORS----------------------------------------------------- //
65 
66 public:
67  Class_Octant(){
68  x = y = z = 0;
69  level = 0;
70  marker = 0;
71  for (int i=0; i<global3D.nfaces; i++){
72  info[i] = true;
73  }
74  };
75 
76  Class_Octant(int8_t level, int32_t x, int32_t y, int32_t z){
77  this->x = x;
78  this->y = y;
79  this->z = z;
80  this->level = level;
81  marker = 0;
82  if (level==0){
83  for (int i=0; i<global3D.nfaces; i++){
84  info[i] = true;
85  }
86  }
87 
88  };
89 
90  Class_Octant(int8_t level, int32_t x, int32_t y, int32_t z, bool bound){
91  this->x = x;
92  this->y = y;
93  this->z = z;
94  this->level = level;
95  marker = 0;
96  if (level==0){
97  for (int i=0; i<global3D.nfaces; i++){
98  info[i] = bound;
99  }
100  }
101 
102  };
103 
104 
105  Class_Octant(const Class_Octant<3> &octant){
106  x = octant.x;
107  y = octant.y;
108  z = octant.z;
109  level = octant.level;
110  marker = octant.marker;
111  info = octant.info;
112  };
113 
116  bool operator ==(const Class_Octant<3> & oct2){
117  bool check = true;
118  check = check && (x == oct2.x);
119  check = check && (y == oct2.y);
120  check = check && (z == oct2.z);
121  check = check && (level == oct2.level);
122  return check;
123  }
124 
125  // ------------------------------------------------------------------------------- //
126  // METHODS ----------------------------------------------------------------------- //
127 
128  // Basic Get/Set methods --------------------------------------------------------- //
129 
130 public:
134  uint32_t getX() const{return x;};
135 
139  uint32_t getY() const{return y;};
140 
144  uint32_t getZ() const{return z;};
145 
149  uint8_t getLevel() const{return level;};
150 
154  int8_t getMarker() const{return marker;};
155 
160  bool getBound(uint8_t face) const{
161  return info[face];
162  };
163 
164 private:
165  void setBound(uint8_t face) {
166  info[face] = true;
167  };
168 
169 public:
174  bool getPbound(uint8_t face) const{
175  return info[global3D.nfaces+face];
176  };
177 
181  bool getIsNewR() const{return info[12];};
182 
186  bool getIsNewC() const{return info[13];};
187 
191  bool getNotBalance() const{return info[14];};
192 
196  bool getBalance() const{return (!info[10]);};
197 
201  void setMarker(int8_t marker){
202  this->marker = marker;
203  };
204 
208  void setBalance(bool balance){
209  info[14] = balance;
210  };
211 
212 private:
213  void setLevel(uint8_t level){
214  this->level = level;
215  };
216 
217  void setPbound(uint8_t face, bool flag){
218  info[global3D.nfaces+face] = flag;
219  };
220 
221  //-------------------------------------------------------------------------------- //
222  // Other Get/Set methods --------------------------------------------------------- //
223 
224 public:
228  uint32_t getSize() const{
229  uint32_t size = uint32_t(pow(double(2),double(MAX_LEVEL_3D-level)));
230  return size;
231  };
232 
236  uint64_t getArea() const{
237  uint64_t area = uint64_t(pow(double(getSize()),2.0));
238  return area;
239  };
240 
244  uint64_t getVolume() const{
245  uint64_t volume = uint64_t(pow(double(getSize()),3.0));
246  return volume;
247  };
248 
249  // ------------------------------------------------------------------------------- //
250 
254  dvector getCenter(){
255  double dh;
256 
257  dh = double(getSize())/2.0;
258  vector<double> center(3);
259 
260  center[0] = (double)x + dh;
261  center[1] = (double)y + dh;
262  center[2] = (double)z + dh;
263  return center;
264  };
265 
266  // ------------------------------------------------------------------------------- //
267 
271  dvector getFaceCenter(uint8_t iface){
272  double dh_2;
273 
274  int A[6][3] = { {0,1,1} , {2,1,1} , {1,0,1} , {1,2,1} , {1,1,0} , {1,1,2} };
275 
276  dh_2 = double(getSize())/2.0;
277  vector<double> center(3);
278 
279  if (iface < global3D.nfaces){
280  center[0] = (double)x + (double)A[iface][0] * dh_2;
281  center[1] = (double)y + (double)A[iface][1] * dh_2;
282  center[2] = (double)z + (double)A[iface][2] * dh_2;
283  }
284  return center;
285  };
286 
287  // ------------------------------------------------------------------------------- //
288 
292  dvector getEdgeCenter(uint8_t iedge){
293  double dh_2;
294 
295  int A[12][3] = { {0,1,0},{2,1,0},{1,0,0},{1,2,0},{0,0,1},{2,0,1},{0,2,1},{2,2,1},{0,1,2},{2,1,2},{1,0,2},{1,2,2} };//{ {0,1,1} , {2,1,1} , {1,0,1} , {1,2,1} , {1,1,0} , {1,1,2} };
296 
297  dh_2 = double(getSize())/2.0;
298  vector<double> center(3);
299 
300  if (iedge < global3D.nedges){
301  center[0] = (double)x + (double)A[iedge][0] * dh_2;
302  center[1] = (double)y + (double)A[iedge][1] * dh_2;
303  center[2] = (double)z + (double)A[iedge][2] * dh_2;
304  }
305  return center;
306  };
307 
308  // ------------------------------------------------------------------------------- //
309 
313  void getNodes(u32vector2D & nodes){
314  uint8_t i, cx, cy, cz;
315  uint32_t dh;
316 
317  dh = getSize();
318  nodes.clear();
319  nodes.resize(global3D.nnodes);
320 
321  for (i = 0; i < global3D.nnodes; i++){
322  nodes[i].resize(3);
323  cx = uint8_t(i%2);
324  cy = uint8_t((i-4*(i/4))/2);
325  cz = uint8_t(i/4);
326  nodes[i][0] = x + cx*dh;
327  nodes[i][1] = y + cy*dh;
328  nodes[i][2] = z + cz*dh;
329 #if defined(__INTEL_COMPILER) || defined(__ICC)
330 #else
331  nodes[i].shrink_to_fit();
332 #endif
333  }
334 #if defined(__INTEL_COMPILER) || defined(__ICC)
335 #else
336  nodes.shrink_to_fit();
337 #endif
338  };
339 
344  void getNode(u32vector & node, uint8_t inode){
345  uint8_t cx, cy, cz;
346  uint32_t dh;
347 
348  dh = getSize();
349  node.clear();
350  node.resize(3);
351  cx = inode%2;
352  cy = (inode-4*(inode/4))/2;
353  cz = inode/4;
354  node[0] = x + cx*dh;
355  node[1] = y + cy*dh;
356  node[2] = z + cz*dh;
357 
358  };
359 
364  u32vector getNode(uint8_t inode){
365  u32vector node;
366  uint8_t cx, cy, cz;
367  uint32_t dh;
368 
369  dh = getSize();
370  node.clear();
371  node.resize(3);
372  cx = inode%2;
373  cy = (inode-4*(inode/4))/2;
374  cz = inode/4;
375  node[0] = x + cx*dh;
376  node[1] = y + cy*dh;
377  node[2] = z + cz*dh;
378  return node;
379  };
380 
381  // ------------------------------------------------------------------------------- //
382 
387  void getNormal(uint8_t & iface,
388  vector<int8_t> & normal){
389  uint8_t i;
390 
391  normal.clear();
392  normal.resize(3);
393  for (i = 0; i < 3; i++){
394  normal[i] = global3D.normals[iface][i];
395  }
396 #if defined(__INTEL_COMPILER) || defined(__ICC)
397 #else
398  normal.shrink_to_fit();
399 #endif
400  };
401 
402  // ------------------------------------------------------------------------------- //
403 
407  uint64_t computeMorton() const{ uint64_t morton = 0;
408  morton = mortonEncode_magicbits(this->x,this->y,this->z);
409  return morton;
410  };
411 
412  // ------------------------------------------------------------------------------- //
413 
417  uint64_t computeMorton(){
418  uint64_t morton = 0;
419  morton = mortonEncode_magicbits(this->x,this->y,this->z);
420  return morton;
421  };
422 
423  // =================================================================================== //
424  // Other methods //
425  // =================================================================================== //
426 
427 private:
428  Class_Octant<3> buildLastDesc(){ // Build last descendant of octant and return the last descendant octant (no info update)
429  uint32_t delta = (uint32_t)pow(2.0,(double)((uint8_t)MAX_LEVEL_3D - level)) - 1;
430  Class_Octant<3> last_desc(MAX_LEVEL_3D,x+delta,y+delta,z+delta);
431  return last_desc;
432  };
433 
434  // =================================================================================== //
435 
436  Class_Octant<3> buildFather(){ // Build father of octant and return the father octant (no info update)
437  uint32_t deltax = x%(uint32_t(pow(2.0,(double)((uint8_t)MAX_LEVEL_3D - (level-1)))));
438  uint32_t deltay = y%(uint32_t(pow(2.0,(double)((uint8_t)MAX_LEVEL_3D - (level-1)))));
439  uint32_t deltaz = z%(uint32_t(pow(2.0,(double)((uint8_t)MAX_LEVEL_3D - (level-1)))));
440  Class_Octant<3> father(level-1, x-deltax, y-deltay, z-deltaz);
441  return father;
442  };
443 
444  // =================================================================================== //
445 
446  // ------------------------------------------------------------------------------- //
447 
451  vector< Class_Octant<3> > buildChildren(){
452  uint8_t xf,yf,zf;
453 
454  if (this->level < MAX_LEVEL_3D){
455  vector< Class_Octant<3> > children(global3D.nchildren);
456  for (int i=0; i<global3D.nchildren; i++){
457  switch (i) {
458  case 0 :
459  {
460  Class_Octant<3> oct(*this);
461  oct.setMarker(max(0,oct.marker-1));
462  oct.setLevel(oct.level+1);
463  oct.info[12]=true;
464  // Update interior face bound and pbound
465  xf=1; yf=3; zf=5;
466  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
467  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
468  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
469  children[0] = oct;
470  }
471  break;
472  case 1 :
473  {
474  Class_Octant<3> oct(*this);
475  oct.setMarker(max(0,oct.marker-1));
476  oct.setLevel(oct.level+1);
477  oct.info[12]=true;
478  uint32_t dh = oct.getSize();
479  oct.x += dh;
480  // Update interior face bound and pbound
481  xf=0; yf=3; zf=5;
482  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
483  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
484  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
485  children[1] = oct;
486  }
487  break;
488  case 2 :
489  {
490  Class_Octant<3> oct(*this);
491  oct.setMarker(max(0,oct.marker-1));
492  oct.setLevel(oct.level+1);
493  oct.info[12]=true;
494  uint32_t dh = oct.getSize();
495  oct.y += dh;
496  // Update interior face bound and pbound
497  xf=1; yf=2; zf=5;
498  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
499  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
500  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
501  children[2] = oct;
502  }
503  break;
504  case 3 :
505  {
506  Class_Octant<3> oct(*this);
507  oct.setMarker(max(0,oct.marker-1));
508  oct.setLevel(oct.level+1);
509  oct.info[12]=true;
510  uint32_t dh = oct.getSize();
511  oct.x += dh;
512  oct.y += dh;
513  // Update interior face bound and pbound
514  xf=0; yf=2; zf=5;
515  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
516  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
517  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
518  children[3] = oct;
519  }
520  break;
521  case 4 :
522  {
523  Class_Octant<3> oct(*this);
524  oct.setMarker(max(0,oct.marker-1));
525  oct.setLevel(oct.level+1);
526  oct.info[12]=true;
527  uint32_t dh = oct.getSize();
528  oct.z += dh;
529  // Update interior face bound and pbound
530  xf=1; yf=3; zf=4;
531  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
532  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
533  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
534  children[4] = oct;
535  }
536  break;
537  case 5 :
538  {
539  Class_Octant<3> oct(*this);
540  oct.setMarker(max(0,oct.marker-1));
541  oct.setLevel(oct.level+1);
542  oct.info[12]=true;
543  uint32_t dh = oct.getSize();
544  oct.x += dh;
545  oct.z += dh;
546  // Update interior face bound and pbound
547  xf=0; yf=3; zf=4;
548  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
549  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
550  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
551  children[5] = oct;
552  }
553  break;
554  case 6 :
555  {
556  Class_Octant<3> oct(*this);
557  oct.setMarker(max(0,oct.marker-1));
558  oct.setLevel(oct.level+1);
559  oct.info[12]=true;
560  uint32_t dh = oct.getSize();
561  oct.y += dh;
562  oct.z += dh;
563  // Update interior face bound and pbound
564  xf=1; yf=2; zf=4;
565  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
566  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
567  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
568  children[6] = oct;
569  }
570  break;
571  case 7 :
572  {
573  Class_Octant<3> oct(*this);
574  oct.setMarker(max(0,oct.marker-1));
575  oct.setLevel(oct.level+1);
576  oct.info[12]=true;
577  uint32_t dh = oct.getSize();
578  oct.x += dh;
579  oct.y += dh;
580  oct.z += dh;
581  // Update interior face bound and pbound
582  xf=0; yf=2; zf=4;
583  oct.info[xf] = oct.info[xf+global3D.nfaces] = false;
584  oct.info[yf] = oct.info[yf+global3D.nfaces] = false;
585  oct.info[zf] = oct.info[zf+global3D.nfaces] = false;
586  children[7] = oct;
587  }
588  break;
589  }
590  }
591  return children;
592  }
593  else{
594  vector< Class_Octant<3> > children(0);
595  //writeLog("Max level reached ---> No Children Built");
596  return children;
597  }
598  };
599 
600  // ------------------------------------------------------------------------------- //
601 
602  vector<uint64_t> computeHalfSizeMorton(uint8_t iface, // Computes Morton index (without level) of "n=sizehf" half-size (or same size if level=maxlevel)
603  uint32_t & sizehf){ // possible neighbours of octant throught face iface (sizehf=0 if boundary octant)
604  uint32_t dh,dh2;
605  uint32_t nneigh;
606  uint32_t i,cx,cy,cz;
607 
608  nneigh = (level < MAX_LEVEL_3D) ? global3D.nchildren/2 : 1;
609  dh = (level < MAX_LEVEL_3D) ? getSize()/2 : getSize();
610  dh2 = getSize();
611 
612  if (info[iface]){
613  sizehf = 0;
614  vector<uint64_t> Morton(0);
615  return Morton;
616  }
617  else{
618  vector<uint64_t> Morton(nneigh);
619  switch (iface) {
620  case 0 :
621  {
622  for (i=0; i<nneigh; i++){
623  cy = (i==1)||(i==3);
624  cz = (i==2)||(i==3);
625  Morton[i] = mortonEncode_magicbits(this->x-dh,this->y+dh*cy,this->z+dh*cz);
626  }
627  }
628  break;
629  case 1 :
630  {
631  for (i=0; i<nneigh; i++){
632  cy = (i==1)||(i==3);
633  cz = (i==2)||(i==3);
634  Morton[i] = mortonEncode_magicbits(this->x+dh2,this->y+dh*cy,this->z+dh*cz);
635  }
636  }
637  break;
638  case 2 :
639  {
640  for (i=0; i<nneigh; i++){
641  cx = (i==1)||(i==3);
642  cz = (i==2)||(i==3);
643  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y-dh,this->z+dh*cz);
644  }
645  }
646  break;
647  case 3 :
648  {
649  for (i=0; i<nneigh; i++){
650  cx = (i==1)||(i==3);
651  cz = (i==2)||(i==3);
652  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2,this->z+dh*cz);
653  }
654  }
655  break;
656  case 4 :
657  {
658  for (i=0; i<nneigh; i++){
659  cx = (i==1)||(i==3);
660  cy = (i==2)||(i==3);
661  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z-dh);
662  }
663  }
664  break;
665  case 5 :
666  {
667  for (i=0; i<nneigh; i++){
668  cx = (i==1)||(i==3);
669  cy = (i==2)||(i==3);
670  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2);
671  }
672  }
673  break;
674  }
675  sizehf = nneigh;
676  return Morton;
677  }
678 
679 
680  };
681 
682  // ------------------------------------------------------------------------------- //
683 
684  vector<uint64_t> computeMinSizeMorton(uint8_t iface, // Computes Morton index (without level) of "n=sizem" min-size (or same size if level=maxlevel)
685  const uint8_t & maxdepth, // possible neighbours of octant throught face iface (sizem=0 if boundary octant)
686  uint32_t & sizem){
687  uint32_t dh,dh2;
688  uint32_t nneigh, nline;
689  uint32_t i,cx,cy,cz;
690 
691  nneigh = (level < MAX_LEVEL_3D) ? uint32_t(pow(2.0,double(2*(maxdepth-level)))) : 1;
692  dh = (level < MAX_LEVEL_3D) ? uint32_t(pow(2.0,double(MAX_LEVEL_3D - maxdepth))) : getSize();
693  dh2 = getSize();
694  nline = uint32_t(pow(2.0,double((maxdepth-level))));
695 
696  if (info[iface]){
697  // uint64_t* Morton = new uint64_t[0];
698  sizem = 0;
699  // return Morton;
700  vector<uint64_t> Morton(0);
701  return Morton;
702  }
703  else{
704  vector<uint64_t> Morton(nneigh);
705  switch (iface) {
706  case 0 :
707  {
708  for (i=0; i<nneigh; i++){
709  cy = (i/nline);
710  cz = (i%nline);
711  Morton[i] = mortonEncode_magicbits(this->x-dh,this->y+dh*cy,this->z+dh*cz);
712  }
713  }
714  break;
715  case 1 :
716  {
717  for (i=0; i<nneigh; i++){
718  cy = (i/nline);
719  cz = (i%nline);
720  Morton[i] = mortonEncode_magicbits(this->x+dh2,this->y+dh*cy,this->z+dh*cz);
721  }
722  }
723  break;
724  case 2 :
725  {
726  for (i=0; i<nneigh; i++){
727  cx = (i/nline);
728  cz = (i%nline);
729  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y-dh,this->z+dh*cz);
730  }
731  }
732  break;
733  case 3 :
734  {
735  for (i=0; i<nneigh; i++){
736  cx = (i/nline);
737  cz = (i%nline);
738  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2,this->z+dh*cz);
739  }
740  }
741  break;
742  case 4 :
743  {
744  for (i=0; i<nneigh; i++){
745  cx = (i/nline);
746  cy = (i%nline);
747  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z-dh);
748  }
749  }
750  break;
751  case 5 :
752  {
753  for (i=0; i<nneigh; i++){
754  cx = (i/nline);
755  cy = (i%nline);
756  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2);
757  }
758  }
759  break;
760  }
761  sizem = nneigh;
762  // sort(Morton,Morton+nneigh);
763  sort(Morton.begin(), Morton.end());
764  return Morton;
765  }
766 
767  };
768 
769  // ------------------------------------------------------------------------------- //
770 
771  vector<uint64_t> computeVirtualMorton(uint8_t iface, // Computes Morton index (without level) of possible (virtual) neighbours of octant throught iface
772  const uint8_t & maxdepth, // Checks if balanced or not and uses half-size or min-size method (sizeneigh=0 if boundary octant)
773  uint32_t & sizeneigh){
774  vector<uint64_t> Morton;
775  if (getNotBalance()){
776  return computeMinSizeMorton(iface,
777  maxdepth,
778  sizeneigh);
779  }
780  else{
781  return computeHalfSizeMorton(iface,
782  sizeneigh);
783  }
784  };
785 
786  // ------------------------------------------------------------------------------- //
787 
788  vector<uint64_t> computeEdgeHalfSizeMorton(uint8_t iedge, // Computes Morton index (without level) of "n=sizehf" half-size (or same size if level=maxlevel)
789  uint32_t & sizehf){ // possible neighbours of octant throught face iface (sizehf=0 if boundary octant)
790  uint32_t dh,dh2;
791  uint32_t nneigh;
792  uint32_t i,cx,cy,cz;
793  uint8_t iface1, iface2;
794 
795  nneigh = (level < MAX_LEVEL_3D) ? 2 : 1;
796  dh = (level < MAX_LEVEL_3D) ? getSize()/2 : getSize();
797  dh2 = getSize();
798  iface1 = global3D.edgeface[iedge][0];
799  iface2 = global3D.edgeface[iedge][1];
800 
801  if (info[iface1] || info[iface2]){
802  sizehf = 0;
803  vector<uint64_t> Morton(0);
804  return Morton;
805  }
806  else{
807  vector<uint64_t> Morton(nneigh);
808  switch (iedge) {
809  case 0 :
810  {
811  for (i=0; i<nneigh; i++){
812  cx = -1;
813  cy = (i==1);
814  cz = -1;
815  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
816  }
817  }
818  break;
819  case 1 :
820  {
821  for (i=0; i<nneigh; i++){
822  cx = 1;
823  cy = (i==1);
824  cz = -1;
825  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh*cz);
826  }
827  }
828  break;
829  case 2 :
830  {
831  for (i=0; i<nneigh; i++){
832  cx = (i==1);
833  cy = -1;
834  cz = -1;
835  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
836  }
837  }
838  break;
839  case 3 :
840  {
841  for (i=0; i<nneigh; i++){
842  cx = (i==1);
843  cy = 1;
844  cz = -1;
845  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh*cz);
846  }
847  }
848  break;
849  case 4 :
850  {
851  for (i=0; i<nneigh; i++){
852  cx = -1;
853  cy = -1;
854  cz = (i==1);
855  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
856  }
857  }
858  break;
859  case 5 :
860  {
861  for (i=0; i<nneigh; i++){
862  cx = 1;
863  cy = -1;
864  cz = (i==1);
865  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh*cz);
866  }
867  }
868  break;
869  case 6 :
870  {
871  for (i=0; i<nneigh; i++){
872  cx = -1;
873  cy = 1;
874  cz = (i==1);
875  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh*cz);
876  }
877  }
878  break;
879  case 7 :
880  {
881  for (i=0; i<nneigh; i++){
882  cx = 1;
883  cy = 1;
884  cz = (i==1);
885  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy,this->z+dh*cz);
886  }
887  }
888  break;
889  case 8 :
890  {
891  for (i=0; i<nneigh; i++){
892  cx = -1;
893  cy = (i==1);
894  cz = 1;
895  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2*cz);
896  }
897  }
898  break;
899  case 9 :
900  {
901  for (i=0; i<nneigh; i++){
902  cx = 1;
903  cy = (i==1);
904  cz = 1;
905  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh2*cz);
906  }
907  }
908  break;
909  case 10 :
910  {
911  for (i=0; i<nneigh; i++){
912  cx = (i==1);
913  cy = -1;
914  cz = 1;
915  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2*cz);
916  }
917  }
918  break;
919  case 11 :
920  {
921  for (i=0; i<nneigh; i++){
922  cx = (i==1);
923  cy = 1;
924  cz = 1;
925  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh2*cz);
926  }
927  }
928  break;
929  }
930  sizehf = nneigh;
931  return Morton;
932  }
933 
934 
935  };
936 
937  // ------------------------------------------------------------------------------- //
938 
939  vector<uint64_t> computeEdgeMinSizeMorton(uint8_t iedge, // Computes Morton index (without level) of "n=sizem" min-size (or same size if level=maxlevel)
940  const uint8_t & maxdepth, // possible neighbours of octant throught edge iedge (sizem=0 if boundary octant)
941  uint32_t & sizem){
942  uint32_t dh,dh2;
943  uint32_t nneigh, nline;
944  uint32_t i,cx,cy,cz;
945  uint8_t iface1, iface2;
946 
947 
948  nneigh = (level < MAX_LEVEL_3D) ? uint32_t(pow(2.0,double((maxdepth-level)))) : 1;
949  dh = (level < MAX_LEVEL_3D) ? uint32_t(pow(2.0,double(MAX_LEVEL_3D - maxdepth))) : getSize();
950  dh2 = getSize();
951  nline = uint32_t(pow(2.0,double((maxdepth-level))));
952  iface1 = global3D.edgeface[iedge][0];
953  iface2 = global3D.edgeface[iedge][1];
954 
955  if (info[iface1] || info[iface2]){
956  sizem = 0;
957  vector<uint64_t> Morton(0);
958  return Morton;
959  }
960  else{
961  vector<uint64_t> Morton(nneigh);
962  switch (iedge) {
963  case 0 :
964  {
965  for (i=0; i<nneigh; i++){
966  cx = -1;
967  cy = i;
968  cz = -1;
969  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
970  }
971  }
972  break;
973  case 1 :
974  {
975  for (i=0; i<nneigh; i++){
976  cx = 1;
977  cy = i;
978  cz = -1;
979  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh*cz);
980  }
981  }
982  break;
983  case 2 :
984  {
985  for (i=0; i<nneigh; i++){
986  cx = i;
987  cy = -1;
988  cz = -1;
989  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
990  }
991  }
992  break;
993  case 3 :
994  {
995  for (i=0; i<nneigh; i++){
996  cx = i;
997  cy = 1;
998  cz = -1;
999  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh*cz);
1000  }
1001  }
1002  break;
1003  case 4 :
1004  {
1005  for (i=0; i<nneigh; i++){
1006  cx = -1;
1007  cy = -1;
1008  cz = i;
1009  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
1010  }
1011  }
1012  break;
1013  case 5 :
1014  {
1015  for (i=0; i<nneigh; i++){
1016  cx = 1;
1017  cy = -1;
1018  cz = i;
1019  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh*cz);
1020  }
1021  }
1022  break;
1023  case 6 :
1024  {
1025  for (i=0; i<nneigh; i++){
1026  cx = -1;
1027  cy = 1;
1028  cz = i;
1029  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh*cz);
1030  }
1031  }
1032  break;
1033  case 7 :
1034  {
1035  for (i=0; i<nneigh; i++){
1036  cx = 1;
1037  cy = 1;
1038  cz = i;
1039  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy,this->z+dh*cz);
1040  }
1041  }
1042  break;
1043  case 8 :
1044  {
1045  for (i=0; i<nneigh; i++){
1046  cx = -1;
1047  cy = i;
1048  cz = 1;
1049  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2*cz);
1050  }
1051  }
1052  break;
1053  case 9 :
1054  {
1055  for (i=0; i<nneigh; i++){
1056  cx = 1;
1057  cy = i;
1058  cz = 1;
1059  Morton[i] = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh2*cz);
1060  }
1061  }
1062  break;
1063  case 10 :
1064  {
1065  for (i=0; i<nneigh; i++){
1066  cx = i;
1067  cy = -1;
1068  cz = 1;
1069  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2*cz);
1070  }
1071  }
1072  break;
1073  case 11 :
1074  {
1075  for (i=0; i<nneigh; i++){
1076  cx = i;
1077  cy = 1;
1078  cz = 1;
1079  Morton[i] = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh2*cz);
1080  }
1081  }
1082  break;
1083  }
1084  sizem = nneigh;
1085  //sort(Morton,Morton+nneigh);
1086  sort(Morton.begin(),Morton.end());
1087  return Morton;
1088  }
1089  };
1090 
1091  // ------------------------------------------------------------------------------- //
1092 
1093  vector<uint64_t> computeEdgeVirtualMorton(uint8_t iedge, // Computes Morton index (without level) of possible (virtual) neighbours of octant throught iface
1094  const uint8_t & maxdepth, // Checks if balanced or not and uses half-size or min-size method (sizeneigh=0 if boundary octant)
1095  uint32_t & sizeneigh,
1096  uint8_t balance_codim){
1097 // if (getNotBalance()){
1098  return computeEdgeMinSizeMorton(iedge,
1099  maxdepth,
1100  sizeneigh);
1101 // }
1102 // else{
1103 // return computeEdgeHalfSizeMorton(iedge,
1104 // sizeneigh);
1105 // }
1106  if(!getNotBalance() && balance_codim > 1){
1107  return computeEdgeHalfSizeMorton(iedge,
1108  sizeneigh);
1109  }
1110  else{
1111  return computeEdgeMinSizeMorton(iedge,
1112  maxdepth,
1113  sizeneigh);
1114  }
1115  };
1116 
1117  // ------------------------------------------------------------------------------- //
1118 
1119 
1120  uint64_t computeNodeHalfSizeMorton(uint8_t inode, // Computes Morton index (without level) of "n=sizehf" half-size (or same size if level=maxlevel)
1121  uint32_t & sizehf){ // possible neighbours of octant throught face iface (sizehf=0 if boundary octant)
1122  uint32_t dh,dh2;
1123  uint32_t nneigh;
1124  int8_t cx,cy,cz;
1125  uint8_t iface1, iface2, iface3;
1126  nneigh = 1;
1127  dh = (level < MAX_LEVEL_3D) ? getSize()/2 : getSize();
1128  dh2 = getSize();
1129  iface1 = global3D.nodeface[inode][0];
1130  iface2 = global3D.nodeface[inode][1];
1131  iface3 = global3D.nodeface[inode][2];
1132 
1133  if (info[iface1] || info[iface2] || info[iface3]){
1134  sizehf = 0;
1135  return this->computeMorton();
1136  }
1137  else{
1138  uint64_t Morton;
1139  switch (inode) {
1140  case 0 :
1141  {
1142  cx = -1;
1143  cy = -1;
1144  cz = -1;
1145  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
1146  }
1147  break;
1148  case 1 :
1149  {
1150  cx = 1;
1151  cy = -1;
1152  cz = -1;
1153  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh*cz);
1154  }
1155  break;
1156  case 2 :
1157  {
1158  cx = -1;
1159  cy = 1;
1160  cz = -1;
1161  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh*cz);
1162  }
1163  break;
1164  case 3 :
1165  {
1166  cx = 1;
1167  cy = 1;
1168  cz = -1;
1169  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy,this->z+dh*cz);
1170  }
1171  break;
1172  case 4 :
1173  {
1174  cx = -1;
1175  cy = -1;
1176  cz = 1;
1177  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2*cz);
1178  }
1179  break;
1180  case 5 :
1181  {
1182  cx = 1;
1183  cy = -1;
1184  cz = 1;
1185  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh2*cz);
1186  }
1187  break;
1188  case 6 :
1189  {
1190  cx = -1;
1191  cy = 1;
1192  cz = 1;
1193  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh2*cz);
1194  }
1195  break;
1196  case 7 :
1197  {
1198  cx = 1;
1199  cy = 1;
1200  cz = 1;
1201  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy,this->z+dh2*cz);
1202  }
1203  break;
1204  }
1205  sizehf = nneigh;
1206  return Morton;
1207  }
1208  };
1209 
1210  // ------------------------------------------------------------------------------- //
1211 
1212  uint64_t computeNodeMinSizeMorton(uint8_t inode, // Computes Morton index (without level) of "n=sizem" min-size (or same size if level=maxlevel)
1213  const uint8_t & maxdepth, // possible neighbours of octant throught face iface (sizem=0 if boundary octant)
1214  uint32_t & sizehf){
1215  uint32_t dh,dh2;
1216  uint32_t nneigh;
1217  int8_t cx,cy,cz;
1218  uint8_t iface1, iface2, iface3;
1219 
1220  nneigh = 1;
1221  dh = (level < MAX_LEVEL_3D) ? uint32_t(pow(2.0,double(MAX_LEVEL_3D - maxdepth))) : getSize();
1222  dh2 = getSize();
1223  iface1 = global3D.nodeface[inode][0];
1224  iface2 = global3D.nodeface[inode][1];
1225  iface3 = global3D.nodeface[inode][2];
1226 
1227  if (info[iface1] || info[iface2] || info[iface3]){
1228  sizehf = 0;
1229  return this->computeMorton();
1230  }
1231  else{
1232  uint64_t Morton;
1233  switch (inode) {
1234  case 0 :
1235  {
1236  cx = -1;
1237  cy = -1;
1238  cz = -1;
1239  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh*cz);
1240  }
1241  break;
1242  case 1 :
1243  {
1244  cx = 1;
1245  cy = -1;
1246  cz = -1;
1247  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh*cz);
1248  }
1249  break;
1250  case 2 :
1251  {
1252  cx = -1;
1253  cy = 1;
1254  cz = -1;
1255  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh*cz);
1256  }
1257  break;
1258  case 3 :
1259  {
1260  cx = 1;
1261  cy = 1;
1262  cz = -1;
1263  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy,this->z+dh*cz);
1264  }
1265  break;
1266  case 4 :
1267  {
1268  cx = -1;
1269  cy = -1;
1270  cz = 1;
1271  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh*cy,this->z+dh2*cz);
1272  }
1273  break;
1274  case 5 :
1275  {
1276  cx = 1;
1277  cy = -1;
1278  cz = 1;
1279  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh*cy,this->z+dh2*cz);
1280  }
1281  break;
1282  case 6 :
1283  {
1284  cx = -1;
1285  cy = 1;
1286  cz = 1;
1287  Morton = mortonEncode_magicbits(this->x+dh*cx,this->y+dh2*cy,this->z+dh2*cz);
1288  }
1289  break;
1290  case 7 :
1291  {
1292  cx = 1;
1293  cy = 1;
1294  cz = 1;
1295  Morton = mortonEncode_magicbits(this->x+dh2*cx,this->y+dh2*cy,this->z+dh2*cz);
1296  }
1297  break;
1298  }
1299  sizehf = nneigh;
1300  return Morton;
1301  }
1302 
1303  };
1304 
1305  // ------------------------------------------------------------------------------- //
1306 
1307  uint64_t computeNodeVirtualMorton(uint8_t inode, // Computes Morton index (without level) of possible (virtual) neighbours of octant throught iface
1308  const uint8_t & maxdepth, // Checks if balanced or not and uses half-size or min-size method (sizeneigh=0 if boundary octant)
1309  uint32_t & sizeneigh){
1310 // if (getNotBalance()){
1311  return computeNodeMinSizeMorton(inode,
1312  maxdepth,
1313  sizeneigh);
1314 // }
1315 // else{
1316 // return computeNodeHalfSizeMorton(inode,
1317 // sizeneigh);
1318 // }
1319  };
1320 
1321  // ------------------------------------------------------------------------------- //
1322 
1323 
1324 };//end Class_Octant<3>
1325 //#endif
1326 
1327 
Class_Octant< 3 >::computeMorton
uint64_t computeMorton() const
Definition: Class_Octant_3D.tpp:407
Class_Octant< 3 >::getEdgeCenter
dvector getEdgeCenter(uint8_t iedge)
Definition: Class_Octant_3D.tpp:292
Class_Global< 3 >::nchildren
uint8_t nchildren
Definition: Class_Global_3D.tpp:180
Class_Para_Tree
Parallel Octree Manager Class.
Definition: Class_Para_Tree.hpp:64
Class_Octant< 3 >::computeMorton
uint64_t computeMorton()
Definition: Class_Octant_3D.tpp:417
Class_Octant< 3 >::getPbound
bool getPbound(uint8_t face) const
Definition: Class_Octant_3D.tpp:174
Class_Octant< 3 >::getNormal
void getNormal(uint8_t &iface, vector< int8_t > &normal)
Definition: Class_Octant_3D.tpp:387
Class_Octant< 3 >::getNotBalance
bool getNotBalance() const
Definition: Class_Octant_3D.tpp:191
Class_Octant< 3 >::getX
uint32_t getX() const
Definition: Class_Octant_3D.tpp:134
Class_Global< 3 >::nedges
uint8_t nedges
Definition: Class_Global_3D.tpp:182
Class_Octant< 3 >::getNode
u32vector getNode(uint8_t inode)
Definition: Class_Octant_3D.tpp:364
Class_Global< 3 >::nnodes
uint8_t nnodes
Definition: Class_Global_3D.tpp:183
Class_Octant< 3 >::getSize
uint32_t getSize() const
Definition: Class_Octant_3D.tpp:228
Class_Global< 3 >::nodeface
uint8_t nodeface[8][3]
Definition: Class_Global_3D.tpp:191
Class_Octant< 3 >::getBalance
bool getBalance() const
Definition: Class_Octant_3D.tpp:196
Class_Global< 3 >::nfaces
uint8_t nfaces
Definition: Class_Global_3D.tpp:181
Class_Octant< 3 >::getIsNewR
bool getIsNewR() const
Definition: Class_Octant_3D.tpp:181
Class_Octant< 3 >
Octant class definition - 3D specialization.
Definition: Class_Octant_3D.tpp:32
Class_Octant< 3 >::getY
uint32_t getY() const
Definition: Class_Octant_3D.tpp:139
Class_Octant< 3 >::setBalance
void setBalance(bool balance)
Definition: Class_Octant_3D.tpp:208
Class_Octant< 3 >::getVolume
uint64_t getVolume() const
Definition: Class_Octant_3D.tpp:244
Class_Octant< 3 >::getMarker
int8_t getMarker() const
Definition: Class_Octant_3D.tpp:154
Class_Octant< 3 >::getArea
uint64_t getArea() const
Definition: Class_Octant_3D.tpp:236
Class_Octant< 3 >::getFaceCenter
dvector getFaceCenter(uint8_t iface)
Definition: Class_Octant_3D.tpp:271
Class_Local_Tree
Local octree portion for each process.
Definition: Class_Local_Tree.hpp:55
Class_Octant< 3 >::getBound
bool getBound(uint8_t face) const
Definition: Class_Octant_3D.tpp:160
Class_Octant< 3 >::getLevel
uint8_t getLevel() const
Definition: Class_Octant_3D.tpp:149
Class_Global< 3 >::normals
int8_t normals[6][3]
Definition: Class_Global_3D.tpp:194
Class_Octant< 3 >::getIsNewC
bool getIsNewC() const
Definition: Class_Octant_3D.tpp:186
Class_Octant< 3 >::getNode
void getNode(u32vector &node, uint8_t inode)
Definition: Class_Octant_3D.tpp:344
Class_Octant< 3 >::getZ
uint32_t getZ() const
Definition: Class_Octant_3D.tpp:144
Class_Octant< 3 >::setMarker
void setMarker(int8_t marker)
Definition: Class_Octant_3D.tpp:201
Class_Octant< 3 >::getCenter
dvector getCenter()
Definition: Class_Octant_3D.tpp:254
Class_Global< 3 >::edgeface
uint8_t edgeface[12][2]
Definition: Class_Global_3D.tpp:193
Class_Octant
Octant class definition.
Definition: Class_Octant.hpp:50
Class_Octant< 3 >::getNodes
void getNodes(u32vector2D &nodes)
Definition: Class_Octant_3D.tpp:313
Generated on Mon Oct 5 2015 15:12:11 for PABLO by   doxygen 1.8.8