ParaTree.hpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit is distributed in the hope that it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
 \*---------------------------------------------------------------------------*/
 
#ifndef __BITPIT_PABLO_PARA_TREE_HPP__
#define __BITPIT_PABLO_PARA_TREE_HPP__
 
// =================================================================================== //
// INCLUDES                                                                            //
// =================================================================================== //
 
// Octant header should be included before bitpit_communications header. If
// this is not the case, octants will be streamed in the communication buffer
// using the standard bitpit operator instead of the custom operator defined
// for the octant. GCC will work also if the Octant header is included after
// bitpit_communications header, but that's a bug in GCC.
//
//    https://gcc.gnu.org/bugzilla/show_bug.cgi?id=51577
//    https://bugs.llvm.org/show_bug.cgi?id=47967
//
// For a dependent name used in a template definition, the lookup is postponed
// until the template arguments are known, at which time ADL examines function
// declarations with external linkage (until C++11) that are visible from the
// template definition context as well as in the template instantiation context
// (in other words, ADL will only look in namespaces associated with the
// argument types), while unqualified lookup only examines function declarations
// with external linkage (until C++11) that are visible from the template
// definition context (in other words, adding a new function declaration after
// template definition does not make it visible except via ADL). The stream
// operator template is defined in the global namespace, whereas the octant
// class is defined in the bitpit namespace. Hence, custom stream operator
// defined for the octant class cannot be selected by ADL (ADL will only look
// in the namaspace of its arguments, i.e., bitpit, whereas the custom stream
// operator is defined in the global namespace). The only chance to select
// the custom operator is through unqualified lookup, and for that to happen,
// the stream operator defined for the octant class should be defined before
// the template.
 
#include "Octant.hpp"
 
#if BITPIT_ENABLE_MPI==1
#include <mpi.h>
#include "DataLBInterface.hpp"
#include "DataCommInterface.hpp"
#include "bitpit_communications.hpp"
#endif
#include "tree_constants.hpp"
#include "LocalTree.hpp"
#include "Map.hpp"
#include <map>
#include <unordered_map>
#include <unordered_set>
#include <set>
#include <bitset>
#include <algorithm>
 
#include "bitpit_common.hpp"
 
namespace bitpit {
 
    // =================================================================================== //
    // TYPEDEFS                                                                            //
    // =================================================================================== //
    typedef std::vector<bool>               bvector;
    typedef std::vector<int>                ivector;
    typedef std::bitset<72>                 octantID;
    typedef std::vector<Octant*>            ptroctvector;
    typedef ptroctvector::iterator          octantIterator;
 
    // =================================================================================== //
    // CLASS DEFINITION                                                                    //
    // =================================================================================== //

    class ParaTree{
 
        // =================================================================================== //
        // MEMBERS                                                                             //
        // =================================================================================== //
    public:
        BITPIT_PUBLIC_API static const std::string  DEFAULT_LOG_FILE;           
 
        typedef std::unordered_map<int, std::array<uint32_t, 2>> ExchangeRanges;
 
        enum Operation {
            OP_NONE,
            OP_INIT,
            OP_PRE_ADAPT,
            OP_ADAPT_MAPPED,
            OP_ADAPT_UNMAPPED,
            OP_LOADBALANCE_FIRST,
            OP_LOADBALANCE
        };
 
        struct LoadBalanceRanges {
            enum ExchangeAction {
                ACTION_UNDEFINED = -1,
                ACTION_NONE,
                ACTION_SEND,
                ACTION_DELETE,
                ACTION_RECEIVE
            };
 
            ExchangeAction sendAction;
            ExchangeRanges sendRanges;
 
            ExchangeAction recvAction;
            ExchangeRanges recvRanges;
 
            LoadBalanceRanges();
            LoadBalanceRanges(bool serial, const ExchangeRanges &_sendRanges, const ExchangeRanges &_recvRanges);
 
            void clear();
        };
 
    private:
        typedef std::unordered_map<int, std::array<uint64_t, 2>> PartitionIntersections;
 
        struct AccretionData {
            int targetRank;
            std::unordered_map<uint64_t, int> internalSeeds;
            std::unordered_map<uint64_t, int> foreignSeeds;
            std::unordered_map<uint64_t, int> population;
        };
 
        //undistributed members
        std::vector<uint64_t>   m_partitionFirstDesc;           
        std::vector<uint64_t>   m_partitionLastDesc;            
        std::vector<uint64_t>   m_partitionRangeGlobalIdx;      
        std::vector<uint64_t>   m_partitionRangeGlobalIdx0;     
        uint64_t                m_globalNumOctants;             
        int                     m_nproc;                        
        int8_t                  m_maxDepth;                     
        const TreeConstants    *m_treeConstants;                
        std::size_t             m_nofGhostLayers;               
 
        //distributed members
        int                     m_rank;                         
        LocalTree               m_octree;                       
        std::map<int,u32vector> m_bordersPerProc;               
        ptroctvector            m_internals;                    
        ptroctvector            m_pborders;                     
 
