mesh.h

//
// Created by milinda on 9/2/16.

//

#ifndef SFCSORTBENCH_MESH_H
#define SFCSORTBENCH_MESH_H

#include <vector>
#include "mpi.h"
#include "octUtils.h"
#include "treenode2vtk.h"
#include "TreeNode.h"
#include "sfcSort.h"
#include "sfcSearch.h"
#include "testUtils.h"

#include <memory>
#include <fdCoefficient.h>
#include "refel.h"
#include "block.h"
#include "stencil.h"
#include "key.h"
#include "skey.h"
#include "node.h"
#include "dendro.h"
#include "asyncExchangeContex.h"

#include "wavelet.h"
#include "dendroProfileParams.h" // only need to profile unzip_asyn for bssn. remove this header file later.

extern double t_e2e; // e2e map generation time
extern double t_e2n; // e2n map generation time
extern double t_sm;  // sm map generation time
extern double t_blk;

extern double t_e2e_g[3];
extern double t_e2n_g[3];
extern double t_sm_g[3];
extern double t_blk_g[3];

#include "dendroProfileParams.h"

enum KeyType
{
    NONE,
    SPLITTER,
    UP,
    DOWN,
    FRONT,
    BACK,
    LEFT,
    RIGHT
};

enum NeighbourLevel
{
    COARSE,
    SAME,
    REFINE
};

//#define DEBUG_MESH_GENERATION

#define KS_MAX 31 // maximum number of points in any direction that wavelet DA supports.
#define KEY_DIR_OFFSET 7
#define CHAINED_GHOST_OFFSET 5u

#define OCT_NO_CHANGE 0u
#define OCT_SPLIT 1u
#define OCT_COARSE 2u

namespace ot
{

enum SM_TYPE
{
    FDM = 0, // Finite Difference Method
    FEM_CG,  // Continous Galerkin  methods
    FEM_DG,   // Discontinous Galerkin methods
    E2E_ONLY // only builds the e2e maps. 
};

enum EType
{
    INDEPENDENT, // all the elemental nodal values,  are local.
    W_DEPENDENT, // element is writable but has ghost nodal values as well. (note that the writable would be Independent U W_dependent)
    UNKWON
};

namespace WaveletDA
{
enum LoopType
{
    ALL,
    INDEPENDENT,
    DEPENDENT
};
}

} // namespace ot

namespace ot
{
struct NodeTuple
{
private:
    unsigned int i;
    unsigned int j;
    unsigned int k;
    unsigned int e;

public:
    NodeTuple(unsigned int xi, unsigned int yj, unsigned int zk, unsigned int owner)
    {
        i = xi;
        j = yj;
        k = zk;
        e = owner;
    }
    inline unsigned int getX() const { return i; }
    inline unsigned int getY() const { return j; }
    inline unsigned int getZ() const { return k; }
    inline unsigned int getLevel() const { return e; }
};

} // namespace ot

namespace ot
{

class Mesh
{

private:
    std::vector<unsigned int> m_uiE2EMapping;
    std::vector<unsigned int> m_uiE2NMapping_CG;
    std::vector<unsigned int> m_uiE2NMapping_DG;
    std::vector<unsigned int> m_uiCG2DG;
    std::vector<unsigned int> m_uiDG2CG;

    std::vector<ot::TreeNode> m_uiLocalSplitterElements; // used to spit the keys to the correct nodes.

    ot::TreeNode *m_uiSplitterNodes;

    // Pre  and Post ghost octants.

    std::vector<ot::TreeNode> m_uiPreGhostOctants;
    std::vector<ot::TreeNode> m_uiPostGhostOctants;

    // Keys Related attributes
    std::vector<Key> m_uiKeys;

    std::vector<Key> m_uiGhostKeys;

    std::vector<Key> m_uiKeysDiag;

    std::vector<ot::TreeNode> m_uiEmbeddedOctree;

    std::vector<ot::TreeNode> m_uiGhostOctants;

    std::vector<unsigned int> m_uiGhostElementRound1Index;

    std::vector<ot::TreeNode> m_uiAllElements;
    std::vector<ot::TreeNode> m_uiAllLocalNode;

    std::vector<ot::Block> m_uiLocalBlockList;

    unsigned int m_uiDmin;
    unsigned int m_uiDmax;

    unsigned int m_uiElementPreGhostBegin = 0; // Not mandatory to store this.
    unsigned int m_uiElementPreGhostEnd;
    unsigned int m_uiElementLocalBegin;
    unsigned int m_uiElementLocalEnd;
    unsigned int m_uiElementPostGhostBegin;
    unsigned int m_uiElementPostGhostEnd;
    unsigned int m_uiFElementPreGhostBegin;
    unsigned int m_uiFElementPreGhostEnd;
    unsigned int m_uiFElementLocalBegin;
    unsigned int m_uiFElementLocalEnd;
    unsigned int m_uiFElementPostGhostBegin;
    unsigned int m_uiFElementPostGhostEnd;

    unsigned int m_uiMeshDomain_min;
    unsigned int m_uiMeshDomain_max;

    unsigned int m_uiNumLocalElements;
    unsigned int m_uiNumPreGhostElements;
    unsigned int m_uiNumPostGhostElements;

