oda.h

Go to the documentation of this file.

#ifndef __OCT_DA_H__
#define __OCT_DA_H__

#include "mpi.h"
#include "matRecord.h"
#include "loopCounters.h"
#include "updateCtx.h"
#include "odaUtils.h"
#include "cnumEtypes.h"
#include "Point.h"
#include <vector>
#include "petsc.h"
#include "petscvec.h"
#include "petscmat.h"
#include "dendro.h"

#ifndef iC
#define iC(fun) {CHKERRQ(fun);}
#endif

#ifdef __DEBUG__
#ifndef __DEBUG_DA__
#define __DEBUG_DA__
#endif
#endif

#ifdef PETSC_USE_LOG

#include "petscsys.h"

namespace ot {
  extern int readFromGhostNodesBeginEvent;
  extern int readFromGhostNodesEndEvent;
  extern int readFromGhostElemsBeginEvent;
  extern int readFromGhostElemsEndEvent;
  extern int writeToGhostNodesBeginEvent;
  extern int writeToGhostNodesEndEvent;
  extern int writeToGhostElemsBeginEvent;
  extern int writeToGhostElemsEndEvent;
  extern int buildDaEvent;
  extern int setMatValuesEvent;
  extern int buildNlistEvent;
  extern int buildNlistCommEvent;

  extern int buildDAstage1Event;
  extern int buildDAstage2Event;
  extern int buildDAstage3Event;
  extern int buildDAstage4Event;
  extern int buildDAstage5Event;
  extern int buildDAstage6Event;
  extern int buildDAstage7Event;
  extern int buildDAstage8Event;
  extern int buildDAstage9Event;
}

#define PROF_WRITE_GHOST_ELEMS_BEGIN_BEGIN PetscLogEventBegin(writeToGhostElemsBeginEvent,0,0,0,0);
#define PROF_WRITE_GHOST_ELEMS_BEGIN_END \
  PetscLogEventEnd(writeToGhostElemsBeginEvent,0,0,0,0); \
return 0;

#define PROF_WRITE_GHOST_ELEMS_END_BEGIN PetscLogEventBegin(writeToGhostElemsEndEvent,0,0,0,0);
#define PROF_WRITE_GHOST_ELEMS_END_END \
  PetscLogEventEnd(writeToGhostElemsEndEvent,0,0,0,0); \
return 0;

#define PROF_WRITE_GHOST_NODES_BEGIN_BEGIN PetscLogEventBegin(writeToGhostNodesBeginEvent,0,0,0,0);
#define PROF_WRITE_GHOST_NODES_BEGIN_END \
  PetscLogEventEnd(writeToGhostNodesBeginEvent,0,0,0,0); \
return 0;

#define PROF_WRITE_GHOST_NODES_END_BEGIN PetscLogEventBegin(writeToGhostNodesEndEvent,0,0,0,0);
#define PROF_WRITE_GHOST_NODES_END_END \
  PetscLogEventEnd(writeToGhostNodesEndEvent,0,0,0,0); \
return 0;

#define PROF_READ_GHOST_ELEMS_BEGIN_BEGIN PetscLogEventBegin(readFromGhostElemsBeginEvent,0,0,0,0);
#define PROF_READ_GHOST_ELEMS_BEGIN_END \
  PetscLogEventEnd(readFromGhostElemsBeginEvent,0,0,0,0); \
return 0;

#define PROF_READ_GHOST_ELEMS_END_BEGIN PetscLogEventBegin(readFromGhostElemsEndEvent,0,0,0,0);
#define PROF_READ_GHOST_ELEMS_END_END \
  PetscLogEventEnd(readFromGhostElemsEndEvent,0,0,0,0); \
return 0;

#define PROF_READ_GHOST_NODES_BEGIN_BEGIN PetscLogEventBegin(readFromGhostNodesBeginEvent,0,0,0,0);
#define PROF_READ_GHOST_NODES_BEGIN_END \
  PetscLogEventEnd(readFromGhostNodesBeginEvent,0,0,0,0); \
return 0;

