EMM_Grid3D.h

Go to the documentation of this file.

#ifndef EMM_Grid3D_H
#define EMM_Grid3D_H

#include "EMM_Object.h"
#include "EMM_Array.h"
#include "CORE_MorseArray.h"

class EMM_MatterField;

DEFINE_SPTR(EMM_Grid3D);
class EMM_Grid3D : public EMM_Object  {
  
    SP_OBJECT(EMM_Grid3D);
    // ATTRIBUTES

public:
    static const tUIndex UNMAGNETIZED_ELEMENT;
    
    //limit condition types
    static const tLimitCondition DIRICHLET_LIMIT_CONDITION=1;
    static const tLimitCondition NEUMANN_LIMIT_CONDITION=0;
    static const tLimitCondition NO_LIMIT_CONDITION=2;
    static const tLimitCondition SLAVE_LIMIT_CONDITION=3;


    //points of mesh
    static const tFlag POINT=0;
    
    //element of mesh
    static const tFlag ELEMENT=1;
    
    //indicator for element in the interior of domain pow(2,FACES_NUMBER_PER_ELEMENT+1)-1=127
    //all bits in {0,1,2,3,n4,5} is 1
    static const tCellFlag MAGNETIZED_DOMAIN_INTERIOR_INDICATOR=127;
    
    //min indicator for element in magnetized domain pow(2,FACES_NUMBER_PER_ELEMENT)=64
    static const tCellFlag MAGNETIZED_ELEMENT_MIN_INDICATOR=64;
    
   
    //number of faces per element
    static const tUCInt FACES_NUMBER_PER_ELEMENT=6;
    
    //number of points by element
    static const tUCInt POINTS_NUMBER_PER_ELEMENT=8;

    //number of points by face
    static const tUCInt POINTS_NUMBER_PER_FACE=4;
    
    static const tBoolean ELEMENT_POINTS[];
    
    static const tUCInt FACE_POINTS[];

    static const tCInt EDGE_POINTS[];
    
private:

   
   
    
    //2^iFace iFace in [0,6[
    static const tCellFlag TWO_POWER_FACES[];
    
    tBoolean mIsDirectionPeriodic[3];
    
    tBoolean mIsPeriodic;
    
    tUInteger mSegmentsNumber[3];
    
    tReal mSegmentsSize[3];
    tReal mAdimensionizedSegmentsSize[3];

    tReal mAdimensionizedVolume;
    
     tReal mLmin;
    tReal mLmax;
    tReal mL;
    
    tUIndex mBoundElement[FACES_NUMBER_PER_ELEMENT];
    
    tIndex mNextElementTranslation[FACES_NUMBER_PER_ELEMENT];
    
    tIndex mPeriodicElementTranslation[FACES_NUMBER_PER_ELEMENT];


    SP::EMM_CellFlagArray mNeighborsIndicators;

    tUIndex mMagnetizedElementsNumber;
    
    tUIndex mElementsNumber;


    //temporary index
    mutable tUIndex mWorkIndex; 

protected: 
    // METHODS
  
    // CONSTRUCTORS 
  
    EMM_Grid3D(void);
 
    // DESTRUCTORS 
  
    
    virtual ~EMM_Grid3D(void);   
  
  
public:
    inline static SP::EMM_Grid3D New() {
        SP::EMM_Grid3D p(new EMM_Grid3D(),EMM_Grid3D::Delete());
        p->setThis(p);
        return p;
    };
    
    
    
    //builder methods
    //===============
private:
    void updateElementTranslations();

public:
    void updateMagnetizedElementsNumber();

    
public:

    //Segments number manipulation
    //=============================
    
    inline void setSegmentsNumber(const tUInteger& nx,const tUInteger& ny,const tUInteger& nz)  { 
        mSegmentsNumber[0]=nx;
        mSegmentsNumber[1]=ny;
        mSegmentsNumber[2]=nz;
        mElementsNumber=nx*ny*nz;
        mNeighborsIndicators->setSize(mElementsNumber);
        updateElementTranslations();
    };
    
    inline void getSegmentsNumber(tUInteger& nx,tUInteger& ny,tUInteger& nz) const { 
        nx=mSegmentsNumber[0];
        ny=mSegmentsNumber[1];
        nz=mSegmentsNumber[2];
        
    }
    inline const tUInteger* getSegmentsNumber() const { 
        return mSegmentsNumber;
    }
    