    /*** Total number of actual elements in the mesh (preG+ local + postG ) */
    unsigned int m_uiNumTotalElements;

    unsigned int m_uiNumActualNodes;

    DendroIntL m_uiUnZippedVecSz;

    unsigned int m_uiNumFakeNodes;

    unsigned int m_uiNodePreGhostBegin;
    unsigned int m_uiNodePreGhostEnd;
    unsigned int m_uiNodeLocalBegin;
    unsigned int m_uiNodeLocalEnd;
    unsigned int m_uiNodePostGhostBegin;
    unsigned int m_uiNodePostGhostEnd;

    // Stuff for communication and synchronization ...
    //-------------------------------------------------------------------------

    int m_uiActiveRank;
    int m_uiActiveNpes;

    MPI_Comm m_uiCommActive;

    MPI_Comm m_uiCommGlobal;

    int m_uiGlobalRank;
    int m_uiGlobalNpes;

    bool m_uiIsActive;

    unsigned int *m_uiSendKeyCount;
    unsigned int *m_uiRecvKeyCount;
    unsigned int *m_uiSendKeyOffset;
    unsigned int *m_uiRecvKeyOffset;

    unsigned int *m_uiSendKeyDiagCount;
    unsigned int *m_uiRecvKeyDiagCount;
    unsigned int *m_uiSendKeyDiagOffset;
    unsigned int *m_uiRecvKeyDiagOffset;

    unsigned int *m_uiSendOctCountRound1;
    unsigned int *m_uiRecvOctCountRound1;
    unsigned int *m_uiSendOctOffsetRound1;
    unsigned int *m_uiRecvOctOffsetRound1;

    unsigned int *m_uiSendOctCountRound1Diag;
    unsigned int *m_uiRecvOctCountRound1Diag;
    unsigned int *m_uiSendOctOffsetRound1Diag;
    unsigned int *m_uiRecvOctOffsetRound1Diag;

    unsigned int *m_uiSendOctCountRound2;
    unsigned int *m_uiRecvOctCountRound2;
    unsigned int *m_uiSendOctOffsetRound2;
    unsigned int *m_uiRecvOctOffsetRound2;

    std::vector<unsigned int> m_uiSendNodeCount;
    std::vector<unsigned int> m_uiRecvNodeCount;
    std::vector<unsigned int> m_uiSendNodeOffset;
    std::vector<unsigned int> m_uiRecvNodeOffset;

    std::vector<unsigned int> m_uiSendProcList;

    std::vector<unsigned int> m_uiRecvProcList;

    std::vector<unsigned int> m_uiGhostElementIDsToBeSent;
    std::vector<unsigned int> m_uiGhostElementIDsToBeRecv;

    std::vector<unsigned int> m_uiPreGhostHangingNodeCGID;

    std::vector<unsigned int> m_uiPostGhostHangingNodeCGID;

    std::vector<ot::TreeNode> m_uiSendBufferElement;

    std::vector<double> m_uiSendBufferNodes;
    std::vector<double> m_uiRecvBufferNodes;

    std::vector<unsigned int> m_uiScatterMapElementRound1;

    std::vector<unsigned int> m_uiScatterMapActualNodeSend;

    std::vector<unsigned int> m_uiScatterMapActualNodeRecv;

    // variables to manage loop access over elements.
    unsigned int m_uiEL_i;

    unsigned int m_uiElementOrder;

    unsigned int m_uiNpE;

    unsigned int m_uiStensilSz;

    unsigned int m_uiNumDirections;

    RefElement m_uiRefEl;

    std::vector<unsigned int> m_uiFEMGhostLev1IDs;

    //===== maps needed for DG computations.

    std::vector<unsigned int> m_uiF2EMap;

    bool m_uiIsBlockSetup;

    SM_TYPE m_uiScatterMapType;

    bool m_uiIsF2ESetup;

    std::vector<AsyncExchangeContex> m_uiMPIContexts;

    unsigned int m_uiCommTag=0;

    std::vector<bool> m_uiIsNodalMapValid;

    // --
    //Note : These are special data stored to search the 3rd point in Dendro-GR unzip with 4th order.
    //

    std::vector<ot::Key> m_uiUnzip_3pt_keys;

    std::vector<ot::Key> m_uiUnzip_3pt_ele;

    std::vector<ot::Key> m_uiUnzip_3pt_recv_keys;

    std::vector<unsigned int> m_uiSendCountRePt;

    std::vector<unsigned int> m_uiSendOffsetRePt;

    std::vector<unsigned int> m_uiRecvCountRePt;

    std::vector<unsigned int> m_uiRecvOffsetRePt;

    std::vector<unsigned int> m_uiReqSendProcList;
    
    std::vector<unsigned int> m_uiReqRecvProcList;

    std::vector<unsigned int> m_uiSendNodeReqPtSM;


private:
    void buildFEM_E2N();

    void generateSearchKeys();

    void generateGhostElementSearchKeys();

    void generateBdyElementDiagonalSearchKeys();

    void buildE2EMap(std::vector<ot::TreeNode> &in, MPI_Comm comm);

    void computeElementOwnerRanks(std::vector<unsigned int> &elementOwner);

    void buildE2EMap(std::vector<ot::TreeNode> &in);

    void buildE2NMap();

    void shrinkE2EAndE2N();