#define PROF_READ_GHOST_NODES_END_BEGIN PetscLogEventBegin(readFromGhostNodesEndEvent,0,0,0,0);
#define PROF_READ_GHOST_NODES_END_END \
  PetscLogEventEnd(readFromGhostNodesEndEvent,0,0,0,0); \
return 0;

#define PROF_SET_MAT_VALUES_BEGIN PetscLogEventBegin(setMatValuesEvent,0,0,0,0);
#define PROF_SET_MAT_VALUES_END \
  PetscLogEventEnd(setMatValuesEvent,0,0,0,0); \
return 0;

#define PROF_BUILD_DA_BEGIN PetscLogEventBegin(buildDaEvent,0,0,0,0);
#define PROF_BUILD_DA_END \
  PetscLogEventEnd(buildDaEvent,0,0,0,0); \
return;

#define PROF_BUILD_NLIST_BEGIN PetscLogEventBegin(buildNlistEvent,0,0,0,0);
#define PROF_BUILD_NLIST_END \
  PetscLogEventEnd(buildNlistEvent,0,0,0,0); \
return ;

#define PROF_BUILD_NLIST_COMM_BEGIN PetscLogEventBegin(buildNlistCommEvent,0,0,0,0);
#define PROF_BUILD_NLIST_COMM_END PetscLogEventEnd(buildNlistCommEvent,0,0,0,0); 

#define PROF_BUILD_DA_STAGE1_BEGIN PetscLogEventBegin(buildDAstage1Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE1_END PetscLogEventEnd(buildDAstage1Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE2_BEGIN PetscLogEventBegin(buildDAstage2Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE2_END PetscLogEventEnd(buildDAstage2Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE3_BEGIN PetscLogEventBegin(buildDAstage3Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE3_END PetscLogEventEnd(buildDAstage3Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE4_BEGIN PetscLogEventBegin(buildDAstage4Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE4_END PetscLogEventEnd(buildDAstage4Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE5_BEGIN PetscLogEventBegin(buildDAstage5Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE5_END PetscLogEventEnd(buildDAstage5Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE6_BEGIN PetscLogEventBegin(buildDAstage6Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE6_END PetscLogEventEnd(buildDAstage6Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE7_BEGIN PetscLogEventBegin(buildDAstage7Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE7_END PetscLogEventEnd(buildDAstage7Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE8_BEGIN PetscLogEventBegin(buildDAstage8Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE8_END PetscLogEventEnd(buildDAstage8Event,0,0,0,0); 

#define PROF_BUILD_DA_STAGE9_BEGIN PetscLogEventBegin(buildDAstage9Event,0,0,0,0);
#define PROF_BUILD_DA_STAGE9_END PetscLogEventEnd(buildDAstage9Event,0,0,0,0); 

#else

#define PROF_READ_GHOST_NODES_BEGIN_BEGIN 
#define PROF_READ_GHOST_NODES_BEGIN_END return 0;

#define PROF_READ_GHOST_NODES_END_BEGIN 
#define PROF_READ_GHOST_NODES_END_END return 0;

#define PROF_READ_GHOST_ELEMS_BEGIN_BEGIN 
#define PROF_READ_GHOST_ELEMS_BEGIN_END return 0;

#define PROF_READ_GHOST_ELEMS_END_BEGIN 
#define PROF_READ_GHOST_ELEMS_END_END return 0;

#define PROF_WRITE_GHOST_NODES_BEGIN_BEGIN 
#define PROF_WRITE_GHOST_NODES_BEGIN_END return 0;

#define PROF_WRITE_GHOST_NODES_END_BEGIN 
#define PROF_WRITE_GHOST_NODES_END_END return 0;

#define PROF_WRITE_GHOST_ELEMS_BEGIN_BEGIN 
#define PROF_WRITE_GHOST_ELEMS_BEGIN_END return 0;

#define PROF_WRITE_GHOST_ELEMS_END_BEGIN 
#define PROF_WRITE_GHOST_ELEMS_END_END return 0;

#define PROF_SET_MAT_VALUES_BEGIN 
#define PROF_SET_MAT_VALUES_END return 0;

#define PROF_BUILD_DA_BEGIN 
#define PROF_BUILD_DA_END return;

#define PROF_BUILD_NLIST_BEGIN 
#define PROF_BUILD_NLIST_END return ;