    inline const tUInteger& getSegmentsNumber(const tUSInt& d) const { 
        return mSegmentsNumber[d];
    }
    
    
    inline tUIndex getElementIndex(const tUIndex& i,const tUIndex& j,const tUIndex& k) const {
        return i+mSegmentsNumber[0]*(j+mSegmentsNumber[1]*k);
    }
    inline void getElementSegmentIndices(const tUIndex& index,tUIndex& i,tUIndex& j,tUIndex& k) const {
        i=index%mSegmentsNumber[0];
        mWorkIndex=index/mSegmentsNumber[0];//mWorkIndex=j+k.Ny
        j=mWorkIndex%mSegmentsNumber[1];
        k=mWorkIndex/mSegmentsNumber[1];
    }
    inline tUIndex  getElementSegmentIndex(const tUIndex& index,const tUSInt& d) const {
        switch(d) {
        case 0:
            return index%mSegmentsNumber[0];
        case 1:
            return (index/mSegmentsNumber[0])%mSegmentsNumber[1];
        case 2:
        default:
            return (index/mSegmentsNumber[0])/mSegmentsNumber[1];
        }
    }


    inline tUIndex getElementsNumber()const{
        return mSegmentsNumber[0]*mSegmentsNumber[1]*mSegmentsNumber[2];
    };
    
    inline  tUIndex getVerticesNumber()const{
        return (mSegmentsNumber[0]+1)*(mSegmentsNumber[1]+1)*(mSegmentsNumber[2]+1);
    };

    //step size methods
    //=================

public:
    inline void setStepsSize(const tReal& lx,const tReal& ly,const tReal& lz) {
        mSegmentsSize[0]=lx;
        mSegmentsSize[1]=ly;
        mSegmentsSize[2]=lz;
    };
    
    
    inline void setStepSize(const tReal& l) {
        setStepsSize(l,l,l);
    }; 
    inline const tReal* getStepsSize() const {
        return mSegmentsSize;
    }  
    inline const tReal& getStepSize(const tUSInt& d) const {
        return mSegmentsSize[d];
    }
    inline tReal getMinStepSize() const {
        tReal dx=mSegmentsSize[0];
        for (tUSInt i=1;i<3;i++) dx=(dx<mSegmentsSize[i])?dx:mSegmentsSize[i];
        return dx;
    }
    inline tReal getMaxStepSize() const {
        tReal dx=mSegmentsSize[0];
        for (tUSInt i=1;i<3;i++) dx=(dx>mSegmentsSize[i])?dx:mSegmentsSize[i];
        return dx;
    }
    inline const tReal* getAdimensionizedSegmentsSize() const {
        return mAdimensionizedSegmentsSize;
    }

    inline const tReal & getReferenceLength() const {
        return mL;
    }
    inline const tReal & getMaxElementSize() const {
        return mLmax;
    }
    inline const tReal & getMinElementSize() const {
        return mLmin;
    }
    inline const tReal* getSegmentsSize() const {
        return mSegmentsSize;
    }

    inline tBoolean isCubicVolume() const {
        const tReal *lk=mAdimensionizedSegmentsSize;
        const tReal *l =&mAdimensionizedSegmentsSize[2];
        tBoolean ok=true;
        for (tUSInt k=0;k<2;k++) {
            ok=ok && (fabs(1.-(*lk)/(*l))<CORE_Object::computeEpsilon<tReal>());
            lk++;
        }
        return ok;
    }
    
   
    inline  tReal adimensionize() {
        tUSInt k=0;
        //get the max size of the element
        mLmax=mSegmentsSize[0];
        //get the min size of the element
        mLmin=mSegmentsSize[0];
        for (k=1;k<3;k++) {
            if (mSegmentsSize[k]>mLmax) {
                mLmax=mSegmentsSize[k];
            }
            if (mSegmentsSize[k]<mLmin) {
                mLmin=mSegmentsSize[k];
            }
        }
        
        //reference length is the max cell size
        mL=mLmax;
        //normalize the element
        mAdimensionizedVolume=1.;
        for (k=0;k<3;k++) {
            mAdimensionizedSegmentsSize[k]=mSegmentsSize[k]/mL;
            mAdimensionizedVolume*=mAdimensionizedSegmentsSize[k];
        }
        return mL;
        
    }

    inline const tReal& getAdimensionizedVolume() const {
        return mAdimensionizedVolume;
    }
    