    void generateFakeElementsAndUpdateE2EAndE2N();

    inline bool computeOveralppingNodes(const ot::TreeNode &parent, const ot::TreeNode &child, int *idx, int *idy, int *idz);

    void computeNodeScatterMaps(MPI_Comm comm);

    void computeNodalScatterMap(MPI_Comm comm);

    void computeNodalScatterMap1(MPI_Comm comm);

    void computeNodalScatterMap2(MPI_Comm comm);

    void computeNodalScatterMap3(MPI_Comm comm);

    void computeNodalScatterMap4(MPI_Comm comm);

    /*
         * NOTE: These private functions does not check for m_uiIsActive
         * sicne they are called in buildE2E and buildE2N mappings.
         *
         * */

    inline void OCT_DIR_LEFT_INTERNAL_EDGE_MAP(unsigned int child, unsigned int parent, bool parentChildLevEqual, std::vector<unsigned int> &edgeChildIndex, std::vector<unsigned int> &edgeOwnerIndex);
    inline void OCT_DIR_RIGHT_INTERNAL_EDGE_MAP(unsigned int child, unsigned int parent, bool parentChildLevEqual, std::vector<unsigned int> &edgeChildIndex, std::vector<unsigned int> &edgeOwnerIndex);
    inline void OCT_DIR_DOWN_INTERNAL_EDGE_MAP(unsigned int child, unsigned int parent, bool parentChildLevEqual, std::vector<unsigned int> &edgeChildIndex, std::vector<unsigned int> &edgeOwnerIndex);
    inline void OCT_DIR_UP_INTERNAL_EDGE_MAP(unsigned int child, unsigned int parent, bool parentChildLevEqual, std::vector<unsigned int> &edgeChildIndex, std::vector<unsigned int> &edgeOwnerIndex);
    inline void OCT_DIR_FRONT_INTERNAL_EDGE_MAP(unsigned int child, unsigned int parent, bool parentChildLevEqual, std::vector<unsigned int> &edgeChildIndex, std::vector<unsigned int> &edgeOwnerIndex);
    inline void OCT_DIR_BACK_INTERNAL_EDGE_MAP(unsigned int child, unsigned int parent, bool parentChildLevEqual, std::vector<unsigned int> &edgeChildIndex, std::vector<unsigned int> &edgeOwnerIndex);

    void SET_FACE_TO_ELEMENT_MAP(unsigned int ele, unsigned int dir, unsigned int dirOp, unsigned int dir1, unsigned int dir2);

    inline void CORNER_NODE_MAP(unsigned int child);

    inline void OCT_DIR_DIAGONAL_E2E(unsigned int elementID, unsigned int face1, unsigned int face2, unsigned int &lookUp) const;

    template <typename pKey, typename pNode>
    void searchKeys(std::vector<pKey> &pKeys, std::vector<pNode> &pNodes);

    inline unsigned int getDIROfANode(unsigned int ii_x, unsigned int jj_y, unsigned int kk_z);

    inline void directionToIJK(unsigned int direction, std::vector<unsigned int> &ii_x, std::vector<unsigned int> &jj_y, std::vector<unsigned int> &kk_z);

    template <typename T>
    inline void interpUpWind(const double *upWind, const unsigned int element, const unsigned int lookup, T *vecLookUp, const unsigned int cnum, const T *parentInterpIn, T *parentInterpOut, const unsigned int padDir, const unsigned int padWidth, const T *zippedVec, T *out);

    template <typename T>
    inline void interpDownWind(const double *downWind, const unsigned int element, const unsigned int lookup, T *vecLookUp, const unsigned int cnum, const T *parentInterpIn, T *parentInterpOut, const unsigned int padDir, const unsigned int padWidth, const T *zippedVec, T *out);

    template <typename T>
    void blockDiagonalUnZip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void blockVertexUnZip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_DOWN_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_UP_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_DOWN_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_UP_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_UP_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_UP_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_DOWN_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_DOWN_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_DOWN_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_DOWN_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_UP_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_UP_BACK_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_DOWN_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_DOWN_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_LEFT_UP_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    template <typename T>
    void OCT_DIR_RIGHT_UP_FRONT_Unzip(const ot::Block &blk, const T *zippedVec, T *unzippedVec);

    //---Note: These functions are specifically written for find missing 3rd block unzip points for the GR application

    void computeSMSpecialPts();

    template <typename T>
    void readSpecialPtsBegin(const T *in);

    template <typename T>
    void readSpecialPtsEnd(const T *in, T* out);

    // --- 3rd point exchange function end.

public:
    Mesh(std::vector<ot::TreeNode> &in, unsigned int k_s, unsigned int pOrder, unsigned int activeNpes, MPI_Comm comm, bool pBlockSetup = true, SM_TYPE smType = SM_TYPE::FDM, unsigned int grainSz = DENDRO_DEFAULT_GRAIN_SZ, double ld_tol = DENDRO_DEFAULT_LB_TOL, unsigned int sf_k = DENDRO_DEFAULT_SF_K);

    Mesh(std::vector<ot::TreeNode> &in, unsigned int k_s, unsigned int pOrder, MPI_Comm comm, bool pBlockSetup = true, SM_TYPE smType = SM_TYPE::FDM, unsigned int grainSz = DENDRO_DEFAULT_GRAIN_SZ, double ld_tol = DENDRO_DEFAULT_LB_TOL, unsigned int sf_k = DENDRO_DEFAULT_SF_K);