        //distributed adpapting memebrs
        u32vector               m_mapIdx;                       
        //elements sent during last loadbalance operation
        LoadBalanceRanges         m_loadBalanceRanges;                                            
 
        //auxiliary members
        int                     m_errorFlag;                    
        bool                    m_serial;                       
        double                  m_tol;                          
 
        //map members
        Map                     m_trans;                        
        uint8_t                 m_dim;                          
 
        //boundary conditions members
        bvector                 m_periodic;                     
 
        //info member
        uint64_t                m_status;                       
        Operation               m_lastOp;                       
 
        //log member
        Logger*                 m_log;                          
 
        //communicator
#if BITPIT_ENABLE_MPI==1
        //TODO Duplicate communicator
        MPI_Comm                m_comm;                         
#endif
 
        // =================================================================================== //
        // CONSTRUCTORS AND OPERATORS                                                          //
        // =================================================================================== //
    public:
#if BITPIT_ENABLE_MPI==1
        ParaTree(const std::string &logfile = DEFAULT_LOG_FILE, MPI_Comm comm = MPI_COMM_WORLD);
        ParaTree(uint8_t dim, const std::string &logfile = DEFAULT_LOG_FILE, MPI_Comm comm = MPI_COMM_WORLD);
        ParaTree(std::istream &stream, const std::string &logfile = DEFAULT_LOG_FILE, MPI_Comm comm = MPI_COMM_WORLD);
#else
        ParaTree(const std::string &logfile = DEFAULT_LOG_FILE);
        ParaTree(uint8_t dim, const std::string &logfile = DEFAULT_LOG_FILE);
        ParaTree(std::istream &stream, const std::string &logfile = DEFAULT_LOG_FILE);
#endif
        virtual ~ParaTree();
 
        ParaTree(const ParaTree & other);
 
    private:
        ParaTree & operator=(const ParaTree & other) = delete;
 
        // =================================================================================== //
        // METHODS                                                                             //
        // =================================================================================== //
    public:
        virtual void    reset();
 
        virtual int     getDumpVersion() const;
        virtual void    dump(std::ostream &stream, bool full = true);
        virtual void    restore(std::istream &stream);
 
        void    printHeader();
 
        // =================================================================================== //
        // BASIC GET/SET METHODS                                                               //
        // =================================================================================== //
        uint8_t     getDim() const;
        uint64_t    getGlobalNumOctants() const;
        bool        getSerial() const;
        bool        getParallel() const;
        int         getRank() const;
        int         getNproc() const;
        Logger&     getLog();
        Operation   getLastOperation() const;
#if BITPIT_ENABLE_MPI==1
        void        setComm(MPI_Comm communicator);
        void        replaceComm(MPI_Comm communicator, MPI_Comm *previousCommunicator = nullptr);
        void        freeComm();
        MPI_Comm    getComm() const;
        bool        isCommSet() const;
#endif
        const std::vector<uint64_t> &getPartitionRangeGlobalIdx() const;
        const std::vector<uint64_t> &getPartitionFirstDesc() const;
        const std::vector<uint64_t> &getPartitionLastDesc() const;
        darray3     getOrigin() const;
        double      getX0() const;
        double      getY0() const;
        double      getZ0() const;
        double      getL() const;
        int         getMaxLevel() const;
        uint32_t    getMaxLength() const;
        uint8_t     getNnodes() const;
        uint8_t     getNfaces() const;
        uint8_t     getNedges() const;
        uint8_t     getNchildren() const;
        uint8_t     getNnodesperface() const;
        void        getNormals(int8_t normals[6][3]) const;
        void        getOppface(uint8_t oppface[6]) const;
        void        getFacenode(uint8_t facenode[6][4]) const;
        void        getNodeface(uint8_t nodeface[8][3]) const;
        void        getEdgeface(uint8_t edgeface[12][2]) const;
        void        getNodecoeffs(int8_t nodecoeffs[8][3]) const;
        void        getEdgecoeffs(int8_t edgecoeffs[12][3]) const;
        const int8_t    (*getNormals() const)[3];
        const uint8_t   *getOppface() const;
        const uint8_t   (*getFacenode() const)[4];
        const uint8_t   (*getNodeface() const)[3];
        const uint8_t   (*getEdgeface() const)[2];
        const int8_t    (*getNodecoeffs() const)[3];
        const int8_t    (*getEdgecoeffs() const)[3];
        bvector     getPeriodic() const;
        double      getTol() const;
        bool        getPeriodic(uint8_t i) const;
        void        setPeriodic(uint8_t i);
        void        setTol(double tol = 1.0e-14);
 
