EMM_MatterField.hpp

Go to the documentation of this file.

#ifndef EMM_MatterField_HPP
#define EMM_MatterField_HPP

#include "OMP_Thread.h"

template<class T>
tBoolean EMM_MatterField::getMatterParameterDistribution(const tFlag& k,EMM_Array<T>& vs) const {
    
   
   
   
    //get the indices of matters at each cell if any
    const EMM_UIndexArray &matterIndices=mMatterIndices;
    const SV::EMM_Matter &matters=mMatters;

    //size of the indices
    tUIndex iCell,nCells=matterIndices.getSize();
    if ((nCells==0) || (nCells==1)) {
         //get the value
        const tReal* parameters=&(matters[0]->getParameter(k));
        const tUSInt& dim=matters[0]->getParameterDimension(k);
        
        //uniform
        vs.setSize(dim);
        T *v=&vs[0];
       
        tUSInt d;
        if (mMatters.getSize()>0) {
            //there is at least one matter
            for (d=0;d<dim;d++) {
                *v=(T) (*parameters);
                v++;
                parameters++;
            }
        } else {
            //there is no matter
            memset(v,0,dim*sizeof(T));
        }
        return true;
    }

    const tUSInt& dim=matters[0]->getParameterDimension(k);
    tUSInt d;
    
    //set the size of vs
    vs.setSize(nCells*dim);

    T *v;
        
    OMP_PARALLEL_PRIVATE_SHARED(private(d,iCell,v),shared(cout,dim,nCells,vs,matterIndices,matters,k)) {
            
        
        // thread id
        tUInteger threadId=OMP_GET_THREAD_ID();
        
        // threads number
        tUInteger nThreads=1;
        nThreads=OMP_GET_THREADS_NUMBER();
        
        //start index
        tUIndex start=threadId*nCells/nThreads;
        //end index
        tUIndex  end=(threadId+1)*nCells/nThreads;
        
        const tUIndex *matterIndex=&matterIndices[start];
        
        tUIndex nMatters=matters.getSize();
        
        const tReal* parameters;

        //iterator on value of the parameter k
        v=&vs[start*dim];
      
        for (iCell=start;iCell<end;iCell++) {//loop on cell 
            if ((*matterIndex)==EMM_MatterField::NO_MATTER) {
                //no matter at iCell
                memset(v,0,dim*sizeof(T));
                v+=dim;
            }  else {
                //get the matter at iCell
                if (*matterIndex<nMatters) {
                    //set the value
                    parameters=&matters[*matterIndex]->getParameter(k);
                    for (d=0;d<dim;d++) {
                        *v=(T) (*parameters);
                        v++;
                        parameters++;
                    }
                } else {//no matter so empty
                    memset(v,0,dim*sizeof(T));
                    v+=dim;
                }
            } 
            //next value pointer
            matterIndex++;
        }//end loop on cell
    } //end OMP parallel loop

    
    //test if the field is uniform
    T vs_old[dim];
    T *vold;
    vold=vs_old;
    for (d=0;d<dim;d++) {
        *vold=(T) 0;
        vold++;
    }
    tBoolean isInitialized=false;
    
    tBoolean isUniform=true;
    iCell=0;

    v=&vs[0];
    tReal eps=1.e-12;

    do {
        if (!isInitialized) {
            //initialize not null uniform value
            vold=vs_old;
            for (d=0;d<dim;d++) {
                (*vold)=*v;
                if ((*v)!=0) isInitialized=true;

                vold++;
                v++;
            }
        } else {
            vold=vs_old;
            for (d=0;d<dim;d++) {
                if (fabs(1-(*v)/(*vold))>eps) {
                    //v!=v_old
                    isUniform=false;
                }
                v++;
                vold++;
            }
        }
        //next iter
        iCell++;
    } while ((iCell<nCells) && (isUniform));
    
    if (isUniform) {
        //resize the array to size dim
        vs.setSize(dim);
        v=&vs[0];
        vold=vs_old;
        for (d=0;d<dim;d++) {
            (*v)=(*vold);
            v++;
            vold++;
        }
        //set capacity to size
        vs.fitToSize();
    }
    return isUniform;

}

template<class T>
tBoolean EMM_MatterField::getMatterParameterDistribution(const tFlag& k,CORE_MorseArray<T>& vs) const {
    
   
   
    
    //get the indices of matters at each cell if any
    const EMM_UIndexArray &matterIndices=mMatterIndices;
    const SV::EMM_Matter &matters=mMatters;

    //size of the indices
    tUIndex iCell,nCells=matterIndices.getSize();
    if ((nCells==0) || (nCells==1)) {
        //mater is uniform on all cells
        //dimension of the parameter
        const tUSInt& dim=matters[0]->getParameterDimension(k);
        tUSInt d;
        //the morse array is uniform
        vs.setUniformSize(1,dim);
        //copy the values of the parameter of size dim
        T *v=vs[0];
        const tReal* parameters=&(matters[0]->getParameter(k));
      
        if (mMatters.getSize()>0) {
            //there is at least one matter
            for (d=0;d<dim;d++) {
                *v=(T) (*parameters);
                parameters++;
                v++;
            }
        } else {
            //no matter
            memset(v,0,dim*sizeof(T));
        }
        //the parameter is uniform
        return true;
    }
    
    
    //get the number of threads
    tUInteger threadId,nThreads=1;
    OMP_PARALLEL_SHARED(shared(nThreads)) {
        if( OMP_GET_THREAD_ID()==0)
            nThreads=OMP_GET_THREADS_NUMBER();
    }

      
    //dimension of the parameter
    const tUSInt& dim=matters[0]->getParameterDimension(k);
    tUSInt d;
    
    //set the morse array size
    vs.setSize(nCells,dim,nThreads);
    vs.reset();

    vector<CORE_MorseArrayIterator<T> > iterVs(nThreads);
    
    OMP_PARALLEL_PRIVATE_SHARED(private(iCell,d,threadId),shared(nThreads,cout,dim,nCells,vs,iterVs,matterIndices,matters,k)) {
            
  
        
        // thread id
        threadId=OMP_GET_THREAD_ID();
        
        
        //start index
        tUIndex start=threadId*nCells/nThreads;
        //end index
        tUIndex  end=(threadId+1)*nCells/nThreads;
        
        const tUIndex *matterIndex=&matterIndices[start];
        
        tUIndex nMatters=matters.getSize();

        //matter of the cell i
        const EMM_Matter *matter=null;
        const tReal* parameters=null;
        
        //get the iterator of the morse array for the thread threadId
        CORE_MorseArrayIterator<T> &iterV=iterVs[threadId];
        T *v=null;
        vs.begin(start,threadId,iterV);
      

        //memory size of the value for the cell (1 or 0)
        tUIndex sV;
        
        for (iCell=start;iCell<end;iCell++) {//loop on cell
            //memory size 0 by default
            sV=0;

            //register the parameter in case of magnetized cell
            if ( (*matterIndex)<nMatters) {
                matter=matters.get(*matterIndex);
                if (matter!=null) {
                    //not null matter
                    sV=dim;
                    v=iterV.values();
                    parameters=&matter->getParameter(k);
                    for (d=0;d<dim;d++) {
                        (*v)=(T) (*parameters);
                        v++;
                        parameters++;
                    }
                    
                } 
            } 
            //next value pointer
            iterV.setSize(sV);
            ++iterV;
            //matter of next cell
            matterIndex++;
        }//end loop on cell
    } //end OMP parallel loop

    
    //merge the packs
    vs.merge();

    //uniformize the morse array
    return vs.uniformize();
    
}
#endif