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

Last change on this file since 871 was 871, checked in by mhnguyen, 5 years ago

Correcting a bug which sometimes make transformations between different size grids incorrect

+) Correct inputs of some functions
+) Add some auxilliary functions.
+) Add new test cases for test_remap

Test
+) Basic test pass
+) Single horizontal and vertical interpolation are correct
+) Chained interpolation are not correct

  • Property svn:executable set to *
File size: 31.8 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.hpp"
15#include "context.hpp"
16#include "context_client.hpp"
17#include "axis_algorithm_transformation.hpp"
18#include "distribution_client.hpp"
19#include "mpi_tag.hpp"
20#include "grid.hpp"
21#include <boost/unordered_map.hpp>
22
23namespace xios {
24CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
25: gridSource_(source), gridDestination_(destination), tmpGridDestination_(destination),
26  originalGridSource_(source), algoTypes_(), nbAlgos_(0), tempGridSrcs_(), tempGridDests_(),
27  auxInputs_(), dynamicalTransformation_(false), timeStamp_()
28
29{
30  //Verify the compatibity between two grids
31  int numElement = gridDestination_->axis_domain_order.numElements();
32  if (numElement != gridSource_->axis_domain_order.numElements())
33    ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
34       << "Two grids have different number of elements"
35       << "Number of elements of grid source " <<gridSource_->getId() << " is " << gridSource_->axis_domain_order.numElements()  << std::endl
36       << "Number of elements of grid destination " <<gridDestination_->getId() << " is " << numElement);
37
38  for (int i = 0; i < numElement; ++i)
39  {
40    if (gridDestination_->axis_domain_order(i) != gridSource_->axis_domain_order(i))
41      ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
42         << "Transformed grid and its grid source have incompatible elements"
43         << "Grid source " <<gridSource_->getId() << std::endl
44         << "Grid destination " <<gridDestination_->getId());
45  }
46
47  initializeTransformations();
48}
49
50/*!
51  Initialize the mapping between the first grid source and the original one
52  In a series of transformation, for each step, there is a need to "create" a new grid that plays a role of "temporary" source.
53Because at the end of the series, we need to know about the index mapping between the final grid destination and original grid source,
54for each transformation, we need to make sure that the current "temporary source" maps its global index correctly to the original one.
55*/
56void CGridTransformation::initializeTransformations()
57{
58  CContext* context = CContext::getCurrent();
59  CContextClient* client = context->client;
60
61  // Initialize algorithms
62  initializeAlgorithms();
63
64  ListAlgoType::const_iterator itb = listAlgos_.begin(),
65                               ite = listAlgos_.end(), it;
66
67  for (it = itb; it != ite; ++it)
68  {
69    ETranformationType transType = (it->second).first;
70    if (!isSpecialTransformation(transType)) ++nbAlgos_;
71  }
72}
73
74CGridTransformation::~CGridTransformation()
75{
76  std::vector<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it,
77                                                                ite = algoTransformation_.end();
78  for (it = itb; it != ite; ++it) delete (*it);
79}
80
81/*!
82  Initialize the algorithms (transformations)
83*/
84void CGridTransformation::initializeAlgorithms()
85{
86  std::vector<int> axisPositionInGrid;
87  std::vector<int> domPositionInGrid;
88  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
89  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
90
91  int idx = 0;
92  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
93  {
94    if (false == (gridDestination_->axis_domain_order)(i))
95    {
96      axisPositionInGrid.push_back(i);
97    }
98    else
99    {
100      domPositionInGrid.push_back(i);
101    }
102  }
103
104  for (int i = 0; i < axisListDestP.size(); ++i)
105  {
106    elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
107  }
108
109  for (int i = 0; i < domListDestP.size(); ++i)
110  {
111    elementPosition2DomainPositionInGrid_[domPositionInGrid[i]] = i;
112  }
113
114  idx = 0;
115  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
116  {
117    if (false == (gridDestination_->axis_domain_order)(i))
118    {
119      initializeAxisAlgorithms(i);
120    }
121    else
122    {
123      initializeDomainAlgorithms(i);
124    }
125  }
126}
127
128/*!
129  Initialize the algorithms corresponding to transformation info contained in each axis.
130If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
131In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
132For now, one approach is to do these combinely but maybe this needs changing.
133\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)
134*/
135void CGridTransformation::initializeAxisAlgorithms(int axisPositionInGrid)
136{
137  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
138  if (!axisListDestP.empty())
139  {
140    if (axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->hasTransformation())
141    {
142      CAxis::TransMapTypes trans = axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->getAllTransformations();
143      CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
144                                           ite = trans.end();
145      int transformationOrder = 0;
146      for (it = itb; it != ite; ++it)
147      {
148        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
149        algoTypes_.push_back(false);
150        ++transformationOrder;
151        std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
152        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]);
153      }
154    }
155  }
156}
157
158/*!
159  Initialize the algorithms corresponding to transformation info contained in each domain.
160If a domain has transformations, they will be represented in form of vector of CTransformation pointers
161In general, each domain can have several transformations performed on itself.
162\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)
163*/
164void CGridTransformation::initializeDomainAlgorithms(int domPositionInGrid)
165{
166  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
167  if (!domListDestP.empty())
168  {
169    if (domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->hasTransformation())
170    {
171      CDomain::TransMapTypes trans = domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->getAllTransformations();
172      CDomain::TransMapTypes::const_iterator itb = trans.begin(), it,
173                                             ite = trans.end();
174      int transformationOrder = 0;
175      for (it = itb; it != ite; ++it)
176      {
177        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
178        algoTypes_.push_back(true);
179        ++transformationOrder;
180        std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
181        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]);
182      }
183    }
184  }
185
186}
187
188/*!
189  Select algorithm correspoding to its transformation type and its position in each element
190  \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)
191                                             and position of axis is 2
192  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
193  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
194  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
195*/
196void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo)
197{
198   if (isDomainAlgo) selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
199   else selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
200}
201
202/*!
203  Select algorithm of an axis correspoding to its transformation type and its position in each element
204  \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)
205                                             and position of axis is 2
206  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
207  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
208*/
209void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
210{
211  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
212  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
213
214  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
215  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
216  CAxis::TransMapTypes::const_iterator it = trans.begin();
217
218  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
219
220  CZoomAxis* zoomAxis = 0;
221  CInterpolateAxis* interpAxis = 0;
222  CGenericAlgorithmTransformation* algo = 0;
223  switch (transType)
224  {
225    case TRANS_INTERPOLATE_AXIS:
226      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
227      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
228      break;
229    case TRANS_ZOOM_AXIS:
230      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
231      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
232      break;
233    case TRANS_INVERSE_AXIS:
234      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
235      break;
236    default:
237      break;
238  }
239  algoTransformation_.push_back(algo);
240
241}
242
243/*!
244  Select algorithm of a domain correspoding to its transformation type and its position in each element
245  \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)
246                                             and position of axis is 2
247  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
248  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
249*/
250void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
251{
252  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
253  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
254
255  int domainIndex =  elementPosition2DomainPositionInGrid_[elementPositionInGrid];
256  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
257  CDomain::TransMapTypes::const_iterator it = trans.begin();
258
259  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
260
261  CZoomDomain* zoomDomain = 0;
262  CInterpolateDomain* interpFileDomain = 0;
263  CGenericAlgorithmTransformation* algo = 0;
264  switch (transType)
265  {
266    case TRANS_INTERPOLATE_DOMAIN:
267      interpFileDomain = dynamic_cast<CInterpolateDomain*> (it->second);
268      algo = new CDomainAlgorithmInterpolate(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
269      break;
270    case TRANS_ZOOM_DOMAIN:
271      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
272      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
273      break;
274    default:
275      break;
276  }
277  algoTransformation_.push_back(algo);
278}
279
280/*!
281  If there are more than one transformation, a new temporary grid will be created and it will play the role of grid destination.
282This new created one keeps a pointer to the real transformed element of grid destination and generate new copies of other elements from grid source.
283  \param [in] elementPositionInGrid position of element in grid
284  \param [in] transType transformation type
285*/
286void CGridTransformation::setUpGridDestination(int elementPositionInGrid, ETranformationType transType, int nbTransformation)
287{
288  if (isSpecialTransformation(transType)) return;
289
290  if (!tempGridDests_.empty() && (getNbAlgo() == tempGridDests_.size()))
291  {
292    tempGridDests_.resize(0);
293//    tmpGridDestination_ = tempGridDests_[nbTransformation];
294//    return;
295  }
296  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
297  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisDst;
298
299  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
300  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainDst;
301
302  int axisIndex = -1, domainIndex = -1;
303  switch (transType)
304  {
305    case TRANS_INTERPOLATE_DOMAIN:
306    case TRANS_ZOOM_DOMAIN:
307      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
308      break;
309
310    case TRANS_INTERPOLATE_AXIS:
311    case TRANS_ZOOM_AXIS:
312    case TRANS_INVERSE_AXIS:
313      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
314      break;
315    default:
316      break;
317  }
318
319  for (int idx = 0; idx < axisListSrcP.size(); ++idx)
320  {
321    CAxis* axis = (axisIndex != idx) ? axisListSrcP[idx] : axisListDestP[idx];
322    axisDst.push_back(axis);
323  }
324
325  for (int idx = 0; idx < domListSrcP.size(); ++idx)
326  {
327    CDomain* domain = (domainIndex != idx) ? domListSrcP[idx] : domListDestP[idx];
328    domainDst.push_back(domain);
329  }
330
331  tmpGridDestination_ = CGrid::createGrid(domainDst, axisDst, gridDestination_->axis_domain_order);
332  tmpGridDestination_->computeGridGlobalDimension(domainDst, axisDst, gridDestination_->axis_domain_order);
333  tempGridDests_.push_back(tmpGridDestination_);
334}
335
336/*!
337  Assign the current grid destination to the grid source in the new transformation.
338The current grid destination plays the role of grid source in next transformation (if any).
339Only element on which the transformation is performed is modified
340  \param [in] elementPositionInGrid position of element in grid
341  \param [in] transType transformation type
342*/
343void CGridTransformation::setUpGridSource(int elementPositionInGrid, ETranformationType transType, int nbTransformation)
344{
345  if (!tempGridSrcs_.empty() && (getNbAlgo()-1) == tempGridSrcs_.size())
346  {
347    tempGridSrcs_.resize(0);
348//    gridSource_ = tempGridSrcs_[nbTransformation];
349//    return;
350  }
351
352  std::vector<CAxis*> axisListDestP = tmpGridDestination_->getAxis();
353  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
354
355  std::vector<CDomain*> domListDestP = tmpGridDestination_->getDomains();
356  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
357
358  int axisIndex = -1, domainIndex = -1;
359  switch (transType)
360  {
361    case TRANS_INTERPOLATE_DOMAIN:
362    case TRANS_ZOOM_DOMAIN:
363      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
364      break;
365
366    case TRANS_INTERPOLATE_AXIS:
367    case TRANS_ZOOM_AXIS:
368    case TRANS_INVERSE_AXIS:
369      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
370      break;
371    default:
372      break;
373  }
374
375  for (int idx = 0; idx < axisListSrcP.size(); ++idx)
376  {
377    if (axisIndex != idx)
378    {
379      CAxis* axis = axisListDestP[idx];
380      axisSrc.push_back(axis);
381    }
382    else
383    {
384      CAxis* axis = CAxis::createAxis();
385      axis->axis_ref.setValue(axisListDestP[idx]->getId());
386      axis->solveRefInheritance(true);
387      axis->checkAttributesOnClient();
388      axisSrc.push_back(axis);
389    }
390  }
391
392  for (int idx = 0; idx < domListSrcP.size(); ++idx)
393  {
394    if (domainIndex != idx)
395    {
396      CDomain* domain = domListDestP[idx];
397      domainSrc.push_back(domain);
398    }
399    else
400    {
401      CDomain* domain = CDomain::createDomain();
402      domain->domain_ref.setValue(domListDestP[idx]->getId());
403      domain->solveRefInheritance(true);
404      domain->checkAttributesOnClient();
405      domainSrc.push_back(domain);
406    }
407  }
408
409//  for (int idx = 0; idx < axisListSrcP.size(); ++idx)
410//  {
411//    CAxis* axis = CAxis::createAxis();
412//    if (axisIndex != idx) axis->axis_ref.setValue(axisListSrcP[idx]->getId());
413//    else axis->axis_ref.setValue(axisListDestP[idx]->getId());
414//    axis->solveRefInheritance(true);
415//    axis->checkAttributesOnClient();
416//    axisSrc.push_back(axis);
417//  }
418//
419//  for (int idx = 0; idx < domListSrcP.size(); ++idx)
420//  {
421//    CDomain* domain = CDomain::createDomain();
422//    if (domainIndex != idx) domain->domain_ref.setValue(domListSrcP[idx]->getId());
423//    else domain->domain_ref.setValue(domListDestP[idx]->getId());
424//    domain->solveRefInheritance(true);
425//    domain->checkAttributesOnClient();
426//    domainSrc.push_back(domain);
427//  }
428
429  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, tmpGridDestination_->axis_domain_order);
430  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, tmpGridDestination_->axis_domain_order);
431
432  tempGridSrcs_.push_back(gridSource_);
433}
434
435/*!
436  Perform all transformations
437  For each transformation, there are some things to do:
438  -) Chose the correct algorithm by transformation type and position of element
439  -) Calculate the mapping of global index between the current grid source and grid destination
440  -) Calculate the mapping of global index between current grid DESTINATION and grid SOURCE
441  -) Make current grid destination become grid source in the next transformation
442*/
443void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
444{
445  if (nbAlgos_ < 1) return;
446  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
447  if (dynamicalTransformation_)
448  {
449    if (timeStamp_.insert(timeStamp).second)   //Reset map
450    {
451      std::list<SendingIndexGridSourceMap>().swap(localIndexToSendFromGridSource_);
452      std::list<RecvIndexGridDestinationMap>().swap(localIndexToReceiveOnGridDest_);
453      std::list<size_t>().swap(nbLocalIndexOnGridDest_);
454    }
455    else
456      return;
457  }
458
459  CContext* context = CContext::getCurrent();
460  CContextClient* client = context->client;
461
462  ListAlgoType::const_iterator itb = listAlgos_.begin(),
463                               ite = listAlgos_.end(), it;
464
465  CGenericAlgorithmTransformation* algo = 0;
466  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
467  for (it = itb; it != ite; ++it)
468  {
469    int elementPositionInGrid = it->first;
470    ETranformationType transType = (it->second).first;
471    int transformationOrder = (it->second).second;
472    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
473
474//    if (1 < nbAlgos_)
475//    {
476//      // Create a temporary grid destination which contains transformed element of grid destination and
477//      // non-transformed elements fo grid source
478////      if (nbAgloTransformation != (nbAlgos_-1)) setUpGridDestination(elementPositionInGrid, transType, nbAgloTransformation);
479//    }
480    setUpGridDestination(elementPositionInGrid, transType, nbAgloTransformation);
481
482    // First of all, select an algorithm
483    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
484    {
485      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
486      algo = algoTransformation_.back();
487    }
488    else
489      algo = algoTransformation_[std::distance(itb, it)];
490
491    if (0 != algo) // Only registered transformation can be executed
492    {
493      algo->computeIndexSourceMapping(dataAuxInputs);
494
495      // ComputeTransformation of global index of each element
496      int elementPosition = it->first;
497      algo->computeGlobalSourceIndex(elementPosition,
498                                     gridSource_,
499                                     tmpGridDestination_,
500                                     globaIndexWeightFromSrcToDst);
501
502      // Compute transformation of global indexes among grids
503      computeTransformationMapping(globaIndexWeightFromSrcToDst);
504
505      if (1 < nbAlgos_)
506      {
507        // Now grid destination becomes grid source in a new transformation
508        if (nbAgloTransformation != (nbAlgos_-1)) setUpGridSource(elementPositionInGrid, transType, nbAgloTransformation);
509      }
510      ++nbAgloTransformation;
511    }
512  }
513}
514
515/*!
516  Compute exchange index between grid source and grid destination
517  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
518*/
519void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
520{
521  CContext* context = CContext::getCurrent();
522  CContextClient* client = context->client;
523  int nbClient = client->clientSize;
524  int clientRank = client->clientRank;
525
526  // Recalculate the distribution of grid destination
527  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
528//  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
529  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
530  const std::vector<int>& localMask = distributionClientDest.getLocalMaskIndexOnClient();
531
532  // Update number of local index on each transformation
533  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
534  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
535  localMaskOnGridDest_.push_back(std::vector<bool>());
536  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
537  tmpMask.resize(nbLocalIndex,false);
538//  for (int idx = 0; idx < nbLocalIndex; ++idx)
539//  {
540//    tmpMask[localMask[idx]] = true;
541//  }
542
543  // Find out number of index sent from grid source and number of index received on grid destination
544  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
545                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
546  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
547  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
548  int connectedClient = globaIndexWeightFromSrcToDst.size();
549  int* recvCount=new int[nbClient];
550  int* displ=new int[nbClient];
551  int* sendRankBuff=new int[connectedClient];
552  int* sendSizeBuff=new int[connectedClient];
553  int n = 0;
554  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
555  {
556    sendRankBuff[n] = itIndex->first;
557    const SendIndexMap& sendIndexMap = itIndex->second;
558    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
559    int sendSize = 0;
560    for (itSend = itbSend; itSend != iteSend; ++itSend)
561    {
562      sendSize += itSend->second.size();
563    }
564    sendSizeBuff[n] = sendSize;
565    sendRankSizeMap[itIndex->first] = sendSize;
566  }
567  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
568
569  displ[0]=0 ;
570  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
571  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
572  int* recvRankBuff=new int[recvSize];
573  int* recvSizeBuff=new int[recvSize];
574  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
575  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
576  for (int i = 0; i < nbClient; ++i)
577  {
578    int currentPos = displ[i];
579    for (int j = 0; j < recvCount[i]; ++j)
580      if (recvRankBuff[currentPos+j] == clientRank)
581      {
582        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
583      }
584  }
585
586  // Sending global index of grid source to corresponding process as well as the corresponding mask
587  std::vector<MPI_Request> requests;
588  std::vector<MPI_Status> status;
589  boost::unordered_map<int, unsigned char* > recvMaskDst;
590  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
591  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
592  {
593    int recvRank = itRecv->first;
594    int recvSize = itRecv->second;
595    recvMaskDst[recvRank] = new unsigned char [recvSize];
596    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
597
598    requests.push_back(MPI_Request());
599    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
600    requests.push_back(MPI_Request());
601    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
602  }
603
604  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
605  boost::unordered_map<int, CArray<double,1> > weightDst;
606  boost::unordered_map<int, unsigned char* > sendMaskDst;
607  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
608  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
609  {
610    int sendRank = itIndex->first;
611    int sendSize = sendRankSizeMap[sendRank];
612    const SendIndexMap& sendIndexMap = itIndex->second;
613    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
614    globalIndexDst[sendRank].resize(sendSize);
615    weightDst[sendRank].resize(sendSize);
616    sendMaskDst[sendRank] = new unsigned char [sendSize];
617    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
618    int countIndex = 0;
619    for (itSend = itbSend; itSend != iteSend; ++itSend)
620    {
621      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
622      for (int idx = 0; idx < dstWeight.size(); ++idx)
623      {
624        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
625        weightDst[sendRank](countIndex) = dstWeight[idx].second;
626        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
627          sendMaskDst[sendRank][countIndex] = 1;
628        else
629          sendMaskDst[sendRank][countIndex] = 0;
630        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
631        ++countIndex;
632      }
633    }
634
635    // Send global index source and mask
636    requests.push_back(MPI_Request());
637    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
638    requests.push_back(MPI_Request());
639    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
640  }
641
642  status.resize(requests.size());
643  MPI_Waitall(requests.size(), &requests[0], &status[0]);
644
645  // Okie, now use the mask to identify which index source we need to send, then also signal the destination which masked index we will return
646  std::vector<MPI_Request>().swap(requests);
647  std::vector<MPI_Status>().swap(status);
648  // Okie, on destination side, we will wait for information of masked index of source
649  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
650  {
651    int recvRank = itSend->first;
652    int recvSize = itSend->second;
653
654    requests.push_back(MPI_Request());
655    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
656  }
657
658  // Ok, now we fill in local index of grid source (we even count for masked index)
659  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
660  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
661  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
662  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
663  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
664  {
665    int recvRank = itRecv->first;
666    int recvSize = itRecv->second;
667    unsigned char* recvMask = recvMaskDst[recvRank];
668    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
669    int realSendSize = 0;
670    for (int idx = 0; idx < recvSize; ++idx)
671    {
672      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
673        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
674         ++realSendSize;
675        else // inform the destination that this index is masked
676         *(recvMask+idx) = 0;
677    }
678
679    tmpSend[recvRank].resize(realSendSize);
680    realSendSize = 0;
681    for (int idx = 0; idx < recvSize; ++idx)
682    {
683      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
684      {
685        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
686         ++realSendSize;
687      }
688    }
689
690    // Okie, now inform the destination which source index are masked
691    requests.push_back(MPI_Request());
692    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
693  }
694  status.resize(requests.size());
695  MPI_Waitall(requests.size(), &requests[0], &status[0]);
696
697  // Cool, now we can fill in local index of grid destination (counted for masked index)
698  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
699  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
700  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
701  {
702    int recvRank = itSend->first;
703    int recvSize = itSend->second;
704    unsigned char* recvMask = sendMaskDst[recvRank];
705
706    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
707    CArray<double,1>& recvWeightDst = weightDst[recvRank];
708    int realRecvSize = 0;
709    for (int idx = 0; idx < recvSize; ++idx)
710    {
711      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
712         ++realRecvSize;
713    }
714
715    int localIndexDst;
716    recvTmp[recvRank].resize(realRecvSize);
717    realRecvSize = 0;
718    for (int idx = 0; idx < recvSize; ++idx)
719    {
720      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
721      {
722        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
723        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
724         ++realRecvSize;
725      }
726    }
727  }
728
729  delete [] recvCount;
730  delete [] displ;
731  delete [] sendRankBuff;
732  delete [] recvRankBuff;
733  delete [] sendSizeBuff;
734  delete [] recvSizeBuff;
735
736  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
737  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
738    delete [] itChar->second;
739  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
740    delete [] itChar->second;
741  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
742  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
743    delete [] itLong->second;
744  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
745    delete [] itLong->second;
746
747}
748
749bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
750{
751  bool res;
752  switch (transType)
753  {
754    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
755     res = true;
756     break;
757    default:
758     res = false;
759     break;
760  }
761
762  return res;
763}
764
765/*!
766  Local index of data which need sending from the grid source
767  \return local index of data
768*/
769const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
770{
771  return localIndexToSendFromGridSource_;
772}
773
774/*!
775  Local index of data which will be received on the grid destination
776  \return local index of data
777*/
778const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
779{
780  return localIndexToReceiveOnGridDest_;
781}
782
783/*!
784 Number of index will be received on the grid destination
785  \return number of index of data
786*/
787const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
788{
789  return nbLocalIndexOnGridDest_;
790}
791
792/*!
793  Local mask of data which will be received on the grid destination
794  \return local mask of data
795*/
796const std::list<std::vector<bool> >& CGridTransformation::getLocalMaskIndexOnGridDest() const
797{
798  return localMaskOnGridDest_;
799}
800
801}
Note: See TracBrowser for help on using the repository browser.