        // =================================================================================== //
        // INDEX BASED METHODS                                                                 //
        // =================================================================================== //
        darray3     getCoordinates(uint32_t idx) const;
        double      getX(uint32_t idx) const;
        double      getY(uint32_t idx) const;
        double      getZ(uint32_t idx) const;
        double      getSize(uint32_t idx) const;
        double      getArea(uint32_t idx) const;
        double      getVolume(uint32_t idx) const;
        void        getCenter(uint32_t idx, darray3& centerCoords) const;
        darray3     getCenter(uint32_t idx) const;
        darray3     getFaceCenter(uint32_t idx, uint8_t face) const;
        void        getFaceCenter(uint32_t idx, uint8_t face, darray3& centerCoords) const;
        darray3     getNode(uint32_t idx, uint8_t node) const;
        void        getNode(uint32_t idx, uint8_t node, darray3& nodeCoords) const;
        void        getNodes(uint32_t idx, darr3vector & nodesCoords) const;
        darr3vector getNodes(uint32_t idx) const;
        void        getNormal(uint32_t idx, uint8_t face, darray3 & normal) const;
        darray3     getNormal(uint32_t idx, uint8_t face) const;
        int8_t      getMarker(uint32_t idx) const;
        uint8_t     getLevel(uint32_t idx) const;
        uint64_t    getMorton(uint32_t idx) const;
        uint64_t    computeNodePersistentKey(uint32_t idx, uint8_t node) const;
        bool        getBalance(uint32_t idx) const;
        bool        getBound(uint32_t idx, uint8_t face) const;
        bool        getBound(uint32_t idx) const;
        bool        getPbound(uint32_t idx, uint8_t face) const;
        bool        getPbound(uint32_t idx) const;
        bool        getIsNewR(uint32_t idx) const;
        bool        getIsNewC(uint32_t idx) const;
        uint64_t    getGlobalIdx(uint32_t idx) const;
        uint64_t    getGhostGlobalIdx(uint32_t idx) const;
        uint32_t    getLocalIdx(uint64_t gidx) const;
        uint32_t    getLocalIdx(uint64_t gidx,int rank) const;
        void        getLocalIdx(uint64_t gidx,uint32_t & lidx,int & rank) const;
        uint32_t    getGhostLocalIdx(uint64_t gidx) const;
        octantID    getPersistentIdx(uint32_t idx) const;
        int8_t      getPreMarker(uint32_t idx);
        void        setMarker(uint32_t idx, int8_t marker);
        void        setBalance(uint32_t idx, bool balance);
 
        // =================================================================================== //
        // POINTER BASED METHODS                                                               //
        // =================================================================================== //
        darray3     getCoordinates(const Octant* oct) const;
        double      getX(const Octant* oct) const;
        double      getY(const Octant* oct) const;
        double      getZ(const Octant* oct) const;
        double      getSize(const Octant* oct) const;
        double      getArea(const Octant* oct) const;
        double      getVolume(const Octant* oct) const;
        void        getCenter(const Octant* oct, darray3& centerCoords) const;
        darray3     getCenter(const Octant* oct) const;
        darray3     getFaceCenter(const Octant* oct, uint8_t face) const;
        void        getFaceCenter(const Octant* oct, uint8_t face, darray3& centerCoords) const;
        darray3     getNode(const Octant* oct, uint8_t node) const;
        void        getNode(const Octant* oct, uint8_t node, darray3& nodeCoords) const;
        void        getNodes(const Octant* oct, darr3vector & nodesCoords) const;
        darr3vector getNodes(const Octant* oct) const;
        void        getNormal(const Octant* oct, uint8_t face, darray3 & normal) const;
        darray3     getNormal(const Octant* oct, uint8_t face) const;
        int8_t      getMarker(const Octant* oct) const;
        uint8_t     getLevel(const Octant* oct) const;
        uint64_t    getMorton(const Octant* oct) const;
        uint64_t    getLastDescMorton(const Octant* oct) const;
        uint64_t    computeNodePersistentKey(const Octant* oct, uint8_t node) const;
        bool        getBalance(const Octant* oct) const;
        bool        getBound(const Octant* oct, uint8_t face) const;
        bool        getBound(const Octant* oct) const;
        bool        getPbound(const Octant* oct, uint8_t face) const;
        bool        getPbound(const Octant* oct) const;
        bool        getIsNewR(const Octant* oct) const;
        bool        getIsNewC(const Octant* oct) const;
        uint32_t    getIdx(const Octant* oct) const;
        uint64_t    getGlobalIdx(const Octant* oct) const;
        octantID    getPersistentIdx(const Octant* oct) const;
        int8_t      getPreMarker(Octant* oct);
        void        setMarker(Octant* oct, int8_t marker);
        void        setBalance(Octant* oct, bool balance);
 
        // =================================================================================== //
        // LOCAL TREE GET/SET METHODS                                                          //
        // =================================================================================== //
        uint64_t    getStatus() const;
        uint32_t    getNumOctants() const;
        uint32_t    getNumGhosts() const;
        uint32_t    getNumNodes() const;
        uint8_t     getLocalMaxDepth() const;
        double      getLocalMaxSize() const;
        double      getLocalMinSize() const;
        uint8_t     getBalanceCodimension() const;
        uint64_t    getFirstDescMorton() const;
        uint64_t    getLastDescMorton() const;
        uint64_t    getLastDescMorton(uint32_t idx) const;
        octantIterator  getInternalOctantsBegin();
        octantIterator  getInternalOctantsEnd();
        octantIterator  getPboundOctantsBegin();
        octantIterator  getPboundOctantsEnd();
        void        setBalanceCodimension(uint8_t b21codim);
 
