oda.h

Go to the documentation of this file.

#ifndef DENDRO_5_ODA_H
#define DENDRO_5_ODA_H

#include "dendro.h"
#include <iostream>
#include <vector>
#include "mpi.h"
#include "mesh.h"
#include "parUtils.h"
#include "testUtils.h"
#include "asyncExchangeContex.h"
#include "oct2vtk.h"
#include "odaUtils.h"
#include "matRecord.h"
#include <assert.h>

#ifdef BUILD_WITH_PETSC
    #include "petsc.h"
    #include "petscvec.h"
    #include "petscmat.h"
    #include "petscdmda.h"
#endif

#define VECType DendroScalar
#define MATType DendroScalar

/***
 * Note: All the Dendro vectors are stored like as shown below in the memory.
 * |------- pre ghost nodes----- |-----------local nodes---------------|--------post ghost nodes--------|
 *
 *
 *
 * */

// ODA flag information.
namespace ot
{
    namespace DA_FLAGS
    {
        enum LoopType {ALL,WRITABLE,INDEPENDENT,W_DEPENDENT,W_BOUNDARY,LOCAL_ELEMENTS,PREGHOST_ELEMENTS,POSTGHOST_ELEMENTS};

        enum BdyType {X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX};


        enum Refine {DA_NO_CHANGE,DA_REFINE,DA_COARSEN};

        enum WriteMode{SET_VALUES=0,ADD_VALUES};



    }

}




namespace ot
{
    class Mesh;
    class TreeNode;

    enum DAType {OCT,PETSC};

    enum MVECType {MESH_BASED,MESH_FREE};


    struct LoopCounter
    {   
        unsigned int currentIndex;
        unsigned int indexBegin;
        unsigned int indexEnd;
    };

    class DA
    {

    private:
        LoopCounter m_uiLoopInfo;

        double m_uiOctreeLowerBound[m_uiDim];

        double m_uiOctreeUpperBound[m_uiDim];

        std::vector<unsigned int> m_uiOctantFlags;

        ot::Mesh* m_uiMesh;


        unsigned int m_uiTotalNodalSz;
        unsigned int m_uiLocalNodalSz;

        unsigned int m_uiTotalElementSz;

        unsigned int m_uiLocalElementSz;

        std::vector<ot::AsyncExchangeContex> m_uiMPIContexts;

        unsigned int m_uiCommTag;

        std::vector<DendroIntL> m_uiLocalToGlobalNodalMap;

        DendroIntL m_uiGlobalNodeSz;




    public:
        std::vector<ot::TreeNode> getLocalOctants();

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

        DA(std::vector<ot::TreeNode> &balOct,MPI_Comm comm,unsigned int order, unsigned int grainSz=100,double sfc_tol=0.3, SM_TYPE smType=SM_TYPE::FEM_CG);

        DA(ot::Mesh* pMesh);


        template<typename T>
        DA(std::function<void(T,T,T,T*)>func,unsigned int dofSz,MPI_Comm comm,unsigned int order,double interp_tol, unsigned int grainSz=100,double sfc_tol=0.3, SM_TYPE smType=SM_TYPE::FEM_CG);

        DA(const Point* pts, unsigned int numPts,Point pt_min,Point pt_max, MPI_Comm comm,unsigned int order,unsigned int grainSz=100,double sfc_tol=0.3, SM_TYPE smType=SM_TYPE::FEM_CG);



        ~DA();

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

        inline unsigned int getPreNodalSz() const {return m_uiMesh->getNumPreMeshNodes();}

        inline unsigned int getPostNodalSz() const {return m_uiMesh->getNumPostMeshNodes();}

        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_uiMesh->isActive();}

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

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

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

        inline const ot::Mesh* getMesh() const {return m_uiMesh;}

        inline const std::vector<unsigned int> getOctFlags() const {return m_uiOctantFlags; }

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

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

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

        inline unsigned int getRankActive() const
        {
            if(m_uiMesh->isActive())
                return m_uiMesh->getMPIRank();
            else
                return m_uiMesh->getMPIRankGlobal();
        }

        bool isBoundaryOctant(unsigned int eleID) const;

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


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


        inline unsigned int getElementSize() const {return m_uiMesh->getNumLocalMeshElements();}

        inline unsigned int getPreGhostElementSize() const {return m_uiMesh->getNumPreGhostElements();}

        inline unsigned int getPostGhostElementSize() const {return m_uiMesh->getNumPostGhostElements(); }

        inline unsigned int getPreAndPostGhostNodeSize() const {return (m_uiMesh->getNumPreMeshNodes() + m_uiMesh->getNumPostMeshNodes());}

        inline unsigned int getNumLayer1GhostEleSz() const {return m_uiMesh->getLevel1GhostElementIndices().size();}

        unsigned int getGhostedElementSize() const {return (m_uiMesh->getNumPreGhostElements() + m_uiMesh->getNumPostGhostElements());}


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

        inline ot::TreeNode getMinTreeNode() const { 

            if(m_uiMesh->isActive())
            {
                const ot::TreeNode* allElements = &(*(m_uiMesh->getAllElements().begin()));
                return allElements[m_uiMesh->getElementLocalBegin()];
            }else
            {
                return ot::TreeNode(0,0,0,0,m_uiDim,m_uiMaxDepth);                
            }
            
        }
        
        inline ot::TreeNode getMaxTreeNode() const { 
            
            if(m_uiMesh->isActive())
            {
                const ot::TreeNode* allElements = &(*(m_uiMesh->getAllElements().begin()));
                return allElements[m_uiMesh->getElementLocalEnd()];
            }else
            {
                return ot::TreeNode(0,0,0,0,m_uiDim,m_uiMaxDepth);                
            }

            
        }

        void computeTreeNodeOwnerProc(const ot::TreeNode * pNodes, unsigned int n, int* ownerranks);

        
         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, DA_FLAGS::WriteMode mode,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_uiMesh->getAllElements()[ele].getLevel();}

        inline ot::TreeNode getOctant(unsigned int ele) const {return m_uiMesh->getAllElements()[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, unsigned int (*getWeight)(const ot::TreeNode *)=NULL) 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,unsigned int dof) const;

        inline unsigned int getMortonChildNum(unsigned int eleID) const
        {
           return m_uiMesh->getMortonchildNum(eleID);
        }

         // 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 "oda.tcc"


#endif