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

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

Implementing some code factoring

+) Replace some slow searching function by faster ones

Test
+) On Curie
+) test_client and test_complete are correct

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