#define PROF_BUILD_NLIST_COMM_BEGIN 
#define PROF_BUILD_NLIST_COMM_END 

#define PROF_BUILD_DA_STAGE1_BEGIN 
#define PROF_BUILD_DA_STAGE1_END 

#define PROF_BUILD_DA_STAGE2_BEGIN 
#define PROF_BUILD_DA_STAGE2_END 

#define PROF_BUILD_DA_STAGE3_BEGIN 
#define PROF_BUILD_DA_STAGE3_END 

#define PROF_BUILD_DA_STAGE4_BEGIN 
#define PROF_BUILD_DA_STAGE4_END 

#define PROF_BUILD_DA_STAGE5_BEGIN 
#define PROF_BUILD_DA_STAGE5_END 

#define PROF_BUILD_DA_STAGE6_BEGIN 
#define PROF_BUILD_DA_STAGE6_END 

#define PROF_BUILD_DA_STAGE7_BEGIN 
#define PROF_BUILD_DA_STAGE7_END 

#define PROF_BUILD_DA_STAGE8_BEGIN 
#define PROF_BUILD_DA_STAGE8_END 

#define PROF_BUILD_DA_STAGE9_BEGIN 
#define PROF_BUILD_DA_STAGE9_END 

#endif

namespace ot {

  //Forward Declaration
  class TreeNode;

    class DA {

      protected:

        LoopCounters m_lcLoopInfo;

        //The compressed Octree. Only the levels of the octants are stored.
        unsigned char*          m_ucpOctLevels;
        bool                    m_bCompressLut;
        bool                    m_bComputedLocalToGlobal;
        bool                    m_bComputedLocalToGlobalElems;
        bool                    m_bIamActive;
        unsigned int            m_uiInputSize;
        std::vector<ot::TreeNode>          m_tnBlocks;
        std::vector<ot::TreeNode>          m_tnMinAllBlocks;
        std::vector<unsigned int>          m_uiNlist;  
        unsigned int*                      m_uiNlistPtr;
        DendroIntL*                      m_dilpLocalToGlobal;
        DendroIntL*                      m_dilpLocalToGlobalElems;

        // The remainders from the Golomb-Rice encoding of the neighbour
        // look-up-table. 
        std::vector<unsigned char>          m_ucpLutRemainders;
        unsigned char*                      m_ucpLutRemaindersPtr;
        std::vector<unsigned char>          m_ucpSortOrders;
        unsigned char*                      m_ucpSortOrdersPtr;

        // The quotient from the Golomb-Rice encoding of the neighbour
        // look-up-table. 
        std::vector<unsigned short>         m_uspLutQuotients;
        unsigned short*                     m_uspLutQuotientsPtr;

        // The masks for the look-up table ...
        // the size of the mask will be 2*nelem, the first character storing the 
        // compression mask and the second the hanging node mask.
        std::vector<unsigned char>      m_ucpLutMasks;
        unsigned char*                  m_ucpLutMasksPtr;

        std::vector<unsigned char>      m_ucpPreGhostConnectivity;
        unsigned char*                  m_ucpPreGhostConnectivityPtr;

        std::vector<Point>              m_ptsPreGhostOffsets;
        Point*                          m_ptsPreGhostOffsetsPtr;

        std::vector<Point>              PreGhostAnchors;

        // The counter that maintains the current
        //quotient index while looping through the DA.

        unsigned int                    m_uiQuotientCounter;
        unsigned int                    m_uiPreGhostQuotientCnt;

        // The quotients for element begin and independent begin.
        unsigned int                    m_uiElementQuotient;
        unsigned int                    m_uiIndependentElementQuotient;

        // The offset (x,y,z) for the region spanned by the
        //octants on the current processor, including ghost octants.
        Point                           m_ptGhostedOffset;

        // The offset (x,y,z) for the region spanned by the octants on the current
        // processor.
        Point                           m_ptOffset;

        Point                           m_ptIndependentOffset;
        Point                           m_ptCurrentOffset;

        //------------------------------------------------------------------------

        // The number of nodes owned by the current processor.
        unsigned int                    m_uiNodeSize;
        unsigned int                    m_uiBoundaryNodeSize;
        // The number of elements owned by the current processor.
        unsigned int                    m_uiElementSize;

