Loading...
Searching...
No Matches
LocalTree.cpp
1/*---------------------------------------------------------------------------*\
2 *
3 * bitpit
4 *
5 * Copyright (C) 2015-2021 OPTIMAD engineering Srl
6 *
7 * -------------------------------------------------------------------------
8 * License
9 * This file is part of bitpit.
10 *
11 * bitpit is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU Lesser General Public License v3 (LGPL)
13 * as published by the Free Software Foundation.
14 *
15 * bitpit is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public License
21 * along with bitpit. If not, see <http://www.gnu.org/licenses/>.
22 *
23 \*---------------------------------------------------------------------------*/
24
25// =================================================================================== //
26// INCLUDES //
27// =================================================================================== //
28#include "LocalTree.hpp"
29#include "morton.hpp"
30
31#include <map>
32#include <unordered_map>
33
34namespace bitpit {
35
36 // =================================================================================== //
37 // NAME SPACES //
38 // =================================================================================== //
39 using namespace std;
40
41 // =================================================================================== //
42 // CLASS IMPLEMENTATION //
43 // =================================================================================== //
44
45 // =================================================================================== //
46 // CONSTRUCTORS AND OPERATORS
47 // =================================================================================== //
48
51 LocalTree::LocalTree() {
52 initialize();
53 reset(false);
54 };
55
59 LocalTree::LocalTree(uint8_t dim){
60 initialize(dim);
61 reset(true);
62 };
63
64 // =================================================================================== //
65 // METHODS
66 // =================================================================================== //
67
68 // =================================================================================== //
69 // BASIC GET/SET METHODS
70 // =================================================================================== //
71
75 uint64_t
76 LocalTree::getFirstDescMorton() const{
77 return m_firstDescMorton;
78 };
79
83 uint64_t
84 LocalTree::getLastDescMorton() const{
85 return m_lastDescMorton;
86 };
87
91 uint32_t
92 LocalTree::getNumGhosts() const{
93 return m_ghosts.size();
94 };
95
99 uint32_t
100 LocalTree::getNumOctants() const{
101 return m_octants.size();
102 };
103
108 int8_t
109 LocalTree::getLocalMaxDepth() const{
110 return m_localMaxDepth;
111 };
112
117 int8_t
118 LocalTree::getMarker(int32_t idx) const {
119 return m_octants[idx].getMarker();
120 };
121
126 uint8_t
127 LocalTree::getLevel(int32_t idx) const {
128 return m_octants[idx].getLevel();
129 };
130
135 uint64_t
136 LocalTree::getMorton(int32_t idx) const {
137 return m_octants[idx].getMorton();
138 };
139
146 uint64_t
147 LocalTree::computeNodePersistentKey(int32_t idx, uint8_t inode) const {
148 return m_octants[idx].computeNodePersistentKey(inode);
149 };
150
155 uint8_t
156 LocalTree::getGhostLevel(int32_t idx) const {
157 return m_ghosts[idx].getLevel();
158 };
159
164 uint64_t
165 LocalTree::computeGhostMorton(int32_t idx) const {
166 return m_ghosts[idx].getMorton();
167 };
168
175 uint64_t
176 LocalTree::computeGhostNodePersistentKey(int32_t idx, uint8_t inode) const {
177 return m_ghosts[idx].computeNodePersistentKey(inode);
178 };
179
184 bool
185 LocalTree::getBalance(int32_t idx) const{
186 return m_octants[idx].getBalance();
187 };
188
192 uint8_t
193 LocalTree::getBalanceCodim() const{
194 return m_balanceCodim;
195 };
196
201 void
202 LocalTree::setMarker(int32_t idx, int8_t marker){
203 m_octants[idx].setMarker(marker);
204 };
205
210 void
211 LocalTree::setBalance(int32_t idx, bool balance){
212 m_octants[idx].setBalance(balance);
213 };
214
218 void
219 LocalTree::setBalanceCodim(uint8_t b21codim){
220 m_balanceCodim = b21codim;
221 };
222
225 void
226 LocalTree::setFirstDescMorton(){
227 if(getNumOctants() > 0){
228 octvector::const_iterator firstOctant = m_octants.begin();
229 m_firstDescMorton = firstOctant->getMorton();
230 } else {
231 m_firstDescMorton = std::numeric_limits<uint64_t>::max();
232 }
233 };
234
237 void
238 LocalTree::setLastDescMorton(){
239 if(getNumOctants() > 0){
240 octvector::const_iterator lastOctant = m_octants.end() - 1;
241 uint32_t x,y,z,delta;
242 delta = (uint32_t)(1<<((uint8_t)m_treeConstants->maxLevel - lastOctant->m_level)) - 1;
243 x = lastOctant->getLogicalCoordinates(0) + delta;
244 y = lastOctant->getLogicalCoordinates(1) + delta;
245 z = lastOctant->getLogicalCoordinates(2) + (m_dim-2)*delta;
246 Octant lastDesc = Octant(m_dim, m_treeConstants->maxLevel,x,y,z);
247 m_lastDescMorton = lastDesc.getMorton();
248 } else {
249 m_lastDescMorton = 0;
250 }
251 };
252
256 void
257 LocalTree::setPeriodic(bvector & periodic){
258 m_periodic = periodic;
259 };
260
261 // =================================================================================== //
262 // OTHER GET/SET METHODS
263 // =================================================================================== //
264
271 bool
272 LocalTree::isPeriodic(const Octant* oct, uint8_t iface) const{
273 // Check if face is on a boundary
274 if (!oct->getBound(iface)) {
275 return false;
276 }
277
278 // Check if boundary is periodic
279 return m_periodic[iface];
280 };
281
288 bool
289 LocalTree::isEdgePeriodic(const Octant* oct, uint8_t iedge) const{
290 // The edge needs to be on a border
291 if (!oct->getEdgeBound(iedge)) {
292 return false;
293 }
294
295 // If one of the edge faces is on a non-periodic border the edge is
296 // non-periodic
297 for (int face : m_treeConstants->edgeFace[iedge]) {
298 if (oct->getBound(face) && !isPeriodic(oct, face)) {
299 return false;
300 }
301 }
302
303 // The edge is periodic
304 return true;
305 };
306
313 bool
314 LocalTree::isNodePeriodic(const Octant* oct, uint8_t inode) const{
315 // The node needs to be on a border
316 if (!oct->getNodeBound(inode)) {
317 return false;
318 }
319
320 // If one of the node faces is on a non-periodic border the node is
321 // non-periodic
322 for (int face : m_treeConstants->nodeFace[inode]) {
323 if (oct->getBound(face) && !isPeriodic(oct, face)) {
324 return false;
325 }
326 }
327
328 // The edge is periodic
329 return true;
330 };
331
332 // =================================================================================== //
333 // OTHER METHODS
334 // =================================================================================== //
335
338 void
339 LocalTree::initialize() {
340 initialize(0);
341 }
342
346 void
347 LocalTree::initialize(uint8_t dim) {
348 m_dim = dim;
349 m_balanceCodim = 1;
350
351 m_periodic.resize(m_dim*2);
352
353 if (m_dim > 0) {
354 m_treeConstants = &(TreeConstants::instance(m_dim));
355 } else {
356 m_treeConstants = nullptr;
357 }
358 }
359
362 void
363 LocalTree::reset(bool createRoot){
364 m_octants.clear();
365 m_ghosts.clear();
366 m_globalIdxGhosts.clear();
367
368 m_lastGhostBros.clear();
369 m_firstGhostBros.clear();
370
371 clearConnectivity();
372 intervector().swap(m_intersections);
373
374 std::fill(m_periodic.begin(), m_periodic.end(), false);
375
376 if (m_dim > 0 && createRoot) {
377 m_localMaxDepth = 0;
378
379 m_octants.push_back(Octant(m_dim));
380 } else {
381 m_localMaxDepth = -1;
382 }
383
384 setFirstDescMorton();
385 setLastDescMorton();
386
387 };
388
393 Octant&
394 LocalTree::extractOctant(uint32_t idx){
395 return m_octants[idx];
396 };
397
402 const Octant&
403 LocalTree::extractOctant(uint32_t idx) const{
404 return m_octants[idx];
405 };
406
411 Octant&
412 LocalTree::extractGhostOctant(uint32_t idx) {
413 return m_ghosts[idx];
414 };
415
420 const Octant&
421 LocalTree::extractGhostOctant(uint32_t idx) const{
422 return m_ghosts[idx];
423 };
424
425 // =================================================================================== //
426
437 bool
438 LocalTree::refine(u32vector & mapidx){
439 // Current number of octants
440 uint32_t nOctants = m_octants.size();
441 if (nOctants == 0) {
442 return false;
443 }
444
445 // Validate markers
446 //
447 // Not all octants marked for refinement can be really refined, for
448 // example octants cannot be refined further than the maximum level.
449 uint32_t firstRefinedIdx = 0;
450 uint32_t nValidRefinements = 0;
451 for (uint32_t idx=0; idx<nOctants; idx++) {
452 // Skip octants not marked for refinement
453 Octant &octant = m_octants[idx];
454 if(octant.getMarker()<= 0){
455 continue;
456 }
457
458 // Octants cannot be refined further than the maximum level
459 if(octant.getLevel()>=m_treeConstants->maxLevel){
460 octant.setMarker(0);
461 continue;
462 }
463
464 // The octant will be refined
465 ++nValidRefinements;
466 if (nValidRefinements == 1) {
467 firstRefinedIdx = idx;
468 }
469 }
470
471 // Early return if no octants need to be refined
472 if (nValidRefinements == 0) {
473 return false;
474 }
475
476 // Resize octant container
477 //
478 // We want to be sure the container capacity is equal to its size.
479 uint32_t nFutureOctants = nOctants + (m_treeConstants->nChildren - 1) * nValidRefinements;
480
481 m_octants.reserve(nFutureOctants);
482 m_octants.resize(nFutureOctants, Octant(m_dim));
483 m_octants.shrink_to_fit();
484
485 // Initialize mapping
486 if(!mapidx.empty()){
487 mapidx.resize(nFutureOctants);
488 }
489
490 // Refine the octants
491 //
492 // We process the octants backwards (starting from the back of their
493 // container):
494 // - if an octant doesn't need refinement, it will only be moved to
495 // its new position;
496 // - if an octant needs to be refined, it will be replaced with its
497 // children.
498 bool refinementCompleted = true;
499
500 uint32_t futureIdx = nFutureOctants - 1;
501 for (uint32_t n=0; n<(nOctants - firstRefinedIdx); ++n) {
502 uint32_t idx = nOctants - n - 1;
503 Octant &octant = m_octants[idx];
504 if(octant.getMarker()<=0){
505 // The octant is not refiend, we need to move it to its new
506 // position in the container.
507 if (futureIdx != idx) {
508 m_octants[futureIdx] = std::move(octant);
509 }
510
511 // Update the mapping
512 if(!mapidx.empty()){
513 mapidx[futureIdx] = mapidx[idx];
514 }
515
516 // Update future octant index
517 --futureIdx;
518 } else {
519 // Move original octant out of the container
520 Octant fatherOctant = std::move(octant);
521
522 // Create children
523 uint8_t nChildren = fatherOctant.countChildren();
524 uint32_t firstChildIdx = futureIdx - (nChildren - 1);
525 fatherOctant.buildChildren(m_octants.data() + firstChildIdx);
526
527 // Update the mapping
528 if(!mapidx.empty()){
529 for (uint8_t i=0; i<nChildren; i++){
530 mapidx[firstChildIdx + i] = mapidx[idx];
531 }
532 }
533
534 // Check if more refinement is needed to satisfy the markers
535 if (refinementCompleted) {
536 refinementCompleted = (m_octants[firstChildIdx].getMarker() <= 0);
537 }
538
539 // Update local max depth
540 uint8_t childrenLevel = m_octants[firstChildIdx].getLevel();
541 if (childrenLevel > m_localMaxDepth){
542 m_localMaxDepth = childrenLevel;
543 }
544
545 // Update future octant index
546 futureIdx -= nChildren;
547 }
548 }
549
550 return (!refinementCompleted);
551
552 };
553
554 // =================================================================================== //
560 bool
561 LocalTree::coarse(u32vector & mapidx){
562
563 u32vector first_child_index;
564 Octant father(m_dim);
565 uint32_t idx, idx2;
566 uint32_t offset;
567 uint32_t idx1_gh;
568 uint32_t idx2_gh;
569 uint32_t nidx;
570 uint32_t mapsize = mapidx.size();
571 int8_t markerfather, marker;
572 uint8_t nbro, nend, nstart;
573 uint8_t nchm1 = m_treeConstants->nChildren-1;
574 bool docoarse = false;
575 bool wstop = false;
576
577 //------------------------------------------ //
578 // Initialization
579
580 uint32_t nInitialOctants = getNumOctants();
581 if (nInitialOctants == 0) {
582 return false;
583 }
584
585 uint32_t nInitialGhosts = getNumGhosts();
586
587 nbro = nend = nstart = 0;
588 nidx = offset = 0;
589
590 idx2_gh = 0;
591 idx1_gh = nInitialGhosts - 1;
592
593 //------------------------------------------ //
594
595 // Set index for start and end check for ghosts
596 if (m_ghosts.size()){
597 bool check = true;
598 while(check){
599 check = idx1_gh < nInitialGhosts;
600 if (check){
601 check = m_ghosts[idx1_gh].getMorton() > m_firstDescMorton;
602 }
603 if (check) idx1_gh--;
604 }
605
606 check = true;
607 while(check){
608 check = idx2_gh < nInitialGhosts;
609 if (check){
610 check = m_ghosts[idx2_gh].getMorton() < m_lastDescMorton;
611 }
612 if (check) idx2_gh++;
613 }
614 }
615
616 // Check and coarse internal octants
617 for (idx=0; idx<nInitialOctants; idx++){
618 if(m_octants[idx].getMarker() < 0 && m_octants[idx].getLevel() > 0){
619 nbro = 0;
620 father = m_octants[idx].buildFather();
621 // Check if family is to be refined
622 for (idx2=idx; idx2<idx+m_treeConstants->nChildren; idx2++){
623 if (idx2<nInitialOctants){
624 if(m_octants[idx2].getMarker() < 0 && m_octants[idx2].buildFather() == father){
625 nbro++;
626 }
627 }
628 }
629 if (nbro == m_treeConstants->nChildren){
630 nidx++;
631 first_child_index.push_back(idx);
632 idx = idx2-1;
633 }
634 }
635 }
636 uint32_t nblock = nInitialOctants;
637 uint32_t nfchild = first_child_index.size();
638 if (nidx!=0){
639 nblock = nInitialOctants - nidx*nchm1;
640 nidx = 0;
641 for (idx=0; idx<nblock; idx++){
642 if (idx+offset < nInitialOctants){
643 if (nidx < nfchild){
644 if (idx+offset == first_child_index[nidx]){
645 markerfather = -m_treeConstants->maxLevel;
646 father = m_octants[idx+offset].buildFather();
647 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT; iii++){
648 father.m_info[iii] = false;
649 }
650 father.setGhostLayer(-1);
651 for(idx2=0; idx2<m_treeConstants->nChildren; idx2++){
652 if (idx2 < nInitialOctants){
653 if (markerfather < m_octants[idx+offset+idx2].getMarker()+1){
654 markerfather = m_octants[idx+offset+idx2].getMarker()+1;
655 }
656 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT; iii++){
657 father.m_info[iii] = father.m_info[iii] || m_octants[idx+offset+idx2].m_info[iii];
658 }
659 }
660 }
661 father.m_info[Octant::INFO_NEW4COARSENING] = true;
662 father.setMarker(markerfather);
663 //Impossible in this version
664// if (markerfather < 0 && mapsize == 0){
665// docoarse = true;
666// }
667 m_octants[idx] = father;
668 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
669 offset += nchm1;
670 nidx++;
671 }
672 else{
673 m_octants[idx] = m_octants[idx+offset];
674 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
675 }
676 }
677 else{
678 m_octants[idx] = m_octants[idx+offset];
679 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
680 }
681 }
682 }
683 }
684 m_octants.resize(nblock, Octant(m_dim));
685 m_octants.shrink_to_fit();
686 nInitialOctants = m_octants.size();
687 if(mapsize > 0){
688 mapidx.resize(nInitialOctants);
689 }
690
691 //Check ghosts
692 if (m_ghosts.size()){
693 // Start on ghosts
694 if (nInitialOctants > 0 && idx1_gh < nInitialGhosts){
695 if (m_ghosts[idx1_gh].buildFather() == m_octants[0].buildFather()){
696 father = m_ghosts[idx1_gh].buildFather();
697 nbro = 0;
698 idx = idx1_gh;
699 marker = m_ghosts[idx].getMarker();
700 bool waszero = (idx == 0);
701 while(marker < 0 && m_ghosts[idx].buildFather() == father){
702 nbro++;
703 if(waszero){
704 break;
705 }
706 idx--;
707 if(idx == 0){
708 waszero = true;
709 }
710 marker = m_ghosts[idx].getMarker();
711 }
712 nstart = 0;
713 idx = 0;
714 marker = m_octants[idx].getMarker();
715 if (idx==nInitialOctants-1) wstop = true;
716 while(marker < 0 && m_octants[idx].buildFather() == father){
717 nbro++;
718 nstart++;
719 if (wstop){
720 break;
721 }
722 idx++;
723 if (idx==nInitialOctants-1){
724 wstop = true;
725 }
726 marker = m_octants[idx].getMarker();
727 }
728 if (nbro == m_treeConstants->nChildren){
729 offset = nstart;
730 }
731 else{
732 nstart = 0;
733 }
734 }
735 if (nstart != 0){
736 for (idx=0; idx<nInitialOctants; idx++){
737 if (idx+offset < nInitialOctants){
738 m_octants[idx] = m_octants[idx+offset];
739 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
740 }
741 }
742 m_octants.resize(nInitialOctants-offset, Octant(m_dim));
743 m_octants.shrink_to_fit();
744 nInitialOctants = m_octants.size();
745 if(mapsize > 0){
746 mapidx.resize(nInitialOctants);
747 }
748 }
749 }
750
751
752 //Verify family between more then two processes
753 if (nInitialOctants > 0 && idx2_gh < nInitialGhosts){
754
755 if (m_ghosts[idx2_gh].buildFather() == father){
756
757 if (m_ghosts[idx2_gh].buildFather() == m_octants[nInitialOctants-1].buildFather()){
758
759 uint64_t idx22_gh = idx2_gh;
760 marker = m_ghosts[idx22_gh].getMarker();
761 while(marker < 0 && m_ghosts[idx22_gh].buildFather() == father){
762 nbro++;
763 idx22_gh++;
764 if(idx22_gh == nInitialGhosts){
765 break;
766 }
767 marker = m_ghosts[idx22_gh].getMarker();
768 }
769
770 if (nbro == m_treeConstants->nChildren){
771 m_octants.clear();
772 nInitialOctants = 0;
773 if(mapsize > 0){
774 mapidx.clear();
775 }
776 }
777 }
778
779 }
780
781 }
782
783
784 // End on ghosts
785 if (nInitialOctants > 0 && idx2_gh < nInitialGhosts){
786 if (m_ghosts[idx2_gh].buildFather() == m_octants[nInitialOctants-1].buildFather()){
787 father = m_ghosts[idx2_gh].buildFather();
788 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT; iii++){
789 father.m_info[iii] = false;
790 }
791 father.setGhostLayer(-1);
792 markerfather = m_ghosts[idx2_gh].getMarker()+1;
793 nbro = 0;
794 idx = idx2_gh;
795 marker = m_ghosts[idx].getMarker();
796 while(marker < 0 && m_ghosts[idx].buildFather() == father){
797 nbro++;
798 if (markerfather < m_ghosts[idx].getMarker()+1){
799 markerfather = m_ghosts[idx].getMarker()+1;
800 }
801 for (uint32_t iii=0; iii<m_treeConstants->nFaces; iii++){
802 father.m_info[iii] = father.m_info[iii] || m_ghosts[idx].m_info[iii];
803 }
804 father.m_info[Octant::INFO_BALANCED] = father.m_info[Octant::INFO_BALANCED] || m_ghosts[idx].m_info[Octant::INFO_BALANCED];
805 idx++;
806 if(idx == nInitialGhosts){
807 break;
808 }
809 marker = m_ghosts[idx].getMarker();
810 }
811 nend = 0;
812 idx = nInitialOctants-1;
813 marker = m_octants[idx].getMarker();
814 if (idx==0) wstop = true;
815 while(marker < 0 && m_octants[idx].buildFather() == father){
816 nbro++;
817 nend++;
818 if (markerfather < m_octants[idx].getMarker()+1){
819 markerfather = m_octants[idx].getMarker()+1;
820 }
821 idx--;
822 if (wstop){
823 break;
824 }
825 if (idx==0){
826 wstop = true;
827 }
828 marker = m_octants[idx].getMarker();
829 }
830 if (nbro == m_treeConstants->nChildren){
831 offset = nend;
832 }
833 else{
834 nend = 0;
835 }
836 }
837 if (nend != 0){
838 for (idx=0; idx < nend; idx++){
839 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT - 1; iii++){
840 father.m_info[iii] = father.m_info[iii] || m_octants[nInitialOctants-idx-1].m_info[iii];
841 }
842 }
843 father.m_info[Octant::INFO_NEW4COARSENING] = true;
844 father.setGhostLayer(-1);
845 //Impossible in this version
846 // if (markerfather < 0 && mapsize == 0){
847 // docoarse = true;
848 // }
849 father.setMarker(markerfather);
850 m_octants.resize(nInitialOctants-offset, Octant(m_dim));
851 m_octants.push_back(father);
852 m_octants.shrink_to_fit();
853 nInitialOctants = m_octants.size();
854 if(mapsize > 0){
855 mapidx.resize(nInitialOctants);
856 }
857 }
858 }
859
860 }//end if ghosts size
861
862 // Update maximum depth
863 updateLocalMaxDepth();
864
865 // Set final last desc
866 setFirstDescMorton();
867 setLastDescMorton();
868
869 return docoarse;
870
871 };
872
873 // =================================================================================== //
874
879 bool
880 LocalTree::globalRefine(u32vector & mapidx){
881
882 for (Octant &octant : m_octants){
883 octant.setMarker(1);
884 }
885
886 return refine(mapidx);
887
888 };
889
890 // =================================================================================== //
891
896 bool
897 LocalTree::globalCoarse(u32vector & mapidx){
898
899 for (Octant &octant : m_octants){
900 octant.setMarker(-1);
901 }
902 for (Octant &octant : m_ghosts){
903 octant.setMarker(-1);
904 }
905
906 return coarse(mapidx);
907
908 };
909
915 void
916 LocalTree::checkCoarse(uint64_t partLastDesc, u32vector & mapidx){
917
918 uint32_t idx;
919 uint32_t mapsize = mapidx.size();
920 uint8_t toDelete = 0;
921
922 uint32_t nInitialOctants = getNumOctants();
923 if (nInitialOctants>0){
924
925 idx = 0;
926 if (m_octants[idx].getMorton() < partLastDesc){
927
928 Octant father0 = m_octants[idx].buildFather();
929 Octant father = father0;
930
931 while(father == father0 && idx < nInitialOctants){
932 toDelete++;
933 idx++;
934 if (idx<nInitialOctants) father = m_octants[idx].buildFather();
935 }
936
937 if (nInitialOctants>toDelete){
938 for(idx=0; idx<nInitialOctants-toDelete; idx++){
939 m_octants[idx] = m_octants[idx+toDelete];
940 if (mapsize>0) mapidx[idx] = mapidx[idx+toDelete];
941 }
942 m_octants.resize(nInitialOctants-toDelete, Octant(m_dim));
943 if (mapsize>0){
944 mapidx.resize(m_octants.size());
945 }
946 }
947 else{
948 m_octants.clear();
949 mapidx.clear();
950 }
951 setFirstDescMorton();
952 }
953
954 }
955 };
956
957 // =================================================================================== //
960 void
961 LocalTree::updateLocalMaxDepth(){
962
963 uint32_t noctants = getNumOctants();
964 uint32_t i;
965
966 m_localMaxDepth = 0;
967 for(i = 0; i < noctants; i++){
968 if(m_octants[i].getLevel() > m_localMaxDepth){
969 m_localMaxDepth = m_octants[i].getLevel();
970 }
971 }
972 };
973
974
986 void
987 LocalTree::findNeighbours(const Octant* oct, uint8_t iface, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
988
989 if (!append) {
990 isghost.clear();
991 neighbours.clear();
992 }
993
994 // Default if iface is nface<iface<0
995 if (iface >= m_treeConstants->nFaces){
996 return;
997 }
998
999 // If a face is a boundary, it can have neighbours only if periodic
1000 bool isperiodic = false;
1001 if (oct->getBound(iface)) {
1002 isperiodic = isPeriodic(oct, iface);
1003 if (!isperiodic) {
1004 return;
1005 }
1006 }
1007
1008 // Exit if is the root octant
1009 uint8_t level = oct->getLevel();
1010 if (level == 0){
1011 // if periodic face return itself
1012 if (isperiodic){
1013 neighbours.push_back(0);
1014 isghost.push_back(false);
1015 }
1016 return;
1017 }
1018
1019 // Initialize search
1020 uint32_t candidateIdx = 0;
1021 uint64_t candidateMorton = 0;
1022
1023 uint64_t neighArea = 0;
1024 uint64_t faceArea = oct->getLogicalArea();
1025
1026 uint32_t size = oct->getLogicalSize();
1027
1028 std::array<uint32_t, 3> coord = oct->getLogicalCoordinates();
1029 if (isperiodic) {
1030 std::array<int64_t, 3> periodicOffset = getPeriodicOffset(*oct, iface);
1031 coord[0] = static_cast<uint32_t>(coord[0] + periodicOffset[0]);
1032 coord[1] = static_cast<uint32_t>(coord[1] + periodicOffset[1]);
1033 coord[2] = static_cast<uint32_t>(coord[2] + periodicOffset[2]);
1034 }
1035
1036 const int8_t (&cxyz)[3] = m_treeConstants->normals[iface];
1037
1038 // Compute same-size virtual neighbour information
1039 std::array<int64_t, 3> sameSizeVirtualNeighOffset = computeFirstVirtualNeighOffset(level, iface, level);
1040 std::array<uint32_t, 3> sameSizeVirtualNeighCoord = coord;
1041 sameSizeVirtualNeighCoord[0] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[0] + sameSizeVirtualNeighOffset[0]);
1042 sameSizeVirtualNeighCoord[1] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[1] + sameSizeVirtualNeighOffset[1]);
1043 sameSizeVirtualNeighCoord[2] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[2] + sameSizeVirtualNeighOffset[2]);
1044 uint64_t sameSizeVirtualNeighMorton = PABLO::computeMorton(m_dim, sameSizeVirtualNeighCoord[0], sameSizeVirtualNeighCoord[1], sameSizeVirtualNeighCoord[2]);
1045
1046 //
1047 // Search in the internal octants
1048 //
1049
1050 // Identify the index of the first neighbour candidate
1051 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_octants, &candidateIdx, &candidateMorton);
1052
1053 // Early return if a neighbour of the same size has been found
1054 if(candidateMorton == sameSizeVirtualNeighMorton && m_octants[candidateIdx].m_level == level){
1055 neighbours.push_back(candidateIdx);
1056 isghost.push_back(false);
1057 return;
1058 }
1059
1060 // Compute the Morton number of the last candidate
1061 //
1062 // This is the Morton number of the last descendant of the same-size
1063 // virtual neighbour.
1064 uint8_t maxNeighLevel = getMaxNeighLevel(*oct);
1065 std::array<int64_t, 3> lastCandidateOffset = computeLastVirtualNeighOffset(level, iface, maxNeighLevel);
1066 std::array<uint32_t, 3> lastCandidateCoord = coord;
1067 lastCandidateCoord[0] = static_cast<uint32_t>(lastCandidateCoord[0] + lastCandidateOffset[0]);
1068 lastCandidateCoord[1] = static_cast<uint32_t>(lastCandidateCoord[1] + lastCandidateOffset[1]);
1069 lastCandidateCoord[2] = static_cast<uint32_t>(lastCandidateCoord[2] + lastCandidateOffset[2]);
1070 uint64_t lastCandidateMorton = PABLO::computeMorton(m_dim, lastCandidateCoord[0], lastCandidateCoord[1], lastCandidateCoord[2]);
1071
1072 // Search for neighbours of different sizes
1073 if (candidateIdx < getNumOctants()){
1074 while(true){
1075 // Detect if the candidate is a neighbour
1076 u32array3 coordtry = {{0, 0, 0}};
1077 bool isNeighbourCandidate = true;
1078 for (int8_t idim=0; idim<m_dim; idim++){
1079 coordtry[idim] = m_octants[candidateIdx].getLogicalCoordinates(idim);
1080
1081 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(coordtry[idim] - coord[idim]));
1082 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_octants[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1083 if (Dx != Dxstar){
1084 isNeighbourCandidate = false;
1085 break;
1086 }
1087 }
1088
1089 if (isNeighbourCandidate){
1090 bool isNeighbour = false;
1091 uint8_t leveltry = m_octants[candidateIdx].getLevel();
1092 if (leveltry > level){
1093 array<int64_t,3> coord1 ={{1, 1, 1}} ;
1094 for (int8_t idim=0; idim<m_dim; idim++){
1095 coord1[idim] = coord[idim] + size;
1096 }
1097
1098 if((abs(cxyz[0])*((coordtry[1]>=coord[1])*(coordtry[1]<coord1[1]))*((coordtry[2]>=coord[2])*(coordtry[2]<coord1[2]))) + (abs(cxyz[1])*((coordtry[0]>=coord[0])*(coordtry[0]<coord1[0]))*((coordtry[2]>=coord[2])*(coordtry[2]<coord1[2]))) + (abs(cxyz[2])*((coordtry[0]>=coord[0])*(coordtry[0]<coord1[0]))*((coordtry[1]>=coord[1])*(coordtry[1]<coord1[1])))){
1099 isNeighbour = true;
1100 }
1101 }
1102 else if (leveltry < level){
1103 u32array3 coordtry1 = {{1, 1, 1}};
1104 for (int8_t idim=0; idim<m_dim; idim++){
1105 coordtry1[idim] = coordtry[idim] + m_octants[candidateIdx].getLogicalSize();
1106 }
1107
1108 if((abs(cxyz[0])*((coord[1]>=coordtry[1])*(coord[1]<coordtry1[1]))*((coord[2]>=coordtry[2])*(coord[2]<coordtry1[2]))) + (abs(cxyz[1])*((coord[0]>=coordtry[0])*(coord[0]<coordtry1[0]))*((coord[2]>=coordtry[2])*(coord[2]<coordtry1[2]))) + (abs(cxyz[2])*((coord[0]>=coordtry[0])*(coord[0]<coordtry1[0]))*((coord[1]>=coordtry[1])*(coord[1]<coordtry1[1])))){
1109 isNeighbour = true;
1110 }
1111 }
1112
1113 if (isNeighbour){
1114 neighbours.push_back(candidateIdx);
1115 isghost.push_back(false);
1116
1117 // If the neighbours already cover the whole face, we have
1118 // found all the neighbours and we can exit.
1119 neighArea += m_octants[candidateIdx].getLogicalArea();
1120 if (neighArea == faceArea){
1121 return;
1122 }
1123 }
1124 }
1125
1126 // Advance to the next candidate
1127 candidateIdx++;
1128 if (candidateIdx > getNumOctants() - 1){
1129 break;
1130 }
1131
1132 candidateMorton = m_octants[candidateIdx].getMorton();
1133 if (candidateMorton > lastCandidateMorton){
1134 break;
1135 }
1136 }
1137 }
1138
1139 //
1140 // Search in the ghost octants
1141 //
1142 // If we want also the neighbours that are ghosts, we always need to
1143 // search in the ghosts, the only exception is for faces of internal
1144 // octants that are not process boundaries.
1145 bool ghostSearch = !onlyinternal && (getNumGhosts() > 0);
1146 if (ghostSearch){
1147 if (!oct->getIsGhost() && !oct->getPbound(iface)){
1148 ghostSearch = false;
1149 }
1150 }
1151
1152 // Search in ghosts
1153 if(ghostSearch){
1154 // Identify the index of the first neighbour candidate
1155 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_ghosts, &candidateIdx, &candidateMorton);
1156
1157 // Early return if a neighbour of the same size has been found
1158 if(candidateMorton == sameSizeVirtualNeighMorton && m_ghosts[candidateIdx].getLevel() == level){
1159 neighbours.push_back(candidateIdx);
1160 isghost.push_back(true);
1161 return;
1162 }
1163
1164 // Search for neighbours of different sizes
1165 if (candidateIdx < getNumGhosts()){
1166 while(true){
1167 // Detect if the candidate is a neighbour
1168 u32array3 coordtry = {{0, 0, 0}};
1169 bool isNeighbourCandidate = true;
1170 for (int8_t idim=0; idim<m_dim; idim++){
1171 coordtry[idim] = m_ghosts[candidateIdx].getLogicalCoordinates(idim);
1172
1173 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(coordtry[idim] - coord[idim]));
1174 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_ghosts[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1175 if (Dx != Dxstar){
1176 isNeighbourCandidate = false;
1177 break;
1178 }
1179 }
1180
1181 if (isNeighbourCandidate){
1182 bool isNeighbour = false;
1183 uint8_t leveltry = m_ghosts[candidateIdx].getLevel();
1184 if (leveltry > level){
1185 array<int64_t, 3> coord1 = {{1, 1, 1}};
1186 for (int8_t idim=0; idim<m_dim; idim++){
1187 coord1[idim] = coord[idim] + size;
1188 }
1189
1190 if((abs(cxyz[0])*((coordtry[1]>=coord[1])*(coordtry[1]<coord1[1]))*((coordtry[2]>=coord[2])*(coordtry[2]<coord1[2]))) + (abs(cxyz[1])*((coordtry[0]>=coord[0])*(coordtry[0]<coord1[0]))*((coordtry[2]>=coord[2])*(coordtry[2]<coord1[2]))) + (abs(cxyz[2])*((coordtry[0]>=coord[0])*(coordtry[0]<coord1[0]))*((coordtry[1]>=coord[1])*(coordtry[1]<coord1[1])))){
1191 isNeighbour = true;
1192 }
1193 }
1194 else if (leveltry < level){
1195 u32array3 coordtry1 = {{1, 1, 1}};
1196 for (int8_t idim=0; idim<m_dim; idim++){
1197 coordtry1[idim] = coordtry[idim] + m_ghosts[candidateIdx].getLogicalSize();
1198 }
1199
1200 if((abs(cxyz[0])*((coord[1]>=coordtry[1])*(coord[1]<coordtry1[1]))*((coord[2]>=coordtry[2])*(coord[2]<coordtry1[2]))) + (abs(cxyz[1])*((coord[0]>=coordtry[0])*(coord[0]<coordtry1[0]))*((coord[2]>=coordtry[2])*(coord[2]<coordtry1[2]))) + (abs(cxyz[2])*((coord[0]>=coordtry[0])*(coord[0]<coordtry1[0]))*((coord[1]>=coordtry[1])*(coord[1]<coordtry1[1])))){
1201 isNeighbour = true;
1202 }
1203 }
1204
1205 if (isNeighbour){
1206 neighbours.push_back(candidateIdx);
1207 isghost.push_back(true);
1208
1209 // If the neighbours already cover the whole face, we have
1210 // found all the neighbours and we can exit.
1211 neighArea += m_ghosts[candidateIdx].getLogicalArea();
1212 if (neighArea == faceArea){
1213 return;
1214 }
1215 }
1216 }
1217
1218 candidateIdx++;
1219 if (candidateIdx > getNumGhosts() - 1){
1220 break;
1221 }
1222
1223 candidateMorton = m_ghosts[candidateIdx].getMorton();
1224 if (candidateMorton > lastCandidateMorton){
1225 break;
1226 }
1227 }
1228 }
1229 }
1230 };
1231
1243 void
1244 LocalTree::findEdgeNeighbours(const Octant* oct, uint8_t iedge, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
1245
1246 if (!append) {
1247 isghost.clear();
1248 neighbours.clear();
1249 }
1250
1251 // Default if iedge is nface<iedge<0
1252 if (iedge >= m_treeConstants->nEdges){
1253 return;
1254 }
1255
1256 // If a edge is a on boundary, it can have neighbours only if periodic
1257 bool isperiodic = false;
1258 if (oct->getEdgeBound(iedge)) {
1259 isperiodic = isEdgePeriodic(oct, iedge);
1260 if (!isperiodic) {
1261 return;
1262 }
1263 }
1264
1265 // Exit if is the root octant
1266 uint8_t level = oct->getLevel();
1267 if (level == 0){
1268 // if periodic face return itself
1269 if (isperiodic){
1270 neighbours.push_back(0);
1271 isghost.push_back(false);
1272 }
1273 return;
1274 }
1275
1276 // Search in the octants
1277 uint32_t candidateIdx = 0;
1278 uint64_t candidateMorton = 0;
1279
1280 uint32_t neighSize = 0;
1281 uint32_t edgeSize = oct->getLogicalSize();
1282
1283 uint32_t size = oct->getLogicalSize();
1284
1285 std::array<uint32_t, 3> coord = oct->getLogicalCoordinates();
1286 if (isperiodic) {
1287 std::array<int64_t, 3> periodicOffset = getEdgePeriodicOffset(*oct, iedge);
1288 coord[0] = static_cast<uint32_t>(coord[0] + periodicOffset[0]);
1289 coord[1] = static_cast<uint32_t>(coord[1] + periodicOffset[1]);
1290 coord[2] = static_cast<uint32_t>(coord[2] + periodicOffset[2]);
1291 }
1292
1293 const int8_t (&cxyz)[3] = m_treeConstants->edgeCoeffs[iedge];
1294
1295 // Compute same-size virtual neighbour information
1296 std::array<int64_t, 3> sameSizeVirtualNeighOffset = computeFirstVirtualEdgeNeighOffset(level, iedge, level);
1297 std::array<uint32_t, 3> sameSizeVirtualNeighCoord = coord;
1298 sameSizeVirtualNeighCoord[0] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[0] + sameSizeVirtualNeighOffset[0]);
1299 sameSizeVirtualNeighCoord[1] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[1] + sameSizeVirtualNeighOffset[1]);
1300 sameSizeVirtualNeighCoord[2] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[2] + sameSizeVirtualNeighOffset[2]);
1301 uint64_t sameSizeVirtualNeighMorton = PABLO::computeMorton(m_dim, sameSizeVirtualNeighCoord[0], sameSizeVirtualNeighCoord[1], sameSizeVirtualNeighCoord[2]);
1302
1303 //
1304 // Search in the internal octants
1305 //
1306
1307 // Identify the index of the first neighbour candidate
1308 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_octants, &candidateIdx, &candidateMorton);
1309
1310 // Early return if a neighbour of the same size has been found
1311 if(candidateMorton == sameSizeVirtualNeighMorton && m_octants[candidateIdx].m_level == level){
1312 isghost.push_back(false);
1313 neighbours.push_back(candidateIdx);
1314 return;
1315 }
1316
1317 // Compute the Morton number of the last candidate
1318 //
1319 // This is the Morton number of the last descendant of the same-size
1320 // virtual neighbour.
1321 uint8_t maxEdgeNeighLevel = getMaxEdgeNeighLevel(*oct);
1322 std::array<int64_t, 3> lastCandidateOffset = computeLastVirtualEdgeNeighOffset(level, iedge, maxEdgeNeighLevel);
1323 std::array<uint32_t, 3> lastCandidateCoord = coord;
1324 lastCandidateCoord[0] = static_cast<uint32_t>(lastCandidateCoord[0] + lastCandidateOffset[0]);
1325 lastCandidateCoord[1] = static_cast<uint32_t>(lastCandidateCoord[1] + lastCandidateOffset[1]);
1326 lastCandidateCoord[2] = static_cast<uint32_t>(lastCandidateCoord[2] + lastCandidateOffset[2]);
1327 uint64_t lastCandidateMorton = PABLO::computeMorton(m_dim, lastCandidateCoord[0], lastCandidateCoord[1], lastCandidateCoord[2]);
1328
1329 // Search for neighbours of different sizes
1330 if (candidateIdx < getNumOctants()) {
1331 while(true){
1332 // Detect if the candidate is a neighbour
1333 u32array3 coordtry = {{0, 0, 0}};
1334 bool isNeighbourCandidate = true;
1335 for (int8_t idim=0; idim<m_dim; idim++){
1336 coordtry[idim] = m_octants[candidateIdx].getLogicalCoordinates(idim);
1337
1338 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1339 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_octants[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1340
1341 if (Dx != Dxstar) {
1342 isNeighbourCandidate = false;
1343 break;
1344 }
1345 }
1346
1347 if (isNeighbourCandidate){
1348 bool isNeighbour = false;
1349 uint8_t leveltry = m_octants[candidateIdx].getLevel();
1350 if (leveltry > level){
1351 u32array3 coord1 = {{1, 1, 1}};
1352 for (int8_t idim=0; idim<m_dim; idim++){
1353 coord1[idim] = coord[idim] + size;
1354 }
1355
1356 if((abs(cxyz[0])*abs(cxyz[2])*((coordtry[1]>=coord[1])*(coordtry[1]<coord1[1]))) + (abs(cxyz[1])*abs(cxyz[2])*((coordtry[0]>=coord[0])*(coordtry[0]<coord1[0]))) + (abs(cxyz[0])*abs(cxyz[1])*((coordtry[2]>=coord[2])*(coordtry[2]<coord1[2])))){
1357 isNeighbour = true;
1358 }
1359 }
1360 else if (leveltry < level){
1361 u32array3 coordtry1 = {{1, 1, 1}};
1362 for (int8_t idim=0; idim<m_dim; idim++){
1363 coordtry1[idim] = coordtry[idim] + m_octants[candidateIdx].getLogicalSize();
1364 }
1365
1366 if((abs(cxyz[0])*abs(cxyz[2])*((coord[1]>=coordtry[1])*(coord[1]<coordtry1[1]))) + (abs(cxyz[1])*abs(cxyz[2])*((coord[0]>=coordtry[0])*(coord[0]<coordtry1[0]))) + (abs(cxyz[0])*abs(cxyz[1])*((coord[2]>=coordtry[2])*(coord[2]<coordtry1[2])))){
1367 isNeighbour = true;
1368 }
1369 }
1370
1371 if (isNeighbour) {
1372 neighbours.push_back(candidateIdx);
1373 isghost.push_back(false);
1374
1375 // If the neighbours already cover the whole edge, we have
1376 // found all the neighbours and we can exit.
1377 neighSize += m_octants[candidateIdx].getLogicalSize();
1378 if (neighSize == edgeSize){
1379 return;
1380 }
1381 }
1382 }
1383
1384 // Advance to the next candidate
1385 candidateIdx++;
1386 if (candidateIdx > getNumOctants()-1){
1387 break;
1388 }
1389
1390 candidateMorton = m_octants[candidateIdx].getMorton();
1391 if (candidateMorton > lastCandidateMorton){
1392 break;
1393 }
1394 }
1395 }
1396
1397 //
1398 // Search in the ghost octants
1399 //
1400 if (getNumGhosts() > 0 && !onlyinternal){
1401 // Identify the index of the first neighbour candidate
1402 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_ghosts, &candidateIdx, &candidateMorton);
1403
1404 // Early return if a neighbour of the same size has been found
1405 if(candidateMorton == sameSizeVirtualNeighMorton && m_ghosts[candidateIdx].m_level == level){
1406 isghost.push_back(true);
1407 neighbours.push_back(candidateIdx);
1408 return;
1409 }
1410
1411 // Search for neighbours of different sizes
1412 if (candidateIdx < getNumGhosts()){
1413 while(true){
1414 // Detect if the candidate is a neighbour
1415 u32array3 coordtry = {{0, 0, 0}};
1416 bool isNeighbourCandidate = true;
1417 for (int8_t idim=0; idim<m_dim; idim++){
1418 coordtry[idim] = m_ghosts[candidateIdx].getLogicalCoordinates(idim);
1419
1420 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1421 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_ghosts[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1422
1423 if (Dx != Dxstar) {
1424 isNeighbourCandidate = false;
1425 break;
1426 }
1427 }
1428
1429 if (isNeighbourCandidate){
1430 bool isNeighbour = false;
1431 uint8_t leveltry = m_ghosts[candidateIdx].getLevel();
1432 if (leveltry > level){
1433 u32array3 coord1 = {{1, 1, 1}};
1434 for (int8_t idim=0; idim<m_dim; idim++){
1435 coord1[idim] = int32_t(coord[idim]) + size;
1436 }
1437
1438 if((abs(cxyz[0])*abs(cxyz[2])*((coordtry[1]>=coord[1])*(coordtry[1]<coord1[1]))) + (abs(cxyz[1])*abs(cxyz[2])*((coordtry[0]>=coord[0])*(coordtry[0]<coord1[0]))) + (abs(cxyz[0])*abs(cxyz[1])*((coordtry[2]>=coord[2])*(coordtry[2]<coord1[2])))){
1439 isNeighbour = true;
1440 }
1441 }
1442 else if (leveltry < level){
1443 u32array3 coordtry1 = {{1, 1, 1}};
1444 for (int8_t idim=0; idim<m_dim; idim++){
1445 coordtry1[idim] = coordtry[idim] + m_ghosts[candidateIdx].getLogicalSize();
1446 }
1447
1448 if((abs(cxyz[0])*abs(cxyz[2])*((coord[1]>=coordtry[1])*(coord[1]<coordtry1[1]))) + (abs(cxyz[1])*abs(cxyz[2])*((coord[0]>=coordtry[0])*(coord[0]<coordtry1[0]))) + (abs(cxyz[0])*abs(cxyz[1])*((coord[2]>=coordtry[2])*(coord[2]<coordtry1[2])))){
1449 isNeighbour = true;
1450 }
1451 }
1452
1453 if (isNeighbour) {
1454 neighbours.push_back(candidateIdx);
1455 isghost.push_back(true);
1456
1457 // If the neighbours already cover the whole edge, we have
1458 // found all the neighbours and we can exit.
1459 neighSize += m_ghosts[candidateIdx].getLogicalSize();
1460 if (neighSize == edgeSize){
1461 return;
1462 }
1463 }
1464 }
1465
1466 // Advance to the next candidate
1467 candidateIdx++;
1468 if (candidateIdx > getNumGhosts()-1){
1469 break;
1470 }
1471
1472 candidateMorton = m_ghosts[candidateIdx].getMorton();
1473 if (candidateMorton > lastCandidateMorton){
1474 break;
1475 }
1476 }
1477 }
1478 }
1479 };
1480
1492 void
1493 LocalTree::findNodeNeighbours(const Octant* oct, uint8_t inode, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
1494
1495 if (!append) {
1496 isghost.clear();
1497 neighbours.clear();
1498 }
1499
1500 // Default if inode is nnodes<inode<0
1501 if (inode >= m_treeConstants->nNodes){
1502 return;
1503 }
1504
1505 // If a node is a on boundary, it can have neighbours only if periodic
1506 bool isperiodic = false;
1507 if (oct->getNodeBound(inode)) {
1508 isperiodic = isNodePeriodic(oct, inode);
1509 if (!isperiodic) {
1510 return;
1511 }
1512 }
1513
1514 // Exit if is the root octant
1515 uint8_t level = oct->getLevel();
1516 if (level == 0){
1517 // if periodic face return itself
1518 if (isperiodic){
1519 neighbours.push_back(0);
1520 isghost.push_back(false);
1521 }
1522 return;
1523 }
1524
1525 // Initialize search
1526 uint32_t candidateIdx = 0;
1527 uint64_t candidateMorton = 0;
1528
1529 uint32_t size = oct->getLogicalSize();
1530
1531 std::array<uint32_t, 3> coord = oct->getLogicalCoordinates();
1532 if (isperiodic) {
1533 std::array<int64_t, 3> periodicOffset = getNodePeriodicOffset(*oct, inode);
1534 coord[0] = static_cast<uint32_t>(coord[0] + periodicOffset[0]);
1535 coord[1] = static_cast<uint32_t>(coord[1] + periodicOffset[1]);
1536 coord[2] = static_cast<uint32_t>(coord[2] + periodicOffset[2]);
1537 }
1538
1539 const int8_t (&cxyz)[3] = m_treeConstants->nodeCoeffs[inode];
1540
1541 // Compute same-size virtual neighbour information
1542 std::array<int64_t, 3> sameSizeVirtualNeighOffset = computeFirstVirtualNodeNeighOffset(level, inode, level);
1543 std::array<uint32_t, 3> sameSizeVirtualNeighCoord = coord;
1544 sameSizeVirtualNeighCoord[0] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[0] + sameSizeVirtualNeighOffset[0]);
1545 sameSizeVirtualNeighCoord[1] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[1] + sameSizeVirtualNeighOffset[1]);
1546 sameSizeVirtualNeighCoord[2] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[2] + sameSizeVirtualNeighOffset[2]);
1547 uint64_t sameSizeVirtualNeighMorton = PABLO::computeMorton(m_dim, sameSizeVirtualNeighCoord[0], sameSizeVirtualNeighCoord[1], sameSizeVirtualNeighCoord[2]);
1548
1549 //
1550 // Search in the internal octants
1551 //
1552
1553 // Identify the index of the first neighbour candidate
1554 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_octants, &candidateIdx, &candidateMorton);
1555
1556 // Early return if a neighbour of the same size has been found
1557 if(candidateMorton == sameSizeVirtualNeighMorton && m_octants[candidateIdx].m_level == oct->m_level){
1558 isghost.push_back(false);
1559 neighbours.push_back(candidateIdx);
1560 return;
1561 }
1562
1563 // Compute the Morton number of the last candidate
1564 //
1565 // This is the Morton number of the last descendant of the same-size
1566 // virtual neighbour.
1567 uint8_t maxNodeNeighLevel = getMaxNodeNeighLevel(*oct);
1568 std::array<int64_t, 3> lastCandidateOffset = computeLastVirtualNodeNeighOffset(level, inode, maxNodeNeighLevel);
1569 std::array<uint32_t, 3> lastCandidateCoord = coord;
1570 lastCandidateCoord[0] = static_cast<uint32_t>(lastCandidateCoord[0] + lastCandidateOffset[0]);
1571 lastCandidateCoord[1] = static_cast<uint32_t>(lastCandidateCoord[1] + lastCandidateOffset[1]);
1572 lastCandidateCoord[2] = static_cast<uint32_t>(lastCandidateCoord[2] + lastCandidateOffset[2]);
1573 uint64_t lastCandidateMorton = PABLO::computeMorton(m_dim, lastCandidateCoord[0], lastCandidateCoord[1], lastCandidateCoord[2]);
1574
1575 // Search for neighbours of different sizes
1576 if (candidateIdx < getNumOctants()) {
1577 while(true){
1578 // Detect if the candidate is a neighbour
1579 u32array3 coordtry = {{0, 0, 0}};
1580 bool isNeighbour = true;
1581 for (int8_t idim=0; idim<m_dim; idim++){
1582 coordtry[idim] = m_octants[candidateIdx].getLogicalCoordinates(idim);
1583
1584 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1585 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_octants[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1586
1587 if (Dx != Dxstar) {
1588 isNeighbour = false;
1589 break;
1590 }
1591 }
1592
1593 // Advance to the next candidate
1594 if (isNeighbour){
1595 neighbours.push_back(candidateIdx);
1596 isghost.push_back(false);
1597 return;
1598 }
1599
1600 candidateIdx++;
1601 if (candidateIdx > getNumOctants()-1){
1602 break;
1603 }
1604
1605 candidateMorton = m_octants[candidateIdx].getMorton();
1606 if (candidateMorton > lastCandidateMorton){
1607 break;
1608 }
1609 }
1610 }
1611
1612 //
1613 // Search in the ghost octants
1614 //
1615
1616 if (getNumGhosts() > 0 && !onlyinternal){
1617 // Identify the index of the first neighbour candidate
1618 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_ghosts, &candidateIdx, &candidateMorton);
1619
1620 // Early return if a neighbour of the same size has been found
1621 if(candidateMorton == sameSizeVirtualNeighMorton && m_ghosts[candidateIdx].m_level == oct->m_level){
1622 isghost.push_back(true);
1623 neighbours.push_back(candidateIdx);
1624 return;
1625 }
1626
1627 // Search for neighbours of different sizes
1628 if (candidateIdx < getNumGhosts()) {
1629 while(true){
1630 // Detect if the candidate is a neighbour
1631 u32array3 coordtry = {{0, 0, 0}};
1632 bool isNeighbour = true;
1633 for (int8_t idim=0; idim<m_dim; idim++){
1634 coordtry[idim] = m_ghosts[candidateIdx].getLogicalCoordinates(idim);
1635
1636 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1637 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_ghosts[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1638
1639 if (Dx != Dxstar) {
1640 isNeighbour = false;
1641 break;
1642 }
1643 }
1644
1645 // Advance to the next candidate
1646 if (isNeighbour){
1647 neighbours.push_back(candidateIdx);
1648 isghost.push_back(true);
1649 return;
1650 }
1651
1652 candidateIdx++;
1653 if (candidateIdx > getNumGhosts()-1){
1654 break;
1655 }
1656
1657 candidateMorton = m_ghosts[candidateIdx].getMorton();
1658 if (candidateMorton > lastCandidateMorton){
1659 break;
1660 }
1661 }
1662 }
1663 }
1664 };
1665
1666 // =================================================================================== //
1683 void
1684 LocalTree::computeNeighSearchBegin(uint64_t sameSizeVirtualNeighMorton, const octvector &octants, uint32_t *searchBeginIdx, uint64_t *searchBeginMorton) const {
1685
1686 // Early return if there are no octants
1687 if (octants.empty()) {
1688 *searchBeginIdx = 0;
1689 *searchBeginMorton = PABLO::INVALID_MORTON;
1690 return;
1691 }
1692
1693 // The search should start from the lower bound if it points to the
1694 // first octant or to an octant whose Morton number is equal to the
1695 // the same-size virtual neighbour Morton number. Otherwise, the
1696 // search should start form the octant preceding the lower bound.
1697 uint32_t lowerBoundIdx;
1698 uint64_t lowerBoundMorton;
1699 findMortonLowerBound(sameSizeVirtualNeighMorton, octants, &lowerBoundIdx, &lowerBoundMorton);
1700
1701 if (lowerBoundMorton == sameSizeVirtualNeighMorton || lowerBoundIdx == 0) {
1702 *searchBeginIdx = lowerBoundIdx;
1703 *searchBeginMorton = lowerBoundMorton;
1704 } else {
1705 *searchBeginIdx = lowerBoundIdx - 1;
1706 *searchBeginMorton = octants[*searchBeginIdx].getMorton();
1707 }
1708
1709 }
1710
1711 // =================================================================================== //
1712
1720 bool
1721 LocalTree::fixBrokenFamiliesMarkers(std::vector<Octant *> *updatedOctants, std::vector<bool> *updatedGhostFlags){
1722
1723 // Initialization
1724 bool updated = false;
1725
1726 uint32_t nocts = m_octants.size();
1727 if (nocts == 0) {
1728 return updated;
1729 }
1730
1731 uint32_t nghosts = m_ghosts.size();
1732
1733 uint32_t idx = 0;
1734 uint32_t idx0 = 0;
1735 uint32_t idx2 = 0;
1736 uint32_t idx1_gh = 0;
1737 uint32_t idx2_gh = 0;
1738 uint32_t last_idx =nocts-1;
1739
1740 uint8_t nbro = 0;
1741 Octant father(m_dim);
1742
1743 //Clean index of ghost brothers in case of coarsening a broken family
1744 m_lastGhostBros.clear();
1745 m_firstGhostBros.clear();
1746
1747 // Process ghosts
1748 bool checkend = true;
1749 bool checkstart = true;
1750 if (nghosts > 0){
1751 // Set index for start and end check for ghosts
1752 while(m_ghosts[idx2_gh].getMorton() <= m_lastDescMorton){
1753 idx2_gh++;
1754 if (idx2_gh > nghosts-1) break;
1755 }
1756 if (idx2_gh > nghosts-1) checkend = false;
1757
1758 while(m_ghosts[idx1_gh].getMorton() <= m_octants[0].getMorton()){
1759 idx1_gh++;
1760 if (idx1_gh > nghosts-1) break;
1761 }
1762 if (idx1_gh == 0) checkstart = false;
1763 idx1_gh-=1;
1764
1765 // Start and End on ghosts
1766 if (checkstart){
1767 if (m_ghosts[idx1_gh].buildFather()==m_octants[0].buildFather()){
1768 father = m_ghosts[idx1_gh].buildFather();
1769 nbro = 0;
1770 idx = idx1_gh;
1771 int8_t marker = m_ghosts[idx].getMarker();
1772 while(marker < 0 && m_ghosts[idx].buildFather() == father){
1773
1774 //Add ghost index to structure for mapper in case of coarsening a broken family
1775 m_firstGhostBros.push_back(idx);
1776
1777 nbro++;
1778 if (idx==0)
1779 break;
1780 idx--;
1781 marker = m_ghosts[idx].getMarker();
1782 }
1783 idx = 0;
1784 while(idx<nocts && m_octants[idx].buildFather() == father){
1785 if (m_octants[idx].getMarker()<0)
1786 nbro++;
1787 idx++;
1788 if(idx==nocts)
1789 break;
1790 }
1791 if (nbro != m_treeConstants->nChildren && idx!=nocts){
1792 for(uint32_t ii=0; ii<idx; ii++){
1793 Octant &octant = m_octants[ii];
1794 if (octant.getMarker()<0){
1795 octant.setMarker(0);
1796 if (updatedOctants) {
1797 updatedOctants->push_back(&octant);
1798 }
1799 if (updatedGhostFlags) {
1800 updatedGhostFlags->push_back(false);
1801 }
1802 updated = true;
1803 }
1804 }
1805 //Clean index of ghost brothers in case of coarsening a broken family
1806 m_firstGhostBros.clear();
1807 }
1808 }
1809 }
1810
1811 if (checkend){
1812 bool checklocal = false;
1813 if (m_ghosts[idx2_gh].buildFather()==m_octants[nocts-1].buildFather()){
1814 father = m_ghosts[idx2_gh].buildFather();
1815 if (idx!=nocts){
1816 nbro = 0;
1817 checklocal = true;
1818 }
1819 idx = idx2_gh;
1820 int8_t marker = m_ghosts[idx].getMarker();
1821 while(marker < 0 && m_ghosts[idx].buildFather() == father){
1822
1823 //Add ghost index to structure for mapper in case of coarsening a broken family
1824 m_lastGhostBros.push_back(idx);
1825
1826 nbro++;
1827 idx++;
1828 if(idx == nghosts){
1829 break;
1830 }
1831 marker = m_ghosts[idx].getMarker();
1832 }
1833 idx = nocts-1;
1834 if (checklocal){
1835 while(m_octants[idx].buildFather() == father){
1836 if (m_octants[idx].getMarker()<0)
1837 nbro++;
1838 if (idx==0)
1839 break;
1840 idx--;
1841 }
1842 }
1843 last_idx=idx;
1844 if (nbro != m_treeConstants->nChildren && idx!=nocts-1){
1845 for(uint32_t ii=idx+1; ii<nocts; ii++){
1846 Octant &octant = m_octants[ii];
1847 if (octant.getMarker()<0){
1848 octant.setMarker(0);
1849 if (updatedOctants) {
1850 updatedOctants->push_back(&octant);
1851 }
1852 if (updatedGhostFlags) {
1853 updatedGhostFlags->push_back(false);
1854 }
1855 updated = true;
1856 }
1857 }
1858 //Clean index of ghost brothers in case of coarsening a broken family
1859 m_lastGhostBros.clear();
1860 }
1861 }
1862 }
1863 }
1864
1865 // Check first internal octants
1866 father = m_octants[0].buildFather();
1867 uint64_t mortonld = father.computeLastDescMorton();
1868 nbro = 0;
1869 for (idx=0; idx<m_treeConstants->nChildren; idx++){
1870 // Check if family is complete or to be checked in the internal loop (some brother refined)
1871 if (idx<nocts){
1872 if (m_octants[idx].getMorton() <= mortonld){
1873 nbro++;
1874 }
1875 }
1876 }
1877
1878 if (nbro != m_treeConstants->nChildren) {
1879 idx0 = nbro;
1880 }
1881
1882 // Check and coarse internal octants
1883 for (idx=idx0; idx<nocts; idx++){
1884 Octant &octant = m_octants[idx];
1885 if(octant.getMarker() < 0 && octant.getLevel() > 0){
1886 nbro = 0;
1887 father = octant.buildFather();
1888 // Check if family is to be coarsened
1889 for (idx2=idx; idx2<idx+m_treeConstants->nChildren; idx2++){
1890 if (idx2<nocts){
1891 if(m_octants[idx2].getMarker() < 0 && m_octants[idx2].buildFather() == father){
1892 nbro++;
1893 }
1894 }
1895 }
1896 if (nbro == m_treeConstants->nChildren){
1897 idx = idx2-1;
1898 }
1899 else{
1900 if (idx<=last_idx){
1901 octant.setMarker(0);
1902 if (updatedOctants) {
1903 updatedOctants->push_back(&octant);
1904 }
1905 if (updatedGhostFlags) {
1906 updatedGhostFlags->push_back(false);
1907 }
1908 updated = true;
1909 }
1910 }
1911 }
1912 }
1913
1914 return updated;
1915 }
1916
1917 // =================================================================================== //
1918
1925 bool
1926 LocalTree::localBalance(bool doNew, bool checkInterior, bool checkGhost){
1927
1928 bool balanceEdges = ((m_balanceCodim>1) && (m_dim==3));
1929 bool balanceNodes = (m_balanceCodim==m_dim);
1930
1931 std::vector<Octant *> processOctants;
1932 std::vector<bool> processGhostFlags;
1933
1934 // Identify internal octants that will be processed
1935 if(checkInterior){
1936 for (Octant &octant : m_octants){
1937 // Skip octants that doesn't need balancing
1938 bool balanceOctant = octant.getBalance();
1939 if (balanceOctant) {
1940 if (doNew) {
1941 balanceOctant = (octant.getMarker() != 0) || octant.getIsNewC() || octant.getIsNewR();
1942 } else {
1943 balanceOctant = (octant.getMarker() != 0);
1944 }
1945 }
1946
1947 if (!balanceOctant) {
1948 continue;
1949 }
1950
1951 // Add octant to the process list
1952 processOctants.push_back(&octant);
1953 processGhostFlags.push_back(false);
1954 }
1955 }
1956
1957 // Identify ghost octants that will be processed
1958 //
1959 // Ghost octants will be balanced by the process that owns them, however they may
1960 // affect balancing of local octants. If ghost octnts are processed, we process all
1961 // ghost octants of the first layer, not only the ones that need balancing. That's
1962 // because it's faster to process all the ghost octants that may affect balacning
1963 // of internal octants, rather than find the ones that actually affect balacing of
1964 // internal octants.
1965 if (checkGhost) {
1966 for (Octant &octant : m_ghosts){
1967 // Only ghosts of the first layer can affect load balance
1968 if (octant.getGhostLayer() > 0) {
1969 continue;
1970 }
1971
1972 // Add octant to the process list
1973 processOctants.push_back(&octant);
1974 processGhostFlags.push_back(true);
1975 }
1976 }
1977
1978 // Iterative balacing
1979 std::vector<uint32_t> neighs;
1980 std::vector<bool> neighGhostFlags;
1981
1982 bool updated = false;
1983 std::vector<Octant *> updatedProcessOctants;
1984 std::vector<bool> updatedProcessGhostFlags;
1985 while (!processOctants.empty()) {
1986 // Fix broken families markers
1987 bool brokenFamilieisFixed = fixBrokenFamiliesMarkers(&processOctants, &processGhostFlags);
1988 if (brokenFamilieisFixed) {
1989 updated = true;
1990 }
1991
1992 // Balance
1993 std::size_t processSize = processOctants.size();
1994 for (std::size_t n = 0; n < processSize; ++n) {
1995 Octant &octant = *(processOctants[n]);
1996 bool ghostFlag = processGhostFlags[n];
1997
1998 // Identify neighbours that will be used for balancing
1999 neighs.clear();
2000 neighGhostFlags.clear();
2001
2002 for (int iface=0; iface < m_treeConstants->nFaces; iface++) {
2003 findNeighbours(&octant, iface, neighs, neighGhostFlags, ghostFlag, true);
2004 }
2005
2006 if (balanceNodes) {
2007 for (int inode=0; inode < m_treeConstants->nNodes; inode++) {
2008 findNodeNeighbours(&octant, inode, neighs, neighGhostFlags, ghostFlag, true);
2009 }
2010 }
2011
2012 if (balanceEdges) {
2013 for (int iedge=0; iedge < m_treeConstants->nEdges; iedge++) {
2014 findEdgeNeighbours(&octant, iedge, neighs, neighGhostFlags, ghostFlag, true);
2015 }
2016 }
2017
2018 // Balance octant
2019 int8_t level = octant.getLevel();
2020 int8_t futureLevel = std::min(m_treeConstants->maxLevel, static_cast<int8_t>(level + octant.getMarker()));
2021
2022 std::size_t nNeighs = neighs.size();
2023 for(std::size_t i = 0; i < nNeighs; i++){
2024 bool neighGhostFlag = neighGhostFlags[i];
2025 Octant &neighOctant = (neighGhostFlag ? m_ghosts[neighs[i]]: m_octants[neighs[i]]);
2026
2027 int8_t neighLevel = neighOctant.getLevel();
2028 int8_t neighFutureLevel = std::min(m_treeConstants->maxLevel, int8_t(neighLevel + neighOctant.getMarker()));
2029 if(!ghostFlag && futureLevel < neighFutureLevel - 1) {
2030 futureLevel = neighFutureLevel - 1;
2031 octant.setMarker(futureLevel - level);
2032 updatedProcessOctants.push_back(&octant);
2033 updatedProcessGhostFlags.push_back(ghostFlag);
2034 updated = true;
2035 }
2036 else if(!neighGhostFlag && neighFutureLevel < futureLevel - 1) {
2037 neighOctant.setMarker((futureLevel - 1) - neighLevel);
2038 if (neighOctant.getBalance()) {
2039 updatedProcessOctants.push_back(&neighOctant);
2040 updatedProcessGhostFlags.push_back(neighGhostFlag);
2041 }
2042 updated = true;
2043 }
2044 }
2045 }
2046
2047 // Update process list
2048 updatedProcessOctants.swap(processOctants);
2049 updatedProcessGhostFlags.swap(processGhostFlags);
2050
2051 updatedProcessOctants.clear();
2052 updatedProcessGhostFlags.clear();
2053 }
2054
2055 return updated;
2056 };
2057
2058 // =================================================================================== //
2059
2068 std::array<int64_t, 3> LocalTree::computeFirstVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const {
2069
2070 // Get normalized offsets
2071 const int8_t (&normalizedOffsets)[3] = m_treeConstants->normals[iface];
2072
2073 // Get octant sizes
2074 uint32_t size = m_treeConstants->lengths[level];
2075 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2076
2077 // Compute the coordinates of the virtual neighbour
2078 std::array<int64_t, 3> neighOffsets;
2079 for (int i = 0; i < 3; ++i) {
2080 neighOffsets[i] = normalizedOffsets[i];
2081 if (neighOffsets[i] > 0) {
2082 neighOffsets[i] *= size;
2083 } else {
2084 neighOffsets[i] *= neighSize;
2085 }
2086 }
2087
2088 return neighOffsets;
2089 }
2090
2099 std::array<int64_t, 3> LocalTree::computeLastVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const {
2100
2101 // Get octant sizes
2102 uint32_t size = m_treeConstants->lengths[level];
2103 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2104
2105 // Compute the coordinates of the virtual neighbour
2106 std::array<int64_t, 3> neighOffsets = computeFirstVirtualNeighOffset(level, iface, neighLevel);
2107
2108 uint32_t delta = size - neighSize;
2109 switch (iface) {
2110 case 0 :
2111 case 1 :
2112 neighOffsets[1] += delta;
2113 neighOffsets[2] += delta;
2114 break;
2115
2116 case 2 :
2117 case 3 :
2118 neighOffsets[0] += delta;
2119 neighOffsets[2] += delta;
2120 break;
2121
2122 case 4 :
2123 case 5 :
2124 neighOffsets[0] += delta;
2125 neighOffsets[1] += delta;
2126 break;
2127 }
2128
2129 return neighOffsets;
2130 }
2131
2140 void LocalTree::computeVirtualNeighOffsets(uint8_t level, uint8_t iface, uint8_t neighLevel, std::vector<std::array<int64_t, 3>> *neighOffsets) const {
2141
2142 // Get octant sizes
2143 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2144
2145 // Direction of the increments
2146 int direction1;
2147 int direction2;
2148 switch (iface) {
2149 case 0:
2150 case 1:
2151 direction1 = 1;
2152 direction2 = 2;
2153 break;
2154
2155 case 2:
2156 case 3:
2157 direction1 = 0;
2158 direction2 = 2;
2159 break;
2160
2161 default:
2162 direction1 = 0;
2163 direction2 = 1;
2164 break;
2165 }
2166
2167 // Compute the coordinates of the virtual neighbour
2168 int nNeighsDirection1 = 1 << (neighLevel - level);
2169 int nNeighsDirection2 = (m_dim > 2) ? nNeighsDirection1 : 1;
2170 int nNeighs = nNeighsDirection1 * nNeighsDirection2;
2171
2172 neighOffsets->assign(nNeighs, computeFirstVirtualNeighOffset(level, iface, neighLevel));
2173
2174 int k = 0;
2175 for (int j = 0; j < nNeighsDirection2; ++j) {
2176 for (int i = 0; i < nNeighsDirection1; ++i) {
2177 (*neighOffsets)[k][direction1] += i * neighSize;
2178 (*neighOffsets)[k][direction2] += j * neighSize;
2179 ++k;
2180 }
2181 }
2182 }
2183
2192 std::array<int64_t, 3> LocalTree::computeFirstVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const {
2193
2194 // Get normalized offsets
2195 const int8_t (&normalizedOffsets)[3] = m_treeConstants->nodeCoeffs[inode];
2196
2197 // Get octant sizes
2198 uint32_t size = m_treeConstants->lengths[level];
2199 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2200
2201 // Compute the coordinates of the virtual neighbour
2202 std::array<int64_t, 3> neighOffsets;
2203 for (int i = 0; i < 3; ++i) {
2204 neighOffsets[i] = normalizedOffsets[i];
2205 if (neighOffsets[i] > 0) {
2206 neighOffsets[i] *= size;
2207 } else {
2208 neighOffsets[i] *= neighSize;
2209 }
2210 }
2211
2212 return neighOffsets;
2213 }
2214
2223 std::array<int64_t, 3> LocalTree::computeLastVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const {
2224
2225 return computeFirstVirtualNodeNeighOffset(level, inode, neighLevel);
2226 }
2227
2236 void LocalTree::computeVirtualNodeNeighOffsets(uint8_t level, uint8_t inode, uint8_t neighLevel, std::vector<std::array<int64_t, 3>> *neighOffsets) const {
2237
2238 neighOffsets->assign(1, computeFirstVirtualNodeNeighOffset(level, inode, neighLevel));
2239 }
2240
2249 std::array<int64_t, 3> LocalTree::computeFirstVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const {
2250
2251 // Get normalized offsets
2252 const int8_t (&normalizedOffsets)[3] = m_treeConstants->edgeCoeffs[iedge];
2253
2254 // Get octant sizes
2255 uint32_t size = m_treeConstants->lengths[level];
2256 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2257
2258 // Compute the coordinates of the virtual neighbour
2259 std::array<int64_t, 3> neighOffsets;
2260 for (int i = 0; i < 3; ++i) {
2261 neighOffsets[i] = normalizedOffsets[i];
2262 if (neighOffsets[i] > 0) {
2263 neighOffsets[i] *= size;
2264 } else {
2265 neighOffsets[i] *= neighSize;
2266 }
2267 }
2268
2269 return neighOffsets;
2270 }
2271
2280 std::array<int64_t, 3> LocalTree::computeLastVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const {
2281
2282 // Get octant sizes
2283 uint32_t size = m_treeConstants->lengths[level];
2284 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2285
2286 // Compute the coordinates of the virtual neighbour
2287 std::array<int64_t, 3> neighOffsets = computeFirstVirtualEdgeNeighOffset(level, iedge, neighLevel);
2288
2289 uint32_t offset = size - neighSize;
2290 switch (iedge) {
2291 case 0:
2292 case 1:
2293 case 8:
2294 case 9:
2295 neighOffsets[1] += offset;
2296 break;
2297
2298 case 2:
2299 case 3:
2300 case 10:
2301 case 11:
2302 neighOffsets[0] += offset;
2303 break;
2304
2305 case 4:
2306 case 5:
2307 case 6:
2308 case 7:
2309 neighOffsets[2] += offset;
2310 break;
2311 }
2312
2313 return neighOffsets;
2314 }
2315
2324 void LocalTree::computeVirtualEdgeNeighOffsets(uint8_t level, uint8_t iedge, uint8_t neighLevel, std::vector<std::array<int64_t, 3>> *neighOffsets) const {
2325
2326 // Get octant sizes
2327 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2328
2329 // Direction of the increments
2330 int direction;
2331 switch (iedge) {
2332 case 0:
2333 case 1:
2334 case 8:
2335 case 9:
2336 direction = 1;
2337 break;
2338
2339 case 2:
2340 case 3:
2341 case 10:
2342 case 11:
2343 direction = 0;
2344 break;
2345
2346 default:
2347 direction = 2;
2348 break;
2349 }
2350
2351 // Compute the coordinates of the virtual neighbours
2352 int nNeighs = 1 << (neighLevel - level);
2353 neighOffsets->assign(nNeighs, computeFirstVirtualEdgeNeighOffset(level, iedge, neighLevel));
2354 for (int i = 1; i < nNeighs; ++i) {
2355 (*neighOffsets)[i][direction] += i * neighSize;
2356 }
2357 }
2358
2359 // =================================================================================== //
2360
2368 std::array<int64_t, 3> LocalTree::getPeriodicOffset(const Octant &octant, uint8_t iface) const {
2369
2370 std::array<int64_t, 3> offset = {{0, 0, 0}};
2371 if (!isPeriodic(&octant, iface)) {
2372 return offset;
2373 }
2374
2375 int64_t maxLength = m_treeConstants->lengths[0];
2376 const int8_t (&referenceOffset)[3] = m_treeConstants->normals[iface];
2377
2378 for (int d = 0; d < m_dim; ++d) {
2379 offset[d] = - maxLength * static_cast<int64_t>(referenceOffset[d]);
2380 }
2381
2382 return offset;
2383 }
2384
2392 std::array<int64_t, 3> LocalTree::getNodePeriodicOffset(const Octant &octant, uint8_t inode) const {
2393
2394 std::array<int64_t, 3> offset = {{0, 0, 0}};
2395 for (int i = 0; i < m_dim; ++i) {
2396 uint8_t face = m_treeConstants->nodeFace[inode][i];
2397 if (!isPeriodic(&octant, face)) {
2398 continue;
2399 }
2400
2401 std::array<int64_t, 3> faceOffset = getPeriodicOffset(octant, face);
2402 for (int d = 0; d < m_dim; ++d) {
2403 offset[d] += faceOffset[d];
2404 }
2405 }
2406
2407 return offset;
2408 }
2409
2417 std::array<int64_t, 3> LocalTree::getEdgePeriodicOffset(const Octant &octant, uint8_t iedge) const {
2418
2419 std::array<int64_t, 3> offset = {{0, 0, 0}};
2420 for (uint8_t face : m_treeConstants->edgeFace[iedge]) {
2421 if (!isPeriodic(&octant, face)) {
2422 continue;
2423 }
2424
2425 std::array<int64_t, 3> faceOffset = getPeriodicOffset(octant, face);
2426 for (int d = 0; d < m_dim; ++d) {
2427 offset[d] += faceOffset[d];
2428 }
2429 }
2430
2431 return offset;
2432 }
2433
2434 // =================================================================================== //
2435
2442 uint8_t LocalTree::getMaxNeighLevel(const Octant &octant) const {
2443
2444 if (octant.getBalance()) {
2445 return octant.getLevel() + 1;
2446 } else {
2447 return m_treeConstants->maxLevel;
2448 }
2449 }
2450
2457 uint8_t LocalTree::getMaxNodeNeighLevel(const Octant &octant) const {
2458
2459 if (octant.getBalance() && m_balanceCodim == m_dim) {
2460 return octant.getLevel() + 1;
2461 } else {
2462 return m_treeConstants->maxLevel;
2463 }
2464 }
2465
2472 uint8_t LocalTree::getMaxEdgeNeighLevel(const Octant &octant) const {
2473
2474 if (octant.getBalance() && m_balanceCodim > 1) {
2475 return octant.getLevel() + 1;
2476 } else {
2477 return m_treeConstants->maxLevel;
2478 }
2479 }
2480
2481 // =================================================================================== //
2482
2485 void
2486 LocalTree::computeIntersections() {
2487
2488 octvector::iterator it, obegin, oend;
2489 Intersection intersection;
2490 u32vector neighbours;
2491 vector<bool> isghost;
2492 uint32_t idx;
2493 uint32_t i, nsize;
2494 uint8_t iface, iface2;
2495
2496 m_intersections.clear();
2497 m_intersections.reserve(2*3*m_octants.size());
2498
2499 idx = 0;
2500
2501 // Loop on ghosts
2502 obegin = m_ghosts.begin();
2503 oend = m_ghosts.end();
2504 for (it = obegin; it != oend; ++it){
2505 for (iface = 0; iface < m_dim; iface++){
2506 iface2 = iface*2;
2507 findNeighbours(m_ghosts.data() + idx, iface2, neighbours, isghost, true, false);
2508 nsize = neighbours.size();
2509 if (!(it->m_info[iface2])){
2510 //Internal intersection
2511 for (i = 0; i < nsize; i++){
2512 intersection.m_dim = m_dim;
2513 intersection.m_finer = getGhostLevel(idx) >= getLevel((int)neighbours[i]);
2514 intersection.m_out = intersection.m_finer;
2515 intersection.m_outisghost = intersection.m_finer;
2516 intersection.m_owners[0] = neighbours[i];
2517 intersection.m_owners[1] = idx;
2518 intersection.m_iface = m_treeConstants->oppositeFace[iface2] - (getGhostLevel(idx) >= getLevel((int)neighbours[i]));
2519 intersection.m_isnew = false;
2520 intersection.m_isghost = true;
2521 intersection.m_bound = false;
2522 intersection.m_pbound = true;
2523 m_intersections.push_back(intersection);
2524 }
2525 }
2526 else{
2527 //Periodic intersection
2528 for (i = 0; i < nsize; i++){
2529 intersection.m_dim = m_dim;
2530 intersection.m_finer = getGhostLevel(idx) >= getLevel((int)neighbours[i]);
2531 intersection.m_out = intersection.m_finer;
2532 intersection.m_outisghost = intersection.m_finer;
2533 intersection.m_owners[0] = neighbours[i];
2534 intersection.m_owners[1] = idx;
2535 intersection.m_iface = m_treeConstants->oppositeFace[iface2] - (getGhostLevel(idx) >= getLevel((int)neighbours[i]));
2536 intersection.m_isnew = false;
2537 intersection.m_isghost = true;
2538 intersection.m_bound = true;
2539 intersection.m_pbound = true;
2540 m_intersections.push_back(intersection);
2541 }
2542 }
2543 }
2544 idx++;
2545 }
2546
2547 // Loop on octants
2548 idx=0;
2549 obegin = m_octants.begin();
2550 oend = m_octants.end();
2551 for (it = obegin; it != oend; ++it){
2552 for (iface = 0; iface < m_dim; iface++){
2553 iface2 = iface*2;
2554 findNeighbours(m_octants.data() + idx, iface2, neighbours, isghost, false, false);
2555 nsize = neighbours.size();
2556 if (nsize) {
2557 if (!(it->m_info[iface2])){
2558 //Internal intersection
2559 for (i = 0; i < nsize; i++){
2560 if (isghost[i]){
2561 intersection.m_dim = m_dim;
2562 intersection.m_owners[0] = idx;
2563 intersection.m_owners[1] = neighbours[i];
2564 intersection.m_finer = (nsize>1);
2565 intersection.m_out = (nsize>1);
2566 intersection.m_outisghost = (nsize>1);
2567 intersection.m_iface = iface2 + (nsize>1);
2568 intersection.m_isnew = false;
2569 intersection.m_isghost = true;
2570 intersection.m_bound = false;
2571 intersection.m_pbound = true;
2572 m_intersections.push_back(intersection);
2573 }
2574 else{
2575 intersection.m_dim = m_dim;
2576 intersection.m_owners[0] = idx;
2577 intersection.m_owners[1] = neighbours[i];
2578 intersection.m_finer = (nsize>1);
2579 intersection.m_out = (nsize>1);
2580 intersection.m_outisghost = false;
2581 intersection.m_iface = iface2 + (nsize>1);
2582 intersection.m_isnew = false;
2583 intersection.m_isghost = false;
2584 intersection.m_bound = false;
2585 intersection.m_pbound = false;
2586 m_intersections.push_back(intersection);
2587 }
2588 }
2589 }
2590 else{
2591 //Periodic intersection
2592 for (i = 0; i < nsize; i++){
2593 if (isghost[i]){
2594 intersection.m_dim = m_dim;
2595 intersection.m_owners[0] = idx;
2596 intersection.m_owners[1] = neighbours[i];
2597 intersection.m_finer = (nsize>1);
2598 intersection.m_out = intersection.m_finer;
2599 intersection.m_outisghost = intersection.m_finer;
2600 intersection.m_iface = iface2 + (nsize>1);
2601 intersection.m_isnew = false;
2602 intersection.m_isghost = true;
2603 intersection.m_bound = true;
2604 intersection.m_pbound = true;
2605 m_intersections.push_back(intersection);
2606 }
2607 else{
2608 intersection.m_dim = m_dim;
2609 intersection.m_owners[0] = idx;
2610 intersection.m_owners[1] = neighbours[i];
2611 intersection.m_finer = (nsize>1);
2612 intersection.m_out = intersection.m_finer;
2613 intersection.m_outisghost = false;
2614 intersection.m_iface = iface2 + (nsize>1);
2615 intersection.m_isnew = false;
2616 intersection.m_isghost = false;
2617 intersection.m_bound = true;
2618 intersection.m_pbound = false;
2619 m_intersections.push_back(intersection);
2620 }
2621 }
2622 }
2623 }
2624 else{
2625 //Boundary intersection
2626 intersection.m_dim = m_dim;
2627 intersection.m_owners[0] = idx;
2628 intersection.m_owners[1] = idx;
2629 intersection.m_finer = 0;
2630 intersection.m_out = 0;
2631 intersection.m_outisghost = false;
2632 intersection.m_iface = iface2;
2633 intersection.m_isnew = false;
2634 intersection.m_isghost = false;
2635 intersection.m_bound = true;
2636 intersection.m_pbound = false;
2637 m_intersections.push_back(intersection);
2638 }
2639 if (it->m_info[iface2+1]){
2640 if (!(m_periodic[iface2+1])){
2641 //Boundary intersection
2642 intersection.m_dim = m_dim;
2643 intersection.m_owners[0] = idx;
2644 intersection.m_owners[1] = idx;
2645 intersection.m_finer = 0;
2646 intersection.m_out = 0;
2647 intersection.m_outisghost = false;
2648 intersection.m_iface = iface2+1;
2649 intersection.m_isnew = false;
2650 intersection.m_isghost = false;
2651 intersection.m_bound = true;
2652 intersection.m_pbound = false;
2653 m_intersections.push_back(intersection);
2654 }
2655 else{
2656 //Periodic intersection
2657 findNeighbours(m_octants.data() + idx, iface2+1, neighbours, isghost, false, false);
2658 nsize = neighbours.size();
2659 for (i = 0; i < nsize; i++){
2660 if (isghost[i]){
2661 intersection.m_dim = m_dim;
2662 intersection.m_owners[0] = idx;
2663 intersection.m_owners[1] = neighbours[i];
2664 intersection.m_finer = (nsize>1);
2665 intersection.m_out = intersection.m_finer;
2666 intersection.m_outisghost = intersection.m_finer;
2667 intersection.m_iface = iface2 + (nsize>1);
2668 intersection.m_isnew = false;
2669 intersection.m_isghost = true;
2670 intersection.m_bound = true;
2671 intersection.m_pbound = true;
2672 m_intersections.push_back(intersection);
2673 }
2674 else{
2675 intersection.m_dim = m_dim;
2676 intersection.m_owners[0] = idx;
2677 intersection.m_owners[1] = neighbours[i];
2678 intersection.m_finer = (nsize>1);
2679 intersection.m_out = intersection.m_finer;
2680 intersection.m_outisghost = false;
2681 intersection.m_iface = iface2 + (nsize>1);
2682 intersection.m_isnew = false;
2683 intersection.m_isghost = false;
2684 intersection.m_bound = true;
2685 intersection.m_pbound = false;
2686 m_intersections.push_back(intersection);
2687 }
2688 }
2689 }
2690 }
2691 }
2692 idx++;
2693 }
2694 intervector(m_intersections).swap(m_intersections);
2695 }
2696
2697 // =================================================================================== //
2702 uint32_t
2703 LocalTree::findMorton(uint64_t targetMorton) const {
2704 return findMorton(targetMorton, m_octants);
2705 };
2706
2707 // =================================================================================== //
2712 uint32_t
2713 LocalTree::findGhostMorton(uint64_t targetMorton) const {
2714 return findMorton(targetMorton, m_ghosts);
2715 };
2716
2717 // =================================================================================== //
2726 uint32_t
2727 LocalTree::findMorton(uint64_t targetMorton, const octvector &octants) const {
2728
2729 uint32_t lowerBoundIdx;
2730 uint64_t lowerBoundMorton;
2731 findMortonLowerBound(targetMorton, octants, &lowerBoundIdx, &lowerBoundMorton);
2732
2733 uint32_t targetIdx;
2734 if (lowerBoundMorton == targetMorton) {
2735 targetIdx = lowerBoundIdx;
2736 } else {
2737 targetIdx = octants.size();
2738 }
2739
2740 return targetIdx;
2741 };
2742
2743 // =================================================================================== //
2761 void
2762 LocalTree::findMortonLowerBound(uint64_t targetMorton, const octvector &octants, uint32_t *lowerBoundIdx, uint64_t *lowerBoundMorton) const {
2763
2764 uint32_t nOctants = octants.size();
2765
2766 uint32_t lowIndex = 0;
2767 uint32_t highIndex = nOctants;
2768 uint32_t midIndex = nOctants;
2769 uint64_t midMorton = PABLO::INVALID_MORTON;
2770 while (lowIndex < highIndex) {
2771 midIndex = lowIndex + (highIndex - lowIndex) / 2;
2772 midMorton = octants[midIndex].getMorton();
2773 if (targetMorton < midMorton) {
2774 highIndex = midIndex;
2775 }
2776 else if (targetMorton > midMorton) {
2777 lowIndex = midIndex + 1;
2778 }
2779 else {
2780 *lowerBoundIdx = midIndex;
2781 *lowerBoundMorton = midMorton;
2782
2783 return;
2784 }
2785 }
2786
2787 *lowerBoundIdx = lowIndex;
2788 if (*lowerBoundIdx == midIndex) {
2789 *lowerBoundMorton = midMorton;
2790 }
2791 else if (*lowerBoundIdx < nOctants) {
2792 *lowerBoundMorton = octants[*lowerBoundIdx].getMorton();
2793 }
2794 else {
2795 *lowerBoundMorton = PABLO::INVALID_MORTON;
2796 }
2797
2798 }
2799
2800 // =================================================================================== //
2816 void
2817 LocalTree::findMortonUpperBound(uint64_t targetMorton, const octvector &octants, uint32_t *upperBoundIdx, uint64_t *upperBoundMorton) const {
2818
2819 uint32_t nOctants = octants.size();
2820
2821 uint32_t lowIndex = 0;
2822 uint32_t highIndex = nOctants;
2823 uint32_t midIndex = nOctants;
2824 uint64_t midMorton = PABLO::INVALID_MORTON;
2825 while (lowIndex < highIndex) {
2826 midIndex = lowIndex + (highIndex - lowIndex) / 2;
2827 midMorton = octants[midIndex].getMorton();
2828 if (targetMorton < midMorton) {
2829 highIndex = midIndex;
2830 }
2831 else {
2832 lowIndex = midIndex + 1;
2833 }
2834 }
2835
2836 *upperBoundIdx = lowIndex;
2837 if (*upperBoundIdx == midIndex) {
2838 *upperBoundMorton = midMorton;
2839 }
2840 else if (*upperBoundIdx < nOctants) {
2841 *upperBoundMorton = octants[*upperBoundIdx].getMorton();
2842 }
2843 else {
2844 *upperBoundMorton = PABLO::INVALID_MORTON;
2845 }
2846
2847 }
2848
2849 // =================================================================================== //
2850
2853 void
2854 LocalTree::computeConnectivity(){
2855 vector<uint64_t> nodeKeys;
2856 unordered_map<uint64_t, array<uint32_t, 3> > nodeCoords;
2857 unordered_map<uint64_t, vector<uint64_t> > nodeOctants;
2858 uint32_t noctants = getNumOctants();
2859 uint32_t nghosts = getNumGhosts();
2860
2861
2862
2863 // Gather node information
2864 nodeKeys.reserve(noctants);
2865 nodeCoords.reserve(noctants);
2866 nodeOctants.reserve(noctants);
2867
2868 for (uint64_t n = 0; n < (noctants + nghosts); n++){
2869 const Octant *octant;
2870 if (n < noctants) {
2871 uint32_t octantId = static_cast<uint32_t>(n);
2872 octant = &(m_octants[octantId]);
2873 } else {
2874 uint32_t octantId = static_cast<uint32_t>(n - noctants);
2875 octant = &(m_ghosts[octantId]);
2876 }
2877
2878 for (uint8_t i = 0; i < m_treeConstants->nNodes; ++i){
2879 u32array3 node;
2880 octant->getLogicalNode(node, i);
2881
2882 uint64_t nodeKey = octant->computeNodePersistentKey(node);
2883 if (nodeCoords.count(nodeKey) == 0) {
2884 nodeKeys.push_back(nodeKey);
2885 nodeCoords.insert({{nodeKey, std::move(node)}});
2886 nodeOctants[nodeKey].reserve(m_treeConstants->nNodes);
2887 }
2888
2889 nodeOctants[nodeKey].push_back(n);
2890 }
2891 }
2892 std::sort(nodeKeys.begin(), nodeKeys.end());
2893
2894 // Build node list and connectivity
2895 m_nodes.clear();
2896 m_connectivity.clear();
2897 m_ghostsConnectivity.clear();
2898
2899 m_nodes.reserve(nodeKeys.size());
2900 m_connectivity.resize(noctants);
2901 m_ghostsConnectivity.resize(nghosts);
2902
2903 uint32_t nodeId = 0;
2904 for (uint64_t nodeKey : nodeKeys) {
2905 m_nodes.emplace_back(std::move(nodeCoords.at(nodeKey)));
2906 for (uint64_t n : nodeOctants.at(nodeKey)) {
2907 std::vector<uint32_t> *octantConnect;
2908 if (n < noctants) {
2909 uint32_t octantId = static_cast<uint32_t>(n);
2910 octantConnect = &(m_connectivity[octantId]);
2911 } else {
2912 uint32_t octantId = static_cast<uint32_t>(n - noctants);
2913 octantConnect = &(m_ghostsConnectivity[octantId]);
2914 }
2915
2916 if (octantConnect->size() == 0) {
2917 octantConnect->reserve(m_treeConstants->nNodes);
2918 }
2919 octantConnect->push_back(nodeId);
2920 }
2921 nodeId++;
2922 }
2923
2924 m_nodes.shrink_to_fit();
2925 };
2926
2932 void
2933 LocalTree::clearConnectivity(bool release){
2934 if (release) {
2935 u32arr3vector().swap(m_nodes);
2936 u32vector2D().swap(m_connectivity);
2937 u32vector2D().swap(m_ghostsConnectivity);
2938 } else {
2939 m_nodes.clear();
2940 m_connectivity.clear();
2941 m_ghostsConnectivity.clear();
2942 }
2943 };
2944
2947 void
2948 LocalTree::updateConnectivity(){
2949 clearConnectivity(false);
2950 computeConnectivity();
2951 };
2952
2953 // =================================================================================== //
2954
2955}
bitpit::Intersection
Intersection class definition.
Definition Intersection.hpp:62
bitpit::LocalTree::computeVirtualNodeNeighOffsets
void computeVirtualNodeNeighOffsets(uint8_t level, uint8_t inode, uint8_t neighLevel, std::vector< std::array< int64_t, 3 > > *neighOffsets) const
Definition LocalTree.cpp:2236
bitpit::LocalTree::getMaxEdgeNeighLevel
uint8_t getMaxEdgeNeighLevel(const Octant &octant) const
Definition LocalTree.cpp:2472
bitpit::LocalTree::computeLastVirtualEdgeNeighOffset
std::array< int64_t, 3 > computeLastVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const
Definition LocalTree.cpp:2280
bitpit::LocalTree::intervector
std::vector< Intersection > intervector
Definition LocalTree.hpp:90
bitpit::LocalTree::getMaxNeighLevel
uint8_t getMaxNeighLevel(const Octant &octant) const
Definition LocalTree.cpp:2442
bitpit::LocalTree::getMaxNodeNeighLevel
uint8_t getMaxNodeNeighLevel(const Octant &octant) const
Definition LocalTree.cpp:2457
bitpit::LocalTree::getEdgePeriodicOffset
std::array< int64_t, 3 > getEdgePeriodicOffset(const Octant &octant, uint8_t iedge) const
Definition LocalTree.cpp:2417
bitpit::LocalTree::getPeriodicOffset
std::array< int64_t, 3 > getPeriodicOffset(const Octant &octant, uint8_t iface) const
Definition LocalTree.cpp:2368
bitpit::LocalTree::computeLastVirtualNeighOffset
std::array< int64_t, 3 > computeLastVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const
Definition LocalTree.cpp:2099
bitpit::LocalTree::u32vector
std::vector< uint32_t > u32vector
Definition LocalTree.hpp:102
bitpit::LocalTree::computeVirtualEdgeNeighOffsets
void computeVirtualEdgeNeighOffsets(uint8_t level, uint8_t iedge, uint8_t neighLevel, std::vector< std::array< int64_t, 3 > > *neighOffsets) const
Definition LocalTree.cpp:2324
bitpit::LocalTree::u32arr3vector
std::vector< u32array3 > u32arr3vector
Definition LocalTree.hpp:110
bitpit::LocalTree::computeLastVirtualNodeNeighOffset
std::array< int64_t, 3 > computeLastVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const
Definition LocalTree.cpp:2223
bitpit::LocalTree::computeFirstVirtualNeighOffset
std::array< int64_t, 3 > computeFirstVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const
Definition LocalTree.cpp:2068
bitpit::LocalTree::computeFirstVirtualEdgeNeighOffset
std::array< int64_t, 3 > computeFirstVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const
Definition LocalTree.cpp:2249
bitpit::LocalTree::computeVirtualNeighOffsets
void computeVirtualNeighOffsets(uint8_t level, uint8_t iface, uint8_t neighLevel, std::vector< std::array< int64_t, 3 > > *neighOffsets) const
Definition LocalTree.cpp:2140
bitpit::LocalTree::computeFirstVirtualNodeNeighOffset
std::array< int64_t, 3 > computeFirstVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const
Definition LocalTree.cpp:2192
bitpit::LocalTree::getNodePeriodicOffset
std::array< int64_t, 3 > getNodePeriodicOffset(const Octant &octant, uint8_t inode) const
Definition LocalTree.cpp:2392
bitpit::Octant
Octant class definition.
Definition Octant.hpp:89
bitpit::Octant::getLevel
uint8_t getLevel() const
Definition Octant.cpp:259
bitpit::Octant::getBalance
bool getBalance() const
Definition Octant.cpp:379
abs
std::array< T, d > abs(const std::array< T, d > &x)
Definition MathOperators_array.tpp:591
bitpit::dgf::check
unsigned int check(std::ifstream &, std::vector< std::vector< int > > &)
Definition DGF.cpp:1169
bitpit::TreeConstants::lengths
std::vector< uint32_t > lengths
Definition tree_constants.hpp:95
bitpit::TreeConstants::instance
static BITPIT_PUBLIC_API const TreeConstants & instance(uint8_t dim)
Definition tree_constants.cpp:41
--- layout: doxygen_footer ---
Generated on Wed Jun 19 2024 10:21:05 for bitpit by   doxygen 1.10.0