sub_oda.h

/***
 * @file        sub_oda.h
 * @author      Hari Sundar   hsundar@gmail.com
 * @author      Milinda Fernando milinda@cs.utah.edu
 * @date        13 June 2019
 * @brief       Subdomain for Octree DA (ODA).
 ***/

#pragma once

#include <functional>
#include <iostream>

#include "oda.h"

namespace ot {
  
  class subDA {

      protected:

        ot::DA*       m_da;
        std::vector<unsigned char>  m_ucpSkipList;
        std::vector<unsigned char>  m_ucpSkipNodeList;
        unsigned int  m_uiLocalElementSize;
        unsigned int  m_uiPreGhostElementSize;
        unsigned int  m_uiPostGhostElementSize;
        unsigned int  m_uiLocalNodeSize;
        unsigned int  m_uiPreGhostNodeSize;
        unsigned int  m_uiPostGhostNodeSize;
        unsigned int m_uiCommTag;

        std::vector<AsyncExchangeContex> m_uiMPIContexts;

        LoopCounter m_uiLoopInfo;
        
        std::vector<AsyncExchangeContex>  m_mpiContexts;
        std::vector<unsigned int> m_uip_sub2DA_ElemMap;      
        std::vector<unsigned int> m_uip_DA2sub_ElemMap;      
        std::vector<unsigned int> m_uip_sub2DA_NodeMap;      
        std::vector<unsigned int> m_uip_DA2sub_NodeMap;

        double   m_dMinBB[3];
        double   m_dMaxBB[3];

        unsigned int m_uiTotalNodalSz;
        unsigned int m_uiLocalNodalSz;

        unsigned int m_uiTotalElementSz;

        unsigned int m_uiLocalElementSz;

        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_uiNodePreGhostBegin=0;
        unsigned int m_uiNodePreGhostEnd;
        unsigned int m_uiNodeLocalBegin;
        unsigned int m_uiNodeLocalEnd;
        unsigned int m_uiNodePostGhostBegin;
        unsigned int m_uiNodePostGhostEnd;

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

        unsigned int m_uiNumPreGhostNodes;
        unsigned int m_uiNumLocalNodes;
        unsigned int m_uiNumPostGhostNodes;

        std::vector<unsigned int > m_uiSendScatterMap;
        std::vector<unsigned int > m_uiRecvScatterMap;
        
        std::vector<unsigned int> m_uiSendCounts;
        std::vector<unsigned int> m_uiRecvCounts;
        std::vector<unsigned int> m_uiSendOffsets;
        std::vector<unsigned int> m_uiRecvOffsets;
        std::vector<unsigned int> m_uiSendProcList;
        std::vector<unsigned int> m_uiRecvProcList;

        bool m_uiIsActive;

        MPI_Comm m_uiCommGlobal;

        MPI_Comm m_uiCommActive;

        int m_uiRankActive;
        
        int m_uiRankGlobal;
        
        int m_uiNpesActive;

        int m_uiNpesGlobal;

       
      public:
        subDA(DA* da, std::function<double ( double, double, double ) > fx_retain, double* gSize);

        ~subDA();
        
        ot::DA* global_domain() { return m_da; }   

        inline unsigned int getLocalNodalSz() const {return m_uiLocalNodalSz;}

        inline unsigned int getPreNodalSz() const { return m_uiNumPreGhostNodes; }

        inline unsigned int getPostNodalSz() const {return m_uiNumPostGhostNodes;}

        inline unsigned int getTotalNodalSz() const {return m_uiTotalNodalSz;}

        inline unsigned int getLocalElemSz() const {return m_uiLocalElementSz;}

        inline unsigned int getTotalElemSz() const {return m_uiTotalElementSz;}

        inline bool isActive() const {return m_da->isActive();}

        inline unsigned int getNumNodesPerElement() const {return m_da->getNumNodesPerElement();}

        inline unsigned int getElementOrder() const {return m_da->getElementOrder();}

        inline MPI_Comm getGlobalComm() const { return m_da->getGlobalComm();}
        
        //inline const std::vector<DendroIntL> getNodeLocalToGlobalMap() const {return m_uiLocalToGlobalNodalMap;}

        inline const ot::Mesh* getMesh() const {return m_da->getMesh();}

        inline MPI_Comm getCommActive() const
        {
           if(m_da->isActive())
               return m_da->getCommActive();
           else
               return MPI_COMM_NULL;
        }

        inline unsigned int getNpesAll() const { return  m_da->getNpesAll();} ;
        
        inline unsigned int getNpesActive() const
        {
            if(m_da->isActive())
                return m_da->getNpesActive();
            else
                return 0;
        }

        inline unsigned int getRankAll() const {return m_da->getRankAll();};

        inline unsigned int getRankActive() const
        {
            if(m_da->isActive())
                return m_da->getRankActive();
            else
                return m_da->getRankAll();
        }

        bool isBoundaryOctant(unsigned int eleID) const;

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


        inline const RefElement* getReferenceElement()const { return  m_da->getReferenceElement();}


        inline unsigned int getElementSize() const {return m_uiNumLocalElements;}

        inline unsigned int getPreGhostElementSize() const {return m_uiNumPreGhostElements;}

        inline unsigned int getPostGhostElementSize() const {return m_uiNumPostGhostElements; }

        inline unsigned int getPreAndPostGhostNodeSize() const {return (m_uiNumPreGhostNodes + m_uiNumPostGhostNodes);}
        