        unsigned int                    m_uiIndependentElementSize;
        // The number of pre-ghosts elements, i.e., ghost elements with morton ids less than
        //  those belonging to this processor.
        unsigned int                    m_uiPreGhostElementSize;

        // The number of pre-ghosts nodes, i.e., ghost nodes with morton ids less than
        //  those belonging to this processor.
        unsigned int                    m_uiPreGhostNodeSize;
        unsigned int                    m_uiPreGhostBoundaryNodeSize;
        // The number of post-ghosts nodes, i.e., ghost nodes with morton ids greater
        //  than those belonging to this processor.
        unsigned int                    m_uiPostGhostNodeSize;    

        unsigned int                    m_uiLocalBufferSize;

        // indices ...
        unsigned int                           m_uiElementBegin;
        unsigned int                           m_uiElementEnd;
        unsigned int                           m_uiPostGhostBegin;

        unsigned int                           m_uiIndependentElementBegin;
        unsigned int                           m_uiIndependentElementEnd;

        unsigned int                           m_uiCurrent;
        //------------------------------------------------------------------------

        // The dimensionality of the underlying space. Only tested in 3D.
        unsigned int                    m_uiDimension;
        // The maximum depth of the octree discretization. Is usually obtained from
        // the octree generation and balancing procedure.
        unsigned int                    m_uiMaxDepth;

        // The MPI Communication context. 
        MPI_Comm                        m_mpiCommAll;
        int                             m_iRankAll;
        int                             m_iNpesAll;
        MPI_Comm                        m_mpiCommActive;
        int                             m_iRankActive;
        int                             m_iNpesActive;
        // Stuff for communication and synchronization ...
        //-------------------------------------------------------------------------
        std::vector<unsigned int>               m_uipScatterMap;
        std::vector<unsigned int>               m_uipSendOffsets;
        std::vector<unsigned int>               m_uipSendProcs;
        std::vector<unsigned int>               m_uipSendCounts;

        std::vector<unsigned int>               m_uipElemScatterMap;
        std::vector<unsigned int>               m_uipElemSendOffsets;
        std::vector<unsigned int>               m_uipElemSendProcs;
        std::vector<unsigned int>               m_uipElemSendCounts;

        std::vector<unsigned int>               m_uipRecvOffsets;
        std::vector<unsigned int>               m_uipRecvProcs;
        std::vector<unsigned int>               m_uipRecvCounts;

        std::vector<unsigned int>               m_uipElemRecvOffsets;
        std::vector<unsigned int>               m_uipElemRecvProcs;
        std::vector<unsigned int>               m_uipElemRecvCounts;

        // tags and contexts ...
        std::vector<updateContext>              m_mpiContexts;
        unsigned int                            m_uiCommTag;

      public:

        std::vector<ot::TreeNode> getBlocks();

        unsigned int getNumBlocks();

        std::vector<ot::TreeNode> getMinAllBlocks();

        unsigned int getInputSize();


        DA(std::vector<ot::TreeNode> &in, MPI_Comm comm, MPI_Comm inputActiveComm,  
            bool compressLut = false, const std::vector<ot::TreeNode> *blocksPtr = NULL,
            bool *iAmActive = NULL);

        DA(unsigned int dummy, std::vector<ot::TreeNode> &in, MPI_Comm comm, 
            MPI_Comm inputActiveComm,  bool compressLut = false,
            const std::vector<ot::TreeNode> *blocksPtr = NULL, bool *iAmActive = NULL);

        ~DA();


        MPI_Comm getComm();

        MPI_Comm getCommActive();

        int getNpesAll();

        int getNpesActive();

        int getRankAll();

        int getRankActive();

        Point getCurrentOffset();

        Point getGhostedOffset();

        Point getOffset();

        Point getNextOffset(Point p, unsigned char d);

        void incrementCurrentOffset();

        void incrementPreGhostOffset();

        bool isBoundaryOctant(unsigned char *flags=NULL);

        DendroIntL* getLocalToGlobalMap();

        DendroIntL* getLocalToGlobalElemsMap();

        unsigned int getNodeSize();