    ~Mesh();

    void performBlocksSetup();

    void buildF2EMap();

    /*** @brief: perform block dependancy flags */
    void flagBlocks();

    inline bool isBlockSetep() { return m_uiIsBlockSetup; }

    inline SM_TYPE getScatterMapType() { return m_uiScatterMapType; }

    // Setters and getters.
    inline unsigned int getNumLocalMeshElements() const { return m_uiNumLocalElements; }
    inline unsigned int getNumPreGhostElements() const { return m_uiNumPreGhostElements; }
    inline unsigned int getNumPostGhostElements() const { return m_uiNumPostGhostElements; }
    // get nodal information.

    inline unsigned int getNumLocalMeshNodes() const { return (m_uiNodeLocalEnd - m_uiNodeLocalBegin); }
    inline unsigned int getNumPreMeshNodes() const { return (m_uiNodePreGhostEnd - m_uiNodePreGhostBegin); }
    inline unsigned int getNumPostMeshNodes() const { return (m_uiNodePostGhostEnd - m_uiNodePostGhostBegin); }
    inline unsigned int getElementPreGhostBegin() const { return m_uiElementPreGhostBegin; }
    inline unsigned int getElementPreGhostEnd() const { return m_uiElementPreGhostEnd; }
    inline unsigned int getElementLocalBegin() const { return m_uiElementLocalBegin; }
    inline unsigned int getElementLocalEnd() const { return m_uiElementLocalEnd; }
    inline unsigned int getElementPostGhostBegin() const { return m_uiElementPostGhostBegin; }
    inline unsigned int getElementPostGhostEnd() const { return m_uiElementPostGhostEnd; }

    inline unsigned int getNodePreGhostBegin() const { return m_uiNodePreGhostBegin; }
    inline unsigned int getNodePreGhostEnd() const { return m_uiNodePreGhostEnd; }
    inline unsigned int getNodeLocalBegin() const { return m_uiNodeLocalBegin; }
    inline unsigned int getNodeLocalEnd() const { return m_uiNodeLocalEnd; }
    inline unsigned int getNodePostGhostBegin() const { return m_uiNodePostGhostBegin; }
    inline unsigned int getNodePostGhostEnd() const { return m_uiNodePostGhostEnd; }

    inline unsigned int getDegOfFreedom() const { return m_uiNumActualNodes; }

    inline unsigned int getDegOfFreedomUnZip() const { return m_uiUnZippedVecSz; }

    inline const std::vector<ot::TreeNode> &getAllElements() const { return m_uiAllElements; }

    inline const std::vector<unsigned int> &getLevel1GhostElementIndices() const { return m_uiGhostElementRound1Index; }

    inline const std::vector<ot::TreeNode> &getSplitterElements() const { return m_uiLocalSplitterElements; }

    inline const std::vector<ot::TreeNode> &getAllLocalNodes() const { return m_uiAllLocalNode; }

    inline const std::vector<unsigned int> &getE2EMapping() const { return m_uiE2EMapping; }

    inline const std::vector<unsigned int> &getE2NMapping() const { return m_uiE2NMapping_CG; }

    inline const std::vector<unsigned int> &getDG2CGMap() const { return m_uiDG2CG; }

    inline const std::vector<unsigned int> &getCG2DGMap() const { return m_uiCG2DG; }

    inline const std::vector<unsigned int> &getE2NMapping_DG() const { return m_uiE2NMapping_DG; }
    inline const std::vector<ot::Block> &getLocalBlockList() const { return m_uiLocalBlockList; }

    inline unsigned int getNumDirections() const { return m_uiNumDirections; }

    inline unsigned int getNumNodesPerElement() const { return m_uiNpE; }

    inline unsigned int getElementOrder() const { return m_uiElementOrder; }

    inline MPI_Comm getMPICommunicator() const { return m_uiCommActive; }

    inline MPI_Comm getMPIGlobalCommunicator() const { return m_uiCommGlobal; }

    inline unsigned int getMPIRankGlobal() const { return m_uiGlobalRank; }

    inline unsigned int getMPICommSizeGlobal() const { return m_uiGlobalNpes; }

    inline unsigned int getMPIRank() const { return m_uiActiveRank; }

    inline unsigned int getMPICommSize() const { return m_uiActiveNpes; }

    inline const RefElement *getReferenceElement() const { return &m_uiRefEl; }

    inline const unsigned int getSendProcListSize() const { return m_uiSendProcList.size(); }

    inline const unsigned int getRecvProcListSize() const { return m_uiRecvProcList.size(); }

    inline const std::vector<unsigned int> &getNodalSendCounts() const { return m_uiSendNodeCount; }

    inline const std::vector<unsigned int> &getNodalSendOffsets() const { return m_uiSendNodeOffset; }

    inline const std::vector<unsigned int> &getNodalRecvCounts() const { return m_uiRecvNodeCount; }

    inline const std::vector<unsigned int> &getNodalRecvOffsets() const { return m_uiRecvNodeOffset; }

    inline const std::vector<unsigned int> &getSendProcList() const { return m_uiSendProcList; }

