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 // Set children information
528 for (int i = 0; i < nChildren; ++i) {
529 m_octants[firstChildIdx + i].m_info[Octant::INFO_NEW4REFINEMENT] = true;
530 }
531
532 // Update the mapping
533 if(!mapidx.empty()){
534 for (uint8_t i=0; i<nChildren; i++){
535 mapidx[firstChildIdx + i] = mapidx[idx];
536 }
537 }
538
539 // Check if more refinement is needed to satisfy the markers
540 if (refinementCompleted) {
541 refinementCompleted = (m_octants[firstChildIdx].getMarker() <= 0);
542 }
543
544 // Update local max depth
545 uint8_t childrenLevel = m_octants[firstChildIdx].getLevel();
546 if (childrenLevel > m_localMaxDepth){
547 m_localMaxDepth = childrenLevel;
548 }
549
550 // Update future octant index
551 futureIdx -= nChildren;
552 }
553 }
554
555 return (!refinementCompleted);
556
557 };
558
559 // =================================================================================== //
565 bool
566 LocalTree::coarse(u32vector & mapidx){
567
568 u32vector first_child_index;
569 Octant father(m_dim);
570 uint32_t idx, idx2;
571 uint32_t offset;
572 uint32_t idx1_gh;
573 uint32_t idx2_gh;
574 uint32_t nidx;
575 uint32_t mapsize = mapidx.size();
576 int8_t markerfather, marker;
577 uint8_t nbro, nend, nstart;
578 uint8_t nchm1 = m_treeConstants->nChildren-1;
579 bool docoarse = false;
580 bool wstop = false;
581
582 //------------------------------------------ //
583 // Initialization
584
585 uint32_t nInitialOctants = getNumOctants();
586 if (nInitialOctants == 0) {
587 return false;
588 }
589
590 uint32_t nInitialGhosts = getNumGhosts();
591
592 nbro = nend = nstart = 0;
593 nidx = offset = 0;
594
595 idx2_gh = 0;
596 idx1_gh = nInitialGhosts - 1;
597
598 //------------------------------------------ //
599
600 // Set index for start and end check for ghosts
601 if (m_ghosts.size()){
602 bool check = true;
603 while(check){
604 check = idx1_gh < nInitialGhosts;
605 if (check){
606 check = m_ghosts[idx1_gh].getMorton() > m_firstDescMorton;
607 }
608 if (check) idx1_gh--;
609 }
610
611 check = true;
612 while(check){
613 check = idx2_gh < nInitialGhosts;
614 if (check){
615 check = m_ghosts[idx2_gh].getMorton() < m_lastDescMorton;
616 }
617 if (check) idx2_gh++;
618 }
619 }
620
621 // Check and coarse internal octants
622 for (idx=0; idx<nInitialOctants; idx++){
623 if(m_octants[idx].getMarker() < 0 && m_octants[idx].getLevel() > 0){
624 nbro = 0;
625 father = m_octants[idx].buildFather();
626 // Check if family is to be refined
627 for (idx2=idx; idx2<idx+m_treeConstants->nChildren; idx2++){
628 if (idx2<nInitialOctants){
629 if(m_octants[idx2].getMarker() < 0 && m_octants[idx2].buildFather() == father){
630 nbro++;
631 }
632 }
633 }
634 if (nbro == m_treeConstants->nChildren){
635 nidx++;
636 first_child_index.push_back(idx);
637 idx = idx2-1;
638 }
639 }
640 }
641 uint32_t nblock = nInitialOctants;
642 uint32_t nfchild = first_child_index.size();
643 if (nidx!=0){
644 nblock = nInitialOctants - nidx*nchm1;
645 nidx = 0;
646 for (idx=0; idx<nblock; idx++){
647 if (idx+offset < nInitialOctants){
648 if (nidx < nfchild){
649 if (idx+offset == first_child_index[nidx]){
650 markerfather = -m_treeConstants->maxLevel;
651 father = m_octants[idx+offset].buildFather();
652 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT; iii++){
653 father.m_info[iii] = false;
654 }
655 father.setGhostLayer(-1);
656 for(idx2=0; idx2<m_treeConstants->nChildren; idx2++){
657 if (idx2 < nInitialOctants){
658 if (markerfather < m_octants[idx+offset+idx2].getMarker()+1){
659 markerfather = m_octants[idx+offset+idx2].getMarker()+1;
660 }
661 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT; iii++){
662 father.m_info[iii] = father.m_info[iii] || m_octants[idx+offset+idx2].m_info[iii];
663 }
664 }
665 }
666 father.m_info[Octant::INFO_NEW4COARSENING] = true;
667 father.setMarker(markerfather);
668 //Impossible in this version
669// if (markerfather < 0 && mapsize == 0){
670// docoarse = true;
671// }
672 m_octants[idx] = father;
673 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
674 offset += nchm1;
675 nidx++;
676 }
677 else{
678 m_octants[idx] = m_octants[idx+offset];
679 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
680 }
681 }
682 else{
683 m_octants[idx] = m_octants[idx+offset];
684 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
685 }
686 }
687 }
688 }
689 m_octants.resize(nblock, Octant(m_dim));
690 m_octants.shrink_to_fit();
691 nInitialOctants = m_octants.size();
692 if(mapsize > 0){
693 mapidx.resize(nInitialOctants);
694 }
695
696 //Check ghosts
697 if (m_ghosts.size()){
698 // Start on ghosts
699 if (nInitialOctants > 0 && idx1_gh < nInitialGhosts){
700 if (m_ghosts[idx1_gh].buildFather() == m_octants[0].buildFather()){
701 father = m_ghosts[idx1_gh].buildFather();
702 nbro = 0;
703 idx = idx1_gh;
704 marker = m_ghosts[idx].getMarker();
705 bool waszero = (idx == 0);
706 while(marker < 0 && m_ghosts[idx].buildFather() == father){
707 nbro++;
708 if(waszero){
709 break;
710 }
711 idx--;
712 if(idx == 0){
713 waszero = true;
714 }
715 marker = m_ghosts[idx].getMarker();
716 }
717 nstart = 0;
718 idx = 0;
719 marker = m_octants[idx].getMarker();
720 if (idx==nInitialOctants-1) wstop = true;
721 while(marker < 0 && m_octants[idx].buildFather() == father){
722 nbro++;
723 nstart++;
724 if (wstop){
725 break;
726 }
727 idx++;
728 if (idx==nInitialOctants-1){
729 wstop = true;
730 }
731 marker = m_octants[idx].getMarker();
732 }
733 if (nbro == m_treeConstants->nChildren){
734 offset = nstart;
735 }
736 else{
737 nstart = 0;
738 }
739 }
740 if (nstart != 0){
741 for (idx=0; idx<nInitialOctants; idx++){
742 if (idx+offset < nInitialOctants){
743 m_octants[idx] = m_octants[idx+offset];
744 if(mapsize > 0) mapidx[idx] = mapidx[idx+offset];
745 }
746 }
747 m_octants.resize(nInitialOctants-offset, Octant(m_dim));
748 m_octants.shrink_to_fit();
749 nInitialOctants = m_octants.size();
750 if(mapsize > 0){
751 mapidx.resize(nInitialOctants);
752 }
753 }
754 }
755
756
757 //Verify family between more then two processes
758 if (nInitialOctants > 0 && idx2_gh < nInitialGhosts){
759
760 if (m_ghosts[idx2_gh].buildFather() == father){
761
762 if (m_ghosts[idx2_gh].buildFather() == m_octants[nInitialOctants-1].buildFather()){
763
764 uint64_t idx22_gh = idx2_gh;
765 marker = m_ghosts[idx22_gh].getMarker();
766 while(marker < 0 && m_ghosts[idx22_gh].buildFather() == father){
767 nbro++;
768 idx22_gh++;
769 if(idx22_gh == nInitialGhosts){
770 break;
771 }
772 marker = m_ghosts[idx22_gh].getMarker();
773 }
774
775 if (nbro == m_treeConstants->nChildren){
776 m_octants.clear();
777 nInitialOctants = 0;
778 if(mapsize > 0){
779 mapidx.clear();
780 }
781 }
782 }
783
784 }
785
786 }
787
788
789 // End on ghosts
790 if (nInitialOctants > 0 && idx2_gh < nInitialGhosts){
791 if (m_ghosts[idx2_gh].buildFather() == m_octants[nInitialOctants-1].buildFather()){
792 father = m_ghosts[idx2_gh].buildFather();
793 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT; iii++){
794 father.m_info[iii] = false;
795 }
796 father.setGhostLayer(-1);
797 markerfather = m_ghosts[idx2_gh].getMarker()+1;
798 nbro = 0;
799 idx = idx2_gh;
800 marker = m_ghosts[idx].getMarker();
801 while(marker < 0 && m_ghosts[idx].buildFather() == father){
802 nbro++;
803 if (markerfather < m_ghosts[idx].getMarker()+1){
804 markerfather = m_ghosts[idx].getMarker()+1;
805 }
806 for (uint32_t iii=0; iii<m_treeConstants->nFaces; iii++){
807 father.m_info[iii] = father.m_info[iii] || m_ghosts[idx].m_info[iii];
808 }
809 father.m_info[Octant::INFO_BALANCED] = father.m_info[Octant::INFO_BALANCED] || m_ghosts[idx].m_info[Octant::INFO_BALANCED];
810 idx++;
811 if(idx == nInitialGhosts){
812 break;
813 }
814 marker = m_ghosts[idx].getMarker();
815 }
816 nend = 0;
817 idx = nInitialOctants-1;
818 marker = m_octants[idx].getMarker();
819 if (idx==0) wstop = true;
820 while(marker < 0 && m_octants[idx].buildFather() == father){
821 nbro++;
822 nend++;
823 if (markerfather < m_octants[idx].getMarker()+1){
824 markerfather = m_octants[idx].getMarker()+1;
825 }
826 idx--;
827 if (wstop){
828 break;
829 }
830 if (idx==0){
831 wstop = true;
832 }
833 marker = m_octants[idx].getMarker();
834 }
835 if (nbro == m_treeConstants->nChildren){
836 offset = nend;
837 }
838 else{
839 nend = 0;
840 }
841 }
842 if (nend != 0){
843 for (idx=0; idx < nend; idx++){
844 for (uint32_t iii=0; iii<Octant::INFO_ITEM_COUNT - 1; iii++){
845 father.m_info[iii] = father.m_info[iii] || m_octants[nInitialOctants-idx-1].m_info[iii];
846 }
847 }
848 father.m_info[Octant::INFO_NEW4COARSENING] = true;
849 father.setGhostLayer(-1);
850 //Impossible in this version
851 // if (markerfather < 0 && mapsize == 0){
852 // docoarse = true;
853 // }
854 father.setMarker(markerfather);
855 m_octants.resize(nInitialOctants-offset, Octant(m_dim));
856 m_octants.push_back(father);
857 m_octants.shrink_to_fit();
858 nInitialOctants = m_octants.size();
859 if(mapsize > 0){
860 mapidx.resize(nInitialOctants);
861 }
862 }
863 }
864
865 }//end if ghosts size
866
867 // Update maximum depth
868 updateLocalMaxDepth();
869
870 // Set final last desc
871 setFirstDescMorton();
872 setLastDescMorton();
873
874 return docoarse;
875
876 };
877
878 // =================================================================================== //
879
884 bool
885 LocalTree::globalRefine(u32vector & mapidx){
886
887 for (Octant &octant : m_octants){
888 octant.setMarker(1);
889 }
890
891 return refine(mapidx);
892
893 };
894
895 // =================================================================================== //
896
901 bool
902 LocalTree::globalCoarse(u32vector & mapidx){
903
904 for (Octant &octant : m_octants){
905 octant.setMarker(-1);
906 }
907 for (Octant &octant : m_ghosts){
908 octant.setMarker(-1);
909 }
910
911 return coarse(mapidx);
912
913 };
914
920 void
921 LocalTree::checkCoarse(uint64_t partLastDesc, u32vector & mapidx){
922
923 uint32_t idx;
924 uint32_t mapsize = mapidx.size();
925 uint8_t toDelete = 0;
926
927 uint32_t nInitialOctants = getNumOctants();
928 if (nInitialOctants>0){
929
930 idx = 0;
931 if (m_octants[idx].getMorton() < partLastDesc){
932
933 Octant father0 = m_octants[idx].buildFather();
934 Octant father = father0;
935
936 while(father == father0 && idx < nInitialOctants){
937 toDelete++;
938 idx++;
939 if (idx<nInitialOctants) father = m_octants[idx].buildFather();
940 }
941
942 if (nInitialOctants>toDelete){
943 for(idx=0; idx<nInitialOctants-toDelete; idx++){
944 m_octants[idx] = m_octants[idx+toDelete];
945 if (mapsize>0) mapidx[idx] = mapidx[idx+toDelete];
946 }
947 m_octants.resize(nInitialOctants-toDelete, Octant(m_dim));
948 if (mapsize>0){
949 mapidx.resize(m_octants.size());
950 }
951 }
952 else{
953 m_octants.clear();
954 mapidx.clear();
955 }
956 setFirstDescMorton();
957 }
958
959 }
960 };
961
962 // =================================================================================== //
965 void
966 LocalTree::updateLocalMaxDepth(){
967
968 uint32_t noctants = getNumOctants();
969 uint32_t i;
970
971 m_localMaxDepth = 0;
972 for(i = 0; i < noctants; i++){
973 if(m_octants[i].getLevel() > m_localMaxDepth){
974 m_localMaxDepth = m_octants[i].getLevel();
975 }
976 }
977 };
978
979
991 void
992 LocalTree::findNeighbours(const Octant* oct, uint8_t iface, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
993
994 if (!append) {
995 isghost.clear();
996 neighbours.clear();
997 }
998
999 // Default if iface is nface<iface<0
1000 if (iface >= m_treeConstants->nFaces){
1001 return;
1002 }
1003
1004 // If a face is a boundary, it can have neighbours only if periodic
1005 bool isperiodic = false;
1006 if (oct->getBound(iface)) {
1007 isperiodic = isPeriodic(oct, iface);
1008 if (!isperiodic) {
1009 return;
1010 }
1011 }
1012
1013 // Exit if is the root octant
1014 uint8_t level = oct->getLevel();
1015 if (level == 0){
1016 // if periodic face return itself
1017 if (isperiodic){
1018 neighbours.push_back(0);
1019 isghost.push_back(false);
1020 }
1021 return;
1022 }
1023
1024 // Initialize search
1025 uint32_t candidateIdx = 0;
1026 uint64_t candidateMorton = 0;
1027
1028 uint64_t neighArea = 0;
1029 uint64_t faceArea = oct->getLogicalArea();
1030
1031 uint32_t size = oct->getLogicalSize();
1032
1033 std::array<uint32_t, 3> coord = oct->getLogicalCoordinates();
1034 if (isperiodic) {
1035 std::array<int64_t, 3> periodicOffset = getPeriodicOffset(*oct, iface);
1036 coord[0] = static_cast<uint32_t>(coord[0] + periodicOffset[0]);
1037 coord[1] = static_cast<uint32_t>(coord[1] + periodicOffset[1]);
1038 coord[2] = static_cast<uint32_t>(coord[2] + periodicOffset[2]);
1039 }
1040
1041 const int8_t (&cxyz)[3] = m_treeConstants->normals[iface];
1042
1043 // Compute same-size virtual neighbour information
1044 std::array<int64_t, 3> sameSizeVirtualNeighOffset = computeFirstVirtualNeighOffset(level, iface, level);
1045 std::array<uint32_t, 3> sameSizeVirtualNeighCoord = coord;
1046 sameSizeVirtualNeighCoord[0] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[0] + sameSizeVirtualNeighOffset[0]);
1047 sameSizeVirtualNeighCoord[1] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[1] + sameSizeVirtualNeighOffset[1]);
1048 sameSizeVirtualNeighCoord[2] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[2] + sameSizeVirtualNeighOffset[2]);
1049 uint64_t sameSizeVirtualNeighMorton = PABLO::computeMorton(m_dim, sameSizeVirtualNeighCoord[0], sameSizeVirtualNeighCoord[1], sameSizeVirtualNeighCoord[2]);
1050
1051 //
1052 // Search in the internal octants
1053 //
1054
1055 // Identify the index of the first neighbour candidate
1056 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_octants, &candidateIdx, &candidateMorton);
1057
1058 // Early return if a neighbour of the same size has been found
1059 if(candidateMorton == sameSizeVirtualNeighMorton && m_octants[candidateIdx].m_level == level){
1060 neighbours.push_back(candidateIdx);
1061 isghost.push_back(false);
1062 return;
1063 }
1064
1065 // Compute the Morton number of the last candidate
1066 //
1067 // This is the Morton number of the last descendant of the same-size
1068 // virtual neighbour.
1069 uint8_t maxNeighLevel = getMaxNeighLevel(*oct);
1070 std::array<int64_t, 3> lastCandidateOffset = computeLastVirtualNeighOffset(level, iface, maxNeighLevel);
1071 std::array<uint32_t, 3> lastCandidateCoord = coord;
1072 lastCandidateCoord[0] = static_cast<uint32_t>(lastCandidateCoord[0] + lastCandidateOffset[0]);
1073 lastCandidateCoord[1] = static_cast<uint32_t>(lastCandidateCoord[1] + lastCandidateOffset[1]);
1074 lastCandidateCoord[2] = static_cast<uint32_t>(lastCandidateCoord[2] + lastCandidateOffset[2]);
1075 uint64_t lastCandidateMorton = PABLO::computeMorton(m_dim, lastCandidateCoord[0], lastCandidateCoord[1], lastCandidateCoord[2]);
1076
1077 // Search for neighbours of different sizes
1078 if (candidateIdx < getNumOctants()){
1079 while(true){
1080 // Detect if the candidate is a neighbour
1081 u32array3 coordtry = {{0, 0, 0}};
1082 bool isNeighbourCandidate = true;
1083 for (int8_t idim=0; idim<m_dim; idim++){
1084 coordtry[idim] = m_octants[candidateIdx].getLogicalCoordinates(idim);
1085
1086 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(coordtry[idim] - coord[idim]));
1087 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_octants[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1088 if (Dx != Dxstar){
1089 isNeighbourCandidate = false;
1090 break;
1091 }
1092 }
1093
1094 if (isNeighbourCandidate){
1095 bool isNeighbour = false;
1096 uint8_t leveltry = m_octants[candidateIdx].getLevel();
1097 if (leveltry > level){
1098 array<int64_t,3> coord1 ={{1, 1, 1}} ;
1099 for (int8_t idim=0; idim<m_dim; idim++){
1100 coord1[idim] = coord[idim] + size;
1101 }
1102
1103 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])))){
1104 isNeighbour = true;
1105 }
1106 }
1107 else if (leveltry < level){
1108 u32array3 coordtry1 = {{1, 1, 1}};
1109 for (int8_t idim=0; idim<m_dim; idim++){
1110 coordtry1[idim] = coordtry[idim] + m_octants[candidateIdx].getLogicalSize();
1111 }
1112
1113 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])))){
1114 isNeighbour = true;
1115 }
1116 }
1117
1118 if (isNeighbour){
1119 neighbours.push_back(candidateIdx);
1120 isghost.push_back(false);
1121
1122 // If the neighbours already cover the whole face, we have
1123 // found all the neighbours and we can exit.
1124 neighArea += m_octants[candidateIdx].getLogicalArea();
1125 if (neighArea == faceArea){
1126 return;
1127 }
1128 }
1129 }
1130
1131 // Advance to the next candidate
1132 candidateIdx++;
1133 if (candidateIdx > getNumOctants() - 1){
1134 break;
1135 }
1136
1137 candidateMorton = m_octants[candidateIdx].getMorton();
1138 if (candidateMorton > lastCandidateMorton){
1139 break;
1140 }
1141 }
1142 }
1143
1144 //
1145 // Search in the ghost octants
1146 //
1147 // If we want also the neighbours that are ghosts, we always need to
1148 // search in the ghosts, the only exception is for faces of internal
1149 // octants that are not process boundaries.
1150 bool ghostSearch = !onlyinternal && (getNumGhosts() > 0);
1151 if (ghostSearch){
1152 if (!oct->getIsGhost() && !oct->getPbound(iface)){
1153 ghostSearch = false;
1154 }
1155 }
1156
1157 // Search in ghosts
1158 if(ghostSearch){
1159 // Identify the index of the first neighbour candidate
1160 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_ghosts, &candidateIdx, &candidateMorton);
1161
1162 // Early return if a neighbour of the same size has been found
1163 if(candidateMorton == sameSizeVirtualNeighMorton && m_ghosts[candidateIdx].getLevel() == level){
1164 neighbours.push_back(candidateIdx);
1165 isghost.push_back(true);
1166 return;
1167 }
1168
1169 // Search for neighbours of different sizes
1170 if (candidateIdx < getNumGhosts()){
1171 while(true){
1172 // Detect if the candidate is a neighbour
1173 u32array3 coordtry = {{0, 0, 0}};
1174 bool isNeighbourCandidate = true;
1175 for (int8_t idim=0; idim<m_dim; idim++){
1176 coordtry[idim] = m_ghosts[candidateIdx].getLogicalCoordinates(idim);
1177
1178 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(coordtry[idim] - coord[idim]));
1179 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_ghosts[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1180 if (Dx != Dxstar){
1181 isNeighbourCandidate = false;
1182 break;
1183 }
1184 }
1185
1186 if (isNeighbourCandidate){
1187 bool isNeighbour = false;
1188 uint8_t leveltry = m_ghosts[candidateIdx].getLevel();
1189 if (leveltry > level){
1190 array<int64_t, 3> coord1 = {{1, 1, 1}};
1191 for (int8_t idim=0; idim<m_dim; idim++){
1192 coord1[idim] = coord[idim] + size;
1193 }
1194
1195 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])))){
1196 isNeighbour = true;
1197 }
1198 }
1199 else if (leveltry < level){
1200 u32array3 coordtry1 = {{1, 1, 1}};
1201 for (int8_t idim=0; idim<m_dim; idim++){
1202 coordtry1[idim] = coordtry[idim] + m_ghosts[candidateIdx].getLogicalSize();
1203 }
1204
1205 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])))){
1206 isNeighbour = true;
1207 }
1208 }
1209
1210 if (isNeighbour){
1211 neighbours.push_back(candidateIdx);
1212 isghost.push_back(true);
1213
1214 // If the neighbours already cover the whole face, we have
1215 // found all the neighbours and we can exit.
1216 neighArea += m_ghosts[candidateIdx].getLogicalArea();
1217 if (neighArea == faceArea){
1218 return;
1219 }
1220 }
1221 }
1222
1223 candidateIdx++;
1224 if (candidateIdx > getNumGhosts() - 1){
1225 break;
1226 }
1227
1228 candidateMorton = m_ghosts[candidateIdx].getMorton();
1229 if (candidateMorton > lastCandidateMorton){
1230 break;
1231 }
1232 }
1233 }
1234 }
1235 };
1236
1248 void
1249 LocalTree::findEdgeNeighbours(const Octant* oct, uint8_t iedge, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
1250
1251 if (!append) {
1252 isghost.clear();
1253 neighbours.clear();
1254 }
1255
1256 // Default if iedge is nface<iedge<0
1257 if (iedge >= m_treeConstants->nEdges){
1258 return;
1259 }
1260
1261 // If a edge is a on boundary, it can have neighbours only if periodic
1262 bool isperiodic = false;
1263 if (oct->getEdgeBound(iedge)) {
1264 isperiodic = isEdgePeriodic(oct, iedge);
1265 if (!isperiodic) {
1266 return;
1267 }
1268 }
1269
1270 // Exit if is the root octant
1271 uint8_t level = oct->getLevel();
1272 if (level == 0){
1273 // if periodic face return itself
1274 if (isperiodic){
1275 neighbours.push_back(0);
1276 isghost.push_back(false);
1277 }
1278 return;
1279 }
1280
1281 // Search in the octants
1282 uint32_t candidateIdx = 0;
1283 uint64_t candidateMorton = 0;
1284
1285 uint32_t neighSize = 0;
1286 uint32_t edgeSize = oct->getLogicalSize();
1287
1288 uint32_t size = oct->getLogicalSize();
1289
1290 std::array<uint32_t, 3> coord = oct->getLogicalCoordinates();
1291 if (isperiodic) {
1292 std::array<int64_t, 3> periodicOffset = getEdgePeriodicOffset(*oct, iedge);
1293 coord[0] = static_cast<uint32_t>(coord[0] + periodicOffset[0]);
1294 coord[1] = static_cast<uint32_t>(coord[1] + periodicOffset[1]);
1295 coord[2] = static_cast<uint32_t>(coord[2] + periodicOffset[2]);
1296 }
1297
1298 const int8_t (&cxyz)[3] = m_treeConstants->edgeCoeffs[iedge];
1299
1300 // Compute same-size virtual neighbour information
1301 std::array<int64_t, 3> sameSizeVirtualNeighOffset = computeFirstVirtualEdgeNeighOffset(level, iedge, level);
1302 std::array<uint32_t, 3> sameSizeVirtualNeighCoord = coord;
1303 sameSizeVirtualNeighCoord[0] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[0] + sameSizeVirtualNeighOffset[0]);
1304 sameSizeVirtualNeighCoord[1] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[1] + sameSizeVirtualNeighOffset[1]);
1305 sameSizeVirtualNeighCoord[2] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[2] + sameSizeVirtualNeighOffset[2]);
1306 uint64_t sameSizeVirtualNeighMorton = PABLO::computeMorton(m_dim, sameSizeVirtualNeighCoord[0], sameSizeVirtualNeighCoord[1], sameSizeVirtualNeighCoord[2]);
1307
1308 //
1309 // Search in the internal octants
1310 //
1311
1312 // Identify the index of the first neighbour candidate
1313 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_octants, &candidateIdx, &candidateMorton);
1314
1315 // Early return if a neighbour of the same size has been found
1316 if(candidateMorton == sameSizeVirtualNeighMorton && m_octants[candidateIdx].m_level == level){
1317 isghost.push_back(false);
1318 neighbours.push_back(candidateIdx);
1319 return;
1320 }
1321
1322 // Compute the Morton number of the last candidate
1323 //
1324 // This is the Morton number of the last descendant of the same-size
1325 // virtual neighbour.
1326 uint8_t maxEdgeNeighLevel = getMaxEdgeNeighLevel(*oct);
1327 std::array<int64_t, 3> lastCandidateOffset = computeLastVirtualEdgeNeighOffset(level, iedge, maxEdgeNeighLevel);
1328 std::array<uint32_t, 3> lastCandidateCoord = coord;
1329 lastCandidateCoord[0] = static_cast<uint32_t>(lastCandidateCoord[0] + lastCandidateOffset[0]);
1330 lastCandidateCoord[1] = static_cast<uint32_t>(lastCandidateCoord[1] + lastCandidateOffset[1]);
1331 lastCandidateCoord[2] = static_cast<uint32_t>(lastCandidateCoord[2] + lastCandidateOffset[2]);
1332 uint64_t lastCandidateMorton = PABLO::computeMorton(m_dim, lastCandidateCoord[0], lastCandidateCoord[1], lastCandidateCoord[2]);
1333
1334 // Search for neighbours of different sizes
1335 if (candidateIdx < getNumOctants()) {
1336 while(true){
1337 // Detect if the candidate is a neighbour
1338 u32array3 coordtry = {{0, 0, 0}};
1339 bool isNeighbourCandidate = true;
1340 for (int8_t idim=0; idim<m_dim; idim++){
1341 coordtry[idim] = m_octants[candidateIdx].getLogicalCoordinates(idim);
1342
1343 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1344 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_octants[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1345
1346 if (Dx != Dxstar) {
1347 isNeighbourCandidate = false;
1348 break;
1349 }
1350 }
1351
1352 if (isNeighbourCandidate){
1353 bool isNeighbour = false;
1354 uint8_t leveltry = m_octants[candidateIdx].getLevel();
1355 if (leveltry > level){
1356 u32array3 coord1 = {{1, 1, 1}};
1357 for (int8_t idim=0; idim<m_dim; idim++){
1358 coord1[idim] = coord[idim] + size;
1359 }
1360
1361 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])))){
1362 isNeighbour = true;
1363 }
1364 }
1365 else if (leveltry < level){
1366 u32array3 coordtry1 = {{1, 1, 1}};
1367 for (int8_t idim=0; idim<m_dim; idim++){
1368 coordtry1[idim] = coordtry[idim] + m_octants[candidateIdx].getLogicalSize();
1369 }
1370
1371 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])))){
1372 isNeighbour = true;
1373 }
1374 }
1375
1376 if (isNeighbour) {
1377 neighbours.push_back(candidateIdx);
1378 isghost.push_back(false);
1379
1380 // If the neighbours already cover the whole edge, we have
1381 // found all the neighbours and we can exit.
1382 neighSize += m_octants[candidateIdx].getLogicalSize();
1383 if (neighSize == edgeSize){
1384 return;
1385 }
1386 }
1387 }
1388
1389 // Advance to the next candidate
1390 candidateIdx++;
1391 if (candidateIdx > getNumOctants()-1){
1392 break;
1393 }
1394
1395 candidateMorton = m_octants[candidateIdx].getMorton();
1396 if (candidateMorton > lastCandidateMorton){
1397 break;
1398 }
1399 }
1400 }
1401
1402 //
1403 // Search in the ghost octants
1404 //
1405 if (getNumGhosts() > 0 && !onlyinternal){
1406 // Identify the index of the first neighbour candidate
1407 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_ghosts, &candidateIdx, &candidateMorton);
1408
1409 // Early return if a neighbour of the same size has been found
1410 if(candidateMorton == sameSizeVirtualNeighMorton && m_ghosts[candidateIdx].m_level == level){
1411 isghost.push_back(true);
1412 neighbours.push_back(candidateIdx);
1413 return;
1414 }
1415
1416 // Search for neighbours of different sizes
1417 if (candidateIdx < getNumGhosts()){
1418 while(true){
1419 // Detect if the candidate is a neighbour
1420 u32array3 coordtry = {{0, 0, 0}};
1421 bool isNeighbourCandidate = true;
1422 for (int8_t idim=0; idim<m_dim; idim++){
1423 coordtry[idim] = m_ghosts[candidateIdx].getLogicalCoordinates(idim);
1424
1425 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1426 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_ghosts[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1427
1428 if (Dx != Dxstar) {
1429 isNeighbourCandidate = false;
1430 break;
1431 }
1432 }
1433
1434 if (isNeighbourCandidate){
1435 bool isNeighbour = false;
1436 uint8_t leveltry = m_ghosts[candidateIdx].getLevel();
1437 if (leveltry > level){
1438 u32array3 coord1 = {{1, 1, 1}};
1439 for (int8_t idim=0; idim<m_dim; idim++){
1440 coord1[idim] = int32_t(coord[idim]) + size;
1441 }
1442
1443 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])))){
1444 isNeighbour = true;
1445 }
1446 }
1447 else if (leveltry < level){
1448 u32array3 coordtry1 = {{1, 1, 1}};
1449 for (int8_t idim=0; idim<m_dim; idim++){
1450 coordtry1[idim] = coordtry[idim] + m_ghosts[candidateIdx].getLogicalSize();
1451 }
1452
1453 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])))){
1454 isNeighbour = true;
1455 }
1456 }
1457
1458 if (isNeighbour) {
1459 neighbours.push_back(candidateIdx);
1460 isghost.push_back(true);
1461
1462 // If the neighbours already cover the whole edge, we have
1463 // found all the neighbours and we can exit.
1464 neighSize += m_ghosts[candidateIdx].getLogicalSize();
1465 if (neighSize == edgeSize){
1466 return;
1467 }
1468 }
1469 }
1470
1471 // Advance to the next candidate
1472 candidateIdx++;
1473 if (candidateIdx > getNumGhosts()-1){
1474 break;
1475 }
1476
1477 candidateMorton = m_ghosts[candidateIdx].getMorton();
1478 if (candidateMorton > lastCandidateMorton){
1479 break;
1480 }
1481 }
1482 }
1483 }
1484 };
1485
1497 void
1498 LocalTree::findNodeNeighbours(const Octant* oct, uint8_t inode, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
1499
1500 if (!append) {
1501 isghost.clear();
1502 neighbours.clear();
1503 }
1504
1505 // Default if inode is nnodes<inode<0
1506 if (inode >= m_treeConstants->nNodes){
1507 return;
1508 }
1509
1510 // If a node is a on boundary, it can have neighbours only if periodic
1511 bool isperiodic = false;
1512 if (oct->getNodeBound(inode)) {
1513 isperiodic = isNodePeriodic(oct, inode);
1514 if (!isperiodic) {
1515 return;
1516 }
1517 }
1518
1519 // Exit if is the root octant
1520 uint8_t level = oct->getLevel();
1521 if (level == 0){
1522 // if periodic face return itself
1523 if (isperiodic){
1524 neighbours.push_back(0);
1525 isghost.push_back(false);
1526 }
1527 return;
1528 }
1529
1530 // Initialize search
1531 uint32_t candidateIdx = 0;
1532 uint64_t candidateMorton = 0;
1533
1534 uint32_t size = oct->getLogicalSize();
1535
1536 std::array<uint32_t, 3> coord = oct->getLogicalCoordinates();
1537 if (isperiodic) {
1538 std::array<int64_t, 3> periodicOffset = getNodePeriodicOffset(*oct, inode);
1539 coord[0] = static_cast<uint32_t>(coord[0] + periodicOffset[0]);
1540 coord[1] = static_cast<uint32_t>(coord[1] + periodicOffset[1]);
1541 coord[2] = static_cast<uint32_t>(coord[2] + periodicOffset[2]);
1542 }
1543
1544 const int8_t (&cxyz)[3] = m_treeConstants->nodeCoeffs[inode];
1545
1546 // Compute same-size virtual neighbour information
1547 std::array<int64_t, 3> sameSizeVirtualNeighOffset = computeFirstVirtualNodeNeighOffset(level, inode, level);
1548 std::array<uint32_t, 3> sameSizeVirtualNeighCoord = coord;
1549 sameSizeVirtualNeighCoord[0] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[0] + sameSizeVirtualNeighOffset[0]);
1550 sameSizeVirtualNeighCoord[1] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[1] + sameSizeVirtualNeighOffset[1]);
1551 sameSizeVirtualNeighCoord[2] = static_cast<uint32_t>(sameSizeVirtualNeighCoord[2] + sameSizeVirtualNeighOffset[2]);
1552 uint64_t sameSizeVirtualNeighMorton = PABLO::computeMorton(m_dim, sameSizeVirtualNeighCoord[0], sameSizeVirtualNeighCoord[1], sameSizeVirtualNeighCoord[2]);
1553
1554 //
1555 // Search in the internal octants
1556 //
1557
1558 // Identify the index of the first neighbour candidate
1559 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_octants, &candidateIdx, &candidateMorton);
1560
1561 // Early return if a neighbour of the same size has been found
1562 if(candidateMorton == sameSizeVirtualNeighMorton && m_octants[candidateIdx].m_level == oct->m_level){
1563 isghost.push_back(false);
1564 neighbours.push_back(candidateIdx);
1565 return;
1566 }
1567
1568 // Compute the Morton number of the last candidate
1569 //
1570 // This is the Morton number of the last descendant of the same-size
1571 // virtual neighbour.
1572 uint8_t maxNodeNeighLevel = getMaxNodeNeighLevel(*oct);
1573 std::array<int64_t, 3> lastCandidateOffset = computeLastVirtualNodeNeighOffset(level, inode, maxNodeNeighLevel);
1574 std::array<uint32_t, 3> lastCandidateCoord = coord;
1575 lastCandidateCoord[0] = static_cast<uint32_t>(lastCandidateCoord[0] + lastCandidateOffset[0]);
1576 lastCandidateCoord[1] = static_cast<uint32_t>(lastCandidateCoord[1] + lastCandidateOffset[1]);
1577 lastCandidateCoord[2] = static_cast<uint32_t>(lastCandidateCoord[2] + lastCandidateOffset[2]);
1578 uint64_t lastCandidateMorton = PABLO::computeMorton(m_dim, lastCandidateCoord[0], lastCandidateCoord[1], lastCandidateCoord[2]);
1579
1580 // Search for neighbours of different sizes
1581 if (candidateIdx < getNumOctants()) {
1582 while(true){
1583 // Detect if the candidate is a neighbour
1584 u32array3 coordtry = {{0, 0, 0}};
1585 bool isNeighbour = true;
1586 for (int8_t idim=0; idim<m_dim; idim++){
1587 coordtry[idim] = m_octants[candidateIdx].getLogicalCoordinates(idim);
1588
1589 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1590 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_octants[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1591
1592 if (Dx != Dxstar) {
1593 isNeighbour = false;
1594 break;
1595 }
1596 }
1597
1598 // Advance to the next candidate
1599 if (isNeighbour){
1600 neighbours.push_back(candidateIdx);
1601 isghost.push_back(false);
1602 return;
1603 }
1604
1605 candidateIdx++;
1606 if (candidateIdx > getNumOctants()-1){
1607 break;
1608 }
1609
1610 candidateMorton = m_octants[candidateIdx].getMorton();
1611 if (candidateMorton > lastCandidateMorton){
1612 break;
1613 }
1614 }
1615 }
1616
1617 //
1618 // Search in the ghost octants
1619 //
1620
1621 if (getNumGhosts() > 0 && !onlyinternal){
1622 // Identify the index of the first neighbour candidate
1623 computeNeighSearchBegin(sameSizeVirtualNeighMorton, m_ghosts, &candidateIdx, &candidateMorton);
1624
1625 // Early return if a neighbour of the same size has been found
1626 if(candidateMorton == sameSizeVirtualNeighMorton && m_ghosts[candidateIdx].m_level == oct->m_level){
1627 isghost.push_back(true);
1628 neighbours.push_back(candidateIdx);
1629 return;
1630 }
1631
1632 // Search for neighbours of different sizes
1633 if (candidateIdx < getNumGhosts()) {
1634 while(true){
1635 // Detect if the candidate is a neighbour
1636 u32array3 coordtry = {{0, 0, 0}};
1637 bool isNeighbour = true;
1638 for (int8_t idim=0; idim<m_dim; idim++){
1639 coordtry[idim] = m_ghosts[candidateIdx].getLogicalCoordinates(idim);
1640
1641 int32_t Dx = int32_t(int32_t(abs(cxyz[idim]))*(-int32_t(coord[idim]) + int32_t(coordtry[idim])));
1642 int32_t Dxstar = int32_t((cxyz[idim]-1)/2)*(m_ghosts[candidateIdx].getLogicalSize()) + int32_t((cxyz[idim]+1)/2)*size;
1643
1644 if (Dx != Dxstar) {
1645 isNeighbour = false;
1646 break;
1647 }
1648 }
1649
1650 // Advance to the next candidate
1651 if (isNeighbour){
1652 neighbours.push_back(candidateIdx);
1653 isghost.push_back(true);
1654 return;
1655 }
1656
1657 candidateIdx++;
1658 if (candidateIdx > getNumGhosts()-1){
1659 break;
1660 }
1661
1662 candidateMorton = m_ghosts[candidateIdx].getMorton();
1663 if (candidateMorton > lastCandidateMorton){
1664 break;
1665 }
1666 }
1667 }
1668 }
1669 };
1670
1671 // =================================================================================== //
1688 void
1689 LocalTree::computeNeighSearchBegin(uint64_t sameSizeVirtualNeighMorton, const octvector &octants, uint32_t *searchBeginIdx, uint64_t *searchBeginMorton) const {
1690
1691 // Early return if there are no octants
1692 if (octants.empty()) {
1693 *searchBeginIdx = 0;
1694 *searchBeginMorton = PABLO::INVALID_MORTON;
1695 return;
1696 }
1697
1698 // The search should start from the lower bound if it points to the
1699 // first octant or to an octant whose Morton number is equal to the
1700 // the same-size virtual neighbour Morton number. Otherwise, the
1701 // search should start form the octant preceding the lower bound.
1702 uint32_t lowerBoundIdx;
1703 uint64_t lowerBoundMorton;
1704 findMortonLowerBound(sameSizeVirtualNeighMorton, octants, &lowerBoundIdx, &lowerBoundMorton);
1705
1706 if (lowerBoundMorton == sameSizeVirtualNeighMorton || lowerBoundIdx == 0) {
1707 *searchBeginIdx = lowerBoundIdx;
1708 *searchBeginMorton = lowerBoundMorton;
1709 } else {
1710 *searchBeginIdx = lowerBoundIdx - 1;
1711 *searchBeginMorton = octants[*searchBeginIdx].getMorton();
1712 }
1713
1714 }
1715
1716 // =================================================================================== //
1717
1725 bool
1726 LocalTree::fixBrokenFamiliesMarkers(std::vector<Octant *> *updatedOctants, std::vector<bool> *updatedGhostFlags){
1727
1728 // Initialization
1729 bool updated = false;
1730
1731 uint32_t nocts = m_octants.size();
1732 if (nocts == 0) {
1733 return updated;
1734 }
1735
1736 uint32_t nghosts = m_ghosts.size();
1737
1738 uint32_t idx = 0;
1739 uint32_t idx0 = 0;
1740 uint32_t idx2 = 0;
1741 uint32_t idx1_gh = 0;
1742 uint32_t idx2_gh = 0;
1743 uint32_t last_idx =nocts-1;
1744
1745 uint8_t nbro = 0;
1746 Octant father(m_dim);
1747
1748 //Clean index of ghost brothers in case of coarsening a broken family
1749 m_lastGhostBros.clear();
1750 m_firstGhostBros.clear();
1751
1752 // Process ghosts
1753 bool checkend = true;
1754 bool checkstart = true;
1755 if (nghosts > 0){
1756 // Set index for start and end check for ghosts
1757 while(m_ghosts[idx2_gh].getMorton() <= m_lastDescMorton){
1758 idx2_gh++;
1759 if (idx2_gh > nghosts-1) break;
1760 }
1761 if (idx2_gh > nghosts-1) checkend = false;
1762
1763 while(m_ghosts[idx1_gh].getMorton() <= m_octants[0].getMorton()){
1764 idx1_gh++;
1765 if (idx1_gh > nghosts-1) break;
1766 }
1767 if (idx1_gh == 0) checkstart = false;
1768 idx1_gh-=1;
1769
1770 // Start and End on ghosts
1771 if (checkstart){
1772 if (m_ghosts[idx1_gh].buildFather()==m_octants[0].buildFather()){
1773 father = m_ghosts[idx1_gh].buildFather();
1774 nbro = 0;
1775 idx = idx1_gh;
1776 int8_t marker = m_ghosts[idx].getMarker();
1777 while(marker < 0 && m_ghosts[idx].buildFather() == father){
1778
1779 //Add ghost index to structure for mapper in case of coarsening a broken family
1780 m_firstGhostBros.push_back(idx);
1781
1782 nbro++;
1783 if (idx==0)
1784 break;
1785 idx--;
1786 marker = m_ghosts[idx].getMarker();
1787 }
1788 idx = 0;
1789 while(idx<nocts && m_octants[idx].buildFather() == father){
1790 if (m_octants[idx].getMarker()<0)
1791 nbro++;
1792 idx++;
1793 if(idx==nocts)
1794 break;
1795 }
1796 if (nbro != m_treeConstants->nChildren && idx!=nocts){
1797 for(uint32_t ii=0; ii<idx; ii++){
1798 Octant &octant = m_octants[ii];
1799 if (octant.getMarker()<0){
1800 octant.setMarker(0);
1801 if (updatedOctants) {
1802 updatedOctants->push_back(&octant);
1803 }
1804 if (updatedGhostFlags) {
1805 updatedGhostFlags->push_back(false);
1806 }
1807 updated = true;
1808 }
1809 }
1810 //Clean index of ghost brothers in case of coarsening a broken family
1811 m_firstGhostBros.clear();
1812 }
1813 }
1814 }
1815
1816 if (checkend){
1817 bool checklocal = false;
1818 if (m_ghosts[idx2_gh].buildFather()==m_octants[nocts-1].buildFather()){
1819 father = m_ghosts[idx2_gh].buildFather();
1820 if (idx!=nocts){
1821 nbro = 0;
1822 checklocal = true;
1823 }
1824 idx = idx2_gh;
1825 int8_t marker = m_ghosts[idx].getMarker();
1826 while(marker < 0 && m_ghosts[idx].buildFather() == father){
1827
1828 //Add ghost index to structure for mapper in case of coarsening a broken family
1829 m_lastGhostBros.push_back(idx);
1830
1831 nbro++;
1832 idx++;
1833 if(idx == nghosts){
1834 break;
1835 }
1836 marker = m_ghosts[idx].getMarker();
1837 }
1838 idx = nocts-1;
1839 if (checklocal){
1840 while(m_octants[idx].buildFather() == father){
1841 if (m_octants[idx].getMarker()<0)
1842 nbro++;
1843 if (idx==0)
1844 break;
1845 idx--;
1846 }
1847 }
1848 last_idx=idx;
1849 if (nbro != m_treeConstants->nChildren && idx!=nocts-1){
1850 for(uint32_t ii=idx+1; ii<nocts; ii++){
1851 Octant &octant = m_octants[ii];
1852 if (octant.getMarker()<0){
1853 octant.setMarker(0);
1854 if (updatedOctants) {
1855 updatedOctants->push_back(&octant);
1856 }
1857 if (updatedGhostFlags) {
1858 updatedGhostFlags->push_back(false);
1859 }
1860 updated = true;
1861 }
1862 }
1863 //Clean index of ghost brothers in case of coarsening a broken family
1864 m_lastGhostBros.clear();
1865 }
1866 }
1867 }
1868 }
1869
1870 // Check first internal octants
1871 father = m_octants[0].buildFather();
1872 uint64_t mortonld = father.computeLastDescMorton();
1873 nbro = 0;
1874 for (idx=0; idx<m_treeConstants->nChildren; idx++){
1875 // Check if family is complete or to be checked in the internal loop (some brother refined)
1876 if (idx<nocts){
1877 if (m_octants[idx].getMorton() <= mortonld){
1878 nbro++;
1879 }
1880 }
1881 }
1882
1883 if (nbro != m_treeConstants->nChildren) {
1884 idx0 = nbro;
1885 }
1886
1887 // Check and coarse internal octants
1888 for (idx=idx0; idx<nocts; idx++){
1889 Octant &octant = m_octants[idx];
1890 if(octant.getMarker() < 0 && octant.getLevel() > 0){
1891 nbro = 0;
1892 father = octant.buildFather();
1893 // Check if family is to be coarsened
1894 for (idx2=idx; idx2<idx+m_treeConstants->nChildren; idx2++){
1895 if (idx2<nocts){
1896 if(m_octants[idx2].getMarker() < 0 && m_octants[idx2].buildFather() == father){
1897 nbro++;
1898 }
1899 }
1900 }
1901 if (nbro == m_treeConstants->nChildren){
1902 idx = idx2-1;
1903 }
1904 else{
1905 if (idx<=last_idx){
1906 octant.setMarker(0);
1907 if (updatedOctants) {
1908 updatedOctants->push_back(&octant);
1909 }
1910 if (updatedGhostFlags) {
1911 updatedGhostFlags->push_back(false);
1912 }
1913 updated = true;
1914 }
1915 }
1916 }
1917 }
1918
1919 return updated;
1920 }
1921
1922 // =================================================================================== //
1923
1930 bool
1931 LocalTree::localBalance(bool doNew, bool checkInterior, bool checkGhost){
1932
1933 bool balanceEdges = ((m_balanceCodim>1) && (m_dim==3));
1934 bool balanceNodes = (m_balanceCodim==m_dim);
1935
1936 std::vector<Octant *> processOctants;
1937 std::vector<bool> processGhostFlags;
1938
1939 // Identify internal octants that will be processed
1940 if(checkInterior){
1941 for (Octant &octant : m_octants){
1942 // Skip octants that doesn't need balancing
1943 bool balanceOctant = octant.getBalance();
1944 if (balanceOctant) {
1945 if (doNew) {
1946 balanceOctant = (octant.getMarker() != 0) || octant.getIsNewC() || octant.getIsNewR();
1947 } else {
1948 balanceOctant = (octant.getMarker() != 0);
1949 }
1950 }
1951
1952 if (!balanceOctant) {
1953 continue;
1954 }
1955
1956 // Add octant to the process list
1957 processOctants.push_back(&octant);
1958 processGhostFlags.push_back(false);
1959 }
1960 }
1961
1962 // Identify ghost octants that will be processed
1963 //
1964 // Ghost octants will be balanced by the process that owns them, however they may
1965 // affect balancing of local octants. If ghost octnts are processed, we process all
1966 // ghost octants of the first layer, not only the ones that need balancing. That's
1967 // because it's faster to process all the ghost octants that may affect balacning
1968 // of internal octants, rather than find the ones that actually affect balacing of
1969 // internal octants.
1970 if (checkGhost) {
1971 for (Octant &octant : m_ghosts){
1972 // Only ghosts of the first layer can affect load balance
1973 if (octant.getGhostLayer() > 0) {
1974 continue;
1975 }
1976
1977 // Add octant to the process list
1978 processOctants.push_back(&octant);
1979 processGhostFlags.push_back(true);
1980 }
1981 }
1982
1983 // Iterative balacing
1984 std::vector<uint32_t> neighs;
1985 std::vector<bool> neighGhostFlags;
1986
1987 bool updated = false;
1988 std::vector<Octant *> updatedProcessOctants;
1989 std::vector<bool> updatedProcessGhostFlags;
1990 while (!processOctants.empty()) {
1991 // Fix broken families markers
1992 bool brokenFamilieisFixed = fixBrokenFamiliesMarkers(&processOctants, &processGhostFlags);
1993 if (brokenFamilieisFixed) {
1994 updated = true;
1995 }
1996
1997 // Balance
1998 std::size_t processSize = processOctants.size();
1999 for (std::size_t n = 0; n < processSize; ++n) {
2000 Octant &octant = *(processOctants[n]);
2001 bool ghostFlag = processGhostFlags[n];
2002
2003 // Identify neighbours that will be used for balancing
2004 neighs.clear();
2005 neighGhostFlags.clear();
2006
2007 for (int iface=0; iface < m_treeConstants->nFaces; iface++) {
2008 findNeighbours(&octant, iface, neighs, neighGhostFlags, ghostFlag, true);
2009 }
2010
2011 if (balanceNodes) {
2012 for (int inode=0; inode < m_treeConstants->nNodes; inode++) {
2013 findNodeNeighbours(&octant, inode, neighs, neighGhostFlags, ghostFlag, true);
2014 }
2015 }
2016
2017 if (balanceEdges) {
2018 for (int iedge=0; iedge < m_treeConstants->nEdges; iedge++) {
2019 findEdgeNeighbours(&octant, iedge, neighs, neighGhostFlags, ghostFlag, true);
2020 }
2021 }
2022
2023 // Balance octant
2024 int8_t level = octant.getLevel();
2025 int8_t futureLevel = std::min(m_treeConstants->maxLevel, static_cast<int8_t>(level + octant.getMarker()));
2026
2027 std::size_t nNeighs = neighs.size();
2028 for(std::size_t i = 0; i < nNeighs; i++){
2029 bool neighGhostFlag = neighGhostFlags[i];
2030 Octant &neighOctant = (neighGhostFlag ? m_ghosts[neighs[i]]: m_octants[neighs[i]]);
2031
2032 int8_t neighLevel = neighOctant.getLevel();
2033 int8_t neighFutureLevel = std::min(m_treeConstants->maxLevel, int8_t(neighLevel + neighOctant.getMarker()));
2034 if(!ghostFlag && futureLevel < neighFutureLevel - 1) {
2035 futureLevel = neighFutureLevel - 1;
2036 octant.setMarker(futureLevel - level);
2037 updatedProcessOctants.push_back(&octant);
2038 updatedProcessGhostFlags.push_back(ghostFlag);
2039 updated = true;
2040 }
2041 else if(!neighGhostFlag && neighFutureLevel < futureLevel - 1) {
2042 neighOctant.setMarker((futureLevel - 1) - neighLevel);
2043 if (neighOctant.getBalance()) {
2044 updatedProcessOctants.push_back(&neighOctant);
2045 updatedProcessGhostFlags.push_back(neighGhostFlag);
2046 }
2047 updated = true;
2048 }
2049 }
2050 }
2051
2052 // Update process list
2053 updatedProcessOctants.swap(processOctants);
2054 updatedProcessGhostFlags.swap(processGhostFlags);
2055
2056 updatedProcessOctants.clear();
2057 updatedProcessGhostFlags.clear();
2058 }
2059
2060 return updated;
2061 };
2062
2063 // =================================================================================== //
2064
2073 std::array<int64_t, 3> LocalTree::computeFirstVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const {
2074
2075 // Get normalized offsets
2076 const int8_t (&normalizedOffsets)[3] = m_treeConstants->normals[iface];
2077
2078 // Get octant sizes
2079 uint32_t size = m_treeConstants->lengths[level];
2080 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2081
2082 // Compute the coordinates of the virtual neighbour
2083 std::array<int64_t, 3> neighOffsets;
2084 for (int i = 0; i < 3; ++i) {
2085 neighOffsets[i] = normalizedOffsets[i];
2086 if (neighOffsets[i] > 0) {
2087 neighOffsets[i] *= size;
2088 } else {
2089 neighOffsets[i] *= neighSize;
2090 }
2091 }
2092
2093 return neighOffsets;
2094 }
2095
2104 std::array<int64_t, 3> LocalTree::computeLastVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const {
2105
2106 // Get octant sizes
2107 uint32_t size = m_treeConstants->lengths[level];
2108 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2109
2110 // Compute the coordinates of the virtual neighbour
2111 std::array<int64_t, 3> neighOffsets = computeFirstVirtualNeighOffset(level, iface, neighLevel);
2112
2113 uint32_t delta = size - neighSize;
2114 switch (iface) {
2115 case 0 :
2116 case 1 :
2117 neighOffsets[1] += delta;
2118 neighOffsets[2] += delta;
2119 break;
2120
2121 case 2 :
2122 case 3 :
2123 neighOffsets[0] += delta;
2124 neighOffsets[2] += delta;
2125 break;
2126
2127 case 4 :
2128 case 5 :
2129 neighOffsets[0] += delta;
2130 neighOffsets[1] += delta;
2131 break;
2132 }
2133
2134 return neighOffsets;
2135 }
2136
2145 void LocalTree::computeVirtualNeighOffsets(uint8_t level, uint8_t iface, uint8_t neighLevel, std::vector<std::array<int64_t, 3>> *neighOffsets) const {
2146
2147 // Get octant sizes
2148 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2149
2150 // Direction of the increments
2151 int direction1;
2152 int direction2;
2153 switch (iface) {
2154 case 0:
2155 case 1:
2156 direction1 = 1;
2157 direction2 = 2;
2158 break;
2159
2160 case 2:
2161 case 3:
2162 direction1 = 0;
2163 direction2 = 2;
2164 break;
2165
2166 default:
2167 direction1 = 0;
2168 direction2 = 1;
2169 break;
2170 }
2171
2172 // Compute the coordinates of the virtual neighbour
2173 int nNeighsDirection1 = 1 << (neighLevel - level);
2174 int nNeighsDirection2 = (m_dim > 2) ? nNeighsDirection1 : 1;
2175 int nNeighs = nNeighsDirection1 * nNeighsDirection2;
2176
2177 neighOffsets->assign(nNeighs, computeFirstVirtualNeighOffset(level, iface, neighLevel));
2178
2179 int k = 0;
2180 for (int j = 0; j < nNeighsDirection2; ++j) {
2181 for (int i = 0; i < nNeighsDirection1; ++i) {
2182 (*neighOffsets)[k][direction1] += i * neighSize;
2183 (*neighOffsets)[k][direction2] += j * neighSize;
2184 ++k;
2185 }
2186 }
2187 }
2188
2197 std::array<int64_t, 3> LocalTree::computeFirstVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const {
2198
2199 // Get normalized offsets
2200 const int8_t (&normalizedOffsets)[3] = m_treeConstants->nodeCoeffs[inode];
2201
2202 // Get octant sizes
2203 uint32_t size = m_treeConstants->lengths[level];
2204 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2205
2206 // Compute the coordinates of the virtual neighbour
2207 std::array<int64_t, 3> neighOffsets;
2208 for (int i = 0; i < 3; ++i) {
2209 neighOffsets[i] = normalizedOffsets[i];
2210 if (neighOffsets[i] > 0) {
2211 neighOffsets[i] *= size;
2212 } else {
2213 neighOffsets[i] *= neighSize;
2214 }
2215 }
2216
2217 return neighOffsets;
2218 }
2219
2228 std::array<int64_t, 3> LocalTree::computeLastVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const {
2229
2230 return computeFirstVirtualNodeNeighOffset(level, inode, neighLevel);
2231 }
2232
2241 void LocalTree::computeVirtualNodeNeighOffsets(uint8_t level, uint8_t inode, uint8_t neighLevel, std::vector<std::array<int64_t, 3>> *neighOffsets) const {
2242
2243 neighOffsets->assign(1, computeFirstVirtualNodeNeighOffset(level, inode, neighLevel));
2244 }
2245
2254 std::array<int64_t, 3> LocalTree::computeFirstVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const {
2255
2256 // Get normalized offsets
2257 const int8_t (&normalizedOffsets)[3] = m_treeConstants->edgeCoeffs[iedge];
2258
2259 // Get octant sizes
2260 uint32_t size = m_treeConstants->lengths[level];
2261 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2262
2263 // Compute the coordinates of the virtual neighbour
2264 std::array<int64_t, 3> neighOffsets;
2265 for (int i = 0; i < 3; ++i) {
2266 neighOffsets[i] = normalizedOffsets[i];
2267 if (neighOffsets[i] > 0) {
2268 neighOffsets[i] *= size;
2269 } else {
2270 neighOffsets[i] *= neighSize;
2271 }
2272 }
2273
2274 return neighOffsets;
2275 }
2276
2285 std::array<int64_t, 3> LocalTree::computeLastVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const {
2286
2287 // Get octant sizes
2288 uint32_t size = m_treeConstants->lengths[level];
2289 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2290
2291 // Compute the coordinates of the virtual neighbour
2292 std::array<int64_t, 3> neighOffsets = computeFirstVirtualEdgeNeighOffset(level, iedge, neighLevel);
2293
2294 uint32_t offset = size - neighSize;
2295 switch (iedge) {
2296 case 0:
2297 case 1:
2298 case 8:
2299 case 9:
2300 neighOffsets[1] += offset;
2301 break;
2302
2303 case 2:
2304 case 3:
2305 case 10:
2306 case 11:
2307 neighOffsets[0] += offset;
2308 break;
2309
2310 case 4:
2311 case 5:
2312 case 6:
2313 case 7:
2314 neighOffsets[2] += offset;
2315 break;
2316 }
2317
2318 return neighOffsets;
2319 }
2320
2329 void LocalTree::computeVirtualEdgeNeighOffsets(uint8_t level, uint8_t iedge, uint8_t neighLevel, std::vector<std::array<int64_t, 3>> *neighOffsets) const {
2330
2331 // Get octant sizes
2332 uint32_t neighSize = m_treeConstants->lengths[neighLevel];
2333
2334 // Direction of the increments
2335 int direction;
2336 switch (iedge) {
2337 case 0:
2338 case 1:
2339 case 8:
2340 case 9:
2341 direction = 1;
2342 break;
2343
2344 case 2:
2345 case 3:
2346 case 10:
2347 case 11:
2348 direction = 0;
2349 break;
2350
2351 default:
2352 direction = 2;
2353 break;
2354 }
2355
2356 // Compute the coordinates of the virtual neighbours
2357 int nNeighs = 1 << (neighLevel - level);
2358 neighOffsets->assign(nNeighs, computeFirstVirtualEdgeNeighOffset(level, iedge, neighLevel));
2359 for (int i = 1; i < nNeighs; ++i) {
2360 (*neighOffsets)[i][direction] += i * neighSize;
2361 }
2362 }
2363
2364 // =================================================================================== //
2365
2373 std::array<int64_t, 3> LocalTree::getPeriodicOffset(const Octant &octant, uint8_t iface) const {
2374
2375 std::array<int64_t, 3> offset = {{0, 0, 0}};
2376 if (!isPeriodic(&octant, iface)) {
2377 return offset;
2378 }
2379
2380 int64_t maxLength = m_treeConstants->lengths[0];
2381 const int8_t (&referenceOffset)[3] = m_treeConstants->normals[iface];
2382
2383 for (int d = 0; d < m_dim; ++d) {
2384 offset[d] = - maxLength * static_cast<int64_t>(referenceOffset[d]);
2385 }
2386
2387 return offset;
2388 }
2389
2397 std::array<int64_t, 3> LocalTree::getNodePeriodicOffset(const Octant &octant, uint8_t inode) const {
2398
2399 std::array<int64_t, 3> offset = {{0, 0, 0}};
2400 for (int i = 0; i < m_dim; ++i) {
2401 uint8_t face = m_treeConstants->nodeFace[inode][i];
2402 if (!isPeriodic(&octant, face)) {
2403 continue;
2404 }
2405
2406 std::array<int64_t, 3> faceOffset = getPeriodicOffset(octant, face);
2407 for (int d = 0; d < m_dim; ++d) {
2408 offset[d] += faceOffset[d];
2409 }
2410 }
2411
2412 return offset;
2413 }
2414
2422 std::array<int64_t, 3> LocalTree::getEdgePeriodicOffset(const Octant &octant, uint8_t iedge) const {
2423
2424 std::array<int64_t, 3> offset = {{0, 0, 0}};
2425 for (uint8_t face : m_treeConstants->edgeFace[iedge]) {
2426 if (!isPeriodic(&octant, face)) {
2427 continue;
2428 }
2429
2430 std::array<int64_t, 3> faceOffset = getPeriodicOffset(octant, face);
2431 for (int d = 0; d < m_dim; ++d) {
2432 offset[d] += faceOffset[d];
2433 }
2434 }
2435
2436 return offset;
2437 }
2438
2439 // =================================================================================== //
2440
2447 uint8_t LocalTree::getMaxNeighLevel(const Octant &octant) const {
2448
2449 if (octant.getBalance()) {
2450 return octant.getLevel() + 1;
2451 } else {
2452 return m_treeConstants->maxLevel;
2453 }
2454 }
2455
2462 uint8_t LocalTree::getMaxNodeNeighLevel(const Octant &octant) const {
2463
2464 if (octant.getBalance() && m_balanceCodim == m_dim) {
2465 return octant.getLevel() + 1;
2466 } else {
2467 return m_treeConstants->maxLevel;
2468 }
2469 }
2470
2477 uint8_t LocalTree::getMaxEdgeNeighLevel(const Octant &octant) const {
2478
2479 if (octant.getBalance() && m_balanceCodim > 1) {
2480 return octant.getLevel() + 1;
2481 } else {
2482 return m_treeConstants->maxLevel;
2483 }
2484 }
2485
2486 // =================================================================================== //
2487
2490 void
2491 LocalTree::computeIntersections() {
2492
2493 octvector::iterator it, obegin, oend;
2494 Intersection intersection;
2495 u32vector neighbours;
2496 vector<bool> isghost;
2497 uint32_t idx;
2498 uint32_t i, nsize;
2499 uint8_t iface, iface2;
2500
2501 m_intersections.clear();
2502 m_intersections.reserve(2*3*m_octants.size());
2503
2504 idx = 0;
2505
2506 // Loop on ghosts
2507 obegin = m_ghosts.begin();
2508 oend = m_ghosts.end();
2509 for (it = obegin; it != oend; ++it){
2510 for (iface = 0; iface < m_dim; iface++){
2511 iface2 = iface*2;
2512 findNeighbours(m_ghosts.data() + idx, iface2, neighbours, isghost, true, false);
2513 nsize = neighbours.size();
2514 if (!(it->m_info[iface2])){
2515 //Internal intersection
2516 for (i = 0; i < nsize; i++){
2517 intersection.m_dim = m_dim;
2518 intersection.m_finer = getGhostLevel(idx) >= getLevel((int)neighbours[i]);
2519 intersection.m_out = intersection.m_finer;
2520 intersection.m_outisghost = intersection.m_finer;
2521 intersection.m_owners[0] = neighbours[i];
2522 intersection.m_owners[1] = idx;
2523 intersection.m_iface = m_treeConstants->oppositeFace[iface2] - (getGhostLevel(idx) >= getLevel((int)neighbours[i]));
2524 intersection.m_isnew = false;
2525 intersection.m_isghost = true;
2526 intersection.m_bound = false;
2527 intersection.m_pbound = true;
2528 m_intersections.push_back(intersection);
2529 }
2530 }
2531 else{
2532 //Periodic intersection
2533 for (i = 0; i < nsize; i++){
2534 intersection.m_dim = m_dim;
2535 intersection.m_finer = getGhostLevel(idx) >= getLevel((int)neighbours[i]);
2536 intersection.m_out = intersection.m_finer;
2537 intersection.m_outisghost = intersection.m_finer;
2538 intersection.m_owners[0] = neighbours[i];
2539 intersection.m_owners[1] = idx;
2540 intersection.m_iface = m_treeConstants->oppositeFace[iface2] - (getGhostLevel(idx) >= getLevel((int)neighbours[i]));
2541 intersection.m_isnew = false;
2542 intersection.m_isghost = true;
2543 intersection.m_bound = true;
2544 intersection.m_pbound = true;
2545 m_intersections.push_back(intersection);
2546 }
2547 }
2548 }
2549 idx++;
2550 }
2551
2552 // Loop on octants
2553 idx=0;
2554 obegin = m_octants.begin();
2555 oend = m_octants.end();
2556 for (it = obegin; it != oend; ++it){
2557 for (iface = 0; iface < m_dim; iface++){
2558 iface2 = iface*2;
2559 findNeighbours(m_octants.data() + idx, iface2, neighbours, isghost, false, false);
2560 nsize = neighbours.size();
2561 if (nsize) {
2562 if (!(it->m_info[iface2])){
2563 //Internal intersection
2564 for (i = 0; i < nsize; i++){
2565 if (isghost[i]){
2566 intersection.m_dim = m_dim;
2567 intersection.m_owners[0] = idx;
2568 intersection.m_owners[1] = neighbours[i];
2569 intersection.m_finer = (nsize>1);
2570 intersection.m_out = (nsize>1);
2571 intersection.m_outisghost = (nsize>1);
2572 intersection.m_iface = iface2 + (nsize>1);
2573 intersection.m_isnew = false;
2574 intersection.m_isghost = true;
2575 intersection.m_bound = false;
2576 intersection.m_pbound = true;
2577 m_intersections.push_back(intersection);
2578 }
2579 else{
2580 intersection.m_dim = m_dim;
2581 intersection.m_owners[0] = idx;
2582 intersection.m_owners[1] = neighbours[i];
2583 intersection.m_finer = (nsize>1);
2584 intersection.m_out = (nsize>1);
2585 intersection.m_outisghost = false;
2586 intersection.m_iface = iface2 + (nsize>1);
2587 intersection.m_isnew = false;
2588 intersection.m_isghost = false;
2589 intersection.m_bound = false;
2590 intersection.m_pbound = false;
2591 m_intersections.push_back(intersection);
2592 }
2593 }
2594 }
2595 else{
2596 //Periodic intersection
2597 for (i = 0; i < nsize; i++){
2598 if (isghost[i]){
2599 intersection.m_dim = m_dim;
2600 intersection.m_owners[0] = idx;
2601 intersection.m_owners[1] = neighbours[i];
2602 intersection.m_finer = (nsize>1);
2603 intersection.m_out = intersection.m_finer;
2604 intersection.m_outisghost = intersection.m_finer;
2605 intersection.m_iface = iface2 + (nsize>1);
2606 intersection.m_isnew = false;
2607 intersection.m_isghost = true;
2608 intersection.m_bound = true;
2609 intersection.m_pbound = true;
2610 m_intersections.push_back(intersection);
2611 }
2612 else{
2613 intersection.m_dim = m_dim;
2614 intersection.m_owners[0] = idx;
2615 intersection.m_owners[1] = neighbours[i];
2616 intersection.m_finer = (nsize>1);
2617 intersection.m_out = intersection.m_finer;
2618 intersection.m_outisghost = false;
2619 intersection.m_iface = iface2 + (nsize>1);
2620 intersection.m_isnew = false;
2621 intersection.m_isghost = false;
2622 intersection.m_bound = true;
2623 intersection.m_pbound = false;
2624 m_intersections.push_back(intersection);
2625 }
2626 }
2627 }
2628 }
2629 else{
2630 //Boundary intersection
2631 intersection.m_dim = m_dim;
2632 intersection.m_owners[0] = idx;
2633 intersection.m_owners[1] = idx;
2634 intersection.m_finer = 0;
2635 intersection.m_out = 0;
2636 intersection.m_outisghost = false;
2637 intersection.m_iface = iface2;
2638 intersection.m_isnew = false;
2639 intersection.m_isghost = false;
2640 intersection.m_bound = true;
2641 intersection.m_pbound = false;
2642 m_intersections.push_back(intersection);
2643 }
2644 if (it->m_info[iface2+1]){
2645 if (!(m_periodic[iface2+1])){
2646 //Boundary intersection
2647 intersection.m_dim = m_dim;
2648 intersection.m_owners[0] = idx;
2649 intersection.m_owners[1] = idx;
2650 intersection.m_finer = 0;
2651 intersection.m_out = 0;
2652 intersection.m_outisghost = false;
2653 intersection.m_iface = iface2+1;
2654 intersection.m_isnew = false;
2655 intersection.m_isghost = false;
2656 intersection.m_bound = true;
2657 intersection.m_pbound = false;
2658 m_intersections.push_back(intersection);
2659 }
2660 else{
2661 //Periodic intersection
2662 findNeighbours(m_octants.data() + idx, iface2+1, neighbours, isghost, false, false);
2663 nsize = neighbours.size();
2664 for (i = 0; i < nsize; i++){
2665 if (isghost[i]){
2666 intersection.m_dim = m_dim;
2667 intersection.m_owners[0] = idx;
2668 intersection.m_owners[1] = neighbours[i];
2669 intersection.m_finer = (nsize>1);
2670 intersection.m_out = intersection.m_finer;
2671 intersection.m_outisghost = intersection.m_finer;
2672 intersection.m_iface = iface2 + (nsize>1);
2673 intersection.m_isnew = false;
2674 intersection.m_isghost = true;
2675 intersection.m_bound = true;
2676 intersection.m_pbound = true;
2677 m_intersections.push_back(intersection);
2678 }
2679 else{
2680 intersection.m_dim = m_dim;
2681 intersection.m_owners[0] = idx;
2682 intersection.m_owners[1] = neighbours[i];
2683 intersection.m_finer = (nsize>1);
2684 intersection.m_out = intersection.m_finer;
2685 intersection.m_outisghost = false;
2686 intersection.m_iface = iface2 + (nsize>1);
2687 intersection.m_isnew = false;
2688 intersection.m_isghost = false;
2689 intersection.m_bound = true;
2690 intersection.m_pbound = false;
2691 m_intersections.push_back(intersection);
2692 }
2693 }
2694 }
2695 }
2696 }
2697 idx++;
2698 }
2699 intervector(m_intersections).swap(m_intersections);
2700 }
2701
2702 // =================================================================================== //
2707 uint32_t
2708 LocalTree::findMorton(uint64_t targetMorton) const {
2709 return findMorton(targetMorton, m_octants);
2710 };
2711
2712 // =================================================================================== //
2717 uint32_t
2718 LocalTree::findGhostMorton(uint64_t targetMorton) const {
2719 return findMorton(targetMorton, m_ghosts);
2720 };
2721
2722 // =================================================================================== //
2731 uint32_t
2732 LocalTree::findMorton(uint64_t targetMorton, const octvector &octants) const {
2733
2734 uint32_t lowerBoundIdx;
2735 uint64_t lowerBoundMorton;
2736 findMortonLowerBound(targetMorton, octants, &lowerBoundIdx, &lowerBoundMorton);
2737
2738 uint32_t targetIdx;
2739 if (lowerBoundMorton == targetMorton) {
2740 targetIdx = lowerBoundIdx;
2741 } else {
2742 targetIdx = octants.size();
2743 }
2744
2745 return targetIdx;
2746 };
2747
2748 // =================================================================================== //
2766 void
2767 LocalTree::findMortonLowerBound(uint64_t targetMorton, const octvector &octants, uint32_t *lowerBoundIdx, uint64_t *lowerBoundMorton) const {
2768
2769 uint32_t nOctants = octants.size();
2770
2771 uint32_t lowIndex = 0;
2772 uint32_t highIndex = nOctants;
2773 uint32_t midIndex = nOctants;
2774 uint64_t midMorton = PABLO::INVALID_MORTON;
2775 while (lowIndex < highIndex) {
2776 midIndex = lowIndex + (highIndex - lowIndex) / 2;
2777 midMorton = octants[midIndex].getMorton();
2778 if (targetMorton < midMorton) {
2779 highIndex = midIndex;
2780 }
2781 else if (targetMorton > midMorton) {
2782 lowIndex = midIndex + 1;
2783 }
2784 else {
2785 *lowerBoundIdx = midIndex;
2786 *lowerBoundMorton = midMorton;
2787
2788 return;
2789 }
2790 }
2791
2792 *lowerBoundIdx = lowIndex;
2793 if (*lowerBoundIdx == midIndex) {
2794 *lowerBoundMorton = midMorton;
2795 }
2796 else if (*lowerBoundIdx < nOctants) {
2797 *lowerBoundMorton = octants[*lowerBoundIdx].getMorton();
2798 }
2799 else {
2800 *lowerBoundMorton = PABLO::INVALID_MORTON;
2801 }
2802
2803 }
2804
2805 // =================================================================================== //
2821 void
2822 LocalTree::findMortonUpperBound(uint64_t targetMorton, const octvector &octants, uint32_t *upperBoundIdx, uint64_t *upperBoundMorton) const {
2823
2824 uint32_t nOctants = octants.size();
2825
2826 uint32_t lowIndex = 0;
2827 uint32_t highIndex = nOctants;
2828 uint32_t midIndex = nOctants;
2829 uint64_t midMorton = PABLO::INVALID_MORTON;
2830 while (lowIndex < highIndex) {
2831 midIndex = lowIndex + (highIndex - lowIndex) / 2;
2832 midMorton = octants[midIndex].getMorton();
2833 if (targetMorton < midMorton) {
2834 highIndex = midIndex;
2835 }
2836 else {
2837 lowIndex = midIndex + 1;
2838 }
2839 }
2840
2841 *upperBoundIdx = lowIndex;
2842 if (*upperBoundIdx == midIndex) {
2843 *upperBoundMorton = midMorton;
2844 }
2845 else if (*upperBoundIdx < nOctants) {
2846 *upperBoundMorton = octants[*upperBoundIdx].getMorton();
2847 }
2848 else {
2849 *upperBoundMorton = PABLO::INVALID_MORTON;
2850 }
2851
2852 }
2853
2854 // =================================================================================== //
2855
2858 void
2859 LocalTree::computeConnectivity(){
2860 vector<uint64_t> nodeKeys;
2861 unordered_map<uint64_t, array<uint32_t, 3> > nodeCoords;
2862 unordered_map<uint64_t, vector<uint64_t> > nodeOctants;
2863 uint32_t noctants = getNumOctants();
2864 uint32_t nghosts = getNumGhosts();
2865
2866
2867
2868 // Gather node information
2869 nodeKeys.reserve(noctants);
2870 nodeCoords.reserve(noctants);
2871 nodeOctants.reserve(noctants);
2872
2873 for (uint64_t n = 0; n < (noctants + nghosts); n++){
2874 const Octant *octant;
2875 if (n < noctants) {
2876 uint32_t octantId = static_cast<uint32_t>(n);
2877 octant = &(m_octants[octantId]);
2878 } else {
2879 uint32_t octantId = static_cast<uint32_t>(n - noctants);
2880 octant = &(m_ghosts[octantId]);
2881 }
2882
2883 for (uint8_t i = 0; i < m_treeConstants->nNodes; ++i){
2884 u32array3 node;
2885 octant->getLogicalNode(node, i);
2886
2887 uint64_t nodeKey = octant->computeNodePersistentKey(node);
2888 if (nodeCoords.count(nodeKey) == 0) {
2889 nodeKeys.push_back(nodeKey);
2890 nodeCoords.insert({{nodeKey, std::move(node)}});
2891 nodeOctants[nodeKey].reserve(m_treeConstants->nNodes);
2892 }
2893
2894 nodeOctants[nodeKey].push_back(n);
2895 }
2896 }
2897 std::sort(nodeKeys.begin(), nodeKeys.end());
2898
2899 // Build node list and connectivity
2900 m_nodes.clear();
2901 m_connectivity.clear();
2902 m_ghostsConnectivity.clear();
2903
2904 m_nodes.reserve(nodeKeys.size());
2905 m_connectivity.resize(noctants);
2906 m_ghostsConnectivity.resize(nghosts);
2907
2908 uint32_t nodeId = 0;
2909 for (uint64_t nodeKey : nodeKeys) {
2910 m_nodes.emplace_back(std::move(nodeCoords.at(nodeKey)));
2911 for (uint64_t n : nodeOctants.at(nodeKey)) {
2912 std::vector<uint32_t> *octantConnect;
2913 if (n < noctants) {
2914 uint32_t octantId = static_cast<uint32_t>(n);
2915 octantConnect = &(m_connectivity[octantId]);
2916 } else {
2917 uint32_t octantId = static_cast<uint32_t>(n - noctants);
2918 octantConnect = &(m_ghostsConnectivity[octantId]);
2919 }
2920
2921 if (octantConnect->size() == 0) {
2922 octantConnect->reserve(m_treeConstants->nNodes);
2923 }
2924 octantConnect->push_back(nodeId);
2925 }
2926 nodeId++;
2927 }
2928
2929 m_nodes.shrink_to_fit();
2930 };
2931
2937 void
2938 LocalTree::clearConnectivity(bool release){
2939 if (release) {
2940 u32arr3vector().swap(m_nodes);
2941 u32vector2D().swap(m_connectivity);
2942 u32vector2D().swap(m_ghostsConnectivity);
2943 } else {
2944 m_nodes.clear();
2945 m_connectivity.clear();
2946 m_ghostsConnectivity.clear();
2947 }
2948 };
2949
2952 void
2953 LocalTree::updateConnectivity(){
2954 clearConnectivity(false);
2955 computeConnectivity();
2956 };
2957
2958 // =================================================================================== //
2959
2960}
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:2241
bitpit::LocalTree::getMaxEdgeNeighLevel
uint8_t getMaxEdgeNeighLevel(const Octant &octant) const
Definition LocalTree.cpp:2477
bitpit::LocalTree::computeLastVirtualEdgeNeighOffset
std::array< int64_t, 3 > computeLastVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const
Definition LocalTree.cpp:2285
bitpit::LocalTree::intervector
std::vector< Intersection > intervector
Definition LocalTree.hpp:90
bitpit::LocalTree::getMaxNeighLevel
uint8_t getMaxNeighLevel(const Octant &octant) const
Definition LocalTree.cpp:2447
bitpit::LocalTree::getMaxNodeNeighLevel
uint8_t getMaxNodeNeighLevel(const Octant &octant) const
Definition LocalTree.cpp:2462
bitpit::LocalTree::getEdgePeriodicOffset
std::array< int64_t, 3 > getEdgePeriodicOffset(const Octant &octant, uint8_t iedge) const
Definition LocalTree.cpp:2422
bitpit::LocalTree::getPeriodicOffset
std::array< int64_t, 3 > getPeriodicOffset(const Octant &octant, uint8_t iface) const
Definition LocalTree.cpp:2373
bitpit::LocalTree::computeLastVirtualNeighOffset
std::array< int64_t, 3 > computeLastVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const
Definition LocalTree.cpp:2104
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:2329
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:2228
bitpit::LocalTree::computeFirstVirtualNeighOffset
std::array< int64_t, 3 > computeFirstVirtualNeighOffset(uint8_t level, uint8_t iface, uint8_t neighLevel) const
Definition LocalTree.cpp:2073
bitpit::LocalTree::computeFirstVirtualEdgeNeighOffset
std::array< int64_t, 3 > computeFirstVirtualEdgeNeighOffset(uint8_t level, uint8_t iedge, uint8_t neighLevel) const
Definition LocalTree.cpp:2254
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:2145
bitpit::LocalTree::computeFirstVirtualNodeNeighOffset
std::array< int64_t, 3 > computeFirstVirtualNodeNeighOffset(uint8_t level, uint8_t inode, uint8_t neighLevel) const
Definition LocalTree.cpp:2197
bitpit::LocalTree::getNodePeriodicOffset
std::array< int64_t, 3 > getNodePeriodicOffset(const Octant &octant, uint8_t inode) const
Definition LocalTree.cpp:2397
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::instance
static BITPIT_PUBLIC_API const TreeConstants & instance(uint8_t dim)
Definition tree_constants.cpp:41
Generated on Mon Apr 28 2025 20:53:28 for bitpit by   doxygen 1.13.2