    void computeNormalFaceAdimensionizedLength(tReal* L2) const;
    
    
  
    
    // PERIODICITY
    //=============
   
    inline void setIsPeriodic(const tBoolean v[3]) {
        mIsPeriodic=false;
        for (tUSInt i=0;i<3;i++) {
            mIsPeriodic=mIsPeriodic || v[i];
            mIsDirectionPeriodic[i]=v[i];
        }
    }
    inline void setIsDirectionPeriodic(const tUSInt& k,const tBoolean& v) {
        mIsDirectionPeriodic[k]=v;
        mIsPeriodic=false;
        for (tUSInt i=0;i<3;i++) {
            mIsPeriodic=mIsPeriodic || mIsDirectionPeriodic[i];
        }
    }
    
    inline const tBoolean* isDirectionPeriodic() const {
        return mIsDirectionPeriodic;
    }
    inline const tBoolean& isDirectionPeriodic(const tFlag& d) const {
        ASSERT_IN(d<3);
        return mIsDirectionPeriodic[d];
    }
    inline const tBoolean& isPeriodic() const {
        return mIsPeriodic;
    }

    inline tCellFlag getPeriodicIndicator() const {
        return getPeriodicIndicator(mIsDirectionPeriodic);
    };
    inline static tCellFlag getPeriodicIndicator(const tBoolean periodicity[]) {
        return 3*periodicity[0]+4*periodicity[1]+16*periodicity[2];
    };
    

  
  
    // cell methods
    //=============

    inline const tUIndex& getMagnetizedElementsNumber() const {
        return mMagnetizedElementsNumber;
    }
    
    inline tBoolean isElementMagnetized(const tUIndex& index)const{
        return  ((*mNeighborsIndicators.get())[index]>=MAGNETIZED_ELEMENT_MIN_INDICATOR);
    };
      
    static inline tBoolean IS_ELEMENT_MAGNETIZED(const tCellFlag& v) {
        return ((64&v)!=0);
    }
    
    inline tBoolean isInMagnetizedDomainInterior(const tUIndex& index)const{
        return  ((*mNeighborsIndicators.get())[index]==MAGNETIZED_DOMAIN_INTERIOR_INDICATOR);
    };

    void buildMagnetizedElements(const EMM_MatterField &matters,
                                 CORE_BooleanMorseArray& magnetizedElements) const;

    
    //neighbor methods
    //=================
    
    inline tBoolean hasNeighborToFace(const tUIndex& iElement,
                                      const tUSInt& iFace) const {
        if (iFace>=FACES_NUMBER_PER_ELEMENT)
            throw EMM_Exception("EMM_Grid3D","hasNeighborToFace",
                                "index of local face is greater than the max faces number in elements");
        tCellFlag p=TWO_POWER_FACES[iFace];
        return (iand((*mNeighborsIndicators.get())[iElement],p)==p);
    };

   
    
    tUIndex getNeighborElement(const tUIndex& index,const tUCInt& iFace) const;

    tUIndex getNeighborElement(const tUIndex& iCell,
                               const tCellFlag& N,
                               const tBoolean isPeriodic[3],
                               const tUCInt &f) const;
    

   
    inline SPC::EMM_CellFlagArray getNeighborsIndicatorsByReference() const {
        return mNeighborsIndicators;
    }
    
    inline const EMM_CellFlagArray& getNeighborsIndicators() const {
        return *mNeighborsIndicators.get();
    }
    
    inline EMM_CellFlagArray& getNeighborsIndicators()  {
        return *mNeighborsIndicators.get();
    }

 
    inline const tCellFlag& getNeighborIndicator(const tUIndex& index) const {
        return (*mNeighborsIndicators.get())[index];
    }

    void buildNeighborsIndices(const EMM_MatterField &matters,
                               CORE_UIndexMorseArray& neighbors) const;

    void buildMagnetizedNeighborsIndices(const EMM_MatterField &matters,
                                         CORE_UIndexMorseArray& neighbors) const;
    
    
    //boundary faces methods
    //======================

    
    static tBoolean isFaceOnDirichletBoundary(const tCellFlag& periodicity,
                                              const tUInteger& xCell,
                                              const tUInteger& yCell,
                                              const tUInteger& zCell,
                                              const tUCInt& f,
                                              const tUInteger& nPx,
                                              const tUInteger& nPy,
                                              const tUInteger& nPz,
                                              const tLimitCondition* dirichletPoints);
    
