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

Last change on this file since 933 was 933, checked in by mhnguyen, 6 years ago

Improving transformation selection. Instead of modifying directly grid_transformation
we only need to register a new transformation with the framework

+) Update all transformations with this new method

Test
+) On Curie
+) Basic tests pass

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