        // =================================================================================== //
        // INTERSECTION GET/SET METHODS                                                        //
        // =================================================================================== //
        uint32_t    getNumIntersections() const;
        Intersection* getIntersection(uint32_t idx);
        uint8_t     getLevel(const Intersection* inter) const;
        bool        getFiner(const Intersection* inter) const;
        bool        getBound(const Intersection* inter) const;
        bool        getIsGhost(const Intersection* inter) const;
        bool        getPbound(const Intersection* inter) const;
        uint8_t     getFace(const Intersection* inter) const;
        u32vector   getOwners(const Intersection* inter) const;
        uint32_t    getIn(const Intersection* inter) const;
        uint32_t    getOut(const Intersection* inter) const;
        bool        getOutIsGhost(const Intersection* inter) const;
        double      getSize(const Intersection* inter) const;
        double      getArea(const Intersection* inter) const;
        darray3     getCenter(const Intersection* inter) const;
        darr3vector getNodes(const Intersection* inter) const;
        darray3     getNormal(const Intersection* inter) const;
 
        // =================================================================================== //
        // OTHER GET/SET METHODS                                                               //
        // =================================================================================== //
        Octant* getOctant(uint32_t idx);
        const Octant* getOctant(uint32_t idx) const;
        Octant* getGhostOctant(uint32_t idx);
        const Octant* getGhostOctant(uint32_t idx) const;
        bool getIsGhost(const Octant* oct) const;
        int getGhostLayer(const Octant* oct) const;
        const LoadBalanceRanges & getLoadBalanceRanges() const;
        std::size_t getNofGhostLayers() const;
        void setNofGhostLayers(std::size_t nofGhostLayers);
        const std::map<int, std::vector<uint32_t>> & getBordersPerProc() const;
 
        // =================================================================================== //
        // PRIVATE GET/SET METHODS                                                             //
        // =================================================================================== //
    private:
        void        setDim(uint8_t dim);
#if BITPIT_ENABLE_MPI==1
        void        updateGlobalFirstDescMorton();
        void        updateGlobalLasttDescMorton();
#endif
 
        // =================================================================================== //
        // OTHER METHODS                                                                       //
        // =================================================================================== //
 
#if BITPIT_ENABLE_MPI==1
        void    _initializeCommunicator(MPI_Comm comm);
#else
        void    _initializeSerialCommunicator();
#endif
        void    _initializePartitions();
        void    _initialize(uint8_t dim, const std::string &logfile);
 
#if BITPIT_ENABLE_MPI==1
        void    initialize(const std::string &logfile, MPI_Comm comm);
        void    initialize(uint8_t dim, const std::string &logfile, MPI_Comm comm);
#else
        void    initialize(const std::string &logfile);
        void    initialize(uint8_t dim, const std::string &logfile);
#endif
        void    initializeLogger(const std::string &logfile);
        void    reinitialize(uint8_t dim, const std::string &logfile);
 
        void    reset(bool createRoot);
 
        uint64_t    evalPointAnchorMorton(const double * point) const;
 
        // =================================================================================== //
        // OTHER OCTANT BASED METHODS                                                          //
        // =================================================================================== //
 
