Loading...
Searching...
No Matches
ParaTree.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 "bitpit_operators.hpp"
29
30#include "morton.hpp"
31#include "ParaTree.hpp"
32
33#include <climits>
34#include <sstream>
35#include <iomanip>
36#include <fstream>
37#include <iterator>
38
39namespace bitpit {
40
41 // =================================================================================== //
42 // NAME SPACES //
43 // =================================================================================== //
44 using namespace std;
45
46 // =================================================================================== //
47 // AUXILIARY IMPLEMENTATIONS //
48 // =================================================================================== //
49
84 ParaTree::LoadBalanceRanges::LoadBalanceRanges()
85 : sendAction(ACTION_UNDEFINED), recvAction(ACTION_UNDEFINED)
86 {
87 }
88
96 ParaTree::LoadBalanceRanges::LoadBalanceRanges(bool serial, const ExchangeRanges &_sendRanges, const ExchangeRanges &_recvRanges)
97 : sendRanges(_sendRanges), recvRanges(_recvRanges)
98 {
99 if (serial) {
100 sendAction = ACTION_DELETE;
101 recvAction = ACTION_NONE;
102 } else {
103 sendAction = ACTION_SEND;
104 recvAction = ACTION_RECEIVE;
105 }
106 }
107
110 void ParaTree::LoadBalanceRanges::clear()
111 {
112 sendAction = ACTION_UNDEFINED;
113 sendRanges.clear();
114
115 recvAction = ACTION_UNDEFINED;
116 recvRanges.clear();
117 }
118
119 // =================================================================================== //
120 // CLASS IMPLEMENTATION //
121 // =================================================================================== //
122
123 const std::string ParaTree::DEFAULT_LOG_FILE = "PABLO";
124
125 // =================================================================================== //
126 // CONSTRUCTORS AND OPERATORS //
127 // =================================================================================== //
128
132#if BITPIT_ENABLE_MPI==1
136 ParaTree::ParaTree(const std::string &logfile, MPI_Comm comm )
137#else
138 ParaTree::ParaTree(const std::string &logfile )
139#endif
140#if BITPIT_ENABLE_MPI==1
141 : m_comm(MPI_COMM_NULL)
142#endif
143 {
144#if BITPIT_ENABLE_MPI==1
145 initialize(logfile, comm);
146#else
147 initialize(logfile);
148#endif
149 reset(true);
150 }
151
157#if BITPIT_ENABLE_MPI==1
161 ParaTree::ParaTree(uint8_t dim, const std::string &logfile, MPI_Comm comm )
162#else
163 ParaTree::ParaTree(uint8_t dim, const std::string &logfile )
164#endif
165 : m_octree(dim), m_trans(dim)
166#if BITPIT_ENABLE_MPI==1
167 , m_comm(MPI_COMM_NULL)
168#endif
169 {
170#if BITPIT_ENABLE_MPI==1
171 initialize(dim, logfile, comm);
172#else
173 initialize(dim, logfile);
174#endif
175 reset(true);
176
177 printHeader();
178 };
179
180 // =============================================================================== //
181
187#if BITPIT_ENABLE_MPI==1
191 ParaTree::ParaTree(std::istream &stream, const std::string &logfile, MPI_Comm comm)
192#else
193 ParaTree::ParaTree(std::istream &stream, const std::string &logfile)
194#endif
195#if BITPIT_ENABLE_MPI==1
196 : m_comm(MPI_COMM_NULL)
197#endif
198 {
199#if BITPIT_ENABLE_MPI==1
200 initialize(logfile, comm);
201#else
202 initialize(logfile);
203#endif
204 restore(stream);
205 }
206
207 // =============================================================================== //
208
211 ParaTree::~ParaTree(){
212#if BITPIT_ENABLE_MPI==1
213 freeComm();
214#endif
215 };
216
217 // =============================================================================== //
218
223 ParaTree::ParaTree(const ParaTree & other)
224 : m_partitionFirstDesc(other.m_partitionFirstDesc),
225 m_partitionLastDesc(other.m_partitionLastDesc),
226 m_partitionRangeGlobalIdx(other.m_partitionRangeGlobalIdx),
227 m_partitionRangeGlobalIdx0(other.m_partitionRangeGlobalIdx0),
228 m_globalNumOctants(other.m_globalNumOctants),
229 m_maxDepth(other.m_maxDepth),
230 m_treeConstants(other.m_treeConstants),
231 m_nofGhostLayers(other.m_nofGhostLayers),
232 m_octree(other.m_octree),
233 m_bordersPerProc(other.m_bordersPerProc),
234 m_internals(other.m_internals),
235 m_pborders(other.m_pborders),
236 m_mapIdx(other.m_mapIdx),
237 m_loadBalanceRanges(other.m_loadBalanceRanges),
238 m_errorFlag(other.m_errorFlag),
239 m_serial(other.m_serial),
240 m_tol(other.m_tol),
241 m_trans(other.m_trans),
242 m_dim(other.m_dim),
243 m_periodic(other.m_periodic),
244 m_status(other.m_status),
245 m_lastOp(other.m_lastOp),
246 m_log(other.m_log)
247#if BITPIT_ENABLE_MPI==1
248 , m_comm(MPI_COMM_NULL)
249#endif
250 {
251#if BITPIT_ENABLE_MPI==1
252 _initializeCommunicator(other.m_comm);
253#else
254 _initializeSerialCommunicator();
255#endif
256 }
257
258 // =================================================================================== //
259 // METHODS
260 // =================================================================================== //
261
262#if BITPIT_ENABLE_MPI==1
270 void
271 ParaTree::_initializeCommunicator(MPI_Comm communicator)
272 {
273 // Communicator can be set just once
274 if (isCommSet()) {
275 throw std::runtime_error ("PABLO communicator can be set just once");
276 }
277
278 // Early return if the communicator is a null communicator
279 if (communicator == MPI_COMM_NULL) {
280 m_comm = MPI_COMM_NULL;
281
282 m_nproc = 1;
283 m_rank = 0;
284
285 return;
286 }
287
288 // Create a copy of the user-specified communicator
289 //
290 // No library routine should use MPI_COMM_WORLD as the communicator;
291 // instead, a duplicate of a user-specified communicator should always
292 // be used.
293 MPI_Comm_dup(communicator, &m_comm);
294
295 // Get MPI information
296 MPI_Comm_size(m_comm, &m_nproc);
297 MPI_Comm_rank(m_comm, &m_rank);
298 }
299#else
303 void
304 ParaTree::_initializeSerialCommunicator()
305 {
306 m_nproc = 1;
307 m_rank = 0;
308 }
309#endif
310
315 void
316 ParaTree::_initializePartitions() {
317 // Create the data structures for storing partion information
318 m_partitionFirstDesc.resize(m_nproc);
319 m_partitionLastDesc.resize(m_nproc);
320 m_partitionRangeGlobalIdx.resize(m_nproc);
321 m_partitionRangeGlobalIdx0.resize(m_nproc);
322 uint64_t lastDescMorton = m_octree.getLastDescMorton();
323 uint64_t firstDescMorton = m_octree.getFirstDescMorton();
324 for(int p = 0; p < m_nproc; ++p){
325 m_partitionRangeGlobalIdx0[p] = 0;
326 m_partitionRangeGlobalIdx[p] = m_globalNumOctants - 1;
327 m_partitionLastDesc[p] = lastDescMorton;
328 m_partitionFirstDesc[p] = firstDescMorton;
329 }
330 }
331
336 void
337 ParaTree::_initialize(uint8_t dim, const std::string &logfile) {
338 // The octree is serial
339 m_serial = true;
340
341 // Initialize the status
342 m_status = 0;
343
344 // Initialize the logger
345 initializeLogger(logfile);
346
347 // Set the dimension to a dummy value
348 setDim(dim);
349
350 // Initialize the global number of octants
351 m_globalNumOctants = 0;
352
353 // Initialize the number of ghost layers
354 m_nofGhostLayers = 1;
355 }
356
361#if BITPIT_ENABLE_MPI==1
365 void
366 ParaTree::initialize(const std::string &logfile, MPI_Comm comm) {
367#else
368 void
369 ParaTree::initialize(const std::string &logfile) {
370#endif
371#if BITPIT_ENABLE_MPI==1
372 // Initialize the communicator
373 _initializeCommunicator(comm);
374#else
375 // Initialize serial communicator
376 _initializeSerialCommunicator();
377#endif
378
379 // Initialized the tree
380 _initialize(0, logfile);
381
382 // Initialize partitions
383 _initializePartitions();
384 }
385
390#if BITPIT_ENABLE_MPI==1
394#endif
395 void
396#if BITPIT_ENABLE_MPI==1
397 ParaTree::initialize(uint8_t dim, const std::string &logfile, MPI_Comm comm) {
398#else
399 ParaTree::initialize(uint8_t dim, const std::string &logfile) {
400#endif
401#if BITPIT_ENABLE_MPI==1
402 // Initialize the communicator
403 _initializeCommunicator(comm);
404#else
405 // Initialize serial communicator
406 _initializeSerialCommunicator();
407#endif
408
409 // Initialized the tree
410 if (dim < 2 || dim > 3) {
411 throw std::runtime_error ("Invalid value for the dimension");
412 }
413
414 _initialize(dim, logfile);
415
416 // Initialize partitions
417 _initializePartitions();
418 }
419
424 void
425 ParaTree::reinitialize(uint8_t dim, const std::string &logfile) {
426 // Initialized the tree
427 if (dim < 2 || dim > 3) {
428 throw std::runtime_error ("Invalid value for the dimension");
429 }
430
431 _initialize(dim, logfile);
432
433 // Initialize partitions
434 _initializePartitions();
435 }
436
439 void
440 ParaTree::initializeLogger(const std::string &logfile){
441 log::manager().create(logfile, false, m_nproc, m_rank);
442 m_log = &log::cout(logfile);
443 }
444
447 void
448 ParaTree::reset(){
449 reset(true);
450 }
451
454 void
455 ParaTree::reset(bool createRoot){
456 m_tol = 1.0e-14;
457 m_serial = true;
458 m_errorFlag = 0;
459
460 if (createRoot) {
461 m_maxDepth = 0;
462 } else {
463 m_maxDepth = -1;
464 }
465
466 m_octree.reset(createRoot);
467 m_globalNumOctants = getNumOctants();
468
469 m_lastOp = OP_INIT;
470
471 m_bordersPerProc.clear();
472 m_internals.clear();
473 m_pborders.clear();
474
475 m_loadBalanceRanges.clear();
476
477 std::fill(m_periodic.begin(), m_periodic.end(), false);
478
479 _initializePartitions();
480 }
481
482 // =============================================================================== //
483
488 int ParaTree::getDumpVersion() const
489 {
490 const int DUMP_VERSION = 1;
491
492 return DUMP_VERSION;
493 }
494
495 // =============================================================================== //
496
502 void ParaTree::dump(std::ostream &stream, bool full)
503 {
504 // Version
505 utils::binary::write(stream, getDumpVersion());
506
507 // Tree data
508 utils::binary::write(stream, getNproc());
509
510 utils::binary::write(stream, getDim());
511
512 utils::binary::write(stream, getSerial());
513 utils::binary::write(stream, getNofGhostLayers());
514 utils::binary::write(stream, getMaxDepth());
515 utils::binary::write(stream, getStatus());
516 utils::binary::write(stream, getBalanceCodimension());
517
518 for (int i = 0; i < m_treeConstants->nFaces; i++) {
519 utils::binary::write(stream, getPeriodic(i));
520 }
521
522 // Octant data
523 uint32_t nOctants = getNumOctants();
524 utils::binary::write(stream, nOctants);
525
526 uint64_t nGlobalOctants = getGlobalNumOctants();
527 utils::binary::write(stream, nGlobalOctants);
528
529 for (uint32_t i = 0; i < nOctants; i++) {
530 const Octant &octant = m_octree.m_octants[i];
531
532 utils::binary::write(stream, octant.getLevel());
533 utils::binary::write(stream, octant.getLogicalCoordinates(0));
534 utils::binary::write(stream, octant.getLogicalCoordinates(1));
535 utils::binary::write(stream, octant.getLogicalCoordinates(2));
536 utils::binary::write(stream, octant.getGhostLayer());
537
538 for (size_t k = 0; k < Octant::INFO_ITEM_COUNT; ++k) {
539 utils::binary::write(stream, (bool) octant.m_info[k]);
540 }
541
542 utils::binary::write(stream, octant.getBalance());
543 utils::binary::write(stream, octant.getMarker());
544 }
545
546 // Information about partitioning
547 for (int k = 0; k < m_nproc; ++k) {
548 utils::binary::write(stream, m_partitionFirstDesc[k]);
549 }
550
551 for (int k = 0; k < m_nproc; ++k) {
552 utils::binary::write(stream, m_partitionLastDesc[k]);
553 }
554
555 for (int k = 0; k < m_nproc; ++k) {
556 utils::binary::write(stream, m_partitionRangeGlobalIdx[k]);
557 }
558
559 // Extended information (mapping, ...)
560 utils::binary::write(stream, full);
561 if (full) {
562 utils::binary::write(stream, m_lastOp);
563 if (m_lastOp == OP_ADAPT_MAPPED){
564 for (auto idx : m_mapIdx) {
565 utils::binary::write(stream, idx);
566 }
567
568 utils::binary::write(stream, m_octree.m_lastGhostBros.size());
569 for (auto lastGhostBrother : m_octree.m_lastGhostBros) {
570 utils::binary::write(stream, lastGhostBrother);
571 }
572 }
573 else if (m_lastOp == OP_LOADBALANCE || m_lastOp == OP_LOADBALANCE_FIRST){
574 for (int i = 0; i < m_nproc; ++i) {
575 utils::binary::write(stream, m_partitionRangeGlobalIdx0[i]);
576 }
577 }
578 }
579
580 }
581
582 // =============================================================================== //
583
588 void ParaTree::restore(std::istream &stream)
589 {
590 // Version
591 int version;
592 utils::binary::read(stream, version);
593 if (version != getDumpVersion()) {
594 throw std::runtime_error ("The version of the file does not match the required version");
595 }
596
597 // Check if the number of processes matches
598 int nProcs;
599 utils::binary::read(stream, nProcs);
600 if (nProcs != m_nproc) {
601 throw std::runtime_error ("The restart was saved with a different number of processes.");
602 }
603
604 // Initialize the tree
605 uint8_t dimension;
606 utils::binary::read(stream, dimension);
607
608 m_octree.initialize(dimension);
609 m_trans.initialize(dimension);
610 reinitialize(dimension, m_log->getName());
611 reset(false);
612
613 // Set tree properties
614 utils::binary::read(stream, m_serial);
615 utils::binary::read(stream, m_nofGhostLayers);
616 utils::binary::read(stream, m_maxDepth);
617 utils::binary::read(stream, m_status);
618
619 bool balanceCodimension;
620 utils::binary::read(stream, balanceCodimension);
621 setBalanceCodimension(balanceCodimension);
622
623 for (int i = 0; i < m_treeConstants->nFaces; i++) {
624 bool periodicBorder;
625 utils::binary::read(stream, periodicBorder);
626 if (periodicBorder){
627 setPeriodic(i);
628 }
629 }
630
631 // Restore octants
632 uint32_t nOctants;
633 utils::binary::read(stream, nOctants);
634
635 uint64_t nGlobalOctants;
636 utils::binary::read(stream, nGlobalOctants);
637 m_globalNumOctants = nGlobalOctants;
638
639 m_octree.m_octants.clear();
640 m_octree.m_octants.reserve(nOctants);
641 for (uint32_t i = 0; i < nOctants; i++) {
642 // Create octant
643 uint8_t level;
644 utils::binary::read(stream, level);
645
646 uint32_t x;
647 utils::binary::read(stream, x);
648
649 uint32_t y;
650 utils::binary::read(stream, y);
651
652 uint32_t z;
653 utils::binary::read(stream, z);
654
655 Octant octant(false, m_dim, level, x, y, z);
656
657 int ghost;
658 utils::binary::read(stream, ghost);
659 octant.setGhostLayer(ghost);
660
661 // Set octant info
662 for (size_t k = 0; k < Octant::INFO_ITEM_COUNT; ++k) {
663 bool bit;
664 utils::binary::read(stream, bit);
665 octant.m_info.set(k, bit);
666 }
667
668 // Set octant 2:1 balance
669 bool balance21;
670 utils::binary::read(stream, balance21);
671 octant.setBalance(balance21);
672
673 // Set marker
674 int8_t marker;
675 utils::binary::read(stream, marker);
676 octant.setMarker(marker);
677
678 // Add octant to the list
679 m_octree.m_octants.push_back(std::move(octant));
680 }
681
682 m_octree.updateLocalMaxDepth();
683
684 // Set first last descendant
685 m_partitionFirstDesc.resize(m_nproc);
686 for (int k = 0; k < m_nproc; ++k) {
687 uint64_t descendant;
688 utils::binary::read(stream, descendant);
689 m_partitionFirstDesc[k] = descendant;
690 }
691 m_octree.m_firstDescMorton = m_partitionFirstDesc[m_rank];
692
693 m_partitionLastDesc.resize(m_nproc);
694 for (int k = 0; k < m_nproc; ++k) {
695 uint64_t descendant;
696 utils::binary::read(stream, descendant);
697 m_partitionLastDesc[k] = descendant;
698 }
699 m_octree.m_lastDescMorton = m_partitionLastDesc[m_rank];
700
701 // Set partitions and parallel information
702 m_partitionRangeGlobalIdx.resize(m_nproc);
703 for (int k = 0; k < m_nproc; ++k) {
704 uint64_t rangeGlobalIdx;
705 utils::binary::read(stream, rangeGlobalIdx);
706 m_partitionRangeGlobalIdx[k] = rangeGlobalIdx;
707 }
708
709#if BITPIT_ENABLE_MPI==1
710 if (!m_serial) {
711 computeGhostHalo();
712 }
713#endif
714
715 // Full restore (i.e. restore with mapper of last operation)
716 m_mapIdx.clear();
717
718 for (int i = 0; i < m_nproc; ++i) {
719 m_partitionRangeGlobalIdx0[i] = 0;
720 }
721
722 bool full;
723 utils::binary::read(stream, full);
724 if (full) {
725 utils::binary::read(stream, m_lastOp);
726 if (m_lastOp == OP_ADAPT_MAPPED) {
727 m_mapIdx.resize(m_octree.m_octants.size());
728 for (size_t i = 0; i < m_octree.m_octants.size(); ++i) {
729 utils::binary::read(stream, m_mapIdx[i]);
730 }
731
732 size_t lastGhostBrosSize;
733 utils::binary::read(stream, lastGhostBrosSize);
734 m_octree.m_lastGhostBros.resize(lastGhostBrosSize);
735 for (size_t i = 0; i < lastGhostBrosSize; ++i) {
736 utils::binary::read(stream, m_octree.m_lastGhostBros[i]);
737 }
738 }
739 else if (m_lastOp == OP_LOADBALANCE || m_lastOp == OP_LOADBALANCE_FIRST){
740 for (int i = 0; i < m_nproc; ++i) {
741 utils::binary::read(stream, m_partitionRangeGlobalIdx0[i]);
742 }
743 }
744 }
745 else {
746 m_lastOp = OP_INIT;
747 }
748 }
749
750 // =============================================================================== //
751
754 void
755 ParaTree::printHeader(){
756 (*m_log) << log::context("PABLO");
757 (*m_log) << "---------------------------------------------" << endl;
758 (*m_log) << "- PABLO PArallel Balanced Linear Octree -" << endl;
759 (*m_log) << "---------------------------------------------" << endl;
760 (*m_log) << " " << endl;
761 (*m_log) << "---------------------------------------------" << endl;
762 (*m_log) << "- PABLO restart -" << endl;
763 (*m_log) << "---------------------------------------------" << endl;
764 (*m_log) << " Number of proc : " + to_string(static_cast<unsigned long long>(m_nproc)) << endl;
765 (*m_log) << " Dimension : " + to_string(static_cast<unsigned long long>(m_dim)) << endl;
766 (*m_log) << " Max allowed level : " + to_string(static_cast<unsigned long long>(m_treeConstants->maxLevel)) << endl;
767 (*m_log) << " Number of octants : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
768 (*m_log) << "---------------------------------------------" << endl;
769 (*m_log) << " " << endl;
770 }
771
772 // =================================================================================== //
773 // BASIC GET/SET METHODS
774 // =================================================================================== //
775
779 uint8_t
780 ParaTree::getDim() const {
781 return m_dim;
782 };
783
787 uint64_t
788 ParaTree::getGlobalNumOctants() const {
789 return m_globalNumOctants;
790 };
791
795 bool
796 ParaTree::getSerial() const {
797 return m_serial;
798 };
799
803 bool
804 ParaTree::getParallel() const {
805 return (!m_serial);
806 };
807
811 int
812 ParaTree::getRank() const {
813 return m_rank;
814 };
815
819 int
820 ParaTree::getNproc() const {
821 return m_nproc;
822 };
823
827 Logger&
828 ParaTree::getLog(){
829 return (*m_log);
830 }
831
835 ParaTree::Operation
836 ParaTree::getLastOperation() const {
837 return m_lastOp;
838 }
839
840#if BITPIT_ENABLE_MPI==1
847 void
848 ParaTree::setComm(MPI_Comm communicator)
849 {
850 // The communicator has to be valid
851 if (communicator == MPI_COMM_NULL) {
852 throw std::runtime_error ("PABLO communicator is not valid");
853 }
854
855 // Initialize the communicator
856 _initializeCommunicator(communicator);
857
858 // Initialize partition data
859 _initializePartitions();
860 }
861
877 void
878 ParaTree::replaceComm(MPI_Comm communicator, MPI_Comm *previousCommunicator)
879 {
880 // The communicator has to be already set
881 if (!isCommSet()) {
882 throw std::runtime_error ("PABLO communicator is not set");
883 }
884
885 // Check if the communicator is valid
886 //
887 // The communicator should be equaivalent to the one currently set,
888 // i.e., the rank of the processes have to be the same in both
889 // communicators.
890 int updatedRank;
891 MPI_Comm_rank(communicator, &updatedRank);
892 int currentRank = getRank();
893
894 int isValid = (updatedRank == currentRank) ? 1 : 0;
895 MPI_Allreduce(MPI_IN_PLACE, &isValid, 1, MPI_INT, MPI_LAND, m_comm);
896 if (!isValid) {
897 throw std::runtime_error ("New communicator is not valid");
898 }
899
900 // Handle previous communicator
901 if (previousCommunicator) {
902 *previousCommunicator = m_comm;
903 m_comm = MPI_COMM_NULL;
904 } else {
905 freeComm();
906 }
907
908 // Set the communicator
909 setComm(communicator);
910 }
911
915 void
916 ParaTree::freeComm() {
917 if (!isCommSet()) {
918 return;
919 }
920
921 // Free MPI communicator
922 int finalizedCalled;
923 MPI_Finalized(&finalizedCalled);
924 if (finalizedCalled) {
925 return;
926 }
927
928 MPI_Comm_free(&m_comm);
929 }
930
935 bool
936 ParaTree::isCommSet() const {
937 return (getComm() != MPI_COMM_NULL);
938 }
939
943 MPI_Comm
944 ParaTree::getComm() const {
945 return m_comm;
946 };
947#endif
948
954 const std::vector<uint64_t> &
955 ParaTree::getPartitionRangeGlobalIdx() const {
956 return m_partitionRangeGlobalIdx;
957 };
958
964 const std::vector<uint64_t> &
965 ParaTree::getPartitionFirstDesc() const {
966 return m_partitionFirstDesc;
967 };
968
974 const std::vector<uint64_t> &
975 ParaTree::getPartitionLastDesc() const {
976 return m_partitionLastDesc;
977 };
978
982 darray3
983 ParaTree::getOrigin() const {
984 return m_trans.m_origin;
985 };
986
990 double
991 ParaTree::getX0() const {
992 return m_trans.m_origin[0];
993 };
994
998 double
999 ParaTree::getY0() const {
1000 return m_trans.m_origin[1];
1001 };
1002
1006 double
1007 ParaTree::getZ0() const {
1008 return m_trans.m_origin[2];
1009 };
1010
1014 double
1015 ParaTree::getL() const {
1016 return m_trans.m_L;
1017 };
1018
1022 int
1023 ParaTree::getMaxLevel() const {
1024 return m_treeConstants->maxLevel;
1025 };
1026
1030 uint32_t
1031 ParaTree::getMaxLength() const {
1032 return m_treeConstants->lengths[0];
1033 }
1034
1038 uint8_t
1039 ParaTree::getNnodes() const {
1040 return m_treeConstants->nNodes;
1041 }
1042
1046 uint8_t
1047 ParaTree::getNfaces() const {
1048 return m_treeConstants->nFaces;
1049 }
1050
1054 uint8_t
1055 ParaTree::getNedges() const {
1056 return m_treeConstants->nEdges;
1057 }
1058
1062 uint8_t
1063 ParaTree::getNchildren() const {
1064 return m_treeConstants->nChildren;
1065 }
1066
1070 uint8_t
1071 ParaTree::getNnodesperface() const {
1072 return m_treeConstants->nNodesPerFace;
1073 }
1074
1078 void
1079 ParaTree::getNormals(int8_t normals[6][3]) const {
1080 for (int i=0; i<6; i++){
1081 for (int j=0; j<3; j++){
1082 normals[i][j] = m_treeConstants->normals[i][j];
1083 }
1084 }
1085 }
1086
1092 void
1093 ParaTree::getOppface(uint8_t oppface[6]) const {
1094 for (int j=0; j<6; j++){
1095 oppface[j] = m_treeConstants->oppositeFace[j];
1096 }
1097 }
1098
1103 void
1104 ParaTree::getFacenode(uint8_t facenode[6][4]) const {
1105 for (int i=0; i<6; i++){
1106 for (int j=0; j<4; j++){
1107 facenode[i][j] = m_treeConstants->faceNode[i][j];
1108 }
1109 }
1110 }
1111
1116 void
1117 ParaTree::getNodeface(uint8_t nodeface[8][3]) const {
1118 for (int i=0; i<8; i++){
1119 for (int j=0; j<3; j++){
1120 nodeface[i][j] = m_treeConstants->nodeFace[i][j];
1121 }
1122 }
1123 }
1124
1129 void
1130 ParaTree::getEdgeface(uint8_t edgeface[12][2]) const {
1131 for (int i=0; i<12; i++){
1132 for (int j=0; j<2; j++){
1133 edgeface[i][j] = m_treeConstants->edgeFace[i][j];
1134 }
1135 }
1136 }
1137
1141 void
1142 ParaTree::getNodecoeffs(int8_t nodecoeffs[8][3]) const {
1143 for (int i=0; i<8; i++){
1144 nodecoeffs[i][2] = 0;
1145 for (int j=0; j<m_dim; j++){
1146 nodecoeffs[i][j] = m_treeConstants->nodeCoeffs[i][j];
1147 }
1148 }
1149 }
1150
1154 void
1155 ParaTree::getEdgecoeffs(int8_t edgecoeffs[12][3]) const {
1156 for (int i=0; i<12; i++){
1157 for (int j=0; j<3; j++){
1158 edgecoeffs[i][j] = m_treeConstants->edgeCoeffs[i][j];
1159 }
1160 }
1161 }
1162
1166 const int8_t
1167 (*ParaTree::getNormals() const) [3] {
1168 return m_treeConstants->normals;
1169 }
1170
1176 const uint8_t
1177 *ParaTree::getOppface() const {
1178 return m_treeConstants->oppositeFace;
1179 }
1180
1185 const uint8_t
1186 (*ParaTree::getFacenode() const)[4] {
1187 return m_treeConstants->faceNode;
1188 }
1189
1194 const uint8_t
1195 (*ParaTree::getNodeface() const)[3] {
1196 return m_treeConstants->nodeFace;
1197 }
1198
1203 const uint8_t
1204 (*ParaTree::getEdgeface() const)[2] {
1205 return m_treeConstants->edgeFace;
1206 }
1207
1211 const int8_t
1212 (*ParaTree::getNodecoeffs() const)[3] {
1213 return m_treeConstants->nodeCoeffs;
1214 };
1215
1219 const int8_t
1220 (*ParaTree::getEdgecoeffs() const)[3] {
1221 return m_treeConstants->edgeCoeffs;
1222 };
1223
1227 bvector
1228 ParaTree::getPeriodic() const {
1229 return m_periodic;
1230 };
1231
1236 bool
1237 ParaTree::getPeriodic(uint8_t i) const {
1238 return m_periodic[i];
1239 };
1240
1243 double
1244 ParaTree::getTol() const {
1245 return m_tol;
1246 };
1247
1251 void
1252 ParaTree::setPeriodic(uint8_t i){
1253 m_periodic[i] = true;
1254 m_periodic[m_treeConstants->oppositeFace[i]] = true;
1255 m_octree.setPeriodic(m_periodic);
1256 };
1257
1261 void
1262 ParaTree::setTol(double tol){
1263 m_tol = tol;
1264 };
1265
1266 // =================================================================================== //
1267 // INDEX BASED METHODS
1268 // =================================================================================== //
1269
1274 darray3
1275 ParaTree::getCoordinates(uint32_t idx) const {
1276 return m_trans.mapCoordinates(m_octree.m_octants[idx].getLogicalCoordinates());
1277 }
1278
1283 double
1284 ParaTree::getX(uint32_t idx) const {
1285 return m_trans.mapX(m_octree.m_octants[idx].getLogicalCoordinates(0));
1286 }
1287
1292 double
1293 ParaTree::getY(uint32_t idx) const {
1294 return m_trans.mapY(m_octree.m_octants[idx].getLogicalCoordinates(1));
1295 }
1296
1301 double
1302 ParaTree::getZ(uint32_t idx) const {
1303 return m_trans.mapZ(m_octree.m_octants[idx].getLogicalCoordinates(2));
1304 }
1305
1310 double
1311 ParaTree::getSize(uint32_t idx) const {
1312 return m_trans.mapSize(m_octree.m_octants[idx].getLogicalSize());
1313 }
1314
1319 double
1320 ParaTree::getArea(uint32_t idx) const {
1321 return m_trans.mapArea(m_octree.m_octants[idx].getLogicalArea());
1322 }
1323
1328 double
1329 ParaTree::getVolume(uint32_t idx) const {
1330 return m_trans.mapVolume(m_octree.m_octants[idx].getLogicalVolume());
1331 }
1332
1337 void
1338 ParaTree::getCenter(uint32_t idx, darray3& centerCoords) const {
1339 darray3 centerCoords_ = m_octree.m_octants[idx].getLogicalCenter();
1340 m_trans.mapCenter(centerCoords_, centerCoords);
1341 }
1342
1347 darray3
1348 ParaTree::getCenter(uint32_t idx) const {
1349 darray3 centerCoords = {{0., 0., 0.}};
1350 darray3 centerCoords_ = m_octree.m_octants[idx].getLogicalCenter();
1351 m_trans.mapCenter(centerCoords_, centerCoords);
1352 return centerCoords;
1353 }
1354
1360 darray3
1361 ParaTree::getFaceCenter(uint32_t idx, uint8_t face) const {
1362 darray3 centerCoords = {{0., 0., 0.}};
1363 darray3 centerCoords_ = m_octree.m_octants[idx].getLogicalFaceCenter(face);
1364 m_trans.mapCenter(centerCoords_, centerCoords);
1365 return centerCoords;
1366 }
1367
1373 void
1374 ParaTree::getFaceCenter(uint32_t idx, uint8_t face, darray3& centerCoords) const {
1375 darray3 centerCoords_ = m_octree.m_octants[idx].getLogicalFaceCenter(face);
1376 m_trans.mapCenter(centerCoords_, centerCoords);
1377 }
1378
1384 darray3
1385 ParaTree::getNode(uint32_t idx, uint8_t node) const {
1386 darray3 nodeCoords = {{0., 0., 0.}};
1387 u32array3 nodeCoords_ = m_octree.m_octants[idx].getLogicalNode(node);
1388 m_trans.mapNode(nodeCoords_, nodeCoords);
1389 return nodeCoords;
1390 }
1391
1397 void
1398 ParaTree::getNode(uint32_t idx, uint8_t node, darray3& nodeCoords) const {
1399 u32array3 nodeCoords_ = m_octree.m_octants[idx].getLogicalNode(node);
1400 m_trans.mapNode(nodeCoords_, nodeCoords);
1401 }
1402
1407 void
1408 ParaTree::getNodes(uint32_t idx, darr3vector & nodesCoords) const {
1409 u32arr3vector nodesCoords_;
1410 m_octree.m_octants[idx].getLogicalNodes(nodesCoords_);
1411 m_trans.mapNodes(nodesCoords_, nodesCoords);
1412 }
1413
1418 darr3vector
1419 ParaTree::getNodes(uint32_t idx) const{
1420 darr3vector nodesCoords;
1421 u32arr3vector nodesCoords_;
1422 m_octree.m_octants[idx].getLogicalNodes(nodesCoords_);
1423 m_trans.mapNodes(nodesCoords_, nodesCoords);
1424 return nodesCoords;
1425 }
1426
1432 void
1433 ParaTree::getNormal(uint32_t idx, uint8_t face, darray3 & normal) const {
1434 i8array3 normal_;
1435 m_octree.m_octants[idx].getNormal(face, normal_, m_treeConstants->normals);
1436 m_trans.mapNormals(normal_, normal);
1437 }
1438
1444 darray3
1445 ParaTree::getNormal(uint32_t idx, uint8_t face) const {
1446 darray3 normal = {{0., 0., 0.}};
1447 i8array3 normal_ = {{0, 0, 0}};
1448 m_octree.m_octants[idx].getNormal(face, normal_, m_treeConstants->normals);
1449 m_trans.mapNormals(normal_, normal);
1450 return normal;
1451 }
1452
1457 int8_t
1458 ParaTree::getMarker(uint32_t idx) const {
1459 return m_octree.getMarker(idx);
1460 };
1461
1467 int8_t
1468 ParaTree::getPreMarker(uint32_t idx){
1469 if (m_lastOp != OP_PRE_ADAPT) {
1470 throw std::runtime_error("Last operation different from preadapt, unable to call getPreMarker function");
1471 }
1472
1473 return m_octree.getMarker(idx);
1474 };
1475
1480 uint8_t
1481 ParaTree::getLevel(uint32_t idx) const {
1482 return m_octree.getLevel(idx);
1483 };
1484
1489 uint64_t
1490 ParaTree::getMorton(uint32_t idx) const {
1491 return m_octree.getMorton(idx);
1492 };
1493
1500 uint64_t
1501 ParaTree::computeNodePersistentKey(uint32_t idx, uint8_t node) const {
1502 return m_octree.computeNodePersistentKey(idx, node);
1503 };
1504
1509 bool
1510 ParaTree::getBalance(uint32_t idx) const {
1511 return m_octree.getBalance(idx);
1512 };
1513
1519 bool
1520 ParaTree::getBound(uint32_t idx, uint8_t face) const {
1521 return m_octree.m_octants[idx].getBound(face);
1522 }
1523
1528 bool
1529 ParaTree::getBound(uint32_t idx) const {
1530 return m_octree.m_octants[idx].getBound();
1531 }
1532
1538 bool
1539 ParaTree::getPbound(uint32_t idx, uint8_t face) const {
1540 return m_octree.m_octants[idx].getPbound(face);
1541 }
1542
1547 bool
1548 ParaTree::getPbound(uint32_t idx) const {
1549 return m_octree.m_octants[idx].getPbound();
1550 }
1551
1556 bool
1557 ParaTree::getIsNewR(uint32_t idx) const {
1558 return m_octree.m_octants[idx].getIsNewR();
1559 };
1560
1565 bool
1566 ParaTree::getIsNewC(uint32_t idx) const {
1567 return m_octree.m_octants[idx].getIsNewC();
1568 };
1569
1574 uint64_t
1575 ParaTree::getGlobalIdx(uint32_t idx) const {
1576 if (m_rank){
1577 return m_partitionRangeGlobalIdx[m_rank-1] + uint64_t(idx + 1);
1578 }
1579 else{
1580 return uint64_t(idx);
1581 };
1582 };
1583
1589 uint32_t
1590 ParaTree::getLocalIdx(uint64_t gidx,int rank) const {
1591 if (rank){
1592 return uint32_t(gidx - 1 - m_partitionRangeGlobalIdx[rank-1]);
1593 }
1594 else{
1595 return uint32_t(gidx);
1596 };
1597 };
1598
1604 void
1605 ParaTree::getLocalIdx(uint64_t gidx, uint32_t & lidx,int & rank) const {
1606 rank = getOwnerRank(gidx);
1607 lidx = getLocalIdx(gidx,rank);
1608 };
1609
1610
1615 uint32_t
1616 ParaTree::getLocalIdx(uint64_t gidx) const {
1617 if (m_rank){
1618 return uint32_t(gidx - 1 - m_partitionRangeGlobalIdx[m_rank-1]);
1619 }
1620 else{
1621 return uint32_t(gidx);
1622 };
1623 };
1624
1629 uint32_t
1630 ParaTree::getGhostLocalIdx(uint64_t gidx) const {
1631
1632 uint32_t index;
1633 typename u64vector::const_iterator findResult;
1634 findResult = std::find(m_octree.m_globalIdxGhosts.begin(),m_octree.m_globalIdxGhosts.end(),gidx);
1635 if(findResult != m_octree.m_globalIdxGhosts.end()){
1636 index = static_cast<uint32_t>(std::distance(m_octree.m_globalIdxGhosts.begin(),findResult));
1637 }
1638 else{
1639 index = std::numeric_limits<uint32_t>::max();
1640 }
1641 return index;
1642 };
1643
1644
1649 uint64_t
1650 ParaTree::getGhostGlobalIdx(uint32_t idx) const {
1651 uint32_t nGhosts = m_octree.getNumGhosts();
1652 if (idx<nGhosts){
1653 return m_octree.m_globalIdxGhosts[idx];
1654 };
1655 return uint64_t(nGhosts);
1656 };
1657
1664 bool
1665 ParaTree::isInternal(uint64_t gidx) const {
1666 if (gidx > m_partitionRangeGlobalIdx[m_rank]){
1667 return false;
1668 }
1669
1670 if (m_rank != 0 && gidx <= m_partitionRangeGlobalIdx[m_rank - 1]){
1671 return false;
1672 }
1673
1674 return true;
1675 };
1676
1682 bitset<72>
1683 ParaTree::getPersistentIdx(uint32_t idx) const {
1684 bitset<72> persistent = getMorton(idx);
1685 bitset<72> level = getLevel(idx);
1686 persistent <<= 8;
1687 persistent |= level;
1688 return persistent;
1689 };
1690
1695 void
1696 ParaTree::setMarker(uint32_t idx, int8_t marker){
1697 if (m_lastOp == OP_PRE_ADAPT) {
1698 throw std::runtime_error("It is not possible to update the tree until the adaption is completed");
1699 }
1700
1701 m_octree.setMarker(idx, marker);
1702 };
1703
1708 void
1709 ParaTree::setBalance(uint32_t idx, bool balance){
1710 if (m_lastOp == OP_PRE_ADAPT) {
1711 throw std::runtime_error("It is not possible to update the tree until the adaption is completed");
1712 }
1713
1714 m_octree.setBalance(idx, balance);
1715 };
1716
1717 // =================================================================================== //
1718 // POINTER BASED METHODS
1719 // =================================================================================== //
1720
1725 darray3
1726 ParaTree::getCoordinates(const Octant* oct) const {
1727 return m_trans.mapCoordinates(oct->getLogicalCoordinates());
1728 }
1729
1734 double
1735 ParaTree::getX(const Octant* oct) const {
1736 return m_trans.mapX(oct->getLogicalCoordinates(0));
1737 }
1738
1743 double
1744 ParaTree::getY(const Octant* oct) const {
1745 return m_trans.mapY(oct->getLogicalCoordinates(1));
1746 }
1747
1752 double
1753 ParaTree::getZ(const Octant* oct) const {
1754 return m_trans.mapZ(oct->getLogicalCoordinates(2));
1755 }
1756
1761 double
1762 ParaTree::getSize(const Octant* oct) const {
1763 return m_trans.mapSize(oct->getLogicalSize());
1764 }
1765
1770 double
1771 ParaTree::getArea(const Octant* oct) const {
1772 return m_trans.mapArea(oct->getLogicalArea());
1773 }
1774
1779 double
1780 ParaTree::getVolume(const Octant* oct) const {
1781 return m_trans.mapVolume(oct->getLogicalVolume());
1782 }
1783
1788 void
1789 ParaTree::getCenter(const Octant* oct, darray3& centerCoords) const {
1790 darray3 centerCoords_ = oct->getLogicalCenter();
1791 m_trans.mapCenter(centerCoords_, centerCoords);
1792 }
1793
1798 darray3
1799 ParaTree::getCenter(const Octant* oct) const {
1800 darray3 centerCoords = {{0., 0., 0.}};
1801 darray3 centerCoords_ = oct->getLogicalCenter();
1802 m_trans.mapCenter(centerCoords_, centerCoords);
1803 return centerCoords;
1804 }
1805
1811 darray3
1812 ParaTree::getFaceCenter(const Octant* oct, uint8_t face) const {
1813 darray3 centerCoords = {{0., 0., 0.}};
1814 darray3 centerCoords_ = oct->getLogicalFaceCenter(face);
1815 m_trans.mapCenter(centerCoords_, centerCoords);
1816 return centerCoords;
1817 }
1818
1824 void
1825 ParaTree::getFaceCenter(const Octant* oct, uint8_t face, darray3& centerCoords) const {
1826 darray3 centerCoords_ = oct->getLogicalFaceCenter(face);
1827 m_trans.mapCenter(centerCoords_, centerCoords);
1828 }
1829
1835 darray3
1836 ParaTree::getNode(const Octant* oct, uint8_t node) const {
1837 darray3 nodeCoords = {{0., 0., 0.}};
1838 u32array3 nodeCoords_ = oct->getLogicalNode(node);
1839 m_trans.mapNode(nodeCoords_, nodeCoords);
1840 return nodeCoords;
1841 }
1842
1848 void
1849 ParaTree::getNode(const Octant* oct, uint8_t node, darray3& nodeCoords) const {
1850 u32array3 nodeCoords_ = oct->getLogicalNode(node);
1851 m_trans.mapNode(nodeCoords_, nodeCoords);
1852 }
1853
1858 void
1859 ParaTree::getNodes(const Octant* oct, darr3vector & nodes) const {
1860 u32arr3vector nodesCoords_;
1861 oct->getLogicalNodes(nodesCoords_);
1862 m_trans.mapNodes(nodesCoords_, nodes);
1863 }
1864
1869 darr3vector
1870 ParaTree::getNodes(const Octant* oct) const {
1871 darr3vector nodes;
1872 u32arr3vector nodesCoords_;
1873 oct->getLogicalNodes(nodesCoords_);
1874 m_trans.mapNodes(nodesCoords_, nodes);
1875 return nodes;
1876 }
1877
1883 void
1884 ParaTree::getNormal(const Octant* oct, uint8_t face, darray3 & normal) const {
1885 i8array3 normal_;
1886 oct->getNormal(face, normal_, m_treeConstants->normals);
1887 m_trans.mapNormals(normal_, normal);
1888 }
1889
1895 darray3
1896 ParaTree::getNormal(const Octant* oct, uint8_t face) const {
1897 darray3 normal = {{0., 0., 0.}};
1898 i8array3 normal_ = {{0, 0, 0}};
1899 oct->getNormal(face, normal_, m_treeConstants->normals);
1900 m_trans.mapNormals(normal_, normal);
1901 return normal;
1902 }
1903
1908 int8_t
1909 ParaTree::getMarker(const Octant* oct) const {
1910 return oct->getMarker();
1911 };
1912
1918 int8_t
1919 ParaTree::getPreMarker(Octant* oct){
1920 if (m_lastOp != OP_PRE_ADAPT) {
1921 throw std::runtime_error("Last operation different from preadapt, unable to call getPreMarker function");
1922 }
1923
1924 return oct->getMarker();
1925 };
1926
1931 uint8_t
1932 ParaTree::getLevel(const Octant* oct) const {
1933 return oct->getLevel();
1934 };
1935
1940 uint64_t
1941 ParaTree::getMorton(const Octant* oct) const {
1942 return oct->getMorton();
1943 };
1944
1949 uint64_t
1950 ParaTree::getLastDescMorton(const Octant* oct) const {
1951 return oct->computeLastDescMorton();
1952 };
1953
1960 uint64_t
1961 ParaTree::computeNodePersistentKey(const Octant* oct, uint8_t node) const {
1962 return oct->computeNodePersistentKey(node);
1963 };
1964
1969 bool
1970 ParaTree::getBalance(const Octant* oct) const {
1971 return oct->getBalance();
1972 };
1973
1979 bool
1980 ParaTree::getBound(const Octant* oct, uint8_t face) const {
1981 return oct->getBound(face);
1982 }
1983
1988 bool
1989 ParaTree::getBound(const Octant* oct) const {
1990 return oct->getBound();
1991 }
1992
1998 bool
1999 ParaTree::getPbound(const Octant* oct, uint8_t face) const {
2000 return oct->getPbound(face);
2001 }
2002
2007 bool
2008 ParaTree::getPbound(const Octant* oct) const {
2009 return oct->getPbound();
2010 }
2011
2016 bool
2017 ParaTree::getIsNewR(const Octant* oct) const {
2018 return oct->getIsNewR();
2019 };
2020
2025 bool
2026 ParaTree::getIsNewC(const Octant* oct) const {
2027 return oct->getIsNewC();
2028 };
2029
2034 uint32_t
2035 ParaTree::getIdx(const Octant* oct) const {
2036#if BITPIT_ENABLE_MPI==1
2037 if (getIsGhost(oct)){
2038 return m_octree.findGhostMorton(oct->getMorton());
2039 }
2040#endif
2041 return m_octree.findMorton(oct->getMorton());
2042 };
2043
2048 uint64_t
2049 ParaTree::getGlobalIdx(const Octant* oct) const {
2050 uint32_t idx = getIdx(oct);
2051#if BITPIT_ENABLE_MPI==1
2052 if (getIsGhost(oct)){
2053 return m_octree.m_globalIdxGhosts[idx];
2054 }
2055#endif
2056 if (m_rank){
2057 return m_partitionRangeGlobalIdx[m_rank-1] + uint64_t(idx + 1);
2058 }
2059 return uint64_t(idx);
2060 };
2061
2067 bitset<72>
2068 ParaTree::getPersistentIdx(const Octant* oct) const {
2069 bitset<72> persistent = getMorton(oct);
2070 bitset<72> level = getLevel(oct);
2071 persistent <<= 8;
2072 persistent |= level;
2073 return persistent;
2074 };
2075
2080 void
2081 ParaTree::setMarker(Octant* oct, int8_t marker){
2082 if (m_lastOp == OP_PRE_ADAPT) {
2083 throw std::runtime_error("It is not possible to update the tree until the adaption is completed");
2084 }
2085
2086 oct->setMarker(marker);
2087 };
2088
2093 void
2094 ParaTree::setBalance(Octant* oct, bool balance){
2095 if (m_lastOp == OP_PRE_ADAPT) {
2096 throw std::runtime_error("It is not possible to update the tree until the adaption is completed");
2097 }
2098
2099 oct->setBalance(balance);
2100 };
2101
2102 // =================================================================================== //
2103 // LOCAL TREE GET/SET METHODS
2104 // =================================================================================== //
2105
2109 uint64_t
2110 ParaTree::getStatus()const {
2111 return m_status;
2112 }
2113
2117 uint32_t
2118 ParaTree::getNumOctants() const{
2119 return m_octree.getNumOctants();
2120 };
2121
2125 uint32_t
2126 ParaTree::getNumGhosts() const{
2127 return m_octree.getNumGhosts();
2128 };
2129
2133 uint32_t
2134 ParaTree::getNumNodes() const{
2135 return m_octree.m_nodes.size();
2136 }
2137
2141 uint8_t
2142 ParaTree::getLocalMaxDepth() const{
2143 return m_octree.getLocalMaxDepth();
2144 };
2145
2149 double
2150 ParaTree::getLocalMinSize() const {
2151 uint32_t size = uint32_t(1)<<(m_treeConstants->maxLevel-m_octree.getLocalMaxDepth());
2152 return m_trans.mapSize(size);
2153 };
2154
2158 double
2159 ParaTree::getLocalMaxSize() const {
2160 uint32_t nocts = getNumOctants();
2161 double octSize = 0;
2162 double size = 0;
2163 for (uint32_t idx = 0; idx < nocts; idx++){
2164 octSize = getSize(idx);
2165 if (octSize > size) size = octSize;
2166 }
2167 return octSize;
2168 };
2169
2173 uint8_t
2174 ParaTree::getBalanceCodimension() const{
2175 return m_octree.getBalanceCodim();
2176 };
2177
2181 uint64_t ParaTree::getFirstDescMorton() const {
2182 return m_octree.getFirstDescMorton();
2183 };
2184
2188 uint64_t ParaTree::getLastDescMorton() const {
2189 return m_octree.getLastDescMorton();
2190 };
2191
2196 uint64_t
2197 ParaTree::getLastDescMorton(uint32_t idx) const {
2198 return m_octree.m_octants[idx].computeLastDescMorton();
2199 };
2200
2204 octantIterator
2205 ParaTree::getInternalOctantsBegin() {
2206 return m_internals.begin();
2207 }
2208
2212 octantIterator
2213 ParaTree::getInternalOctantsEnd() {
2214 return m_internals.end();
2215 }
2216
2220 octantIterator
2221 ParaTree::getPboundOctantsBegin() {
2222 return m_pborders.begin();
2223 }
2224
2228 octantIterator
2229 ParaTree::getPboundOctantsEnd() {
2230 return m_pborders.end();
2231 }
2232
2236 void
2237 ParaTree::setBalanceCodimension(uint8_t b21codim){
2238 if (m_lastOp == OP_PRE_ADAPT) {
2239 throw std::runtime_error("It is not possible to update the tree until the adaption is completed");
2240 }
2241
2242 m_octree.setBalanceCodim(b21codim);
2243 };
2244
2245 // =================================================================================== //
2246 // INTERSECTION GET/SET METHODS
2247 // =================================================================================== //
2248
2252 uint32_t
2253 ParaTree::getNumIntersections() const {
2254 return m_octree.m_intersections.size();
2255 }
2256
2261 Intersection*
2262 ParaTree::getIntersection(uint32_t idx) {
2263 if (idx < m_octree.m_intersections.size()){
2264 return &m_octree.m_intersections[idx];
2265 }
2266 return NULL;
2267 }
2268
2273 uint8_t
2274 ParaTree::getLevel(const Intersection* inter) const {
2275 if(inter->m_finer && inter->m_isghost)
2276 return m_octree.extractGhostOctant(inter->m_owners[inter->m_finer]).getLevel();
2277 else
2278 return m_octree.extractOctant(inter->m_owners[inter->m_finer]).getLevel();
2279 }
2280
2285 bool
2286 ParaTree::getFiner(const Intersection* inter) const {
2287 return inter->m_finer;
2288 }
2289
2294 bool
2295 ParaTree::getBound(const Intersection* inter) const {
2296 return inter->getBound();
2297 }
2298
2303 bool
2304 ParaTree::getIsGhost(const Intersection* inter) const {
2305 return inter->getIsGhost();
2306 }
2307
2312 bool
2313 ParaTree::getPbound(const Intersection* inter) const {
2314 return inter->getPbound();
2315 }
2316
2321 uint8_t
2322 ParaTree::getFace(const Intersection* inter) const {
2323 return inter->m_iface;
2324 }
2325
2330 u32vector
2331 ParaTree::getOwners(const Intersection* inter) const {
2332 u32vector owners(2);
2333 owners[0] = inter->m_owners[0];
2334 owners[1] = inter->m_owners[1];
2335 return owners;
2336 }
2337
2342 uint32_t
2343 ParaTree::getIn(const Intersection* inter) const {
2344 return inter->getIn();
2345 }
2346
2351 uint32_t
2352 ParaTree::getOut(const Intersection* inter) const {
2353 return inter->getOut();
2354 }
2355
2360 bool
2361 ParaTree::getOutIsGhost(const Intersection* inter) const {
2362 return inter->getOutIsGhost();
2363 }
2364
2369 double
2370 ParaTree::getSize(const Intersection* inter) const {
2371 uint32_t Size;
2372 if(inter->m_finer && inter->m_isghost)
2373 Size = m_octree.extractGhostOctant(inter->m_owners[inter->m_finer]).getLogicalSize();
2374 else
2375 Size = m_octree.extractOctant(inter->m_owners[inter->m_finer]).getLogicalSize();
2376 return m_trans.mapSize(Size);
2377 }
2378
2383 double
2384 ParaTree::getArea(const Intersection* inter) const {
2385 uint64_t Area;
2386 if(inter->m_finer && inter->m_isghost)
2387 Area = m_octree.extractGhostOctant(inter->m_owners[1]).getLogicalArea();
2388 else
2389 Area = m_octree.extractOctant(inter->m_owners[inter->m_finer]).getLogicalArea();
2390 return m_trans.mapArea(Area);
2391 }
2392
2397 darray3
2398 ParaTree::getCenter(const Intersection* inter) const {
2399 Octant oct(m_dim);
2400 if(inter->m_finer && inter->m_isghost)
2401 oct = m_octree.extractGhostOctant(inter->m_owners[inter->m_finer]);
2402 else
2403 oct = m_octree.extractOctant(inter->m_owners[inter->m_finer]);
2404
2405 darray3 center = {{0., 0., 0.}};
2406 darray3 centerCoords_ = oct.getLogicalCenter();
2407 int sign = ( int(2*((inter->m_iface)%2)) - 1);
2408 double deplace = double (sign * int(oct.getLogicalSize())) / 2;
2409 centerCoords_[inter->m_iface/2] = uint32_t(int(centerCoords_[inter->m_iface/2]) + deplace);
2410 m_trans.mapCenter(centerCoords_, center);
2411 return center;
2412 }
2413
2418 darr3vector
2419 ParaTree::getNodes(const Intersection* inter) const {
2420 darr3vector nodes;
2421 Octant oct(m_dim);
2422 if(inter->m_finer && inter->m_isghost)
2423 oct = m_octree.extractGhostOctant(inter->m_owners[inter->m_finer]);
2424 else
2425 oct = m_octree.extractOctant(inter->m_owners[inter->m_finer]);
2426 uint8_t face = inter->m_iface;
2427 u32arr3vector nodesCoords_all;
2428 oct.getLogicalNodes(nodesCoords_all);
2429 u32arr3vector nodesCoords_(m_treeConstants->nNodesPerFace);
2430 for (int i=0; i<m_treeConstants->nNodesPerFace; i++){
2431 for (int j=0; j<3; j++){
2432 nodesCoords_[i][j] = nodesCoords_all[m_treeConstants->faceNode[face][i]][j];
2433 }
2434 }
2435 m_trans.mapNodesIntersection(nodesCoords_, nodes);
2436 return nodes;
2437 }
2438
2443 darray3
2444 ParaTree::getNormal(const Intersection* inter) const {
2445 Octant oct(m_dim);
2446 if(inter->m_finer && inter->m_isghost)
2447 oct = m_octree.extractGhostOctant(inter->m_owners[inter->m_finer]);
2448 else
2449 oct = m_octree.extractOctant(inter->m_owners[inter->m_finer]);
2450
2451 uint8_t face = inter->m_iface;
2452 darray3 normal = {{0., 0., 0.}};
2453 i8array3 normal_ = {{0, 0, 0}};
2454 oct.getNormal(face, normal_, m_treeConstants->normals);
2455 m_trans.mapNormals(normal_, normal);
2456 return normal;
2457 }
2458
2459 // =================================================================================== //
2460 // OTHER GET/SET METHODS
2461 // =================================================================================== //
2462
2468 Octant*
2469 ParaTree::getOctant(uint32_t idx) {
2470 return (m_octree.m_octants.data() + idx);
2471 };
2472
2478 const Octant*
2479 ParaTree::getOctant(uint32_t idx) const {
2480 return (m_octree.m_octants.data() + idx);
2481 };
2482
2488 Octant*
2489 ParaTree::getGhostOctant(uint32_t idx) {
2490 return (m_octree.m_ghosts.data() + idx);
2491 };
2492
2498 const Octant*
2499 ParaTree::getGhostOctant(uint32_t idx) const {
2500 return (m_octree.m_ghosts.data() + idx);
2501 };
2502
2507 bool
2508 ParaTree::getIsGhost(const Octant* oct) const {
2509 return oct->getIsGhost();
2510 };
2511
2517 int
2518 ParaTree::getGhostLayer(const Octant* oct) const {
2519 return oct->getGhostLayer();
2520 };
2521
2525 const ParaTree::LoadBalanceRanges &
2526 ParaTree::getLoadBalanceRanges() const {
2527 return m_loadBalanceRanges;
2528 };
2529
2533 std::size_t
2534 ParaTree::getNofGhostLayers() const {
2535 return m_nofGhostLayers;
2536 };
2537
2541 void
2542 ParaTree::setNofGhostLayers(std::size_t nofGhostLayers) {
2543 if (m_nofGhostLayers == 0) {
2544 throw std::runtime_error ("It is not possible to disable the ghost halo!");
2545 }
2546
2547 typedef std::result_of<decltype(&Octant::getGhostLayer)(Octant)>::type layer_t;
2548 typedef std::make_unsigned<layer_t>::type ulayer_t;
2549 ulayer_t maxNofGhostLayers = std::numeric_limits<layer_t>::max() + (ulayer_t) 1;
2550 if (nofGhostLayers > maxNofGhostLayers) {
2551 throw std::runtime_error ("Halo size exceeds the maximum allowed value.");
2552 }
2553
2554 m_nofGhostLayers = nofGhostLayers;
2555 };
2556
2560 const std::map<int, u32vector> & ParaTree::getBordersPerProc() const {
2561 return m_bordersPerProc;
2562 }
2563
2564 // =================================================================================== //
2565 // PRIVATE GET/SET METHODS
2566 // =================================================================================== //
2567
2570 void
2571 ParaTree::setDim(uint8_t dim){
2572 m_dim = dim;
2573 if (m_dim != 0) {
2574 m_treeConstants = &(TreeConstants::instance(m_dim));
2575 m_periodic.resize(m_treeConstants->nFaces, false);
2576 } else {
2577 m_treeConstants = nullptr;
2578 m_periodic.clear();
2579 }
2580 }
2581
2582#if BITPIT_ENABLE_MPI==1
2585 void
2586 ParaTree::updateGlobalFirstDescMorton(){
2587
2588 // Exchange first descendant information
2589 uint64_t firstDescMorton = m_octree.getFirstDescMorton();
2590 m_errorFlag = MPI_Allgather(&firstDescMorton,1,MPI_UINT64_T,m_partitionFirstDesc.data(),1,MPI_UINT64_T,m_comm);
2591
2592 // Fix first descendant for empty partitions
2593 int pp = m_nproc-1;
2594 if (m_partitionRangeGlobalIdx[pp] == m_partitionRangeGlobalIdx[pp-1]){
2595 m_partitionFirstDesc[pp] = std::numeric_limits<uint64_t>::max();
2596 if (m_rank == pp){
2597 m_octree.m_firstDescMorton = std::numeric_limits<uint64_t>::max();
2598 }
2599 }
2600 for (int p = pp-1; p > 0; --p){
2601 if (m_partitionRangeGlobalIdx[p] == m_partitionRangeGlobalIdx[p-1]){
2602 m_partitionFirstDesc[p] = m_partitionFirstDesc[p+1];
2603 if (m_rank == p){
2604 m_octree.m_firstDescMorton = m_partitionFirstDesc[p+1];
2605 }
2606 }
2607 }
2608 };
2609
2612 void
2613 ParaTree::updateGlobalLasttDescMorton(){
2614
2615 // Exchange last descendant information
2616 uint64_t lastDescMorton = m_octree.getLastDescMorton();
2617 m_errorFlag = MPI_Allgather(&lastDescMorton,1,MPI_UINT64_T,m_partitionLastDesc.data(),1,MPI_UINT64_T,m_comm);
2618
2619 // Fix first descendant for empty partitions
2620 //
2621 // Attention rank = 0 can't be empty
2622 for (int p = 1; p < m_nproc; ++p){
2623 if (m_partitionRangeGlobalIdx[p] == m_partitionRangeGlobalIdx[p-1]){
2624 m_partitionLastDesc[p] = m_partitionLastDesc[p-1];
2625 if (m_rank == p){
2626 m_octree.m_lastDescMorton = m_partitionLastDesc[p-1];
2627 }
2628 }
2629 }
2630 };
2631#endif
2632
2633 // =================================================================================== //
2634 // OTHER METHODS //
2635 // =================================================================================== //
2636
2637 // =================================================================================== //
2638 // OTHER OCTANT BASED METHODS //
2639 // =================================================================================== //
2640
2651 void
2652 ParaTree::findNeighbours(const Octant* oct, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours, bvector & isghost, bool onlyinternal, bool append) const{
2653
2654 if (entityCodim == 1){
2655 m_octree.findNeighbours(oct, entityIdx, neighbours, isghost, onlyinternal, append);
2656 }
2657 else if (entityCodim == 2 && m_dim == 3){
2658 m_octree.findEdgeNeighbours(oct, entityIdx, neighbours, isghost, onlyinternal, append);
2659 }
2660 else if (entityCodim == m_dim){
2661 m_octree.findNodeNeighbours(oct, entityIdx, neighbours, isghost, onlyinternal, append);
2662 }
2663 else {
2664 neighbours.clear();
2665 isghost.clear();
2666 }
2667
2668 };
2669
2678 void
2679 ParaTree::findNeighbours(uint32_t idx, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours, bvector & isghost) const {
2680
2681 const Octant* oct = &m_octree.m_octants[idx];
2682
2683 findNeighbours(oct, entityIdx, entityCodim, neighbours, isghost, false, false);
2684
2685 };
2686
2693 void
2694 ParaTree::findNeighbours(uint32_t idx, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours) const {
2695
2696 const Octant* oct = &m_octree.m_octants[idx];
2697 bvector isghost;
2698 findNeighbours(oct, entityIdx, entityCodim, neighbours, isghost, true, false);
2699
2700 };
2701
2710 void
2711 ParaTree::findNeighbours(const Octant* oct, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours, bvector & isghost) const {
2712
2713 findNeighbours(oct, entityIdx, entityCodim, neighbours, isghost, false, false);
2714
2715 };
2716
2724 void
2725 ParaTree::findGhostNeighbours(uint32_t idx, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours) const {
2726
2727 const Octant* oct = &m_octree.m_ghosts[idx];
2728 bvector isghost;
2729 findNeighbours(oct, entityIdx, entityCodim, neighbours, isghost, true, false);
2730
2731 };
2732
2742 void
2743 ParaTree::findGhostNeighbours(uint32_t idx, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours, bvector & isghost) const {
2744
2745 const Octant* oct = &m_octree.m_ghosts[idx];
2746 findNeighbours(oct, entityIdx, entityCodim, neighbours, isghost, false, false);
2747
2748 };
2749
2759 void
2760 ParaTree::findGhostNeighbours(const Octant* oct, uint8_t entityIdx, uint8_t entityCodim, u32vector & neighbours, bvector & isghost) const {
2761
2762 findNeighbours(oct, entityIdx, entityCodim, neighbours, isghost, false, false);
2763
2764 };
2765
2772 void
2773 ParaTree::findAllNodeNeighbours(const Octant* oct, uint32_t node, u32vector & neighbours, bvector & isghost) const {
2774
2775 int dim = getDim();
2776 int octantLevel = getLevel(oct);
2777
2778 u32vector codimNeighnours;
2779 bvector codimIsGhost;
2780
2781 // Get vertex neighbours
2782 findNeighbours(oct, node, m_dim, neighbours, isghost, false, false);
2783
2784 // Get edge neighbours
2785 //
2786 // On non uniform trees the vertex can be inside the edge of the
2787 // neighbour (hanging nodes). To correctly consider these neighbours,
2788 // the following logic can be used:
2789 // - if an edge neighbour has the same level or a lower level than
2790 // the current cell, then it certainly is also a vertex neighbour;
2791 // - if a edge neighbour has a higher level than the current cell,
2792 // it is necessary to check if the neighbour actually contains the
2793 // vertex.
2794 if (dim == 3) {
2795 for (int edge : m_treeConstants->nodeEdge[node]) {
2796 findNeighbours(oct, edge, 2, codimNeighnours, codimIsGhost, false, false);
2797 for (std::size_t i = 0; i < codimNeighnours.size(); ++i) {
2798 const Octant* neighOctant;
2799 if (codimIsGhost[i]==0) {
2800 neighOctant = &m_octree.m_octants[codimNeighnours[i]];
2801 }
2802 else {
2803 neighOctant = &m_octree.m_ghosts[codimNeighnours[i]];
2804 }
2805 int neighOctantLevel = getLevel(neighOctant);
2806 if (neighOctantLevel <= octantLevel) {
2807 neighbours.push_back(codimNeighnours[i]);
2808 isghost.push_back(codimIsGhost[i]);
2809 } else if (isNodeOnOctant(oct, node, neighOctant)) {
2810 neighbours.push_back(codimNeighnours[i]);
2811 isghost.push_back(codimIsGhost[i]);
2812 }
2813 }
2814 }
2815 }
2816
2817 // Get face neighbours
2818 //
2819 // On non uniform trees the vertex can be inside the face of the
2820 // neighbour (hanging nodes). To correctly consider these neighbours,
2821 // the following logic can be used:
2822 // - if an face neighbour has the same level or a lower level than
2823 // the current cell, then it certainly is also a vertex neighbour;
2824 // - if a face neighbour has a higher level than the current cell,
2825 // it is necessary to check if the neighbour actually contains the
2826 // vertex.
2827 for (int j = 0; j < dim; ++j) {
2828 int face = m_treeConstants->nodeFace[node][j];
2829 findNeighbours(oct, face, 1, codimNeighnours, codimIsGhost, false, false);
2830 for (std::size_t i = 0; i < codimNeighnours.size(); ++i) {
2831 const Octant* neighOctant;
2832 if (codimIsGhost[i]==0) {
2833 neighOctant = &m_octree.m_octants[codimNeighnours[i]];
2834 }
2835 else {
2836 neighOctant = &m_octree.m_ghosts[codimNeighnours[i]];
2837 }
2838 int neighOctantLevel = getLevel(neighOctant);
2839 if (neighOctantLevel <= octantLevel) {
2840 neighbours.push_back(codimNeighnours[i]);
2841 isghost.push_back(codimIsGhost[i]);
2842 } else if (isNodeOnOctant(oct, node, neighOctant)) {
2843 neighbours.push_back(codimNeighnours[i]);
2844 isghost.push_back(codimIsGhost[i]);
2845 }
2846 }
2847 }
2848 };
2849
2856 void
2857 ParaTree::findAllNodeNeighbours(uint32_t idx, uint32_t node, u32vector & neighbours, bvector & isghost) {
2858
2859 Octant* oct = getOctant(idx);
2860 findAllNodeNeighbours(oct, node, neighbours, isghost);
2861 };
2862
2871 void
2872 ParaTree::findAllCodimensionNeighbours(uint32_t idx, u32vector & neighbours, bvector & isghost){
2873 Octant* oct = getOctant(idx);
2874 findAllCodimensionNeighbours(oct,neighbours,isghost);
2875 }
2876
2885 void
2886 ParaTree::findAllCodimensionNeighbours(const Octant* oct, u32vector & neighbours, bvector & isghost){
2887 std::array<uint8_t, 4> nCodimensionItems;
2888 nCodimensionItems[0] = 0;
2889 nCodimensionItems[1] = getNfaces();
2890 if(m_dim == 3){
2891 nCodimensionItems[2] = getNedges();
2892 }
2893 nCodimensionItems[m_dim] = getNnodes();
2894
2895 neighbours.clear();
2896 isghost.clear();
2897
2898 int initalCapacity = uipow(3, m_dim) - 1;
2899 neighbours.reserve(initalCapacity);
2900 isghost.reserve(initalCapacity);
2901
2902 for(uint8_t codim = 1; codim <= m_dim; ++codim){
2903 for(int item = 0; item < nCodimensionItems[codim]; ++item){
2904 findNeighbours(oct,item,codim,neighbours,isghost,false,true);
2905 }
2906 }
2907 }
2908
2917 void
2918 ParaTree::findGhostAllCodimensionNeighbours(uint32_t idx, u32vector & neighbours, bvector & isghost){
2919 Octant* oct = getGhostOctant(idx);
2920 findGhostAllCodimensionNeighbours(oct,neighbours,isghost);
2921 }
2922
2931 void
2932 ParaTree::findGhostAllCodimensionNeighbours(Octant* oct, u32vector & neighbours, bvector & isghost){
2933 findAllCodimensionNeighbours(oct, neighbours, isghost);
2934 }
2935
2940 Octant*
2941 ParaTree::getPointOwner(const dvector &point){
2942 uint32_t idx = getPointOwnerIdx(point);
2943 if(idx < numeric_limits<uint32_t>::max())
2944 return &m_octree.m_octants[idx];
2945 else
2946 return NULL;
2947 };
2948
2954 Octant*
2955 ParaTree::getPointOwner(const dvector &point, bool & isghost){
2956 uint32_t idx = getPointOwnerIdx(point,isghost);
2957 if(idx < numeric_limits<uint32_t>::max())
2958 if(isghost)
2959 return &m_octree.m_ghosts[idx];
2960 else
2961 return &m_octree.m_octants[idx];
2962 else
2963 return NULL;
2964 };
2965
2966
2971 Octant*
2972 ParaTree::getPointOwner(const darray3 &point){
2973 uint32_t idx = getPointOwnerIdx(point);
2974 if(idx < numeric_limits<uint32_t>::max())
2975 return &m_octree.m_octants[idx];
2976 else
2977 return NULL;
2978 };
2979
2985 Octant*
2986 ParaTree::getPointOwner(const darray3 &point, bool & isghost){
2987 uint32_t idx = getPointOwnerIdx(point,isghost);
2988 if(idx < numeric_limits<uint32_t>::max())
2989 if(isghost)
2990 return &m_octree.m_ghosts[idx];
2991 else
2992 return &m_octree.m_octants[idx];
2993 else
2994 return NULL;
2995 };
2996
3001 uint32_t
3002 ParaTree::getPointOwnerIdx(const darray3 &point) const {
3003 return getPointOwnerIdx(point.data());
3004 }
3005
3011 uint32_t
3012 ParaTree::getPointOwnerIdx(const darray3 &point, bool & isghost) const {
3013 return getPointOwnerIdx(point.data(),isghost);
3014 }
3015
3020 uint32_t
3021 ParaTree::getPointOwnerIdx(const dvector &point) const {
3022 assert(point.size() >= 3);
3023 return getPointOwnerIdx(point.data());
3024 };
3025
3031 uint32_t
3032 ParaTree::getPointOwnerIdx(const dvector &point, bool & isghost) const {
3033 assert(point.size() >= 3);
3034 return getPointOwnerIdx(point.data(),isghost);
3035 };
3036
3037
3042 uint32_t
3043 ParaTree::getPointOwnerIdx(const double * point) const {
3044 // Evaluate the Morton associated with the point
3045 uint64_t morton = evalPointAnchorMorton(point);
3046 if (morton == PABLO::INVALID_MORTON) {
3047 return numeric_limits<uint32_t>::max();
3048 }
3049
3050 // Early return if the point belongs to a ghost
3051 if (findOwner(morton) != m_rank && (!m_serial)) {
3052 return numeric_limits<uint32_t>::max();
3053 }
3054
3055 // Identify the octant that contains the point
3056 //
3057 // If the owner of the point is an internal octant, decrementing the Morton upper bound
3058 // by one will give us the requested octant.
3059 uint32_t pointOwnerIdx;
3060 uint64_t pointOwnerMorton;
3061 m_octree.findMortonUpperBound(morton, m_octree.m_octants, &pointOwnerIdx, &pointOwnerMorton);
3062 --pointOwnerIdx;
3063
3064 return pointOwnerIdx;
3065 };
3066
3072 uint32_t
3073 ParaTree::getPointOwnerIdx(const double * point, bool & isghost) const {
3074 uint64_t morton = evalPointAnchorMorton(point);
3075 if (morton == PABLO::INVALID_MORTON) {
3076 isghost = false;
3077
3078 return numeric_limits<uint32_t>::max();
3079 }
3080
3081 // Check if the point belongs to a ghost
3082 isghost = (findOwner(morton) != m_rank);
3083
3084 // Identify the octant that contains the point
3085 //
3086 // If the owner of the point is an internal octant, decrementing the Morton upper bound
3087 // by one will give us the requested octant. If the owner of the point is a ghost octant,
3088 // we need to test all the octants before the Morton upper bound.
3089 uint32_t pointOwnerIdx;
3090 uint64_t pointOwnerMorton;
3091 if (!isghost) {
3092 m_octree.findMortonUpperBound(morton, m_octree.m_octants, &pointOwnerIdx, &pointOwnerMorton);
3093 --pointOwnerIdx;
3094 } else {
3095 m_octree.findMortonUpperBound(morton, m_octree.m_ghosts, &pointOwnerIdx, &pointOwnerMorton);
3096 while (pointOwnerIdx > 0) {
3097 --pointOwnerIdx;
3098
3099 const Octant *octant = getGhostOctant(pointOwnerIdx);
3100 darray3 octantAnchor = getCoordinates(octant);
3101 double octantSize = getSize(octant);
3102
3103 bool found = true;
3104 for (int i = 0; i < m_dim; ++i) {
3105 if (point[i] < octantAnchor[i] - m_tol || point[i] > (octantAnchor[i] + octantSize + m_tol)) {
3106 found = false;
3107 break;
3108 }
3109 }
3110
3111 if (found) {
3112 return pointOwnerIdx;
3113 }
3114 }
3115
3116 return numeric_limits<uint32_t>::max();
3117 }
3118
3119 return pointOwnerIdx;
3120 };
3121
3126 int
3127 ParaTree::getPointOwnerRank(const darray3 &point){
3128 // Early return if the point is not associatd with a valid Morton
3129 uint64_t morton = evalPointAnchorMorton(point.data());
3130 if (morton == PABLO::INVALID_MORTON) {
3131 return -1;
3132 }
3133
3134 // Identify that partition that contains the point
3135 return findOwner(morton);
3136 };
3137
3144 uint64_t
3145 ParaTree::evalPointAnchorMorton(const double * point) const {
3146 // Early return if the tree is empty
3147 if (m_octree.m_octants.empty()) {
3148 return PABLO::INVALID_MORTON;
3149 }
3150
3151 // Early return if the point is outside the domain
3152 if (point[0] > 1+m_tol || point[1] > 1+m_tol || point[2] > 1+m_tol) {
3153 return PABLO::INVALID_MORTON;
3154 } else if (point[0] < -m_tol || point[1] < -m_tol || point[2] < -m_tol) {
3155 return PABLO::INVALID_MORTON;
3156 }
3157
3158 // Evaluate the Morton associated to the point
3159 uint32_t x = m_trans.mapX(std::min(std::max(point[0], 0.0), 1.0));
3160 uint32_t y = m_trans.mapY(std::min(std::max(point[1], 0.0), 1.0));
3161 uint32_t z = m_trans.mapZ(std::min(std::max(point[2], 0.0), 1.0));
3162
3163 uint32_t maxLength = getMaxLength();
3164 if (x == maxLength) x = x - 1;
3165 if (y == maxLength) y = y - 1;
3166 if (z == maxLength) z = z - 1;
3167
3168 return PABLO::computeMorton(m_dim, x, y, z);
3169 }
3170
3176 uint8_t
3177 ParaTree::getFamilySplittingNode(const Octant* oct) const {
3178 return oct->getFamilySplittingNode();
3179 };
3180
3186 void
3187 ParaTree::expectedOctantAdapt(const Octant* oct, int8_t marker, octvector* result) const {
3188
3189 if (oct == nullptr || result == nullptr){
3190 return;
3191 }
3192
3193 if (marker > 0){
3194 uint8_t nChildren = oct->countChildren();
3195 result->resize(nChildren);
3196 oct->buildChildren(result->data());
3197 }
3198 else if (marker < 0){
3199 result->clear();
3200 result->push_back(oct->buildFather());
3201 }
3202
3203 result->push_back(*oct);
3204 };
3205
3214 void
3215 ParaTree::getMapping(uint32_t idx, u32vector & mapper, bvector & isghost) const {
3216
3217 if (idx >= m_mapIdx.size()){
3218 throw std::runtime_error ("Invalid value for input index in getMapping");
3219 }
3220
3221 // Coarsening has to be handled separately, all other changes can just
3222 // return the value stored in the mapper.
3223 if (getIsNewC(idx)){
3224 // Count the children
3225 int nChildren = m_treeConstants->nChildren;
3226
3227 int nInternalChildren = nChildren;
3228 if (idx == (getNumOctants() - 1)) {
3229 nInternalChildren -= m_octree.m_lastGhostBros.size();
3230 }
3231
3232 // Fill the mapper
3233 mapper.resize(nChildren);
3234 isghost.resize(nChildren);
3235
3236 for (int i = 0; i < nInternalChildren; i++){
3237 mapper[i] = m_mapIdx[idx] + i;
3238 isghost[i] = false;
3239 }
3240
3241 for (int i = nInternalChildren; i < nChildren; i++){
3242 mapper[i] = m_octree.m_lastGhostBros[i - nInternalChildren];
3243 isghost[i] = true;
3244 }
3245 } else {
3246 mapper.resize(1);
3247 isghost.resize(1);
3248
3249 mapper[0] = m_mapIdx[idx];
3250 isghost[0] = false;
3251 }
3252
3253 };
3254
3264 void
3265 ParaTree::getMapping(uint32_t idx, u32vector & mapper, bvector & isghost, ivector & rank) const {
3266
3267 if (m_lastOp == OP_ADAPT_MAPPED){
3268 getMapping(idx, mapper, isghost);
3269 int n = isghost.size();
3270 rank.resize(n);
3271 for (int i=0; i<n; i++){
3272 rank[i] = m_rank;
3273 }
3274 }
3275 else if (m_lastOp == OP_LOADBALANCE || m_lastOp == OP_LOADBALANCE_FIRST){
3276 mapper.resize(1);
3277 isghost.resize(1);
3278 rank.resize(1);
3279 uint64_t gidx = getGlobalIdx(idx);
3280 uint64_t previousIdx = gidx;
3281 for (int iproc=0; iproc<m_nproc; ++iproc){
3282 if (m_partitionRangeGlobalIdx0[iproc]>=gidx){
3283 if (iproc > 0)
3284 previousIdx -= m_partitionRangeGlobalIdx0[iproc-1] + 1;
3285 rank[0] = (m_lastOp == OP_LOADBALANCE_FIRST ? m_rank : iproc);
3286 isghost[0] = false;
3287 break;
3288 }
3289 }
3290 mapper[0] = static_cast<uint32_t>(previousIdx);
3291 }
3292 };
3293
3300 void
3301 ParaTree::getPreMapping(u32vector & idx, vector<int8_t> & markers, vector<bool> & isghost)
3302 {
3303 if (m_lastOp != OP_PRE_ADAPT) {
3304 throw std::runtime_error("Last operation different from preadapt, unable to call getPreMarker function");
3305 }
3306
3307 idx.clear();
3308 markers.clear();
3309 isghost.clear();
3310
3311 std::size_t firstGsize = m_octree.m_firstGhostBros.size();
3312 std::size_t lastGsize = m_octree.m_lastGhostBros.size();
3313
3314 idx.reserve(getNumOctants() + firstGsize + lastGsize);
3315 markers.reserve(getNumOctants() + firstGsize + lastGsize);
3316 isghost.reserve(getNumOctants() + firstGsize + lastGsize);
3317
3318 //insert first ghost brothers if present
3319 {
3320 for (uint32_t id : m_octree.m_firstGhostBros){
3321 int8_t marker = m_octree.m_ghosts[id].getMarker();
3322 idx.push_back(id);
3323 markers.push_back(marker);
3324 isghost.push_back(true);
3325 }
3326 }
3327
3328 {
3329 for (uint32_t count = 0; count < m_octree.getNumOctants(); count++){
3330 int8_t marker = m_octree.m_octants[count].getMarker();
3331 if (marker != 0){
3332 idx.push_back(count);
3333 markers.push_back(marker);
3334 isghost.push_back(false);
3335 }
3336 }
3337 }
3338
3339 //insert last ghost brothers if present
3340 {
3341 for (uint32_t id : m_octree.m_lastGhostBros){
3342 int8_t marker = m_octree.m_ghosts[id].getMarker();
3343 idx.push_back(id);
3344 markers.push_back(marker);
3345 isghost.push_back(true);
3346 }
3347 }
3348 };
3349
3356 bool
3357 ParaTree::isNodeOnOctant(const Octant *nodeOctant, uint8_t nodeIndex, const Octant *octant) const {
3358
3359 int dim = octant->getDim();
3360
3361 // Get the coordinates of the node
3362 std::array<uint32_t, 3> nodeCoords = nodeOctant->getLogicalNode(nodeIndex);
3363
3364 // Get minimum/maximum coordinates of the contant
3365 std::array<uint32_t, 3> minOctantCoords = octant->getLogicalNode(0);
3366 std::array<uint32_t, 3> maxOctantCoords = octant->getLogicalNode(3 + 4 * (dim - 2));
3367
3368 // Check if the node intersects the bounding box octant
3369 //
3370 // NOTE: since the octants are cubes, the bounding box coincides with
3371 // the octant.
3372 for (int i = 0; i < dim; ++i) {
3373 // I-th coordinate of the node
3374 uint32_t nodeCoord = nodeCoords[i];
3375
3376 // Minimum/maximum i-th coordinate of the octant
3377 uint32_t minOctantBBCoord = minOctantCoords[i];
3378 uint32_t maxOctantBBCoord = maxOctantCoords[i];
3379
3380 // Check if the node intersects the octant
3381 if (nodeCoord < minOctantBBCoord) {
3382 return false;
3383 } else if (nodeCoord > maxOctantBBCoord) {
3384 return false;
3385 }
3386 }
3387
3388 return true;
3389 }
3390
3397 bool
3398 ParaTree::isEdgeOnOctant(const Octant* edgeOctant, uint8_t edgeIndex, const Octant* octant) const {
3399
3400 // Edges are only defined on three-dimensional trees.
3401 int dim = octant->getDim();
3402 assert(dim == 3);
3403
3404 // Get the coordinates of the edge
3405 const uint8_t (*edgeNodes)[2] = &m_treeConstants->edgeNode[edgeIndex];
3406 std::array<uint32_t, 3> minEdgeCoords = edgeOctant->getLogicalNode((*edgeNodes)[0]);
3407 std::array<uint32_t, 3> maxEdgeCoords = edgeOctant->getLogicalNode((*edgeNodes)[1]);
3408
3409 // Get minimum/maximum coordinates of the contant
3410 std::array<uint32_t, 3> minOctantCoords = octant->getLogicalNode(0);
3411 std::array<uint32_t, 3> maxOctantCoords = octant->getLogicalNode(7);
3412
3413 // Check if the edge intersects the bounding box octant
3414 //
3415 // NOTE: since the octants are cubes, the bounding box coincides with
3416 // the octant.
3417 for (int i = 0; i < dim; ++i) {
3418 // Minimum/maximum i-th coordinate of the edge
3419 uint32_t minEdgeCoord = minEdgeCoords[i];
3420 uint32_t maxEdgeCoord = maxEdgeCoords[i];
3421
3422 // Minimum/maximum i-th coordinate of the octant
3423 uint32_t minOctantBBCoord = minOctantCoords[i];
3424 uint32_t maxOctantBBCoord = maxOctantCoords[i];
3425
3426 // Check if the edge intersects the octant
3427 if (minEdgeCoord < minOctantBBCoord && maxEdgeCoord < minOctantBBCoord) {
3428 return false;
3429 } else if (minEdgeCoord > maxOctantBBCoord && maxEdgeCoord > maxOctantBBCoord) {
3430 return false;
3431 }
3432 }
3433
3434 return true;
3435 }
3436
3443 bool
3444 ParaTree::isFaceOnOctant(const Octant* faceOctant, uint8_t faceIndex, const Octant* octant) const {
3445
3446 int dim = octant->getDim();
3447
3448 // Get minimum/maximum coordinates of the face
3449 const uint8_t (*faceNodes)[6][4] = &m_treeConstants->faceNode;
3450 std::array<uint32_t, 3> minFaceCoords = faceOctant->getLogicalNode((*faceNodes)[faceIndex][0]);
3451 std::array<uint32_t, 3> maxFaceCoords = faceOctant->getLogicalNode((*faceNodes)[faceIndex][2 * dim - 1]);
3452
3453 // Get minimum/maximum coordinates of the contant
3454 std::array<uint32_t, 3> minOctantCoords = octant->getLogicalNode(0);
3455 std::array<uint32_t, 3> maxOctantCoords = octant->getLogicalNode(3 + 4 * (dim - 2));
3456
3457 // Check if the face intersects the bounding box octant
3458 for (int i = 0; i < dim; ++i) {
3459 // Minimum/maximum i-th coordinate of the face
3460 uint32_t minFaceCoord = minFaceCoords[i];
3461 uint32_t maxFaceCoord = maxFaceCoords[i];
3462
3463 // Minimum/maximum i-th coordinate of the octant
3464 uint32_t minOctantBBCoord = minOctantCoords[i];
3465 uint32_t maxOctantBBCoord = maxOctantCoords[i];
3466
3467 // Check if the face intersects the octant
3468 if (minFaceCoord < minOctantBBCoord && maxFaceCoord < minOctantBBCoord) {
3469 return false;
3470 } else if (minFaceCoord > maxOctantBBCoord && maxFaceCoord > maxOctantBBCoord) {
3471 return false;
3472 }
3473 }
3474
3475 return true;
3476 }
3477
3478 // =================================================================================== //
3479 // OTHER PARATREE BASED METHODS //
3480 // =================================================================================== //
3481
3482
3490 void
3491 ParaTree::settleMarkers(){
3492
3493 (*m_log) << "---------------------------------------------" << endl;
3494 (*m_log) << " SETTLE MARKERS " << endl;
3495
3496 balance21(true, false);
3497
3498 (*m_log) << " " << endl;
3499 (*m_log) << "---------------------------------------------" << endl;
3500
3501 };
3502
3509 void
3510 ParaTree::preadapt(){
3511
3512 balance21(true, false);
3513
3514 m_lastOp = OP_PRE_ADAPT;
3515
3516 (*m_log) << "---------------------------------------------" << endl;
3517 (*m_log) << " PRE-ADAPT " << endl;
3518
3519 (*m_log) << " " << endl;
3520 (*m_log) << "---------------------------------------------" << endl;
3521
3522 };
3523
3527 bool
3528 ParaTree::checkToAdapt(){
3529
3530 bool lcheck = false;
3531 bool gcheck = false;
3532 octvector::iterator it = m_octree.m_octants.begin();
3533 octvector::iterator itend = m_octree.m_octants.end();
3534 while(!lcheck && it != itend){
3535 lcheck = (it->getMarker() != 0);
3536 ++it;
3537 }
3538 if (m_nproc == 1){
3539 gcheck = lcheck;
3540 }
3541 else{
3542#if BITPIT_ENABLE_MPI==1
3543 m_errorFlag = MPI_Allreduce(&lcheck,&gcheck,1,MPI_C_BOOL,MPI_LOR,m_comm);
3544#endif
3545 }
3546 return gcheck;
3547 };
3548
3557 bool
3558 ParaTree::adapt(bool mapper_flag){
3559
3560 bool done = false;
3561
3562 done = private_adapt_mapidx(mapper_flag);
3563 m_status += done;
3564 return done;
3565
3566 };
3567
3572 bool
3573 ParaTree::adaptGlobalRefine(bool mapper_flag) {
3574 //TODO recoding for adapting with abs(marker) > 1
3575 uint32_t nocts0 = getNumOctants();
3576 vector<Octant>::iterator iter, iterend = m_octree.m_octants.end();
3577
3578 for (iter = m_octree.m_octants.begin(); iter != iterend; ++iter){
3579 iter->m_info[Octant::INFO_NEW4REFINEMENT] = false;
3580 iter->m_info[Octant::INFO_NEW4COARSENING] = false;
3581 }
3582
3583 // Initialize mapping
3584 m_mapIdx.resize(nocts0);
3585 m_mapIdx.shrink_to_fit();
3586 for (uint32_t i=0; i<nocts0; i++){
3587 m_mapIdx[i] = i;
3588 }
3589
3590 // Update tree
3591 bool globalDone = false;
3592#if BITPIT_ENABLE_MPI==1
3593 if(m_serial){
3594#endif
3595 (*m_log) << "---------------------------------------------" << endl;
3596 (*m_log) << " ADAPT (Global Refine)" << endl;
3597 (*m_log) << " " << endl;
3598
3599 (*m_log) << " " << endl;
3600 (*m_log) << " Initial Number of octants : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
3601
3602 // Refine
3603 while(m_octree.globalRefine(m_mapIdx));
3604
3605 if (getNumOctants() > nocts0)
3606 globalDone = true;
3607 (*m_log) << " Number of octants after Refine : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
3608 updateAdapt();
3609
3610 (*m_log) << " " << endl;
3611 (*m_log) << "---------------------------------------------" << endl;
3612#if BITPIT_ENABLE_MPI==1
3613 }
3614 else{
3615 (*m_log) << "---------------------------------------------" << endl;
3616 (*m_log) << " ADAPT (Global Refine)" << endl;
3617 (*m_log) << " " << endl;
3618
3619 (*m_log) << " " << endl;
3620 (*m_log) << " Initial Number of octants : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
3621
3622 // Refine
3623 while(m_octree.globalRefine(m_mapIdx));
3624
3625 bool localDone = false;
3626 if (getNumOctants() > nocts0)
3627 localDone = true;
3628 updateAdapt();
3629 computeGhostHalo();
3630 (*m_log) << " Number of octants after Refine : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
3631
3632 m_errorFlag = MPI_Allreduce(&localDone,&globalDone,1,MPI_C_BOOL,MPI_LOR,m_comm);
3633 (*m_log) << " " << endl;
3634 (*m_log) << "---------------------------------------------" << endl;
3635 }
3636#endif
3637
3638 // Update last operation
3639 if (mapper_flag) {
3640 m_lastOp = OP_ADAPT_MAPPED;
3641 }
3642 else{
3643 m_lastOp = OP_ADAPT_UNMAPPED;
3644 }
3645
3646 return globalDone;
3647 }
3648
3653 bool
3654 ParaTree::adaptGlobalCoarse(bool mapper_flag) {
3655 //TODO recoding for adapting with abs(marker) > 1
3656 uint32_t nocts0 = getNumOctants();
3657 vector<Octant>::iterator iter, iterend = m_octree.m_octants.end();
3658
3659 for (iter = m_octree.m_octants.begin(); iter != iterend; ++iter){
3660 iter->m_info[Octant::INFO_NEW4REFINEMENT] = false;
3661 iter->m_info[Octant::INFO_NEW4COARSENING] = false;
3662 }
3663
3664 if (mapper_flag){
3665 m_mapIdx.resize(nocts0);
3666 for (uint32_t i=0; i<nocts0; i++){
3667 m_mapIdx[i] = i;
3668 }
3669 } else {
3670 m_mapIdx.clear();
3671 }
3672 m_mapIdx.shrink_to_fit();
3673
3674 bool globalDone = false;
3675#if BITPIT_ENABLE_MPI==1
3676 if(m_serial){
3677#endif
3678 (*m_log) << "---------------------------------------------" << endl;
3679 (*m_log) << " ADAPT (Global Coarse)" << endl;
3680 (*m_log) << " " << endl;
3681
3682 // 2:1 Balance
3683 balance21(true, false);
3684
3685 (*m_log) << " " << endl;
3686 (*m_log) << " Initial Number of octants : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
3687
3688 // Coarse
3689 while(m_octree.globalCoarse(m_mapIdx));
3690 updateAfterCoarse();
3691 balance21(false, true);
3692 while(m_octree.refine(m_mapIdx));
3693 updateAdapt();
3694
3695 if (getNumOctants() < nocts0){
3696 globalDone = true;
3697 }
3698
3699 (*m_log) << " Number of octants after Coarse : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
3700 (*m_log) << " " << endl;
3701 (*m_log) << "---------------------------------------------" << endl;
3702#if BITPIT_ENABLE_MPI==1
3703 }
3704 else{
3705 (*m_log) << "---------------------------------------------" << endl;
3706 (*m_log) << " ADAPT (Global Coarse)" << endl;
3707 (*m_log) << " " << endl;
3708
3709 // 2:1 Balance
3710 balance21(true, false);
3711
3712 (*m_log) << " " << endl;
3713 (*m_log) << " Initial Number of octants : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
3714
3715 // Coarse
3716 while(m_octree.globalCoarse(m_mapIdx));
3717 updateAfterCoarse();
3718 computeGhostHalo();
3719 balance21(false, true);
3720 while(m_octree.refine(m_mapIdx));
3721 updateAdapt();
3722
3723 computeGhostHalo();
3724 bool localDone = false;
3725 if (getNumOctants() < nocts0){
3726 localDone = true;
3727 }
3728
3729 m_errorFlag = MPI_Allreduce(&localDone,&globalDone,1,MPI_C_BOOL,MPI_LOR,m_comm);
3730 (*m_log) << " Number of octants after Coarse : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
3731 (*m_log) << " " << endl;
3732 (*m_log) << "---------------------------------------------" << endl;
3733 }
3734#endif
3735 return globalDone;
3736 }
3737
3742 int8_t
3743 ParaTree::getMaxDepth() const{
3744 return m_maxDepth;
3745 };
3746
3752 int
3753 ParaTree::findOwner(uint64_t morton) const {
3754 // Early return if the requested morton is on first partition
3755 if (morton <= m_partitionLastDesc[0]) {
3756 return 0;
3757 }
3758 if (morton > m_partitionLastDesc[m_nproc-1]) {
3759 return -1;
3760 }
3761
3762 // Find the partition using a bisection method
3763 int p = -1;
3764 int beg = 0;
3765 int end = m_nproc -1;
3766 int seed = m_nproc/2;
3767 while(beg != end){
3768 if(morton <= m_partitionLastDesc[seed]){
3769 end = seed;
3770 if(morton > m_partitionLastDesc[seed-1])
3771 beg = seed;
3772 }
3773 else{
3774 beg = seed;
3775 if(morton <= m_partitionLastDesc[seed+1])
3776 beg = seed + 1;
3777 }
3778 seed = beg + (end - beg) / 2;
3779 }
3780 if(beg!=0){
3781 while(m_partitionLastDesc[beg] == m_partitionLastDesc[beg-1]){
3782 --beg;
3783 if(beg==0)
3784 break;
3785 }
3786 }
3787 p = beg;
3788 return p;
3789 }
3790
3796 int
3797 ParaTree::getOwnerRank(uint64_t globalIndex) const {
3798 // Get the iterator point to the onwer rank
3799 std::vector<uint64_t>::const_iterator rankItr = std::lower_bound (m_partitionRangeGlobalIdx.begin(), m_partitionRangeGlobalIdx.end(), globalIndex);
3800
3801 // Get the onwer rank
3802 int ownerRank;
3803 if (rankItr == m_partitionRangeGlobalIdx.end()) {
3804 ownerRank = -1;
3805 } else {
3806 ownerRank = static_cast<int>(std::distance(m_partitionRangeGlobalIdx.begin(), rankItr));
3807 }
3808
3809 return ownerRank;
3810 }
3811
3814 void
3815 ParaTree::computeConnectivity() {
3816 m_octree.computeConnectivity();
3817 }
3818
3824 void
3825 ParaTree::clearConnectivity(bool release) {
3826 m_octree.clearConnectivity(release);
3827 }
3828
3831 void
3832 ParaTree::updateConnectivity() {
3833 m_octree.updateConnectivity();
3834 }
3835
3840 const u32vector2D &
3841 ParaTree::getConnectivity() const {
3842 return m_octree.m_connectivity;
3843 }
3844
3850 const u32vector &
3851 ParaTree::getConnectivity(uint32_t idx) const {
3852 return m_octree.m_connectivity[idx];
3853 }
3854
3859 const u32vector &
3860 ParaTree::getConnectivity(Octant* oct) const {
3861 return m_octree.m_connectivity[getIdx(oct)];
3862 }
3863
3867 const u32arr3vector &
3868 ParaTree::getNodes() const {
3869 return m_octree.m_nodes;
3870 }
3871
3876 const u32array3 &
3877 ParaTree::getNodeLogicalCoordinates(uint32_t node) const {
3878 return m_octree.m_nodes[node];
3879 }
3880
3885 darray3
3886 ParaTree::getNodeCoordinates(uint32_t node) const {
3887 return m_trans.mapCoordinates(m_octree.m_nodes[node]);
3888 }
3889
3894 const u32vector2D &
3895 ParaTree::getGhostConnectivity() const {
3896 return m_octree.m_ghostsConnectivity;
3897 }
3898
3904 const u32vector &
3905 ParaTree::getGhostConnectivity(uint32_t idx) const {
3906 return m_octree.m_ghostsConnectivity[idx];
3907 }
3908
3914 const u32vector &
3915 ParaTree::getGhostConnectivity(const Octant* oct) const {
3916 return m_octree.m_ghostsConnectivity[getIdx(oct)];
3917 }
3918
3919
3920
3924 bool
3925 ParaTree::check21Balance(){
3926
3927 bool balanced=true;
3928 uint32_t nocts = getNumOctants();
3929 uint8_t balanceCodim = getBalanceCodimension();
3930 uint8_t maxCodim = balanceCodim;
3931 u32vector neighs;
3932 bvector isghost;
3933 uint8_t levelDiff;
3934 for(uint8_t c = 1; c <= maxCodim; ++c){
3935 uint8_t nCodimSubelements = c==1 ? getNfaces() : (c==m_dim ? getNnodes() : getNedges());
3936 for(uint32_t i = 0; i < nocts; ++i){
3937 uint8_t level = getLevel(i);
3938 for(uint8_t f = 0; f < nCodimSubelements; ++f){
3939 findNeighbours(i,f,c,neighs,isghost);
3940 size_t nsize = neighs.size();
3941 for(size_t n = 0; n < nsize; ++n){
3942 const Octant* noct = isghost[n] ? getGhostOctant(neighs[n]) : getOctant(neighs[n]);
3943 uint8_t nlevel = noct->getLevel();
3944 levelDiff = uint8_t(abs(int(nlevel)-int(level)));
3945 if(levelDiff > 1){
3946 log::Visibility visi = m_log->getDefaultVisibility();
3947 m_log->setDefaultVisibility(log::VISIBILITY_GLOBAL);
3948 (*m_log) << "---------------------------------------------" << std::endl;
3949 (*m_log) << "LOCALLY 2:1 UNBALANCED OCTREE" << std::endl;
3950 (*m_log) << "I'm " << getRank() << ": element " << i << " is not 2:1 balanced across " << int(f) << " subentity of codim " << int(c) << ", relative to " << (isghost[n] ? "ghost" : "internal") << "neighbour " << neighs[n] << std::endl;
3951 (*m_log) << "---------------------------------------------" << std::endl;
3952 m_log->setDefaultVisibility(visi);
3953 balanced = false;
3954 break;
3955 }
3956 }
3957 if(!balanced)
3958 break;
3959 }
3960 if(!balanced)
3961 break;
3962 }
3963 if(!balanced)
3964 break;
3965 }
3966
3967 bool gBalanced;
3968#if BITPIT_ENABLE_MPI==1
3969 MPI_Allreduce(&balanced,&gBalanced,1,MPI_C_BOOL,MPI_LAND,getComm());
3970#else
3971 gBalanced = balanced;
3972#endif
3973
3974 if(gBalanced){
3975 (*m_log) << "---------------------------------------------" << std::endl;
3976 (*m_log) << "CORRECTLY GLOBAL 2:1 BALANCED OCTREE" << std::endl;
3977 (*m_log) << "---------------------------------------------" << std::endl;
3978 }
3979 else{
3980 (*m_log) << "---------------------------------------------" << std::endl;
3981 (*m_log) << "UNCORRECTLY GLOBAL 2:1 BALANCED OCTREE" << std::endl;
3982 (*m_log) << "---------------------------------------------" << std::endl;
3983 }
3984 return gBalanced;
3985 }
3986
3987
3988
3989#if BITPIT_ENABLE_MPI==1
3990
3996 void
3997 ParaTree::loadBalance(const dvector* weight){
3998
3999 //Write info on log
4000 (*m_log) << "---------------------------------------------" << endl;
4001 (*m_log) << " LOAD BALANCE " << endl;
4002
4003 m_lastOp = OP_LOADBALANCE;
4004 if (m_nproc>1){
4005
4006 std::vector<uint32_t> partition(m_nproc);
4007 if (weight == NULL)
4008 computePartition(partition.data());
4009 else
4010 computePartition(weight, partition.data());
4011
4012 weight = NULL;
4013
4014 privateLoadBalance<DummyDataLBImpl>(partition.data(), nullptr);
4015
4016 //Write info of final partition on log
4017 (*m_log) << " " << endl;
4018 (*m_log) << " Final Parallel partition : " << endl;
4019 (*m_log) << " Octants for proc "+ to_string(static_cast<unsigned long long>(0))+" : " + to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << endl;
4020 for(int ii=1; ii<m_nproc; ii++){
4021 (*m_log) << " Octants for proc "+ to_string(static_cast<unsigned long long>(ii))+" : " + to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[ii]-m_partitionRangeGlobalIdx[ii-1])) << endl;
4022 }
4023 (*m_log) << " " << endl;
4024 (*m_log) << "---------------------------------------------" << endl;
4025
4026 }
4027 else{
4028 (*m_log) << " " << endl;
4029 (*m_log) << " Serial partition : " << endl;
4030 (*m_log) << " Octants for proc "+ to_string(static_cast<unsigned long long>(0))+" : " + to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << endl;
4031 (*m_log) << " " << endl;
4032 (*m_log) << "---------------------------------------------" << endl;
4033 }
4034
4035 }
4036
4043 void
4044 ParaTree::loadBalance(uint8_t & level, const dvector* weight){
4045
4046 //Write info on log
4047 (*m_log) << "---------------------------------------------" << endl;
4048 (*m_log) << " LOAD BALANCE " << endl;
4049
4050 m_lastOp = OP_LOADBALANCE;
4051 if (m_nproc>1){
4052
4053 std::vector<uint32_t> partition(m_nproc);
4054 computePartition(level, weight, partition.data());
4055
4056 privateLoadBalance<DummyDataLBImpl>(partition.data(), nullptr);
4057
4058 //Write info of final partition on log
4059 (*m_log) << " " << endl;
4060 (*m_log) << " Final Parallel partition : " << endl;
4061 (*m_log) << " Octants for proc "+ to_string(static_cast<unsigned long long>(0))+" : " + to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << endl;
4062 for(int ii=1; ii<m_nproc; ii++){
4063 (*m_log) << " Octants for proc "+ to_string(static_cast<unsigned long long>(ii))+" : " + to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[ii]-m_partitionRangeGlobalIdx[ii-1])) << endl;
4064 }
4065 (*m_log) << " " << endl;
4066 (*m_log) << "---------------------------------------------" << endl;
4067
4068 }
4069 else{
4070 (*m_log) << " " << endl;
4071 (*m_log) << " Serial partition : " << endl;
4072 (*m_log) << " Octants for proc "+ to_string(static_cast<unsigned long long>(0))+" : " + to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << endl;
4073 (*m_log) << " " << endl;
4074 (*m_log) << "---------------------------------------------" << endl;
4075 }
4076
4077 };
4078
4088 ParaTree::LoadBalanceRanges
4089 ParaTree::evalLoadBalanceRanges(dvector *weights){
4090
4091 // If there is only one process no octants can be exchanged
4092 if (m_nproc == 1) {
4093 LoadBalanceRanges loadBalanceInfo;
4094 loadBalanceInfo.sendAction = LoadBalanceRanges::ACTION_NONE;
4095 loadBalanceInfo.recvAction = LoadBalanceRanges::ACTION_NONE;
4096
4097 return loadBalanceInfo;
4098 }
4099
4100 // Compute updated partition
4101 std::vector<uint32_t> updatedPartition(m_nproc);
4102 if (weights) {
4103 computePartition(weights, updatedPartition.data());
4104 } else {
4105 computePartition(updatedPartition.data());
4106 }
4107
4108 // Evaluate send ranges
4109 return evalLoadBalanceRanges(updatedPartition.data());
4110 }
4111
4126 ParaTree::LoadBalanceRanges
4127 ParaTree::evalLoadBalanceRanges(uint8_t level, dvector *weights){
4128
4129 // If there is only one process no octants can be sent
4130 if (m_nproc == 1) {
4131 LoadBalanceRanges loadBalanceInfo;
4132 loadBalanceInfo.sendAction = LoadBalanceRanges::ACTION_NONE;
4133 loadBalanceInfo.recvAction = LoadBalanceRanges::ACTION_NONE;
4134
4135 return loadBalanceInfo;
4136 }
4137
4138 // Compute updated partition
4139 std::vector<uint32_t> updatedPartition(m_nproc);
4140 computePartition(level, weights, updatedPartition.data());
4141
4142 // Evaluate send ranges
4143 return evalLoadBalanceRanges(updatedPartition.data());
4144 }
4145
4154 ParaTree::LoadBalanceRanges
4155 ParaTree::evalLoadBalanceRanges(const uint32_t *updatedPartition){
4156
4157 // If there is only one process no octants can be sent
4158 if (m_nproc == 1) {
4159 LoadBalanceRanges loadBalanceInfo;
4160 loadBalanceInfo.sendAction = LoadBalanceRanges::ACTION_NONE;
4161 loadBalanceInfo.recvAction = LoadBalanceRanges::ACTION_NONE;
4162
4163 return loadBalanceInfo;
4164 }
4165
4166 ExchangeRanges sendRanges = evalLoadBalanceSendRanges(updatedPartition);
4167 ExchangeRanges recvRanges = evalLoadBalanceRecvRanges(updatedPartition);
4168
4169 return LoadBalanceRanges(m_serial, std::move(sendRanges), std::move(recvRanges));
4170 }
4171
4179 ParaTree::ExchangeRanges
4180 ParaTree::evalLoadBalanceSendRanges(const uint32_t *updatedPartition){
4181
4182 ExchangeRanges sendRanges;
4183
4184 // If there is only one process no octants can be sent
4185 if (m_nproc == 1) {
4186 return sendRanges;
4187 }
4188
4189 // Compute current partition schema
4190 std::vector<uint32_t> currentPartition(m_nproc, 0);
4191 if (!m_serial) {
4192 currentPartition[0] = m_partitionRangeGlobalIdx[0] + 1;
4193 for (int i = 1; i < m_nproc; ++i) {
4194 currentPartition[i] = m_partitionRangeGlobalIdx[i] - m_partitionRangeGlobalIdx[i - 1];
4195 }
4196 } else {
4197 currentPartition[m_rank] = getNumOctants();
4198 }
4199
4200 // Get the intersections
4201 PartitionIntersections globalIntersections = evalPartitionIntersections(currentPartition.data(), m_rank, updatedPartition);
4202
4203 // Evaluate the send local indexes
4204 uint64_t offset = 0;
4205 for (int i = 0; i < m_rank; ++i) {
4206 offset += currentPartition[i];
4207 }
4208
4209 for (const auto &intersectionEntry : globalIntersections) {
4210 int rank = intersectionEntry.first;
4211 if (rank == m_rank) {
4212 continue;
4213 }
4214
4215 const std::array<uint64_t, 2> &intersection = intersectionEntry.second;
4216
4217 std::array<uint32_t, 2> &sendRange = sendRanges[rank];
4218 sendRange[0] = intersection[0] - offset;
4219 sendRange[1] = intersection[1] - offset;
4220 }
4221
4222 return sendRanges;
4223 }
4224
4233 ParaTree::ExchangeRanges
4234 ParaTree::evalLoadBalanceRecvRanges(const uint32_t *updatedPartition){
4235
4236 ExchangeRanges recvRanges;
4237
4238 // If there is only one process no octants can be received
4239 if (m_nproc == 1) {
4240 return recvRanges;
4241 }
4242
4243 // Compute current partition schema
4244 std::vector<uint32_t> currentPartition(m_nproc, 0);
4245 if (!m_serial) {
4246 currentPartition[0] = m_partitionRangeGlobalIdx[0] + 1;
4247 for (int i = 1; i < m_nproc; ++i) {
4248 currentPartition[i] = m_partitionRangeGlobalIdx[i] - m_partitionRangeGlobalIdx[i - 1];
4249 }
4250 } else {
4251 currentPartition[m_rank] = getNumOctants();
4252 }
4253
4254 // Get the intersections
4255 PartitionIntersections globalIntersections = evalPartitionIntersections(updatedPartition, m_rank, currentPartition.data());
4256
4257 // Evaluate the receive local indexes
4258 uint64_t offset = 0;
4259 for (int i = 0; i < m_rank; ++i) {
4260 offset += updatedPartition[i];
4261 }
4262
4263 for (const auto &intersectionEntry : globalIntersections) {
4264 int rank = intersectionEntry.first;
4265 if (rank == m_rank) {
4266 continue;
4267 }
4268
4269 const std::array<uint64_t, 2> &intersection = intersectionEntry.second;
4270
4271 std::array<uint32_t, 2> &recvRange = recvRanges[rank];
4272 recvRange[0] = intersection[0] - offset;
4273 recvRange[1] = intersection[1] - offset;
4274 }
4275
4276 return recvRanges;
4277 }
4278
4296 ParaTree::PartitionIntersections
4297 ParaTree::evalPartitionIntersections(const uint32_t *partition_A, int rank_A, const uint32_t *partition_B){
4298
4299 PartitionIntersections intersections;
4300
4301 // If the partition is empty there are no intersections.
4302 if (partition_A[rank_A] == 0) {
4303 return intersections;
4304 }
4305
4306 // Calculate partition offsets
4307 std::vector<uint64_t> offsets_A(m_nproc + 1, 0);
4308 std::vector<uint64_t> offsets_B(m_nproc + 1, 0);
4309 for (int i = 0; i < m_nproc; ++i) {
4310 offsets_A[i + 1] = offsets_A[i] + partition_A[i];
4311 offsets_B[i + 1] = offsets_B[i] + partition_B[i];
4312 }
4313
4314 uint64_t beginGlobalId_A = offsets_A[m_rank];
4315 uint64_t endGlobalId_A = offsets_A[m_rank + 1];
4316
4317 auto firstRankItr = std::upper_bound(offsets_B.begin(), offsets_B.end(), beginGlobalId_A);
4318 assert(firstRankItr != offsets_B.begin());
4319 firstRankItr--;
4320
4321 for (auto itr = firstRankItr; itr != offsets_B.end(); ++itr) {
4322 int rank_B = static_cast<int>(std::distance(offsets_B.begin(), itr));
4323 uint64_t beginGlobalId_B = offsets_B[rank_B];
4324 uint64_t endGlobalId_B = offsets_B[rank_B + 1];
4325
4326 std::array<uint64_t, 2> &intersection = intersections[rank_B];
4327 intersection[0] = std::max(beginGlobalId_A, beginGlobalId_B);
4328 intersection[1] = std::min(endGlobalId_A, endGlobalId_B);
4329
4330 if (endGlobalId_B >= endGlobalId_A) {
4331 break;
4332 }
4333 }
4334
4335 return intersections;
4336 }
4337#endif
4338
4343 double
4344 ParaTree::levelToSize(uint8_t level) const {
4345 uint32_t size = uint32_t(1)<<(m_treeConstants->maxLevel-level);
4346 return m_trans.mapSize(size);
4347 }
4348
4349 // =================================================================================== //
4350 // OTHER INTERSECTION BASED METHODS //
4351 // =================================================================================== //
4352
4355 void
4356 ParaTree::computeIntersections(){
4357 m_octree.computeIntersections();
4358 }
4359
4360 // =================================================================================== //
4361 // OTHER PRIVATE METHODS //
4362 // =================================================================================== //
4363
4368 Octant&
4369 ParaTree::extractOctant(uint32_t idx) {
4370 return m_octree.extractOctant(idx) ;
4371 };
4372
4376 bool
4377 ParaTree::private_adapt_mapidx(bool mapflag) {
4378 //TODO recoding for adapting with abs(marker) > 1
4379
4380 m_loadBalanceRanges.clear();
4381 uint32_t nocts0 = getNumOctants();
4382 vector<Octant >::iterator iter, iterend = m_octree.m_octants.end();
4383
4384 for (iter = m_octree.m_octants.begin(); iter != iterend; ++iter){
4385 iter->m_info[Octant::INFO_NEW4REFINEMENT] = false;
4386 iter->m_info[Octant::INFO_NEW4COARSENING] = false;
4387 }
4388
4389 // m_mapIdx init
4390 if (mapflag) {
4391 m_mapIdx.resize(nocts0);
4392 for (uint32_t i=0; i<nocts0; i++){
4393 m_mapIdx[i] = i;
4394 }
4395 } else {
4396 m_mapIdx.clear();
4397 }
4398 m_mapIdx.shrink_to_fit();
4399
4400 bool globalDone = false;
4401#if BITPIT_ENABLE_MPI==1
4402 if(m_serial){
4403#endif
4404 (*m_log) << "---------------------------------------------" << endl;
4405 (*m_log) << " ADAPT (Refine/Coarse)" << endl;
4406 (*m_log) << " " << endl;
4407
4408 // 2:1 Balance
4409 if (m_lastOp != OP_PRE_ADAPT) {
4410 balance21(true, false);
4411 }
4412
4413 (*m_log) << " " << endl;
4414 (*m_log) << " Initial Number of octants : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
4415
4416 // Refine
4417 while(m_octree.refine(m_mapIdx));
4418 if (getNumOctants() > nocts0)
4419 globalDone = true;
4420 (*m_log) << " Number of octants after Refine : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
4421 nocts0 = getNumOctants();
4422 updateAdapt();
4423
4424 // Coarse
4425 while(m_octree.coarse(m_mapIdx));
4426 updateAfterCoarse();
4427 if (getNumOctants() < nocts0){
4428 globalDone = true;
4429 }
4430
4431 (*m_log) << " Number of octants after Coarse : " + to_string(static_cast<unsigned long long>(getNumOctants())) << endl;
4432 (*m_log) << " " << endl;
4433 (*m_log) << "---------------------------------------------" << endl;
4434#if BITPIT_ENABLE_MPI==1
4435 }
4436 else{
4437 (*m_log) << "---------------------------------------------" << endl;
4438 (*m_log) << " ADAPT (Refine/Coarse)" << endl;
4439 (*m_log) << " " << endl;
4440
4441 // 2:1 Balance
4442 if (m_lastOp != OP_PRE_ADAPT) {
4443 balance21(true, false);
4444 }
4445
4446 (*m_log) << " " << endl;
4447 (*m_log) << " Initial Number of octants : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
4448
4449 // Refine
4450 while(m_octree.refine(m_mapIdx));
4451 bool localDone = false;
4452 if (getNumOctants() > nocts0)
4453 localDone = true;
4454 updateAdapt();
4455 computeGhostHalo();
4456 (*m_log) << " Number of octants after Refine : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
4457 nocts0 = getNumOctants();
4458
4459
4460 // Coarse
4461 while(m_octree.coarse(m_mapIdx));
4462 updateAfterCoarse();
4463 computeGhostHalo();
4464 if (getNumOctants() < nocts0){
4465 localDone = true;
4466 }
4467
4468 m_errorFlag = MPI_Allreduce(&localDone,&globalDone,1,MPI_C_BOOL,MPI_LOR,m_comm);
4469 (*m_log) << " Number of octants after Coarse : " + to_string(static_cast<unsigned long long>(m_globalNumOctants)) << endl;
4470 (*m_log) << " " << endl;
4471 (*m_log) << "---------------------------------------------" << endl;
4472 }
4473#endif
4474
4475 // Update last operation
4476 if (mapflag) {
4477 m_lastOp = OP_ADAPT_MAPPED;
4478 }
4479 else{
4480 m_lastOp = OP_ADAPT_UNMAPPED;
4481 }
4482
4483 return globalDone;
4484 }
4485
4488 void
4489 ParaTree::updateAdapt(){
4490#if BITPIT_ENABLE_MPI==1
4491 if(m_serial)
4492 {
4493#endif
4494 for (int iproc=0; iproc<m_nproc; iproc++){
4495 m_partitionRangeGlobalIdx0[iproc] = m_partitionRangeGlobalIdx[iproc];
4496 }
4497 m_maxDepth = m_octree.m_localMaxDepth;
4498 m_globalNumOctants = getNumOctants();
4499 for(int p = 0; p < m_nproc; ++p){
4500 m_partitionRangeGlobalIdx[p] = m_globalNumOctants - 1;
4501 }
4502 m_internals.resize(getNumOctants());
4503 int i = 0;
4504 octvector::iterator itend = m_octree.m_octants.end();
4505 for (octvector::iterator it = m_octree.m_octants.begin(); it != itend; ++it){
4506 m_internals[i] = &(*it);
4507 i++;
4508 }
4509#if BITPIT_ENABLE_MPI==1
4510 }
4511 else
4512 {
4513 for (int iproc=0; iproc<m_nproc; iproc++){
4514 m_partitionRangeGlobalIdx0[iproc] = m_partitionRangeGlobalIdx[iproc];
4515 }
4516 //update m_maxDepth
4517 m_errorFlag = MPI_Allreduce(&m_octree.m_localMaxDepth,&m_maxDepth,1,MPI_INT8_T,MPI_MAX,m_comm);
4518 //update m_globalNumOctants
4519 uint64_t local_num_octants = getNumOctants();
4520 m_errorFlag = MPI_Allreduce(&local_num_octants,&m_globalNumOctants,1,MPI_UINT64_T,MPI_SUM,m_comm);
4521 //update m_partitionRangeGlobalIdx
4522 uint64_t* rbuff = new uint64_t[m_nproc];
4523 m_errorFlag = MPI_Allgather(&local_num_octants,1,MPI_UINT64_T,rbuff,1,MPI_UINT64_T,m_comm);
4524 for(int p = 0; p < m_nproc; ++p){
4525 m_partitionRangeGlobalIdx[p] = 0;
4526 for(int pp = 0; pp <=p; ++pp)
4527 m_partitionRangeGlobalIdx[p] += rbuff[pp];
4528 --m_partitionRangeGlobalIdx[p];
4529 }
4530 delete [] rbuff; rbuff = NULL;
4531 }
4532#endif
4533 }
4534
4535#if BITPIT_ENABLE_MPI==1
4541 void
4542 ParaTree::computePartition(uint32_t* partition){
4543
4544 uint32_t division_result = 0;
4545 uint32_t remind = 0;
4546
4547 division_result = uint32_t(m_globalNumOctants/(uint64_t)m_nproc);
4548 remind = (uint32_t)(m_globalNumOctants%(uint64_t)m_nproc);
4549
4550 for(uint32_t i = 0; i < (uint32_t)m_nproc; ++i)
4551 if(i<remind)
4552 partition[i] = division_result + 1;
4553 else
4554 partition[i] = division_result;
4555
4556 }
4557
4564 void
4565 ParaTree::computePartition(const dvector *weight, uint32_t *partition){
4566
4567 assert(weight->size() >= m_octree.getNumOctants());
4568
4569 // Evaluate global weights
4570 //
4571 // If the tree is serial, all process have all the octants, hence
4572 // global weights and local weights are the same.
4573 std::vector<double> globalWeightsStorage;
4574
4575 const double *globalWeights;
4576 if (m_serial) {
4577 globalWeights = weight->data();
4578 } else {
4579 // Information about current partitioning and displacements should
4580 // be stored using uint32_t, however the maximum number of items
4581 // that can be exchanged using MPI funcitons is limited to INT_MAX.
4582 assert(m_globalNumOctants <= INT_MAX);
4583 std::vector<int> currentPartition(m_nproc);
4584 int nOctants = m_octree.getNumOctants();
4585 MPI_Allgather(&nOctants, 1, MPI_INT, currentPartition.data(), 1, MPI_INT, m_comm);
4586
4587 std::vector<int> displacements(m_nproc);
4588 displacements[0] = 0;
4589 for(int i = 1; i < m_nproc; ++i){
4590 displacements[i] = displacements[i - 1] + currentPartition[i - 1];
4591 }
4592
4593 globalWeightsStorage.resize(m_globalNumOctants);
4594 globalWeights = globalWeightsStorage.data();
4595 MPI_Allgatherv(weight->data(), nOctants, MPI_DOUBLE, globalWeightsStorage.data(),
4596 currentPartition.data(), displacements.data(), MPI_DOUBLE, m_comm);
4597 }
4598
4599 // Initialize partitioning
4600 for (int i = 0; i < m_nproc; ++i) {
4601 partition[i] = 0;
4602 }
4603
4604 // Assign octants to partitions
4605 //
4606 // After evaluationg the target weight of a partition, octants will
4607 // be added to that partition until the weigth of the partition is
4608 // greater or equal the target weigth or until all the octant are
4609 // assigned
4610 uint32_t nAssigendOctants = 0;
4611 for (int i = 0; i < m_nproc - 1; ++i) {
4612 double unassignedWeight = 0.;
4613 for (uint32_t n=nAssigendOctants; n<m_globalNumOctants; n++){
4614 unassignedWeight += globalWeights[n];
4615 }
4616 double targetWeight = unassignedWeight / (m_nproc - i);
4617
4618 double partitionWeight = 0.;
4619 while(partitionWeight < targetWeight){
4620 partitionWeight += globalWeights[nAssigendOctants];
4621 partition[i]++;
4622
4623 nAssigendOctants++;
4624 if (nAssigendOctants == m_globalNumOctants) {
4625 break;
4626 }
4627 }
4628
4629 if (nAssigendOctants == m_globalNumOctants) {
4630 break;
4631 }
4632 }
4633 partition[m_nproc-1] = static_cast<uint32_t>(m_globalNumOctants - nAssigendOctants);
4634 };
4635
4648 void
4649 ParaTree::computePartition(uint8_t level_, const dvector *weight, uint32_t *partition) {
4650
4651 // Compute partitioning without family constrains
4652 uint32_t* partition_temp = new uint32_t[m_nproc];
4653 if (weight==NULL){
4654 computePartition(partition_temp);
4655 }
4656 else{
4657 computePartition(weight, partition_temp);
4658 }
4659
4660 // Modify partitioning to take into account family constrains
4661 //
4662 // Partitioning is modified to guarantee that families of octants at
4663 // the desired level above the maximum depth reached in the tree
4664 // are retained compact on the same process.
4665 uint8_t level = uint8_t(min(int(max(int(m_maxDepth) - int(level_), int(1))) , int(m_treeConstants->maxLevel)));
4666 uint8_t* new_boundary_owner = new uint8_t[m_nproc-1];
4667 uint8_t new_interfaces_count;
4668 uint8_t first_new_interface_rank_index;
4669 uint8_t* new_interfaces_count_per_rank = new uint8_t[m_nproc];
4670 uint8_t* first_new_interface_rank_index_per_rank = new uint8_t[m_nproc];
4671
4672 uint32_t Dh = uint32_t(pow(double(2),double(m_treeConstants->maxLevel-level)));
4673 uint32_t istart, nocts, rest;
4674 uint32_t forw = 0, backw = 0;
4675 uint32_t i = 0, iproc, j;
4676 uint64_t sum;
4677 int32_t* pointercomm;
4678 int32_t* deplace = new int32_t[m_nproc-1];
4679
4680 // Find processes currently owning the new incoming process boundaries
4681 // boundary_proc[i] = j means that the new process interface between i-th and (i+1)-th
4682 // processes is falling currently on j-th process
4683 j = 0;
4684 sum = 0;
4685 for (iproc=0; iproc<(uint32_t)(m_nproc-1); iproc++){
4686 sum += partition_temp[iproc];
4687 while(sum > m_partitionRangeGlobalIdx[j]){
4688 j++;
4689 }
4690 new_boundary_owner[iproc] = j;
4691 }
4692
4693 nocts = getNumOctants();
4694 sum = 0;
4695
4696 // Store how many process interfaces fall in the current rank. For these interfaces
4697 // the current rank has to communicate the info about the correction to the partition
4698 // strcture aimed to maintain the families compact.
4699 new_interfaces_count = 0;
4700
4701 // Store the index of the first process owner of the new process interface
4702 // It will be the starting index for the communication of the corrections on
4703 // the partition structure.
4704 first_new_interface_rank_index = 0;
4705
4706 for (iproc=0; iproc<(uint32_t)(m_nproc-1); iproc++){
4707 deplace[iproc] = 1;
4708 sum += partition_temp[iproc];
4709
4710 // If the current process owns a new process interface, check if
4711 // the family at the interface is compact and store the correction on the
4712 // temporary partition structure.
4713 if (new_boundary_owner[iproc] == m_rank){
4714 if (new_interfaces_count == 0){
4715 first_new_interface_rank_index = iproc;
4716 }
4717 new_interfaces_count++;
4718
4719 // Place istart at index of the last octant at new incoming process interface
4720 if (m_rank!=0)
4721 istart = static_cast<uint32_t>(sum - m_partitionRangeGlobalIdx[m_rank-1] - 1);
4722 else
4723 istart = static_cast<uint32_t>(sum);
4724
4725 // Compute the rest w.r.t. the size of the octant of the same size
4726 // of a compact family of the desired level
4727 i = istart;
4728 rest = m_octree.m_octants[i].getLogicalCoordinates(0)%Dh + m_octree.m_octants[i].getLogicalCoordinates(1)%Dh;
4729 if (m_dim == 3){
4730 rest += m_octree.m_octants[i].getLogicalCoordinates(2)%Dh;
4731 }
4732
4733 // Deplace forward the index of the octant i defining the process boundary until the previous
4734 // defined rest is zero, i.e. the families are compact up to the target level.
4735 while(rest!=0){
4736 if (i==nocts-1){
4737 i = istart + nocts;
4738 break;
4739 }
4740 i++;
4741 rest = m_octree.m_octants[i].getLogicalCoordinates(0)%Dh + m_octree.m_octants[i].getLogicalCoordinates(1)%Dh;
4742 if (m_dim == 3){
4743 rest += m_octree.m_octants[i].getLogicalCoordinates(2)%Dh;
4744 }
4745 }
4746
4747 // Store the computed forward correction (= local number of octants if not found)
4748 forw = i - istart;
4749
4750 // Do the same for a backward correction try
4751 i = istart;
4752 rest = m_octree.m_octants[i].getLogicalCoordinates(0)%Dh + m_octree.m_octants[i].getLogicalCoordinates(1)%Dh;
4753 if (m_dim == 3){
4754 rest += m_octree.m_octants[i].getLogicalCoordinates(2)%Dh;
4755 }
4756 while(rest!=0){
4757 if (i==0){
4758 i = istart - nocts;
4759 break;
4760 }
4761 i--;
4762 rest = m_octree.m_octants[i].getLogicalCoordinates(0)%Dh + m_octree.m_octants[i].getLogicalCoordinates(1)%Dh;
4763 if (m_dim == 3){
4764 rest += m_octree.m_octants[i].getLogicalCoordinates(2)%Dh;
4765 }
4766 }
4767
4768 // Store the backward correction previously computed
4769 backw = istart - i;
4770
4771 // If no correction are needed forward and backward corrections are zero,
4772 // so the deplace term will be zero and, even if communicated it will no correct
4773 // the partition structure.
4774 // If correction is needed the lower between forward and backward corrections
4775 // is the correct one, beeing the other one overdimensioned with local number of octants value.
4776 if (forw<backw)
4777 deplace[iproc] = forw;
4778 else
4779 deplace[iproc] = -(int32_t)backw;
4780 }
4781 }
4782
4783 // Communicate the right corrections to other processes
4784 m_errorFlag = MPI_Allgather(&new_interfaces_count,1,MPI_UINT8_T,new_interfaces_count_per_rank,1,MPI_UINT8_T,m_comm);
4785 m_errorFlag = MPI_Allgather(&first_new_interface_rank_index,1,MPI_UINT8_T,first_new_interface_rank_index_per_rank,1,MPI_UINT8_T,m_comm);
4786 for (iproc=0; iproc<(uint32_t)(m_nproc); iproc++){
4787 pointercomm = &deplace[first_new_interface_rank_index_per_rank[iproc]];
4788 m_errorFlag = MPI_Bcast(pointercomm, new_interfaces_count_per_rank[iproc], MPI_INT32_T, iproc, m_comm);
4789 }
4790
4791 // Apply the corrections stored in deplace container to the termporary
4792 // computed partition structure.
4793 // The new number of octants for each processors are increased by the
4794 // (signed) deplace value related to itself and decreased by the (signed)
4795 // deplace value of the previous process.
4796 for (iproc=0; iproc<(uint32_t)(m_nproc); iproc++){
4797 if (iproc < (uint32_t)(m_nproc-1))
4798 partition[iproc] = partition_temp[iproc] + deplace[iproc];
4799 else
4800 partition[iproc] = partition_temp[iproc];
4801 if (iproc !=0)
4802 partition[iproc] = partition[iproc] - deplace[iproc-1];
4803 }
4804
4805 delete [] partition_temp; partition_temp = NULL;
4806 delete [] new_boundary_owner; new_boundary_owner = NULL;
4807 delete [] new_interfaces_count_per_rank; new_interfaces_count_per_rank = NULL;
4808 delete [] first_new_interface_rank_index_per_rank; first_new_interface_rank_index_per_rank = NULL;
4809 delete [] deplace; deplace = NULL;
4810 }
4811
4814 void
4815 ParaTree::updateLoadBalance() {
4816 //update sizes
4817 m_octree.updateLocalMaxDepth();
4818 uint64_t* rbuff = new uint64_t[m_nproc];
4819
4820 uint64_t local_num_octants = getNumOctants();
4821 m_errorFlag = MPI_Allgather(&local_num_octants,1,MPI_UINT64_T,rbuff,1,MPI_UINT64_T,m_comm);
4822 for (int iproc=0; iproc<m_nproc; iproc++){
4823 m_partitionRangeGlobalIdx0[iproc] = m_partitionRangeGlobalIdx[iproc];
4824 }
4825 for(int p = 0; p < m_nproc; ++p){
4826 m_partitionRangeGlobalIdx[p] = 0;
4827 for(int pp = 0; pp <=p; ++pp)
4828 m_partitionRangeGlobalIdx[p] += rbuff[pp];
4829 --m_partitionRangeGlobalIdx[p];
4830 }
4831
4832 m_serial = false;
4833
4834 delete [] rbuff; rbuff = NULL;
4835
4836 //update first and last descendant
4837 m_octree.setFirstDescMorton();
4838 m_octree.setLastDescMorton();
4839 updateGlobalFirstDescMorton();
4840 updateGlobalLasttDescMorton();
4841 }
4842
4846 void
4847 ParaTree::setPboundGhosts() {
4848 //BUILD BORDER OCTANT INDECES VECTOR (map value) TO BE SENT TO THE RIGHT PROCESS (map key)
4849 //find local octants to be sent as ghost to the right processes
4850 //it visits the local octants building virtual neighbors on each octant face
4851 //find the owner of these virtual neighbor and build a map (process,border octants)
4852 //this map contains the local octants as ghosts for neighbor processes
4853
4854 // NO PBORDERS !
4855 //
4856 // TODO: provide an estimate of the border octants in order to reserve
4857 // the vectors that will contain them.
4858 m_bordersPerProc.clear();
4859 m_internals.resize(getNumOctants());
4860 m_pborders.resize(getNumOctants());
4861
4862 int countpbd = 0;
4863 int countint = 0;
4864 std::set<int> neighProcs;
4865 std::vector<std::array<int64_t, 3>> virtualNeighOffsets;
4866 for (uint32_t idx = 0; idx < m_octree.getNumOctants(); ++idx) {
4867 neighProcs.clear();
4868 Octant &octant = m_octree.m_octants[idx];
4869 u32array3 octantCoord = octant.getLogicalCoordinates();
4870 uint8_t octantLevel = octant.getLevel();
4871
4872 // Virtual Face Neighbors
4873 uint8_t maxFaceNeighLevel = std::min(m_octree.getMaxNeighLevel(octant), static_cast<uint8_t>(m_maxDepth));
4874 for(uint8_t i = 0; i < m_treeConstants->nFaces; ++i){
4875 bool isFacePeriodic = m_octree.isPeriodic(&octant, i);
4876 if (!isFacePeriodic) {
4877 bool isFaceBoundary = octant.getBound(i);
4878 if (isFaceBoundary) {
4879 continue;
4880 }
4881 }
4882
4883 std::array<uint32_t, 3> virtualOctantOrigin = octantCoord;
4884 if (isFacePeriodic) {
4885 std::array<int64_t, 3> periodicOffset = m_octree.getPeriodicOffset(octant, i);
4886 virtualOctantOrigin[0] = static_cast<uint32_t>(virtualOctantOrigin[0] + periodicOffset[0]);
4887 virtualOctantOrigin[1] = static_cast<uint32_t>(virtualOctantOrigin[1] + periodicOffset[1]);
4888 virtualOctantOrigin[2] = static_cast<uint32_t>(virtualOctantOrigin[2] + periodicOffset[2]);
4889 }
4890
4891 m_octree.computeVirtualNeighOffsets(octantLevel, i, maxFaceNeighLevel, &virtualNeighOffsets);
4892
4893 bool isFacePbound = false;
4894 for(const std::array<int64_t, 3> &virtualNeighOffset : virtualNeighOffsets){
4895 std::array<uint32_t, 3> virtualNeighCoord = virtualOctantOrigin;
4896 virtualNeighCoord[0] = static_cast<uint32_t>(virtualNeighCoord[0] + virtualNeighOffset[0]);
4897 virtualNeighCoord[1] = static_cast<uint32_t>(virtualNeighCoord[1] + virtualNeighOffset[1]);
4898 virtualNeighCoord[2] = static_cast<uint32_t>(virtualNeighCoord[2] + virtualNeighOffset[2]);
4899 uint64_t virtualNeighMorton = PABLO::computeMorton(m_dim, virtualNeighCoord[0], virtualNeighCoord[1], virtualNeighCoord[2]);
4900
4901 int virtualNeighProc = findOwner(virtualNeighMorton);
4902 if (virtualNeighProc != m_rank) {
4903 neighProcs.insert(virtualNeighProc);
4904 isFacePbound = true;
4905 }
4906 }
4907
4908 octant.setPbound(i, isFacePbound);
4909 }
4910
4911 // Virtual Edge Neighbors
4912 uint8_t maxEdgeNeighLevel = std::min(m_octree.getMaxEdgeNeighLevel(octant), static_cast<uint8_t>(m_maxDepth));
4913 for(uint8_t e = 0; e < m_treeConstants->nEdges; ++e){
4914 bool isEdgePeriodic = m_octree.isEdgePeriodic(&octant, e);
4915 if (!isEdgePeriodic) {
4916 bool isEdgeBoundary = octant.getEdgeBound(e);
4917 if (isEdgeBoundary) {
4918 continue;
4919 }
4920 }
4921
4922 std::array<uint32_t, 3> virtualOctantOrigin = octantCoord;
4923 if (isEdgePeriodic) {
4924 std::array<int64_t, 3> periodicOffset = m_octree.getEdgePeriodicOffset(octant, e);
4925 virtualOctantOrigin[0] = static_cast<uint32_t>(virtualOctantOrigin[0] + periodicOffset[0]);
4926 virtualOctantOrigin[1] = static_cast<uint32_t>(virtualOctantOrigin[1] + periodicOffset[1]);
4927 virtualOctantOrigin[2] = static_cast<uint32_t>(virtualOctantOrigin[2] + periodicOffset[2]);
4928 }
4929
4930 m_octree.computeVirtualEdgeNeighOffsets(octantLevel, e, maxEdgeNeighLevel, &virtualNeighOffsets);
4931
4932 for(const std::array<int64_t, 3> &virtualNeighOffset : virtualNeighOffsets){
4933 std::array<uint32_t, 3> virtualNeighCoord = virtualOctantOrigin;
4934 virtualNeighCoord[0] = static_cast<uint32_t>(virtualNeighCoord[0] + virtualNeighOffset[0]);
4935 virtualNeighCoord[1] = static_cast<uint32_t>(virtualNeighCoord[1] + virtualNeighOffset[1]);
4936 virtualNeighCoord[2] = static_cast<uint32_t>(virtualNeighCoord[2] + virtualNeighOffset[2]);
4937 uint64_t virtualNeighMorton = PABLO::computeMorton(m_dim, virtualNeighCoord[0], virtualNeighCoord[1], virtualNeighCoord[2]);
4938
4939 int virtualNeighProc = findOwner(virtualNeighMorton);
4940 if (virtualNeighProc != m_rank) {
4941 neighProcs.insert(virtualNeighProc);
4942 }
4943 }
4944 }
4945
4946 // Virtual Corner Neighbors
4947 uint8_t maxNodeNeighLevel = std::min(m_octree.getMaxNodeNeighLevel(octant), static_cast<uint8_t>(m_maxDepth));
4948 for(uint8_t c = 0; c < m_treeConstants->nNodes; ++c){
4949 bool isNodePeriodic = m_octree.isNodePeriodic(&octant, c);
4950 if (!isNodePeriodic) {
4951 bool isNodeBoundary = octant.getNodeBound(c);
4952 if (isNodeBoundary) {
4953 continue;
4954 }
4955 }
4956
4957 std::array<uint32_t, 3> virtualOctantOrigin = octantCoord;
4958 if (isNodePeriodic) {
4959 std::array<int64_t, 3> periodicOffset = m_octree.getNodePeriodicOffset(octant, c);
4960 virtualOctantOrigin[0] = static_cast<uint32_t>(virtualOctantOrigin[0] + periodicOffset[0]);
4961 virtualOctantOrigin[1] = static_cast<uint32_t>(virtualOctantOrigin[1] + periodicOffset[1]);
4962 virtualOctantOrigin[2] = static_cast<uint32_t>(virtualOctantOrigin[2] + periodicOffset[2]);
4963 }
4964
4965 m_octree.computeVirtualNodeNeighOffsets(octantLevel, c, maxNodeNeighLevel, &virtualNeighOffsets);
4966
4967 for(const std::array<int64_t, 3> &virtualNeighOffset : virtualNeighOffsets){
4968 std::array<uint32_t, 3> virtualNeighCoord = virtualOctantOrigin;
4969 virtualNeighCoord[0] = static_cast<uint32_t>(virtualNeighCoord[0] + virtualNeighOffset[0]);
4970 virtualNeighCoord[1] = static_cast<uint32_t>(virtualNeighCoord[1] + virtualNeighOffset[1]);
4971 virtualNeighCoord[2] = static_cast<uint32_t>(virtualNeighCoord[2] + virtualNeighOffset[2]);
4972 uint64_t virtualNeighMorton = PABLO::computeMorton(m_dim, virtualNeighCoord[0], virtualNeighCoord[1], virtualNeighCoord[2]);
4973
4974 int virtualNeighProc = findOwner(virtualNeighMorton);
4975 if (virtualNeighProc != m_rank) {
4976 neighProcs.insert(virtualNeighProc);
4977 }
4978 }
4979 }
4980
4981 // Build list of internal and process borders octants
4982 if (neighProcs.empty()){
4983 m_internals[countint] = &octant;
4984 countint++;
4985 } else {
4986 m_pborders[countpbd] = &octant;
4987 countpbd++;
4988
4989 for (int neighProc : neighProcs) {
4990 assert(neighProc != m_rank);
4991 m_bordersPerProc[neighProc].push_back(idx);
4992 }
4993 }
4994 }
4995 m_pborders.resize(countpbd);
4996 m_pborders.shrink_to_fit();
4997 m_internals.resize(countint);
4998 m_internals.shrink_to_fit();
4999
5000 // Build ghosts
5001 buildGhostOctants(m_bordersPerProc, std::vector<AccretionData>());
5002 }
5003
5007 void
5008 ParaTree::computeGhostHalo(){
5009 // Build first layer of ghosts
5010 setPboundGhosts();
5011
5012 // Early return if we need to build only one layer
5013 if(m_nofGhostLayers <= 1){
5014 return;
5015 }
5016
5017 //
5018 // Accrete sources
5019 //
5020 // We don't build ghost layers directly, instead we identify the
5021 // internal cells that are ghosts for the neighboring process and
5022 // we use this list to create the ghosts. We use the term "sources"
5023 // to identify internal cells that are ghosts for the neighboring
5024 // process. For each layer of ghosts, a corresponding layer of
5025 // sources exists.
5026 //
5027 // Sources are identified one layer at a time. The first layer is
5028 // already known: the process-border octants. The neighbors of
5029 // process-border octants are the second layer of sources; the
5030 // neighbors of the second layer of sources are the third layer,
5031 // and so on an so forth.
5032 //
5033 // To identify the sources, an auxiliary data structure is used. This
5034 // data structure is called accretion and contains the list of sources
5035 // currently identified (population), a list of octants to be used for
5036 // building the next layer of sources (seeds) and the rank on which
5037 // the sources gathered by the accretion will be ghosts.
5038 //
5039 // The identification of the sources starts creating one accretions for
5040 // each of the neighboring processes. The accretions are initialized
5041 // using the process-border octants already build: those octants are
5042 // the first layer of sources and the seeds for the generation of the
5043 // second layer. Adding the internal neighbors of the internal seeds to
5044 // the population, accretions are grown one layer at a time. When an
5045 // accretion reaches a neighboring processes (i.e., when a first-layer
5046 // ghost enters in the list of foreign seeds), we communicate to the
5047 // owner of the ghost to create a new accretion and continue the search
5048 // for the sources. At the end of the procedure, the population of the
5049 // accretions on each process will contain the desired sources.
5050
5051 // Initialize cache for 1-rings of the internal octants
5052 std::unordered_map<uint32_t, std::vector<uint64_t>> oneRingsCache;
5053 oneRingsCache.reserve(getNumOctants());
5054
5055 // Initialize data communicator
5056 DataCommunicator accretionDataCommunicator(m_comm);
5057
5058 // Initialize accretions
5059 std::vector<AccretionData> accretions;
5060 initializeGhostHaloAccretions(&accretions);
5061
5062 // Grow the accretions
5063 for(std::size_t layer = 1; layer < m_nofGhostLayers; ++layer){
5064 // Exchange accretions
5065 //
5066 // When a ghost is incorporated in the seeds, the accretion
5067 // needs to continue on the process that owns the ghost.
5068 exchangeGhostHaloAccretions(&accretionDataCommunicator, &accretions);
5069
5070 // Grow accretions
5071 growGhostHaloAccretions(&oneRingsCache, &accretions);
5072 }
5073
5074 // To correctly identify the population of the last layer of sources,
5075 // we need to exchange the accretions one more time. This allows to
5076 // communicate the foreign seeds found during the last growth to the
5077 // processes that own them.
5078 exchangeGhostHaloAccretions(&accretionDataCommunicator, &accretions);
5079
5080 //
5081 // Extract list of sources
5082 //
5083 // Sources are internal octants that are ghosts for other processes,
5084 // i.e., internal octants on processes borders (pborder octants).
5085 // The population of the accretions is exaclty made of the sources
5086 // for the target rank of the accretion.
5087 for(const AccretionData &accretion : accretions){
5088 int targetRank = accretion.targetRank;
5089
5090 std::vector<uint32_t> &rankBordersPerProc = m_bordersPerProc[targetRank];
5091 rankBordersPerProc.resize(accretion.population.size());
5092
5093 std::size_t index = 0;
5094 for (const auto &populationEntry : accretion.population) {
5095 uint64_t globalIdx = populationEntry.first;
5096 uint32_t localIdx = getLocalIdx(globalIdx);
5097 rankBordersPerProc[index] = localIdx;
5098 ++index;
5099 }
5100
5101 std::sort(rankBordersPerProc.begin(), rankBordersPerProc.end());
5102 }
5103
5104 //
5105 // Build the ghosts
5106 //
5107 buildGhostOctants(m_bordersPerProc, accretions);
5108 }
5109
5118 void
5119 ParaTree::initializeGhostHaloAccretions(std::vector<AccretionData> *accretions) {
5120
5121 const int FIRST_LAYER = 0;
5122
5123 accretions->reserve(m_bordersPerProc.size());
5124 for(const auto &bordersPerProcEntry : m_bordersPerProc){
5125 accretions->emplace_back();
5126 AccretionData &accretion = accretions->back();
5127
5128 // Rank for which the accretion will gather sources
5129 int targetRank = bordersPerProcEntry.first;
5130 accretion.targetRank = targetRank;
5131
5132 // Initialize the first layer
5133 //
5134 // The population and the seeds of the first layer are the octants
5135 // on processes borders.
5136 const std::vector<uint32_t> &rankBordersPerProc = bordersPerProcEntry.second;
5137 const std::size_t nRankBordersPerProc = rankBordersPerProc.size();
5138
5139 accretion.population.reserve(m_nofGhostLayers * nRankBordersPerProc);
5140 accretion.internalSeeds.reserve(nRankBordersPerProc);
5141 accretion.foreignSeeds.reserve(nRankBordersPerProc);
5142 for(uint32_t pborderLocalIdx : rankBordersPerProc){
5143 uint64_t pborderGlobalIdx = getGlobalIdx(pborderLocalIdx);
5144 if (isInternal(pborderGlobalIdx)) {
5145 accretion.population.insert({pborderGlobalIdx, FIRST_LAYER});
5146 accretion.internalSeeds.insert({pborderGlobalIdx, FIRST_LAYER});
5147 } else {
5148 accretion.foreignSeeds.insert({pborderGlobalIdx, FIRST_LAYER});
5149 }
5150 }
5151 }
5152 }
5153
5163 void
5164 ParaTree::growGhostHaloAccretions(std::unordered_map<uint32_t, std::vector<uint64_t>> *oneRingsCache,
5165 std::vector<AccretionData> *accretions) {
5166
5167 // The neighbour of the internal seeds are the next layer of sources.
5168 u32vector seedNeighLocalIds;
5169 bvector seedNeighGhostFlag;
5170 for(AccretionData &accretion : *accretions){
5171 // If the accretion doesn't have internal seeds we can skip it
5172 std::size_t nSeeds = accretion.internalSeeds.size();
5173 if (nSeeds == 0) {
5174 continue;
5175 }
5176
5177 // Rank for which accretion is gathering data
5178 const int targetRank = accretion.targetRank;
5179
5180 // Seeds
5181 //
5182 // We make a copy of the seeds and then we clear the original
5183 // list in order to generate the seeds for the next layer.
5184 std::unordered_map<uint64_t, int> currentInternalSeeds;
5185 currentInternalSeeds.reserve(accretion.internalSeeds.size());
5186 accretion.internalSeeds.swap(currentInternalSeeds);
5187
5188 // The next layer is obtained adding the 1-ring neighbours of
5189 // the internal octants of the previous layer.
5190 for(const auto &seedEntry : currentInternalSeeds){
5191 // Get seed information
5192 int seedLayer = seedEntry.second;
5193 uint64_t seedGlobalIdx = seedEntry.first;
5194 uint32_t seedLocalIdx = getLocalIdx(seedGlobalIdx);
5195
5196 // Find the 1-ring of the source
5197 auto oneRingsCacheItr = oneRingsCache->find(seedLocalIdx);
5198 if (oneRingsCacheItr == oneRingsCache->end()) {
5199 const Octant *seedOctant = getOctant(seedLocalIdx);
5200 findAllCodimensionNeighbours(seedOctant, seedNeighLocalIds, seedNeighGhostFlag);
5201 size_t nSeedNeighs = seedNeighLocalIds.size();
5202
5203 oneRingsCacheItr = oneRingsCache->insert({seedLocalIdx, std::vector<uint64_t>()}).first;
5204 oneRingsCacheItr->second.resize(nSeedNeighs + 1);
5205 for(size_t n = 0; n < nSeedNeighs; ++n){
5206 if (!seedNeighGhostFlag[n]) {
5207 oneRingsCacheItr->second[n] = getGlobalIdx(seedNeighLocalIds[n]);
5208 } else {
5209 oneRingsCacheItr->second[n] = getGhostGlobalIdx(seedNeighLocalIds[n]);
5210 }
5211 }
5212 oneRingsCacheItr->second[nSeedNeighs] = getGlobalIdx(seedLocalIdx);
5213 }
5214 const std::vector<uint64_t> &seedOneRing = oneRingsCacheItr->second;
5215
5216 // Add the 1-ring of the octant to the sources
5217 for(uint64_t neighGlobalIdx : seedOneRing){
5218 // Discard octants already in the population
5219 if (accretion.population.count(neighGlobalIdx) > 0) {
5220 continue;
5221 }
5222
5223 // Get neighbour information
5224 bool isNeighInternal = isInternal(neighGlobalIdx);
5225
5226 int neighRank;
5227 if (isNeighInternal) {
5228 neighRank = getRank();
5229 } else {
5230 neighRank = getOwnerRank(neighGlobalIdx);
5231 }
5232
5233 // Add the neighbour to the population
5234 //
5235 // Population should only contain internal octants.
5236 if (isNeighInternal) {
5237 accretion.population.insert({neighGlobalIdx, seedLayer + 1});
5238 }
5239
5240 // Add the neighbour to the seeds
5241 //
5242 // Internal seeds contains only internal octants, foreign
5243 // seeds containt ghost octants that are not owned by the
5244 // target rank (if an octant is owned by the target rank,
5245 // by definition it will not be a source for that rank).
5246 if (isNeighInternal) {
5247 accretion.internalSeeds.insert({neighGlobalIdx, seedLayer + 1});
5248 } else if (neighRank != targetRank) {
5249 accretion.foreignSeeds.insert({neighGlobalIdx, seedLayer + 1});
5250 }
5251 }
5252 }
5253 }
5254 }
5255
5266 void
5267 ParaTree::exchangeGhostHaloAccretions(DataCommunicator *dataCommunicator,
5268 std::vector<AccretionData> *accretions) {
5269
5270 // Generate accretions that has to be sent to other processes
5271 //
5272 // When the accretion reaches a foreign process (i.e., a ghost is
5273 // added to the seeds), the ranks that owns the ghost seed have to
5274 // continue the propagation of those seeds (which will be local
5275 // seeds for that foreign rank) locally. Therefore we need to
5276 // communicate to the ranks that own ghosts seeds the information
5277 // about the accretion and the list of seeds.
5278 std::unordered_map<int, std::vector<AccretionData>> foreignAccretions;
5279 for(const AccretionData &accretion : *accretions){
5280 for(const auto &seedEntry : accretion.foreignSeeds){
5281 // Find the foreign accretion the ghost seed belogns to
5282 int seedRank = getOwnerRank(seedEntry.first);
5283 std::vector<AccretionData> &foreignRankAccretions = foreignAccretions[seedRank];
5284
5285 auto foreignRankAccretionsItr = foreignRankAccretions.begin();
5286 for (; foreignRankAccretionsItr != foreignRankAccretions.end(); ++foreignRankAccretionsItr) {
5287 if (foreignRankAccretionsItr->targetRank!= accretion.targetRank) {
5288 continue;
5289 }
5290
5291 break;
5292 }
5293
5294 if (foreignRankAccretionsItr == foreignRankAccretions.end()) {
5295 foreignRankAccretions.emplace_back();
5296 foreignRankAccretionsItr = foreignRankAccretions.begin() + foreignRankAccretions.size() - 1;
5297 foreignRankAccretionsItr->targetRank = accretion.targetRank;
5298 }
5299
5300 AccretionData &foreignAccretion = *foreignRankAccretionsItr;
5301
5302 // Add the seed
5303 foreignAccretion.internalSeeds.insert(seedEntry);
5304 }
5305 }
5306
5307 // Early return if no communications are needed
5308 bool foreignAccrationExchangeNeeded = (foreignAccretions.size() > 0);
5309 MPI_Allreduce(MPI_IN_PLACE, &foreignAccrationExchangeNeeded, 1, MPI_C_BOOL, MPI_LOR, m_comm);
5310 if (!foreignAccrationExchangeNeeded) {
5311 return;
5312 }
5313
5314 // Clear previous communications
5315 dataCommunicator->clearAllSends();
5316 dataCommunicator->clearAllRecvs();
5317
5318 // Fill send buffers with accretions data
5319 for(const auto &foreignAccretionEntry : foreignAccretions){
5320 int receiverRank = foreignAccretionEntry.first;
5321
5322 // Evaluate buffer size
5323 std::size_t buffSize = 0;
5324 buffSize += sizeof(std::size_t);
5325 for(const auto &foreignAccretion: foreignAccretionEntry.second){
5326 std::size_t nSeeds = foreignAccretion.internalSeeds.size();
5327
5328 buffSize += sizeof(int);
5329 buffSize += sizeof(std::size_t);
5330 buffSize += nSeeds * (sizeof(uint64_t) + sizeof(int));
5331 }
5332
5333 dataCommunicator->setSend(receiverRank, buffSize);
5334
5335 // Fill buffer
5336 SendBuffer &sendBuffer = dataCommunicator->getSendBuffer(receiverRank);
5337
5338 const std::vector<AccretionData> &foreignRankAccretions = foreignAccretionEntry.second;
5339
5340 sendBuffer << foreignRankAccretions.size();
5341 for(const auto &foreignAccretion: foreignRankAccretions){
5342 sendBuffer << foreignAccretion.targetRank;
5343 sendBuffer << foreignAccretion.internalSeeds.size();
5344 for(const auto &seedEntry : foreignAccretion.internalSeeds){
5345 sendBuffer << seedEntry.first;
5346 sendBuffer << seedEntry.second;
5347 }
5348 }
5349 }
5350
5351 // Start communications
5352 dataCommunicator->discoverRecvs();
5353 dataCommunicator->startAllRecvs();
5354 dataCommunicator->startAllSends();
5355
5356 // Receive the accretiation to grow on behalf of neighbour processes
5357 int nCompletedRecvs = 0;
5358 while (nCompletedRecvs < dataCommunicator->getRecvCount()) {
5359 int senderRank = dataCommunicator->waitAnyRecv();
5360 RecvBuffer &recvBuffer = dataCommunicator->getRecvBuffer(senderRank);
5361
5362 std::size_t nForeignAccretions;
5363 recvBuffer >> nForeignAccretions;
5364
5365 for (std::size_t k = 0; k < nForeignAccretions; ++k) {
5366 // Target rank
5367 int targetRank;
5368 recvBuffer >> targetRank;
5369
5370 // Get the accretion to update
5371 //
5372 // If an accretions with the requeste owner and target
5373 // ranks doesn't exist create a new one.
5374 auto accretionsItr = accretions->begin();
5375 for(; accretionsItr != accretions->end(); ++accretionsItr){
5376 if (accretionsItr->targetRank!= targetRank) {
5377 continue;
5378 }
5379
5380 break;
5381 }
5382
5383 if (accretionsItr == accretions->end()) {
5384 accretions->emplace_back();
5385 accretionsItr = accretions->begin() + accretions->size() - 1;
5386 accretionsItr->targetRank = targetRank;
5387 }
5388
5389 AccretionData &accretion = *accretionsItr;
5390
5391 // Initialize accretion seeds and population
5392 //
5393 // We are receiving only internal octants, there is no need
5394 // to explicitly check if and octant is internal before add
5395 // it to the population and to the internal seeds.
5396 std::size_t nSeeds;
5397 recvBuffer >> nSeeds;
5398
5399 for(std::size_t n = 0; n < nSeeds; ++n){
5400 uint64_t globalIdx;
5401 recvBuffer >> globalIdx;
5402
5403 int layer;
5404 recvBuffer >> layer;
5405
5406 assert(isInternal(globalIdx));
5407 accretion.population.insert({globalIdx, layer});
5408 accretion.internalSeeds.insert({globalIdx, layer});
5409 }
5410 }
5411
5412 ++nCompletedRecvs;
5413 }
5414
5415 // Wait until all exchanges are completed
5416 dataCommunicator->waitAllSends();
5417 }
5418
5424 void
5425 ParaTree::buildGhostOctants(const std::map<int, u32vector> &bordersPerProc,
5426 const std::vector<AccretionData> &accretions){
5427
5428 DataCommunicator ghostDataCommunicator(m_comm);
5429
5430 // Binary size of a ghost entry in the communication buffer
5431 const std::size_t GHOST_ENTRY_BINARY_SIZE = sizeof(uint64_t) + Octant::getBinarySize() + sizeof(int);
5432
5433 // Fill the send buffers with source octants
5434 //
5435 // A source octant is an internal octants that is ghosts on other
5436 // process.
5437 for(const auto &bordersPerProcEntry : bordersPerProc){
5438 int rank = bordersPerProcEntry.first;
5439 const std::vector<uint32_t> &rankBordersPerProc = bordersPerProcEntry.second;
5440 std::size_t nRankBordersPerProc = rankBordersPerProc.size();
5441
5442 // Get the accretion associated with this rank
5443 //
5444 // If there are no accretions it means we are building the first
5445 // layer of ghosts.
5446 std::vector<AccretionData>::const_iterator accretionsItr;
5447 if (accretions.size() > 0) {
5448 bool accretionFound = false;
5449 BITPIT_UNUSED(accretionFound);
5450 accretionsItr = accretions.begin();
5451 for (; accretionsItr != accretions.end(); ++accretionsItr) {
5452 if (accretionsItr->targetRank == rank) {
5453 accretionFound = true;
5454 break;
5455 }
5456 }
5457 assert(accretionFound);
5458 } else {
5459 accretionsItr = accretions.end();
5460 }
5461
5462 // Initialize the send
5463 std::size_t buffSize = GHOST_ENTRY_BINARY_SIZE * nRankBordersPerProc;
5464 ghostDataCommunicator.setSend(rank, buffSize);
5465
5466 // Fill the buffer
5467 SendBuffer &sendBuffer = ghostDataCommunicator.getSendBuffer(rank);
5468 for(uint32_t sourceLocalIdx : rankBordersPerProc){
5469 // Global index
5470 uint64_t sourceGlobalIdx = getGlobalIdx(sourceLocalIdx);
5471 sendBuffer << sourceGlobalIdx;
5472
5473 // Source data
5474 sendBuffer << m_octree.m_octants[sourceLocalIdx];
5475
5476 // Layer information
5477 //
5478 // If there are no accretions it means we are building the
5479 // first layer of ghosts.
5480 int layer;
5481 if (accretionsItr != accretions.end()) {
5482 layer = accretionsItr->population.at(sourceGlobalIdx);
5483 } else {
5484 layer = 0;
5485 }
5486
5487 sendBuffer << layer;
5488 }
5489 }
5490
5491 // Discover the receives
5492 ghostDataCommunicator.discoverRecvs();
5493 ghostDataCommunicator.startAllRecvs();
5494 ghostDataCommunicator.startAllSends();
5495
5496 // Get the ranks from which ghosts will be received
5497 std::vector<int> ghostCommunicatorRecvsRanks = ghostDataCommunicator.getRecvRanks();
5498 std::sort(ghostCommunicatorRecvsRanks.begin(), ghostCommunicatorRecvsRanks.end());
5499
5500 // Prepare ghost data structures
5501 uint32_t nGhosts = 0;
5502 for (int rank : ghostCommunicatorRecvsRanks) {
5503 RecvBuffer &recvBuffer = ghostDataCommunicator.getRecvBuffer(rank);
5504 std::size_t nRankGhosts = recvBuffer.getSize() / GHOST_ENTRY_BINARY_SIZE;
5505 nGhosts += nRankGhosts;
5506 }
5507
5508 m_octree.m_ghosts.resize(nGhosts);
5509 m_octree.m_globalIdxGhosts.resize(nGhosts);
5510
5511 // Receive the ghosts
5512 //
5513 // Ghosts have to be received following the rank order.
5514 uint32_t ghostLocalIdx = 0;
5515 for (int rank : ghostCommunicatorRecvsRanks) {
5516 ghostDataCommunicator.waitRecv(rank);
5517 RecvBuffer &recvBuffer = ghostDataCommunicator.getRecvBuffer(rank);
5518
5519 std::size_t nRankGhosts = recvBuffer.getSize() / GHOST_ENTRY_BINARY_SIZE;
5520 for(std::size_t n = 0; n < nRankGhosts; ++n){
5521 // Assign the global index
5522 uint64_t ghostGlobalIdx;
5523 recvBuffer >> ghostGlobalIdx;
5524 m_octree.m_globalIdxGhosts[ghostLocalIdx] = ghostGlobalIdx;
5525
5526 // Build the ghosts
5527 //
5528 // The layer of the received octant will be overwritten with
5529 // the actual ghost layer.
5530 Octant &ghostOctant = m_octree.m_ghosts[ghostLocalIdx];
5531 recvBuffer >> ghostOctant;
5532
5533 // Set the layer of the ghost
5534 int ghostLayer;
5535 recvBuffer >> ghostLayer;
5536 ghostOctant.setGhostLayer(ghostLayer);
5537
5538 // Increase the ghost index
5539 ++ghostLocalIdx;
5540 }
5541 }
5542
5543 // Wait for the communications to complete
5544 ghostDataCommunicator.waitAllSends();
5545 }
5546
5551 bool
5552 ParaTree::commMarker() {
5553 // If the tree is not partitioned, there is nothing to communicate.
5554 if (m_serial) {
5555 return false;
5556 }
5557
5558 // Binary size of a marker entry in the communication buffer
5559 const std::size_t MARKER_ENTRY_BINARY_SIZE = sizeof(Octant::m_marker);
5560
5561 // Fill communication buffer with level and marker
5562 //
5563 // It visits every element in m_bordersPerProc (one for every neighbor proc)
5564 // for every element it visits the border octants it contains and write them in the bitpit communication structure, DataCommunicator
5565 // this structure has a buffer for every proc containing the octants to be sent to that proc written in a char* buffer
5566 DataCommunicator markerCommunicator(m_comm);
5567
5568 for(const auto &bordersPerProcEntry : m_bordersPerProc){
5569 int rank = bordersPerProcEntry.first;
5570 const std::vector<uint32_t> &rankBordersPerProc = m_bordersPerProc.at(rank);
5571 const std::size_t nRankBorders = rankBordersPerProc.size();
5572
5573 std::size_t buffSize = nRankBorders * MARKER_ENTRY_BINARY_SIZE;
5574 markerCommunicator.setSend(rank, buffSize);
5575
5576 SendBuffer &sendBuffer = markerCommunicator.getSendBuffer(rank);
5577 for(std::size_t i = 0; i < nRankBorders; ++i){
5578 const Octant &octant = m_octree.m_octants[rankBordersPerProc[i]];
5579 sendBuffer << octant.getMarker();
5580 }
5581 }
5582
5583 markerCommunicator.discoverRecvs();
5584 markerCommunicator.startAllRecvs();
5585 markerCommunicator.startAllSends();
5586
5587 // Read level and marker from communication buffer
5588 //
5589 // every receive buffer is visited, and read octant by octant.
5590 // every ghost octant level and marker are updated
5591 std::vector<int> recvRanks = markerCommunicator.getRecvRanks();
5592 std::sort(recvRanks.begin(), recvRanks.end());
5593
5594 bool updated = false;
5595 uint32_t ghostIdx = 0;
5596 for(int rank : recvRanks){
5597 markerCommunicator.waitRecv(rank);
5598 RecvBuffer &recvBuffer = markerCommunicator.getRecvBuffer(rank);
5599
5600 const std::size_t nRankGhosts = recvBuffer.getSize() / MARKER_ENTRY_BINARY_SIZE;
5601 for(std::size_t i = 0; i < nRankGhosts; ++i){
5602 int8_t marker;
5603 recvBuffer >> marker;
5604
5605 Octant &octant = m_octree.m_ghosts[ghostIdx];
5606 if (octant.getMarker() != marker) {
5607 octant.setMarker(marker);
5608 updated = true;
5609 }
5610
5611 ++ghostIdx;
5612 }
5613 }
5614
5615 markerCommunicator.waitAllSends();
5616
5617 return updated;
5618 }
5619#endif
5620
5624 void
5625 ParaTree::updateAfterCoarse(){
5626
5627 updateAdapt();
5628
5629#if BITPIT_ENABLE_MPI==1
5630 if(!m_serial){
5631 updateGlobalFirstDescMorton();
5632 updateGlobalLasttDescMorton();
5633 }
5634#endif
5635 }
5636
5643 void
5644 ParaTree::balance21(bool verbose, bool balanceNewOctants){
5645
5646 // Print header
5647 if (verbose){
5648 (*m_log) << "---------------------------------------------" << endl;
5649 (*m_log) << " 2:1 BALANCE (balancing Marker before Adapt)" << endl;
5650 (*m_log) << " " << endl;
5651 (*m_log) << " Iterative procedure " << endl;
5652 (*m_log) << " " << endl;
5653 }
5654
5655 bool initialStep = true;
5656 while (true) {
5657#if BITPIT_ENABLE_MPI==1
5658 // Exchange markers across processes
5659 if (!m_serial) {
5660 bool markersUpdated = commMarker();
5661 if (!initialStep) {
5662 MPI_Allreduce(MPI_IN_PLACE, &markersUpdated, 1, MPI_C_BOOL, MPI_LOR, m_comm);
5663 if (!markersUpdated) {
5664 break;
5665 }
5666 }
5667 }
5668#endif
5669
5670 // Balance 2:1 the tree
5671 //
5672 // After the first step only ghost octants are checked for balance, that's because
5673 // starting from the second step we only need to propagate marker information across
5674 // the processes.
5675 bool checkInterior = initialStep;
5676 bool checkGhosts = true;
5677
5678 bool balanceUpdated = m_octree.localBalance(balanceNewOctants, checkInterior, checkGhosts);
5679#if BITPIT_ENABLE_MPI==1
5680 if (!m_serial) {
5681 MPI_Allreduce(MPI_IN_PLACE, &balanceUpdated, 1, MPI_C_BOOL, MPI_LOR, m_comm);
5682 if (!balanceUpdated) {
5683 break;
5684 }
5685 } else {
5686 break;
5687 }
5688#else
5689 BITPIT_UNUSED(balanceUpdated);
5690 break;
5691#endif
5692
5693 // This is not anymore the first balance step
5694 initialStep = false;
5695
5696 }
5697
5698 // Print footer
5699 if (verbose){
5700 (*m_log) << " 2:1 Balancing reached " << endl;
5701 (*m_log) << " " << endl;
5702 (*m_log) << "---------------------------------------------" << endl;
5703 }
5704 };
5705
5706 // =================================================================================== //
5707 // TESTING OUTPUT METHODS //
5708 // =================================================================================== //
5709
5715 void
5716 ParaTree::write(const std::string &filename) {
5717
5718 // Update connectivity
5719 m_octree.updateConnectivity();
5720
5721 // Write output file
5722 stringstream name;
5723 name << "s" << std::setfill('0') << std::setw(4) << m_nproc << "-p" << std::setfill('0') << std::setw(4) << m_rank << "-" << filename << ".vtu";
5724
5725 ofstream out(name.str().c_str());
5726 if(!out.is_open()){
5727 stringstream ss;
5728 ss << filename << "*.vtu cannot be opened and it won't be written." << endl;
5729 (*m_log) << ss.str();
5730 return;
5731 }
5732 int nofNodes = m_octree.m_nodes.size();
5733 int nofOctants = m_octree.m_connectivity.size();
5734 int nofGhosts = m_octree.m_ghostsConnectivity.size();
5735 int nofAll = nofGhosts + nofOctants;
5736 out << "<?xml version=\"1.0\"?>" << endl
5737 << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
5738 << " <UnstructuredGrid>" << endl
5739 << " <Piece NumberOfCells=\"" << m_octree.m_connectivity.size() + m_octree.m_ghostsConnectivity.size() << "\" NumberOfPoints=\"" << m_octree.m_nodes.size() << "\">" << endl;
5740 out << " <Points>" << endl
5741 << " <DataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\""<< 3 <<"\" format=\"ascii\">" << endl
5742 << " " << std::fixed;
5743 for(int i = 0; i < nofNodes; i++)
5744 {
5745 for(int j = 0; j < 3; ++j){
5746 if (j==0) out << std::setprecision(6) << m_trans.mapX(m_octree.m_nodes[i][j]) << " ";
5747 if (j==1) out << std::setprecision(6) << m_trans.mapY(m_octree.m_nodes[i][j]) << " ";
5748 if (j==2) out << std::setprecision(6) << m_trans.mapZ(m_octree.m_nodes[i][j]) << " ";
5749 }
5750 if((i+1)%4==0 && i!=nofNodes-1)
5751 out << endl << " ";
5752 }
5753 out << endl << " </DataArray>" << endl
5754 << " </Points>" << endl
5755 << " <Cells>" << endl
5756 << " <DataArray type=\"UInt64\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5757 << " ";
5758 for(int i = 0; i < nofOctants; i++)
5759 {
5760 for(int j = 0; j < m_treeConstants->nNodes; j++)
5761 {
5762 int jj = j;
5763 if (m_dim==2){
5764 if (j<2){
5765 jj = j;
5766 }
5767 else if(j==2){
5768 jj = 3;
5769 }
5770 else if(j==3){
5771 jj = 2;
5772 }
5773 }
5774 out << m_octree.m_connectivity[i][jj] << " ";
5775 }
5776 if((i+1)%3==0 && i!=nofOctants-1)
5777 out << endl << " ";
5778 }
5779 for(int i = 0; i < nofGhosts; i++)
5780 {
5781 for(int j = 0; j < m_treeConstants->nNodes; j++)
5782 {
5783 int jj = j;
5784 if (m_dim==2){
5785 if (j<2){
5786 jj = j;
5787 }
5788 else if(j==2){
5789 jj = 3;
5790 }
5791 else if(j==3){
5792 jj = 2;
5793 }
5794 }
5795 out << m_octree.m_ghostsConnectivity[i][jj] << " ";
5796 }
5797 if((i+1)%3==0 && i!=nofGhosts-1)
5798 out << endl << " ";
5799 }
5800 out << endl << " </DataArray>" << endl
5801 << " <DataArray type=\"UInt64\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5802 << " ";
5803 for(int i = 0; i < nofAll; i++)
5804 {
5805 out << (i+1)*m_treeConstants->nNodes << " ";
5806 if((i+1)%12==0 && i!=nofAll-1)
5807 out << endl << " ";
5808 }
5809 out << endl << " </DataArray>" << endl
5810 << " <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5811 << " ";
5812 for(int i = 0; i < nofAll; i++)
5813 {
5814 int type;
5815 type = 5 + (m_dim*2);
5816 out << type << " ";
5817 if((i+1)%12==0 && i!=nofAll-1)
5818 out << endl << " ";
5819 }
5820 out << endl << " </DataArray>" << endl
5821 << " </Cells>" << endl
5822 << " </Piece>" << endl
5823 << " </UnstructuredGrid>" << endl
5824 << "</VTKFile>" << endl;
5825
5826
5827 if(m_rank == 0){
5828 name.str("");
5829 name << "s" << std::setfill('0') << std::setw(4) << m_nproc << "-" << filename << ".pvtu";
5830 ofstream pout(name.str().c_str());
5831 if(!pout.is_open()){
5832 stringstream ss;
5833 ss << filename << "*.pvtu cannot be opened and it won't be written." << endl;
5834 (*m_log) << ss.str();
5835 return;
5836 }
5837
5838 pout << "<?xml version=\"1.0\"?>" << endl
5839 << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
5840 << " <PUnstructuredGrid GhostLevel=\"0\">" << endl
5841 << " <PPointData>" << endl
5842 << " </PPointData>" << endl
5843 << " <PCellData Scalars=\"\">" << endl;
5844 pout << " </PCellData>" << endl
5845 << " <PPoints>" << endl
5846 << " <PDataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\"3\"/>" << endl
5847 << " </PPoints>" << endl;
5848 for(int i = 0; i < m_nproc; i++)
5849 pout << " <Piece Source=\"s" << std::setw(4) << std::setfill('0') << m_nproc << "-p" << std::setw(4) << std::setfill('0') << i << "-" << filename << ".vtu\"/>" << endl;
5850 pout << " </PUnstructuredGrid>" << endl
5851 << "</VTKFile>";
5852
5853 pout.close();
5854
5855 }
5856#if BITPIT_ENABLE_MPI==1
5857 if (isCommSet()) {
5858 MPI_Barrier(m_comm);
5859 }
5860#endif
5861
5862 }
5863
5870 void
5871 ParaTree::writeTest(const std::string &filename, vector<double> data) {
5872
5873 // Update connectivity
5874 m_octree.updateConnectivity();
5875
5876 // Write output file
5877 stringstream name;
5878 name << "s" << std::setfill('0') << std::setw(4) << m_nproc << "-p" << std::setfill('0') << std::setw(4) << m_rank << "-" << filename << ".vtu";
5879
5880 ofstream out(name.str().c_str());
5881 if(!out.is_open()){
5882 stringstream ss;
5883 ss << filename << "*.vtu cannot be opened and it won't be written.";
5884 (*m_log) << ss.str();
5885 return;
5886 }
5887 int nofNodes = m_octree.m_nodes.size();
5888 int nofOctants = m_octree.m_connectivity.size();
5889 int nofAll = nofOctants;
5890 out << "<?xml version=\"1.0\"?>" << endl
5891 << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
5892 << " <UnstructuredGrid>" << endl
5893 << " <Piece NumberOfCells=\"" << m_octree.m_connectivity.size() << "\" NumberOfPoints=\"" << m_octree.m_nodes.size() << "\">" << endl;
5894 out << " <CellData Scalars=\"Data\">" << endl;
5895 out << " <DataArray type=\"Float64\" Name=\"Data\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5896 << " " << std::fixed;
5897 int ndata = m_octree.m_connectivity.size();
5898 for(int i = 0; i < ndata; i++)
5899 {
5900 out << std::setprecision(6) << data[i] << " ";
5901 if((i+1)%4==0 && i!=ndata-1)
5902 out << endl << " ";
5903 }
5904 out << endl << " </DataArray>" << endl
5905 << " </CellData>" << endl
5906 << " <Points>" << endl
5907 << " <DataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\""<< 3 <<"\" format=\"ascii\">" << endl
5908 << " " << std::fixed;
5909 for(int i = 0; i < nofNodes; i++)
5910 {
5911 for(int j = 0; j < 3; ++j){
5912 if (j==0) out << std::setprecision(6) << m_trans.mapX(m_octree.m_nodes[i][j]) << " ";
5913 if (j==1) out << std::setprecision(6) << m_trans.mapY(m_octree.m_nodes[i][j]) << " ";
5914 if (j==2) out << std::setprecision(6) << m_trans.mapZ(m_octree.m_nodes[i][j]) << " ";
5915 }
5916 if((i+1)%4==0 && i!=nofNodes-1)
5917 out << endl << " ";
5918 }
5919 out << endl << " </DataArray>" << endl
5920 << " </Points>" << endl
5921 << " <Cells>" << endl
5922 << " <DataArray type=\"UInt64\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5923 << " ";
5924 for(int i = 0; i < nofOctants; i++)
5925 {
5926 for(int j = 0; j < m_treeConstants->nNodes; j++)
5927 {
5928 int jj = j;
5929 if (m_dim==2){
5930 if (j<2){
5931 jj = j;
5932 }
5933 else if(j==2){
5934 jj = 3;
5935 }
5936 else if(j==3){
5937 jj = 2;
5938 }
5939 }
5940 out << m_octree.m_connectivity[i][jj] << " ";
5941 }
5942 if((i+1)%3==0 && i!=nofOctants-1)
5943 out << endl << " ";
5944 }
5945 out << endl << " </DataArray>" << endl
5946 << " <DataArray type=\"UInt64\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5947 << " ";
5948 for(int i = 0; i < nofAll; i++)
5949 {
5950 out << (i+1)*m_treeConstants->nNodes << " ";
5951 if((i+1)%12==0 && i!=nofAll-1)
5952 out << endl << " ";
5953 }
5954 out << endl << " </DataArray>" << endl
5955 << " <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
5956 << " ";
5957 for(int i = 0; i < nofAll; i++)
5958 {
5959 int type;
5960 type = 5 + (m_dim*2);
5961 out << type << " ";
5962 if((i+1)%12==0 && i!=nofAll-1)
5963 out << endl << " ";
5964 }
5965 out << endl << " </DataArray>" << endl
5966 << " </Cells>" << endl
5967 << " </Piece>" << endl
5968 << " </UnstructuredGrid>" << endl
5969 << "</VTKFile>" << endl;
5970
5971
5972 if(m_rank == 0){
5973 name.str("");
5974 name << "s" << std::setfill('0') << std::setw(4) << m_nproc << "-" << filename << ".pvtu";
5975 ofstream pout(name.str().c_str());
5976 if(!pout.is_open()){
5977 stringstream ss;
5978 ss << filename << "*.pvtu cannot be opened and it won't be written." << endl;
5979 (*m_log) << ss.str();
5980 return;
5981 }
5982
5983 pout << "<?xml version=\"1.0\"?>" << endl
5984 << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
5985 << " <PUnstructuredGrid GhostLevel=\"0\">" << endl
5986 << " <PPointData>" << endl
5987 << " </PPointData>" << endl
5988 << " <PCellData Scalars=\"Data\">" << endl
5989 << " <PDataArray type=\"Float64\" Name=\"Data\" NumberOfComponents=\"1\"/>" << endl
5990 << " </PCellData>" << endl
5991 << " <PPoints>" << endl
5992 << " <PDataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\"3\"/>" << endl
5993 << " </PPoints>" << endl;
5994 for(int i = 0; i < m_nproc; i++)
5995 pout << " <Piece Source=\"s" << std::setw(4) << std::setfill('0') << m_nproc << "-p" << std::setw(4) << std::setfill('0') << i << "-" << filename << ".vtu\"/>" << endl;
5996 pout << " </PUnstructuredGrid>" << endl
5997 << "</VTKFile>";
5998
5999 pout.close();
6000
6001 }
6002#if BITPIT_ENABLE_MPI==1
6003 if (isCommSet()) {
6004 MPI_Barrier(m_comm);
6005 }
6006#endif
6007
6008 }
6009
6010 // =============================================================================== //
6011
6012}
6013
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::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::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::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::getNodePeriodicOffset
std::array< int64_t, 3 > getNodePeriodicOffset(const Octant &octant, uint8_t inode) const
Definition LocalTree.cpp:2392
bitpit::LoggerManager::create
void create(const std::string &name, bool reset=false, int nProcesses=1, int rank=0)
Definition logger.cpp:1327
bitpit::Logger
Message logger.
Definition logger.hpp:143
bitpit::Logger::setDefaultVisibility
void setDefaultVisibility(log::Visibility visibility)
Definition logger.cpp:579
bitpit::Logger::getDefaultVisibility
log::Visibility getDefaultVisibility()
Definition logger.cpp:591
bitpit::Logger::getName
std::string getName() const
Definition logger.cpp:842
bitpit::Octant
Octant class definition.
Definition Octant.hpp:89
bitpit::Octant::computeNodePersistentKey
uint64_t computeNodePersistentKey(uint8_t inode) const
Definition Octant.cpp:635
bitpit::Octant::getGhostLayer
int getGhostLayer() const
Definition Octant.cpp:373
bitpit::Octant::getPbound
bool getPbound(uint8_t face) const
Definition Octant.cpp:333
bitpit::Octant::countChildren
uint8_t countChildren() const
Definition Octant.cpp:700
bitpit::Octant::getLogicalNode
void getLogicalNode(u32array3 &node, uint8_t inode) const
Definition Octant.cpp:541
bitpit::Octant::getLogicalNodes
void getLogicalNodes(u32arr3vector &nodes) const
Definition Octant.cpp:510
bitpit::Octant::getLogicalCoordinates
u32array3 getLogicalCoordinates() const
Definition Octant.cpp:243
bitpit::Octant::getBinarySize
static unsigned int getBinarySize()
Definition Octant.cpp:658
bitpit::Octant::getDim
uint32_t getDim() const
Definition Octant.cpp:237
bitpit::Octant::getNormal
void getNormal(uint8_t iface, i8array3 &normal, const int8_t(&normals)[6][3]) const
Definition Octant.cpp:569
bitpit::Octant::getIsNewR
bool getIsNewR() const
Definition Octant.cpp:355
bitpit::Octant::buildFather
Octant buildFather() const
Definition Octant.cpp:689
bitpit::Octant::getMorton
uint64_t getMorton() const
Definition Octant.cpp:580
bitpit::Octant::setMarker
void setMarker(int8_t marker)
Definition Octant.cpp:385
bitpit::Octant::getBound
bool getBound(uint8_t face) const
Definition Octant.cpp:272
bitpit::Octant::getLogicalSize
uint32_t getLogicalSize() const
Definition Octant.cpp:432
bitpit::Octant::setBalance
void setBalance(bool balance)
Definition Octant.cpp:393
bitpit::Octant::getLogicalArea
uint64_t getLogicalArea() const
Definition Octant.cpp:440
bitpit::Octant::getLogicalCenter
darray3 getLogicalCenter() const
Definition Octant.cpp:458
bitpit::Octant::buildChildren
std::vector< Octant > buildChildren() const
Definition Octant.cpp:713
bitpit::Octant::getLevel
uint8_t getLevel() const
Definition Octant.cpp:259
bitpit::Octant::setGhostLayer
void setGhostLayer(int ghostLayer)
Definition Octant.cpp:419
bitpit::Octant::getIsNewC
bool getIsNewC() const
Definition Octant.cpp:361
bitpit::Octant::getLogicalFaceCenter
darray3 getLogicalFaceCenter(uint8_t iface) const
Definition Octant.cpp:474
bitpit::Octant::getMarker
int8_t getMarker() const
Definition Octant.cpp:265
bitpit::Octant::getFamilySplittingNode
uint8_t getFamilySplittingNode() const
Definition Octant.cpp:863
bitpit::Octant::getIsGhost
bool getIsGhost() const
Definition Octant.cpp:367
bitpit::Octant::getLogicalVolume
uint64_t getLogicalVolume() const
Definition Octant.cpp:448
bitpit::Octant::computeLastDescMorton
uint64_t computeLastDescMorton() const
Definition Octant.cpp:587
bitpit::Octant::getBalance
bool getBalance() const
Definition Octant.cpp:379
bitpit::ParaTree
Para Tree is the user interface class.
Definition ParaTree.hpp:113
bitpit::ParaTree::getNnodesperface
uint8_t getNnodesperface() const
Definition ParaTree.cpp:1071
bitpit::ParaTree::getPointOwner
Octant * getPointOwner(const dvector &point)
Definition ParaTree.cpp:2941
bitpit::ParaTree::write
void write(const std::string &filename)
Definition ParaTree.cpp:5716
bitpit::ParaTree::getLocalIdx
uint32_t getLocalIdx(uint64_t gidx) const
Definition ParaTree.cpp:1616
bitpit::ParaTree::levelToSize
double levelToSize(uint8_t level) const
Definition ParaTree.cpp:4344
bitpit::ParaTree::getTol
double getTol() const
Definition ParaTree.cpp:1244
bitpit::ParaTree::getGhostLocalIdx
uint32_t getGhostLocalIdx(uint64_t gidx) const
Definition ParaTree.cpp:1630
bitpit::ParaTree::getDim
uint8_t getDim() const
Definition ParaTree.cpp:780
bitpit::ParaTree::setMarker
void setMarker(uint32_t idx, int8_t marker)
Definition ParaTree.cpp:1696
bitpit::ParaTree::isFaceOnOctant
bool isFaceOnOctant(const Octant *faceOctant, uint8_t faceIndex, const Octant *octant) const
Definition ParaTree.cpp:3444
bitpit::ParaTree::getNodeLogicalCoordinates
const u32array3 & getNodeLogicalCoordinates(uint32_t node) const
Definition ParaTree.cpp:3877
bitpit::ParaTree::getNumIntersections
uint32_t getNumIntersections() const
Definition ParaTree.cpp:2253
bitpit::ParaTree::getMaxDepth
int8_t getMaxDepth() const
Definition ParaTree.cpp:3743
bitpit::ParaTree::getNumOctants
uint32_t getNumOctants() const
Definition ParaTree.cpp:2118
bitpit::ParaTree::getNchildren
uint8_t getNchildren() const
Definition ParaTree.cpp:1063
bitpit::ParaTree::getPreMapping
void getPreMapping(u32vector &idx, std::vector< int8_t > &markers, std::vector< bool > &isghost)
Definition ParaTree.cpp:3301
bitpit::ParaTree::getFacenode
const uint8_t(* getFacenode() const)[4]
Definition ParaTree.cpp:1186
bitpit::ParaTree::getPbound
bool getPbound(uint32_t idx, uint8_t face) const
Definition ParaTree.cpp:1539
bitpit::ParaTree::expectedOctantAdapt
void expectedOctantAdapt(const Octant *octant, int8_t marker, octvector *result) const
Definition ParaTree.cpp:3187
bitpit::ParaTree::getParallel
bool getParallel() const
Definition ParaTree.cpp:804
bitpit::ParaTree::getNodes
const u32arr3vector & getNodes() const
Definition ParaTree.cpp:3868
bitpit::ParaTree::adaptGlobalCoarse
bool adaptGlobalCoarse(bool mapper_flag=false)
Definition ParaTree.cpp:3654
bitpit::ParaTree::getFamilySplittingNode
uint8_t getFamilySplittingNode(const Octant *) const
Definition ParaTree.cpp:3177
bitpit::ParaTree::adapt
bool adapt(bool mapper_flag=false)
Definition ParaTree.cpp:3558
bitpit::ParaTree::getNormals
const int8_t(* getNormals() const)[3]
Definition ParaTree.cpp:1167
bitpit::ParaTree::getLog
Logger & getLog()
Definition ParaTree.cpp:828
bitpit::ParaTree::findAllNodeNeighbours
void findAllNodeNeighbours(uint32_t idx, uint32_t node, u32vector &neighbours, bvector &isghost)
Definition ParaTree.cpp:2857
bitpit::ParaTree::findGhostAllCodimensionNeighbours
void findGhostAllCodimensionNeighbours(uint32_t idx, u32vector &neighbours, bvector &isghost)
Definition ParaTree.cpp:2918
bitpit::ParaTree::getOctant
Octant * getOctant(uint32_t idx)
Definition ParaTree.cpp:2469
bitpit::ParaTree::getMapping
void getMapping(uint32_t idx, u32vector &mapper, bvector &isghost) const
Definition ParaTree.cpp:3215
bitpit::ParaTree::getMaxLevel
int getMaxLevel() const
Definition ParaTree.cpp:1023
bitpit::ParaTree::getPersistentIdx
octantID getPersistentIdx(uint32_t idx) const
Definition ParaTree.cpp:1683
bitpit::ParaTree::restore
virtual void restore(std::istream &stream)
Definition ParaTree.cpp:588
bitpit::ParaTree::getBalance
bool getBalance(uint32_t idx) const
Definition ParaTree.cpp:1510
bitpit::ParaTree::getSerial
bool getSerial() const
Definition ParaTree.cpp:796
bitpit::ParaTree::getLocalMinSize
double getLocalMinSize() const
Definition ParaTree.cpp:2150
bitpit::ParaTree::loadBalance
void loadBalance(const dvector *weight=NULL)
Definition ParaTree.cpp:3997
bitpit::ParaTree::setNofGhostLayers
void setNofGhostLayers(std::size_t nofGhostLayers)
Definition ParaTree.cpp:2542
bitpit::ParaTree::getMorton
uint64_t getMorton(uint32_t idx) const
Definition ParaTree.cpp:1490
bitpit::ParaTree::computeConnectivity
void computeConnectivity()
Definition ParaTree.cpp:3815
bitpit::ParaTree::ParaTree
ParaTree(const std::string &logfile=DEFAULT_LOG_FILE, MPI_Comm comm=MPI_COMM_WORLD)
Definition ParaTree.cpp:136
bitpit::ParaTree::getFaceCenter
darray3 getFaceCenter(uint32_t idx, uint8_t face) const
Definition ParaTree.cpp:1361
bitpit::ParaTree::getNnodes
uint8_t getNnodes() const
Definition ParaTree.cpp:1039
bitpit::ParaTree::getLastDescMorton
uint64_t getLastDescMorton() const
Definition ParaTree.cpp:2188
bitpit::ParaTree::getPointOwnerRank
int getPointOwnerRank(const darray3 &point)
Definition ParaTree.cpp:3127
bitpit::ParaTree::getBound
bool getBound(uint32_t idx, uint8_t face) const
Definition ParaTree.cpp:1520
bitpit::ParaTree::findOwner
int findOwner(uint64_t morton) const
Definition ParaTree.cpp:3753
bitpit::ParaTree::~ParaTree
virtual ~ParaTree()
Definition ParaTree.cpp:211
bitpit::ParaTree::getRank
int getRank() const
Definition ParaTree.cpp:812
bitpit::ParaTree::DEFAULT_LOG_FILE
static BITPIT_PUBLIC_API const std::string DEFAULT_LOG_FILE
Definition ParaTree.hpp:119
bitpit::ParaTree::getIn
uint32_t getIn(const Intersection *inter) const
Definition ParaTree.cpp:2343
bitpit::ParaTree::getEdgeface
const uint8_t(* getEdgeface() const)[2]
Definition ParaTree.cpp:1204
bitpit::ParaTree::getLastOperation
Operation getLastOperation() const
Definition ParaTree.cpp:836
bitpit::ParaTree::getPartitionLastDesc
const std::vector< uint64_t > & getPartitionLastDesc() const
Definition ParaTree.cpp:975
bitpit::ParaTree::getMaxLength
uint32_t getMaxLength() const
Definition ParaTree.cpp:1031
bitpit::ParaTree::adaptGlobalRefine
bool adaptGlobalRefine(bool mapper_flag=false)
Definition ParaTree.cpp:3573
bitpit::ParaTree::getInternalOctantsBegin
octantIterator getInternalOctantsBegin()
Definition ParaTree.cpp:2205
bitpit::ParaTree::getSize
double getSize(uint32_t idx) const
Definition ParaTree.cpp:1311
bitpit::ParaTree::isCommSet
bool isCommSet() const
Definition ParaTree.cpp:936
bitpit::ParaTree::getPartitionFirstDesc
const std::vector< uint64_t > & getPartitionFirstDesc() const
Definition ParaTree.cpp:965
bitpit::ParaTree::evalLoadBalanceRanges
LoadBalanceRanges evalLoadBalanceRanges(dvector *weights)
Definition ParaTree.cpp:4089
bitpit::ParaTree::getY
double getY(uint32_t idx) const
Definition ParaTree.cpp:1293
bitpit::ParaTree::getEdgecoeffs
const int8_t(* getEdgecoeffs() const)[3]
Definition ParaTree.cpp:1220
bitpit::ParaTree::getPartitionRangeGlobalIdx
const std::vector< uint64_t > & getPartitionRangeGlobalIdx() const
Definition ParaTree.cpp:955
bitpit::ParaTree::reset
virtual void reset()
Definition ParaTree.cpp:448
bitpit::ParaTree::getPeriodic
bvector getPeriodic() const
Definition ParaTree.cpp:1228
bitpit::ParaTree::getZ0
double getZ0() const
Definition ParaTree.cpp:1007
bitpit::ParaTree::getPointOwnerIdx
uint32_t getPointOwnerIdx(const double *point) const
Definition ParaTree.cpp:3043
bitpit::ParaTree::setBalanceCodimension
void setBalanceCodimension(uint8_t b21codim)
Definition ParaTree.cpp:2237
bitpit::ParaTree::getNode
darray3 getNode(uint32_t idx, uint8_t node) const
Definition ParaTree.cpp:1385
bitpit::ParaTree::writeTest
void writeTest(const std::string &filename, dvector data)
Definition ParaTree.cpp:5871
bitpit::ParaTree::isNodeOnOctant
bool isNodeOnOctant(const Octant *nodeOctant, uint8_t nodeIndex, const Octant *octant) const
Definition ParaTree.cpp:3357
bitpit::ParaTree::getZ
double getZ(uint32_t idx) const
Definition ParaTree.cpp:1302
bitpit::ParaTree::getBordersPerProc
const std::map< int, std::vector< uint32_t > > & getBordersPerProc() const
Definition ParaTree.cpp:2560
bitpit::ParaTree::getY0
double getY0() const
Definition ParaTree.cpp:999
bitpit::ParaTree::getLoadBalanceRanges
const LoadBalanceRanges & getLoadBalanceRanges() const
Definition ParaTree.cpp:2526
bitpit::ParaTree::getComm
MPI_Comm getComm() const
Definition ParaTree.cpp:944
bitpit::ParaTree::setTol
void setTol(double tol=1.0e-14)
Definition ParaTree.cpp:1262
bitpit::ParaTree::getFiner
bool getFiner(const Intersection *inter) const
Definition ParaTree.cpp:2286
bitpit::ParaTree::computeNodePersistentKey
uint64_t computeNodePersistentKey(uint32_t idx, uint8_t node) const
Definition ParaTree.cpp:1501
bitpit::ParaTree::getArea
double getArea(uint32_t idx) const
Definition ParaTree.cpp:1320
bitpit::ParaTree::getL
double getL() const
Definition ParaTree.cpp:1015
bitpit::ParaTree::getIsGhost
bool getIsGhost(const Intersection *inter) const
Definition ParaTree.cpp:2304
bitpit::ParaTree::getIdx
uint32_t getIdx(const Octant *oct) const
Definition ParaTree.cpp:2035
bitpit::ParaTree::isEdgeOnOctant
bool isEdgeOnOctant(const Octant *edgeOctant, uint8_t edgeIndex, const Octant *octant) const
Definition ParaTree.cpp:3398
bitpit::ParaTree::getNofGhostLayers
std::size_t getNofGhostLayers() const
Definition ParaTree.cpp:2534
bitpit::ParaTree::findAllCodimensionNeighbours
void findAllCodimensionNeighbours(uint32_t idx, u32vector &neighbours, bvector &isghost)
Definition ParaTree.cpp:2872
bitpit::ParaTree::getFirstDescMorton
uint64_t getFirstDescMorton() const
Definition ParaTree.cpp:2181
bitpit::ParaTree::findGhostNeighbours
void findGhostNeighbours(uint32_t idx, uint8_t face, uint8_t codim, u32vector &neighbours) const
Definition ParaTree.cpp:2725
bitpit::ParaTree::getCoordinates
darray3 getCoordinates(uint32_t idx) const
Definition ParaTree.cpp:1275
bitpit::ParaTree::getPboundOctantsBegin
octantIterator getPboundOctantsBegin()
Definition ParaTree.cpp:2221
bitpit::ParaTree::getFace
uint8_t getFace(const Intersection *inter) const
Definition ParaTree.cpp:2322
bitpit::ParaTree::getMarker
int8_t getMarker(uint32_t idx) const
Definition ParaTree.cpp:1458
bitpit::ParaTree::getNumGhosts
uint32_t getNumGhosts() const
Definition ParaTree.cpp:2126
bitpit::ParaTree::getGhostConnectivity
const u32vector2D & getGhostConnectivity() const
Definition ParaTree.cpp:3895
bitpit::ParaTree::getNedges
uint8_t getNedges() const
Definition ParaTree.cpp:1055
bitpit::ParaTree::getNodeface
const uint8_t(* getNodeface() const)[3]
Definition ParaTree.cpp:1195
bitpit::ParaTree::getOwnerRank
int getOwnerRank(uint64_t globalIdx) const
Definition ParaTree.cpp:3797
bitpit::ParaTree::getInternalOctantsEnd
octantIterator getInternalOctantsEnd()
Definition ParaTree.cpp:2213
bitpit::ParaTree::getGlobalNumOctants
uint64_t getGlobalNumOctants() const
Definition ParaTree.cpp:788
bitpit::ParaTree::getNumNodes
uint32_t getNumNodes() const
Definition ParaTree.cpp:2134
bitpit::ParaTree::getOut
uint32_t getOut(const Intersection *inter) const
Definition ParaTree.cpp:2352
bitpit::ParaTree::updateConnectivity
void updateConnectivity()
Definition ParaTree.cpp:3832
bitpit::ParaTree::getNproc
int getNproc() const
Definition ParaTree.cpp:820
bitpit::ParaTree::getNfaces
uint8_t getNfaces() const
Definition ParaTree.cpp:1047
bitpit::ParaTree::getOutIsGhost
bool getOutIsGhost(const Intersection *inter) const
Definition ParaTree.cpp:2361
bitpit::ParaTree::getLocalMaxDepth
uint8_t getLocalMaxDepth() const
Definition ParaTree.cpp:2142
bitpit::ParaTree::getGhostLayer
int getGhostLayer(const Octant *oct) const
Definition ParaTree.cpp:2518
bitpit::ParaTree::getVolume
double getVolume(uint32_t idx) const
Definition ParaTree.cpp:1329
bitpit::ParaTree::getDumpVersion
virtual int getDumpVersion() const
Definition ParaTree.cpp:488
bitpit::ParaTree::getNodecoeffs
const int8_t(* getNodecoeffs() const)[3]
Definition ParaTree.cpp:1212
bitpit::ParaTree::getLevel
uint8_t getLevel(uint32_t idx) const
Definition ParaTree.cpp:1481
bitpit::ParaTree::getOppface
const uint8_t * getOppface() const
Definition ParaTree.cpp:1177
bitpit::ParaTree::clearConnectivity
void clearConnectivity(bool release=true)
Definition ParaTree.cpp:3825
bitpit::ParaTree::getX
double getX(uint32_t idx) const
Definition ParaTree.cpp:1284
bitpit::ParaTree::computeIntersections
void computeIntersections()
Definition ParaTree.cpp:4356
bitpit::ParaTree::replaceComm
void replaceComm(MPI_Comm communicator, MPI_Comm *previousCommunicator=nullptr)
Definition ParaTree.cpp:878
bitpit::ParaTree::getIsNewC
bool getIsNewC(uint32_t idx) const
Definition ParaTree.cpp:1566
bitpit::ParaTree::getGlobalIdx
uint64_t getGlobalIdx(uint32_t idx) const
Definition ParaTree.cpp:1575
bitpit::ParaTree::getCenter
void getCenter(uint32_t idx, darray3 &centerCoords) const
Definition ParaTree.cpp:1338
bitpit::ParaTree::getNodeCoordinates
darray3 getNodeCoordinates(uint32_t node) const
Definition ParaTree.cpp:3886
bitpit::ParaTree::getOwners
u32vector getOwners(const Intersection *inter) const
Definition ParaTree.cpp:2331
bitpit::ParaTree::setPeriodic
void setPeriodic(uint8_t i)
Definition ParaTree.cpp:1252
bitpit::ParaTree::getConnectivity
const u32vector2D & getConnectivity() const
Definition ParaTree.cpp:3841
bitpit::ParaTree::setBalance
void setBalance(uint32_t idx, bool balance)
Definition ParaTree.cpp:1709
bitpit::ParaTree::getPreMarker
int8_t getPreMarker(uint32_t idx)
Definition ParaTree.cpp:1468
bitpit::ParaTree::dump
virtual void dump(std::ostream &stream, bool full=true)
Definition ParaTree.cpp:502
bitpit::ParaTree::getGhostGlobalIdx
uint64_t getGhostGlobalIdx(uint32_t idx) const
Definition ParaTree.cpp:1650
bitpit::ParaTree::setComm
void setComm(MPI_Comm communicator)
Definition ParaTree.cpp:848
bitpit::ParaTree::getBalanceCodimension
uint8_t getBalanceCodimension() const
Definition ParaTree.cpp:2174
bitpit::ParaTree::getNormal
void getNormal(uint32_t idx, uint8_t face, darray3 &normal) const
Definition ParaTree.cpp:1433
bitpit::ParaTree::getPboundOctantsEnd
octantIterator getPboundOctantsEnd()
Definition ParaTree.cpp:2229
bitpit::ParaTree::getX0
double getX0() const
Definition ParaTree.cpp:991
bitpit::ParaTree::getStatus
uint64_t getStatus() const
Definition ParaTree.cpp:2110
bitpit::ParaTree::getIntersection
Intersection * getIntersection(uint32_t idx)
Definition ParaTree.cpp:2262
bitpit::ParaTree::getGhostOctant
Octant * getGhostOctant(uint32_t idx)
Definition ParaTree.cpp:2489
bitpit::ParaTree::getLocalMaxSize
double getLocalMaxSize() const
Definition ParaTree.cpp:2159
bitpit::ParaTree::getOrigin
darray3 getOrigin() const
Definition ParaTree.cpp:983
bitpit::ParaTree::getIsNewR
bool getIsNewR(uint32_t idx) const
Definition ParaTree.cpp:1557
abs
std::array< T, d > abs(const std::array< T, d > &x)
Definition MathOperators_array.tpp:591
sum
void sum(const std::array< T, d > &x, T1 &s)
Definition MathOperators_array.tpp:544
sign
T sign(const T &)
Definition Operators.tpp:38
max
std::array< T, d > max(const std::array< T, d > &x, const std::array< T, d > &y)
Definition MathOperators_array.tpp:306
pow
std::array< T, d > pow(std::array< T, d > &x, double p)
Definition MathOperators_array.tpp:637
uipow
T uipow(const T &, unsigned int)
Definition Operators.tpp:55
min
std::array< T, d > min(const std::array< T, d > &x, const std::array< T, d > &y)
Definition MathOperators_array.tpp:67
bitpit::ParaTree::ExchangeRanges
std::unordered_map< int, std::array< uint32_t, 2 > > ExchangeRanges
Definition ParaTree.hpp:121
bitpit::utils::binary::write
void write(std::ostream &stream, const std::vector< bool > &container)
Definition binaryUtils.cpp:41
bitpit::utils::binary::read
void read(std::istream &stream, std::vector< bool > &container)
Definition binaryUtils.cpp:72
BITPIT_UNUSED
#define BITPIT_UNUSED(variable)
Definition compiler.hpp:63
bitpit::log::cout
Logger & cout(log::Level defaultSeverity, log::Visibility defaultVisibility)
Definition logger.cpp:1705
bitpit::log::manager
LoggerManager & manager()
Definition logger.cpp:1685
bitpit::log::context
LoggerManipulator< std::string > context(const std::string &context)
Definition logger.cpp:1929
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