source: XIOS/trunk/src/transformation/grid_transformation.cpp @ 687

Last change on this file since 687 was 687, checked in by mhnguyen, 7 years ago

Implementing auto-generate rectilinear domain

+) Add a new special transformation to generate (complete) rectilinear domain

Test
+) On Curie
+) test_new_feature passed

File size: 28.3 KB
Line 
1/*!
2   \file grid_transformation.cpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 02 Jul 2015
6
7   \brief Interface for all transformations.
8 */
9#include "grid_transformation.hpp"
10#include "axis_algorithm_inverse.hpp"
11#include "axis_algorithm_zoom.hpp"
12#include "axis_algorithm_interpolate.hpp"
13#include "domain_algorithm_zoom.hpp"
14#include "domain_algorithm_interpolate_from_file.hpp"
15#include "context.hpp"
16#include "context_client.hpp"
17#include "transformation_mapping.hpp"
18#include "axis_algorithm_transformation.hpp"
19#include "distribution_client.hpp"
20
21namespace xios {
22CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
23: gridSource_(source), gridDestination_(destination), originalGridSource_(source),
24  globalIndexOfCurrentGridSource_(), globalIndexOfOriginalGridSource_(), weightOfGlobalIndexOfOriginalGridSource_(0), algoTypes_()
25{
26  //Verify the compatibity between two grids
27  int numElement = gridDestination_->axis_domain_order.numElements();
28  if (numElement != gridSource_->axis_domain_order.numElements())
29    ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
30       << "Two grids have different number of elements"
31       << "Number of elements of grid source " <<gridSource_->getId() << " is " << gridSource_->axis_domain_order.numElements()  << std::endl
32       << "Number of elements of grid destination " <<gridDestination_->getId() << " is " << numElement);
33
34  for (int i = 0; i < numElement; ++i)
35  {
36    if (gridDestination_->axis_domain_order(i) != gridSource_->axis_domain_order(i))
37      ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
38         << "Transformed grid and its grid source have incompatible elements"
39         << "Grid source " <<gridSource_->getId() << std::endl
40         << "Grid destination " <<gridDestination_->getId());
41  }
42
43  std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
44  std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
45  for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
46  {
47    CAxis* axis = CAxis::createAxis();
48    axis->setAttributes(axisSrcTmp[idx]);
49    axisSrc.push_back(axis);
50  }
51
52  for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
53  {
54    CDomain* domain = CDomain::createDomain();
55    domain->setAttributes(domainSrcTmp[idx]);
56    domainSrc.push_back(domain);
57  }
58
59  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
60
61  initializeMappingOfOriginalGridSource();
62  initializeAlgorithms();
63}
64
65/*!
66  Initialize the mapping between the first grid source and the original one
67  In a series of transformation, for each step, there is a need to "create" a new grid that plays a role of "temporary" source.
68Because at the end of the series, we need to know about the index mapping between the final grid destination and original grid source,
69for each transformation, we need to make sure that the current "temporary source" maps its global index correctly to the original one.
70*/
71void CGridTransformation::initializeMappingOfOriginalGridSource()
72{
73  CContext* context = CContext::getCurrent();
74  CContextClient* client = context->client;
75
76  CDistributionClient distribution(client->clientRank, originalGridSource_);
77  const std::vector<size_t>& globalIndexGridSrcSendToServer = distribution.getGlobalDataIndexSendToServer();
78
79  weightOfGlobalIndexOfOriginalGridSource_.resize(globalIndexGridSrcSendToServer.size());
80  globalIndexOfCurrentGridSource_  = globalIndexGridSrcSendToServer;
81  globalIndexOfOriginalGridSource_ = globalIndexGridSrcSendToServer;
82  weightOfGlobalIndexOfOriginalGridSource_ = 1.0;
83}
84
85CGridTransformation::~CGridTransformation()
86{
87  std::list<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it,
88                                                              ite = algoTransformation_.end();
89  for (it = itb; it != ite; ++it) delete (*it);
90}
91
92/*!
93  Initialize the algorithms (transformations)
94*/
95void CGridTransformation::initializeAlgorithms()
96{
97  std::vector<int> axisPositionInGrid;
98  std::vector<int> domPositionInGrid;
99  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
100  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
101
102  int idx = 0;
103  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
104  {
105    if (false == (gridDestination_->axis_domain_order)(i))
106    {
107      axisPositionInGrid.push_back(idx);
108      ++idx;
109    }
110    else
111    {
112      ++idx;
113      domPositionInGrid.push_back(idx);
114      ++idx;
115    }
116  }
117
118  for (int i = 0; i < axisListDestP.size(); ++i)
119  {
120    elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
121  }
122
123  for (int i = 0; i < domListDestP.size(); ++i)
124  {
125    elementPosition2DomainPositionInGrid_[domPositionInGrid[i]] = i;
126  }
127
128  idx = 0;
129  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
130  {
131    if (false == (gridDestination_->axis_domain_order)(i))
132    {
133      initializeAxisAlgorithms(idx);
134      ++idx;
135    }
136    else
137    {
138      ++idx;
139      initializeDomainAlgorithms(idx);
140      ++idx;
141    }
142  }
143}
144
145
146
147/*!
148  Initialize the algorithms corresponding to transformation info contained in each axis.
149If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
150In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
151For now, one approach is to do these combinely but maybe this needs changing.
152\param [in] axisPositionInGrid position of an axis in grid. (for example: a grid with one domain and one axis, position of domain is 1, position of axis is 2)
153*/
154void CGridTransformation::initializeAxisAlgorithms(int axisPositionInGrid)
155{
156  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
157  if (!axisListDestP.empty())
158  {
159    if (axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->hasTransformation())
160    {
161      CAxis::TransMapTypes trans = axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->getAllTransformations();
162      CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
163                                           ite = trans.end();
164      int transformationOrder = 0;
165      for (it = itb; it != ite; ++it)
166      {
167        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
168        algoTypes_.push_back(false);
169        ++transformationOrder;
170      }
171    }
172  }
173}
174
175/*!
176  Initialize the algorithms corresponding to transformation info contained in each domain.
177If a domain has transformations, they will be represented in form of vector of CTransformation pointers
178In general, each domain can have several transformations performed on itself.
179\param [in] domPositionInGrid position of a domain in grid. (for example: a grid with one domain and one axis, position of domain is 1, position of axis is 2)
180*/
181void CGridTransformation::initializeDomainAlgorithms(int domPositionInGrid)
182{
183  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
184  if (!domListDestP.empty())
185  {
186    if (domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->hasTransformation())
187    {
188      CDomain::TransMapTypes trans = domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->getAllTransformations();
189      CDomain::TransMapTypes::const_iterator itb = trans.begin(), it,
190                                             ite = trans.end();
191      int transformationOrder = 0;
192      for (it = itb; it != ite; ++it)
193      {
194        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
195        algoTypes_.push_back(true);
196        ++transformationOrder;
197      }
198    }
199  }
200
201}
202
203/*!
204  Select algorithm correspoding to its transformation type and its position in each element
205  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
206                                             and position of axis is 2
207  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
208  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
209  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
210*/
211void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo)
212{
213   if (isDomainAlgo) selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
214   else selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
215}
216
217/*!
218  Select algorithm of an axis correspoding to its transformation type and its position in each element
219  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
220                                             and position of axis is 2
221  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
222  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
223*/
224void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
225{
226  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
227  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
228
229  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
230  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
231  CAxis::TransMapTypes::const_iterator it = trans.begin();
232
233  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
234
235  CZoomAxis* zoomAxis = 0;
236  CInterpolateAxis* interpAxis = 0;
237  CGenericAlgorithmTransformation* algo = 0;
238  switch (transType)
239  {
240    case TRANS_INTERPOLATE_AXIS:
241      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
242      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
243      break;
244    case TRANS_ZOOM_AXIS:
245      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
246      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
247      break;
248    case TRANS_INVERSE_AXIS:
249      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
250      break;
251    default:
252      break;
253  }
254  algoTransformation_.push_back(algo);
255
256}
257
258/*!
259  Select algorithm of a domain correspoding to its transformation type and its position in each element
260  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
261                                             and position of axis is 2
262  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
263  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
264*/
265void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
266{
267  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
268  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
269
270  int domainIndex =  elementPosition2DomainPositionInGrid_[elementPositionInGrid];
271  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
272  CDomain::TransMapTypes::const_iterator it = trans.begin();
273
274  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
275
276  CZoomDomain* zoomDomain = 0;
277  CInterpolateFromFileDomain* interpFileDomain = 0;
278  CGenericAlgorithmTransformation* algo = 0;
279  switch (transType)
280  {
281    case TRANS_INTERPOLATE_DOMAIN_FROM_FILE:
282      interpFileDomain = dynamic_cast<CInterpolateFromFileDomain*> (it->second);
283      algo = new CDomainAlgorithmInterpolateFromFile(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
284      break;
285    case TRANS_ZOOM_DOMAIN:
286      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
287      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
288      break;
289    default:
290      break;
291  }
292  algoTransformation_.push_back(algo);
293}
294
295/*!
296  Assign the current grid destination to the grid source in the new transformation.
297The current grid destination plays the role of grid source in next transformation (if any).
298Only element on which the transformation is performed is modified
299  \param [in] elementPositionInGrid position of element in grid
300  \param [in] transType transformation type
301*/
302void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType)
303{
304  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
305  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
306
307  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
308  std::vector<CDomain*> domListSrcP = gridSource_->getDomains();
309
310  int axisIndex, domainIndex;
311  switch (transType)
312  {
313    case TRANS_INTERPOLATE_DOMAIN_FROM_FILE:
314    case TRANS_ZOOM_DOMAIN:
315      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
316      domListSrcP[domainIndex]->duplicateAttributes(domListDestP[domainIndex]);
317      break;
318
319    case TRANS_INTERPOLATE_AXIS:
320    case TRANS_ZOOM_AXIS:
321    case TRANS_INVERSE_AXIS:
322      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
323      axisListSrcP[axisIndex]->duplicateAttributes(axisListDestP[axisIndex]);
324      break;
325    default:
326      break;
327  }
328}
329
330/*!
331  Perform all transformations
332  For each transformation, there are some things to do:
333  -) Chose the correct algorithm by transformation type and position of element
334  -) Calculate the mapping of global index between the current grid source and grid destination
335  -) Calculate the mapping of global index between current grid DESTINATION and ORIGINAL grid SOURCE
336  -) Make current grid destination become grid source in the next transformation
337*/
338void CGridTransformation::computeAll()
339{
340  CContext* context = CContext::getCurrent();
341  CContextClient* client = context->client;
342
343  ListAlgoType::const_iterator itb = listAlgos_.begin(),
344                               ite = listAlgos_.end(), it;
345  CGenericAlgorithmTransformation* algo = 0;
346  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
347  for (it = itb; it != ite; ++it)
348  {
349    int elementPositionInGrid = it->first;
350    ETranformationType transType = (it->second).first;
351    int transformationOrder = (it->second).second;
352    std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource;
353
354    // First of all, select an algorithm
355    selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
356    algo = algoTransformation_.back();
357
358    if (0 != algo) // Only registered transformation can be executed
359    {
360      // Recalculate the distribution of grid destination
361      CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
362      const std::vector<size_t>& globalIndexGridDestSendToServer = distributionClientDest.getGlobalDataIndexSendToServer();
363
364      // ComputeTransformation of global index of each element
365      std::vector<int> gridDestinationDimensionSize = gridDestination_->getGlobalDimension();
366      std::vector<int> gridSrcDimensionSize = gridSource_->getGlobalDimension();
367      int elementPosition = it->first;
368      algo->computeGlobalSourceIndex(elementPosition,
369                                     gridDestinationDimensionSize,
370                                     gridSrcDimensionSize,
371                                     globalIndexGridDestSendToServer,
372                                     globaIndexWeightFromDestToSource);
373
374      // Compute transformation of global indexes among grids
375      computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource);
376
377      // Now grid destination becomes grid source in a new transformation
378      setUpGrid(elementPositionInGrid, transType);
379      ++nbAgloTransformation;
380    }
381  }
382
383  if (0 != nbAgloTransformation)
384  {
385    updateFinalGridDestination();
386    computeFinalTransformationMapping();
387  }
388}
389
390
391/*!
392  After applying the algorithms, there are some informations on grid destination needing change, for now, there are:
393   +) mask
394*/
395void CGridTransformation::updateFinalGridDestination()
396{
397  CContext* context = CContext::getCurrent();
398  CContextClient* client = context->client;
399
400  //First of all, retrieve info of local mask of grid destination
401  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
402  const std::vector<int>& localMaskIndexOnClientDest = distributionClientDest.getLocalMaskIndexOnClient();
403  const std::vector<size_t>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer();
404
405  std::vector<size_t>::const_iterator itbArr, itArr, iteArr;
406  itbArr = globalIndexOnClientDest.begin();
407  iteArr = globalIndexOnClientDest.end();
408
409  // Then find out which index became invalid (become masked after being applied the algorithms, or demande some masked points from grid source)
410  int num = globalIndexOfOriginalGridSource_.size();
411  const size_t sfmax = NumTraits<unsigned long>::sfmax();
412  int maskIndexNum = 0;
413  for (int idx = 0; idx < num; ++idx)
414  {
415    if (sfmax == globalIndexOfOriginalGridSource_[idx])
416    {
417      size_t maskedGlobalIndex = globalIndexOfCurrentGridSource_[idx];
418      itArr = std::find(itbArr, iteArr, maskedGlobalIndex);
419      if (iteArr != itArr) ++maskIndexNum;
420    }
421  }
422
423  CArray<int,1> maskIndexToModify(maskIndexNum);
424  maskIndexNum = 0;
425  for (int idx = 0; idx < num; ++idx)
426  {
427    if (sfmax == globalIndexOfOriginalGridSource_[idx])
428    {
429      size_t maskedGlobalIndex = globalIndexOfCurrentGridSource_[idx];
430      itArr = std::find(itbArr, iteArr, maskedGlobalIndex);
431      if (iteArr != itArr)
432      {
433        int localIdx = std::distance(itbArr, itArr);
434        maskIndexToModify(maskIndexNum) = localMaskIndexOnClientDest[localIdx];
435        ++maskIndexNum;
436      }
437    }
438  }
439
440  gridDestination_->modifyMask(maskIndexToModify);
441}
442
443/*!
444  A transformation from a grid source to grid destination often passes several intermediate grids, which play a role of
445temporary grid source and/or grid destination. This function makes sure that global index of original grid source are mapped correctly to
446the final grid destination
447*/
448void CGridTransformation::computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource)
449{
450  CContext* context = CContext::getCurrent();
451  CContextClient* client = context->client;
452
453  CTransformationMapping transformationMap(gridDestination_, gridSource_);
454
455    // Then compute transformation mapping among clients
456  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
457
458  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
459  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
460
461 // Sending global index of original grid source
462  std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend,
463                                                     iteSend = globalIndexToSend.end();
464  std::vector<size_t>::const_iterator itbArr = globalIndexOfCurrentGridSource_.begin(), itArr,
465                                      iteArr = globalIndexOfCurrentGridSource_.end();
466 int sendBuffSize = 0;
467 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
468
469 typedef unsigned long Scalar;
470 unsigned long* sendBuff, *currentSendBuff;
471 if (0 != sendBuffSize) sendBuff = new unsigned long [sendBuffSize];
472 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
473
474 std::map<int, MPI_Request> requests;
475
476 std::vector<int> permutIndex(globalIndexOfCurrentGridSource_.size());
477 typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
478 XIOSAlgorithms::fillInIndex(globalIndexOfCurrentGridSource_.size(), permutIndex);
479 XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOfCurrentGridSource_, permutIndex);
480 BinarySearch searchCurrentSrc(globalIndexOfCurrentGridSource_);
481 std::vector<int>::iterator itbIndex = permutIndex.begin(), itIndex,
482                            iteIndex = permutIndex.end();
483
484  // Find out local index on grid destination (received)
485 int currentBuffPosition = 0;
486 for (itSend = itbSend; itSend != iteSend; ++itSend)
487 {
488   int destRank = itSend->first;
489   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
490   int countSize = globalIndexOfCurrentGridSourceToSend.size();
491   for (int idx = 0; idx < (countSize); ++idx)
492   {
493     if (searchCurrentSrc.search(itbIndex, iteIndex, globalIndexOfCurrentGridSourceToSend[idx], itIndex))
494     {
495//       int index = std::distance(itbArr, itArr);
496       sendBuff[idx+currentBuffPosition] = globalIndexOfOriginalGridSource_[*itIndex]; //[index];
497     }
498   }
499   currentSendBuff = sendBuff + currentBuffPosition;
500   MPI_Isend(currentSendBuff, countSize, MPI_UNSIGNED_LONG, destRank, 14, client->intraComm, &requests[destRank]);
501   currentBuffPosition += countSize;
502 }
503
504 // Receiving global index of grid source sending from current grid source
505 std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
506                                                                                     iteRecv = globalIndexToReceive.end();
507 int recvBuffSize = 0;
508 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
509
510 unsigned long* recvBuff, *currentRecvBuff;
511 if (0 != recvBuffSize) recvBuff = new unsigned long [recvBuffSize];
512 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
513
514 int currentRecvBuffPosition = 0;
515 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
516 {
517   MPI_Status status;
518   int srcRank = itRecv->first;
519   int countSize = (itRecv->second).size();
520   currentRecvBuff = recvBuff + currentRecvBuffPosition;
521   MPI_Recv(currentRecvBuff, countSize, MPI_UNSIGNED_LONG, srcRank, 14, client->intraComm, &status);
522   currentRecvBuffPosition += countSize;
523 }
524
525 int nbCurrentGridSource = 0;
526 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
527 {
528   int ssize = (itRecv->second).size();
529   for (int idx = 0; idx < ssize; ++idx)
530   {
531     nbCurrentGridSource += (itRecv->second)[idx].size();
532   }
533 }
534
535 if (globalIndexOfCurrentGridSource_.size() != nbCurrentGridSource)
536 {
537   globalIndexOfCurrentGridSource_.resize(nbCurrentGridSource);
538   globalIndexOfOriginalGridSource_.resize(nbCurrentGridSource);
539   weightOfGlobalIndexOfOriginalGridSource_.resize(nbCurrentGridSource);
540 }
541
542 int k = 0;
543 currentRecvBuff = recvBuff;
544 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
545 {
546   int countSize = (itRecv->second).size();
547   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff)
548   {
549     int ssize = (itRecv->second)[idx].size();
550     for (int i = 0; i < ssize; ++i)
551     {
552       globalIndexOfCurrentGridSource_[k] = ((itRecv->second)[idx][i]).first;
553       weightOfGlobalIndexOfOriginalGridSource_(k) = ((itRecv->second)[idx][i]).second;
554       globalIndexOfOriginalGridSource_[k] = *currentRecvBuff;
555       ++k;
556     }
557   }
558 }
559
560 std::map<int, MPI_Request>::iterator itRequest;
561 for (itRequest = requests.begin(); itRequest != requests.end(); ++itRequest)
562   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
563
564 if (0 != sendBuffSize) delete [] sendBuff;
565 if (0 != recvBuffSize) delete [] recvBuff;
566}
567
568/*!
569  Compute transformation mapping between grid source and grid destination
570  The transformation between grid source and grid destination is represented in form of mapping between global index
571of two grids. Then local index mapping between data on each grid will be found out thanks to these global indexes
572*/
573void CGridTransformation::computeFinalTransformationMapping()
574{
575  CContext* context = CContext::getCurrent();
576  CContextClient* client = context->client;
577
578  CTransformationMapping transformationMap(gridDestination_, originalGridSource_);
579
580  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource;
581  int nb = globalIndexOfCurrentGridSource_.size();
582  const size_t sfmax = NumTraits<unsigned long>::sfmax();
583  for (int idx = 0; idx < nb; ++idx)
584  {
585    if (sfmax != globalIndexOfOriginalGridSource_[idx])
586      globaIndexWeightFromDestToSource[globalIndexOfCurrentGridSource_[idx]].push_back(make_pair(globalIndexOfOriginalGridSource_[idx], weightOfGlobalIndexOfOriginalGridSource_(idx))) ;
587  }
588
589  // Then compute transformation mapping among clients
590  transformationMap.computeTransformationMapping(globaIndexWeightFromDestToSource);
591
592  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
593  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
594
595  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
596  CDistributionClient distributionClientSrc(client->clientRank, originalGridSource_);
597
598  const std::vector<size_t>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer();
599  const std::vector<size_t>& globalIndexOnClientSrc = distributionClientSrc.getGlobalDataIndexSendToServer();
600
601  std::vector<size_t>::const_iterator itbArr, itArr, iteArr;
602  std::vector<int>::const_iterator itIndex, itbIndex, iteIndex;
603  std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
604
605  std::vector<int> permutIndex;
606  typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
607
608  // Find out local index on grid destination (received)
609  XIOSAlgorithms::fillInIndex(globalIndexOnClientDest.size(), permutIndex);
610  XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOnClientDest, permutIndex);
611  itbIndex = permutIndex.begin();
612  iteIndex = permutIndex.end();
613  BinarySearch searchClientDest(globalIndexOnClientDest);
614  itbMapRecv = globalIndexToReceive.begin();
615  iteMapRecv = globalIndexToReceive.end();
616  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
617  {
618    int sourceRank = itMapRecv->first;
619    int numGlobalIndex = (itMapRecv->second).size();
620    for (int i = 0; i < numGlobalIndex; ++i)
621    {
622      int vecSize = ((itMapRecv->second)[i]).size();
623      std::vector<std::pair<int,double> > tmpVec;
624      for (int idx = 0; idx < vecSize; ++idx)
625      {
626        size_t globalIndex = (itMapRecv->second)[i][idx].first;
627        double weight = (itMapRecv->second)[i][idx].second;
628        if (searchClientDest.search(itbIndex, iteIndex, globalIndex, itIndex))
629        {
630          tmpVec.push_back(make_pair(*itIndex, weight));
631        }
632      }
633      localIndexToReceiveOnGridDest_[sourceRank].push_back(tmpVec);
634    }
635  }
636
637  // Find out local index on grid source (to send)
638  std::map<int,std::vector<size_t> >::const_iterator itbMap, itMap, iteMap;
639  XIOSAlgorithms::fillInIndex(globalIndexOnClientSrc.size(), permutIndex);
640  XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOnClientSrc, permutIndex);
641  itbIndex = permutIndex.begin();
642  iteIndex = permutIndex.end();
643  BinarySearch searchClientSrc(globalIndexOnClientSrc);
644  itbMap = globalIndexToSend.begin();
645  iteMap = globalIndexToSend.end();
646  for (itMap = itbMap; itMap != iteMap; ++itMap)
647  {
648    int destRank = itMap->first;
649    int vecSize = itMap->second.size();
650    localIndexToSendFromGridSource_[destRank].resize(vecSize);
651    for (int idx = 0; idx < vecSize; ++idx)
652    {
653      if (searchClientSrc.search(itbIndex, iteIndex, itMap->second[idx], itIndex))
654      {
655        localIndexToSendFromGridSource_[destRank](idx) = *itIndex;
656      }
657    }
658  }
659}
660
661/*!
662  Local index of data which need sending from the grid source
663  \return local index of data
664*/
665const std::map<int, CArray<int,1> >& CGridTransformation::getLocalIndexToSendFromGridSource() const
666{
667  return localIndexToSendFromGridSource_;
668}
669
670/*!
671  Local index of data which will be received on the grid destination
672  \return local index of data
673*/
674const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
675{
676  return localIndexToReceiveOnGridDest_;
677}
678
679}
Note: See TracBrowser for help on using the repository browser.