Ignore:
Timestamp:
01/22/21 12:00:29 (3 years ago)
Author:
yushan
Message:

Graph intermedia commit to a tmp branch

File:
1 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_trunk_graph/src/transformation/generic_algorithm_transformation.cpp

    r1612 r2019  
    1414#include "timer.hpp" 
    1515#include "mpi.hpp" 
     16#include "transform_connector.hpp" 
     17#include "weight_transform_connector.hpp" 
     18#include "grid_algorithm_generic.hpp" 
     19#include "transform_filter.hpp" 
    1620 
    17 namespace xios { 
     21namespace xios  
     22{ 
    1823 
    19 CGenericAlgorithmTransformation::CGenericAlgorithmTransformation() 
    20  : transformationMapping_(), transformationWeight_(), transformationPosition_(), 
    21    idAuxInputs_(), type_(ELEMENT_NO_MODIFICATION_WITH_DATA), indexElementSrc_(), 
    22    computedProcSrcNonTransformedElement_(false), eliminateRedondantSrc_(true), isDistributedComputed_(false) 
     24CGenericAlgorithmTransformation::CGenericAlgorithmTransformation(bool isSource) 
     25 : isSource_(isSource) 
    2326{ 
    2427} 
    2528 
    26 void CGenericAlgorithmTransformation::updateData(CArray<double,1>& dataOut) 
     29 
     30 
     31/////////////////////////////////////////////////////////////// 
     32////////// new algorithm for new method               ///////// 
     33/////////////////////////////////////////////////////////////// 
     34 
     35CGridAlgorithm* CGenericAlgorithmTransformation::createGridAlgorithm(CGrid* gridSrc, CGrid* gridDst, int pos) 
    2736{ 
    28  
     37  return new CGridAlgorithmGeneric(gridSrc, gridDst, pos, this) ; 
    2938} 
    30  
    31 void CGenericAlgorithmTransformation::apply(const std::vector<std::pair<int,double> >& localIndex, 
    32                                             const double* dataInput, 
    33                                             CArray<double,1>& dataOut, 
    34                                             std::vector<bool>& flagInitial, 
    35                                             bool ignoreMissingValue, bool firstPass  ) 
    36 TRY 
    37 { 
    38   int nbLocalIndex = localIndex.size();    
    39   double defaultValue = std::numeric_limits<double>::quiet_NaN(); 
    40      
    41   if (ignoreMissingValue) 
    42   { 
    43     if (firstPass) dataOut=defaultValue ; 
    44      
    45     for (int idx = 0; idx < nbLocalIndex; ++idx) 
    46     { 
    47       if (! NumTraits<double>::isNan(*(dataInput + idx))) 
    48       { 
    49         if (flagInitial[localIndex[idx].first]) dataOut(localIndex[idx].first) = *(dataInput + idx) * localIndex[idx].second; 
    50         else dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second; 
    51         flagInitial[localIndex[idx].first] = false; // Reset flag to indicate not all data source are nan 
    52       } 
    53     } 
    54  
    55   } 
    56   else 
    57   { 
    58     for (int idx = 0; idx < nbLocalIndex; ++idx) 
    59     { 
    60       dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second; 
    61     } 
    62   } 
    63 } 
    64 CATCH 
    65  
    66 void CGenericAlgorithmTransformation::computePositionElements(CGrid* dst, CGrid* src) 
    67 TRY 
    68 { 
    69   int idxScalar = 0, idxAxis = 0, idxDomain = 0; 
    70   CArray<int,1> axisDomainOrderDst = dst->axis_domain_order; 
    71   for (int i = 0; i < axisDomainOrderDst.numElements(); ++i) 
    72   { 
    73     int dimElement = axisDomainOrderDst(i); 
    74     if (2 == dimElement) 
    75     { 
    76       elementPositionInGridDst2DomainPosition_[i] = idxDomain; 
    77       ++idxDomain; 
    78     } 
    79     else if (1 == dimElement) 
    80     { 
    81       elementPositionInGridDst2AxisPosition_[i] = idxAxis; 
    82       ++idxAxis; 
    83     } 
    84     else 
    85     { 
    86       elementPositionInGridDst2ScalarPosition_[i] = idxScalar; 
    87       ++idxScalar; 
    88     } 
    89   } 
    90  
    91   idxScalar = idxAxis = idxDomain = 0; 
    92   CArray<int,1> axisDomainOrderSrc = src->axis_domain_order; 
    93   for (int i = 0; i < axisDomainOrderSrc.numElements(); ++i) 
    94   { 
    95     int dimElement = axisDomainOrderSrc(i); 
    96     if (2 == dimElement) 
    97     { 
    98       elementPositionInGridSrc2DomainPosition_[i] = idxDomain; 
    99       ++idxDomain; 
    100     } 
    101     else if (1 == dimElement) 
    102     { 
    103       elementPositionInGridSrc2AxisPosition_[i] = idxAxis; 
    104       ++idxAxis; 
    105     } 
    106     else 
    107     { 
    108       elementPositionInGridSrc2ScalarPosition_[i] = idxScalar; 
    109       ++idxScalar; 
    110     } 
    111   } 
    112 } 
    113 CATCH 
    114  
    115 bool CGenericAlgorithmTransformation::isDistributedTransformation(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst) 
    116 TRY 
    117 { 
    118  
    119   if (!isDistributedComputed_) 
    120   { 
    121     isDistributedComputed_=true ; 
    122     if (!eliminateRedondantSrc_) isDistributed_=true ; 
    123     else 
    124     { 
    125       CContext* context = CContext::getCurrent(); 
    126       CContextClient* client = context->client; 
    127    
    128       computePositionElements(gridSrc, gridDst); 
    129       std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars(); 
    130       std::vector<CAxis*> axisListSrcP = gridSrc->getAxis(); 
    131       std::vector<CDomain*> domainListSrcP = gridSrc->getDomains(); 
    132       int distributed, distributed_glo ; 
    133    
    134       CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order; 
    135       if (2 == axisDomainSrcOrder(elementPositionInGrid)) // It's domain 
    136       { 
    137         distributed=domainListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]]->isDistributed() ; 
    138         MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, client->intraComm) ; 
    139      
    140       } 
    141       else if (1 == axisDomainSrcOrder(elementPositionInGrid))//it's an axis 
    142       { 
    143         distributed=axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]]->isDistributed() ; 
    144         MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, client->intraComm) ; 
    145       } 
    146       else //it's a scalar 
    147       { 
    148         distributed_glo=false ; 
    149       }  
    150       isDistributed_=distributed_glo ; 
    151     } 
    152   } 
    153   return isDistributed_ ; 
    154 } 
    155 CATCH 
    156  
    157 /*! 
    158   This function computes the global indexes of grid source, which the grid destination is in demand. 
    159   \param[in] elementPositionInGrid position of an element in a grid .E.g: if grid is composed of domain and axis (in order), 
    160                 then position of axis in grid is 1 and domain is positioned at 0. 
    161   \param[in] gridSrc Grid source 
    162   \param[in] gridDst Grid destination 
    163   \param[in\out] globaIndexWeightFromSrcToDst mapping of each global index source and weight to index destination 
    164 */ 
    165 void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid, 
    166                                                                CGrid* gridSrc, 
    167                                                                CGrid* gridDst, 
    168                                                                SourceDestinationIndexMap& globaIndexWeightFromSrcToDst) 
    169 TRY 
    170  { 
    171   CContext* context = CContext::getCurrent(); 
    172   CContextClient* client = context->client; 
    173   int nbClient = client->clientSize; 
    174  
    175   typedef std::unordered_map<int, std::vector<std::pair<int,double> > > SrcToDstMap; 
    176   int idx; 
    177  
    178   // compute position of elements on grids 
    179   computePositionElements(gridDst, gridSrc); 
    180   std::vector<CScalar*> scalarListDestP = gridDst->getScalars(); 
    181   std::vector<CAxis*> axisListDestP = gridDst->getAxis(); 
    182   std::vector<CDomain*> domainListDestP = gridDst->getDomains(); 
    183   CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order; 
    184   std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars(); 
    185   std::vector<CAxis*> axisListSrcP = gridSrc->getAxis(); 
    186   std::vector<CDomain*> domainListSrcP = gridSrc->getDomains(); 
    187   CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;   
    188    
    189   bool isTransPosEmpty = transformationPosition_.empty(); 
    190   CArray<size_t,1> transPos; 
    191   if (!isTransPosEmpty) transPos.resize(transformationMapping_.size()); 
    192   std::set<size_t> allIndexSrc; // All index of source, which can be scattered among processes, need for doing transformation  
    193    
    194   for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans) 
    195   { 
    196     TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap, 
    197                                            iteTransMap = transformationMapping_[idxTrans].end(); 
    198     TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight; 
    199  
    200     // Build mapping between global source element index and global destination element index. 
    201     itTransWeight = itbTransWeight; 
    202     for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight) 
    203     { 
    204       const std::vector<int>& srcIndex = itTransMap->second; 
    205       for (idx = 0; idx < srcIndex.size(); ++idx) 
    206         allIndexSrc.insert(srcIndex[idx]); 
    207     } 
    208  
    209     if (!isTransPosEmpty) 
    210     { 
    211       TransformationPositionMap::const_iterator itPosMap = transformationPosition_[idxTrans].begin(); 
    212       transPos(idxTrans) = itPosMap->second[0]; 
    213     } 
    214   } 
    215  
    216   size_t indexSrcSize = 0; 
    217   CArray<size_t,1> indexSrc(allIndexSrc.size()); 
    218   for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it, ++indexSrcSize) 
    219     indexSrc(indexSrcSize) = *it; 
    220  
    221   // Flag to indicate whether we will recompute distribution of source global index  on processes 
    222   bool computeGlobalIndexOnProc = false;  
    223   if (indexElementSrc_.size() != allIndexSrc.size()) 
    224     computeGlobalIndexOnProc = true;   
    225   else 
    226   { 
    227     for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it) 
    228       if (0 == indexElementSrc_.count(*it)) 
    229       { 
    230         computeGlobalIndexOnProc = true; 
    231         break;         
    232       } 
    233   } 
    234  
    235   if (computeGlobalIndexOnProc) 
    236     indexElementSrc_.swap(allIndexSrc); 
    237        
    238   int sendValue = (computeGlobalIndexOnProc) ? 1 : 0; 
    239   int recvValue = 0; 
    240   MPI_Allreduce(&sendValue, &recvValue, 1, MPI_INT, MPI_SUM, client->intraComm); 
    241   computeGlobalIndexOnProc = (0 < recvValue); 
    242  
    243 //  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc; 
    244  
    245   if (computeGlobalIndexOnProc || !computedProcSrcNonTransformedElement_) 
    246   {     
    247     { 
    248       CClientClientDHTInt::Index2VectorInfoTypeMap tmp ; 
    249       globalIndexOfTransformedElementOnProc_.swap(tmp) ; 
    250     } 
    251     // Find out global index source of transformed element on corresponding process.     
    252     if (globalElementIndexOnProc_.empty()) 
    253       globalElementIndexOnProc_.resize(axisDomainDstOrder.numElements()); 
    254      
    255     for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx) 
    256     { 
    257        
    258       if (idx == elementPositionInGrid) 
    259         computeExchangeGlobalIndex(indexSrc, axisDomainSrcOrder(idx), globalIndexOfTransformedElementOnProc_); //globalElementIndexOnProc[idx]); 
    260       if (!computedProcSrcNonTransformedElement_) 
    261       { 
    262         if (2 == axisDomainDstOrder(idx)) // It's domain 
    263         { 
    264           if (idx != elementPositionInGrid) 
    265             computeExchangeDomainIndex(domainListDestP[elementPositionInGridDst2DomainPosition_[idx]], 
    266                                        domainListSrcP[elementPositionInGridSrc2DomainPosition_[idx]], 
    267                                        transPos, 
    268                                        globalElementIndexOnProc_[idx]);       
    269  
    270         } 
    271         else if (1 == axisDomainDstOrder(idx))//it's an axis 
    272         { 
    273           if (idx != elementPositionInGrid) 
    274             computeExchangeAxisIndex(axisListDestP[elementPositionInGridDst2AxisPosition_[idx]], 
    275                                      axisListSrcP[elementPositionInGridSrc2AxisPosition_[idx]], 
    276                                      transPos, 
    277                                      globalElementIndexOnProc_[idx]); 
    278         } 
    279         else //it's a scalar 
    280         { 
    281           if (idx != elementPositionInGrid) 
    282             computeExchangeScalarIndex(scalarListDestP[elementPositionInGridDst2ScalarPosition_[idx]], 
    283                                        scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[idx]], 
    284                                        transPos, 
    285                                        globalElementIndexOnProc_[idx]); 
    286  
    287         } 
    288       } 
    289     } 
    290  
    291     if (!isTransPosEmpty && !computedProcSrcNonTransformedElement_) 
    292     { 
    293       for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx) 
    294       { 
    295         if (idx != elementPositionInGrid) 
    296         { 
    297           std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it, 
    298                                                                    ite = globalElementIndexOnProc_[idx].end(); 
    299           for (it = itb; it != ite; ++it) it->second.resize(1); 
    300         } 
    301       } 
    302     } 
    303  
    304 /*       
    305     if (!computedProcSrcNonTransformedElement_) 
    306     { 
    307       for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx) 
    308       { 
    309         if (idx != elementPositionInGrid) 
    310         { 
    311           std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it, 
    312                                                                    ite = globalElementIndexOnProc_[idx].end(); 
    313           for (it = itb; it != ite; ++it) procOfNonTransformedElements_.insert(it->first); 
    314           if (procOfNonTransformedElements_.size() == nbClient) 
    315             break; 
    316         } 
    317       } 
    318     } 
    319     
    320     // Processes contain the source index of transformed element 
    321     std::set<int> procOfTransformedElement;  
    322     CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc.begin(), 
    323                                                            itIdxe = globalIndexOfTransformedElementOnProc.end(), itIdx; 
    324     for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx) 
    325     { 
    326       std::vector<int>& tmp = itIdx->second; 
    327       for (int i = 0; i < tmp.size(); ++i) 
    328         procOfTransformedElement.insert(tmp[i]); 
    329       if (tmp.size() == nbClient) 
    330         break; 
    331     }                                                            
    332      
    333     std::set<int>& commonProc = (procOfTransformedElement.size() < procOfNonTransformedElements_.size()) ? procOfTransformedElement  
    334                               : (!procOfNonTransformedElements_.empty() ? procOfNonTransformedElements_ : procOfTransformedElement); 
    335      
    336     std::vector<int> procContainSrcElementIdx(commonProc.size()); 
    337     int count = 0; 
    338     for (std::set<int>::iterator it = commonProc.begin(); it != commonProc.end(); ++it) 
    339       procContainSrcElementIdx[count++] = *it; 
    340  
    341     procContainSrcElementIdx_.swap(procContainSrcElementIdx);     
    342 */ 
    343      
    344     if (procElementList_.empty()) procElementList_.resize(axisDomainDstOrder.numElements()) ; 
    345     for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx) 
    346     { 
    347       std::set<int>& procList=procElementList_[idx] ; 
    348       std::set<int> commonTmp ; 
    349       if (idx == elementPositionInGrid) 
    350       { 
    351           set<int> tmpSet ;  
    352           procList.swap(tmpSet) ; 
    353           CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc_.begin(), 
    354                                                                  itIdxe = globalIndexOfTransformedElementOnProc_.end(), itIdx; 
    355           for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx) 
    356           { 
    357              std::vector<int>& tmp = itIdx->second; 
    358              for (int i = 0; i < tmp.size(); ++i) procList.insert(tmp[i]); 
    359              if (tmp.size() == nbClient) 
    360              break; 
    361           } 
    362       } 
    363       else 
    364       { 
    365         if (!computedProcSrcNonTransformedElement_) 
    366         { 
    367           set<int> tmpSet ;  
    368           procList.swap(tmpSet) ; 
    369           std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it, 
    370                                                                    ite = globalElementIndexOnProc_[idx].end(); 
    371           for (it = itb; it != ite; ++it) 
    372           { 
    373             procList.insert(it->first); 
    374             if (procList.size() == nbClient)  break; 
    375           } 
    376         } 
    377       } 
    378  
    379       if (idx==0) commonProc_= procList ; 
    380       else 
    381       { 
    382         for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it) 
    383           if (procList.count(*it)==1) commonTmp.insert(*it) ; 
    384         commonProc_.swap(commonTmp) ; 
    385       } 
    386     } 
    387     std::vector<int> procContainSrcElementIdx(commonProc_.size()); 
    388     int count = 0; 
    389     for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it) procContainSrcElementIdx[count++] = *it; 
    390     procContainSrcElementIdx_.swap(procContainSrcElementIdx); 
    391      
    392         // For the first time, surely we calculate proc containing non transformed elements 
    393     if (!computedProcSrcNonTransformedElement_) computedProcSrcNonTransformedElement_ = true; 
    394   } 
    395    
    396   for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans) 
    397   { 
    398     TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap, 
    399                                            iteTransMap = transformationMapping_[idxTrans].end(); 
    400     TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight; 
    401     SrcToDstMap src2DstMap; 
    402     src2DstMap.rehash(std::ceil(transformationMapping_[idxTrans].size()/src2DstMap.max_load_factor())); 
    403  
    404     // Build mapping between global source element index and global destination element index. 
    405     std::unordered_map<int,std::vector<size_t> >().swap(globalElementIndexOnProc_[elementPositionInGrid]); 
    406     std::set<int> tmpCounter; 
    407     itTransWeight = itbTransWeight; 
    408     for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight) 
    409     { 
    410       const std::vector<int>& srcIndex = itTransMap->second; 
    411       const std::vector<double>& weight = itTransWeight->second; 
    412       for (idx = 0; idx < srcIndex.size(); ++idx) 
    413       { 
    414         src2DstMap[srcIndex[idx]].push_back(make_pair(itTransMap->first, weight[idx])); 
    415         if (0 == tmpCounter.count(srcIndex[idx])) 
    416         {           
    417           tmpCounter.insert(srcIndex[idx]); 
    418         
    419           vector<int>& rankSrc = globalIndexOfTransformedElementOnProc_[srcIndex[idx]] ; 
    420           for (int n=0;n<rankSrc.size();++n) 
    421           { 
    422             if (commonProc_.count(rankSrc[n])==1) globalElementIndexOnProc_[elementPositionInGrid][rankSrc[n]].push_back(srcIndex[idx]); 
    423           } 
    424 //          for (int j = 0; j < procContainSrcElementIdx_.size(); ++j) 
    425 //            globalElementIndexOnProc_[elementPositionInGrid][procContainSrcElementIdx_[j]].push_back(srcIndex[idx]); 
    426         } 
    427       } 
    428     } 
    429    
    430     if (!isTransPosEmpty) 
    431     { 
    432       for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx) 
    433       { 
    434         if (idx != elementPositionInGrid) 
    435         { 
    436           std::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it, 
    437                                                                    ite = globalElementIndexOnProc_[idx].end(); 
    438           for (it = itb; it != ite; ++it) it->second[0] = transPos(idxTrans); 
    439         } 
    440       } 
    441     } 
    442  
    443     // Ok, now compute global index of grid source and ones of grid destination 
    444     computeGlobalGridIndexMapping(elementPositionInGrid, 
    445                                   procContainSrcElementIdx_, //srcRank, 
    446                                   src2DstMap, 
    447                                   gridSrc, 
    448                                   gridDst, 
    449                                   globalElementIndexOnProc_, 
    450                                   globaIndexWeightFromSrcToDst); 
    451   }   
    452  } 
    453 CATCH 
    454  
    455 /*! 
    456   Compute mapping of global index of grid source and grid destination 
    457   \param [in] elementPositionInGrid position of element in grid. E.x: grid composed of domain and axis, domain has position 0 and axis 1. 
    458   \param [in] srcRank rank of client from which we demand global index of element source 
    459   \param [in] src2DstMap mapping of global index of element source and global index of element destination 
    460   \param [in] gridSrc Grid source 
    461   \param [in] gridDst Grid destination 
    462   \param [in] globalElementIndexOnProc Global index of element source on different client rank 
    463   \param [out] globaIndexWeightFromSrcToDst Mapping of global index of grid source and grid destination 
    464 */ 
    465 void CGenericAlgorithmTransformation::computeGlobalGridIndexMapping(int elementPositionInGrid, 
    466                                                                    const std::vector<int>& srcRank, 
    467                                                                    std::unordered_map<int, std::vector<std::pair<int,double> > >& src2DstMap, 
    468                                                                    CGrid* gridSrc, 
    469                                                                    CGrid* gridDst, 
    470                                                                    std::vector<std::unordered_map<int,std::vector<size_t> > >& globalElementIndexOnProc, 
    471                                                                    SourceDestinationIndexMap& globaIndexWeightFromSrcToDst) 
    472 TRY 
    473 { 
    474   SourceDestinationIndexMap globaIndexWeightFromSrcToDst_tmp ; 
    475    
    476   CContext* context = CContext::getCurrent(); 
    477   CContextClient* client=context->client; 
    478   int clientRank = client->clientRank; 
    479    
    480   std::vector<CDomain*> domainListSrcP = gridSrc->getDomains(); 
    481   std::vector<CAxis*> axisListSrcP = gridSrc->getAxis(); 
    482   std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars(); 
    483   CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order; 
    484  
    485   size_t nbElement = axisDomainSrcOrder.numElements(); 
    486   std::vector<size_t> nGlobSrc(nbElement); 
    487   size_t globalSrcSize = 1; 
    488   int domainIndex = 0, axisIndex = 0, scalarIndex = 0; 
    489   for (int idx = 0; idx < nbElement; ++idx) 
    490   { 
    491     nGlobSrc[idx] = globalSrcSize; 
    492     int elementDimension = axisDomainSrcOrder(idx); 
    493  
    494     // If this is a domain 
    495     if (2 == elementDimension) 
    496     { 
    497       globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue(); 
    498       ++domainIndex; 
    499     } 
    500     else if (1 == elementDimension) // So it's an axis 
    501     { 
    502       globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue(); 
    503       ++axisIndex; 
    504     } 
    505     else 
    506     { 
    507       globalSrcSize *= 1; 
    508       ++scalarIndex; 
    509     } 
    510   } 
    511  
    512   std::vector<CDomain*> domainListDestP = gridDst->getDomains(); 
    513   std::vector<CAxis*> axisListDestP = gridDst->getAxis(); 
    514   std::vector<CScalar*> scalarListDestP = gridDst->getScalars(); 
    515   CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order; 
    516  
    517   std::vector<size_t> nGlobDst(nbElement); 
    518   size_t globalDstSize = 1; 
    519   domainIndex = axisIndex = scalarIndex = 0; 
    520   set<size_t>  globalSrcStoreIndex ; 
    521    
    522   for (int idx = 0; idx < nbElement; ++idx) 
    523   { 
    524     nGlobDst[idx] = globalDstSize; 
    525     int elementDimension = axisDomainDstOrder(idx); 
    526  
    527     // If this is a domain 
    528     if (2 == elementDimension) 
    529     { 
    530       globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue(); 
    531       ++domainIndex; 
    532     } 
    533     else if (1 == elementDimension) // So it's an axis 
    534     { 
    535       globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue(); 
    536       ++axisIndex; 
    537     } 
    538     else 
    539     { 
    540       globalDstSize *= 1; 
    541       ++scalarIndex; 
    542     } 
    543   } 
    544  
    545   std::map< std::pair<size_t,size_t>, int > rankMap ; 
    546   std::map< std::pair<size_t,size_t>, int >:: iterator rankMapIt ; 
    547    
    548   for (int i = 0; i < srcRank.size(); ++i) 
    549   { 
    550     size_t ssize = 1; 
    551     int rankSrc = srcRank[i]; 
    552     for (int idx = 0; idx < nbElement; ++idx) 
    553     { 
    554       ssize *= (globalElementIndexOnProc[idx][rankSrc]).size(); 
    555     } 
    556  
    557     std::vector<int> idxLoop(nbElement,0); 
    558     std::vector<int> currentIndexSrc(nbElement, 0); 
    559     std::vector<int> currentIndexDst(nbElement, 0); 
    560     int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size(); 
    561     size_t idx = 0; 
    562     while (idx < ssize) 
    563     { 
    564       for (int ind = 0; ind < nbElement; ++ind) 
    565       { 
    566         if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size()) 
    567         { 
    568           idxLoop[ind] = 0; 
    569           ++idxLoop[ind+1]; 
    570         } 
    571  
    572         currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]]; 
    573       } 
    574  
    575       for (int ind = 0; ind < innnerLoopSize; ++ind) 
    576       { 
    577         currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind]; 
    578         int globalElementDstIndexSize = 0; 
    579         if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid])) 
    580         { 
    581           globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size(); 
    582         } 
    583  
    584         std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0); 
    585         size_t globalSrcIndex = 0; 
    586         for (int idxElement = 0; idxElement < nbElement; ++idxElement) 
    587         { 
    588           if (idxElement == elementPositionInGrid) 
    589           { 
    590             for (int k = 0; k < globalElementDstIndexSize; ++k) 
    591             { 
    592               globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement]; 
    593             } 
    594           } 
    595           else 
    596           { 
    597             for (int k = 0; k < globalElementDstIndexSize; ++k) 
    598             { 
    599               globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement]; 
    600             } 
    601           } 
    602           globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement]; 
    603         } 
    604  
    605         for (int k = 0; k < globalElementDstIndexSize; ++k) 
    606         { 
    607            
    608           globaIndexWeightFromSrcToDst_tmp[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second )); 
    609           rankMapIt=rankMap.find(make_pair(globalSrcIndex,globalDstVecIndex[k])) ; 
    610           if (rankMapIt==rankMap.end()) rankMap[make_pair(globalSrcIndex,globalDstVecIndex[k])] = rankSrc ; 
    611           else if (rankSrc==clientRank) rankMapIt->second = rankSrc ; 
    612         } 
    613         ++idxLoop[0]; 
    614       } 
    615       idx += innnerLoopSize; 
    616     } 
    617   } 
    618  
    619   // eliminate redondant global src point owned by differrent processes. 
    620   // Avoid as possible to tranfer data from an other process if the src point is also owned by current process  
    621       int rankSrc ; 
    622       size_t globalSrcIndex ; 
    623       size_t globalDstIndex ; 
    624       double weight ; 
    625   
    626       SourceDestinationIndexMap::iterator rankIt,rankIte ; 
    627       std::unordered_map<size_t, std::vector<std::pair<size_t,double> > >::iterator globalSrcIndexIt, globalSrcIndexIte ; 
    628       std::vector<std::pair<size_t,double> >::iterator vectIt,vectIte ; 
    629     
    630       rankIt=globaIndexWeightFromSrcToDst_tmp.begin() ; rankIte=globaIndexWeightFromSrcToDst_tmp.end() ; 
    631       for(;rankIt!=rankIte;++rankIt) 
    632       { 
    633         rankSrc = rankIt->first ; 
    634         globalSrcIndexIt = rankIt->second.begin() ; globalSrcIndexIte = rankIt->second.end() ; 
    635         for(;globalSrcIndexIt!=globalSrcIndexIte;++globalSrcIndexIt) 
    636         {  
    637           globalSrcIndex = globalSrcIndexIt->first ; 
    638           vectIt = globalSrcIndexIt->second.begin() ; vectIte = globalSrcIndexIt->second.end() ; 
    639           for(vectIt; vectIt!=vectIte; vectIt++) 
    640           { 
    641             globalDstIndex = vectIt->first ; 
    642             weight = vectIt->second ; 
    643             if (eliminateRedondantSrc_) 
    644             { 
    645               if (rankMap[make_pair(globalSrcIndex,globalDstIndex)] == rankSrc)   
    646                 globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstIndex,weight)) ; 
    647             } 
    648             else globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstIndex,weight)) ; 
    649          } 
    650        } 
    651      } 
    652 } 
    653 CATCH 
    654  
    655 /*! 
    656   Find out proc and global index of axis source which axis destination is on demande 
    657   \param[in] scalar Scalar destination 
    658   \param[in] scalar Scalar source 
    659   \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid. 
    660   \param[out] globalScalarIndexOnProc Global index of axis source on different procs 
    661 */ 
    662 void CGenericAlgorithmTransformation::computeExchangeScalarIndex(CScalar* scalarDst, 
    663                                                                  CScalar* scalarSrc, 
    664                                                                  CArray<size_t,1>& destGlobalIndexPositionInGrid, 
    665                                                                  std::unordered_map<int,std::vector<size_t> >& globalScalarIndexOnProc) 
    666 TRY 
    667 { 
    668   CContext* context = CContext::getCurrent(); 
    669   CContextClient* client=context->client; 
    670   int clientRank = client->clientRank; 
    671   int clientSize = client->clientSize; 
    672  
    673   globalScalarIndexOnProc.rehash(std::ceil(clientSize/globalScalarIndexOnProc.max_load_factor())); 
    674   for (int idx = 0; idx < clientSize; ++idx) 
    675   { 
    676     globalScalarIndexOnProc[idx].push_back(0); 
    677   } 
    678 } 
    679 CATCH 
    680  
    681 /*! 
    682   Find out proc and global index of axis source which axis destination is on demande 
    683   \param[in] axisDst Axis destination 
    684   \param[in] axisSrc Axis source 
    685   \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid. 
    686   \param[out] globalAxisIndexOnProc Global index of axis source on different procs 
    687 */ 
    688 void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst, 
    689                                                                CAxis* axisSrc, 
    690                                                                CArray<size_t,1>& destGlobalIndexPositionInGrid, 
    691                                                                std::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc) 
    692 TRY 
    693 { 
    694   CContext* context = CContext::getCurrent(); 
    695   CContextClient* client=context->client; 
    696   int clientRank = client->clientRank; 
    697   int clientSize = client->clientSize; 
    698  
    699   size_t globalIndex; 
    700   int nIndexSize = axisSrc->index.numElements(); 
    701   CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank; 
    702   globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor())); 
    703   for (int idx = 0; idx < nIndexSize; ++idx) 
    704   { 
    705     if (axisSrc->mask(idx)) 
    706     { 
    707       globalIndex = axisSrc->index(idx); 
    708       globalIndex2ProcRank[globalIndex].push_back(clientRank); 
    709     } 
    710   } 
    711  
    712   CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm); 
    713   CArray<size_t,1> globalAxisIndex(axisDst->index.numElements()); 
    714   for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx) 
    715   { 
    716     globalAxisIndex(idx) = axisDst->index(idx); 
    717   } 
    718   dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex); 
    719  
    720   std::vector<int> countIndex(clientSize,0); 
    721   const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap(); 
    722   CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it, 
    723                                                                ite = computedGlobalIndexOnProc.end(); 
    724   for (it = itb; it != ite; ++it) 
    725   { 
    726     const std::vector<int>& procList = it->second; 
    727     for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]]; 
    728   } 
    729  
    730   globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor())); 
    731   for (int idx = 0; idx < clientSize; ++idx) 
    732   { 
    733     if (0 != countIndex[idx]) 
    734     { 
    735       globalAxisIndexOnProc[idx].resize(countIndex[idx]); 
    736       countIndex[idx] = 0; 
    737     } 
    738   } 
    739  
    740   for (it = itb; it != ite; ++it) 
    741   { 
    742     const std::vector<int>& procList = it->second; 
    743     for (int idx = 0; idx < procList.size(); ++idx) 
    744     { 
    745       globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first; 
    746       ++countIndex[procList[idx]]; 
    747     } 
    748   } 
    749 } 
    750 CATCH 
    751  
    752 /*! 
    753   Find out proc and global index of domain source which domain destination is on demande 
    754   \param[in] domainDst Domain destination 
    755   \param[in] domainSrc Domain source 
    756   \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid. 
    757   \param[out] globalDomainIndexOnProc Global index of domain source on different procs 
    758 */ 
    759 void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst, 
    760                                                                  CDomain* domainSrc, 
    761                                                                  CArray<size_t,1>& destGlobalIndexPositionInGrid, 
    762                                                                  std::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc) 
    763 TRY 
    764 { 
    765   CContext* context = CContext::getCurrent(); 
    766   CContextClient* client=context->client; 
    767   int clientRank = client->clientRank; 
    768   int clientSize = client->clientSize; 
    769  
    770   int niGlobSrc = domainSrc->ni_glo.getValue(); 
    771   size_t globalIndex; 
    772   int i_ind, j_ind; 
    773   int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements() 
    774                                                              : destGlobalIndexPositionInGrid.numElements(); 
    775   CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank; 
    776   globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor())); 
    777    
    778   if (destGlobalIndexPositionInGrid.isEmpty()) 
    779   { 
    780     for (int idx = 0; idx < nIndexSize; ++idx) 
    781     { 
    782       i_ind=domainSrc->i_index(idx) ; 
    783       j_ind=domainSrc->j_index(idx) ; 
    784  
    785       if (domainSrc->localMask(idx)) 
    786       { 
    787         globalIndex = i_ind + j_ind * niGlobSrc; 
    788         globalIndex2ProcRank[globalIndex].resize(1); 
    789         globalIndex2ProcRank[globalIndex][0] = clientRank; 
    790       } 
    791     } 
    792   } 
    793   else 
    794   { 
    795     for (int idx = 0; idx < nIndexSize; ++idx) 
    796     { 
    797 //      if (domainSrc->localMask(idx)) -> not necessairy, mask seem to be included in  destGlobalIndexPositionInGrid(idx)    (ym) 
    798         globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank); 
    799     } 
    800   } 
    801  
    802   CArray<size_t,1> globalDomainIndex; 
    803   if (destGlobalIndexPositionInGrid.isEmpty()) 
    804   { 
    805     int niGlobDst = domainDst->ni_glo.getValue(); 
    806     globalDomainIndex.resize(domainDst->i_index.numElements()); 
    807     nIndexSize = domainDst->i_index.numElements(); 
    808  
    809     for (int idx = 0; idx < nIndexSize; ++idx) 
    810     { 
    811       i_ind=domainDst->i_index(idx) ; 
    812       j_ind=domainDst->j_index(idx) ; 
    813  
    814       globalDomainIndex(idx) = i_ind + j_ind * niGlobDst; 
    815     } 
    816   } 
    817   else 
    818   { 
    819     globalDomainIndex.reference(destGlobalIndexPositionInGrid); 
    820   } 
    821  
    822   CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm); 
    823   dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex); 
    824  
    825   std::vector<int> countIndex(clientSize,0); 
    826   const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap(); 
    827   CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it, 
    828                                                                ite = computedGlobalIndexOnProc.end(); 
    829   for (it = itb; it != ite; ++it) 
    830   { 
    831     const std::vector<int>& procList = it->second; 
    832     for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]]; 
    833   } 
    834  
    835   globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor())); 
    836   for (int idx = 0; idx < clientSize; ++idx) 
    837   { 
    838     if (0 != countIndex[idx]) 
    839     { 
    840       globalDomainIndexOnProc[idx].resize(countIndex[idx]); 
    841       countIndex[idx] = 0; 
    842     } 
    843   } 
    844  
    845   for (it = itb; it != ite; ++it) 
    846   { 
    847     const std::vector<int>& procList = it->second; 
    848     for (int idx = 0; idx < procList.size(); ++idx) 
    849     { 
    850       globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first; 
    851       ++countIndex[procList[idx]]; 
    852     } 
    853   } 
    854 } 
    855 CATCH 
    856  
    857 void CGenericAlgorithmTransformation::computeTransformationMappingNonDistributed(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst, 
    858                                                                                  vector<int>& localSrc, vector<int>& localDst, vector<double>& weight, 
    859                                                                                  int& nlocalIndexDest) 
    860 TRY 
    861 { 
    862  
    863   CContext* context = CContext::getCurrent(); 
    864   CContextClient* client = context->client; 
    865   int nbClient = client->clientSize; 
    866  
    867   computePositionElements(gridDst, gridSrc); 
    868   std::vector<CScalar*> scalarListDstP = gridDst->getScalars(); 
    869   std::vector<CAxis*> axisListDstP = gridDst->getAxis(); 
    870   std::vector<CDomain*> domainListDstP = gridDst->getDomains(); 
    871   CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order; 
    872   std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars(); 
    873   std::vector<CAxis*> axisListSrcP = gridSrc->getAxis(); 
    874   std::vector<CDomain*> domainListSrcP = gridSrc->getDomains(); 
    875   CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;   
    876  
    877   int nElement=axisDomainSrcOrder.numElements() ; 
    878   int indSrc=1 ; 
    879   int indDst=1 ; 
    880   vector<int> nIndexSrc(nElement) ; 
    881   vector<int> nIndexDst(nElement) ; 
    882   vector< CArray<bool,1>* > maskSrc(nElement) ; 
    883   vector< CArray<bool,1>* > maskDst(nElement) ; 
    884      
    885   int nlocalIndexSrc=1 ; 
    886 //  int nlocalIndexDest=1 ; 
    887   nlocalIndexDest=1 ; 
    888   CArray<bool,1> maskScalar(1) ; 
    889   maskScalar  = true ; 
    890    
    891    
    892   for(int i=0 ; i<nElement; i++) 
    893   { 
    894     int dimElement = axisDomainSrcOrder(i); 
    895     if (2 == dimElement) //domain 
    896     { 
    897       CDomain* domain=domainListSrcP[elementPositionInGridSrc2DomainPosition_[i]] ; 
    898       nIndexSrc[i] = domain->i_index.numElements() ; 
    899       maskSrc[i]=&domain->localMask ; 
    900     } 
    901     else if (1 == dimElement) //axis 
    902     { 
    903       CAxis* axis=axisListSrcP[elementPositionInGridSrc2AxisPosition_[i]] ; 
    904       nIndexSrc[i] = axis->index.numElements() ; 
    905       maskSrc[i]=&axis->mask ; 
    906     } 
    907     else  //scalar 
    908     { 
    909       nIndexSrc[i]=1 ; 
    910       maskSrc[i]=&maskScalar ; 
    911     } 
    912     nlocalIndexSrc=nlocalIndexSrc*nIndexSrc[i] ; 
    913   } 
    91439 
    91540 
    91641 
    917   int offset=1 ; 
    918   for(int i=0 ; i<nElement; i++) 
    919   { 
    920     int dimElement = axisDomainDstOrder(i); 
    921     if (2 == dimElement) //domain 
    922     { 
    923       CDomain* domain=domainListDstP[elementPositionInGridDst2DomainPosition_[i]] ; 
    924       int nIndex=domain->i_index.numElements() ; 
    925       CArray<bool,1>& localMask=domain->localMask ; 
    926       int nbInd=0 ; 
    927       for(int j=0;j<nIndex;j++) if (localMask(j)) nbInd++ ; 
    928       nIndexDst[i] = nbInd ; 
    929       maskDst[i]=&domain->localMask ; 
    930     } 
    931     else if (1 == dimElement) //axis 
    932     { 
    933       CAxis* axis = axisListDstP[elementPositionInGridDst2AxisPosition_[i]] ; 
    934       int nIndex=axis->index.numElements() ; 
    935       CArray<bool,1>& localMask=axis->mask ; 
    936       int nbInd=0 ; 
    937       for(int j=0;j<nIndex;j++) if (localMask(j)) nbInd++ ; 
    938       nIndexDst[i] = nbInd ; 
    939       maskDst[i]=&axis->mask ; 
    940     } 
    941     else  //scalar 
    942     { 
    943       nIndexDst[i]=1 ; 
    944       maskDst[i]=&maskScalar ; 
    945     } 
    946     if (i<elementPositionInGrid) offset=offset*nIndexDst[i] ; 
    947     nlocalIndexDest=nlocalIndexDest*nIndexDst[i] ; 
    948   } 
     42CTransformFilter* CGenericAlgorithmTransformation::createTransformFilter(CGarbageCollector& gc, CGridAlgorithm* algo, bool detectMissingValues, double defaultValue) 
     43{ 
     44  return new CTransformFilter(gc, 1, algo, detectMissingValues, defaultValue) ; 
     45} 
    94946 
    950   vector<int> dstLocalInd ; 
    951   int dimElement = axisDomainDstOrder(elementPositionInGrid); 
    952   if (2 == dimElement) //domain 
    953   { 
    954     CDomain* domain = domainListDstP[elementPositionInGridDst2DomainPosition_[elementPositionInGrid]] ; 
    955     int ni_glo=domain->ni_glo ; 
    956     int nj_glo=domain->nj_glo ; 
    957     int nindex_glo=ni_glo*nj_glo ; 
    958     dstLocalInd.resize(nindex_glo,-1) ; 
    959     int nIndex=domain->i_index.numElements() ; 
    960     CArray<bool,1>& localMask=domain->localMask ; 
    961     int unmaskedInd=0 ; 
    962     int globIndex ; 
    963     for(int i=0;i<nIndex;i++) 
    964     { 
    965       if (localMask(i)) 
    966       { 
    967         globIndex=domain->j_index(i)*ni_glo+domain->i_index(i) ; 
    968         dstLocalInd[globIndex]=unmaskedInd ; 
    969         unmaskedInd++ ; 
    970       } 
    971     } 
    972   } 
    973   else if (1 == dimElement) //axis 
    974   { 
    975     CAxis* axis = axisListDstP[elementPositionInGridDst2AxisPosition_[elementPositionInGrid]] ; 
    976     int nindex_glo=axis->n_glo ; 
    977     dstLocalInd.resize(nindex_glo,-1) ; 
    978     int nIndex=axis->index.numElements() ; 
    979     CArray<bool,1>& localMask=axis->mask ; // axis mask must include later data_index 
    980     int unmaskedInd=0 ; 
    981     for(int i=0;i<nIndex;i++) 
    982     { 
    983       if (localMask(i)) 
    984       { 
    985         dstLocalInd[axis->index(i)]=unmaskedInd ; 
    986         unmaskedInd++ ; 
    987       } 
    988     } 
    989   } 
    990   else  //scalar 
    991   { 
    992     dstLocalInd.resize(1) ;  
    993     dstLocalInd[0]=0 ;  
    994   } 
    995  
    996   vector<vector<vector<pair<int,double> > > > dstIndWeight(transformationMapping_.size()) ; 
    997     
    998   for(int t=0;t<transformationMapping_.size();++t) 
    999   { 
    1000     TransformationIndexMap::const_iterator   itTransMap = transformationMapping_[t].begin(), 
    1001                                              iteTransMap = transformationMapping_[t].end(); 
    1002     TransformationWeightMap::const_iterator itTransWeight = transformationWeight_[t].begin(); 
    1003     dstIndWeight[t].resize(nIndexSrc[elementPositionInGrid]) ; 
    1004      
    1005     for(;itTransMap!=iteTransMap;++itTransMap,++itTransWeight) 
    1006     { 
    1007       int dst=dstLocalInd[itTransMap->first] ; 
    1008       if (dst!=-1) 
    1009       { 
    1010         const vector<int>& srcs=itTransMap->second; 
    1011         const vector<double>& weights=itTransWeight->second; 
    1012         for(int i=0;i<srcs.size() ;i++) dstIndWeight[t][srcs[i]].push_back(pair<int,double>(dst*offset+t,weights[i])) ; 
    1013       } 
    1014     } 
    1015   } 
    1016   int srcInd=0 ;   
    1017   int currentInd ; 
    1018   int t=0 ;   
    1019   int srcIndCompressed=0 ; 
    1020    
    1021   nonDistributedrecursiveFunct(nElement-1,true,elementPositionInGrid,maskSrc,maskDst, srcInd, srcIndCompressed, nIndexSrc, t, dstIndWeight,   
    1022                                currentInd,localSrc,localDst,weight); 
    1023                 
     47vector<string> CGenericAlgorithmTransformation::getAuxFieldId(void)  
     48{ 
     49  return vector<string>() ; 
    102450} 
    1025 CATCH 
    1026  
    1027  
    1028 void CGenericAlgorithmTransformation::nonDistributedrecursiveFunct(int currentPos, bool masked, int elementPositionInGrid, 
    1029                                                                    vector< CArray<bool,1>* >& maskSrc, vector< CArray<bool,1>* >& maskDst, 
    1030                                                                    int& srcInd, int& srcIndCompressed, vector<int>& nIndexSrc, 
    1031                                                                    int& t, vector<vector<vector<pair<int,double> > > >& dstIndWeight, int currentInd, 
    1032                                                                    vector<int>& localSrc, vector<int>& localDst, vector<double>& weight) 
    1033 TRY 
    1034 { 
    1035   int masked_ ; 
    1036   if (currentPos!=elementPositionInGrid) 
    1037   { 
    1038     if (currentPos!=0) 
    1039     { 
    1040       CArray<bool,1>& mask = *maskSrc[currentPos] ; 
    1041       
    1042       for(int i=0;i<nIndexSrc[currentPos];i++) 
    1043       { 
    1044         masked_=masked ; 
    1045         if (!mask(i)) masked_=false ; 
    1046         nonDistributedrecursiveFunct(currentPos-1, masked_, elementPositionInGrid, maskSrc, maskDst, srcInd, srcIndCompressed, nIndexSrc, t, 
    1047                                      dstIndWeight, currentInd, localSrc, localDst, weight); 
    1048       } 
    1049     } 
    1050     else 
    1051     { 
    1052       CArray<bool,1>& mask = *maskSrc[currentPos] ; 
    1053       for(int i=0;i<nIndexSrc[currentPos];i++) 
    1054       { 
    1055         if (masked && mask(i)) 
    1056         { 
    1057           if (dstIndWeight[t][currentInd].size()>0) 
    1058           { 
    1059             for(vector<pair<int,double> >::iterator it = dstIndWeight[t][currentInd].begin(); it!=dstIndWeight[t][currentInd].end(); ++it) 
    1060             { 
    1061               localSrc.push_back(srcIndCompressed) ; 
    1062               localDst.push_back(it->first) ; 
    1063               weight.push_back(it->second) ; 
    1064               (it->first)++ ; 
    1065             } 
    1066           } 
    1067           if (t < dstIndWeight.size()-1) t++ ; 
    1068             srcIndCompressed ++ ; 
    1069         } 
    1070         srcInd++ ; 
    1071       } 
    1072     } 
    1073   } 
    1074   else 
    1075   { 
    1076   
    1077     if (currentPos!=0) 
    1078     { 
    1079  
    1080       CArray<bool,1>& mask = *maskSrc[currentPos] ; 
    1081       for(int i=0;i<nIndexSrc[currentPos];i++) 
    1082       { 
    1083         t=0 ; 
    1084         masked_=masked ; 
    1085         if (!mask(i)) masked_=false ;  
    1086         nonDistributedrecursiveFunct(currentPos-1, masked_, elementPositionInGrid, maskSrc, maskDst, srcInd, 
    1087                                      srcIndCompressed, nIndexSrc, t, dstIndWeight , i,  localSrc, localDst, weight); 
    1088       } 
    1089     } 
    1090     else 
    1091     { 
    1092       for(int i=0;i<nIndexSrc[currentPos];i++) 
    1093       { 
    1094         if (masked) 
    1095         { 
    1096           t=0 ;         
    1097           if (dstIndWeight[t][i].size()>0) 
    1098           { 
    1099             for(vector<pair<int,double> >::iterator it = dstIndWeight[t][i].begin(); it!=dstIndWeight[t][i].end(); ++it) 
    1100             { 
    1101               localSrc.push_back(srcIndCompressed) ; 
    1102               localDst.push_back(it->first) ; 
    1103               weight.push_back(it->second) ; 
    1104               (it->first)++ ; 
    1105             } 
    1106            } 
    1107           if (t < dstIndWeight.size()-1) t++ ; 
    1108           srcIndCompressed ++ ; 
    1109         } 
    1110         srcInd++ ; 
    1111       } 
    1112     } 
    1113   } 
    111451 
    111552} 
    1116 CATCH 
    1117  
    1118 /*! 
    1119   Compute index mapping between element source and element destination with an auxiliary inputs which determine 
    1120 position of each mapped index in global index of grid destination. 
    1121   \param [in] dataAuxInputs auxiliary inputs 
    1122 */ 
    1123 void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    1124 TRY 
    1125 { 
    1126   computeIndexSourceMapping_(dataAuxInputs); 
    1127 } 
    1128 CATCH 
    1129  
    1130 std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs() 
    1131 TRY 
    1132 { 
    1133   return idAuxInputs_; 
    1134 } 
    1135 CATCH 
    1136  
    1137 CGenericAlgorithmTransformation::AlgoTransType CGenericAlgorithmTransformation::type() 
    1138 TRY 
    1139 { 
    1140   return type_; 
    1141 } 
    1142 CATCH 
    1143  
    1144 } 
Note: See TracChangeset for help on using the changeset viewer.