        void        findNeighbours(const Octant* oct, uint8_t face, uint8_t codim, u32vector & neighbours, bvector & isghost, bool onlyinternals, bool append) const;
    public:
        void        findNeighbours(uint32_t idx, uint8_t face, uint8_t codim, u32vector & neighbours, bvector & isghost) const;
        void        findNeighbours(uint32_t idx, uint8_t face, uint8_t codim, u32vector & neighbours) const;
        void        findNeighbours(const Octant* oct, uint8_t face, uint8_t codim, u32vector & neighbours, bvector & isghost) const ;
        void        findGhostNeighbours(uint32_t idx, uint8_t face, uint8_t codim, u32vector & neighbours) const;
        void        findGhostNeighbours(uint32_t idx, uint8_t face, uint8_t codim, u32vector & neighbours, bvector & isghost) const;
        void        findGhostNeighbours(const Octant* oct, uint8_t face, uint8_t codim, u32vector & neighbours, bvector & isghost) const;
        void        findAllNodeNeighbours(uint32_t idx, uint32_t node, u32vector & neighbours, bvector & isghost);
        void        findAllNodeNeighbours(const Octant* oct, uint32_t node, u32vector & neighbours, bvector & isghost) const;
        void        findAllCodimensionNeighbours(uint32_t idx, u32vector & neighbours, bvector & isghost);
        void        findAllCodimensionNeighbours(const Octant* oct, u32vector & neighbours, bvector & isghost);
        void        findGhostAllCodimensionNeighbours(uint32_t idx, u32vector & neighbours, bvector & isghost);
        void        findGhostAllCodimensionNeighbours(Octant* oct, u32vector & neighbours, bvector & isghost);
        Octant*     getPointOwner(const dvector &point);
        Octant*     getPointOwner(const dvector &point, bool & isghost);
        Octant*     getPointOwner(const darray3 &point);
        Octant*     getPointOwner(const darray3 &point, bool & isghost);
        uint32_t    getPointOwnerIdx(const double * point) const;
        uint32_t    getPointOwnerIdx(const double * point, bool & isghost) const;
        uint32_t    getPointOwnerIdx(const dvector &point) const;
        uint32_t    getPointOwnerIdx(const dvector &point, bool & isghost) const;
        uint32_t    getPointOwnerIdx(const darray3 &point) const;
        uint32_t    getPointOwnerIdx(const darray3 &point, bool & isghost) const;
        void        getMapping(uint32_t idx, u32vector & mapper, bvector & isghost) const;
        void        getMapping(uint32_t idx, u32vector & mapper, bvector & isghost, ivector & rank) const;
        void        getPreMapping(u32vector & idx, std::vector<int8_t> & markers, std::vector<bool> & isghost);
        bool        isNodeOnOctant(const Octant* nodeOctant, uint8_t nodeIndex, const Octant* octant) const;
        bool        isEdgeOnOctant(const Octant* edgeOctant, uint8_t edgeIndex, const Octant* octant) const;
        bool        isFaceOnOctant(const Octant* faceOctant, uint8_t faceIndex, const Octant* octant) const;
        int         getPointOwnerRank(const darray3 &point);
        uint8_t     getFamilySplittingNode(const Octant*) const;
        void        expectedOctantAdapt(const Octant* octant, int8_t marker, octvector* result) const;
 
        // =================================================================================== //
        // OTHER PARATREE BASED METHODS                                                        //
        // =================================================================================== //
        int8_t      getMaxDepth() const;
        int         findOwner(uint64_t morton) const;
        int         getOwnerRank(uint64_t globalIdx) const;
        void        settleMarkers();
        void        preadapt();
        bool        checkToAdapt();
        bool        adapt(bool mapper_flag = false);
        bool        adaptGlobalRefine(bool mapper_flag = false);
        bool        adaptGlobalCoarse(bool mapper_flag = false);
        void        computeConnectivity();
        void        clearConnectivity(bool release = true);
        void        updateConnectivity();
        const u32vector2D & getConnectivity() const;
        const u32vector & getConnectivity(uint32_t idx) const;
        const u32vector & getConnectivity(Octant* oct) const;
        const u32arr3vector & getNodes() const;
        const u32array3 & getNodeLogicalCoordinates(uint32_t node) const;
        darray3     getNodeCoordinates(uint32_t node) const;
        const u32vector2D & getGhostConnectivity() const;
        const u32vector & getGhostConnectivity(uint32_t idx) const;
        const u32vector & getGhostConnectivity(const Octant* oct) const;
        bool        check21Balance();
#if BITPIT_ENABLE_MPI==1
        void        loadBalance(const dvector* weight = NULL);
        void        loadBalance(uint8_t & level, const dvector* weight = NULL);
 
        LoadBalanceRanges evalLoadBalanceRanges(dvector *weights);
        LoadBalanceRanges evalLoadBalanceRanges(uint8_t level, dvector *weights);
    private:
        LoadBalanceRanges evalLoadBalanceRanges(const uint32_t *updatedPartition);
 
        ExchangeRanges evalLoadBalanceSendRanges(const uint32_t *updatedPartition);
        ExchangeRanges evalLoadBalanceRecvRanges(const uint32_t *updatedPartition);
 
        PartitionIntersections evalPartitionIntersections(const uint32_t *schema_A, int rank_A, const uint32_t *schema_B);
#endif
    public:
        double      levelToSize(uint8_t level) const;
 
        // =================================================================================== //
        // OTHER INTERSECTION BASED METHODS                                                    //
        // =================================================================================== //
        void        computeIntersections();
 
        // =================================================================================== //
        // OTHER PRIVATE METHODS                                                               //
        // =================================================================================== //
    private:
        Octant& extractOctant(uint32_t idx);
        bool        private_adapt_mapidx(bool mapflag);
        void        updateAdapt();
#if BITPIT_ENABLE_MPI==1
        void        computePartition(uint32_t *partition);
        void        computePartition(const dvector *weight, uint32_t *partition);
        void        computePartition(uint8_t level_, const dvector *weight, uint32_t *partition);
        void        updateLoadBalance();
        void        setPboundGhosts();
        void        buildGhostOctants(const std::map<int, u32vector> &bordersPerProc, const std::vector<AccretionData> &accretions);
        void        initializeGhostHaloAccretions(std::vector<AccretionData> *accretions);
        void        growGhostHaloAccretions(std::unordered_map<uint32_t, std::vector<uint64_t>> *cachedOneRings, std::vector<AccretionData> *accretions);
        void        exchangeGhostHaloAccretions(DataCommunicator *dataCommunicator, std::vector<AccretionData> *accretions);
 
