PABLO  0.1
PArallel Balanced Linear Octree
 All Classes Functions Variables Pages
Class_Local_Tree_3D.tpp
1 
27 // =================================================================================== //
28 // CLASS SPECIALIZATION //
29 // =================================================================================== //
30 template<>
31 class Class_Local_Tree<3>{
32  // ------------------------------------------------------------------------------- //
33  // FRIENDSHIPS ------------------------------------------------------------------- //
34  template<int dim> friend class Class_Para_Tree;
35 
36  // ------------------------------------------------------------------------------- //
37  // TYPEDEFS ----------------------------------------------------------------------- //
38 public:
39  typedef vector<Class_Octant<3> > OctantsType;
40  typedef vector<Class_Intersection<3> > IntersectionsType;
41  typedef vector<uint32_t> u32vector;
42  typedef vector<uint64_t> u64vector;
43  typedef vector<vector<uint32_t> > u32vector2D;
44  typedef vector<vector<uint64_t> > u64vector2D;
45 
46 
47  // ------------------------------------------------------------------------------- //
48  // MEMBERS ----------------------------------------------------------------------- //
49 
50 private:
51  OctantsType octants;
52  OctantsType ghosts;
53  IntersectionsType intersections;
54  u64vector globalidx_ghosts;
55  Class_Octant<3> first_desc;
56  Class_Octant<3> last_desc;
57  uint32_t size_ghosts;
58  uint8_t local_max_depth;
59  uint8_t balance_codim;
62  u32vector last_ghost_bros;
64  // connectivity
65  u32vector2D nodes;
66  u32vector2D connectivity;
68  u32vector2D ghostsnodes;
69  u32vector2D ghostsconnectivity;
72  // ------------------------------------------------------------------------------- //
73  // CONSTRUCTORS ------------------------------------------------------------------ //
74 
75 public:
76  Class_Local_Tree(){
77  Class_Octant<3> oct0;
78  Class_Octant<3> octf(MAX_LEVEL_3D,0,0,0);
79  Class_Octant<3> octl(MAX_LEVEL_3D,global3D.max_length-1,global3D.max_length-1,global3D.max_length-1);
80  octants.resize(1);
81  octants[0] = oct0;
82  first_desc = octf;
83  last_desc = octl;
84  size_ghosts = 0;
85  local_max_depth = 0;
86  balance_codim = 1;
87  };
88 
89  ~Class_Local_Tree(){};
90 
91  // ------------------------------------------------------------------------------- //
92  // METHODS ----------------------------------------------------------------------- //
93 
94  // Basic Get/Set methods --------------------------------------------------------- //
95 
96 private:
97  const Class_Octant<3> & getFirstDesc() const{
98  return first_desc;
99  };
100  const Class_Octant<3> & getLastDesc() const{
101  return last_desc;
102  };
103  uint32_t getSizeGhost() const{
104  return size_ghosts;
105  };
106  uint32_t getNumOctants() const{
107  return octants.size();
108  };
109  uint8_t getLocalMaxDepth() const{ // Get max depth reached in local tree
110  return local_max_depth;
111  };
112  int8_t getMarker(int32_t idx){ // Get refinement/coarsening marker for idx-th octant
113  return octants[idx].getMarker();
114  };
115  uint8_t getLevel(int32_t idx){ // Get refinement/coarsening marker for idx-th octant
116  return octants[idx].getLevel();
117  };
118  uint8_t getGhostLevel(int32_t idx){ // Get refinement/coarsening marker for idx-th ghost octant
119  return ghosts[idx].getLevel();
120  };
121  bool getBalance(int32_t idx){ // Get if balancing-blocked idx-th octant
122  return octants[idx].getNotBalance();
123  };
124 
128  uint8_t getBalanceCodim() const{
129  return balance_codim;
130  };
131 
132  void setMarker(int32_t idx, int8_t marker){ // Set refinement/coarsening marker for idx-th octant
133  octants[idx].setMarker(marker);
134  };
135  void setBalance(int32_t idx, bool balance){ // Set if balancing-blocked idx-th octant
136  octants[idx].setBalance(balance);
137  };
138 
142  void setBalanceCodim(uint8_t b21codim){
143  balance_codim = b21codim;
144  };
145 
146  void setFirstDesc(){
147  OctantsType::const_iterator firstOctant = octants.begin();
148  first_desc = Class_Octant<3>(MAX_LEVEL_3D,firstOctant->x,firstOctant->y,firstOctant->z);
149  };
150  void setLastDesc(){
151  OctantsType::const_iterator lastOctant = octants.end() - 1;
152  uint32_t x,y,z,delta;
153  delta = (uint32_t)pow(2.0,(double)((uint8_t)MAX_LEVEL_3D - lastOctant->level)) - 1;
154  x = lastOctant->x + delta;
155  y = lastOctant->y + delta;
156  z = lastOctant->z + delta;
157  last_desc = Class_Octant<3>(MAX_LEVEL_3D,x,y,z);
158 
159  };
160 
161  //-------------------------------------------------------------------------------- //
162  // Other methods ----------------------------------------------------------------- //
163 
164  Class_Octant<3>& extractOctant(uint32_t idx){
165  return octants[idx];
166  };
167 
168  const Class_Octant<3>& extractOctant(uint32_t idx) const{
169  return octants[idx];
170  };
171 
172  Class_Octant<3>& extractGhostOctant(uint32_t idx) {
173  return ghosts[idx];
174  };
175 
176  const Class_Octant<3>& extractGhostOctant(uint32_t idx) const{
177  return ghosts[idx];
178  };
179 
180  // =================================================================================== //
181 
182  bool refine(){ // Refine local tree: refine one time octants with marker >0
183 
184  // Local variables
185  vector<uint32_t> last_child_index;
186  vector< Class_Octant<3> > children;
187  uint32_t idx, nocts, ilastch;
188  uint32_t offset = 0, blockidx;
189  uint8_t nchm1 = global3D.nchildren-1, ich;
190  bool dorefine = false;
191 
192  nocts = octants.size();
193  for (idx=0; idx<nocts; idx++){
194  if(octants[idx].getMarker() > 0 && octants[idx].getLevel() < MAX_LEVEL_3D){
195  last_child_index.push_back(idx+nchm1+offset);
196  offset += nchm1;
197  }
198  else{
199  // octants[idx].info[12] = false;
200  if (octants[idx].marker > 0)
201  octants[idx].marker = 0;
202  octants[idx].info[15] = true;
203  }
204  }
205  if (offset > 0){
206  octants.resize(octants.size()+offset);
207  blockidx = last_child_index[0]-nchm1;
208  idx = octants.size();
209  ilastch = last_child_index.size()-1;
210  while (idx>blockidx){
211  idx--;
212  if(idx == last_child_index[ilastch]){
213  // if(octants[idx-offset].getMarker() > 0 && octants[idx-offset].getLevel() < MAX_LEVEL_3D){
214  children = octants[idx-offset].buildChildren();
215  for (ich=0; ich<global3D.nchildren; ich++){
216  octants[idx-ich] = (children[nchm1-ich]);
217  }
218  offset -= nchm1;
219  idx -= nchm1;
220  //Update local max depth
221  if (children[0].getLevel() > local_max_depth){
222  local_max_depth = children[0].getLevel();
223  }
224  if (children[0].getMarker() > 0){
225  //More Refinement to do
226  dorefine = true;
227  }
228  //delete []children;
229  if (ilastch != 0){
230  ilastch--;
231  }
232  }
233  else {
234  octants[idx] = octants[idx-offset];
235  }
236  }
237  }
238 #if defined(__INTEL_COMPILER) || defined(__ICC)
239 #else
240  octants.shrink_to_fit();
241 #endif
242  nocts = octants.size();
243 
244  setFirstDesc();
245  setLastDesc();
246 
247  return dorefine;
248 
249  };
250 
251  // =================================================================================== //
252 
253  bool coarse(){ // Coarse local tree: coarse one time family of octants with marker <0
254  // Local variables // (if at least one octant of family has marker>=0 set marker=0 for the entire family)
255  vector<uint32_t> first_child_index;
256  Class_Octant<3> father;
257  uint32_t nocts;
258  uint32_t idx, idx2;
259  uint32_t offset;
260  uint32_t idx2_gh;
261  uint32_t nidx;
262  int8_t markerfather, marker;
263  uint8_t nbro, nend;
264  uint8_t nchm1 = global3D.nchildren-1;
265  bool docoarse = false;
266  bool wstop = false;
267 
268  //------------------------------------------ //
269  // Initialization
270 
271  nbro = nend = 0;
272  nidx = offset = 0;
273 
274  idx2_gh = 0;
275 
276  nocts = octants.size();
277  size_ghosts = ghosts.size();
278 
279  first_child_index.clear();
280 
281  // Init first and last desc (even if already calculated)
282  setFirstDesc();
283  setLastDesc();
284 
285  //------------------------------------------ //
286 
287  // Set index for start and end check for ghosts
288  if (ghosts.size()){
289  while(idx2_gh < size_ghosts && ghosts[idx2_gh].computeMorton() <= last_desc.computeMorton()){
290  idx2_gh++;
291  }
292  idx2_gh = min((size_ghosts-1), idx2_gh);
293  }
294 
295  // Check and coarse internal octants
296  for (idx=0; idx<nocts; idx++){
297  if(octants[idx].getMarker() < 0 && octants[idx].getLevel() > 0){
298  nbro = 0;
299  father = octants[idx].buildFather();
300  // Check if family is to be refined
301  for (idx2=idx; idx2<idx+global3D.nchildren; idx2++){
302  if (idx2<nocts){
303  if(octants[idx2].getMarker() < 0 && octants[idx2].buildFather() == father){
304  nbro++;
305  }
306  }
307  }
308  if (nbro == global3D.nchildren){
309  nidx++;
310  first_child_index.push_back(idx);
311  idx = idx2-1;
312  }
313 // else{
314 // if (idx < (nocts>global3D.nchildren)*(nocts-global3D.nchildren)){
315 // octants[idx].setMarker(0);
316 // octants[idx].info[15] = true;
317 // }
318 // }
319  }
320  // else{
321  // // octants[idx].info[13] = false;
322  // }
323  }
324  uint32_t nblock = nocts;
325  uint32_t nfchild = first_child_index.size();
326  if (nidx!=0){
327  nblock = nocts - nidx*nchm1;
328  nidx = 0;
329  for (idx=0; idx<nblock; idx++){
330  if (nidx < nfchild){
331  if (idx+offset == first_child_index[nidx]){
332  markerfather = -MAX_LEVEL_3D;
333  father = octants[idx+offset].buildFather();
334  for(idx2=0; idx2<global3D.nchildren; idx2++){
335  if (markerfather < octants[idx+offset+idx2].getMarker()+1){
336  markerfather = octants[idx+offset+idx2].getMarker()+1;
337  }
338  for (int iii=0; iii<16; iii++){
339  father.info[iii] = father.info[iii] || octants[idx+offset+idx2].info[iii];
340  }
341  }
342  father.info[13] = true;
343  father.info[15] = true;
344  if (markerfather < 0){
345  docoarse = true;
346  }
347  father.setMarker(markerfather);
348  octants[idx] = father;
349  offset += nchm1;
350  nidx++;
351  }
352  else{
353  octants[idx] = octants[idx+offset];
354  }
355  }
356  else{
357  octants[idx] = octants[idx+offset];
358  }
359  }
360  }
361  // octants.resize(nocts-offset);
362  octants.resize(nblock);
363 #if defined(__INTEL_COMPILER) || defined(__ICC)
364 #else
365  octants.shrink_to_fit();
366 #endif
367  nocts = octants.size();
368 
369  // End on ghosts
370  if (ghosts.size() && nocts > 0){
371 // if ((ghosts[idx2_gh].getMarker() < 0) && (octants[nocts-1].getMarker() < 0)){
372  if (ghosts[idx2_gh].buildFather() == octants[nocts-1].buildFather()){
373  father = ghosts[idx2_gh].buildFather();
374  for (uint32_t iii=0; iii<16; iii++){
375  father.info[iii] = false;
376  }
377  markerfather = ghosts[idx2_gh].getMarker()+1;
378  nbro = 0;
379  idx = idx2_gh;
380  marker = ghosts[idx].getMarker();
381  while(marker < 0 && ghosts[idx].buildFather() == father){
382  nbro++;
383  if (markerfather < ghosts[idx].getMarker()+1){
384  markerfather = ghosts[idx].getMarker()+1;
385  }
386  for (uint32_t iii=0; iii<global3D.nfaces; iii++){
387  father.info[iii] = father.info[iii] || ghosts[idx].info[iii];
388  }
389  father.info[14] = father.info[14] || ghosts[idx].info[14];
390  idx++;
391  if(idx == size_ghosts){
392  break;
393  }
394  marker = ghosts[idx].getMarker();
395 // for (int iii=0; iii<16; iii++){
396 // father.info[iii] = father.info[iii] || ghosts[idx].info[iii];
397 // }
398  }
399  nend = 0;
400  idx = nocts-1;
401  marker = octants[idx].getMarker();
402  while(marker < 0 && octants[idx].buildFather() == father && idx >= 0){
403  nbro++;
404  nend++;
405  if (markerfather < octants[idx].getMarker()+1){
406  markerfather = octants[idx].getMarker()+1;
407  }
408  idx--;
409  marker = octants[idx].getMarker();
410  if (wstop){
411  break;
412  }
413  if (idx==0){
414  wstop = true;
415  }
416  }
417  if (nbro == global3D.nchildren){
418  offset = nend;
419  }
420  else{
421  nend = 0;
422 // for(uint32_t ii=nocts-global3D.nchildren; ii<nocts; ii++){
423 // octants[ii].setMarker(0);
424 // octants[ii].info[15] = true;
425 // }
426  }
427  }
428  if (nend != 0){
429  for (idx=0; idx < nend; idx++){
430 // for (int iii=0; iii<16; iii++){
431 // father.info[iii] = father.info[iii] || octants[nocts-idx-1].info[iii];
432 // }
433  }
434  father.info[13] = true;
435  father.info[15] = true;
436  if (markerfather < 0){
437  docoarse = true;
438  }
439  father.setMarker(markerfather);
440  octants.resize(nocts-offset);
441  octants.push_back(father);
442 #if defined(__INTEL_COMPILER) || defined(__ICC)
443 #else
444  octants.shrink_to_fit();
445 #endif
446  nocts = octants.size();
447  }
448 
449  }
450 
451  // Set final first and last desc
452  if(nocts>0){
453  setFirstDesc();
454  setLastDesc();
455  }
456  return docoarse;
457 
458  };
459 
460  // =================================================================================== //
461 
462  bool refine(u32vector & mapidx){ // Refine local tree: refine one time octants with marker >0
463  // mapidx[i] = index in old octants vector of the i-th octant (index of father if octant is new after)
464  // Local variables
465  vector<uint32_t> last_child_index;
466  vector< Class_Octant<3> > children;
467  uint32_t idx, nocts, ilastch;
468  uint32_t offset = 0, blockidx;
469  uint8_t nchm1 = global3D.nchildren-1, ich;
470  bool dorefine = false;
471 
472  nocts = octants.size();
473  for (idx=0; idx<nocts; idx++){
474  if(octants[idx].getMarker() > 0 && octants[idx].getLevel() < MAX_LEVEL_3D){
475  last_child_index.push_back(idx+nchm1+offset);
476  offset += nchm1;
477  }
478  else{
479  // octants[idx].info[12] = false;
480  if (octants[idx].marker > 0){
481  octants[idx].marker = 0;
482  octants[idx].info[15] = false;
483  }
484  }
485  }
486  if (offset > 0){
487  mapidx.resize(octants.size()+offset);
488 #if defined(__INTEL_COMPILER) || defined(__ICC)
489 #else
490  mapidx.shrink_to_fit();
491 #endif
492  octants.resize(octants.size()+offset);
493  blockidx = last_child_index[0]-nchm1;
494  idx = octants.size();
495  ilastch = last_child_index.size()-1;
496  while (idx>blockidx){
497  // while (idx>0){
498  idx--;
499  // if(octants[idx-offset].getMarker() > 0 && octants[idx-offset].getLevel() < MAX_LEVEL_3D){
500  if(idx == last_child_index[ilastch]){
501  children = octants[idx-offset].buildChildren();
502  for (ich=0; ich<global3D.nchildren; ich++){
503  octants[idx-ich] = (children[nchm1-ich]);
504  mapidx[idx-ich] = mapidx[idx-offset];
505  }
506  offset -= nchm1;
507  idx -= nchm1;
508  //Update local max depth
509  if (children[0].getLevel() > local_max_depth){
510  local_max_depth = children[0].getLevel();
511  }
512  if (children[0].getMarker() > 0){
513  //More Refinement to do
514  dorefine = true;
515  }
516  //delete []children;
517  if (ilastch != 0){
518  ilastch--;
519  }
520  }
521  else {
522  octants[idx] = octants[idx-offset];
523  mapidx[idx] = mapidx[idx-offset];
524  }
525  }
526  }
527 #if defined(__INTEL_COMPILER) || defined(__ICC)
528 #else
529  octants.shrink_to_fit();
530 #endif
531  nocts = octants.size();
532 
533  setFirstDesc();
534  setLastDesc();
535 
536  return dorefine;
537 
538  };
539 
540  // =================================================================================== //
541 
542  bool coarse(u32vector & mapidx){ // Coarse local tree: coarse one time family of octants with marker <0
543  // (if at least one octant of family has marker>=0 set marker=0 for the entire family)
544  // mapidx[i] = index in old octants vector of the i-th octant (index of father if octant is new after)
545  // Local variables
546  vector<uint32_t> first_child_index;
547  Class_Octant<3> father;
548  uint32_t nocts, nocts0;
549  uint32_t idx, idx2;
550  uint32_t offset;
551  uint32_t idx2_gh;
552  uint32_t nidx;
553  int8_t markerfather, marker;
554  uint8_t nbro, nend;
555  uint8_t nchm1 = global3D.nchildren-1;
556  bool docoarse = false;
557  bool wstop = false;
558 
559  //------------------------------------------ //
560  // Initialization
561 
562  nbro = nend = 0;
563  nidx = offset = 0;
564 
565  idx2_gh = 0;
566 
567  nocts = nocts0 = octants.size();
568  size_ghosts = ghosts.size();
569 
570 
571  // Init first and last desc (even if already calculated)
572  setFirstDesc();
573  setLastDesc();
574 
575  //------------------------------------------ //
576 
577  // Set index for start and end check for ghosts
578  if (ghosts.size()){
579  while(idx2_gh < size_ghosts && ghosts[idx2_gh].computeMorton() < last_desc.computeMorton()){
580  idx2_gh++;
581  }
582  idx2_gh = min((size_ghosts-1), idx2_gh);
583  }
584 
585  // Check and coarse internal octants
586  for (idx=0; idx<nocts; idx++){
587  if(octants[idx].getMarker() < 0 && octants[idx].getLevel() > 0){
588  nbro = 0;
589  father = octants[idx].buildFather();
590  // Check if family is to be refined
591  for (idx2=idx; idx2<idx+global3D.nchildren; idx2++){
592  if (idx2<nocts){
593  if(octants[idx2].getMarker() < 0 && octants[idx2].buildFather() == father){
594  nbro++;
595  }
596  }
597  }
598  if (nbro == global3D.nchildren){
599  nidx++;
600  first_child_index.push_back(idx);
601  idx = idx2-1;
602  }
603 // else{
604 // if (idx < (nocts>global3D.nchildren)*(nocts-global3D.nchildren)){
605 // octants[idx].setMarker(0);
606 // octants[idx].info[15] = true;
607 // }
608 // }
609  }
610  // else{
611  // // octants[idx].info[13] = false;
612  // }
613  }
614  uint32_t nblock = nocts;
615  uint32_t nfchild = first_child_index.size();
616  if (nidx!=0){
617  nblock = nocts - nidx*nchm1;
618  nidx = 0;
619  //for (idx=0; idx<nblock; idx++){
620  for (idx=0; idx<nblock; idx++){
621  if (nidx < nfchild){
622  if (idx+offset == first_child_index[nidx]){
623  markerfather = -MAX_LEVEL_3D;
624  father = octants[idx+offset].buildFather();
625  for (uint32_t iii=0; iii<16; iii++){
626  father.info[iii] = false;
627  }
628  for(idx2=0; idx2<global3D.nchildren; idx2++){
629  if (markerfather < octants[idx+offset+idx2].getMarker()+1){
630  markerfather = octants[idx+offset+idx2].getMarker()+1;
631  }
632  for (uint32_t iii=0; iii<16; iii++){
633  father.info[iii] = father.info[iii] || octants[idx+offset+idx2].info[iii];
634  }
635  }
636  father.info[13] = true;
637  father.info[15] = true;
638  father.setMarker(markerfather);
639  if (markerfather < 0){
640  docoarse = true;
641  }
642  octants[idx] = father;
643  mapidx[idx] = mapidx[idx+offset];
644  offset += nchm1;
645  nidx++;
646  }
647  else{
648  octants[idx] = octants[idx+offset];
649  mapidx[idx] = mapidx[idx+offset];
650  }
651  }
652  else{
653  octants[idx] = octants[idx+offset];
654  mapidx[idx] = mapidx[idx+offset];
655  }
656  }
657  }
658  octants.resize(nblock);
659 #if defined(__INTEL_COMPILER) || defined(__ICC)
660 #else
661  octants.shrink_to_fit();
662 #endif
663  nocts = octants.size();
664  mapidx.resize(nocts);
665 #if defined(__INTEL_COMPILER) || defined(__ICC)
666 #else
667  mapidx.shrink_to_fit();
668 #endif
669  // End on ghosts
670  if (ghosts.size() && nocts > 0){
671 // if ((ghosts[idx2_gh].getMarker() < 0) && (octants[nocts-1].getMarker() < 0)){
672  if (ghosts[idx2_gh].buildFather() == octants[nocts-1].buildFather()){
673  father = ghosts[idx2_gh].buildFather();
674  for (uint32_t iii=0; iii<16; iii++){
675  father.info[iii] = false;
676  }
677  markerfather = ghosts[idx2_gh].getMarker()+1;
678  nbro = 0;
679  idx = idx2_gh;
680  marker = ghosts[idx].getMarker();
681  while(marker < 0 && ghosts[idx].buildFather() == father){
682  nbro++;
683  if (markerfather < ghosts[idx].getMarker()+1){
684  markerfather = ghosts[idx].getMarker()+1;
685  }
686  for (uint32_t iii=0; iii<global3D.nfaces; iii++){
687  father.info[iii] = father.info[iii] || ghosts[idx].info[iii];
688  }
689  father.info[14] = father.info[14] || ghosts[idx].info[14];
690  idx++;
691  if(idx == size_ghosts){
692  break;
693  }
694  marker = ghosts[idx].getMarker();
695 // for (uint32_t iii=0; iii<16; iii++){
696 // father.info[iii] = father.info[iii] || ghosts[idx].info[iii];
697 // }
698  }
699  nend = 0;
700  idx = nocts-1;
701  marker = octants[idx].getMarker();
702  while(marker < 0 && octants[idx].buildFather() == father && idx >= 0){
703  nbro++;
704  nend++;
705  if (markerfather < octants[idx].getMarker()+1){
706  markerfather = octants[idx].getMarker()+1;
707  }
708  idx--;
709  marker = octants[idx].getMarker();
710  if (wstop){
711  break;
712  }
713  if (idx==0){
714  wstop = true;
715  }
716  }
717  if (nbro == global3D.nchildren){
718  offset = nend;
719  }
720  else{
721  nend = 0;
722 // for(uint32_t ii=nocts-global3D.nchildren; ii<nocts; ii++){
723 // octants[ii].setMarker(0);
724 // octants[ii].info[15] = true;
725 // }
726  }
727  }
728  if (nend != 0){
729  for (idx=0; idx < nend; idx++){
730  for (uint32_t iii=0; iii<16; iii++){
731  father.info[iii] = father.info[iii] || octants[nocts-idx-1].info[iii];
732  }
733  }
734  father.info[13] = true;
735  father.info[15] = true;
736  if (markerfather < 0){
737  docoarse = true;
738  }
739  father.setMarker(markerfather);
740  octants.resize(nocts-offset);
741  octants.push_back(father);
742 #if defined(__INTEL_COMPILER) || defined(__ICC)
743 #else
744  octants.shrink_to_fit();
745 #endif
746  nocts = octants.size();
747  mapidx.resize(nocts);
748 #if defined(__INTEL_COMPILER) || defined(__ICC)
749 #else
750  mapidx.shrink_to_fit();
751 #endif
752  }
753 
754  }
755 
756  // Set final first and last desc
757  if(nocts>0){
758  setFirstDesc();
759  setLastDesc();
760  }
761  return docoarse;
762 
763  };
764 
765  // =================================================================================== //
766 
767  // Global refine of octree (one level every element)
768  bool globalRefine(){
769 
770  // Local variables
771  vector<uint32_t> last_child_index;
772  vector< Class_Octant<3> > children;
773  uint32_t idx, nocts, ilastch;
774  uint32_t offset = 0, blockidx;
775  uint8_t nchm1 = global3D.nchildren-1, ich;
776  bool dorefine = false;
777 
778  nocts = octants.size();
779  for (idx=0; idx<nocts; idx++){
780  octants[idx].setMarker(1);
781  if(octants[idx].getMarker() > 0 && octants[idx].getLevel() < MAX_LEVEL_3D){
782  last_child_index.push_back(idx+nchm1+offset);
783  offset += nchm1;
784  }
785  else{
786  if (octants[idx].marker > 0)
787  octants[idx].marker = 0;
788  octants[idx].info[15] = true;
789  }
790  }
791  if (offset > 0){
792  octants.resize(octants.size()+offset);
793  blockidx = last_child_index[0]-nchm1;
794  idx = octants.size();
795  ilastch = last_child_index.size()-1;
796  while (idx>blockidx){
797  idx--;
798  if(idx == last_child_index[ilastch]){
799  // if(octants[idx-offset].getMarker() > 0 && octants[idx-offset].getLevel() < MAX_LEVEL_3D){
800  children = octants[idx-offset].buildChildren();
801  for (ich=0; ich<global3D.nchildren; ich++){
802  octants[idx-ich] = (children[nchm1-ich]);
803  }
804  offset -= nchm1;
805  idx -= nchm1;
806  //Update local max depth
807  if (children[0].getLevel() > local_max_depth){
808  local_max_depth = children[0].getLevel();
809  }
810  if (children[0].getMarker() > 0){
811  //More Refinement to do
812  dorefine = true;
813  }
814  if (ilastch != 0){
815  ilastch--;
816  }
817  }
818  else {
819  octants[idx] = octants[idx-offset];
820  }
821  }
822  }
823 #if defined(__INTEL_COMPILER) || defined(__ICC)
824 #else
825  octants.shrink_to_fit();
826 #endif
827  nocts = octants.size();
828 
829  setFirstDesc();
830  setLastDesc();
831 
832  return dorefine;
833 
834  };
835 
836  // =================================================================================== //
837 
838  // Global coarse of octree (every marker set =-1)
839  bool globalCoarse(){
840  //Local Variables
841  vector<uint32_t> first_child_index;
842  Class_Octant<3> father;
843  uint32_t nocts;
844  uint32_t idx, idx2;
845  uint32_t offset;
846  uint32_t idx2_gh;
847  uint32_t nidx;
848  int8_t markerfather, marker;
849  uint8_t nbro, nend;
850  uint8_t nchm1 = global3D.nchildren-1;
851  bool docoarse = false;
852  bool wstop = false;
853 
854  //------------------------------------------ //
855  // Initialization
856 
857  nbro = nend = 0;
858  nidx = offset = 0;
859 
860  idx2_gh = 0;
861 
862  nocts = octants.size();
863  size_ghosts = ghosts.size();
864 
865  // Init first and last desc (even if already calculated)
866  setFirstDesc();
867  setLastDesc();
868 
869  //------------------------------------------ //
870 
871  // Set index for start and end check for ghosts
872  if (ghosts.size()){
873  while(idx2_gh < size_ghosts && ghosts[idx2_gh].computeMorton() <= last_desc.computeMorton()){
874  idx2_gh++;
875  }
876  idx2_gh = min((size_ghosts-1), idx2_gh);
877  }
878 
879  // Check and coarse internal octants
880  for (idx=0; idx<nocts; idx++){
881  octants[idx].setMarker(-1);
882  if(octants[idx].getMarker() < 0 && octants[idx].getLevel() > 0){
883  nbro = 0;
884  father = octants[idx].buildFather();
885  // Check if family is to be refined
886  for (idx2=idx; idx2<idx+global3D.nchildren; idx2++){
887  if (idx2<nocts){
888  octants[idx2].setMarker(-1);
889  if(octants[idx2].getMarker() < 0 && octants[idx2].buildFather() == father){
890  nbro++;
891  }
892  }
893  }
894  if (nbro == global3D.nchildren){
895  nidx++;
896  first_child_index.push_back(idx);
897  idx = idx2-1;
898  }
899  else{
900  if (idx < (nocts>global3D.nchildren)*(nocts-global3D.nchildren)){
901  octants[idx].setMarker(0);
902  octants[idx].info[15] = true;
903  }
904  }
905  }
906  // else{
907  // // octants[idx].info[13] = false;
908  // }
909  }
910  uint32_t nblock = nocts;
911  uint32_t nfchild = first_child_index.size();
912  if (nidx!=0){
913  nblock = nocts - nidx*nchm1;
914  nidx = 0;
915  for (idx=0; idx<nblock; idx++){
916  if (nidx < nfchild){
917  if (idx+offset == first_child_index[nidx]){
918  markerfather = -MAX_LEVEL_3D;
919  father = octants[idx+offset].buildFather();
920  for(idx2=0; idx2<global3D.nchildren; idx2++){
921  if (markerfather < octants[idx+offset+idx2].getMarker()+1){
922  markerfather = octants[idx+offset+idx2].getMarker()+1;
923  }
924  for (uint32_t iii=0; iii<16; iii++){
925  father.info[iii] = father.info[iii] || octants[idx+offset+idx2].info[iii];
926  }
927  }
928  father.info[13] = true;
929  if (markerfather < 0){
930  docoarse = true;
931  }
932  father.setMarker(markerfather);
933  octants[idx] = father;
934  offset += nchm1;
935  nidx++;
936  }
937  else{
938  octants[idx] = octants[idx+offset];
939  }
940  }
941  else{
942  octants[idx] = octants[idx+offset];
943  }
944  }
945  }
946  octants.resize(nblock);
947 #if defined(__INTEL_COMPILER) || defined(__ICC)
948 #else
949  octants.shrink_to_fit();
950 #endif
951  nocts = octants.size();
952 
953  // End on ghosts
954  if (ghosts.size() && nocts > 0){
955  ghosts[idx2_gh].setMarker(-1);
956  if ((ghosts[idx2_gh].getMarker() < 0) && (octants[nocts-1].getMarker() < 0)){
957  father = ghosts[idx2_gh].buildFather();
958  for (uint32_t iii=0; iii<16; iii++){
959  father.info[iii] = false;
960  }
961  markerfather = ghosts[idx2_gh].getMarker()+1;
962  nbro = 0;
963  idx = idx2_gh;
964  ghosts[idx].setMarker(-1);
965  marker = ghosts[idx].getMarker();
966  while(marker < 0 && ghosts[idx].buildFather() == father){
967  nbro++;
968  if (markerfather < ghosts[idx].getMarker()+1){
969  markerfather = ghosts[idx].getMarker()+1;
970  }
971  idx++;
972  if(idx == size_ghosts){
973  break;
974  }
975  ghosts[idx].setMarker(-1);
976  marker = ghosts[idx].getMarker();
977  for (uint32_t iii=0; iii<16; iii++){
978  father.info[iii] = father.info[iii] || ghosts[idx].info[iii];
979  }
980  }
981  nend = 0;
982  idx = nocts-1;
983  octants[idx].setMarker(-1);
984  marker = octants[idx].getMarker();
985  while(marker < 0 && octants[idx].buildFather() == father && idx >= 0){
986  nbro++;
987  nend++;
988  if (markerfather < octants[idx].getMarker()+1){
989  markerfather = octants[idx].getMarker()+1;
990  }
991  idx--;
992  octants[idx].setMarker(-1);
993  marker = octants[idx].getMarker();
994  if (wstop){
995  break;
996  }
997  if (idx==0){
998  wstop = true;
999  }
1000  }
1001  if (nbro == global3D.nchildren){
1002  offset = nend;
1003  }
1004  else{
1005  nend = 0;
1006  for(uint32_t ii=nocts-global3D.nchildren; ii<nocts; ii++){
1007  octants[ii].setMarker(0);
1008  octants[ii].info[15] = true;
1009  }
1010  }
1011  }
1012  if (nend != 0){
1013  for (idx=0; idx < nend; idx++){
1014  for (int iii=0; iii<16; iii++){
1015  father.info[iii] = father.info[iii] || octants[nocts-idx-1].info[iii];
1016  }
1017  }
1018  father.info[13] = true;
1019  if (markerfather < 0){
1020  docoarse = true;
1021  }
1022  father.setMarker(markerfather);
1023  octants.resize(nocts-offset);
1024  octants.push_back(father);
1025 #if defined(__INTEL_COMPILER) || defined(__ICC)
1026 #else
1027  octants.shrink_to_fit();
1028 #endif
1029  nocts = octants.size();
1030  }
1031 
1032  }
1033 
1034  // Set final first and last desc
1035  if(nocts>0){
1036  setFirstDesc();
1037  setLastDesc();
1038  }
1039  return docoarse;
1040 
1041  };
1042 
1043  // =================================================================================== //
1044 
1045  // One level global refine with mapidx
1046  bool globalRefine(u32vector & mapidx){
1047 
1048  vector<uint32_t> last_child_index;
1049  vector< Class_Octant<3> > children;
1050  uint32_t idx, nocts, ilastch;
1051  uint32_t offset = 0, blockidx;
1052  uint8_t nchm1 = global3D.nchildren-1, ich;
1053  bool dorefine = false;
1054 
1055  nocts = octants.size();
1056  for (idx=0; idx<nocts; idx++){
1057  octants[idx].setMarker(1);
1058  if(octants[idx].getMarker() > 0 && octants[idx].getLevel() < MAX_LEVEL_3D){
1059  last_child_index.push_back(idx+nchm1+offset);
1060  offset += nchm1;
1061  }
1062  else{
1063  // octants[idx].info[12] = false;
1064  if (octants[idx].marker > 0){
1065  octants[idx].marker = 0;
1066  octants[idx].info[15] = false;
1067  }
1068  }
1069  }
1070  if (offset > 0){
1071  mapidx.resize(octants.size()+offset);
1072 #if defined(__INTEL_COMPILER) || defined(__ICC)
1073 #else
1074  mapidx.shrink_to_fit();
1075 #endif
1076  octants.resize(octants.size()+offset);
1077  blockidx = last_child_index[0]-nchm1;
1078  idx = octants.size();
1079  ilastch = last_child_index.size()-1;
1080  while (idx>blockidx){
1081  // while (idx>0){
1082  idx--;
1083  // if(octants[idx-offset].getMarker() > 0 && octants[idx-offset].getLevel() < MAX_LEVEL_3D){
1084  if(idx == last_child_index[ilastch]){
1085  children = octants[idx-offset].buildChildren();
1086  for (ich=0; ich<global3D.nchildren; ich++){
1087  octants[idx-ich] = (children[nchm1-ich]);
1088  mapidx[idx-ich] = mapidx[idx-offset];
1089  }
1090  offset -= nchm1;
1091  idx -= nchm1;
1092  //Update local max depth
1093  if (children[0].getLevel() > local_max_depth){
1094  local_max_depth = children[0].getLevel();
1095  }
1096  if (children[0].getMarker() > 0){
1097  //More Refinement to do
1098  dorefine = true;
1099  }
1100  //delete []children;
1101  if (ilastch != 0){
1102  ilastch--;
1103  }
1104  }
1105  else {
1106  octants[idx] = octants[idx-offset];
1107  mapidx[idx] = mapidx[idx-offset];
1108  }
1109  }
1110  }
1111 #if defined(__INTEL_COMPILER) || defined(__ICC)
1112 #else
1113  octants.shrink_to_fit();
1114 #endif
1115  nocts = octants.size();
1116 
1117  setFirstDesc();
1118  setLastDesc();
1119 
1120  return dorefine;
1121 
1122  };
1123 
1124  // =================================================================================== //
1125 
1126  bool globalCoarse(u32vector & mapidx){ // Coarse local tree: coarse one time family of octants with marker <0
1127  // (if at least one octant of family has marker>=0 set marker=0 for the entire family)
1128  // mapidx[i] = index in old octants vector of the i-th octant (index of father if octant is new after)
1129  // Local variables
1130  vector<uint32_t> first_child_index;
1131  Class_Octant<3> father;
1132  uint32_t nocts, nocts0;
1133  uint32_t idx, idx2;
1134  uint32_t offset;
1135  uint32_t idx2_gh;
1136  uint32_t nidx;
1137  int8_t markerfather, marker;
1138  uint8_t nbro, nstart, nend;
1139  uint8_t nchm1 = global3D.nchildren-1;
1140  bool docoarse = false;
1141  bool wstop = false;
1142 
1143  //------------------------------------------ //
1144  // Initialization
1145 
1146  nbro = nstart = nend = 0;
1147  nidx = offset = 0;
1148 
1149  idx2_gh = 0;
1150 
1151  nocts = nocts0 = octants.size();
1152  size_ghosts = ghosts.size();
1153 
1154 
1155  // Init first and last desc (even if already calculated)
1156  setFirstDesc();
1157  setLastDesc();
1158 
1159  //------------------------------------------ //
1160 
1161  // Set index for start and end check for ghosts
1162  if (ghosts.size()){
1163  while(idx2_gh < size_ghosts && ghosts[idx2_gh].computeMorton() < last_desc.computeMorton()){
1164  idx2_gh++;
1165  }
1166  idx2_gh = min((size_ghosts-1), idx2_gh);
1167  }
1168 
1169  // Check and coarse internal octants
1170  for (idx=0; idx<nocts; idx++){
1171  octants[idx].setMarker(-1);
1172  if(octants[idx].getMarker() < 0 && octants[idx].getLevel() > 0){
1173  nbro = 0;
1174  father = octants[idx].buildFather();
1175  // Check if family is to be refined
1176  for (idx2=idx; idx2<idx+global3D.nchildren; idx2++){
1177  if (idx2<nocts){
1178  octants[idx2].setMarker(-1);
1179  if(octants[idx2].getMarker() < 0 && octants[idx2].buildFather() == father){
1180  nbro++;
1181  }
1182  }
1183  }
1184  if (nbro == global3D.nchildren){
1185  nidx++;
1186  first_child_index.push_back(idx);
1187  idx = idx2-1;
1188  }
1189  else{
1190  if (idx < (nocts>global3D.nchildren)*(nocts-global3D.nchildren)){
1191  octants[idx].setMarker(0);
1192  octants[idx].info[15] = true;
1193  }
1194  }
1195  }
1196  // else{
1197  // // octants[idx].info[13] = false;
1198  // }
1199  }
1200  if (nidx!=0){
1201  uint32_t nblock = nocts - nidx*nchm1 - nstart;
1202  uint32_t nfchild = first_child_index.size();
1203  nidx = 0;
1204  //for (idx=0; idx<nblock; idx++){
1205  for (idx=0; idx<nocts-offset; idx++){
1206  if (nidx < nfchild){
1207  if (idx+offset == first_child_index[nidx]){
1208  markerfather = -MAX_LEVEL_3D;
1209  father = octants[idx+offset].buildFather();
1210  for (uint32_t iii=0; iii<16; iii++){
1211  father.info[iii] = false;
1212  }
1213  for(idx2=0; idx2<global3D.nchildren; idx2++){
1214  if (markerfather < octants[idx+offset+idx2].getMarker()+1){
1215  markerfather = octants[idx+offset+idx2].getMarker()+1;
1216  }
1217  for (uint32_t iii=0; iii<16; iii++){
1218  father.info[iii] = father.info[iii] || octants[idx+offset+idx2].info[iii];
1219  }
1220  }
1221  father.info[13] = true;
1222  father.setMarker(markerfather);
1223  if (markerfather < 0){
1224  docoarse = true;
1225  }
1226  octants[idx] = father;
1227  mapidx[idx] = mapidx[idx+offset];
1228  offset += nchm1;
1229  nidx++;
1230  }
1231  else{
1232  octants[idx] = octants[idx+offset];
1233  mapidx[idx] = mapidx[idx+offset];
1234  }
1235  }
1236  else{
1237  octants[idx] = octants[idx+offset];
1238  mapidx[idx] = mapidx[idx+offset];
1239  }
1240  }
1241  }
1242  octants.resize(nocts-offset);
1243 #if defined(__INTEL_COMPILER) || defined(__ICC)
1244 #else
1245  octants.shrink_to_fit();
1246 #endif
1247  nocts = octants.size();
1248  mapidx.resize(nocts);
1249 #if defined(__INTEL_COMPILER) || defined(__ICC)
1250 #else
1251  mapidx.shrink_to_fit();
1252 #endif
1253  // End on ghosts
1254  if (ghosts.size() && nocts > 0){
1255  ghosts[idx2_gh].setMarker(-1);
1256  if ((ghosts[idx2_gh].getMarker() < 0) && (octants[nocts-1].getMarker() < 0)){
1257  father = ghosts[idx2_gh].buildFather();
1258  markerfather = ghosts[idx2_gh].getMarker()+1;
1259  nbro = 0;
1260  idx = idx2_gh;
1261  ghosts[idx].setMarker(-1);
1262  marker = ghosts[idx].getMarker();
1263  while(marker < 0 && ghosts[idx].buildFather() == father){
1264  nbro++;
1265  if (markerfather < ghosts[idx].getMarker()+1){
1266  markerfather = ghosts[idx].getMarker()+1;
1267  }
1268  idx++;
1269  if(idx == size_ghosts){
1270  break;
1271  }
1272  ghosts[idx].setMarker(-1);
1273  marker = ghosts[idx].getMarker();
1274  }
1275  nend = 0;
1276  idx = nocts-1;
1277  octants[idx].setMarker(-1);
1278  marker = octants[idx].getMarker();
1279  while(marker < 0 && octants[idx].buildFather() == father && idx >= 0){
1280  nbro++;
1281  nend++;
1282  if (markerfather < octants[idx].getMarker()+1){
1283  markerfather = octants[idx].getMarker()+1;
1284  }
1285  idx--;
1286  octants[idx].setMarker(-1);
1287  marker = octants[idx].getMarker();
1288  if (wstop){
1289  break;
1290  }
1291  if (idx==0){
1292  wstop=true;
1293  }
1294  }
1295  if (nbro == global3D.nchildren){
1296  offset = nend;
1297  }
1298  else{
1299  nend = 0;
1300  for(uint32_t ii=nocts-global3D.nchildren; ii<nocts; ii++){
1301  octants[ii].setMarker(0);
1302  octants[ii].info[15] = true;
1303  }
1304  }
1305  }
1306  if (nend != 0){
1307  for (uint32_t iii=0; iii<16; iii++){
1308  father.info[iii] = false;
1309  }
1310  for (idx=0; idx < nend; idx++){
1311  for (uint32_t iii=0; iii<16; iii++){
1312  father.info[iii] = father.info[iii] || octants[nocts-idx-1].info[iii];
1313  }
1314  }
1315  father.info[13] = true;
1316  if (markerfather < 0){
1317  docoarse = true;
1318  }
1319  father.setMarker(markerfather);
1320  octants.resize(nocts-offset);
1321  octants.push_back(father);
1322 #if defined(__INTEL_COMPILER) || defined(__ICC)
1323 #else
1324  octants.shrink_to_fit();
1325 #endif
1326  nocts = octants.size();
1327  mapidx.resize(nocts);
1328 #if defined(__INTEL_COMPILER) || defined(__ICC)
1329 #else
1330  mapidx.shrink_to_fit();
1331 #endif
1332  }
1333  }
1334 
1335  // Set final first and last desc
1336  if(nocts>0){
1337  setFirstDesc();
1338  setLastDesc();
1339  }
1340  return docoarse;
1341 
1342  };
1343 
1344  // =================================================================================== //
1345 
1346  void checkCoarse(uint64_t lastDescPre, // Delete overlapping octants after coarse local tree. Check first and last descendants
1347  uint64_t firstDescPost){ // of process before and after the local process
1348  uint32_t idx;
1349  uint32_t nocts;
1350  uint64_t Morton;
1351  uint8_t toDelete = 0;
1352 
1353  nocts = getNumOctants();
1354  idx = 0;
1355  Morton = octants[idx].computeMorton();
1356  while(Morton <= lastDescPre && idx < nocts && Morton != 0){
1357  // To delete, the father is in proc before me
1358  toDelete++;
1359  idx++;
1360  Morton = octants[idx].computeMorton();
1361  }
1362  for(idx=0; idx<nocts-toDelete; idx++){
1363  octants[idx] = octants[idx+toDelete];
1364  }
1365  octants.resize(nocts-toDelete);
1366 #if defined(__INTEL_COMPILER) || defined(__ICC)
1367 #else
1368  octants.shrink_to_fit();
1369 #endif
1370  nocts = getNumOctants();
1371 
1372  setFirstDesc();
1373  setLastDesc();
1374 
1375  };
1376 
1377  // =================================================================================== //
1378 
1379  void checkCoarse(uint64_t lastDescPre, // Delete overlapping octants after coarse local tree. Check first and last descendants
1380  uint64_t firstDescPost,
1381  u32vector & mapidx){ // of process before and after the local process
1382  uint32_t idx;
1383  uint32_t nocts;
1384  uint64_t Morton;
1385  uint8_t toDelete = 0;
1386 
1387  nocts = getNumOctants();
1388  idx = 0;
1389  Morton = octants[idx].computeMorton();
1390  while(Morton <= lastDescPre && idx < nocts && Morton != 0){
1391  // To delete, the father is in proc before me
1392  toDelete++;
1393  idx++;
1394  Morton = octants[idx].computeMorton();
1395  }
1396  for(idx=0; idx<nocts-toDelete; idx++){
1397  octants[idx] = octants[idx+toDelete];
1398  mapidx[idx] = mapidx[idx+toDelete];
1399  }
1400  octants.resize(nocts-toDelete);
1401 #if defined(__INTEL_COMPILER) || defined(__ICC)
1402 #else
1403  octants.shrink_to_fit();
1404 #endif
1405  mapidx.resize(nocts-toDelete);
1406 #if defined(__INTEL_COMPILER) || defined(__ICC)
1407 #else
1408  mapidx.shrink_to_fit();
1409 #endif
1410  nocts = getNumOctants();
1411 
1412  setFirstDesc();
1413  setLastDesc();
1414 
1415  };
1416 
1417  // =================================================================================== //
1418 
1419  void updateLocalMaxDepth(){ // Update max depth reached in local tree
1420  uint32_t noctants = getNumOctants();
1421  uint32_t i;
1422 
1423  local_max_depth = 0;
1424  for(i = 0; i < noctants; i++){
1425  if(octants[i].getLevel() > local_max_depth){
1426  local_max_depth = octants[i].getLevel();
1427  }
1428  }
1429 
1430  };
1431 
1432  // =================================================================================== //
1433 
1434  void findNeighbours(uint32_t idx, // Finds neighbours of idx-th octant through iface in vector octants.
1435  uint8_t iface, // Returns a vector (empty if iface is a bound face) with the index of neighbours
1436  u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
1437  vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
1438 
1439  uint64_t Morton, Mortontry;
1440  uint32_t noctants = getNumOctants();
1441  uint32_t idxtry;
1442  Class_Octant<3>* oct = &octants[idx];
1443  uint32_t size = oct->getSize();
1444 
1445  //Alternative to switch case
1446 // int8_t cx = int8_t((iface<2)*(int8_t(2*iface-1)));
1447 // int8_t cy = int8_t((iface<4)*(int8_t(iface/2))*(int8_t(2*iface-5)));
1448 // int8_t cz = int8_t((int8_t(iface/4))*(int8_t(2*iface-9)));
1449  int8_t cx = global3D.normals[iface][0];
1450  int8_t cy = global3D.normals[iface][1];
1451  int8_t cz = global3D.normals[iface][2];
1452 
1453  isghost.clear();
1454  neighbours.clear();
1455 
1456  // Default if iface is nface<iface<0
1457  if (iface < 0 || iface > global3D.nfaces){
1458  return;
1459  }
1460 
1461  // Check if octants face is a process boundary
1462  if (oct->info[global3D.nfaces+iface] == false){
1463 
1464  // Check if octants face is a boundary
1465  if (oct->info[iface] == false){
1466 
1467  //Build Morton number of virtual neigh of same size
1468  Class_Octant<3> samesizeoct(oct->level, int32_t(oct->x)+int32_t(cx*size), int32_t(oct->y)+int32_t(cy*size), int32_t(oct->z)+int32_t(cz*size));
1469  Morton = samesizeoct.computeMorton();
1470  // Search morton in octants
1471  // If a even face morton is lower than morton of oct, if odd higher
1472  // ---> can i search only before or after idx in octants
1473  int32_t jump = (oct->computeMorton() > Morton) ? int32_t(idx/2+1) : int32_t((noctants -idx)/2+1);
1474  idxtry = uint32_t(idx +((oct->computeMorton()<Morton)-(oct->computeMorton()>Morton))*jump);
1475  if (idxtry > noctants-1) idxtry = noctants-1;
1476  Mortontry = oct->computeMorton();
1477  while(abs(jump) > 0){
1478  Mortontry = octants[idxtry].computeMorton();
1479  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
1480  idxtry += jump;
1481  if (idxtry > noctants-1){
1482  if (jump > 0){
1483  idxtry = noctants - 1;
1484  jump = 0;
1485  }
1486  else if (jump < 0){
1487  idxtry = 0;
1488  jump = 0;
1489  }
1490  }
1491  }
1492  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
1493  //Found neighbour of same size
1494  isghost.push_back(false);
1495  neighbours.push_back(idxtry);
1496  return;
1497  }
1498  else{
1499  // Step until the mortontry lower than morton (one idx of distance)
1500  {
1501  while(octants[idxtry].computeMorton() < Morton){
1502  idxtry++;
1503  if(idxtry > noctants-1){
1504  idxtry = noctants-1;
1505  break;
1506  }
1507  }
1508  while(octants[idxtry].computeMorton() > Morton){
1509  idxtry--;
1510  if(idxtry > noctants-1){
1511  idxtry = 0;
1512  break;
1513  }
1514  }
1515  }
1516  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
1517  //Found neighbour of same size
1518  isghost.push_back(false);
1519  neighbours.push_back(idxtry);
1520  return;
1521  }
1522  // Compute Last discendent of virtual octant of same size
1523  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
1524  uint64_t Mortonlast = last_desc.computeMorton();
1525  Mortontry = octants[idxtry].computeMorton();
1526  int32_t Dx, Dy, Dz;
1527  int32_t Dxstar, Dystar, Dzstar;
1528  while(Mortontry < Mortonlast && idxtry < noctants){
1529  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
1530  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
1531  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
1532  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
1533  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
1534  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
1535 
1536  uint32_t x0 = oct->x;
1537  uint32_t x1 = x0 + size;
1538  uint32_t y0 = oct->y;
1539  uint32_t y1 = y0 + size;
1540  uint32_t z0 = oct->z;
1541  uint32_t z1 = z0 + size;
1542  uint32_t x0try = octants[idxtry].x;
1543  uint32_t x1try = x0try + octants[idxtry].getSize();
1544  uint32_t y0try = octants[idxtry].y;
1545  uint32_t y1try = y0try + octants[idxtry].getSize();
1546  uint32_t z0try = octants[idxtry].z;
1547  uint32_t z1try = z0try + octants[idxtry].getSize();
1548  uint8_t level = oct->level;
1549  uint8_t leveltry = octants[idxtry].getLevel();
1550 
1551  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
1552  if (leveltry > level){
1553  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
1554  neighbours.push_back(idxtry);
1555  isghost.push_back(false);
1556  }
1557  }
1558  if (leveltry < level){
1559  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
1560  neighbours.push_back(idxtry);
1561  isghost.push_back(false);
1562  }
1563  }
1564  }
1565 
1566  idxtry++;
1567  if(idxtry>noctants-1){
1568  break;
1569  }
1570  Mortontry = octants[idxtry].computeMorton();
1571  }
1572  return;
1573  }
1574  }
1575  else{
1576  // Boundary Face
1577  return;
1578  }
1579  }
1580  //--------------------------------------------------------------- //
1581  //--------------------------------------------------------------- //
1582  else{
1583  // Check if octants face is a boundary
1584  if (oct->info[iface] == false){
1585  // IF OCTANT FACE IS A PROCESS BOUNDARY SEARCH ALSO IN GHOSTS
1586 
1587  if (ghosts.size()>0){
1588  // Search in ghosts
1589 
1590  uint32_t idxghost = uint32_t(size_ghosts/2);
1591  Class_Octant<3>* octghost = &ghosts[idxghost];
1592 
1593  //Build Morton number of virtual neigh of same size
1594  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
1595  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
1596  // Search morton in octants
1597  // If a even face morton is lower than morton of oct, if odd higher
1598  // ---> can i search only before or after idx in octants
1599  int32_t jump = (octghost->computeMorton() > Morton) ? int32_t(idxghost/2+1) : int32_t((size_ghosts -idxghost)/2+1);
1600  idxtry = uint32_t(idxghost +((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
1601  if (idxtry > ghosts.size()-1) idxtry = ghosts.size()-1;
1602  while(abs(jump) > 0){
1603  Mortontry = ghosts[idxtry].computeMorton();
1604  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
1605  idxtry += jump;
1606  if (idxtry > ghosts.size()-1){
1607  if (jump > 0){
1608  idxtry = ghosts.size() - 1;
1609  jump = 0;
1610  }
1611  else if (jump < 0){
1612  idxtry = 0;
1613  jump = 0;
1614  }
1615  }
1616  }
1617  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
1618  //Found neighbour of same size
1619  isghost.push_back(true);
1620  neighbours.push_back(idxtry);
1621  return;
1622  }
1623  else{
1624  // Step until the mortontry lower than morton (one idx of distance)
1625  {
1626  while(ghosts[idxtry].computeMorton() < Morton){
1627  idxtry++;
1628  if(idxtry > ghosts.size()-1){
1629  idxtry = ghosts.size()-1;
1630  break;
1631  }
1632  }
1633  while(ghosts[idxtry].computeMorton() > Morton){
1634  idxtry--;
1635  if(idxtry > ghosts.size()-1){
1636  idxtry = 0;
1637  break;
1638  }
1639  }
1640  }
1641  if(idxtry < size_ghosts){
1642  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
1643  //Found neighbour of same size
1644  isghost.push_back(true);
1645  neighbours.push_back(idxtry);
1646  return;
1647  }
1648  // Compute Last discendent of virtual octant of same size
1649  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
1650  uint64_t Mortonlast = last_desc.computeMorton();
1651  Mortontry = ghosts[idxtry].computeMorton();
1652  int32_t Dx, Dy, Dz;
1653  int32_t Dxstar, Dystar, Dzstar;
1654  while(Mortontry < Mortonlast && idxtry < size_ghosts){
1655  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(ghosts[idxtry].x));
1656  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(ghosts[idxtry].y));
1657  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(ghosts[idxtry].z));
1658  Dxstar = int32_t((cx-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cx+1)/2)*size;
1659  Dystar = int32_t((cy-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cy+1)/2)*size;
1660  Dzstar = int32_t((cz-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cz+1)/2)*size;
1661 
1662  uint32_t x0 = oct->x;
1663  uint32_t x1 = x0 + size;
1664  uint32_t y0 = oct->y;
1665  uint32_t y1 = y0 + size;
1666  uint32_t z0 = oct->z;
1667  uint32_t z1 = z0 + size;
1668  uint32_t x0try = ghosts[idxtry].x;
1669  uint32_t x1try = x0try + ghosts[idxtry].getSize();
1670  uint32_t y0try = ghosts[idxtry].y;
1671  uint32_t y1try = y0try + ghosts[idxtry].getSize();
1672  uint32_t z0try = ghosts[idxtry].z;
1673  uint32_t z1try = z0try + ghosts[idxtry].getSize();
1674  uint8_t level = oct->level;
1675  uint8_t leveltry = ghosts[idxtry].getLevel();
1676 
1677  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
1678  if (leveltry > level){
1679  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
1680  neighbours.push_back(idxtry);
1681  isghost.push_back(true);
1682  }
1683  }
1684  if (leveltry < level){
1685  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
1686  neighbours.push_back(idxtry);
1687  isghost.push_back(true);
1688  }
1689  }
1690  }
1691 
1692  idxtry++;
1693  if(idxtry>size_ghosts-1){
1694  break;
1695  }
1696  Mortontry = ghosts[idxtry].computeMorton();
1697  }
1698  }
1699  }
1700 
1701  uint32_t lengthneigh = 0;
1702  uint32_t sizeneigh = neighbours.size();
1703  for (idxtry=0; idxtry<sizeneigh; idxtry++){
1704  lengthneigh += ghosts[neighbours[idxtry]].getArea();
1705  }
1706  if (lengthneigh < oct->getArea()){
1707  // Search in octants
1708 
1709  // Check if octants face is a boundary
1710  if (oct->info[iface] == false){
1711 
1712  //Build Morton number of virtual neigh of same size
1713  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
1714  Morton = samesizeoct.computeMorton();
1715  // Search morton in octants
1716  // If a even face morton is lower than morton of oct, if odd higher
1717  // ---> can i search only before or after idx in octants
1718  int32_t jump = (oct->computeMorton() > Morton) ? int32_t(idx/2+1) : int32_t((noctants -idx)/2+1);
1719  idxtry = uint32_t(idx +((oct->computeMorton()<Morton)-(oct->computeMorton()>Morton))*jump);
1720  if (idxtry > noctants-1) idxtry = noctants-1;
1721  while(abs(jump) > 0){
1722  Mortontry = octants[idxtry].computeMorton();
1723  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
1724  idxtry += jump;
1725  if (idxtry > noctants-1){
1726  if (jump > 0){
1727  idxtry = noctants - 1;
1728  jump = 0;
1729  }
1730  else if (jump < 0){
1731  idxtry = 0;
1732  jump = 0;
1733  }
1734  }
1735  }
1736  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
1737  //Found neighbour of same size
1738  isghost.push_back(false);
1739  neighbours.push_back(idxtry);
1740  return;
1741  }
1742  else{
1743  // Step until the mortontry lower than morton (one idx of distance)
1744  {
1745  while(octants[idxtry].computeMorton() < Morton){
1746  idxtry++;
1747  if(idxtry > noctants-1){
1748  idxtry = noctants-1;
1749  break;
1750  }
1751  }
1752  while(octants[idxtry].computeMorton() > Morton){
1753  idxtry--;
1754  if(idxtry > noctants-1){
1755  idxtry = 0;
1756  break;
1757  }
1758  }
1759  }
1760  if (idxtry < noctants){
1761  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
1762  //Found neighbour of same size
1763  isghost.push_back(false);
1764  neighbours.push_back(idxtry);
1765  return;
1766  }
1767  // Compute Last discendent of virtual octant of same size
1768  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
1769  uint64_t Mortonlast = last_desc.computeMorton();
1770  Mortontry = octants[idxtry].computeMorton();
1771  int32_t Dx, Dy, Dz;
1772  int32_t Dxstar, Dystar, Dzstar;
1773  while(Mortontry < Mortonlast && idxtry <= noctants-1){
1774  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
1775  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
1776  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
1777  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
1778  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
1779  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
1780 
1781  uint32_t x0 = oct->x;
1782  uint32_t x1 = x0 + size;
1783  uint32_t y0 = oct->y;
1784  uint32_t y1 = y0 + size;
1785  uint32_t z0 = oct->z;
1786  uint32_t z1 = z0 + size;
1787  uint32_t x0try = octants[idxtry].x;
1788  uint32_t x1try = x0try + octants[idxtry].getSize();
1789  uint32_t y0try = octants[idxtry].y;
1790  uint32_t y1try = y0try + octants[idxtry].getSize();
1791  uint32_t z0try = octants[idxtry].z;
1792  uint32_t z1try = z0try + octants[idxtry].getSize();
1793  uint8_t level = oct->level;
1794  uint8_t leveltry = octants[idxtry].getLevel();
1795 
1796  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
1797  if (leveltry > level){
1798  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
1799  neighbours.push_back(idxtry);
1800  isghost.push_back(false);
1801  }
1802  }
1803  if (leveltry < level){
1804  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
1805  neighbours.push_back(idxtry);
1806  isghost.push_back(false);
1807  }
1808  }
1809  }
1810 
1811  idxtry++;
1812  if(idxtry>noctants-1){
1813  break;
1814  }
1815  Mortontry = octants[idxtry].computeMorton();
1816  }
1817  }
1818  }
1819  }
1820  }
1821  return;
1822  }
1823  }
1824  else{
1825  // Boundary Face
1826  return;
1827  }
1828  }
1829  };
1830 
1831  // =================================================================================== //
1832 
1833  void findNeighbours(Class_Octant<3>* oct, // Finds neighbours of octant through iface in vector octants.
1834  uint8_t iface, // Returns a vector (empty if iface is a bound face) with the index of neighbours
1835  u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
1836  vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
1837 
1838  uint64_t Morton, Mortontry;
1839  uint32_t noctants = getNumOctants();
1840  uint32_t idxtry;
1841  uint32_t size = oct->getSize();
1842 
1843  // TODO Create a global matrix
1844  //Alternative to switch case
1845  int8_t cx = int8_t((iface<2)*(int8_t(2*iface-1)));
1846  int8_t cy = int8_t((iface<4)*(int8_t(iface/2))*(int8_t(2*iface-5)));
1847  int8_t cz = int8_t((int8_t(iface/4))*(int8_t(2*iface-9)));
1848 
1849  isghost.clear();
1850  neighbours.clear();
1851 
1852  // Default if iface is nface<iface<0
1853  if (iface < 0 || iface > global3D.nfaces){
1854  return;
1855  }
1856 
1857  // Check if octants face is a process boundary
1858  if (oct->info[global3D.nfaces+iface] == false){
1859 
1860  // Check if octants face is a boundary
1861  if (oct->info[iface] == false){
1862 
1863  //Build Morton number of virtual neigh of same size
1864  Class_Octant<3> samesizeoct(oct->level, int32_t(oct->x)+int32_t(cx*size), int32_t(oct->y)+int32_t(cy*size), int32_t(oct->z)+int32_t(cz*size));
1865  Morton = samesizeoct.computeMorton();
1866  // Search morton in octants
1867  // If a even face morton is lower than morton of oct, if odd higher
1868  // ---> can i search only before or after idx in octants
1869  int32_t jump = int32_t((noctants)/2+1);
1870  idxtry = uint32_t(jump);
1871  Mortontry = oct->computeMorton();
1872  while(abs(jump) > 0){
1873  Mortontry = octants[idxtry].computeMorton();
1874  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
1875  idxtry += jump;
1876  if (idxtry > noctants-1){
1877  if (jump > 0){
1878  idxtry = noctants - 1;
1879  jump = 0;
1880  }
1881  else if (jump < 0){
1882  idxtry = 0;
1883  jump = 0;
1884  }
1885  }
1886  }
1887  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
1888  //Found neighbour of same size
1889  isghost.push_back(false);
1890  neighbours.push_back(idxtry);
1891  return;
1892  }
1893  else{
1894  // Step until the mortontry lower than morton (one idx of distance)
1895  {
1896  while(octants[idxtry].computeMorton() < Morton){
1897  idxtry++;
1898  if(idxtry > noctants-1){
1899  idxtry = noctants-1;
1900  break;
1901  }
1902  }
1903  while(octants[idxtry].computeMorton() > Morton){
1904  idxtry--;
1905  if(idxtry > noctants-1){
1906  idxtry = 0;
1907  break;
1908  }
1909  }
1910  }
1911  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
1912  //Found neighbour of same size
1913  isghost.push_back(false);
1914  neighbours.push_back(idxtry);
1915  return;
1916  }
1917  // Compute Last discendent of virtual octant of same size
1918  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
1919  uint64_t Mortonlast = last_desc.computeMorton();
1920  Mortontry = octants[idxtry].computeMorton();
1921  int32_t Dx, Dy, Dz;
1922  int32_t Dxstar, Dystar, Dzstar;
1923  while(Mortontry < Mortonlast && idxtry < noctants){
1924  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
1925  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
1926  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
1927  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
1928  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
1929  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
1930 
1931  uint32_t x0 = oct->x;
1932  uint32_t x1 = x0 + size;
1933  uint32_t y0 = oct->y;
1934  uint32_t y1 = y0 + size;
1935  uint32_t z0 = oct->z;
1936  uint32_t z1 = z0 + size;
1937  uint32_t x0try = octants[idxtry].x;
1938  uint32_t x1try = x0try + octants[idxtry].getSize();
1939  uint32_t y0try = octants[idxtry].y;
1940  uint32_t y1try = y0try + octants[idxtry].getSize();
1941  uint32_t z0try = octants[idxtry].z;
1942  uint32_t z1try = z0try + octants[idxtry].getSize();
1943  uint8_t level = oct->level;
1944  uint8_t leveltry = octants[idxtry].getLevel();
1945 
1946  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
1947  if (leveltry > level){
1948  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
1949  neighbours.push_back(idxtry);
1950  isghost.push_back(false);
1951  }
1952  }
1953  if (leveltry < level){
1954  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
1955  neighbours.push_back(idxtry);
1956  isghost.push_back(false);
1957  }
1958  }
1959  }
1960 
1961  idxtry++;
1962  if(idxtry>noctants-1){
1963  break;
1964  }
1965  Mortontry = octants[idxtry].computeMorton();
1966  }
1967  return;
1968  }
1969  }
1970  else{
1971  // Boundary Face
1972  return;
1973  }
1974  }
1975  //--------------------------------------------------------------- //
1976  //--------------------------------------------------------------- //
1977  else{
1978  // Check if octants face is a boundary
1979  if (oct->info[iface] == false){
1980  // IF OCTANT FACE IS A PROCESS BOUNDARY SEARCH ALSO IN GHOSTS
1981 
1982  if (ghosts.size()>0){
1983  // Search in ghosts
1984 
1985  uint32_t idxghost = uint32_t(size_ghosts/2);
1986  Class_Octant<3>* octghost = &ghosts[idxghost];
1987 
1988  //Build Morton number of virtual neigh of same size
1989  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
1990  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
1991  // Search morton in octants
1992  // If a even face morton is lower than morton of oct, if odd higher
1993  // ---> can i search only before or after idx in octants
1994  int32_t jump = (octghost->computeMorton() > Morton) ? int32_t(idxghost/2+1) : int32_t((size_ghosts -idxghost)/2+1);
1995  idxtry = uint32_t(idxghost +((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
1996  if (idxtry > ghosts.size()-1) idxtry = ghosts.size()-1;
1997  while(abs(jump) > 0){
1998  Mortontry = ghosts[idxtry].computeMorton();
1999  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
2000  idxtry += jump;
2001  if (idxtry > ghosts.size()-1){
2002  if (jump > 0){
2003  idxtry = ghosts.size() - 1;
2004  jump = 0;
2005  }
2006  else if (jump < 0){
2007  idxtry = 0;
2008  jump = 0;
2009  }
2010  }
2011  }
2012  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
2013  //Found neighbour of same size
2014  isghost.push_back(true);
2015  neighbours.push_back(idxtry);
2016  return;
2017  }
2018  else{
2019  // Step until the mortontry lower than morton (one idx of distance)
2020  {
2021  while(ghosts[idxtry].computeMorton() < Morton){
2022  idxtry++;
2023  if(idxtry > ghosts.size()-1){
2024  idxtry = ghosts.size()-1;
2025  break;
2026  }
2027  }
2028  while(ghosts[idxtry].computeMorton() > Morton){
2029  idxtry--;
2030  if(idxtry > ghosts.size()-1){
2031  idxtry = 0;
2032  break;
2033  }
2034  }
2035  }
2036  if(idxtry < size_ghosts){
2037  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
2038  //Found neighbour of same size
2039  isghost.push_back(true);
2040  neighbours.push_back(idxtry);
2041  return;
2042  }
2043  // Compute Last discendent of virtual octant of same size
2044  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
2045  uint64_t Mortonlast = last_desc.computeMorton();
2046  Mortontry = ghosts[idxtry].computeMorton();
2047  int32_t Dx, Dy, Dz;
2048  int32_t Dxstar, Dystar, Dzstar;
2049  while(Mortontry < Mortonlast && idxtry < size_ghosts){
2050  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(ghosts[idxtry].x));
2051  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(ghosts[idxtry].y));
2052  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(ghosts[idxtry].z));
2053  Dxstar = int32_t((cx-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cx+1)/2)*size;
2054  Dystar = int32_t((cy-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cy+1)/2)*size;
2055  Dzstar = int32_t((cz-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cz+1)/2)*size;
2056 
2057  uint32_t x0 = oct->x;
2058  uint32_t x1 = x0 + size;
2059  uint32_t y0 = oct->y;
2060  uint32_t y1 = y0 + size;
2061  uint32_t z0 = oct->z;
2062  uint32_t z1 = z0 + size;
2063  uint32_t x0try = ghosts[idxtry].x;
2064  uint32_t x1try = x0try + ghosts[idxtry].getSize();
2065  uint32_t y0try = ghosts[idxtry].y;
2066  uint32_t y1try = y0try + ghosts[idxtry].getSize();
2067  uint32_t z0try = ghosts[idxtry].z;
2068  uint32_t z1try = z0try + ghosts[idxtry].getSize();
2069  uint8_t level = oct->level;
2070  uint8_t leveltry = ghosts[idxtry].getLevel();
2071 
2072  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
2073  if (leveltry > level){
2074  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
2075  neighbours.push_back(idxtry);
2076  isghost.push_back(true);
2077  }
2078  }
2079  if (leveltry < level){
2080  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
2081  neighbours.push_back(idxtry);
2082  isghost.push_back(true);
2083  }
2084  }
2085  }
2086 
2087  idxtry++;
2088  if(idxtry>size_ghosts-1){
2089  break;
2090  }
2091  Mortontry = ghosts[idxtry].computeMorton();
2092  }
2093  }
2094  }
2095 
2096  uint32_t lengthneigh = 0;
2097  uint32_t sizeneigh = neighbours.size();
2098  for (idxtry=0; idxtry<sizeneigh; idxtry++){
2099  lengthneigh += ghosts[neighbours[idxtry]].getArea();
2100  }
2101  if (lengthneigh < oct->getArea()){
2102  // Search in octants
2103 
2104  // Check if octants face is a boundary
2105  if (oct->info[iface] == false){
2106 
2107  //Build Morton number of virtual neigh of same size
2108  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
2109  Morton = samesizeoct.computeMorton();
2110  // Search morton in octants
2111  // If a even face morton is lower than morton of oct, if odd higher
2112  // ---> can i search only before or after idx in octants
2113  int32_t jump = (int32_t((noctants)/2+1));
2114  idxtry = uint32_t(jump);
2115  while(abs(jump) > 0){
2116  Mortontry = octants[idxtry].computeMorton();
2117  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
2118  idxtry += jump;
2119  if (idxtry > noctants-1){
2120  if (jump > 0){
2121  idxtry = noctants - 1;
2122  jump = 0;
2123  }
2124  else if (jump < 0){
2125  idxtry = 0;
2126  jump = 0;
2127  }
2128  }
2129  }
2130  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
2131  //Found neighbour of same size
2132  isghost.push_back(false);
2133  neighbours.push_back(idxtry);
2134  return;
2135  }
2136  else{
2137  // Step until the mortontry lower than morton (one idx of distance)
2138  {
2139  while(octants[idxtry].computeMorton() < Morton){
2140  idxtry++;
2141  if(idxtry > noctants-1){
2142  idxtry = noctants-1;
2143  break;
2144  }
2145  }
2146  while(octants[idxtry].computeMorton() > Morton){
2147  idxtry--;
2148  if(idxtry > noctants-1){
2149  idxtry = 0;
2150  break;
2151  }
2152  }
2153  }
2154  if (idxtry < noctants){
2155  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
2156  //Found neighbour of same size
2157  isghost.push_back(false);
2158  neighbours.push_back(idxtry);
2159  return;
2160  }
2161  // Compute Last discendent of virtual octant of same size
2162  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
2163  uint64_t Mortonlast = last_desc.computeMorton();
2164  Mortontry = octants[idxtry].computeMorton();
2165  int32_t Dx, Dy, Dz;
2166  int32_t Dxstar, Dystar, Dzstar;
2167  while(Mortontry < Mortonlast && idxtry <= noctants-1){
2168  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
2169  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
2170  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
2171  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
2172  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
2173  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
2174 
2175  uint32_t x0 = oct->x;
2176  uint32_t x1 = x0 + size;
2177  uint32_t y0 = oct->y;
2178  uint32_t y1 = y0 + size;
2179  uint32_t z0 = oct->z;
2180  uint32_t z1 = z0 + size;
2181  uint32_t x0try = octants[idxtry].x;
2182  uint32_t x1try = x0try + octants[idxtry].getSize();
2183  uint32_t y0try = octants[idxtry].y;
2184  uint32_t y1try = y0try + octants[idxtry].getSize();
2185  uint32_t z0try = octants[idxtry].z;
2186  uint32_t z1try = z0try + octants[idxtry].getSize();
2187  uint8_t level = oct->level;
2188  uint8_t leveltry = octants[idxtry].getLevel();
2189 
2190  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
2191  if (leveltry > level){
2192  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
2193  neighbours.push_back(idxtry);
2194  isghost.push_back(false);
2195  }
2196  }
2197  if (leveltry < level){
2198  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
2199  neighbours.push_back(idxtry);
2200  isghost.push_back(false);
2201  }
2202  }
2203  }
2204 
2205  idxtry++;
2206  Mortontry = octants[idxtry].computeMorton();
2207  }
2208  }
2209  }
2210  }
2211  }
2212  return;
2213  }
2214  }
2215  else{
2216  // Boundary Face
2217  return;
2218  }
2219  }
2220  };
2221 
2222  // =================================================================================== //
2223 
2224  void findGhostNeighbours(uint32_t const idx, // Finds neighbours of idx-th ghost octant through iface in vector octants.
2225  uint8_t iface, // Returns a vector (empty if iface is not the pbound face for ghost) with the index of neighbours
2226  u32vector & neighbours){ // in the structure octants
2227 
2228  uint64_t Morton, Mortontry;
2229  uint32_t noctants = getNumOctants();
2230  uint32_t idxtry;
2231  Class_Octant<3>* oct = &ghosts[idx];
2232  uint32_t size = oct->getSize();
2233 
2234  //Alternative to switch case
2235  int8_t cx = int8_t((iface<2)*(int8_t(2*iface-1)));
2236  int8_t cy = int8_t((iface<4)*(int8_t(iface/2))*(int8_t(2*iface-5)));
2237  int8_t cz = int8_t((int8_t(iface/4))*(int8_t(2*iface-9)));
2238 
2239  neighbours.clear();
2240 
2241  // Default if iface is nface<iface<0
2242  if (iface < 0 || iface > global3D.nfaces){
2243  return;
2244  }
2245 
2246  // Check if octants face is a process boundary
2247  if (oct->info[global3D.nfaces+iface] == true){
2248 
2249  //Build Morton number of virtual neigh of same size
2250  Class_Octant<3> samesizeoct(oct->level, int32_t(oct->x)+int32_t(cx*size), int32_t(oct->y)+int32_t(cy*size), int32_t(oct->z)+int32_t(cz*size));
2251  Morton = samesizeoct.computeMorton();
2252  // Search morton in octants
2253  // If a even face morton is lower than morton of oct, if odd higher
2254  // ---> can i search only before or after idx in octants
2255  int32_t jump = getNumOctants()/2;
2256  idxtry = uint32_t(getNumOctants()/2);
2257  Mortontry = octants[idxtry].computeMorton();
2258  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
2259  while(abs(jump) > 0){
2260  Mortontry = octants[idxtry].computeMorton();
2261  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
2262  idxtry += jump;
2263  if (idxtry > noctants-1){
2264  if (jump > 0){
2265  idxtry = noctants - 1;
2266  jump = 0;
2267  }
2268  else if (jump < 0){
2269  idxtry = 0;
2270  jump = 0;
2271  }
2272  }
2273  }
2274  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
2275  //Found neighbour of same size
2276  neighbours.push_back(idxtry);
2277  return;
2278  }
2279  else{
2280  // Step until the mortontry lower than morton (one idx of distance)
2281  {
2282  while(octants[idxtry].computeMorton() < Morton){
2283  idxtry++;
2284  if(idxtry > noctants-1){
2285  idxtry = noctants-1;
2286  break;
2287  }
2288  }
2289  while(octants[idxtry].computeMorton() > Morton){
2290  idxtry--;
2291  if(idxtry > noctants-1){
2292  idxtry = 0;
2293  break;
2294  }
2295  }
2296  }
2297  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
2298  //Found neighbour of same size
2299  neighbours.push_back(idxtry);
2300  return;
2301  }
2302  // Compute Last discendent of virtual octant of same size
2303  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
2304  uint64_t Mortonlast = last_desc.computeMorton();
2305  Mortontry = octants[idxtry].computeMorton();
2306  int32_t Dx, Dy, Dz;
2307  int32_t Dxstar, Dystar, Dzstar;
2308  while(Mortontry < Mortonlast && idxtry < noctants){
2309  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
2310  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
2311  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
2312  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
2313  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
2314  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
2315 
2316  uint32_t x0 = oct->x;
2317  uint32_t x1 = x0 + size;
2318  uint32_t y0 = oct->y;
2319  uint32_t y1 = y0 + size;
2320  uint32_t z0 = oct->z;
2321  uint32_t z1 = z0 + size;
2322  uint32_t x0try = octants[idxtry].x;
2323  uint32_t x1try = x0try + octants[idxtry].getSize();
2324  uint32_t y0try = octants[idxtry].y;
2325  uint32_t y1try = y0try + octants[idxtry].getSize();
2326  uint32_t z0try = octants[idxtry].z;
2327  uint32_t z1try = z0try + octants[idxtry].getSize();
2328  uint8_t level = oct->level;
2329  uint8_t leveltry = octants[idxtry].getLevel();
2330 
2331  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
2332  if (leveltry > level){
2333  if((abs(cx)*((y0try>=y0)*(y0try<y1))*((z0try>=z0)*(z0try<z1))) + (abs(cy)*((x0try>=x0)*(x0try<x1))*((z0try>=z0)*(z0try<z1))) + (abs(cz)*((x0try>=x0)*(x0try<x1))*((y0try>=y0)*(y0try<y1)))){
2334  neighbours.push_back(idxtry);
2335  }
2336  }
2337  if (leveltry < level){
2338  if((abs(cx)*((y0>=y0try)*(y0<y1try))*((z0>=z0try)*(z0<z1try))) + (abs(cy)*((x0>=x0try)*(x0<x1try))*((z0>=z0try)*(z0<z1try))) + (abs(cz)*((x0>=x0try)*(x0<x1try))*((y0>=y0try)*(y0<y1try)))){
2339  neighbours.push_back(idxtry);
2340  }
2341  }
2342  }
2343 
2344  idxtry++;
2345  if(idxtry>noctants-1){
2346  break;
2347  }
2348  Mortontry = octants[idxtry].computeMorton();
2349  }
2350  return;
2351  }
2352  }
2353  //--------------------------------------------------------------- //
2354  //-----Not Pbound face------------- ----------------------------- //
2355  else{
2356  return;
2357  }
2358 
2359  };
2360 
2361  // =================================================================================== //
2362 
2363  void preBalance21(bool internal){
2364  // Local variables
2365  Class_Octant<3> father, lastdesc;
2366  uint64_t mortonld;
2367  uint32_t nocts;
2368  uint32_t idx, idx2, idx0, last_idx;
2369  uint32_t idx1_gh, idx2_gh;
2370  int8_t markerfather, marker;
2371  uint8_t nbro;
2372  uint8_t nchm1 = global3D.nchildren-1;
2373  bool Bdone=false;
2374 
2375  //------------------------------------------ //
2376  // Initialization
2377 
2378  nbro = 0;
2379  idx2_gh = idx0 = 0;
2380  idx1_gh=0;
2381 
2382  nocts = octants.size();
2383  size_ghosts = ghosts.size();
2384  last_idx=nocts-1;
2385 
2386  //Clean index of ghost brothers in case of coarsening a broken family
2387  last_ghost_bros.clear();
2388 
2389 
2390  // Set index for start and end check for ghosts
2391  if (ghosts.size()){
2392  while(idx2_gh < size_ghosts && ghosts[idx2_gh].computeMorton() <= last_desc.computeMorton()){
2393  idx2_gh++;
2394  }
2395  idx2_gh = min((size_ghosts-1), idx2_gh);
2396 
2397  while(idx1_gh < size_ghosts && ghosts[idx1_gh].computeMorton() <= octants[0].computeMorton()){
2398  idx1_gh++;
2399  }
2400  idx1_gh-=1;
2401  if (idx1_gh==-1) idx1_gh=0;
2402  }
2403 
2404  // End on ghosts
2405  if (ghosts.size() && nocts > 0){
2406  if (ghosts[idx1_gh].buildFather()==octants[0].buildFather()){
2407  father = ghosts[idx1_gh].buildFather();
2408  nbro = 0;
2409  idx = idx1_gh;
2410  marker = ghosts[idx].getMarker();
2411  while(marker < 0 && ghosts[idx].buildFather() == father){
2412  nbro++;
2413  if (idx==0)
2414  break;
2415  idx--;
2416  marker = ghosts[idx].getMarker();
2417  }
2418  idx = 0;
2419  //marker = octants[idx].getMarker();
2420  //while(marker<0 && octants[idx].buildFather() == father){
2421  while(idx<nocts && octants[idx].buildFather() == father){
2422  if(octants[idx].getMarker()<0)
2423  nbro++;
2424  idx++;
2425  if(idx==nocts)
2426  break;
2427  //marker = octants[idx].getMarker();
2428  }
2429  if (nbro != global3D.nchildren && idx!=nocts-1){
2430  for(uint32_t ii=0; ii<idx; ii++){
2431  if (octants[ii].getMarker()<0){
2432  octants[ii].setMarker(0);
2433  octants[ii].info[15]=true;
2434  Bdone=true;
2435  }
2436  }
2437  }
2438  }
2439 
2440  //if ((ghosts[idx2_gh].getMarker() < 0) && (octants[nocts-1].getMarker() < 0)){
2441  if (ghosts[idx2_gh].buildFather()==octants[nocts-1].buildFather()){
2442  father = ghosts[idx2_gh].buildFather();
2443  nbro = 0;
2444  idx = idx2_gh;
2445  marker = ghosts[idx].getMarker();
2446  while(marker < 0 && ghosts[idx].buildFather() == father){
2447 
2448  //Add ghost index to structure for mapper in case of coarsening a broken family
2449  last_ghost_bros.push_back(idx);
2450 
2451  nbro++;
2452  idx++;
2453  if(idx == size_ghosts){
2454  break;
2455  }
2456  marker = ghosts[idx].getMarker();
2457  }
2458  idx = nocts-1;
2459  //marker = octants[idx].getMarker();
2460  //while(marker<0 && octants[idx].buildFather() == father && idx >= 0){
2461  while(octants[idx].buildFather() == father ){
2462  if (octants[idx].getMarker()<0)
2463  nbro++;
2464  //marker = octants[idx].getMarker();
2465  if (idx==0)
2466  break;
2467  idx--;
2468  }
2469  last_idx=idx;
2470  if (nbro != global3D.nchildren && idx!=nocts-1){
2471  for(uint32_t ii=idx+1; ii<nocts; ii++){
2472  if (octants[ii].getMarker()<0){
2473  octants[ii].setMarker(0);
2474  octants[ii].info[15]=true;
2475  Bdone=true;
2476  }
2477  }
2478  //Clean ghost index to structure for mapper in case of coarsening a broken family
2479  last_ghost_bros.clear();
2480  }
2481  }
2482  }
2483 
2484  // Check first internal octants
2485  if (internal){
2486  father = octants[0].buildFather();
2487  lastdesc = father.buildLastDesc();
2488  mortonld = lastdesc.computeMorton();
2489  nbro = 0;
2490  for (idx=0; idx<global3D.nchildren; idx++){
2491  // Check if family is complete or to be checked in the internal loop (some brother refined)
2492  if (octants[idx].computeMorton() <= mortonld){
2493  nbro++;
2494  }
2495  }
2496  if (nbro != global3D.nchildren)
2497  idx0 = nbro;
2498 
2499  // Check and coarse internal octants
2500  for (idx=idx0; idx<nocts; idx++){
2501  if(octants[idx].getMarker() < 0 && octants[idx].getLevel() > 0){
2502  nbro = 0;
2503  father = octants[idx].buildFather();
2504  // Check if family is to be coarsened
2505  for (idx2=idx; idx2<idx+global3D.nchildren; idx2++){
2506  if (idx2<nocts){
2507  if(octants[idx2].getMarker() < 0 && octants[idx2].buildFather() == father){
2508  nbro++;
2509  }
2510  }
2511  }
2512  if (nbro == global3D.nchildren){
2513  idx = idx2-1;
2514  }
2515  else{
2516  if (idx<=last_idx){
2517  octants[idx].setMarker(0);
2518  octants[idx].info[15]=true;
2519  Bdone=true;
2520  }
2521  }
2522  }
2523  }
2524  }
2525  };
2526 
2527  // =================================================================================== //
2528 
2529  void preBalance21(u32vector& newmodified){
2530  // Local variables
2531  Class_Octant<3> father, lastdesc;
2532  uint64_t mortonld;
2533  uint32_t nocts;
2534  uint32_t idx, idx2, idx0, last_idx;
2535  uint32_t idx1_gh, idx2_gh;
2536  int8_t markerfather, marker;
2537  uint8_t nbro;
2538  uint8_t nchm1 = global3D.nchildren-1;
2539  bool Bdone=false;
2540 
2541  //------------------------------------------ //
2542  // Initialization
2543 
2544  nbro = 0;
2545  idx2_gh = idx0 = 0;
2546  idx1_gh=0;
2547 
2548  nocts = octants.size();
2549  size_ghosts = ghosts.size();
2550  last_idx=nocts-1;
2551 
2552  //Clean index of ghost brothers in case of coarsening a broken family
2553  last_ghost_bros.clear();
2554 
2555  // Set index for start and end check for ghosts
2556  if (ghosts.size()){
2557  while(idx2_gh < size_ghosts && ghosts[idx2_gh].computeMorton() <= last_desc.computeMorton()){
2558  idx2_gh++;
2559  }
2560  idx2_gh = min((size_ghosts-1), idx2_gh);
2561 
2562  while(idx1_gh < size_ghosts && ghosts[idx1_gh].computeMorton() <= octants[0].computeMorton()){
2563  idx1_gh++;
2564  }
2565  idx1_gh-=1;
2566  if (idx1_gh==-1) idx1_gh=0;
2567  }
2568 
2569  // End on ghosts
2570  if (ghosts.size() && nocts > 0){
2571  if (ghosts[idx1_gh].buildFather()==octants[0].buildFather()){
2572  father = ghosts[idx1_gh].buildFather();
2573  nbro = 0;
2574  idx = idx1_gh;
2575  marker = ghosts[idx].getMarker();
2576  while(marker < 0 && ghosts[idx].buildFather() == father){
2577 
2578  //Add ghost index to structure for mapper in case of coarsening a broken family
2579  last_ghost_bros.push_back(idx);
2580 
2581  nbro++;
2582  if (idx==0)
2583  break;
2584  idx--;
2585  marker = ghosts[idx].getMarker();
2586  }
2587  idx = 0;
2588  //marker = octants[idx].getMarker();
2589  //while(marker<0 && octants[idx].buildFather() == father){
2590  while(idx<nocts && octants[idx].buildFather() == father){
2591  if (octants[idx].getMarker()<0)
2592  nbro++;
2593  idx++;
2594  if(idx==nocts)
2595  break;
2596  }
2597  if (nbro != global3D.nchildren && idx!=nocts-1){
2598  for(uint32_t ii=0; ii<idx; ii++){
2599  if (octants[ii].getMarker()<0){
2600  octants[ii].setMarker(0);
2601  octants[ii].info[15]=true;
2602  Bdone=true;
2603  newmodified.push_back(ii);
2604  }
2605  }
2606  //Clean index of ghost brothers in case of coarsening a broken family
2607  last_ghost_bros.clear();
2608  }
2609  }
2610 
2611  //if ((ghosts[idx2_gh].getMarker() < 0) && (octants[nocts-1].getMarker() < 0)){
2612  if (ghosts[idx2_gh].buildFather()==octants[nocts-1].buildFather()){
2613  father = ghosts[idx2_gh].buildFather();
2614  nbro = 0;
2615  idx = idx2_gh;
2616  marker = ghosts[idx].getMarker();
2617  while(marker < 0 && ghosts[idx].buildFather() == father){
2618  nbro++;
2619  idx++;
2620  if(idx == size_ghosts){
2621  break;
2622  }
2623  marker = ghosts[idx].getMarker();
2624  }
2625  idx = nocts-1;
2626  //marker = octants[idx].getMarker();
2627  //while(marker<0 && octants[idx].buildFather() == father && idx >= 0){
2628  while(octants[idx].buildFather() == father){
2629  if (octants[idx].getMarker()<0)
2630  nbro++;
2631  idx--;
2632  if (idx==0)
2633  break;
2634  }
2635  last_idx=idx;
2636  if (nbro != global3D.nchildren && idx!=nocts-1){
2637  for(uint32_t ii=idx+1; ii<nocts; ii++){
2638  if (octants[ii].getMarker()<0){
2639  octants[ii].setMarker(0);
2640  octants[ii].info[15]=true;
2641  Bdone=true;
2642  newmodified.push_back(ii);
2643  }
2644  }
2645  }
2646  }
2647  }
2648 
2649  // Check first internal octants
2650  father = octants[0].buildFather();
2651  lastdesc = father.buildLastDesc();
2652  mortonld = lastdesc.computeMorton();
2653  nbro = 0;
2654  for (idx=0; idx<global3D.nchildren; idx++){
2655  // Check if family is complete or to be checked in the internal loop (some brother refined)
2656  if (octants[idx].computeMorton() <= mortonld){
2657  nbro++;
2658  }
2659  }
2660  if (nbro != global3D.nchildren)
2661  idx0 = nbro;
2662 
2663  // Check and coarse internal octants
2664  for (idx=idx0; idx<nocts; idx++){
2665  if(octants[idx].getMarker() < 0 && octants[idx].getLevel() > 0){
2666  nbro = 0;
2667  father = octants[idx].buildFather();
2668  // Check if family is to be coarsened
2669  for (idx2=idx; idx2<idx+global3D.nchildren; idx2++){
2670  if (idx2<nocts){
2671  if(octants[idx2].getMarker() < 0 && octants[idx2].buildFather() == father){
2672  nbro++;
2673  }
2674  }
2675  }
2676  if (nbro == global3D.nchildren){
2677  idx = idx2-1;
2678  }
2679  else{
2680  if (idx<=last_idx){
2681  octants[idx].setMarker(0);
2682  octants[idx].info[15]=true;
2683  Bdone=true;
2684  newmodified.push_back(idx);
2685  }
2686  }
2687  }
2688  }
2689  };
2690 
2691  // =================================================================================== //
2692 
2693  bool localBalance(bool doInterior){ // 2:1 balancing on level a local tree already adapted (balance only the octants with info[14] = false) (refinement wins!)
2694  // Return true if balanced done with some markers modification
2695  // Seto doInterior = false if the interior octants are already balanced
2696  // Local variables
2697  uint32_t sizeneigh, modsize;
2698  u32vector neigh;
2699  u32vector modified, newmodified;
2700  uint32_t i, idx;
2701  uint8_t iface, iedge, inode;
2702  int8_t targetmarker;
2703  vector<bool> isghost;
2704  bool Bdone = false;
2705 
2706  OctantsType::iterator obegin, oend, it;
2707  u32vector::iterator ibegin, iend, iit;
2708 
2709  //If interior octants have to be balanced
2710  if(doInterior){
2711  // First loop on the octants
2712  obegin = octants.begin();
2713  oend = octants.end();
2714  idx = 0;
2715  for (it=obegin; it!=oend; it++){
2716  if (!it->getNotBalance() && it->getMarker() != 0){
2717  targetmarker = min(MAX_LEVEL_3D, (octants[idx].getLevel() + octants[idx].getMarker()));
2718 
2719  //Balance through faces
2720  for (iface=0; iface<global3D.nfaces; iface++){
2721  if(!it->getBound(iface)){
2722  findNeighbours(idx, iface, neigh, isghost);
2723  sizeneigh = neigh.size();
2724  for(i=0; i<sizeneigh; i++){
2725  if (!isghost[i]){
2726  {
2727  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1) ){
2728  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-1-octants[idx].getLevel());
2729  octants[idx].info[15] = true;
2730  modified.push_back(idx);
2731  Bdone = true;
2732  }
2733  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2734  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2735  octants[neigh[i]].info[15] = true;
2736  modified.push_back(neigh[i]);
2737  Bdone = true;
2738  }
2739  };
2740  }
2741  else{
2742  {
2743  if((ghosts[neigh[i]].getLevel() + ghosts[neigh[i]].getMarker()) > (targetmarker + 1) ){
2744  octants[idx].setMarker(ghosts[neigh[i]].getLevel()+ghosts[neigh[i]].getMarker()-1-octants[idx].getLevel());
2745  octants[idx].info[15] = true;
2746  modified.push_back(idx);
2747  Bdone = true;
2748  }
2749  };
2750  }
2751  }
2752  }
2753  }
2754 
2755  if (balance_codim>1){
2756  //Balance through edges
2757  for (iedge=0; iedge<global3D.nedges; iedge++){
2758  if(!it->getBound(global3D.edgeface[iedge][0]) && !it->getBound(global3D.edgeface[iedge][1])){
2759  findEdgeNeighbours(idx, iedge, neigh, isghost);
2760  sizeneigh = neigh.size();
2761  for(i=0; i<sizeneigh; i++){
2762  if (!isghost[i]){
2763  {
2764  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1) ){
2765  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-1-octants[idx].getLevel());
2766  octants[idx].info[15] = true;
2767  modified.push_back(idx);
2768  Bdone = true;
2769  }
2770  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2771  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2772  octants[neigh[i]].info[15] = true;
2773  modified.push_back(neigh[i]);
2774  Bdone = true;
2775  }
2776  };
2777  }
2778  else{
2779  if((ghosts[neigh[i]].getLevel() + ghosts[neigh[i]].getMarker()) > (targetmarker + 1) ){
2780  octants[idx].setMarker(ghosts[neigh[i]].getLevel()+ghosts[neigh[i]].getMarker()-1-octants[idx].getLevel());
2781  octants[idx].info[15] = true;
2782  modified.push_back(idx);
2783  Bdone = true;
2784  }
2785  }
2786  }
2787  }
2788  }
2789  }
2790 
2791  if (balance_codim>2){
2792  //Balance through nodes
2793  for (inode=0; inode<global3D.nnodes; inode++){
2794  if(!it->getBound(global3D.nodeface[inode][0]) && !it->getBound(global3D.nodeface[inode][1]) && !it->getBound(global3D.nodeface[inode][2])){
2795  findNodeNeighbours(idx, inode, neigh, isghost);
2796  sizeneigh = neigh.size();
2797  for(i=0; i<sizeneigh; i++){
2798  if (!isghost[i]){
2799  {
2800  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1) ){
2801  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-1-octants[idx].getLevel());
2802  octants[idx].info[15] = true;
2803  modified.push_back(idx);
2804  Bdone = true;
2805  }
2806  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2807  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2808  octants[neigh[i]].info[15] = true;
2809  modified.push_back(neigh[i]);
2810  Bdone = true;
2811  }
2812  };
2813  }
2814  else{
2815  if((ghosts[neigh[i]].getLevel() + ghosts[neigh[i]].getMarker()) > (targetmarker + 1) ){
2816  octants[idx].setMarker(ghosts[neigh[i]].getLevel()+ghosts[neigh[i]].getMarker()-1-octants[idx].getLevel());
2817  octants[idx].info[15] = true;
2818  modified.push_back(idx);
2819  Bdone = true;
2820  }
2821  }
2822  }
2823  }
2824  }
2825  }
2826 
2827  }
2828  idx++;
2829  }
2830  // Loop on ghost octants (influence over interior borders)
2831  obegin = ghosts.begin();
2832  oend = ghosts.end();
2833  idx = 0;
2834  for (it=obegin; it!=oend; it++){
2835  if (!it->getNotBalance() && it->getMarker() != 0){
2836  targetmarker = min(MAX_LEVEL_3D, (it->getLevel()+it->getMarker()));
2837 
2838  //Balance through faces
2839  for (iface=0; iface<global3D.nfaces; iface++){
2840  if(it->getPbound(iface) == true){
2841  neigh.clear();
2842  findGhostNeighbours(idx, iface, neigh);
2843  sizeneigh = neigh.size();
2844  for(i=0; i<sizeneigh; i++){
2845  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2846  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2847  octants[neigh[i]].info[15] = true;
2848  modified.push_back(neigh[i]);
2849  Bdone = true;
2850  }
2851  }
2852  }
2853  }
2854 
2855  if (balance_codim>1){
2856  //Balance through edges
2857  for (iedge=0; iedge<global3D.nedges; iedge++){
2858  if(it->getPbound(global3D.nodeface[iedge][0]) == true || it->getPbound(global3D.nodeface[iedge][1]) == true){
2859  neigh.clear();
2860  findGhostEdgeNeighbours(idx, iedge, neigh);
2861  sizeneigh = neigh.size();
2862  for(i=0; i<sizeneigh; i++){
2863  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2864  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2865  octants[neigh[i]].info[15] = true;
2866  modified.push_back(neigh[i]);
2867  Bdone = true;
2868  }
2869  }
2870  }
2871  }
2872  }
2873 
2874  if (balance_codim>2){
2875  //Balance through nodes
2876  for (inode=0; inode<global3D.nnodes; inode++){
2877  if(it->getPbound(global3D.nodeface[inode][0]) == true || it->getPbound(global3D.nodeface[inode][1]) == true || it->getPbound(global3D.nodeface[inode][2]) == true){
2878  neigh.clear();
2879  findGhostNodeNeighbours(idx, inode, neigh);
2880  sizeneigh = neigh.size();
2881  for(i=0; i<sizeneigh; i++){
2882  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2883  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2884  octants[neigh[i]].info[15] = true;
2885  modified.push_back(neigh[i]);
2886  Bdone = true;
2887  }
2888  }
2889  }
2890  }
2891  }
2892 
2893  }
2894  idx++;
2895  }
2896 
2897  // While loop for iterative balancing
2898  u32vector().swap(newmodified);
2899  modsize = modified.size();
2900  while(modsize!=0){
2901  ibegin = modified.begin();
2902  iend = modified.end();
2903  for (iit=ibegin; iit!=iend; iit++){
2904  idx = *iit;
2905  if (!octants[idx].getNotBalance()){
2906  targetmarker = min(MAX_LEVEL_3D, (octants[idx].getLevel()+octants[idx].getMarker()));
2907 
2908  //Balance through faces
2909  for (iface=0; iface<global3D.nfaces; iface++){
2910  if(!octants[idx].getPbound(iface)){
2911  findNeighbours(idx, iface, neigh, isghost);
2912  sizeneigh = neigh.size();
2913  for(i=0; i<sizeneigh; i++){
2914  if (!isghost[i]){
2915  {
2916  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
2917  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
2918  octants[idx].info[15] = true;
2919  newmodified.push_back(idx);
2920  Bdone = true;
2921  }
2922  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2923  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2924  octants[neigh[i]].info[15] = true;
2925  newmodified.push_back(neigh[i]);
2926  Bdone = true;
2927  }
2928  };
2929  }
2930  }
2931  }
2932  }
2933 
2934  if (balance_codim>1){
2935  //Balance through edges
2936  for (iedge=0; iedge<global3D.nedges; iedge++){
2937  if(!octants[idx].getPbound(global3D.edgeface[iedge][0]) || !octants[idx].getPbound(global3D.edgeface[iedge][1])){
2938  findEdgeNeighbours(idx, iedge, neigh, isghost);
2939  sizeneigh = neigh.size();
2940  for(i=0; i<sizeneigh; i++){
2941  if (!isghost[i]){
2942  {
2943  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
2944  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
2945  octants[idx].info[15] = true;
2946  newmodified.push_back(idx);
2947  Bdone = true;
2948  }
2949  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2950  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2951  octants[neigh[i]].info[15] = true;
2952  newmodified.push_back(neigh[i]);
2953  Bdone = true;
2954  }
2955  };
2956  }
2957  }
2958  }
2959  }
2960  }
2961 
2962  if (balance_codim>2){
2963  //Balance through nodes
2964  for (inode=0; inode<global3D.nnodes; inode++){
2965  if(!octants[idx].getPbound(global3D.nodeface[inode][0]) || !octants[idx].getPbound(global3D.nodeface[inode][1]) || !octants[idx].getPbound(global3D.nodeface[inode][2])){
2966  findNodeNeighbours(idx, inode, neigh, isghost);
2967  sizeneigh = neigh.size();
2968  for(i=0; i<sizeneigh; i++){
2969  if (!isghost[i]){
2970  {
2971  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
2972  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
2973  octants[idx].info[15] = true;
2974  newmodified.push_back(idx);
2975  Bdone = true;
2976  }
2977  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
2978  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
2979  octants[neigh[i]].info[15] = true;
2980  newmodified.push_back(neigh[i]);
2981  Bdone = true;
2982  }
2983  };
2984  }
2985  }
2986  }
2987  }
2988  }
2989 
2990  }
2991  }
2992  preBalance21(newmodified);
2993  u32vector().swap(modified);
2994  swap(modified,newmodified);
2995  modsize = modified.size();
2996  u32vector().swap(newmodified);
2997  }// end while
2998 
2999  }
3000  else{
3001 
3002  // Loop on ghost octants (influence over interior borders)
3003  obegin = ghosts.begin();
3004  oend = ghosts.end();
3005  idx = 0;
3006  for (it=obegin; it!=oend; it++){
3007  if (!it->getNotBalance() && it->info[15]){
3008  targetmarker = min(MAX_LEVEL_3D, (it->getLevel()+it->getMarker()));
3009 
3010  //Balance through faces
3011  for (iface=0; iface<global3D.nfaces; iface++){
3012  if(it->getPbound(iface) == true){
3013  neigh.clear();
3014  findGhostNeighbours(idx, iface, neigh);
3015  sizeneigh = neigh.size();
3016  for(i=0; i<sizeneigh; i++){
3017  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3018  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3019  octants[neigh[i]].info[15] = true;
3020  modified.push_back(neigh[i]);
3021  Bdone = true;
3022  }
3023  }
3024  }
3025  }
3026 
3027  if (balance_codim>1){
3028  //Balance through edges
3029  for (iedge=0; iedge<global3D.nedges; iedge++){
3030  if(it->getPbound(global3D.nodeface[iedge][0]) == true || it->getPbound(global3D.nodeface[iedge][1]) == true){
3031  neigh.clear();
3032  findGhostEdgeNeighbours(idx, iedge, neigh);
3033  sizeneigh = neigh.size();
3034  for(i=0; i<sizeneigh; i++){
3035  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3036  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3037  octants[neigh[i]].info[15] = true;
3038  modified.push_back(neigh[i]);
3039  Bdone = true;
3040  }
3041  }
3042  }
3043  }
3044  }
3045 
3046  if (balance_codim>2){
3047  //Balance through nodes
3048  for (inode=0; inode<global3D.nnodes; inode++){
3049  if(it->getPbound(global3D.nodeface[inode][0]) == true || it->getPbound(global3D.nodeface[inode][1]) == true || it->getPbound(global3D.nodeface[inode][2]) == true){
3050  neigh.clear();
3051  findGhostNodeNeighbours(idx, inode, neigh);
3052  sizeneigh = neigh.size();
3053  for(i=0; i<sizeneigh; i++){
3054  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3055  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3056  octants[neigh[i]].info[15] = true;
3057  modified.push_back(neigh[i]);
3058  Bdone = true;
3059  }
3060  }
3061  }
3062  }
3063  }
3064 
3065  }
3066  idx++;
3067  }
3068 
3069  // While loop for iterative balancing
3070  u32vector().swap(newmodified);
3071  modsize = modified.size();
3072  while(modsize!=0){
3073  ibegin = modified.begin();
3074  iend = modified.end();
3075  for (iit=ibegin; iit!=iend; iit++){
3076  idx = *iit;
3077  if (!octants[idx].getNotBalance()){
3078  targetmarker = min(MAX_LEVEL_3D, (octants[idx].getLevel()+octants[idx].getMarker()));
3079 
3080  //Balance through faces
3081  for (iface=0; iface<global3D.nfaces; iface++){
3082  if(!octants[idx].getPbound(iface)){
3083  findNeighbours(idx, iface, neigh, isghost);
3084  sizeneigh = neigh.size();
3085  for(i=0; i<sizeneigh; i++){
3086  if (!isghost[i]){
3087  {
3088  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3089  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3090  octants[idx].info[15] = true;
3091  newmodified.push_back(idx);
3092  Bdone = true;
3093  }
3094  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3095  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3096  octants[neigh[i]].info[15] = true;
3097  newmodified.push_back(neigh[i]);
3098  Bdone = true;
3099  }
3100  };
3101  }
3102  }
3103  }
3104  }
3105 
3106  if (balance_codim>1){
3107  //Balance through edges
3108  for (iedge=0; iedge<global3D.nedges; iedge++){
3109  if(!octants[idx].getPbound(global3D.edgeface[iedge][0]) || !octants[idx].getPbound(global3D.edgeface[iedge][1])){
3110  findEdgeNeighbours(idx, iedge, neigh, isghost);
3111  sizeneigh = neigh.size();
3112  for(i=0; i<sizeneigh; i++){
3113  if (!isghost[i]){
3114  {
3115  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3116  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3117  octants[idx].info[15] = true;
3118  newmodified.push_back(idx);
3119  Bdone = true;
3120  }
3121  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3122  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3123  octants[neigh[i]].info[15] = true;
3124  newmodified.push_back(neigh[i]);
3125  Bdone = true;
3126  }
3127  };
3128  }
3129  }
3130  }
3131  }
3132  }
3133 
3134  if (balance_codim>2){
3135  //Balance through nodes
3136  for (inode=0; inode<global3D.nnodes; inode++){
3137  if(!octants[idx].getPbound(global3D.nodeface[inode][0]) || !octants[idx].getPbound(global3D.nodeface[inode][1]) || !octants[idx].getPbound(global3D.nodeface[inode][2])){
3138  findNodeNeighbours(idx, inode, neigh, isghost);
3139  sizeneigh = neigh.size();
3140  for(i=0; i<sizeneigh; i++){
3141  if (!isghost[i]){
3142  {
3143  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3144  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3145  octants[idx].info[15] = true;
3146  newmodified.push_back(idx);
3147  Bdone = true;
3148  }
3149  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3150  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3151  octants[neigh[i]].info[15] = true;
3152  newmodified.push_back(neigh[i]);
3153  Bdone = true;
3154  }
3155  };
3156  }
3157  }
3158  }
3159  }
3160  }
3161 
3162  }
3163  }
3164  preBalance21(newmodified);
3165  u32vector().swap(modified);
3166  swap(modified,newmodified);
3167  modsize = modified.size();
3168  u32vector().swap(newmodified);
3169  }// end while
3170  obegin = oend = octants.end();
3171  ibegin = iend = modified.end();
3172  }
3173  return Bdone;
3174  // Pay attention : info[15] may be true after local balance for some octants
3175  };
3176 
3177  // =================================================================================== //
3178 
3179  bool localBalanceAll(bool doInterior){ // 2:1 balancing on level a local tree already adapted (balance only the octants with info[14] = false) (refinement wins!)
3180  // Return true if balanced done with some markers modification
3181  // Seto doInterior = false if the interior octants are already balanced
3182  // Local variables
3183  uint32_t sizeneigh, modsize;
3184  u32vector neigh;
3185  u32vector modified, newmodified;
3186  uint32_t i, idx;
3187  uint8_t iface, iedge, inode;
3188  int8_t targetmarker;
3189  vector<bool> isghost;
3190  bool Bdone = false;
3191 
3192  OctantsType::iterator obegin, oend, it;
3193  u32vector::iterator ibegin, iend, iit;
3194 
3195 
3196  //If interior octants have to be balanced
3197  if(doInterior){
3198  // First loop on the octants
3199  obegin = octants.begin();
3200  oend = octants.end();
3201  idx = 0;
3202  for (it=obegin; it!=oend; it++){
3203  if ((!it->getNotBalance()) && ((it->info[15]) || (it->getMarker()!=0) || ((it->getIsNewC()) || (it->getIsNewR())))){
3204  targetmarker = min(MAX_LEVEL_3D, int(octants[idx].getLevel()) + int(octants[idx].getMarker()));
3205 
3206  //Balance through faces
3207  for (iface=0; iface<global3D.nfaces; iface++){
3208  if(!it->getBound(iface)){
3209  findNeighbours(idx, iface, neigh, isghost);
3210  sizeneigh = neigh.size();
3211  for(i=0; i<sizeneigh; i++){
3212  if (!isghost[i]){
3213  {
3214  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1) ){
3215  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-1-octants[idx].getLevel());
3216  octants[idx].info[15] = true;
3217  modified.push_back(idx);
3218  Bdone = true;
3219  }
3220  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3221  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3222  octants[neigh[i]].info[15] = true;
3223  modified.push_back(neigh[i]);
3224  Bdone = true;
3225  }
3226  };
3227  }
3228  else{
3229  {
3230  if((ghosts[neigh[i]].getLevel() + ghosts[neigh[i]].getMarker()) > (targetmarker + 1) ){
3231  octants[idx].setMarker(ghosts[neigh[i]].getLevel()+ghosts[neigh[i]].getMarker()-1-octants[idx].getLevel());
3232  octants[idx].info[15] = true;
3233  modified.push_back(idx);
3234  Bdone = true;
3235  }
3236  };
3237 
3238  }
3239  }
3240  }
3241  }
3242 
3243  if (balance_codim>1){
3244  //Balance through edges
3245  for (iedge=0; iedge<global3D.nedges; iedge++){
3246  if(!it->getBound(global3D.edgeface[iedge][0]) && !it->getBound(global3D.edgeface[iedge][1])){
3247  findEdgeNeighbours(idx, iedge, neigh, isghost);
3248  sizeneigh = neigh.size();
3249  for(i=0; i<sizeneigh; i++){
3250  if (!isghost[i]){
3251  {
3252  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1) ){
3253  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-1-octants[idx].getLevel());
3254  octants[idx].info[15] = true;
3255  modified.push_back(idx);
3256  Bdone = true;
3257  }
3258  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3259  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3260  octants[neigh[i]].info[15] = true;
3261  modified.push_back(neigh[i]);
3262  Bdone = true;
3263  }
3264  };
3265  }
3266  else{
3267  if((ghosts[neigh[i]].getLevel() + ghosts[neigh[i]].getMarker()) > (targetmarker + 1) ){
3268  octants[idx].setMarker(ghosts[neigh[i]].getLevel()+ghosts[neigh[i]].getMarker()-1-octants[idx].getLevel());
3269  octants[idx].info[15] = true;
3270  modified.push_back(idx);
3271  Bdone = true;
3272  }
3273  }
3274  }
3275  }
3276  }
3277  }
3278 
3279  if (balance_codim>2){
3280  //Balance through nodes
3281  for (inode=0; inode<global3D.nnodes; inode++){
3282  if(!it->getBound(global3D.nodeface[inode][0]) && !it->getBound(global3D.nodeface[inode][1]) && !it->getBound(global3D.nodeface[inode][2])){
3283  findNodeNeighbours(idx, inode, neigh, isghost);
3284  sizeneigh = neigh.size();
3285  for(i=0; i<sizeneigh; i++){
3286  if (!isghost[i]){
3287  {
3288  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1) ){
3289  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-1-octants[idx].getLevel());
3290  octants[idx].info[15] = true;
3291  modified.push_back(idx);
3292  Bdone = true;
3293  }
3294  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3295  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3296  octants[neigh[i]].info[15] = true;
3297  modified.push_back(neigh[i]);
3298  Bdone = true;
3299  }
3300  };
3301  }
3302  else{
3303  if((ghosts[neigh[i]].getLevel() + ghosts[neigh[i]].getMarker()) > (targetmarker + 1) ){
3304  octants[idx].setMarker(ghosts[neigh[i]].getLevel()+ghosts[neigh[i]].getMarker()-1-octants[idx].getLevel());
3305  octants[idx].info[15] = true;
3306  modified.push_back(idx);
3307  Bdone = true;
3308  }
3309  }
3310  }
3311  }
3312  }
3313  }
3314 
3315  }
3316  idx++;
3317  }
3318  // Loop on ghost octants (influence over interior borders)
3319  obegin = ghosts.begin();
3320  oend = ghosts.end();
3321  idx = 0;
3322  for (it=obegin; it!=oend; it++){
3323  if (!it->getNotBalance() && (it->info[15] || (it->getIsNewC() || it->getIsNewR()))){
3324  targetmarker = min(MAX_LEVEL_3D, (it->getLevel()+it->getMarker()));
3325 
3326  //Balance through faces
3327  for (iface=0; iface<global3D.nfaces; iface++){
3328  if(it->getPbound(iface) == true){
3329  neigh.clear();
3330  findGhostNeighbours(idx, iface, neigh);
3331  sizeneigh = neigh.size();
3332  for(i=0; i<sizeneigh; i++){
3333  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3334  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3335  octants[neigh[i]].info[15] = true;
3336  modified.push_back(neigh[i]);
3337  Bdone = true;
3338  }
3339  }
3340  }
3341  }
3342 
3343  if (balance_codim>1){
3344  //Balance through edges
3345  for (iedge=0; iedge<global3D.nedges; iedge++){
3346  if(it->getPbound(global3D.nodeface[iedge][0]) == true || it->getPbound(global3D.nodeface[iedge][1]) == true){
3347  neigh.clear();
3348  findGhostEdgeNeighbours(idx, iedge, neigh);
3349  sizeneigh = neigh.size();
3350  for(i=0; i<sizeneigh; i++){
3351  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3352  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3353  octants[neigh[i]].info[15] = true;
3354  modified.push_back(neigh[i]);
3355  Bdone = true;
3356  }
3357  }
3358  }
3359  }
3360  }
3361 
3362  if (balance_codim>2){
3363  //Balance through nodes
3364  for (inode=0; inode<global3D.nnodes; inode++){
3365  if(it->getPbound(global3D.nodeface[inode][0]) == true || it->getPbound(global3D.nodeface[inode][1]) == true || it->getPbound(global3D.nodeface[inode][2]) == true){
3366  neigh.clear();
3367  findGhostNodeNeighbours(idx, inode, neigh);
3368  sizeneigh = neigh.size();
3369  for(i=0; i<sizeneigh; i++){
3370  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3371  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3372  octants[neigh[i]].info[15] = true;
3373  modified.push_back(neigh[i]);
3374  Bdone = true;
3375  }
3376  }
3377  }
3378  }
3379  }
3380 
3381  }
3382  idx++;
3383  }
3384 
3385  // While loop for iterative balancing
3386  u32vector().swap(newmodified);
3387  modsize = modified.size();
3388  while(modsize!=0){
3389  ibegin = modified.begin();
3390  iend = modified.end();
3391  for (iit=ibegin; iit!=iend; iit++){
3392  idx = *iit;
3393  if (!octants[idx].getNotBalance()){
3394  targetmarker = min(MAX_LEVEL_3D, (octants[idx].getLevel()+octants[idx].getMarker()));
3395 
3396  //Balance through faces
3397  for (iface=0; iface<global3D.nfaces; iface++){
3398  if(!octants[idx].getPbound(iface)){
3399  findNeighbours(idx, iface, neigh, isghost);
3400  sizeneigh = neigh.size();
3401  for(i=0; i<sizeneigh; i++){
3402  if (!isghost[i]){
3403  {
3404  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3405  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3406  octants[idx].info[15] = true;
3407  newmodified.push_back(idx);
3408  Bdone = true;
3409  }
3410  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3411  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3412  octants[neigh[i]].info[15] = true;
3413  newmodified.push_back(neigh[i]);
3414  Bdone = true;
3415  }
3416  };
3417  }
3418  }
3419  }
3420  }
3421 
3422  if (balance_codim>1){
3423  //Balance through edges
3424  for (iedge=0; iedge<global3D.nedges; iedge++){
3425  if(!octants[idx].getPbound(global3D.edgeface[iedge][0]) || !octants[idx].getPbound(global3D.edgeface[iedge][1])){
3426  findEdgeNeighbours(idx, iedge, neigh, isghost);
3427  sizeneigh = neigh.size();
3428  for(i=0; i<sizeneigh; i++){
3429  if (!isghost[i]){
3430  {
3431  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3432  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3433  octants[idx].info[15] = true;
3434  newmodified.push_back(idx);
3435  Bdone = true;
3436  }
3437  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3438  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3439  octants[neigh[i]].info[15] = true;
3440  newmodified.push_back(neigh[i]);
3441  Bdone = true;
3442  }
3443  };
3444  }
3445  }
3446  }
3447  }
3448  }
3449 
3450  if (balance_codim>2){
3451  //Balance through nodes
3452  for (inode=0; inode<global3D.nnodes; inode++){
3453  if(!octants[idx].getPbound(global3D.nodeface[inode][0]) || !octants[idx].getPbound(global3D.nodeface[inode][1]) || !octants[idx].getPbound(global3D.nodeface[inode][2])){
3454  findNodeNeighbours(idx, inode, neigh, isghost);
3455  sizeneigh = neigh.size();
3456  for(i=0; i<sizeneigh; i++){
3457  if (!isghost[i]){
3458  {
3459  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3460  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3461  octants[idx].info[15] = true;
3462  newmodified.push_back(idx);
3463  Bdone = true;
3464  }
3465  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3466  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3467  octants[neigh[i]].info[15] = true;
3468  newmodified.push_back(neigh[i]);
3469  Bdone = true;
3470  }
3471  };
3472  }
3473  }
3474  }
3475  }
3476  }
3477 
3478  }
3479  }
3480  preBalance21(newmodified);
3481  u32vector().swap(modified);
3482  swap(modified,newmodified);
3483  modsize = modified.size();
3484  u32vector().swap(newmodified);
3485  }// end while
3486 
3487  }
3488  else{
3489 
3490  // Loop on ghost octants (influence over interior borders)
3491  obegin = ghosts.begin();
3492  oend = ghosts.end();
3493  idx = 0;
3494  for (it=obegin; it!=oend; it++){
3495  if (!it->getNotBalance() && (it->info[15] || (it->getIsNewC() || it->getIsNewR()))){
3496  targetmarker = min(MAX_LEVEL_3D, (it->getLevel()+it->getMarker()));
3497 
3498  //Balance through faces
3499  for (iface=0; iface<global3D.nfaces; iface++){
3500  if(it->getPbound(iface) == true){
3501  neigh.clear();
3502  findGhostNeighbours(idx, iface, neigh);
3503  sizeneigh = neigh.size();
3504  for(i=0; i<sizeneigh; i++){
3505  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3506  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3507  octants[neigh[i]].info[15] = true;
3508  modified.push_back(neigh[i]);
3509  Bdone = true;
3510  }
3511  }
3512  }
3513  }
3514 
3515  if (balance_codim>1){
3516  //Balance through edges
3517  for (iedge=0; iedge<global3D.nedges; iedge++){
3518  if(it->getPbound(global3D.nodeface[iedge][0]) == true || it->getPbound(global3D.nodeface[iedge][1]) == true){
3519  neigh.clear();
3520  findGhostEdgeNeighbours(idx, iedge, neigh);
3521  sizeneigh = neigh.size();
3522  for(i=0; i<sizeneigh; i++){
3523  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3524  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3525  octants[neigh[i]].info[15] = true;
3526  modified.push_back(neigh[i]);
3527  Bdone = true;
3528  }
3529  }
3530  }
3531  }
3532  }
3533 
3534  if (balance_codim>2){
3535  //Balance through nodes
3536  for (inode=0; inode<global3D.nnodes; inode++){
3537  if(it->getPbound(global3D.nodeface[inode][0]) == true || it->getPbound(global3D.nodeface[inode][1]) == true || it->getPbound(global3D.nodeface[inode][2]) == true){
3538  neigh.clear();
3539  findGhostNodeNeighbours(idx, inode, neigh);
3540  sizeneigh = neigh.size();
3541  for(i=0; i<sizeneigh; i++){
3542  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3543  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3544  octants[neigh[i]].info[15] = true;
3545  modified.push_back(neigh[i]);
3546  Bdone = true;
3547  }
3548  }
3549  }
3550  }
3551  }
3552 
3553  }
3554  idx++;
3555  }
3556 
3557  // While loop for iterative balancing
3558  u32vector().swap(newmodified);
3559  modsize = modified.size();
3560  while(modsize!=0){
3561  ibegin = modified.begin();
3562  iend = modified.end();
3563  for (iit=ibegin; iit!=iend; iit++){
3564  idx = *iit;
3565  if (!octants[idx].getNotBalance()){
3566  targetmarker = min(MAX_LEVEL_3D, (octants[idx].getLevel()+octants[idx].getMarker()));
3567 
3568  //Balance through faces
3569  for (iface=0; iface<global3D.nfaces; iface++){
3570  if(!octants[idx].getPbound(iface)){
3571  findNeighbours(idx, iface, neigh, isghost);
3572  sizeneigh = neigh.size();
3573  for(i=0; i<sizeneigh; i++){
3574  if (!isghost[i]){
3575  {
3576  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3577  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3578  octants[idx].info[15] = true;
3579  newmodified.push_back(idx);
3580  Bdone = true;
3581  }
3582  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3583  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3584  octants[neigh[i]].info[15] = true;
3585  newmodified.push_back(neigh[i]);
3586  Bdone = true;
3587  }
3588  };
3589  }
3590  }
3591  }
3592  }
3593 
3594  if (balance_codim>1){
3595  //Balance through edges
3596  for (iedge=0; iedge<global3D.nedges; iedge++){
3597  if(!octants[idx].getPbound(global3D.edgeface[iedge][0]) || !octants[idx].getPbound(global3D.edgeface[iedge][1])){
3598  findEdgeNeighbours(idx, iedge, neigh, isghost);
3599  sizeneigh = neigh.size();
3600  for(i=0; i<sizeneigh; i++){
3601  if (!isghost[i]){
3602  {
3603  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3604  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3605  octants[idx].info[15] = true;
3606  newmodified.push_back(idx);
3607  Bdone = true;
3608  }
3609  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3610  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3611  octants[neigh[i]].info[15] = true;
3612  newmodified.push_back(neigh[i]);
3613  Bdone = true;
3614  }
3615  };
3616  }
3617  }
3618  }
3619  }
3620  }
3621 
3622  if (balance_codim>2){
3623  //Balance through nodes
3624  for (inode=0; inode<global3D.nnodes; inode++){
3625  if(!octants[idx].getPbound(global3D.nodeface[inode][0]) || !octants[idx].getPbound(global3D.nodeface[inode][1]) || !octants[idx].getPbound(global3D.nodeface[inode][2])){
3626  findNodeNeighbours(idx, inode, neigh, isghost);
3627  sizeneigh = neigh.size();
3628  for(i=0; i<sizeneigh; i++){
3629  if (!isghost[i]){
3630  {
3631  if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) > (targetmarker + 1)){
3632  octants[idx].setMarker(octants[neigh[i]].getLevel()+octants[neigh[i]].getMarker()-octants[idx].getLevel()-1);
3633  octants[idx].info[15] = true;
3634  newmodified.push_back(idx);
3635  Bdone = true;
3636  }
3637  else if((octants[neigh[i]].getLevel() + octants[neigh[i]].getMarker()) < (targetmarker - 1)){
3638  octants[neigh[i]].setMarker(targetmarker-octants[neigh[i]].getLevel()-1);
3639  octants[neigh[i]].info[15] = true;
3640  newmodified.push_back(neigh[i]);
3641  Bdone = true;
3642  }
3643  };
3644  }
3645  }
3646  }
3647  }
3648  }
3649 
3650  }
3651  }
3652  preBalance21(newmodified);
3653  u32vector().swap(modified);
3654  swap(modified,newmodified);
3655  modsize = modified.size();
3656  u32vector().swap(newmodified);
3657  }// end while
3658  obegin = oend = octants.end();
3659  ibegin = iend = modified.end();
3660  }
3661  return Bdone;
3662  // Pay attention : info[15] may be true after local balance for some octants
3663  };
3664 
3665  // =================================================================================== //
3666 
3667  void findEdgeNeighbours(uint32_t idx, // Finds neighbours of idx-th octant through iedge in vector octants.
3668  uint8_t iedge, // Returns a vector (empty if iedge is a bound edge) with the index of neighbours
3669  u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
3670  vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
3671 
3672  uint64_t Morton, Mortontry;
3673  uint32_t noctants = getNumOctants();
3674  uint32_t idxtry;
3675  Class_Octant<3>* oct = &octants[idx];
3676  uint32_t size = oct->getSize();
3677  uint8_t iface1, iface2;
3678  int32_t Dx, Dy, Dz;
3679  int32_t Dxstar,Dystar,Dzstar;
3680 
3681  //Alternative to switch case
3682  // int8_t Cx[12] = {-1,1,0,0,-1,1,-1,1,-1,1,0,0};
3683  // int8_t Cy[12] = {0,0,-1,1,-1,-1,1,1,0,0,-1,1};
3684  // int8_t Cz[12] = {-1,-1,-1,-1,0,0,0,0,1,1,1,1};
3685  int8_t cx = global3D.edgecoeffs[iedge][0];
3686  int8_t cy = global3D.edgecoeffs[iedge][1];
3687  int8_t cz = global3D.edgecoeffs[iedge][2];
3688 
3689  isghost.clear();
3690  neighbours.clear();
3691 
3692  // Default if iedge is nface<iedge<0
3693  if (iedge < 0 || iedge > global3D.nfaces*2){
3694  return;
3695  }
3696 
3697  // Check if octants edge is a process boundary
3698  iface1 = global3D.edgeface[iedge][0];
3699  iface2 = global3D.edgeface[iedge][1];
3700 
3701  // Check if octants edge is a boundary
3702  if (oct->info[iface1] == false && oct->info[iface2] == false){
3703 
3704  //Build Morton number of virtual neigh of same size
3705  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
3706  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
3707 
3708  //SEARCH IN GHOSTS
3709 
3710  if (ghosts.size()>0){
3711  // Search in ghosts
3712  uint32_t idxghost = uint32_t(size_ghosts/2);
3713  Class_Octant<3>* octghost = &ghosts[idxghost];
3714 
3715  // Search morton in octants
3716  // If a even face morton is lower than morton of oct, if odd higher
3717  // ---> can i search only before or after idx in octants
3718  int32_t jump = int32_t(idxghost/2);
3719  idxtry = uint32_t(((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
3720  if (idxtry > ghosts.size()-1) idxtry = ghosts.size()-1;
3721  while(abs(jump) > 0){
3722  Mortontry = ghosts[idxtry].computeMorton();
3723  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
3724  idxtry += jump;
3725  if (idxtry > ghosts.size()-1){
3726  if (jump > 0){
3727  idxtry = ghosts.size() - 1;
3728  jump = 0;
3729  }
3730  else if (jump < 0){
3731  idxtry = 0;
3732  jump = 0;
3733  }
3734  }
3735  }
3736  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
3737  //Found neighbour of same size
3738  isghost.push_back(true);
3739  neighbours.push_back(idxtry);
3740  return;
3741  }
3742  else{
3743  // Step until the mortontry lower than morton (one idx of distance)
3744  {
3745  while(ghosts[idxtry].computeMorton() < Morton){
3746  idxtry++;
3747  if(idxtry > ghosts.size()-1){
3748  idxtry = ghosts.size()-1;
3749  break;
3750  }
3751  }
3752  while(ghosts[idxtry].computeMorton() > Morton){
3753  idxtry--;
3754  if(idxtry > ghosts.size()-1){
3755  idxtry = 0;
3756  break;
3757  }
3758  }
3759  }
3760  if(idxtry < size_ghosts){
3761  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
3762  //Found neighbour of same size
3763  isghost.push_back(true);
3764  neighbours.push_back(idxtry);
3765  return;
3766  }
3767  // Compute Last discendent of virtual octant of same size
3768  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
3769  uint64_t Mortonlast = last_desc.computeMorton();
3770  Mortontry = ghosts[idxtry].computeMorton();
3771  while(Mortontry < Mortonlast && idxtry < ghosts.size()){
3772  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(ghosts[idxtry].x));
3773  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(ghosts[idxtry].y));
3774  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(ghosts[idxtry].z));
3775  Dxstar = int32_t((cx-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cx+1)/2)*size;
3776  Dystar = int32_t((cy-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cy+1)/2)*size;
3777  Dzstar = int32_t((cz-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cz+1)/2)*size;
3778 
3779  uint32_t x0 = oct->x;
3780  uint32_t x1 = x0 + size;
3781  uint32_t y0 = oct->y;
3782  uint32_t y1 = y0 + size;
3783  uint32_t z0 = oct->z;
3784  uint32_t z1 = z0 + size;
3785  uint32_t x0try = ghosts[idxtry].x;
3786  uint32_t x1try = x0try + ghosts[idxtry].getSize();
3787  uint32_t y0try = ghosts[idxtry].y;
3788  uint32_t y1try = y0try + ghosts[idxtry].getSize();
3789  uint32_t z0try = ghosts[idxtry].z;
3790  uint32_t z1try = z0try + ghosts[idxtry].getSize();
3791  uint8_t level = oct->level;
3792  uint8_t leveltry = ghosts[idxtry].getLevel();
3793 
3794  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
3795  if (leveltry > level){
3796  if((abs(cx)*abs(cz)*((y0try>=y0)*(y0try<y1))) + (abs(cy)*abs(cz)*((x0try>=x0)*(x0try<x1))) + (abs(cx)*abs(cy)*((z0try>=z0)*(z0try<z1)))){
3797  neighbours.push_back(idxtry);
3798  isghost.push_back(true);
3799  }
3800  }
3801  if (leveltry < level){
3802  if((abs(cx)*abs(cz)*((y0>=y0try)*(y0<y1try))) + (abs(cy)*abs(cz)*((x0>=x0try)*(x0<x1try))) + (abs(cx)*abs(cy)*((z0>=z0try)*(z0<z1try)))){
3803  neighbours.push_back(idxtry);
3804  isghost.push_back(true);
3805  }
3806  }
3807  }
3808 
3809  idxtry++;
3810  if(idxtry>size_ghosts-1){
3811  break;
3812  }
3813  Mortontry = ghosts[idxtry].computeMorton();
3814  }
3815  }
3816  }
3817  }
3818 
3819  // Search in octants
3820 
3821  //Build Morton number of virtual neigh of same size
3822  // Search morton in octants
3823  // If a even face morton is lower than morton of oct, if odd higher
3824  // ---> can i search only before or after idx in octants
3825 // int32_t jump = int32_t(noctants/2+1);
3826 // idxtry = uint32_t(((oct->computeMorton()<Morton)-(oct->computeMorton()>Morton))*jump);
3827  int32_t jump = (oct->computeMorton() > Morton) ? int32_t(idx/2+1) : int32_t((noctants -idx)/2+1);
3828  idxtry = uint32_t(idx +((oct->computeMorton()<Morton)-(oct->computeMorton()>Morton))*jump);
3829  if (idxtry > noctants-1)
3830  idxtry = noctants-1;
3831  while(abs(jump) > 0){
3832  Mortontry = octants[idxtry].computeMorton();
3833  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
3834  idxtry += jump;
3835  if (idxtry > octants.size()-1){
3836  if (jump > 0){
3837  idxtry = octants.size() - 1;
3838  jump = 0;
3839  }
3840  else if (jump < 0){
3841  idxtry = 0;
3842  jump = 0;
3843  }
3844  }
3845  }
3846  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
3847  //Found neighbour of same size
3848  isghost.push_back(false);
3849  neighbours.push_back(idxtry);
3850  return;
3851  }
3852  else{
3853  // Step until the mortontry lower than morton (one idx of distance)
3854  {
3855  while(octants[idxtry].computeMorton() < Morton){
3856  idxtry++;
3857  if(idxtry > noctants-1){
3858  idxtry = noctants-1;
3859  break;
3860  }
3861  }
3862  while(octants[idxtry].computeMorton() > Morton){
3863  idxtry--;
3864  if(idxtry > noctants-1){
3865  idxtry = 0;
3866  break;
3867  }
3868  }
3869  }
3870  if (idxtry < noctants){
3871  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
3872  //Found neighbour of same size
3873  isghost.push_back(false);
3874  neighbours.push_back(idxtry);
3875  return;
3876  }
3877  // Compute Last discendent of virtual octant of same size
3878  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
3879  uint64_t Mortonlast = last_desc.computeMorton();
3880  Mortontry = octants[idxtry].computeMorton();
3881  while(Mortontry < Mortonlast && idxtry <= noctants-1){
3882  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
3883  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
3884  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
3885  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
3886  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
3887  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
3888 
3889  uint32_t x0 = oct->x;
3890  uint32_t x1 = x0 + size;
3891  uint32_t y0 = oct->y;
3892  uint32_t y1 = y0 + size;
3893  uint32_t z0 = oct->z;
3894  uint32_t z1 = z0 + size;
3895  uint32_t x0try = octants[idxtry].x;
3896  uint32_t x1try = x0try + octants[idxtry].getSize();
3897  uint32_t y0try = octants[idxtry].y;
3898  uint32_t y1try = y0try + octants[idxtry].getSize();
3899  uint32_t z0try = octants[idxtry].z;
3900  uint32_t z1try = z0try + octants[idxtry].getSize();
3901  uint8_t level = oct->level;
3902  uint8_t leveltry = octants[idxtry].getLevel();
3903 
3904  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
3905  if (leveltry > level){
3906  if((abs(cx)*abs(cz)*((y0try>=y0)*(y0try<y1))) + (abs(cy)*abs(cz)*((x0try>=x0)*(x0try<x1))) + (abs(cx)*abs(cy)*((z0try>=z0)*(z0try<z1)))){
3907  neighbours.push_back(idxtry);
3908  isghost.push_back(false);
3909  }
3910  }
3911  if (leveltry < level){
3912  if((abs(cx)*abs(cz)*((y0>=y0try)*(y0<y1try))) + (abs(cy)*abs(cz)*((x0>=x0try)*(x0<x1try))) + (abs(cx)*abs(cy)*((z0>=z0try)*(z0<z1try)))){
3913  neighbours.push_back(idxtry);
3914  isghost.push_back(false);
3915  }
3916  }
3917  }
3918 
3919  idxtry++;
3920  if(idxtry>noctants-1){
3921  break;
3922  }
3923  Mortontry = octants[idxtry].computeMorton();
3924  }
3925  }
3926  }
3927  return;
3928  }
3929  else{
3930  // Boundary Face
3931  return;
3932  }
3933  };
3934 
3935  // =================================================================================== //
3936 
3937  void findEdgeNeighbours(Class_Octant<3>* oct, // Finds neighbours of idx-th octant through iedge in vector octants.
3938  uint8_t iedge, // Returns a vector (empty if iedge is a bound edge) with the index of neighbours
3939  u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
3940  vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
3941 
3942  uint64_t Morton, Mortontry;
3943  uint32_t noctants = getNumOctants();
3944  uint32_t idxtry;
3945  uint32_t size = oct->getSize();
3946  uint8_t iface1, iface2;
3947  int32_t Dx, Dy, Dz;
3948  int32_t Dxstar,Dystar,Dzstar;
3949 
3950  //Alternative to switch case
3951  // int8_t Cx[12] = {-1,1,0,0,-1,1,-1,1,-1,1,0,0};
3952  // int8_t Cy[12] = {0,0,-1,1,-1,-1,1,1,0,0,-1,1};
3953  // int8_t Cz[12] = {-1,-1,-1,-1,0,0,0,0,1,1,1,1};
3954  int8_t cx = global3D.edgecoeffs[iedge][0];
3955  int8_t cy = global3D.edgecoeffs[iedge][1];
3956  int8_t cz = global3D.edgecoeffs[iedge][2];
3957 
3958  isghost.clear();
3959  neighbours.clear();
3960 
3961  // Default if iedge is nface<iedge<0
3962  if (iedge < 0 || iedge > global3D.nfaces*2){
3963  return;
3964  }
3965 
3966  // Check if octants edge is a process boundary
3967  iface1 = global3D.edgeface[iedge][0];
3968  iface2 = global3D.edgeface[iedge][1];
3969 
3970  // Check if octants edge is a boundary
3971  if (oct->info[iface1] == false && oct->info[iface2] == false){
3972 
3973  //Build Morton number of virtual neigh of same size
3974  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
3975  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
3976 
3977  //SEARCH IN GHOSTS
3978 
3979  if (ghosts.size()>0){
3980  // Search in ghosts
3981  uint32_t idxghost = uint32_t(size_ghosts/2);
3982  Class_Octant<3>* octghost = &ghosts[idxghost];
3983 
3984  // Search morton in octants
3985  // If a even face morton is lower than morton of oct, if odd higher
3986  // ---> can i search only before or after idx in octants
3987  int32_t jump = int32_t(idxghost/2+1);
3988  idxtry = uint32_t(idxghost +((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
3989  while(abs(jump) > 0){
3990  Mortontry = ghosts[idxtry].computeMorton();
3991  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
3992  idxtry += jump;
3993  if (idxtry > ghosts.size()-1){
3994  if (jump > 0){
3995  idxtry = ghosts.size() - 1;
3996  jump = 0;
3997  }
3998  else if (jump < 0){
3999  idxtry = 0;
4000  jump = 0;
4001  }
4002  }
4003  }
4004  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4005  //Found neighbour of same size
4006  isghost.push_back(true);
4007  neighbours.push_back(idxtry);
4008  return;
4009  }
4010  else{
4011  // Step until the mortontry lower than morton (one idx of distance)
4012  {
4013  while(ghosts[idxtry].computeMorton() < Morton){
4014  idxtry++;
4015  if(idxtry > ghosts.size()-1){
4016  idxtry = ghosts.size()-1;
4017  break;
4018  }
4019  }
4020  while(ghosts[idxtry].computeMorton() > Morton){
4021  idxtry--;
4022  if(idxtry > ghosts.size()-1){
4023  idxtry = 0;
4024  break;
4025  }
4026  }
4027  }
4028  if(idxtry < size_ghosts){
4029  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4030  //Found neighbour of same size
4031  isghost.push_back(true);
4032  neighbours.push_back(idxtry);
4033  return;
4034  }
4035  // Compute Last discendent of virtual octant of same size
4036  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4037  uint64_t Mortonlast = last_desc.computeMorton();
4038  Mortontry = ghosts[idxtry].computeMorton();
4039  while(Mortontry < Mortonlast && idxtry < ghosts.size()){
4040  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(ghosts[idxtry].x));
4041  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(ghosts[idxtry].y));
4042  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(ghosts[idxtry].z));
4043  Dxstar = int32_t((cx-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cx+1)/2)*size;
4044  Dystar = int32_t((cy-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cy+1)/2)*size;
4045  Dzstar = int32_t((cz-1)/2)*(ghosts[idxtry].getSize()) + int32_t((cz+1)/2)*size;
4046 
4047  uint32_t x0 = oct->x;
4048  uint32_t x1 = x0 + size;
4049  uint32_t y0 = oct->y;
4050  uint32_t y1 = y0 + size;
4051  uint32_t z0 = oct->z;
4052  uint32_t z1 = z0 + size;
4053  uint32_t x0try = ghosts[idxtry].x;
4054  uint32_t x1try = x0try + ghosts[idxtry].getSize();
4055  uint32_t y0try = ghosts[idxtry].y;
4056  uint32_t y1try = y0try + ghosts[idxtry].getSize();
4057  uint32_t z0try = ghosts[idxtry].z;
4058  uint32_t z1try = z0try + ghosts[idxtry].getSize();
4059  uint8_t level = oct->level;
4060  uint8_t leveltry = ghosts[idxtry].getLevel();
4061 
4062  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
4063  if (leveltry > level){
4064  if((abs(cx)*abs(cz)*((y0try>=y0)*(y0try<y1))) + (abs(cy)*abs(cz)*((x0try>=x0)*(x0try<x1))) + (abs(cx)*abs(cy)*((z0try>=z0)*(z0try<z1)))){
4065  neighbours.push_back(idxtry);
4066  isghost.push_back(true);
4067  }
4068  }
4069  if (leveltry < level){
4070  if((abs(cx)*abs(cz)*((y0>=y0try)*(y0<y1try))) + (abs(cy)*abs(cz)*((x0>=x0try)*(x0<x1try))) + (abs(cx)*abs(cy)*((z0>=z0try)*(z0<z1try)))){
4071  neighbours.push_back(idxtry);
4072  isghost.push_back(true);
4073  }
4074  }
4075  }
4076 
4077  idxtry++;
4078  if(idxtry>size_ghosts-1){
4079  break;
4080  }
4081  Mortontry = ghosts[idxtry].computeMorton();
4082  }
4083  }
4084  }
4085  }
4086 
4087  // Search in octants
4088 
4089  //Build Morton number of virtual neigh of same size
4090  // Search morton in octants
4091  // If a even face morton is lower than morton of oct, if odd higher
4092  // ---> can i search only before or after idx in octants
4093  int32_t jump = int32_t((noctants)/2+1);
4094  idxtry = uint32_t(jump);
4095  if (idxtry > noctants-1)
4096  idxtry = noctants-1;
4097  while(abs(jump) > 0){
4098  Mortontry = octants[idxtry].computeMorton();
4099  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4100  idxtry += jump;
4101  if (idxtry > octants.size()-1){
4102  if (jump > 0){
4103  idxtry = octants.size() - 1;
4104  jump = 0;
4105  }
4106  else if (jump < 0){
4107  idxtry = 0;
4108  jump = 0;
4109  }
4110  }
4111  }
4112  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4113  //Found neighbour of same size
4114  isghost.push_back(false);
4115  neighbours.push_back(idxtry);
4116  return;
4117  }
4118  else{
4119  // Step until the mortontry lower than morton (one idx of distance)
4120  {
4121  while(octants[idxtry].computeMorton() < Morton){
4122  idxtry++;
4123  if(idxtry > noctants-1){
4124  idxtry = noctants-1;
4125  break;
4126  }
4127  }
4128  while(octants[idxtry].computeMorton() > Morton){
4129  idxtry--;
4130  if(idxtry > noctants-1){
4131  idxtry = 0;
4132  break;
4133  }
4134  }
4135  }
4136  if (idxtry < noctants){
4137  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4138  //Found neighbour of same size
4139  isghost.push_back(false);
4140  neighbours.push_back(idxtry);
4141  return;
4142  }
4143  // Compute Last discendent of virtual octant of same size
4144  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4145  uint64_t Mortonlast = last_desc.computeMorton();
4146  Mortontry = octants[idxtry].computeMorton();
4147  while(Mortontry < Mortonlast && idxtry <= noctants-1){
4148  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
4149  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
4150  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
4151  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
4152  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
4153  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
4154 
4155  uint32_t x0 = oct->x;
4156  uint32_t x1 = x0 + size;
4157  uint32_t y0 = oct->y;
4158  uint32_t y1 = y0 + size;
4159  uint32_t z0 = oct->z;
4160  uint32_t z1 = z0 + size;
4161  uint32_t x0try = octants[idxtry].x;
4162  uint32_t x1try = x0try + octants[idxtry].getSize();
4163  uint32_t y0try = octants[idxtry].y;
4164  uint32_t y1try = y0try + octants[idxtry].getSize();
4165  uint32_t z0try = octants[idxtry].z;
4166  uint32_t z1try = z0try + octants[idxtry].getSize();
4167  uint8_t level = oct->level;
4168  uint8_t leveltry = octants[idxtry].getLevel();
4169 
4170  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
4171  if (leveltry > level){
4172  if((abs(cx)*abs(cz)*((y0try>=y0)*(y0try<y1))) + (abs(cy)*abs(cz)*((x0try>=x0)*(x0try<x1))) + (abs(cx)*abs(cy)*((z0try>=z0)*(z0try<z1)))){
4173  neighbours.push_back(idxtry);
4174  isghost.push_back(false);
4175  }
4176  }
4177  if (leveltry < level){
4178  if((abs(cx)*abs(cz)*((y0>=y0try)*(y0<y1try))) + (abs(cy)*abs(cz)*((x0>=x0try)*(x0<x1try))) + (abs(cx)*abs(cy)*((z0>=z0try)*(z0<z1try)))){
4179  neighbours.push_back(idxtry);
4180  isghost.push_back(false);
4181  }
4182  }
4183  }
4184 
4185  idxtry++;
4186  if(idxtry>noctants-1){
4187  break;
4188  }
4189  Mortontry = octants[idxtry].computeMorton();
4190  }
4191  }
4192  }
4193  return;
4194  }
4195  else{
4196  // Boundary Face
4197  return;
4198  }
4199 
4200  };
4201 
4202  // =================================================================================== //
4203 
4204  void findGhostEdgeNeighbours(uint32_t idx, // Finds neighbours of idx-th ghost through iedge in vector octants.
4205  uint8_t iedge, // Returns a vector (empty if iedge is not a pbound edge) with the index of neighbours
4206  u32vector & neighbours){ // in the structure octants.
4207 
4208  uint64_t Morton, Mortontry;
4209  uint32_t noctants = getNumOctants();
4210  uint32_t idxtry;
4211  Class_Octant<3>* oct = &ghosts[idx];
4212  uint32_t size = oct->getSize();
4213  uint8_t iface1, iface2;
4214  int32_t Dx, Dy, Dz;
4215  int32_t Dxstar,Dystar,Dzstar;
4216 
4217  //Alternative to switch case
4218  int8_t cx = global3D.edgecoeffs[iedge][0];
4219  int8_t cy = global3D.edgecoeffs[iedge][1];
4220  int8_t cz = global3D.edgecoeffs[iedge][2];
4221 
4222  neighbours.clear();
4223 
4224  // Default if iedge is nface<iedge<0
4225  if (iedge < 0 || iedge > global3D.nfaces*2){
4226  return;
4227  }
4228 
4229  // Check if octants edge is a process boundary
4230  iface1 = global3D.edgeface[iedge][0];
4231  iface2 = global3D.edgeface[iedge][1];
4232 
4233  // Check if octants edge is a pboundary edge
4234  if (oct->info[iface1+global3D.nfaces] == true || oct->info[iface2+global3D.nfaces] == true){
4235 
4236  //Build Morton number of virtual neigh of same size
4237  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
4238  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
4239 
4240  //Build Morton number of virtual neigh of same size
4241  // Search morton in octants
4242  // If a even face morton is lower than morton of oct, if odd higher
4243  // ---> can i search only before or after idx in octants
4244  int32_t jump = getNumOctants()/2;
4245  idxtry = uint32_t(getNumOctants()/2);
4246  while(abs(jump) > 0){
4247  Mortontry = octants[idxtry].computeMorton();
4248  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4249  idxtry += jump;
4250  if (idxtry > octants.size()-1){
4251  if (jump > 0){
4252  idxtry = octants.size() - 1;
4253  jump = 0;
4254  }
4255  else if (jump < 0){
4256  idxtry = 0;
4257  jump = 0;
4258  }
4259  }
4260  }
4261  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4262  //Found neighbour of same size
4263  neighbours.push_back(idxtry);
4264  return;
4265  }
4266  else{
4267  // Step until the mortontry lower than morton (one idx of distance)
4268  {
4269  while(octants[idxtry].computeMorton() < Morton){
4270  idxtry++;
4271  if(idxtry > noctants-1){
4272  idxtry = noctants-1;
4273  break;
4274  }
4275  }
4276  while(octants[idxtry].computeMorton() > Morton){
4277  idxtry--;
4278  if(idxtry > noctants-1){
4279  idxtry = 0;
4280  break;
4281  }
4282  }
4283  }
4284  if (idxtry < noctants){
4285  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4286  //Found neighbour of same size
4287  neighbours.push_back(idxtry);
4288  return;
4289  }
4290  // Compute Last discendent of virtual octant of same size
4291  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4292  uint64_t Mortonlast = last_desc.computeMorton();
4293  Mortontry = octants[idxtry].computeMorton();
4294  while(Mortontry < Mortonlast && idxtry <= noctants-1){
4295  Dx = int32_t(abs(cx))*(-int32_t(oct->x) + int32_t(octants[idxtry].x));
4296  Dy = int32_t(abs(cy))*(-int32_t(oct->y) + int32_t(octants[idxtry].y));
4297  Dz = int32_t(abs(cz))*(-int32_t(oct->z) + int32_t(octants[idxtry].z));
4298  Dxstar = int32_t((cx-1)/2)*(octants[idxtry].getSize()) + int32_t((cx+1)/2)*size;
4299  Dystar = int32_t((cy-1)/2)*(octants[idxtry].getSize()) + int32_t((cy+1)/2)*size;
4300  Dzstar = int32_t((cz-1)/2)*(octants[idxtry].getSize()) + int32_t((cz+1)/2)*size;
4301 
4302  uint32_t x0 = oct->x;
4303  uint32_t x1 = x0 + size;
4304  uint32_t y0 = oct->y;
4305  uint32_t y1 = y0 + size;
4306  uint32_t z0 = oct->z;
4307  uint32_t z1 = z0 + size;
4308  uint32_t x0try = octants[idxtry].x;
4309  uint32_t x1try = x0try + octants[idxtry].getSize();
4310  uint32_t y0try = octants[idxtry].y;
4311  uint32_t y1try = y0try + octants[idxtry].getSize();
4312  uint32_t z0try = octants[idxtry].z;
4313  uint32_t z1try = z0try + octants[idxtry].getSize();
4314  uint8_t level = oct->level;
4315  uint8_t leveltry = octants[idxtry].getLevel();
4316 
4317  if (Dx == Dxstar && Dy == Dystar && Dz == Dzstar){
4318  if (leveltry > level){
4319  if((abs(cx)*abs(cz)*((y0try>=y0)*(y0try<y1))) + (abs(cy)*abs(cz)*((x0try>=x0)*(x0try<x1))) + (abs(cx)*abs(cy)*((z0try>=z0)*(z0try<z1)))){
4320  neighbours.push_back(idxtry);
4321  }
4322  }
4323  if (leveltry < level){
4324  if((abs(cx)*abs(cz)*((y0>=y0try)*(y0<y1try))) + (abs(cy)*abs(cz)*((x0>=x0try)*(x0<x1try))) + (abs(cx)*abs(cy)*((z0>=z0try)*(z0<z1try)))){
4325  neighbours.push_back(idxtry);
4326  }
4327  }
4328  }
4329  idxtry++;
4330  if(idxtry>noctants-1){
4331  break;
4332  }
4333  Mortontry = octants[idxtry].computeMorton();
4334  }
4335  }
4336  }
4337  return;
4338  }
4339  else{
4340  // Boundary Face
4341  return;
4342  }
4343  };
4344 
4345  // =================================================================================== //
4346 
4347 // void findNodeNeighbours(uint32_t idx, // Finds neighbours of idx-th octant through inode in vector octants.
4348 // uint8_t inode, // Returns a vector (empty if inode is a bound node) with the index of neighbours
4349 // u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
4350 // vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
4351 //
4352 // uint64_t Morton, Mortontry;
4353 // uint32_t noctants = getNumOctants();
4354 // uint32_t idxtry;
4355 // Class_Octant<3>* oct = &octants[idx];
4356 // uint32_t size = oct->getSize();
4357 // uint8_t iface1, iface2, iface3;
4358 // int32_t Dhx, Dhy, Dhz;
4359 // int32_t Dhxref, Dhyref, Dhzref;
4360 //
4361 // //Alternative to switch case
4362 // int8_t Cx[8] = {-1,1,-1,1,-1,1,-1,1};
4363 // int8_t Cy[8] = {-1,-1,1,1,-1,-1,1,1};
4364 // int8_t Cz[8] = {-1,-1,-1,-1,1,1,1,1};
4365 // int8_t cx = Cx[inode];
4366 // int8_t cy = Cy[inode];
4367 // int8_t cz = Cz[inode];
4368 //
4369 // isghost.clear();
4370 // neighbours.clear();
4371 //
4372 // // Default if inode is nnodes<inode<0
4373 // if (inode < 0 || inode > global3D.nnodes){
4374 // return;
4375 // }
4376 //
4377 // // Check if octants node is a boundary
4378 // iface1 = global3D.nodeface[inode][0];
4379 // iface2 = global3D.nodeface[inode][1];
4380 // iface3 = global3D.nodeface[inode][2];
4381 //
4382 // // Check if octants node is a boundary
4383 // if (oct->info[iface1] == false && oct->info[iface2] == false && oct->info[iface3] == false){
4384 //
4385 // //Build Morton number of virtual neigh of same size
4386 // Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
4387 // Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
4388 //
4389 // //SEARCH IN GHOSTS
4390 //
4391 // if (ghosts.size()>0){
4392 //
4393 // // Search in ghosts
4394 // uint32_t idxghost = uint32_t(size_ghosts/2);
4395 // Class_Octant<3>* octghost = &ghosts[idxghost];
4396 //
4397 // // Search morton in octants
4398 // // If a even face morton is lower than morton of oct, if odd higher
4399 // // ---> can i search only before or after idx in octants
4400 // int32_t jump;
4401 // if (inode==7 || inode ==0){
4402 // jump = (octghost->computeMorton() > Morton) ? int32_t(idxghost/2+1) : int32_t((size_ghosts -idxghost)/2+1);
4403 // idxtry = uint32_t(idxghost +((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
4404 // if (idxtry > size_ghosts-1)
4405 // idxtry = size_ghosts-1;
4406 // }
4407 // else{
4408 // jump = idxghost;
4409 // idxtry = uint32_t(jump);
4410 // }
4411 // while(abs(jump) > 0){
4412 // Mortontry = ghosts[idxtry].computeMorton();
4413 // jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4414 // idxtry += jump;
4415 // if (idxtry > ghosts.size()-1){
4416 // if (jump > 0){
4417 // idxtry = ghosts.size() - 1;
4418 // jump = 0;
4419 // }
4420 // else if (jump < 0){
4421 // idxtry = 0;
4422 // jump = 0;
4423 // }
4424 // }
4425 // }
4426 // if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4427 // //Found neighbour of same size
4428 // isghost.push_back(true);
4429 // neighbours.push_back(idxtry);
4430 // return;
4431 // }
4432 // else{
4433 // // Step until the mortontry lower than morton (one idx of distance)
4434 // {
4435 // while(ghosts[idxtry].computeMorton() < Morton){
4436 // idxtry++;
4437 // if(idxtry > ghosts.size()-1){
4438 // idxtry = ghosts.size()-1;
4439 // break;
4440 // }
4441 // }
4442 // while(ghosts[idxtry].computeMorton() > Morton){
4443 // idxtry--;
4444 // if(idxtry > ghosts.size()-1){
4445 // idxtry = 0;
4446 // break;
4447 // }
4448 // }
4449 // }
4450 // if(idxtry < size_ghosts){
4451 // if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4452 // //Found neighbour of same size
4453 // isghost.push_back(true);
4454 // neighbours.push_back(idxtry);
4455 // return;
4456 // }
4457 // // Compute Last discendent of virtual octant of same size
4458 // Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4459 // uint64_t Mortonlast = last_desc.computeMorton();
4460 // Mortontry = ghosts[idxtry].computeMorton();
4461 // while(Mortontry < Mortonlast && idxtry < size_ghosts){
4462 // // Dhx = (-int32_t(oct->x) + int32_t(ghosts[idxtry].x));
4463 // // Dhy = (-int32_t(oct->y) + int32_t(ghosts[idxtry].y));
4464 // // Dhz = (-int32_t(oct->z) + int32_t(ghosts[idxtry].z));
4465 // // Dhxref = int32_t(cx<0)*(-int32_t(ghosts[idxtry].getSize())) + int32_t(cx>0)*size;
4466 // // Dhyref = int32_t(cy<0)*(-int32_t(ghosts[idxtry].getSize())) + int32_t(cy>0)*size;
4467 // // Dhzref = int32_t(cz<0)*(-int32_t(ghosts[idxtry].getSize())) + int32_t(cz>0)*size;
4468 // Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(ghosts[idxtry].x));
4469 // Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(ghosts[idxtry].y));
4470 // Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(ghosts[idxtry].z));
4471 // Dhxref = int32_t(cx<0)*ghosts[idxtry].getSize() + int32_t(cx>0)*size;
4472 // Dhyref = int32_t(cy<0)*ghosts[idxtry].getSize() + int32_t(cy>0)*size;
4473 // Dhzref = int32_t(cz<0)*ghosts[idxtry].getSize() + int32_t(cz>0)*size;
4474 // if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
4475 // neighbours.push_back(idxtry);
4476 // isghost.push_back(true);
4477 // return;
4478 // }
4479 // idxtry++;
4480 // if(idxtry>size_ghosts-1){
4481 // break;
4482 // }
4483 // Mortontry = ghosts[idxtry].computeMorton();
4484 // }
4485 // }
4486 // }
4487 // }
4488 //
4489 // // Search in octants
4490 //
4491 // //Build Morton number of virtual neigh of same size
4492 // // Search morton in octants
4493 // // If a even face morton is lower than morton of oct, if odd higher
4494 // // ---> can i search only before or after idx in octants
4495 // int32_t jump;
4496 // if (inode==0 || inode==7 ){
4497 // jump = (oct->computeMorton() > Morton) ? int32_t(idx/2+1) : int32_t((noctants -idx)/2+1);
4498 // idxtry = uint32_t(idx +((oct->computeMorton()<Morton)-(oct->computeMorton()>Morton))*jump);
4499 // if (idxtry > noctants-1)
4500 // idxtry = noctants-1;
4501 // }
4502 // else{
4503 // jump = noctants/2;
4504 // idxtry = jump;
4505 // }
4506 // while(abs(jump) > 0){
4507 // Mortontry = octants[idxtry].computeMorton();
4508 // jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4509 // idxtry += jump;
4510 // if (idxtry > octants.size()-1){
4511 // if (jump > 0){
4512 // idxtry = octants.size() - 1;
4513 // jump = 0;
4514 // }
4515 // else if (jump < 0){
4516 // idxtry = 0;
4517 // jump = 0;
4518 // }
4519 // }
4520 // }
4521 // if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4522 // //Found neighbour of same size
4523 // isghost.push_back(false);
4524 // neighbours.push_back(idxtry);
4525 // return;
4526 // }
4527 // else{
4528 // // Step until the mortontry lower than morton (one idx of distance)
4529 // {
4530 // while(octants[idxtry].computeMorton() < Morton){
4531 // idxtry++;
4532 // if(idxtry > noctants-1){
4533 // idxtry = noctants-1;
4534 // break;
4535 // }
4536 // }
4537 // while(octants[idxtry].computeMorton() > Morton){
4538 // idxtry--;
4539 // if(idxtry > noctants-1){
4540 // idxtry = 0;
4541 // break;
4542 // }
4543 // }
4544 // }
4545 // if (idxtry < noctants){
4546 // if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4547 // //Found neighbour of same size
4548 // isghost.push_back(false);
4549 // neighbours.push_back(idxtry);
4550 // return;
4551 // }
4552 // // Compute Last discendent of virtual octant of same size
4553 // Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4554 // uint64_t Mortonlast = last_desc.computeMorton();
4555 // Mortontry = octants[idxtry].computeMorton();
4556 // while(Mortontry < Mortonlast && idxtry <= noctants-1){
4557 // // Dhx = (-int32_t(oct->x) + int32_t(octants[idxtry].x));
4558 // // Dhy = (-int32_t(oct->y) + int32_t(octants[idxtry].y));
4559 // // Dhz = (-int32_t(oct->z) + int32_t(octants[idxtry].z));
4560 // // Dhxref = int32_t(cx<0)*(-int32_t(octants[idxtry].getSize())) + int32_t(cx>0)*size;
4561 // // Dhyref = int32_t(cy<0)*(-int32_t(octants[idxtry].getSize())) + int32_t(cy>0)*size;
4562 // // Dhzref = int32_t(cz<0)*(-int32_t(octants[idxtry].getSize())) + int32_t(cz>0)*size;
4563 // Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(octants[idxtry].x));
4564 // Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(octants[idxtry].y));
4565 // Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(octants[idxtry].z));
4566 // Dhxref = int32_t(cx<0)*octants[idxtry].getSize() + int32_t(cx>0)*size;
4567 // Dhyref = int32_t(cy<0)*octants[idxtry].getSize() + int32_t(cy>0)*size;
4568 // Dhzref = int32_t(cz<0)*octants[idxtry].getSize() + int32_t(cz>0)*size;
4569 // if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
4570 // neighbours.push_back(idxtry);
4571 // isghost.push_back(false);
4572 // return;
4573 // }
4574 // idxtry++;
4575 // if(idxtry>noctants-1){
4576 // break;
4577 // }
4578 // Mortontry = octants[idxtry].computeMorton();
4579 // }
4580 // }
4581 // }
4582 // return;
4583 // }
4584 // else{
4585 // // Boundary Node
4586 // return;
4587 // }
4588 //
4589 // };
4590 
4591  // =================================================================================== //
4592 
4593  void findNodeNeighbours(Class_Octant<3>* oct, // Finds neighbours of idx-th octant through inode in vector octants.
4594  uint8_t inode, // Returns a vector (empty if inode is a bound node) with the index of neighbours
4595  u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
4596  vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
4597 
4598  uint64_t Morton, Mortontry;
4599  uint32_t noctants = getNumOctants();
4600  uint32_t idxtry;
4601  uint32_t size = oct->getSize();
4602  uint8_t iface1, iface2, iface3;
4603  int32_t Dhx, Dhy, Dhz;
4604  int32_t Dhxref, Dhyref, Dhzref;
4605 
4606  //Alternative to switch case
4607  int8_t Cx[8] = {-1,1,-1,1,-1,1,-1,1};
4608  int8_t Cy[8] = {-1,-1,1,1,-1,-1,1,1};
4609  int8_t Cz[8] = {-1,-1,-1,-1,1,1,1,1};
4610  int8_t cx = Cx[inode];
4611  int8_t cy = Cy[inode];
4612  int8_t cz = Cz[inode];
4613 
4614  isghost.clear();
4615  neighbours.clear();
4616 
4617  // Default if inode is nnodes<inode<0
4618  if (inode < 0 || inode > global3D.nnodes){
4619  return;
4620  }
4621 
4622  // Check if octants node is a boundary
4623  iface1 = global3D.nodeface[inode][0];
4624  iface2 = global3D.nodeface[inode][1];
4625  iface3 = global3D.nodeface[inode][2];
4626 
4627  // Check if octants node is a boundary
4628  if (oct->info[iface1] == false && oct->info[iface2] == false && oct->info[iface3] == false){
4629 
4630  //Build Morton number of virtual neigh of same size
4631  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
4632  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
4633 
4634  //SEARCH IN GHOSTS
4635 
4636  if (ghosts.size()>0){
4637  // Search in ghosts
4638  uint32_t idxghost = uint32_t(size_ghosts/2);
4639  Class_Octant<3>* octghost = &ghosts[idxghost];
4640 
4641  // Search morton in octants
4642  // If a even face morton is lower than morton of oct, if odd higher
4643  // ---> can i search only before or after idx in octants
4644  int32_t jump = (octghost->computeMorton() > Morton) ? int32_t(idxghost/2+1) : int32_t((size_ghosts -idxghost)/2+1);
4645  idxtry = uint32_t(idxghost +((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
4646  if (idxtry > size_ghosts-1)
4647  idxtry = size_ghosts-1;
4648  while(abs(jump) > 0){
4649  Mortontry = ghosts[idxtry].computeMorton();
4650  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4651  idxtry += jump;
4652  if (idxtry > ghosts.size()-1){
4653  if (jump > 0){
4654  idxtry = ghosts.size() - 1;
4655  jump = 0;
4656  }
4657  else if (jump < 0){
4658  idxtry = 0;
4659  jump = 0;
4660  }
4661  }
4662  }
4663  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4664  //Found neighbour of same size
4665  isghost.push_back(true);
4666  neighbours.push_back(idxtry);
4667  return;
4668  }
4669  else{
4670  // Step until the mortontry lower than morton (one idx of distance)
4671  {
4672  while(ghosts[idxtry].computeMorton() < Morton){
4673  idxtry++;
4674  if(idxtry > ghosts.size()-1){
4675  idxtry = ghosts.size()-1;
4676  break;
4677  }
4678  }
4679  while(ghosts[idxtry].computeMorton() > Morton){
4680  idxtry--;
4681  if(idxtry > ghosts.size()-1){
4682  idxtry = 0;
4683  break;
4684  }
4685  }
4686  }
4687  if(idxtry < size_ghosts){
4688  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4689  //Found neighbour of same size
4690  isghost.push_back(true);
4691  neighbours.push_back(idxtry);
4692  return;
4693  }
4694  // Compute Last discendent of virtual octant of same size
4695  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4696  uint64_t Mortonlast = last_desc.computeMorton();
4697  Mortontry = ghosts[idxtry].computeMorton();
4698  while(Mortontry < Mortonlast && idxtry < size_ghosts){
4699  Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(ghosts[idxtry].x));
4700  Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(ghosts[idxtry].y));
4701  Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(ghosts[idxtry].z));
4702  Dhxref = int32_t(cx<0)*ghosts[idxtry].getSize() + int32_t(cx>0)*size;
4703  Dhyref = int32_t(cy<0)*ghosts[idxtry].getSize() + int32_t(cy>0)*size;
4704  Dhzref = int32_t(cz<0)*ghosts[idxtry].getSize() + int32_t(cz>0)*size;
4705  if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
4706  neighbours.push_back(idxtry);
4707  isghost.push_back(true);
4708  return;
4709  }
4710  idxtry++;
4711  if(idxtry>size_ghosts-1){
4712  break;
4713  }
4714  Mortontry = ghosts[idxtry].computeMorton();
4715  }
4716  }
4717  }
4718  }
4719 
4720  // Search in octants
4721 
4722  //Build Morton number of virtual neigh of same size
4723  // Search morton in octants
4724  // If a even face morton is lower than morton of oct, if odd higher
4725  // ---> can i search only before or after idx in octants
4726  int32_t jump = (int32_t((noctants)/2+1));
4727  idxtry = uint32_t(jump);
4728  if (idxtry > noctants-1)
4729  idxtry = noctants-1;
4730  while(abs(jump) > 0){
4731  Mortontry = octants[idxtry].computeMorton();
4732  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4733  idxtry += jump;
4734  if (idxtry > octants.size()-1){
4735  if (jump > 0){
4736  idxtry = octants.size() - 1;
4737  jump = 0;
4738  }
4739  else if (jump < 0){
4740  idxtry = 0;
4741  jump = 0;
4742  }
4743  }
4744  }
4745  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4746  //Found neighbour of same size
4747  isghost.push_back(false);
4748  neighbours.push_back(idxtry);
4749  return;
4750  }
4751  else{
4752  // Step until the mortontry lower than morton (one idx of distance)
4753  {
4754  while(octants[idxtry].computeMorton() < Morton){
4755  idxtry++;
4756  if(idxtry > noctants-1){
4757  idxtry = noctants-1;
4758  break;
4759  }
4760  }
4761  while(octants[idxtry].computeMorton() > Morton){
4762  idxtry--;
4763  if(idxtry > noctants-1){
4764  idxtry = 0;
4765  break;
4766  }
4767  }
4768  }
4769  if (idxtry < noctants){
4770  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4771  //Found neighbour of same size
4772  isghost.push_back(false);
4773  neighbours.push_back(idxtry);
4774  return;
4775  }
4776  // Compute Last discendent of virtual octant of same size
4777  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4778  uint64_t Mortonlast = last_desc.computeMorton();
4779  Mortontry = octants[idxtry].computeMorton();
4780  while(Mortontry < Mortonlast && idxtry <= noctants-1){
4781  Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(octants[idxtry].x));
4782  Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(octants[idxtry].y));
4783  Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(octants[idxtry].z));
4784  Dhxref = int32_t(cx<0)*octants[idxtry].getSize() + int32_t(cx>0)*size;
4785  Dhyref = int32_t(cy<0)*octants[idxtry].getSize() + int32_t(cy>0)*size;
4786  Dhzref = int32_t(cz<0)*octants[idxtry].getSize() + int32_t(cz>0)*size;
4787  if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
4788  neighbours.push_back(idxtry);
4789  isghost.push_back(false);
4790  return;
4791  }
4792  idxtry++;
4793  if(idxtry>noctants-1){
4794  break;
4795  }
4796  Mortontry = octants[idxtry].computeMorton();
4797  }
4798  }
4799  }
4800  return;
4801  }
4802  else{
4803  // Boundary Node
4804  return;
4805  }
4806  };
4807 
4808  // =================================================================================== //
4809 
4810  void findNodeNeighbours(uint32_t idx, // Finds neighbours of idx-th octant through inode in vector octants.
4811  uint8_t inode, // Returns a vector (empty if inode is a bound node) with the index of neighbours
4812  u32vector & neighbours, // in their structure (octants or ghosts) and sets isghost[i] = true if the
4813  vector<bool> & isghost){ // i-th neighbour is ghost in the local tree
4814 
4815  uint64_t Morton, Mortontry;
4816  uint32_t noctants = getNumOctants();
4817  uint32_t idxtry;
4818  Class_Octant<3>* oct = &octants[idx];
4819  uint32_t size = oct->getSize();
4820  uint8_t iface1, iface2, iface3;
4821  int32_t Dhx, Dhy, Dhz;
4822  int32_t Dhxref, Dhyref, Dhzref;
4823 
4824  //Alternative to switch case
4825  int8_t Cx[8] = {-1,1,-1,1,-1,1,-1,1};
4826  int8_t Cy[8] = {-1,-1,1,1,-1,-1,1,1};
4827  int8_t Cz[8] = {-1,-1,-1,-1,1,1,1,1};
4828  int8_t cx = Cx[inode];
4829  int8_t cy = Cy[inode];
4830  int8_t cz = Cz[inode];
4831 
4832  isghost.clear();
4833  neighbours.clear();
4834 
4835  // Default if inode is nnodes<inode<0
4836  if (inode < 0 || inode > global3D.nnodes){
4837  return;
4838  }
4839 
4840  // Check if octants node is a boundary
4841  iface1 = global3D.nodeface[inode][0];
4842  iface2 = global3D.nodeface[inode][1];
4843  iface3 = global3D.nodeface[inode][2];
4844 
4845  // Check if octants node is a boundary
4846  if (oct->info[iface1] == false && oct->info[iface2] == false && oct->info[iface3] == false){
4847 
4848  //Build Morton number of virtual neigh of same size
4849  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
4850  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
4851 
4852  //SEARCH IN GHOSTS
4853 
4854  if (ghosts.size()>0){
4855  // Search in ghosts
4856  uint32_t idxghost = uint32_t(size_ghosts/2);
4857  Class_Octant<3>* octghost = &ghosts[idxghost];
4858 
4859  // Search morton in octants
4860  // If a even face morton is lower than morton of oct, if odd higher
4861  // ---> can i search only before or after idx in octants
4862  int32_t jump = (octghost->computeMorton() > Morton) ? int32_t(idxghost/2+1) : int32_t((size_ghosts -idxghost)/2+1);
4863  idxtry = uint32_t(idxghost +((octghost->computeMorton()<Morton)-(octghost->computeMorton()>Morton))*jump);
4864  if (idxtry > size_ghosts-1)
4865  idxtry = size_ghosts-1;
4866  while(abs(jump) > 0){
4867  Mortontry = ghosts[idxtry].computeMorton();
4868  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4869  idxtry += jump;
4870  if (idxtry > ghosts.size()-1){
4871  if (jump > 0){
4872  idxtry = ghosts.size() - 1;
4873  jump = 0;
4874  }
4875  else if (jump < 0){
4876  idxtry = 0;
4877  jump = 0;
4878  }
4879  }
4880  }
4881  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4882  //Found neighbour of same size
4883  isghost.push_back(true);
4884  neighbours.push_back(idxtry);
4885  return;
4886  }
4887  else{
4888  // Step until the mortontry lower than morton (one idx of distance)
4889  {
4890  while(ghosts[idxtry].computeMorton() < Morton){
4891  idxtry++;
4892  if(idxtry > ghosts.size()-1){
4893  idxtry = ghosts.size()-1;
4894  break;
4895  }
4896  }
4897  while(ghosts[idxtry].computeMorton() > Morton){
4898  idxtry--;
4899  if(idxtry > ghosts.size()-1){
4900  idxtry = 0;
4901  break;
4902  }
4903  }
4904  }
4905  if(idxtry < size_ghosts){
4906  if(ghosts[idxtry].computeMorton() == Morton && ghosts[idxtry].level == oct->level){
4907  //Found neighbour of same size
4908  isghost.push_back(true);
4909  neighbours.push_back(idxtry);
4910  return;
4911  }
4912  // Compute Last discendent of virtual octant of same size
4913  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4914  uint64_t Mortonlast = last_desc.computeMorton();
4915  Mortontry = ghosts[idxtry].computeMorton();
4916  while(Mortontry < Mortonlast && idxtry < size_ghosts){
4917  Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(ghosts[idxtry].x));
4918  Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(ghosts[idxtry].y));
4919  Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(ghosts[idxtry].z));
4920  Dhxref = int32_t(cx<0)*ghosts[idxtry].getSize() + int32_t(cx>0)*size;
4921  Dhyref = int32_t(cy<0)*ghosts[idxtry].getSize() + int32_t(cy>0)*size;
4922  Dhzref = int32_t(cz<0)*ghosts[idxtry].getSize() + int32_t(cz>0)*size;
4923  if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
4924  neighbours.push_back(idxtry);
4925  isghost.push_back(true);
4926  return;
4927  }
4928  idxtry++;
4929  if(idxtry>size_ghosts-1){
4930  break;
4931  }
4932  Mortontry = ghosts[idxtry].computeMorton();
4933  }
4934  }
4935  }
4936  }
4937 
4938  // Search in octants
4939 
4940  //Build Morton number of virtual neigh of same size
4941  // Search morton in octants
4942  // If a even face morton is lower than morton of oct, if odd higher
4943  // ---> can i search only before or after idx in octants
4944  int32_t jump = (int32_t((noctants)/2+1));
4945  idxtry = uint32_t(jump);
4946  if (idxtry > noctants-1)
4947  idxtry = noctants-1;
4948  while(abs(jump) > 0){
4949  Mortontry = octants[idxtry].computeMorton();
4950  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
4951  idxtry += jump;
4952  if (idxtry > octants.size()-1){
4953  if (jump > 0){
4954  idxtry = octants.size() - 1;
4955  jump = 0;
4956  }
4957  else if (jump < 0){
4958  idxtry = 0;
4959  jump = 0;
4960  }
4961  }
4962  }
4963  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4964  //Found neighbour of same size
4965  isghost.push_back(false);
4966  neighbours.push_back(idxtry);
4967  return;
4968  }
4969  else{
4970  // Step until the mortontry lower than morton (one idx of distance)
4971  {
4972  while(octants[idxtry].computeMorton() < Morton){
4973  idxtry++;
4974  if(idxtry > noctants-1){
4975  idxtry = noctants-1;
4976  break;
4977  }
4978  }
4979  while(octants[idxtry].computeMorton() > Morton){
4980  idxtry--;
4981  if(idxtry > noctants-1){
4982  idxtry = 0;
4983  break;
4984  }
4985  }
4986  }
4987  if (idxtry < noctants){
4988  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
4989  //Found neighbour of same size
4990  isghost.push_back(false);
4991  neighbours.push_back(idxtry);
4992  return;
4993  }
4994  // Compute Last discendent of virtual octant of same size
4995  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
4996  uint64_t Mortonlast = last_desc.computeMorton();
4997  Mortontry = octants[idxtry].computeMorton();
4998  while(Mortontry < Mortonlast && idxtry <= noctants-1){
4999  Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(octants[idxtry].x));
5000  Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(octants[idxtry].y));
5001  Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(octants[idxtry].z));
5002  Dhxref = int32_t(cx<0)*octants[idxtry].getSize() + int32_t(cx>0)*size;
5003  Dhyref = int32_t(cy<0)*octants[idxtry].getSize() + int32_t(cy>0)*size;
5004  Dhzref = int32_t(cz<0)*octants[idxtry].getSize() + int32_t(cz>0)*size;
5005  if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
5006  neighbours.push_back(idxtry);
5007  isghost.push_back(false);
5008  return;
5009  }
5010  idxtry++;
5011  if(idxtry>noctants-1){
5012  break;
5013  }
5014  Mortontry = octants[idxtry].computeMorton();
5015  }
5016  }
5017  }
5018  return;
5019  }
5020  else{
5021  // Boundary Node
5022  return;
5023  }
5024  };
5025 
5026  // =================================================================================== //
5027 
5028 void findGhostNodeNeighbours(uint32_t idx, // Finds neighbours of idx-th ghost through inode in vector octants.
5029  uint8_t inode, // Returns a vector (empty if inode is not a pbound node) with the index of neighbours
5030  u32vector & neighbours){ // in the structure octants.
5031 
5032  uint64_t Morton, Mortontry;
5033  uint32_t noctants = getNumOctants();
5034  uint32_t idxtry;
5035  Class_Octant<3>* oct = &ghosts[idx];
5036  uint32_t size = oct->getSize();
5037  uint8_t iface1, iface2, iface3;
5038  int32_t Dhx, Dhy, Dhz;
5039  int32_t Dhxref, Dhyref, Dhzref;
5040 
5041  //Alternative to switch case
5042  int8_t Cx[8] = {-1,1,-1,1,-1,1,-1,1};
5043  int8_t Cy[8] = {-1,-1,1,1,-1,-1,1,1};
5044  int8_t Cz[8] = {-1,-1,-1,-1,1,1,1,1};
5045  int8_t cx = Cx[inode];
5046  int8_t cy = Cy[inode];
5047  int8_t cz = Cz[inode];
5048 
5049  neighbours.clear();
5050 
5051  // Default if inode is nnodes<inode<0
5052  if (inode < 0 || inode > global3D.nnodes){
5053  return;
5054  }
5055 
5056  // Check if octants node is a boundary
5057  iface1 = global3D.nodeface[inode][0];
5058  iface2 = global3D.nodeface[inode][1];
5059  iface3 = global3D.nodeface[inode][2];
5060 
5061 // // Check if octants node is a boundary
5062 // if (oct->info[iface1] == false && oct->info[iface2] == false && oct->info[iface3] == false){
5063  // Check if octants node is a pboundary node
5064  if (oct->info[iface1+global3D.nfaces] == true || oct->info[iface2+global3D.nfaces] == true || oct->info[iface3+global3D.nfaces] == true){
5065 
5066  //Build Morton number of virtual neigh of same size
5067  Class_Octant<3> samesizeoct(oct->level, oct->x+cx*size, oct->y+cy*size, oct->z+cz*size);
5068  Morton = samesizeoct.computeMorton(); //mortonEncode_magicbits(oct->x-size,oct->y,oct->z);
5069  int32_t jump = noctants/2;
5070  idxtry = jump;
5071  while(abs(jump) > 0){
5072  Mortontry = octants[idxtry].computeMorton();
5073  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
5074  idxtry += jump;
5075  }
5076  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
5077  //Found neighbour of same size
5078  neighbours.push_back(idxtry);
5079  return;
5080  }
5081  else{
5082  // Step until the mortontry lower than morton (one idx of distance)
5083  {
5084  while(octants[idxtry].computeMorton() < Morton){
5085  idxtry++;
5086  if(idxtry > noctants-1){
5087  idxtry = noctants-1;
5088  break;
5089  }
5090  }
5091  while(octants[idxtry].computeMorton() > Morton){
5092  idxtry--;
5093  if(idxtry > noctants-1){
5094  idxtry = 0;
5095  break;
5096  }
5097  }
5098  }
5099  if (idxtry < noctants){
5100  if(octants[idxtry].computeMorton() == Morton && octants[idxtry].level == oct->level){
5101  //Found neighbour of same size
5102  neighbours.push_back(idxtry);
5103  return;
5104  }
5105  // Compute Last discendent of virtual octant of same size
5106  Class_Octant<3> last_desc = samesizeoct.buildLastDesc();
5107  uint64_t Mortonlast = last_desc.computeMorton();
5108  Mortontry = octants[idxtry].computeMorton();
5109  while(Mortontry < Mortonlast && idxtry <= noctants-1){
5110  Dhx = int32_t(cx)*(int32_t(oct->x) - int32_t(octants[idxtry].x));
5111  Dhy = int32_t(cy)*(int32_t(oct->y) - int32_t(octants[idxtry].y));
5112  Dhz = int32_t(cz)*(int32_t(oct->z) - int32_t(octants[idxtry].z));
5113  Dhxref = int32_t(cx<0)*octants[idxtry].getSize() + int32_t(cx>0)*size;
5114  Dhyref = int32_t(cy<0)*octants[idxtry].getSize() + int32_t(cy>0)*size;
5115  Dhzref = int32_t(cz<0)*octants[idxtry].getSize() + int32_t(cz>0)*size;
5116  if ((abs(Dhx) == Dhxref) && (abs(Dhy) == Dhyref) && (abs(Dhz) == Dhzref)){
5117  neighbours.push_back(idxtry);
5118  }
5119  idxtry++;
5120  Mortontry = octants[idxtry].computeMorton();
5121  }
5122  }
5123  }
5124  return;
5125  }
5126  else{
5127  // Boundary Node
5128  return;
5129  }
5130 };
5131 
5132 // =================================================================================== //
5133 
5134 void computeIntersections() {
5135 
5136  OctantsType::iterator it, obegin, oend;
5137  Class_Intersection<3> intersection;
5138  u32vector neighbours;
5139  vector<bool> isghost;
5140  uint32_t counter, idx;
5141  uint32_t i, nsize;
5142  uint8_t iface, iface2;
5143 
5144  intersections.clear();
5145  intersections.reserve(2*3*octants.size());
5146 
5147  counter = idx = 0;
5148 
5149  // Loop on ghosts
5150  obegin = ghosts.begin();
5151  oend = ghosts.end();
5152  for (it = obegin; it != oend; it++){
5153  for (iface = 0; iface < 3; iface++){
5154  iface2 = iface*2;
5155  findGhostNeighbours(idx, iface2, neighbours);
5156  nsize = neighbours.size();
5157  for (i = 0; i < nsize; i++){
5158  intersection.finer = getGhostLevel(idx) >= getLevel((int)neighbours[i]);
5159  intersection.owners[0] = neighbours[i];
5160  intersection.owners[1] = idx;
5161  intersection.iface = global3D.oppface[iface2] - (getGhostLevel(idx) >= getLevel((int)neighbours[i]));
5162  intersection.isnew = false;
5163  intersection.isghost = true;
5164  intersection.bound = false;
5165  intersection.pbound = true;
5166  intersections.push_back(intersection);
5167  counter++;
5168  }
5169  }
5170  idx++;
5171  }
5172 
5173  // Loop on octants
5174  idx=0;
5175  obegin = octants.begin();
5176  oend = octants.end();
5177  for (it = obegin; it != oend; it++){
5178  for (iface = 0; iface < 3; iface++){
5179  iface2 = iface*2;
5180  findNeighbours(idx, iface2, neighbours, isghost);
5181  nsize = neighbours.size();
5182  if (nsize) {
5183  for (i = 0; i < nsize; i++){
5184  if (isghost[i]){
5185  intersection.owners[0] = idx;
5186  intersection.owners[1] = neighbours[i];
5187  intersection.finer = (nsize>1);
5188  intersection.iface = iface2 + (nsize>1);
5189  intersection.isnew = false;
5190  intersection.isghost = true;
5191  intersection.bound = false;
5192  intersection.pbound = true;
5193  intersections.push_back(intersection);
5194  counter++;
5195  }
5196  else{
5197  intersection.owners[0] = idx;
5198  intersection.owners[1] = neighbours[i];
5199  intersection.finer = (nsize>1);
5200  intersection.iface = iface2 + (nsize>1);
5201  intersection.isnew = false;
5202  intersection.isghost = false;
5203  intersection.bound = false;
5204  intersection.pbound = false;
5205  intersections.push_back(intersection);
5206  counter++;
5207  }
5208  }
5209  }
5210  else{
5211  intersection.owners[0] = idx;
5212  intersection.owners[1] = idx;
5213  intersection.finer = 0;
5214  intersection.iface = iface2;
5215  intersection.isnew = false;
5216  intersection.isghost = false;
5217  intersection.bound = true;
5218  intersection.pbound = false;
5219  intersections.push_back(intersection);
5220  counter++;
5221  }
5222  if (it->info[iface2+1]){
5223  intersection.owners[0] = idx;
5224  intersection.owners[1] = idx;
5225  intersection.finer = 0;
5226  intersection.iface = iface2+1;
5227  intersection.isnew = false;
5228  intersection.isghost = false;
5229  intersection.bound = true;
5230  intersection.pbound = false;
5231  intersections.push_back(intersection);
5232  counter++;
5233  }
5234  }
5235  idx++;
5236  }
5237 #if defined(__INTEL_COMPILER) || defined(__ICC)
5238 #else
5239  intersections.shrink_to_fit();
5240 #endif
5241  }
5242 
5243  // =================================================================================== //
5244 
5245  uint32_t findMorton(uint64_t Morton){ // Find an input Morton in octants and return the local idx
5246  uint32_t nocts = octants.size();
5247  uint32_t idx = nocts/2;
5248  uint64_t Mortontry = octants[idx].computeMorton();
5249  int32_t jump = nocts/2;
5250  while(abs(jump)>0){
5251  if (Mortontry == Morton){
5252  return idx;
5253  }
5254  Mortontry = octants[idx].computeMorton();
5255  jump = ((Mortontry<Morton)-(Mortontry>Morton))*abs(jump)/2;
5256  idx += jump;
5257  if (idx > nocts){
5258  return nocts-1; // return nocts if not found the Morton
5259  }
5260  }
5261  if (Mortontry<Morton){
5262  for (uint32_t idx2=idx; idx2<nocts; idx2++){
5263  Mortontry = octants[idx2].computeMorton();
5264  if (Mortontry == Morton){
5265  return idx2;
5266  }
5267  }
5268  }
5269  else{
5270  for(uint32_t idx2=0; idx2<idx+1; idx2++){
5271  Mortontry = octants[idx2].computeMorton();
5272  if (Mortontry == Morton){
5273  return idx2;
5274  }
5275  }
5276  }
5277  return nocts-1;
5278  };
5279 
5280  // =================================================================================== //
5281 
5282  uint32_t findGhostMorton(uint64_t Morton){ // Find an input Morton in ghosts and return the local idx
5283  uint32_t nocts = ghosts.size();
5284  uint32_t idx = nocts/2;
5285  uint64_t Mortontry = ghosts[idx].computeMorton();
5286  int32_t jump = nocts/2;
5287  while(abs(jump)>0){
5288  if (Mortontry == Morton){
5289  return idx;
5290  }
5291  Mortontry = ghosts[idx].computeMorton();
5292  jump = (Mortontry<Morton)*jump/4;
5293  idx += jump;
5294  if (idx > nocts){
5295  return nocts; // return nocts if not found the Morton
5296  }
5297  }
5298  if (Mortontry<Morton){
5299  for (uint32_t idx2=idx; idx2<nocts; idx2++){
5300  Mortontry = ghosts[idx2].computeMorton();
5301  if (Mortontry == Morton){
5302  return idx2;
5303  }
5304  }
5305  }
5306  else{
5307  for(uint32_t idx2=0; idx2<idx; idx2++){
5308  Mortontry = ghosts[idx2].computeMorton();
5309  if (Mortontry == Morton){
5310  return idx2;
5311  }
5312  }
5313  }
5314  return nocts;
5315  };
5316 
5317  // =================================================================================== //
5318 
5321  void computeConnectivity(){ // Computes nodes vector and connectivity of octants of local tree
5322  map<uint64_t, vector<uint32_t> > mapnodes;
5323  map<uint64_t, vector<uint32_t> >::iterator iter, iterend;
5324  uint32_t i, k, counter;
5325  uint64_t morton;
5326  uint32_t noctants = getNumOctants();
5327  u32vector2D octnodes;
5328  uint8_t j;
5329 
5330  clearConnectivity();
5331  octnodes.reserve(global3D.nnodes);
5332 
5333  if (nodes.size() == 0){
5334  connectivity.resize(noctants);
5335  for (i = 0; i < noctants; i++){
5336  octants[i].getNodes(octnodes);
5337  for (j = 0; j < global3D.nnodes; j++){
5338 // morton = mortonEncode_magicbits(octnodes[j][0], octnodes[j][1], octnodes[j][2]);
5339 /*debug version*/ morton = keyXYZ(octnodes[j][0], octnodes[j][1], octnodes[j][2]);
5340  if (mapnodes[morton].size()==0){
5341  mapnodes[morton].reserve(16);
5342  for (k = 0; k < 3; k++){
5343  mapnodes[morton].push_back(octnodes[j][k]);
5344  }
5345  }
5346  mapnodes[morton].push_back(i);
5347  }
5348  u32vector2D().swap(octnodes);
5349  }
5350  iter = mapnodes.begin();
5351  iterend = mapnodes.end();
5352  counter = 0;
5353  uint32_t numnodes = mapnodes.size();
5354  nodes.resize(numnodes);
5355  while (iter != iterend){
5356  vector<uint32_t> nodecasting(iter->second.begin(), iter->second.begin()+3);
5357  nodes[counter] = nodecasting;
5358 #if defined(__INTEL_COMPILER) || defined(__ICC)
5359 #else
5360  nodes[counter].shrink_to_fit();
5361 #endif
5362  for(vector<uint32_t>::iterator iter2 = iter->second.begin()+3; iter2 != iter->second.end(); iter2++){
5363  if (connectivity[(*iter2)].size()==0){
5364  connectivity[(*iter2)].reserve(8);
5365  }
5366  connectivity[(*iter2)].push_back(counter);
5367  }
5368  mapnodes.erase(iter++);
5369  counter++;
5370  }
5371 #if defined(__INTEL_COMPILER) || defined(__ICC)
5372 #else
5373  nodes.shrink_to_fit();
5374 #endif
5375  //Lento. Solo per risparmiare memoria
5376  for (uint32_t ii=0; ii<noctants; ii++){
5377 #if defined(__INTEL_COMPILER) || defined(__ICC)
5378 #else
5379  connectivity[ii].shrink_to_fit();
5380 #endif
5381  }
5382 #if defined(__INTEL_COMPILER) || defined(__ICC)
5383 #else
5384  connectivity.shrink_to_fit();
5385 #endif
5386  }
5387  map<uint64_t, vector<uint32_t> >().swap(mapnodes);
5388  iter = mapnodes.end();
5389 
5390  };
5391 
5392  // =================================================================================== //
5393 
5394  void clearConnectivity(){ // Clear nodes vector and connectivity of octants of local tree
5395  u32vector2D().swap(nodes);
5396  u32vector2D().swap(connectivity);
5397  };
5398 
5399  // =================================================================================== //
5400 
5401  void updateConnectivity(){ // Updates nodes vector and connectivity of octants of local tree
5402  clearConnectivity();
5403  computeConnectivity();
5404  };
5405 
5406  // =================================================================================== //
5407 
5408  void computeGhostsConnectivity(){ // Computes ghosts nodes vector and connectivity of ghosts octants of local tree
5409  map<uint64_t, vector<uint32_t> > mapnodes;
5410  map<uint64_t, vector<uint32_t> >::iterator iter, iterend;
5411  uint32_t i, k, counter;
5412  uint64_t morton;
5413  uint32_t noctants = size_ghosts;
5414  u32vector2D octnodes;
5415  uint8_t j;
5416 
5417  octnodes.reserve(global3D.nnodes);
5418  if (ghostsnodes.size() == 0){
5419  ghostsconnectivity.resize(noctants);
5420  for (i = 0; i < noctants; i++){
5421  ghosts[i].getNodes(octnodes);
5422  for (j = 0; j < global3D.nnodes; j++){
5423 // morton = mortonEncode_magicbits(octnodes[j][0], octnodes[j][1], octnodes[j][2]);
5424 /*debug version*/ morton = keyXYZ(octnodes[j][0], octnodes[j][1], octnodes[j][2]);
5425  if (mapnodes[morton].size()==0){
5426  for (k = 0; k < 3; k++){
5427  mapnodes[morton].push_back(octnodes[j][k]);
5428  }
5429  }
5430  mapnodes[morton].push_back(i);
5431  }
5432  u32vector2D().swap(octnodes);
5433  }
5434  iter = mapnodes.begin();
5435  iterend = mapnodes.end();
5436  uint32_t numnodes = mapnodes.size();
5437  ghostsnodes.resize(numnodes);
5438  counter = 0;
5439  while (iter != iterend){
5440  vector<uint32_t> nodecasting(iter->second.begin(), iter->second.begin()+3);
5441  ghostsnodes[counter] = nodecasting;
5442 #if defined(__INTEL_COMPILER) || defined(__ICC)
5443 #else
5444  ghostsnodes[counter].shrink_to_fit();
5445 #endif
5446  for(vector<uint32_t>::iterator iter2 = iter->second.begin()+3; iter2 != iter->second.end(); iter2++){
5447  if (ghostsconnectivity[(*iter2)].size()==0){
5448  ghostsconnectivity[(*iter2)].reserve(8);
5449  }
5450  ghostsconnectivity[(*iter2)].push_back(counter);
5451  }
5452  mapnodes.erase(iter++);
5453  counter++;
5454  }
5455 #if defined(__INTEL_COMPILER) || defined(__ICC)
5456 #else
5457  ghostsnodes.shrink_to_fit();
5458 #endif
5459  //Lento. Solo per risparmiare memoria
5460  for (uint32_t ii=0; ii<noctants; ii++){
5461 #if defined(__INTEL_COMPILER) || defined(__ICC)
5462 #else
5463  ghostsconnectivity[ii].shrink_to_fit();
5464 #endif
5465  }
5466 #if defined(__INTEL_COMPILER) || defined(__ICC)
5467 #else
5468  ghostsconnectivity.shrink_to_fit();
5469 #endif
5470  }
5471  iter = mapnodes.end();
5472 
5473  };
5474 
5475  // =================================================================================== //
5476 
5477  void clearGhostsConnectivity(){ // Clear ghosts nodes vector and connectivity of ghosts octants of local tree
5478  u32vector2D().swap(ghostsnodes);
5479  u32vector2D().swap(ghostsconnectivity);
5480  };
5481 
5482  // =================================================================================== //
5483 
5484  void updateGhostsConnectivity(){ // Update ghosts nodes vector and connectivity of ghosts octants of local tree
5485  clearGhostsConnectivity();
5486  computeGhostsConnectivity();
5487  };
5488 
5489  // =================================================================================== //
5490 
5491 };//end Class_Local_Tree<3> specialization;
5492 
5493 
Class_Octant< 3 >::computeMorton
uint64_t computeMorton() const
Definition: Class_Octant_3D.tpp:407
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_Global< 3 >::oppface
uint8_t oppface[6]
Definition: Class_Global_3D.tpp:190
Class_Global< 3 >::nedges
uint8_t nedges
Definition: Class_Global_3D.tpp:182
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 >::edgecoeffs
int8_t edgecoeffs[12][3]
Definition: Class_Global_3D.tpp:195
Class_Global< 3 >::nodeface
uint8_t nodeface[8][3]
Definition: Class_Global_3D.tpp:191
Class_Global< 3 >::max_length
uint32_t max_length
Definition: Class_Global_3D.tpp:176
Class_Global< 3 >::nfaces
uint8_t nfaces
Definition: Class_Global_3D.tpp:181
Class_Octant< 3 >
Octant class definition - 3D specialization.
Definition: Class_Octant_3D.tpp:32
Class_Intersection< 3 >
Intersection class definition - 3D specialization.
Definition: Class_Intersection_3D.tpp:27
Class_Octant< 3 >::getMarker
int8_t getMarker() const
Definition: Class_Octant_3D.tpp:154
Class_Local_Tree
Local octree portion for each process.
Definition: Class_Local_Tree.hpp:55
Class_Global< 3 >::normals
int8_t normals[6][3]
Definition: Class_Global_3D.tpp:194
Class_Octant< 3 >::setMarker
void setMarker(int8_t marker)
Definition: Class_Octant_3D.tpp:201
Class_Global< 3 >::edgeface
uint8_t edgeface[12][2]
Definition: Class_Global_3D.tpp:193
Generated on Mon Oct 5 2015 15:12:11 for PABLO by   doxygen 1.8.8