source: XIOS/trunk/src/transformation/generic_algorithm_transformation.cpp @ 1018

Last change on this file since 1018 was 1018, checked in by mhnguyen, 4 years ago

Improving missing-value processing
If detect_missing_value is activated, then all missing value will be converted to
NaN (Not-a-number) in input of data flow then they will be reconverted to missing value on output

+) Update SourceFilter?, TemporalFilter? and SpatialTransformFilter? with new processing
+) Update all transformations with new processing

Test
+) On Curie
+) Work

File size: 25.7 KB
Line 
1/*!
2   \file generic_algorithm_transformation.hpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 21 Mars 2016
6
7   \brief Interface for all transformation algorithms.
8 */
9#include "generic_algorithm_transformation.hpp"
10#include "context.hpp"
11#include "context_client.hpp"
12#include "client_client_dht_template.hpp"
13#include "utils.hpp"
14
15namespace xios {
16
17CGenericAlgorithmTransformation::CGenericAlgorithmTransformation()
18 : transformationMapping_(), transformationWeight_(), transformationPosition_(),
19   idAuxInputs_(), type_(ELEMENT_NO_MODIFICATION_WITH_DATA)
20{
21}
22
23void CGenericAlgorithmTransformation::updateData(CArray<double,1>& dataOut)
24{
25
26}
27
28void CGenericAlgorithmTransformation::apply(const std::vector<std::pair<int,double> >& localIndex,
29                                            const double* dataInput,
30                                            CArray<double,1>& dataOut,
31                                            std::vector<bool>& flagInitial,
32                                            const double& defaultValue)
33{
34  int nbLocalIndex = localIndex.size(); 
35  bool hasMissingValue = NumTraits<double>::isnan(defaultValue); 
36  if (hasMissingValue)
37  {
38    for (int idx = 0; idx < nbLocalIndex; ++idx)
39    {
40      if (NumTraits<double>::isnan(*(dataInput + idx)))
41      {
42        flagInitial[localIndex[idx].first] = false;
43      }
44      else
45      {
46        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
47      }
48    }
49
50    for (int idx = 0; idx < nbLocalIndex; ++idx)
51    {
52      if (!flagInitial[localIndex[idx].first])
53        dataOut(localIndex[idx].first) = defaultValue;
54    }
55  }
56  else
57  {
58    for (int idx = 0; idx < nbLocalIndex; ++idx)
59    {
60      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
61    }
62  }
63}
64
65void CGenericAlgorithmTransformation::computePositionElements(CGrid* dst, CGrid* src)
66{
67  int idxScalar = 0, idxAxis = 0, idxDomain = 0;
68  CArray<int,1> axisDomainOrderDst = dst->axis_domain_order;
69  for (int i = 0; i < axisDomainOrderDst.numElements(); ++i)
70  {
71    int dimElement = axisDomainOrderDst(i);
72    if (2 == dimElement)
73    {
74      elementPositionInGridDst2DomainPosition_[i] = idxDomain;
75      ++idxDomain;
76    }
77    else if (1 == dimElement)
78    {
79      elementPositionInGridDst2AxisPosition_[i] = idxAxis;
80      ++idxAxis;
81    }
82    else
83    {
84      elementPositionInGridDst2ScalarPosition_[i] = idxScalar;
85      ++idxScalar;
86    }
87  }
88
89  idxScalar = idxAxis = idxDomain = 0;
90  CArray<int,1> axisDomainOrderSrc = src->axis_domain_order;
91  for (int i = 0; i < axisDomainOrderSrc.numElements(); ++i)
92  {
93    int dimElement = axisDomainOrderSrc(i);
94    if (2 == dimElement)
95    {
96      elementPositionInGridSrc2DomainPosition_[i] = idxDomain;
97      ++idxDomain;
98    }
99    else if (1 == dimElement)
100    {
101      elementPositionInGridSrc2AxisPosition_[i] = idxAxis;
102      ++idxAxis;
103    }
104    else
105    {
106      elementPositionInGridSrc2ScalarPosition_[i] = idxScalar;
107      ++idxScalar;
108    }
109  }
110}
111
112/*!
113  This function computes the global indexes of grid source, which the grid destination is in demand.
114  \param[in] elementPositionInGrid position of an element in a grid .E.g: if grid is composed of domain and axis (in order),
115                then position of axis in grid is 1 and domain is positioned at 0.
116  \param[in] gridSrc Grid source
117  \param[in] gridDst Grid destination
118  \param[in\out] globaIndexWeightFromSrcToDst mapping of each global index source and weight to index destination
119*/
120void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid,
121                                                               CGrid* gridSrc,
122                                                               CGrid* gridDst,
123                                                               SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
124 {
125  CContext* context = CContext::getCurrent();
126  CContextClient* client = context->client;
127  int nbClient = client->clientSize;
128
129  typedef boost::unordered_map<int, std::vector<std::pair<int,double> > > SrcToDstMap;
130
131  size_t indexSrcSize = 0;
132  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
133  {
134    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
135                                           iteTransMap = transformationMapping_[idxTrans].end();
136    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
137
138    itTransWeight = itbTransWeight;
139    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
140    {
141       indexSrcSize += (itTransMap->second).size();
142    }
143  }
144
145  bool isTransPosEmpty = transformationPosition_.empty();
146  CArray<size_t,1> transPos;
147  if (!isTransPosEmpty) transPos.resize(transformationMapping_.size());
148  CArray<size_t,1> indexSrc(indexSrcSize);
149  indexSrcSize = 0;
150  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
151  {
152    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
153                                           iteTransMap = transformationMapping_[idxTrans].end();
154    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
155
156    // Build mapping between global source element index and global destination element index.
157    itTransWeight = itbTransWeight;
158    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
159    {
160      const std::vector<int>& srcIndex = itTransMap->second;
161      for (int idx = 0; idx < srcIndex.size(); ++idx)
162      {
163        indexSrc(indexSrcSize) = srcIndex[idx];
164        ++indexSrcSize;
165      }
166    }
167
168    if (!isTransPosEmpty)
169    {
170      TransformationPositionMap::const_iterator itPosMap = transformationPosition_[idxTrans].begin();
171      transPos(idxTrans) = itPosMap->second[0];
172    }
173  }
174
175  // compute position of elements on grids
176  computePositionElements(gridDst, gridSrc);
177  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
178  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
179  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
180  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
181  std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
182  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
183  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
184  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
185
186  // Find out global index source of transformed element on corresponding process.
187  std::vector<boost::unordered_map<int,std::vector<size_t> > > globalElementIndexOnProc(axisDomainDstOrder.numElements());
188  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc; 
189  for (int idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
190  {
191    if (idx == elementPositionInGrid)
192      computeExchangeGlobalIndex(indexSrc, axisDomainSrcOrder(idx), globalIndexOfTransformedElementOnProc); //globalElementIndexOnProc[idx]);
193    if (2 == axisDomainDstOrder(idx)) // It's domain
194    {
195      if (idx != elementPositionInGrid)
196        computeExchangeDomainIndex(domainListDestP[elementPositionInGridDst2DomainPosition_[idx]],
197                                   domainListSrcP[elementPositionInGridSrc2DomainPosition_[idx]],
198                                   transPos,
199                                   globalElementIndexOnProc[idx]);     
200
201    }
202    else if (1 == axisDomainDstOrder(idx))//it's an axis
203    {
204      if (idx != elementPositionInGrid)
205        computeExchangeAxisIndex(axisListDestP[elementPositionInGridDst2AxisPosition_[idx]],
206                                 axisListSrcP[elementPositionInGridSrc2AxisPosition_[idx]],
207                                 transPos,
208                                 globalElementIndexOnProc[idx]);
209    }
210    else //it's a scalar
211    {
212      if (idx != elementPositionInGrid)
213        computeExchangeScalarIndex(scalarListDestP[elementPositionInGridDst2ScalarPosition_[idx]],
214                                   scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[idx]],
215                                   transPos,
216                                   globalElementIndexOnProc[idx]);
217
218    }
219  }
220
221  if (!isTransPosEmpty)
222  {
223    for (int idx = 0; idx < globalElementIndexOnProc.size(); ++idx)
224    {
225      if (idx != elementPositionInGrid)
226      {
227        boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc[idx].begin(), it,
228                                                                 ite = globalElementIndexOnProc[idx].end();
229        for (it = itb; it != ite; ++it) it->second.resize(1);
230      }
231    }
232  }
233
234  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
235  {
236    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
237                                           iteTransMap = transformationMapping_[idxTrans].end();
238    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
239    SrcToDstMap src2DstMap;
240    src2DstMap.rehash(std::ceil(transformationMapping_[idxTrans].size()/src2DstMap.max_load_factor()));
241
242    // Build mapping between global source element index and global destination element index.
243    boost::unordered_map<int,std::vector<size_t> >().swap(globalElementIndexOnProc[elementPositionInGrid]);
244    boost::unordered_map<int,int> tmpCounter;
245    itTransWeight = itbTransWeight;
246    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
247    {
248      const std::vector<int>& srcIndex = itTransMap->second;
249      const std::vector<double>& weight = itTransWeight->second;
250      for (int idx = 0; idx < srcIndex.size(); ++idx)
251      {
252        src2DstMap[srcIndex[idx]].push_back(make_pair(itTransMap->first, weight[idx]));
253        if (1 == globalIndexOfTransformedElementOnProc.count(srcIndex[idx]) && (0 == tmpCounter.count(srcIndex[idx])))
254        {
255          tmpCounter[srcIndex[idx]] = 1;
256          std::vector<int>& srcProc = globalIndexOfTransformedElementOnProc[srcIndex[idx]];
257          for (int j = 0; j < srcProc.size(); ++j)
258            globalElementIndexOnProc[elementPositionInGrid][srcProc[j]].push_back(srcIndex[idx]);
259        }
260      }
261    }
262
263    if (!isTransPosEmpty)
264    {
265      for (int idx = 0; idx < globalElementIndexOnProc.size(); ++idx)
266      {
267        if (idx != elementPositionInGrid)
268        {
269          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc[idx].begin(), it,
270                                                                   ite = globalElementIndexOnProc[idx].end();
271          for (it = itb; it != ite; ++it) it->second[0] = transPos(idxTrans);
272        }
273      }
274    }
275
276    std::vector<std::vector<bool> > elementOnProc(axisDomainDstOrder.numElements(), std::vector<bool>(nbClient, false));
277    boost::unordered_map<int,std::vector<size_t> >::const_iterator it, itb, ite;
278    for (int idx = 0; idx < globalElementIndexOnProc.size(); ++idx)
279    {
280      itb = globalElementIndexOnProc[idx].begin();
281      ite = globalElementIndexOnProc[idx].end();
282      for (it = itb; it != ite; ++it) elementOnProc[idx][it->first] = true;
283    }
284
285    // Determine procs which contain global source index
286    std::vector<bool> intersectedProc(nbClient, true);
287    for (int idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
288    {
289      std::transform(elementOnProc[idx].begin(), elementOnProc[idx].end(),
290                     intersectedProc.begin(), intersectedProc.begin(),
291                     std::logical_and<bool>());
292    }
293
294    std::vector<int> srcRank;
295    for (int idx = 0; idx < nbClient; ++idx)
296    {
297      if (intersectedProc[idx]) srcRank.push_back(idx);
298    }
299
300    // Ok, now compute global index of grid source and ones of grid destination
301    computeGlobalGridIndexMapping(elementPositionInGrid,
302                                  srcRank,
303                                  src2DstMap,
304                                  gridSrc,
305                                  gridDst,
306                                  globalElementIndexOnProc,
307                                  globaIndexWeightFromSrcToDst);
308  }
309 }
310
311/*!
312  Compute mapping of global index of grid source and grid destination
313  \param [in] elementPositionInGrid position of element in grid. E.x: grid composed of domain and axis, domain has position 0 and axis 1.
314  \param [in] srcRank rank of client from which we demand global index of element source
315  \param [in] src2DstMap mapping of global index of element source and global index of element destination
316  \param[in] gridSrc Grid source
317  \param[in] gridDst Grid destination
318  \param[in] globalElementIndexOnProc Global index of element source on different client rank
319  \param[out] globaIndexWeightFromSrcToDst Mapping of global index of grid source and grid destination
320*/
321void CGenericAlgorithmTransformation::computeGlobalGridIndexMapping(int elementPositionInGrid,
322                                                                   const std::vector<int>& srcRank,
323                                                                   boost::unordered_map<int, std::vector<std::pair<int,double> > >& src2DstMap,
324                                                                   CGrid* gridSrc,
325                                                                   CGrid* gridDst,
326                                                                   std::vector<boost::unordered_map<int,std::vector<size_t> > >& globalElementIndexOnProc,
327                                                                   SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
328{
329  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
330  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
331  std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();
332  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
333
334  size_t nbElement = axisDomainSrcOrder.numElements();
335  std::vector<size_t> nGlobSrc(nbElement);
336  size_t globalSrcSize = 1;
337  int domainIndex = 0, axisIndex = 0, scalarIndex = 0;
338  for (int idx = 0; idx < nbElement; ++idx)
339  {
340    nGlobSrc[idx] = globalSrcSize;
341    int elementDimension = axisDomainSrcOrder(idx);
342
343    // If this is a domain
344    if (2 == elementDimension)
345    {
346      globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue();
347      ++domainIndex;
348    }
349    else if (1 == elementDimension) // So it's an axis
350    {
351      globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue();
352      ++axisIndex;
353    }
354    else
355    {
356      globalSrcSize *= 1;
357      ++scalarIndex;
358    }
359  }
360
361  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
362  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
363  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
364  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
365
366  std::vector<size_t> nGlobDst(nbElement);
367  size_t globalDstSize = 1;
368  domainIndex = axisIndex = scalarIndex = 0;
369  for (int idx = 0; idx < nbElement; ++idx)
370  {
371    nGlobDst[idx] = globalDstSize;
372    int elementDimension = axisDomainDstOrder(idx);
373
374    // If this is a domain
375    if (2 == elementDimension)
376    {
377      globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue();
378      ++domainIndex;
379    }
380    else if (1 == elementDimension) // So it's an axis
381    {
382      globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue();
383      ++axisIndex;
384    }
385    else
386    {
387      globalDstSize *= 1;
388      ++scalarIndex;
389    }
390  }
391
392  for (int i = 0; i < srcRank.size(); ++i)
393  {
394    size_t ssize = 1;
395    int rankSrc = srcRank[i];
396    for (int idx = 0; idx < nbElement; ++idx)
397    {
398      ssize *= (globalElementIndexOnProc[idx][rankSrc]).size();
399    }
400
401    std::vector<int> idxLoop(nbElement,0);
402    std::vector<int> currentIndexSrc(nbElement, 0);
403    std::vector<int> currentIndexDst(nbElement, 0);
404    int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size();
405    size_t idx = 0;
406    while (idx < ssize)
407    {
408      for (int ind = 0; ind < nbElement; ++ind)
409      {
410        if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size())
411        {
412          idxLoop[ind] = 0;
413          ++idxLoop[ind+1];
414        }
415
416        currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]];
417      }
418
419      for (int ind = 0; ind < innnerLoopSize; ++ind)
420      {
421        currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind];
422        int globalElementDstIndexSize = 0;
423        if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid]))
424        {
425          globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size();
426        }
427
428        std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0);
429        size_t globalSrcIndex = 0;
430        for (int idxElement = 0; idxElement < nbElement; ++idxElement)
431        {
432          if (idxElement == elementPositionInGrid)
433          {
434            for (int k = 0; k < globalElementDstIndexSize; ++k)
435            {
436              globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement];
437            }
438          }
439          else
440          {
441            for (int k = 0; k < globalElementDstIndexSize; ++k)
442            {
443              globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement];
444            }
445          }
446          globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement];
447        }
448
449        for (int k = 0; k < globalElementDstIndexSize; ++k)
450        {
451          globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second ));
452        }
453        ++idxLoop[0];
454      }
455      idx += innnerLoopSize;
456    }
457  }
458}
459
460/*!
461  Find out proc and global index of axis source which axis destination is on demande
462  \param[in] scalar Scalar destination
463  \param[in] scalar Scalar source
464  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
465  \param[out] globalScalarIndexOnProc Global index of axis source on different procs
466*/
467void CGenericAlgorithmTransformation::computeExchangeScalarIndex(CScalar* scalarDst,
468                                                                 CScalar* scalarSrc,
469                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
470                                                                 boost::unordered_map<int,std::vector<size_t> >& globalScalarIndexOnProc)
471{
472  CContext* context = CContext::getCurrent();
473  CContextClient* client=context->client;
474  int clientRank = client->clientRank;
475  int clientSize = client->clientSize;
476
477  globalScalarIndexOnProc.rehash(std::ceil(clientSize/globalScalarIndexOnProc.max_load_factor()));
478  for (int idx = 0; idx < clientSize; ++idx)
479  {
480    globalScalarIndexOnProc[idx].push_back(0);
481  }
482}
483
484/*!
485  Find out proc and global index of axis source which axis destination is on demande
486  \param[in] axisDst Axis destination
487  \param[in] axisSrc Axis source
488  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
489  \param[out] globalAxisIndexOnProc Global index of axis source on different procs
490*/
491void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst,
492                                                               CAxis* axisSrc,
493                                                               CArray<size_t,1>& destGlobalIndexPositionInGrid,
494                                                               boost::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc)
495{
496  CContext* context = CContext::getCurrent();
497  CContextClient* client=context->client;
498  int clientRank = client->clientRank;
499  int clientSize = client->clientSize;
500
501  size_t globalIndex;
502  int nIndexSize = axisSrc->index.numElements();
503  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
504  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
505  for (int idx = 0; idx < nIndexSize; ++idx)
506  {
507    globalIndex = axisSrc->index(idx);
508    globalIndex2ProcRank[globalIndex].push_back(clientRank);
509  }
510
511  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
512  CArray<size_t,1> globalAxisIndex(axisDst->index.numElements());
513  for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx)
514  {
515    globalAxisIndex(idx) = axisDst->index(idx);
516  }
517  dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);
518
519  std::vector<int> countIndex(clientSize,0);
520  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
521  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
522                                                               ite = computedGlobalIndexOnProc.end();
523  for (it = itb; it != ite; ++it)
524  {
525    const std::vector<int>& procList = it->second;
526    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
527  }
528
529  globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor()));
530  for (int idx = 0; idx < clientSize; ++idx)
531  {
532    if (0 != countIndex[idx])
533    {
534      globalAxisIndexOnProc[idx].resize(countIndex[idx]);
535      countIndex[idx] = 0;
536    }
537  }
538
539  for (it = itb; it != ite; ++it)
540  {
541    const std::vector<int>& procList = it->second;
542    for (int idx = 0; idx < procList.size(); ++idx)
543    {
544      globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
545      ++countIndex[procList[idx]];
546    }
547  }
548}
549
550/*!
551  Find out proc and global index of domain source which domain destination is on demande
552  \param[in] domainDst Domain destination
553  \param[in] domainSrc Domain source
554  \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid.
555  \param[out] globalDomainIndexOnProc Global index of domain source on different procs
556*/
557void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst,
558                                                                 CDomain* domainSrc,
559                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
560                                                                 boost::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc)
561{
562  CContext* context = CContext::getCurrent();
563  CContextClient* client=context->client;
564  int clientRank = client->clientRank;
565  int clientSize = client->clientSize;
566
567  int niGlobSrc = domainSrc->ni_glo.getValue();
568  size_t globalIndex;
569  int i_ind, j_ind;
570  int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements()
571                                                             : destGlobalIndexPositionInGrid.numElements();
572  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
573  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
574  if (destGlobalIndexPositionInGrid.isEmpty())
575  {
576    for (int idx = 0; idx < nIndexSize; ++idx)
577    {
578      i_ind=domainSrc->i_index(idx) ;
579      j_ind=domainSrc->j_index(idx) ;
580
581      globalIndex = i_ind + j_ind * niGlobSrc;
582      globalIndex2ProcRank[globalIndex].resize(1);
583      globalIndex2ProcRank[globalIndex][0] = clientRank;
584    }
585  }
586  else
587  {
588    for (int idx = 0; idx < nIndexSize; ++idx)
589    {
590      globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank);
591    }
592  }
593
594  CArray<size_t,1> globalDomainIndex;
595  if (destGlobalIndexPositionInGrid.isEmpty())
596  {
597    int niGlobDst = domainDst->ni_glo.getValue();
598    globalDomainIndex.resize(domainDst->i_index.numElements());
599    nIndexSize = domainDst->i_index.numElements();
600
601    for (int idx = 0; idx < nIndexSize; ++idx)
602    {
603      i_ind=domainDst->i_index(idx) ;
604      j_ind=domainDst->j_index(idx) ;
605
606      globalDomainIndex(idx) = i_ind + j_ind * niGlobDst;
607    }
608  }
609  else
610  {
611    globalDomainIndex.reference(destGlobalIndexPositionInGrid);
612  }
613
614  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
615  dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);
616
617  std::vector<int> countIndex(clientSize,0);
618  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
619  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
620                                                               ite = computedGlobalIndexOnProc.end();
621  for (it = itb; it != ite; ++it)
622  {
623    const std::vector<int>& procList = it->second;
624    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
625  }
626
627  globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor()));
628  for (int idx = 0; idx < clientSize; ++idx)
629  {
630    if (0 != countIndex[idx])
631    {
632      globalDomainIndexOnProc[idx].resize(countIndex[idx]);
633      countIndex[idx] = 0;
634    }
635  }
636
637  for (it = itb; it != ite; ++it)
638  {
639    const std::vector<int>& procList = it->second;
640    for (int idx = 0; idx < procList.size(); ++idx)
641    {
642      globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
643      ++countIndex[procList[idx]];
644    }
645  }
646}
647
648/*!
649  Compute index mapping between element source and element destination with an auxiliary inputs which determine
650position of each mapped index in global index of grid destination.
651  \param [in] dataAuxInputs auxiliary inputs
652*/
653void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs)
654{
655  computeIndexSourceMapping_(dataAuxInputs);
656}
657
658std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs()
659{
660  return idAuxInputs_;
661}
662
663CGenericAlgorithmTransformation::AlgoTransType CGenericAlgorithmTransformation::type()
664{
665  return type_;
666}
667
668}
Note: See TracBrowser for help on using the repository browser.