        void        computeGhostHalo();
        bool        commMarker();
#endif
        void        updateAfterCoarse();
        void        balance21(bool verbose, bool balanceNewOctants);
        void        createPartitionInfo();
 
        bool        isInternal(uint64_t gidx) const;
 
        // =================================================================================== //
        // TESTING OUTPUT METHODS                                                                  //
        // =================================================================================== //
    public:
        void        write(const std::string &filename);
        void        writeTest(const std::string &filename, dvector data);
 
        // =================================================================================== //
        // TEMPLATE METHODS                                                                    //
        // =================================================================================== //
#if BITPIT_ENABLE_MPI==1

        template<class Impl>
        void
        communicate(DataCommInterface<Impl> & userData){
            //BUILD SEND BUFFERS
            DataCommunicator communicator(m_comm);
            size_t fixedDataSize = userData.fixedSize();
            std::map<int,u32vector >::iterator bitend = m_bordersPerProc.end();
            std::map<int,u32vector >::iterator bitbegin = m_bordersPerProc.begin();
            for(std::map<int,u32vector >::iterator bit = bitbegin; bit != bitend; ++bit){
                int  key = bit->first;
                const u32vector & pborders = bit->second;
                size_t buffSize = 0;
                size_t nofPbordersPerProc = pborders.size();
                if(fixedDataSize != 0){
                    buffSize = fixedDataSize*nofPbordersPerProc;
                }
                else{
                    for(size_t i = 0; i < nofPbordersPerProc; ++i){
                        buffSize += userData.size(pborders[i]);
                    }
                }
                //enlarge buffer to store number of pborders from this proc
                buffSize += sizeof(size_t);
                //build buffer for this proc
                communicator.setSend(key,buffSize);
                SendBuffer & sendBuffer = communicator.getSendBuffer(key);
                //store number of pborders from this proc at the begining
                sendBuffer << nofPbordersPerProc;
 
                //WRITE SEND BUFFERS
                for(size_t j = 0; j < nofPbordersPerProc; ++j){
                    userData.gather(sendBuffer,pborders[j]);
                }
            }
 
            communicator.discoverRecvs();
            communicator.startAllRecvs();
 
            communicator.startAllSends();
 
            //READ RECEIVE BUFFERS
            int ghostOffset = 0;
            std::vector<int> recvRanks = communicator.getRecvRanks();
            std::sort(recvRanks.begin(),recvRanks.end());
            for(int rank : recvRanks){
                communicator.waitRecv(rank);
                RecvBuffer & recvBuffer = communicator.getRecvBuffer(rank);
                size_t nofGhostFromThisProc = 0;
                recvBuffer >> nofGhostFromThisProc;
                for(size_t k = 0; k < nofGhostFromThisProc; ++k){
                    userData.scatter(recvBuffer, k+ghostOffset);
                }
                ghostOffset += nofGhostFromThisProc;
            }
            communicator.waitAllSends();
        }

        template<class Impl>
        void
        loadBalance(DataLBInterface<Impl> & userData, dvector* weight = NULL){
            //Write info on log
            (*m_log) << "---------------------------------------------" << std::endl;
            (*m_log) << " LOAD BALANCE " << std::endl;
 
            m_lastOp = OP_LOADBALANCE;
            if (m_nproc>1){
 
                std::vector<uint32_t> partition(m_nproc);
                if (weight == NULL)
                    computePartition(partition.data());
                else
                    computePartition(weight, partition.data());
 
                weight = NULL;
 
                privateLoadBalance(partition.data(), &userData);
 
                //Write info of final partition on log
                (*m_log) << " " << std::endl;
                (*m_log) << " Final Parallel partition : " << std::endl;
                (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(0))+" :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << std::endl;
                for(int ii=1; ii<m_nproc; ii++){
                    (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(ii))+"    :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[ii]-m_partitionRangeGlobalIdx[ii-1])) << std::endl;
                }
                (*m_log) << " " << std::endl;
                (*m_log) << "---------------------------------------------" << std::endl;
 
            }
            else{
                m_loadBalanceRanges.clear();
 
                (*m_log) << " " << std::endl;
                (*m_log) << " Serial partition : " << std::endl;
                (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(0))+" :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << std::endl;
                (*m_log) << " " << std::endl;
                (*m_log) << "---------------------------------------------" << std::endl;
            }
 
        }

