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

Last change on this file since 1261 was 1261, checked in by ymipsl, 4 years ago

backport of rev1260 in trunk.

YM

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