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

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

Correcting a bug invovling with the order of transformations

+) Add a check to make sure that, if the grid destination and grid source share a
same element, all transformations associated with this elements will not be accounted
in transformations list

Test
+) On Curie
+) all tests pass

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