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

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

Adding interpolation test_remap

+) Add new test case for domain interpolation
+) Resolve circular dependence
+) Add function to read weigh file

Test
+) On Curie
+) test_remap can print out .nc

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