        template<class Impl>
        void
        loadBalance(DataLBInterface<Impl> & userData, uint8_t & level, dvector* weight = NULL){
 
            //Write info on log
            (*m_log) << "---------------------------------------------" << std::endl;
            (*m_log) << " LOAD BALANCE " << std::endl;
 
            m_lastOp = OP_LOADBALANCE;
            if (m_nproc>1){
 
                std::vector<uint32_t> partition(m_nproc);
                computePartition(level, weight, partition.data());
 
                privateLoadBalance(partition.data(), &userData);
 
                //Write info of final partition on log
                (*m_log) << " " << std::endl;
                (*m_log) << " Final Parallel partition : " << std::endl;
                (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(0))+" :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << std::endl;
                for(int ii=1; ii<m_nproc; ii++){
                    (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(ii))+"    :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[ii]-m_partitionRangeGlobalIdx[ii-1])) << std::endl;
                }
                (*m_log) << " " << std::endl;
                (*m_log) << "---------------------------------------------" << std::endl;
 
            }
            else{
                m_loadBalanceRanges.clear();
 
                (*m_log) << " " << std::endl;
                (*m_log) << " Serial partition : " << std::endl;
                (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(0))+" :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << std::endl;
                (*m_log) << " " << std::endl;
                (*m_log) << "---------------------------------------------" << std::endl;
            }
 
        }