    static tBoolean isFaceOnNeumannBoundary(const tCellFlag& periodicity,
                                            const tUInteger& xCell,
                                            const tUInteger& yCell,
                                            const tUInteger& zCell,
                                            const tUCInt& f,
                                            const tUInteger& nPx,
                                            const tUInteger& nPy,
                                            const tUInteger& nPz,
                                            const tLimitCondition* dirichletPoints);

    
    static inline tBoolean GET_MASTER_PERIODIC_POINT(const tCellFlag& periodicity,
                                                     const tUIndex& Nx,const tUIndex& Ny,const tUIndex& Nz,
                                                     tUIndex& iQ,tUIndex& jQ,tUIndex& kQ) {
        tBoolean moved=false;
        if ((periodicity&3)!=0) {
            //x-direction is periodic
            if (iQ>=Nx) {
                iQ-=Nx;
                moved=true;
            }
           
            
        }
        if ((periodicity&12)!=0) {
            //y-direction is periodic
            if (jQ>=Ny) {
                jQ-=Ny;
                moved=true;
            }
           
            
        }
        if ((periodicity&48)!=0) {
            //z-direction is periodic
            if (kQ>=Nz) {
                kQ-=Nz;
                moved=true;
            } 
            
        }
        return moved;
    }

    
    
    static inline tBoolean IS_POINT_IN_BOUNDARY(const tBoolean x[],
                                                const tCellFlag& N,
                                                const tBoolean periodicity[]) {
        return IS_POINT_IN_BOUNDARY(21+(*x)+4*x[1]+16*x[2],
                                    N,
                                    3*((*periodicity)+4*periodicity[1]+16*periodicity[2]));
    }
    
    static inline tBoolean IS_POINT_IN_BOUNDARY(const tCellFlag& X,
                                                const tCellFlag& N,
                                                const tCellFlag& P) {
        return ((X&((~N)&(~P)))!=0);
    }
    
    static inline tBoolean IS_POINT_IN_PERIODIC_BOUNDARY(const tBoolean x[],
                                                         const tCellFlag& N,
                                                         const tBoolean periodicity[]) {
        //Pf=X0.2 +(1-X0).1+
        //   X1.8 +(1-x1).4+
        //   X2.32+(1-x2).16
        //  = 21+x0+4.x1+16.X2
        return IS_POINT_IN_PERIODIC_BOUNDARY(21+(*x)+4*x[1]+16*x[2],
                                             N,
                                             3*((*periodicity)+4*(periodicity[1]+4*periodicity[2])));
    }
    static inline tBoolean IS_POINT_IN_PERIODIC_BOUNDARY(const tCellFlag& X,
                                                         const tCellFlag& N,
                                                         const tCellFlag& P) {
        return ((X&((~N)&(P)))!=0);
    }
    static inline tBoolean IS_POINT_IN_SLAVE_PERIODIC_BOUNDARY(const tBoolean x[],
                                                               const tCellFlag& N,
                                                               const tBoolean periodicity[]) {
        
        return IS_POINT_IN_SLAVE_PERIODIC_BOUNDARY(21+(*x)+4*x[1]+16*x[2],
                                                   N,
                                                   3*((*periodicity)+4*(periodicity[1]+4*periodicity[2])));
        
    }

    static inline tBoolean IS_POINT_IN_SLAVE_PERIODIC_BOUNDARY(const tCellFlag& X,
                                                               const tCellFlag& N,
                                                               const tCellFlag& P) {
        return ((X&((~N)&(P))&(42))!=0);
    }
    
    //limit conditions
    //=================

    
    void updateLimitConditionOnPoints(EMM_LimitConditionArray& lc) const;
    

 
  
    //  GEO File saver & loader
    // ========================
public:
    tBoolean loadFromFile(const tString& fileName);

    tBoolean saveToFile(const tString& fileName) const;
  
                      
   
private: 
    tBoolean loadFromGEOFile(const tString& fileName);

    tBoolean saveToGEOFile(const tString& fileName) const;
    
    tBoolean loadHeaderFromGEOStream(ifstream& f,tUInteger& nReadLines,tUInteger& nReadProperties);
    tBoolean loadFromGEOStream(ifstream& f,tUInteger& nReadLines);
    
    tBoolean loadFooterFromGEOStream(ifstream& f,tUInteger& nReadLines,tUInteger& nReadProperties);
    

    //STRING REPRESENTATION
public:
    virtual tString toString() const;
  
};   
    

#endif