        unsigned int getGhostedElementSize() const {return (m_uiNumPreGhostElements + m_uiNumPostGhostElements);}

        inline unsigned int getMaxDepth() const {return  m_uiMaxDepth;};
        inline unsigned int getDimension() const {return m_uiDim;};

         template<typename T>
         int createVector(T*& local, bool isElemental=false, bool isGhosted=false, unsigned int dof=1) const;

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

        template <typename T>
        void destroyVector(T*& local) const;

        template <typename T>
        void destroyVector(std::vector<T>& local) const;



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

         unsigned int curr();


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

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

         template<typename T>
         void readFromGhostBegin(T* vec, unsigned int dof=1);

         template<typename T>
         void readFromGhostEnd(T* vec, unsigned int dof=1);


         template<typename T>
         void writeToGhostsBegin(T* vec, unsigned int dof=1) ;

         template<typename T>
         void writeToGhostsEnd(T* vec, unsigned int dof=1) ;


         template<typename T>
         void nodalVecToGhostedNodal(const T* in,T*& out,bool isAllocated=false, unsigned int dof=1)const;


         template<typename T>
         void ghostedNodalToNodalVec(const T* gVec,T*& local,bool isAllocated=false, unsigned int dof=1) const;



         template <typename T>
         void getElementNodalValues(const T*in, T* eleVecOut,unsigned int eleID,unsigned int dof=1)const;

        template <typename T>
        void eleVecToVecAccumilation(T*out, const T* eleVecIn,unsigned int eleID,unsigned int dof=1)const;

        void getOctreeBoundaryNodeIndices(std::vector<unsigned int >& bdyIndex,std::vector<double>& coords,bool isGhosted=false);

        int getNodeIndices(DendroIntL* nodeIdx,unsigned int ele,bool isGhosted) const;
        
        int getGlobalNodeIndices(DendroIntL* nodeIdx,unsigned int ele) const ;

        template <typename T>
        void setVectorByFunction(T* local,std::function<void(T,T,T,T*)>func,bool isElemental=false, bool isGhosted=false, unsigned int dof=1) const;


        template <typename T>
        void setVectorByScalar(T* local,const T* value,bool isElemental=false, bool isGhosted=false, unsigned int dof=1) const;


        template <typename T>
        void vecTopvtu(T* local, const char * fPrefix,char** nodalVarNames=NULL,bool isElemental=false,bool isGhosted=false,unsigned int dof=1);


        inline unsigned int getLevel(unsigned int ele) const { return m_da->getLevel(m_uip_sub2DA_ElemMap[ele]);}

        inline ot::TreeNode getOctant(unsigned int ele) const {return m_da->getOctant(m_uip_sub2DA_ElemMap[ele]);}

        template<typename T>
        T* getVecPointerToDof(T* in ,unsigned int dofInex, bool isElemental=false,bool isGhosted=false) const;


        template<typename T>
        void copyVectors(T* dest,const T* source,bool isElemental=false,bool isGhosted=false,unsigned int dof=1) const;

        template<typename T>
        void copyVector(T* dest,const T* source,bool isElemental=false,bool isGhosted=false) const;


        ot::DA* remesh(const DA_FLAGS::Refine * flags, unsigned int sz,unsigned int grainSz=100,double ld_bal=0.3, unsigned int sfK=2) const;

        template<typename T>
        void intergridTransfer(const T* varIn, T* & varOut, const ot::DA* newDA, bool isElemental=false, bool isGhosted=false, unsigned int dof=1);



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

        
         // all the petsc functionalities goes below with the pre-processor gards.
#ifdef BUILD_WITH_PETSC

        PetscErrorCode petscCreateVector(Vec &local, bool isElemental, bool isGhosted, unsigned int dof) const;


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



        PetscErrorCode petscNodalVecToGhostedNodal(const Vec& in,Vec& out,bool isAllocated=false,unsigned int dof=1) const;


        PetscErrorCode petscGhostedNodalToNodalVec(const Vec& gVec,Vec& local,bool isAllocated=false,unsigned int dof=1) const;

        void petscReadFromGhostBegin(PetscScalar * vecArry, unsigned int dof=1) ;

        void petscReadFromGhostEnd(PetscScalar * vecArry, unsigned int dof=1) ;


        template<typename T>
        void petscSetVectorByFunction(Vec& local,std::function<void(T,T,T,T*)>func,bool isElemental=false, bool isGhosted=false, unsigned int dof=1) const;


        template <typename T>
        void petscSetVectorByScalar(Vec& local,const T* value,bool isElemental=false, bool isGhosted=false, unsigned int dof=1) const;


        void petscVecTopvtu(const Vec& local, const char * fPrefix,char** nodalVarNames=NULL,bool isElemental=false,bool isGhosted=false,unsigned int dof=1) ;


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

        
        
        PetscErrorCode petscChangeVecToMatBased(Vec& v1,bool isElemental,bool isGhosted,unsigned int dof=1) const;
        
        
        PetscErrorCode petscChangeVecToMatFree(Vec& v1,bool isElemental,bool isGhosted,unsigned int dof=1) const;

        template<typename T>
        void petscIntergridTransfer(const Vec &varIn, Vec &varOut, const ot::DA *newDA, bool isElemental = false,bool isGhosted = false, unsigned int dof = 1);

        PetscErrorCode petscDestroyVec(Vec & vec);
      

#endif



    };


}

#include "sub_oda.tcc"