        unsigned int getBoundaryNodeSize();

        unsigned int getInternalNodeSize();

        unsigned int getElementSize();

        unsigned int getPreGhostElementSize();

        unsigned int getIndependentSize();

        unsigned int getGhostedNodeSize();

        unsigned int getGhostedElementSize();

        unsigned int getIdxElementBegin();

        unsigned int getIdxElementEnd();

        unsigned int getIdxPostGhostBegin();

        unsigned int getMaxDepth();

        unsigned int getDimension();



        //Communicating Ghost Nodes
        template <typename T>
          int ReadFromGhostsBegin ( T* arr, unsigned int dof=1);

        template <typename T>
          int ReadFromGhostsEnd(T* arr);


        template <typename T>
          int WriteToGhostsBegin ( T* arr, unsigned int dof=1);

        template <typename T>
          int WriteToGhostsEnd(T* arr, unsigned int dof=1);


        template <typename T>
          int ReadFromGhostElemsBegin ( T* arr, unsigned int dof=1);

        template <typename T>
          int ReadFromGhostElemsEnd(T* arr);

        template <typename T>
          int WriteToGhostElemsBegin ( T* arr, unsigned int dof=1);

        template <typename T>
          int WriteToGhostElemsEnd(T* arr, unsigned int dof=1);



        int createVector(Vec &local, bool isElemental, bool isGhosted, unsigned int dof=1);

        int createActiveVector(Vec &local, bool isElemental, bool isGhosted, unsigned int dof=1);

        int createMatrix(Mat &M, MatType mtype, unsigned int dof=1);

        int createActiveMatrix(Mat &M, MatType mtype, unsigned int dof=1);

        int computeLocalToGlobalMappings();

        int computeLocalToGlobalElemMappings();

        bool computedLocalToGlobal();

        bool computedLocalToGlobalElems();

        int setValuesInMatrix(Mat mat, std::vector<ot::MatRecord>& records,
            unsigned int dof, InsertMode mode);

        template <typename T>
          int  createVector(std::vector<T> &local, bool isElemental,
              bool isGhosted, unsigned int dof=1);

        int vecGetBuffer(Vec in, PetscScalar* &out, bool isElemental,
            bool isGhosted, bool isReadOnly, unsigned int dof=1); 

        template < typename T >
          int vecGetBuffer(std::vector<T> &in, T* &out, bool isElemental,
              bool isGhosted, bool isReadOnly, unsigned int dof=1); 

        int vecRestoreBuffer(Vec in, PetscScalar* out, bool isElemental, 
            bool isGhosted, bool isReadOnly, unsigned int dof=1); 

        template < typename T >
          int vecRestoreBuffer(std::vector<T> &in, T* out, bool isElemental,
              bool isGhosted, bool isReadOnly, unsigned int dof=1); 

        //----------------------------------


        template<ot::DA_FLAGS::loopType type>
          void init();

        unsigned int curr();

        unsigned int currWithInfo();

        template<ot::DA_FLAGS::loopType type>
          unsigned int end();

        template<ot::DA_FLAGS::loopType type>
          unsigned int next();

        unsigned char getChildNumber();

        unsigned char getFlag(unsigned int i);

        bool isHanging(unsigned int i);

        bool isGhost(unsigned int i);

        bool isNode(unsigned int i);

        unsigned char getHangingNodeIndex(unsigned int i);

        unsigned char getTypeMask(unsigned int i);

        unsigned char getLevel(unsigned int i);

        int getNodeIndices(unsigned int* nodes);

        unsigned int getLocalBufferSize(); 


        void updateQuotientCounter();

        bool isLUTcompressed();

        bool iAmActive();

      protected:

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

        void DA_FactoryPart0(std::vector<ot::TreeNode>& in, MPI_Comm comm,
            MPI_Comm activeInputComm, bool compressLut, bool* iAmActive);

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

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

        void DA_FactoryPart3(std::vector<ot::TreeNode>& in, MPI_Comm comm, bool compressLut, 
            const std::vector<ot::TreeNode>* blocksPtr, bool* iAmActive);

    };//end DA class definition

}// end namespace

#include "oda.txx"

#endif

---