    inline const std::vector<unsigned int> &getRecvProcList() const { return m_uiRecvProcList; }

    inline const std::vector<unsigned int> &getSendNodeSM() const { return m_uiScatterMapActualNodeSend; }

    inline const std::vector<unsigned int> &getRecvNodeSM() const { return m_uiScatterMapActualNodeRecv; }

    inline void dg2eijk(unsigned int dg_index, unsigned int &e, unsigned int &i, unsigned int &j, unsigned int &k) const
    {
        e = dg_index / m_uiNpE;
        k = 0;
        j = 0;
        i = 0;

        //std::cout<<"e: "<<e<<std::endl;
        if (dg_index > e * m_uiNpE)
            k = (dg_index - e * m_uiNpE) / ((m_uiElementOrder + 1) * (m_uiElementOrder + 1));
        //std::cout<<"k: "<<k<<std::endl;
        if ((dg_index + k * ((m_uiElementOrder + 1) * (m_uiElementOrder + 1))) > (e * m_uiNpE))
            j = (dg_index - e * m_uiNpE - k * ((m_uiElementOrder + 1) * (m_uiElementOrder + 1))) / (m_uiElementOrder + 1);
        //std::cout<<"j: "<<j<<std::endl;
        if ((dg_index + k * ((m_uiElementOrder + 1) * (m_uiElementOrder + 1)) + j * (m_uiElementOrder + 1)) > (e * m_uiNpE))
            i = (dg_index - e * m_uiNpE - k * ((m_uiElementOrder + 1) * (m_uiElementOrder + 1)) - j * (m_uiElementOrder + 1));
        //std::cout<<"i: "<<i<<std::endl;
    }

    inline unsigned int getMortonchildNum(unsigned int eleID) const
    {
        return m_uiAllElements[eleID].getMortonIndex();
    }

    inline bool isActive() const { return m_uiIsActive; }

    inline void waitAll() const { MPI_Barrier(m_uiCommGlobal); }

    void setOctreeRefineFlags(unsigned int *flags, unsigned int sz);

    void getElementalFaceNeighbors(const unsigned int eID, const unsigned int dir, unsigned int *lookup) const;
    void getElementalEdgeNeighbors(const unsigned int eID, const unsigned int dir, unsigned int *lookup) const;
    void getElementalVertexNeighbors(const unsigned int eID, const unsigned int dir, unsigned int *lookup) const;

    void getElementQMat(unsigned int currentId, double *&qMat, bool isAllocated = true) const;


    template<typename T>
    void getUnzipElementalNodalValues(const T* uzipVec, unsigned int blkID, unsigned int ele, T*out, bool isPadded=true) const;

    // Wavelet Init functions
    template <ot::WaveletDA::LoopType type>
    void init();

    template <ot::WaveletDA::LoopType type>
    bool nextAvailable();

    template <ot::WaveletDA::LoopType type>
    void next();

    inline const ot::TreeNode &currentOctant();

    inline unsigned int currentIndex();

    inline void currentElementNeighbourIndexList(unsigned int *neighList);

    inline void currentElementNodeList(unsigned int *nodeList);

    inline void currentElementNodeList_DG(unsigned int *nodeList);

    // functions to access, faces and edges of a given element.

    inline void faceNodesIndex(unsigned int elementID, unsigned int face, std::vector<unsigned int> &index,
                               bool isInternal) const;

    inline void edgeNodeIndex(unsigned int elementID, unsigned int face1, unsigned int face2, std::vector<unsigned int> &index, bool isInternal) const;

    inline void cornerNodeIndex(unsigned int elementID, unsigned int mortonIndex, unsigned int &index) const;

    inline void elementNodeIndex(unsigned int elementID, std::vector<unsigned int> &index, bool isInternal) const;

    bool isEdgeHanging(unsigned int elementId, unsigned int edgeId, unsigned int &cnum) const;

    bool isFaceHanging(unsigned int elementId, unsigned int faceId, unsigned int &cnum) const;

    bool isNodeHanging(unsigned int eleID, unsigned int ix, unsigned int jy, unsigned int kz) const;

    inline bool isNodeLocal(unsigned int eleID, unsigned int ix, unsigned int jy, unsigned int kz) const
    {
        return ((m_uiE2NMapping_CG[eleID * m_uiNpE + kz * (m_uiElementOrder + 1) * (m_uiElementOrder + 1) + jy * (m_uiElementOrder + 1) + ix] >= m_uiNodeLocalBegin) && (m_uiE2NMapping_CG[eleID * m_uiNpE + kz * (m_uiElementOrder + 1) * (m_uiElementOrder + 1) + jy * (m_uiElementOrder + 1) + ix] < m_uiNodeLocalEnd));
    };

    inline bool isBoundaryOctant(unsigned int ele) const 
    {
        assert(ele < m_uiAllElements.size() );
        return (m_uiAllElements[ele].minX() ==0 || m_uiAllElements[ele].minY() ==0 || m_uiAllElements[ele].minZ()==0 || m_uiAllElements[ele].maxX() == 1u<<(m_uiMaxDepth) || m_uiAllElements[ele].maxY() == 1u<<(m_uiMaxDepth) || m_uiAllElements[ele].maxZ() == 1u<<(m_uiMaxDepth) );
    }

