source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/generic_algorithm_transformation.cpp @ 1438

Last change on this file since 1438 was 1438, checked in by ymipsl, 6 years ago

Non Distributed transformation : fix small typo that may lead to big problems...

YM

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