        template<class Impl>
        void
        privateLoadBalance(const uint32_t *partition, DataLBInterface<Impl> *userData = nullptr){
 
            (*m_log) << " " << std::endl;
            if (m_serial) {
                (*m_log) << " Initial Serial distribution : " << std::endl;
            } else {
                (*m_log) << " Initial Parallel partition : " << std::endl;
            }
            (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(0))+" :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[0]+1)) << std::endl;
            for(int ii=1; ii<m_nproc; ii++){
                (*m_log) << " Octants for proc  "+ std::to_string(static_cast<unsigned long long>(ii))+"    :   " + std::to_string(static_cast<unsigned long long>(m_partitionRangeGlobalIdx[ii]-m_partitionRangeGlobalIdx[ii-1])) << std::endl;
            }
 
            // Compute load balance ranges
            std::unordered_map<int, std::array<uint32_t, 2>> sendRanges = evalLoadBalanceSendRanges(partition);
            std::unordered_map<int, std::array<uint32_t, 2>> recvRanges = evalLoadBalanceRecvRanges(partition);
 
            m_loadBalanceRanges = LoadBalanceRanges(m_serial, sendRanges, recvRanges);
 
            // Compute information about the new partitioning
            assert(m_nproc > 0);
 
            uint32_t newSizeOctants = partition[m_rank];
 
            uint64_t newLastOctantGlobalIdx;
            uint64_t newFirstOctantGlobalIdx;
            if (newSizeOctants > 0) {
                newLastOctantGlobalIdx = partition[0];
                for (int p = 1; p < m_rank + 1; ++p) {
                    newLastOctantGlobalIdx += partition[p];
                }
                --newLastOctantGlobalIdx;
 
                newFirstOctantGlobalIdx = newLastOctantGlobalIdx - newSizeOctants + 1;
            } else {
                newLastOctantGlobalIdx  = 0;
                newFirstOctantGlobalIdx = 1;
            }
 
            // Partition internal octants
            if (m_serial) {
                m_lastOp = OP_LOADBALANCE_FIRST;
 
                if (newFirstOctantGlobalIdx != 0) {
                    for (uint32_t i = 0; i < newSizeOctants; ++i){
                        m_octree.m_octants[i] = m_octree.m_octants[newFirstOctantGlobalIdx + i];
                        if (userData) {
                            userData->move(newFirstOctantGlobalIdx + i, i);
                        }
                    }
                }
 
                m_octree.m_octants.resize(newSizeOctants);
                m_octree.m_octants.shrink_to_fit();
 
                if (userData) {
                    userData->resize(newSizeOctants);
                    userData->shrink();
                }
            } else {
                // Compute information about the current partitioning
                uint64_t lastOctantGlobalIdx  = m_partitionRangeGlobalIdx[m_rank];
                uint64_t firstOctantGlobalIdx = lastOctantGlobalIdx - m_octree.m_octants.size() + 1;
 
                // Initialize communications
                DataCommunicator lbCommunicator(m_comm);
 
                if (!userData || userData->fixedSize()) {
                    for (const auto &entry : recvRanges) {
                        int rank = entry.first;
                        uint32_t beginRecvIdx = entry.second[0];
                        uint32_t endRecvIdx   = entry.second[1];
 
                        uint32_t nOctantsToReceive = endRecvIdx - beginRecvIdx;
                        std::size_t buffSize = nOctantsToReceive * Octant::getBinarySize();
                        if (userData) {
                            buffSize += nOctantsToReceive * userData->fixedSize();
                        }
 
                        lbCommunicator.setRecv(rank, buffSize);
                        lbCommunicator.startRecv(rank);
                    }
                }
 
                for (const auto &entry : sendRanges) {
                    int rank = entry.first;
                    uint32_t beginSendIdx = entry.second[0];
                    uint32_t endSendIdx   = entry.second[1];
 
                    uint32_t nOctantsToSend = endSendIdx - beginSendIdx;
                    std::size_t buffSize = nOctantsToSend * Octant::getBinarySize();
                    if (userData) {
                        if (userData->fixedSize()) {
                            buffSize += nOctantsToSend * userData->fixedSize();
                        }  else {
                            for (uint32_t i = beginSendIdx; i < endSendIdx; ++i) {
                                buffSize += userData->size(i);
                            }
                        }
                    }
                    lbCommunicator.setSend(rank, buffSize);
 
                    SendBuffer &sendBuffer = lbCommunicator.getSendBuffer(rank);
                    for (uint32_t i = beginSendIdx; i < endSendIdx; ++i) {
                        sendBuffer << m_octree.m_octants[i];
                        userData->gather(sendBuffer, i);
                    }
 
                    lbCommunicator.startSend(rank);
                }
 
                if (userData && !userData->fixedSize()) {
                    lbCommunicator.discoverRecvs();
                    lbCommunicator.startAllRecvs();
                }
 
                // Move resident octants into place
                if (newSizeOctants > m_octree.m_octants.size()) {
                    m_octree.m_octants.reserve(newSizeOctants);
                    m_octree.m_octants.resize(newSizeOctants);
 
                    if (userData) {
                        userData->resize(newSizeOctants);
                    }
                }
 
                bool hasResidentOctants;
                if (newSizeOctants > 0) {
                    hasResidentOctants = true;
                    if (newFirstOctantGlobalIdx > lastOctantGlobalIdx) {
                        hasResidentOctants = false;
                    } else if (newLastOctantGlobalIdx < firstOctantGlobalIdx) {
                        hasResidentOctants = false;
                    }
                } else {
                    hasResidentOctants = false;
                }
 
                if (hasResidentOctants) {
                    uint64_t firstResidentGlobalOctantIdx = std::max(firstOctantGlobalIdx, newFirstOctantGlobalIdx);
                    uint64_t lastResidentGlobalOctantIdx  = std::min(lastOctantGlobalIdx, newLastOctantGlobalIdx);
 
                    uint32_t nofResidents              = static_cast<uint32_t>(lastResidentGlobalOctantIdx - firstResidentGlobalOctantIdx + 1);
                    uint32_t newFirstResidentOffsetIdx = static_cast<uint32_t>(firstResidentGlobalOctantIdx - newFirstOctantGlobalIdx);
                    uint32_t firstResidentOffsetIdx    = static_cast<uint32_t>(firstResidentGlobalOctantIdx - firstOctantGlobalIdx);
 
                    // If residents are moved closer to the head, we need to move
                    // them from the first to the last. Otherwise, if resident are
                    // moved closer to the tail, we need to move them in reversed
                    // order, i.e., from the last to the first (otherwise octants
                    // will be overwritten during the relocaiton).
                    if (newFirstResidentOffsetIdx != firstResidentOffsetIdx) {
                        for (uint32_t i = 0; i < nofResidents; ++i) {
                            int residentIdx;
                            if (newFirstResidentOffsetIdx < firstResidentOffsetIdx) {
                                residentIdx = i;
                            } else {
                                residentIdx = nofResidents - i - 1;
                            }
 
                            m_octree.m_octants[newFirstResidentOffsetIdx + residentIdx] = m_octree.m_octants[firstResidentOffsetIdx + residentIdx];
                            if (userData) {
                                userData->move(firstResidentOffsetIdx + residentIdx, newFirstResidentOffsetIdx + residentIdx);
                            }
                        }
                    }
                }
 
                if (newSizeOctants < m_octree.m_octants.size()) {
                    m_octree.m_octants.resize(newSizeOctants);
                    m_octree.m_octants.shrink_to_fit();
 
                    if (userData) {
                        userData->resize(newSizeOctants);
                        userData->shrink();
                    }
                }
 
                // Read buffers and build new octants
                int nCompletedRecvs = 0;
                while (nCompletedRecvs < lbCommunicator.getRecvCount()) {
                    int senderRank = lbCommunicator.waitAnyRecv();
                    RecvBuffer &recvBuffer = lbCommunicator.getRecvBuffer(senderRank);
 
                    const std::array<uint32_t, 2> &recvRange = recvRanges.at(senderRank);
                    uint32_t beginRecvIdx = recvRange[0];
                    uint32_t endRecvIdx   = recvRange[1];
                    assert(userData || ((endRecvIdx - beginRecvIdx) == (recvBuffer.getSize() / (Octant::getBinarySize()))));
                    assert(!userData || !userData->fixedSize() || ((endRecvIdx - beginRecvIdx) == (recvBuffer.getSize() / (Octant::getBinarySize() + userData->fixedSize()))));
 
                    for (uint32_t i = beginRecvIdx; i < endRecvIdx; ++i) {
                        recvBuffer >> m_octree.m_octants[i];
                        if (userData) {
                            userData->scatter(recvBuffer, i);
                        }
                    }
 
                    ++nCompletedRecvs;
                }
 
                lbCommunicator.waitAllSends();
            }
 
            // Update load balance information
            updateLoadBalance();
 
            // Update ghosts
            m_octree.m_ghosts.clear();
 
            computeGhostHalo();
 
            if (userData) {
                userData->resizeGhost(m_octree.getNumGhosts());
            }
        }
#endif
 
        // =============================================================================== //
 
    };
 
}
 
#endif /* __BITPIT_PARA_TREE_HPP__ */