    inline DendroIntL getGhostExcgTotalSendNodeCount() const
    {
        if (m_uiGlobalNpes == 1)
            return 0;

        if (m_uiIsActive)
            return (m_uiSendNodeOffset[m_uiActiveNpes - 1] + m_uiSendNodeCount[m_uiActiveNpes - 1]);
        else
            return 0;
    }

    inline DendroIntL getGhostExcgTotalRecvNodeCount() const
    {
        if (m_uiGlobalNpes == 1)
            return 0;

        if (m_uiIsActive)
            return (m_uiRecvNodeOffset[m_uiActiveNpes - 1] + m_uiRecvNodeCount[m_uiActiveNpes - 1]);
        else
            return 0;
    }

    inline const ot::TreeNode *getNodalSplitterNodes() const { return m_uiSplitterNodes; }

    // Methods needed for PDE & ODE and other solvers.
    template <typename T>
    T *createVector() const;

    template <typename T>
    void createVector(std::vector<T> &vec) const;

    template <typename T>
    T *createVector(const T initValue) const;

    template <typename T>
    T *createVector(std::function<T(T, T, T)> func) const;

    template <typename T>
    void createVector(std::vector<T> &vec, const T initValue) const;

    template <typename T>
    void createVector(std::vector<T> &vec, std::function<T(T, T, T)> func) const;

    template <typename T>
    void createUnZippedVector(std::vector<T> &uvec) const;

    template <typename T>
    T *createUnZippedVector() const;

    template <typename T>
    void createUnZippedVector(std::vector<T> &uvec, const T initValue) const;

    template <typename T>
    T *createUnZippedVector(const T initValue) const;

    inline void parent2ChildInterpolation(const double *in, double *out, unsigned int cnum, unsigned int dim = 3) const;

    inline void child2ParentInterpolation(const double *in, double *out, unsigned int cnum, unsigned int dim = 3) const;

    inline void child2ParentInjection(double *in, double *out, bool *isHanging) const;

    template <typename T>
    void unzip(const T *zippedVec, T *unzippedVec);

    template <typename T>
    void unzip_async(T *zippedVec, T *unzippedVec, MPI_Request *send_reqs, MPI_Request *recv_reqs, MPI_Status *send_sts, MPI_Status *recv_sts);

    template <typename T>
    void zip(const T *unzippedVec, T *zippedVec);

    template <typename T, unsigned int length, unsigned int offsetCentered, unsigned int offsetBackward, unsigned int offsetForward>
    void applyStencil(const std::vector<T> &in, std::vector<T> &out,
                      const Stencil<T, length, offsetCentered> &centered,
                      const Stencil<T, length, offsetBackward> &backward,
                      const Stencil<T, length, offsetForward> &forward);

    template <typename T>
    void performGhostExchange(std::vector<T> &vec);

    template <typename T>
    void performGhostExchange(T *vec);

    template <typename T>
    void ghostExchangeStart(T *vec, T *sendNodeBuffer, T *recvNodeBuffer, MPI_Request *send_reqs, MPI_Request *recv_reqs);

    template <typename T>
    void ghostExchangeRecvSync(T *vec, T *recvNodeBuffer, MPI_Request *recv_reqs, MPI_Status *recv_sts);

    inline void ghostExchangeSendSync(MPI_Request *send_reqs, MPI_Status *send_sts)
    {
        MPI_Waitall(m_uiSendProcList.size(), send_reqs, send_sts);
    }

    template <typename T>
    void vectorToVTK(const std::vector<T> &vec, char *fprefix, double pTime = 0.0, unsigned int nCycle = 0) const;

    template <typename T>
    bool isReMesh(const T **vec, const unsigned int *varIds, const unsigned int numVars, double tol, double amr_coarse_fac = DENDRO_AMR_COARSEN_FAC);

    template <typename T>
    bool isReMeshUnzip(const T **unzippedVec, const unsigned int *varIds, const unsigned int numVars, std::function<double(double, double, double)> wavelet_tol, double amr_coarse_fac = DENDRO_AMR_COARSEN_FAC, double coarsen_hx = DENDRO_REMESH_UNZIP_SCALE_FAC);

    ot::Mesh *ReMesh(unsigned int grainSz = DENDRO_DEFAULT_GRAIN_SZ, double ld_tol = DENDRO_DEFAULT_LB_TOL, unsigned int sfK = DENDRO_DEFAULT_SF_K,unsigned int (*getWeight)(const ot::TreeNode *)=NULL);

    template <typename T>
    void interGridTransfer(std::vector<T> &vec, const ot::Mesh *pMesh);

    template <typename T>
    void interGridTransfer(T *&vec, const ot::Mesh *pMesh);

    template<typename T>
    void interGridTransferUnzip(T*& unzip, T*& vec, const ot::Mesh *pMesh);

    template <typename T>
    void getElementNodalValues(const T *vec, T *nodalValues, unsigned int elementID) const;

    template <typename T>
    void computeElementalContribution(const T *in, T *out, unsigned int elementID) const;

    void getElementCoordinates(unsigned int eleID, double *coords) const;

    template <typename T>
    int getFaceNeighborValues(unsigned int eleID, const T *in, T *out, T *coords, unsigned int *neighID, unsigned int face, NeighbourLevel &level) const;

