Ignore:
Timestamp:
01/12/21 23:05:02 (3 years ago)
Author:
ymipsl
Message:
  • bug fix when createing mask on server side when overlapping grid
  • implement axis interpolation on pressure coordinate
  • big cleaning in transformation

YM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/generic_algorithm_transformation.cpp

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