    EType getElementType(unsigned int eleID);

    int getBlkBdyParentNodeIndices(unsigned int blkId, unsigned int eleId, unsigned int dir, unsigned int* nid, unsigned int* child, unsigned int* fid, unsigned int* cid);

};

template <>
inline void Mesh::init<WaveletDA::LoopType ::ALL>()
{
    m_uiEL_i = m_uiElementPreGhostBegin;
}

template <>
inline void Mesh::init<WaveletDA::INDEPENDENT>()
{
    m_uiEL_i = m_uiElementLocalBegin;
}

template <>
inline void Mesh::init<WaveletDA::DEPENDENT>()
{
    m_uiEL_i = m_uiElementPreGhostBegin;
}

template <>
inline bool Mesh::nextAvailable<WaveletDA::ALL>()
{
    return (m_uiEL_i < m_uiElementPostGhostEnd);
}

template <>
inline bool Mesh::nextAvailable<WaveletDA::INDEPENDENT>()
{
    return (m_uiEL_i < m_uiElementLocalEnd);
}
template <>

inline bool Mesh::nextAvailable<WaveletDA::DEPENDENT>()
{
    return (m_uiEL_i < m_uiElementPreGhostEnd) || ((m_uiEL_i > m_uiElementLocalEnd) && m_uiEL_i < m_uiElementPostGhostEnd);
}

template <>
inline void Mesh::next<WaveletDA::ALL>()
{
    // interpolation should come here.
    m_uiEL_i++;
}

template <>
inline void Mesh::next<WaveletDA::INDEPENDENT>()
{
    // interpolation should come here.
    m_uiEL_i++;
}

template <>
inline void Mesh::next<WaveletDA::DEPENDENT>()
{
    // interpolation should come here.
    m_uiEL_i++;
    if (m_uiEL_i == m_uiElementPreGhostEnd)
        m_uiEL_i = m_uiElementPostGhostBegin;
}

inline const ot::TreeNode &Mesh::currentOctant()
{
    return m_uiAllElements[m_uiEL_i];
}

inline unsigned int Mesh::currentIndex() { return m_uiEL_i; }

inline void Mesh::currentElementNeighbourIndexList(unsigned int *neighList)
{
    if (!m_uiIsActive)
        return;

    for (unsigned int k = 0; k < m_uiNumDirections; k++)
        neighList[k] = m_uiE2EMapping[m_uiEL_i * m_uiNumDirections + k];
}

inline void Mesh::currentElementNodeList(unsigned int *nodeList)
{
    if (!m_uiIsActive)
        return;

    for (unsigned int k = 0; k < m_uiNpE; k++)
    {
        nodeList[k] = m_uiE2NMapping_CG[m_uiEL_i * m_uiNpE + k];
    }
}

inline void Mesh::currentElementNodeList_DG(unsigned int *nodeList)
{
    if (!m_uiIsActive)
        return;

    for (unsigned int k = 0; k < m_uiNpE; k++)
    {
        nodeList[k] = m_uiE2NMapping_DG[m_uiEL_i * m_uiNpE + k];
    }
}

inline void Mesh::parent2ChildInterpolation(const double *in, double *out, unsigned int cnum, unsigned int dim) const
{

    dendro::timer::t_unzip_p2c.start();
    if (dim == 3)
        m_uiRefEl.I3D_Parent2Child(in, out, cnum);
    else if (dim == 2)
        m_uiRefEl.I2D_Parent2Child(in, out, cnum);
    else if (dim == 1)
        m_uiRefEl.I1D_Parent2Child(in, out, cnum);
    dendro::timer::t_unzip_p2c.stop();
}

inline void Mesh::child2ParentInterpolation(const double *in, double *out, unsigned int cnum, unsigned int dim) const
{

    if (dim == 3)
        m_uiRefEl.I3D_Child2Parent(in, out, cnum);
    else if (dim == 2)
        m_uiRefEl.I2D_Child2Parent(in, out, cnum);
    else if (dim == 1)
        m_uiRefEl.I1D_Child2Parent(in, out, cnum);
}

inline void Mesh::child2ParentInjection(double *in, double *out, bool *isHanging) const
{

    assert(m_uiElementOrder % 2 == 0 || m_uiElementOrder == 1);
    if (m_uiElementOrder == 1)
    { // special case
        for (unsigned int child = 0; child < NUM_CHILDREN; child++)
        {
            for (unsigned int k = 0; k < (m_uiElementOrder + 1); k++)
                for (unsigned int j = 0; j < (m_uiElementOrder + 1); j++)
                    for (unsigned int i = 0; i < (m_uiElementOrder + 1); i++)
                        out[k * (m_uiElementOrder + 1) * (m_uiElementOrder + 1) + j * (m_uiElementOrder + 1) + i] = in[child * m_uiNpE + k * (m_uiElementOrder + 1) + j * (m_uiElementOrder + 1) + i];
        }
    }
    else
    {

        for (unsigned int child = 0; child < NUM_CHILDREN; child++)
        {

            for (unsigned int k = 0; k < (m_uiElementOrder + 1); k++)
                for (unsigned int j = 0; j < (m_uiElementOrder + 1); j++)
                    for (unsigned int i = 0; i < (m_uiElementOrder + 1); i++)
                    {
                        if ((i % 2 == 0) && (j % 2 == 0) && (k % 2 == 0) && isHanging[k * (m_uiElementOrder + 1) * (m_uiElementOrder + 1) + j * (m_uiElementOrder + 1) + i])
                        {
                            out[((((child & 4u) >> 2u) * m_uiElementOrder + k) >> 1) * (m_uiElementOrder + 1) * (m_uiElementOrder + 1) + ((((child & 2u) >> 1u) * m_uiElementOrder + j) >> 1) * (m_uiElementOrder + 1) + (((((child & (1u))) * m_uiElementOrder + i) >> 1))] = in[child * m_uiNpE + k * (m_uiElementOrder + 1) * (m_uiElementOrder + 1) + j * (m_uiElementOrder + 1) + i];
                            //if(!m_uiActiveRank) std::cout<<"rank: "<<m_uiActiveRank<<" child: "<<child<<" i: "<<i<<" j: "<<j<<" k: "<<k<<" inserted to node : "<<((((child & 4u)>>2u)*m_uiElementOrder+k)>>1)*(m_uiElementOrder+1)*(m_uiElementOrder+1)+((((child & 2u)>>1u)*m_uiElementOrder+j)>>1)*(m_uiElementOrder+1)+(((((child & (1u)))*m_uiElementOrder+i)>>1))<<std::endl;
                        }
                    }
        }
    }
}

inline bool Mesh::computeOveralppingNodes(const ot::TreeNode &parent, const ot::TreeNode &child, int *idx, int *idy, int *idz)
{

    unsigned int Lp = 1u << (m_uiMaxDepth - parent.getLevel());
    unsigned int Lc = 1u << (m_uiMaxDepth - child.getLevel());
    //intilize the mapping to -1. -1 denotes that mapping is not defined for given k value.

    unsigned int dp, dc;
    dp = (m_uiElementOrder);
    dc = m_uiElementOrder;

    assert(Lp % dp == 0);
    assert(Lc % dc == 0);

    for (unsigned int k = 0; k < (m_uiElementOrder + 1); k++)
    {
        idx[k] = -1;
        idy[k] = -1;
        idz[k] = -1;
    }
    bool state = false;
    bool stateX = false;
    bool stateY = false;
    bool stateZ = false;
    if (parent == child)
    {
        for (unsigned int k = 0; k < (m_uiElementOrder + 1); k++)
        {
            idx[k] = k;
            idy[k] = k;
            idz[k] = k;
        }
        return true;
    }
    else if (parent.isAncestor(child))
    {

        /*if((((child.getX()-parent.getX())*m_uiElementOrder)%Lp) || (((child.getY()-parent.getY())*m_uiElementOrder)%Lp) || (((child.getZ()-parent.getZ())*m_uiElementOrder)%Lp)) return false;
            else*/
        {
            unsigned int index[3];
            for (unsigned int k = 0; k < (m_uiElementOrder + 1); k++)
            {

                index[0] = (m_uiElementOrder + 1);
                index[1] = (m_uiElementOrder + 1);
                index[2] = (m_uiElementOrder + 1);

                if (!(((child.getX() - parent.getX()) * dp * dc + k * Lc * dp) % (Lp * dc)))
                    index[0] = ((child.getX() - parent.getX()) * dp * dc + k * Lc * dp) / (Lp * dc); //((child.getX()-parent.getX())*m_uiElementOrder + k*Lc)/Lp;

                if (!(((child.getY() - parent.getY()) * dp * dc + k * Lc * dp) % (Lp * dc)))
                    index[1] = ((child.getY() - parent.getY()) * dp * dc + k * Lc * dp) / (Lp * dc); //((child.getY()-parent.getY())*m_uiElementOrder + k*Lc)/Lp;

                if (!(((child.getZ() - parent.getZ()) * dp * dc + k * Lc * dp) % (Lp * dc)))
                    index[2] = ((child.getZ() - parent.getZ()) * dp * dc + k * Lc * dp) / (Lp * dc); //((child.getZ()-parent.getZ())*m_uiElementOrder + k*Lc)/Lp;

                if (!stateX && index[0] < (m_uiElementOrder + 1))
                    stateX = true;

                if (!stateY && index[1] < (m_uiElementOrder + 1))
                    stateY = true;

                if (!stateZ && index[2] < (m_uiElementOrder + 1))
                    stateZ = true;

                if (index[0] < (m_uiElementOrder + 1))
                {
                    idx[k] = index[0];
                    assert((parent.getX() + idx[k] * Lp / dp) == (child.getX() + k * Lc / dc));
                }
                if (index[1] < (m_uiElementOrder + 1))
                {
                    idy[k] = index[1];
                    assert((parent.getY() + idy[k] * Lp / dp) == (child.getY() + k * Lc / dc));
                }
                if (index[2] < (m_uiElementOrder + 1))
                {
                    idz[k] = index[2];
                    assert((parent.getZ() + idz[k] * Lp / dp) == (child.getZ() + k * Lc / dc));
                }
            }
            state = stateX & stateY & stateZ;
            return state;
        }
    }
    else
    {
        return false;
    }
}

} // namespace ot

#include "mesh.tcc"
#include "meshE2NUtils.tcc"

#endif //SFCSORTBENCH_MESH_H