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

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

Adding a new type of element into grid: Scalar

+) Add a new node Scalar for xml
+) Make some change on writing scalar value
+) Reorganize some codes
+) Remove some redundant codes

Test
+) On Curie
+) All tests pass

  • Property svn:executable set to *
File size: 25.4 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 : CGridTransformationSelector(destination, source),
26  tmpGridDestination_(destination), originalGridSource_(source),
27  tempGridSrcs_(), tempGridDests_(),
28  dynamicalTransformation_(false), timeStamp_()
29{
30}
31
32CGridTransformation::~CGridTransformation()
33{
34}
35
36/*!
37  Select algorithm of a scalar correspoding to its transformation type and its position in each element
38  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
39  \param [in] transType transformation type, for now we have
40  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
41*/
42void CGridTransformation::selectScalarAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
43{
44
45}
46
47/*!
48  Select algorithm of an axis correspoding to its transformation type and its position in each element
49  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
50  \param [in] transType transformation type, for now we have zoom_axis, inverse_axis, interpolate_axis
51  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
52*/
53void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
54{
55  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
56  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
57
58  int axisIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
59  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
60  CAxis::TransMapTypes::const_iterator it = trans.begin();
61
62  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
63
64  CZoomAxis* zoomAxis = 0;
65  CInterpolateAxis* interpAxis = 0;
66  CGenericAlgorithmTransformation* algo = 0;
67  switch (transType)
68  {
69    case TRANS_INTERPOLATE_AXIS:
70      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
71      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
72      break;
73    case TRANS_ZOOM_AXIS:
74      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
75      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
76      break;
77    case TRANS_INVERSE_AXIS:
78      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
79      break;
80    default:
81      break;
82  }
83  algoTransformation_.push_back(algo);
84}
85
86/*!
87  Select algorithm of a domain correspoding to its transformation type and its position in each element
88  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
89  \param [in] transType transformation type, for now we have zoom_domain, interpolate_domain
90  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
91*/
92void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
93{
94  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
95  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
96
97  int domainIndex =  elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
98  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
99  CDomain::TransMapTypes::const_iterator it = trans.begin();
100
101  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
102
103  CZoomDomain* zoomDomain = 0;
104  CInterpolateDomain* interpFileDomain = 0;
105  CGenericAlgorithmTransformation* algo = 0;
106  switch (transType)
107  {
108    case TRANS_INTERPOLATE_DOMAIN:
109      interpFileDomain = dynamic_cast<CInterpolateDomain*> (it->second);
110      algo = new CDomainAlgorithmInterpolate(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
111      break;
112    case TRANS_ZOOM_DOMAIN:
113      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
114      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
115      break;
116    default:
117      break;
118  }
119  algoTransformation_.push_back(algo);
120}
121
122/*!
123  If there are more than one transformation, a new temporary grid will be created and it will play the role of grid destination.
124This new created one keeps a pointer to the real transformed element of grid destination and generate new copies of other elements from grid source.
125  \param [in] elementPositionInGrid position of element in grid
126  \param [in] transType transformation type
127*/
128void CGridTransformation::setUpGridDestination(int elementPositionInGrid, ETranformationType transType, int nbTransformation)
129{
130  if (isSpecialTransformation(transType)) return;
131
132  if (!tempGridDests_.empty() && (getNbAlgo() == tempGridDests_.size()))
133  {
134    tempGridDests_.resize(0);
135  }
136
137  std::vector<CScalar*> scalarListDestP = gridDestination_->getScalars();
138  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarDst;
139
140  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
141  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisDst;
142
143  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
144  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainDst;
145
146  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
147  CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
148
149  int scalarIndex = -1, axisIndex = -1, domainIndex = -1;
150  switch (transType)
151  {
152    case TRANS_INTERPOLATE_DOMAIN:
153    case TRANS_ZOOM_DOMAIN:
154      domainIndex = elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
155      break;
156
157    case TRANS_INTERPOLATE_AXIS:
158    case TRANS_ZOOM_AXIS:
159    case TRANS_INVERSE_AXIS:
160      axisIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
161      break;
162    default:
163      break;
164  }
165
166  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
167  {
168    int dimElementDst = axisDomainOrderDst(idx);
169    if (2 == dimElementDst)
170    {
171      if (elementPositionInGrid == idx)
172        domainDst.push_back(domListDestP[domainIndex]);
173      else
174        domainDst.push_back(domListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]]);
175    }
176    else if (1 == dimElementDst)
177    {
178      if (elementPositionInGrid == idx)
179        axisDst.push_back(axisListDestP[axisIndex]);
180      else
181        axisDst.push_back(axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]]);
182    }
183    else
184    {
185      if (elementPositionInGrid == idx)
186        scalarDst.push_back(scalarListDestP[scalarIndex]);
187      else
188        scalarDst.push_back(scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[elementPositionInGrid]]);
189    }
190  }
191
192  tmpGridDestination_ = CGrid::createGrid(domainDst, axisDst, scalarDst, gridDestination_->axis_domain_order);
193  tmpGridDestination_->computeGridGlobalDimension(domainDst, axisDst, scalarDst, gridDestination_->axis_domain_order);
194  tempGridDests_.push_back(tmpGridDestination_);
195}
196
197/*!
198  Assign the current grid destination to the grid source in the new transformation.
199The current grid destination plays the role of grid source in next transformation (if any).
200Only element on which the transformation is performed is modified
201  \param [in] elementPositionInGrid position of element in grid
202  \param [in] transType transformation type
203*/
204void CGridTransformation::setUpGridSource(int elementPositionInGrid, ETranformationType transType, int nbTransformation)
205{
206  if (!tempGridSrcs_.empty() && (getNbAlgo()-1) == tempGridSrcs_.size())
207  {
208    tempGridSrcs_.resize(0);
209  }
210
211  std::vector<CScalar*> scalarListDestP = tmpGridDestination_->getScalars();
212  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarSrc;
213
214  std::vector<CAxis*> axisListDestP = tmpGridDestination_->getAxis();
215  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
216
217  std::vector<CDomain*> domListDestP = tmpGridDestination_->getDomains();
218  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
219
220  CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
221
222  int axisIndex = -1, domainIndex = -1, scalarIndex = -1;
223  int axisListIndex = 0, domainListIndex = 0, scalarListIndex = 0;
224  switch (transType)
225  {
226    case TRANS_INTERPOLATE_DOMAIN:
227    case TRANS_ZOOM_DOMAIN:
228      domainIndex = elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
229      break;
230
231    case TRANS_INTERPOLATE_AXIS:
232    case TRANS_ZOOM_AXIS:
233    case TRANS_INVERSE_AXIS:
234      axisIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
235      break;
236    default:
237      break;
238  }
239
240  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
241  {
242    int dimElementDst = axisDomainOrderDst(idx);
243    if (2 == dimElementDst)
244    {
245      if (elementPositionInGrid == idx)
246        domainSrc.push_back(domListDestP[domainIndex]);
247      else
248      {
249        CDomain* domain = CDomain::createDomain();
250        domain->domain_ref.setValue(domListDestP[domainListIndex]->getId());
251        domain->solveRefInheritance(true);
252        domain->checkAttributesOnClient();
253        domainSrc.push_back(domain);
254      }
255      ++domainListIndex;
256    }
257    else if (1 == dimElementDst)
258    {
259      if (elementPositionInGrid == idx)
260        axisSrc.push_back(axisListDestP[axisIndex]);
261      else
262      {
263        CAxis* axis = CAxis::createAxis();
264        axis->axis_ref.setValue(axisListDestP[axisListIndex]->getId());
265        axis->solveRefInheritance(true);
266        axis->checkAttributesOnClient();
267        axisSrc.push_back(axis);
268      }
269      ++axisListIndex;
270    }
271    else
272    {
273      if (elementPositionInGrid == idx)
274        scalarSrc.push_back(scalarListDestP[scalarIndex]);
275      else
276      {
277        CScalar* scalar = CScalar::createScalar();
278        scalar->scalar_ref.setValue(scalarListDestP[scalarListIndex]->getId());
279        scalar->solveRefInheritance(true);
280        scalar->checkAttributesOnClient();
281        scalarSrc.push_back(scalar);
282      }
283      ++scalarListIndex;
284    }
285  }
286
287  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, scalarSrc, tmpGridDestination_->axis_domain_order);
288  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, scalarSrc, tmpGridDestination_->axis_domain_order);
289
290  tempGridSrcs_.push_back(gridSource_);
291}
292
293/*!
294  Perform all transformations
295  For each transformation, there are some things to do:
296  -) Chose the correct algorithm by transformation type and position of element
297  -) Calculate the mapping of global index between the current grid source and grid destination
298  -) Calculate the mapping of global index between current grid DESTINATION and grid SOURCE
299  -) Make current grid destination become grid source in the next transformation
300*/
301void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
302{
303  if (nbNormalAlgos_ < 1) return;
304  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
305  if (dynamicalTransformation_)
306  {
307    if (timeStamp_.insert(timeStamp).second)   //Reset map
308    {
309      std::list<SendingIndexGridSourceMap>().swap(localIndexToSendFromGridSource_);
310      std::list<RecvIndexGridDestinationMap>().swap(localIndexToReceiveOnGridDest_);
311      std::list<size_t>().swap(nbLocalIndexOnGridDest_);
312      std::list<std::vector<bool> >().swap(localMaskOnGridDest_);
313    }
314    else
315      return;
316  }
317
318  CContext* context = CContext::getCurrent();
319  CContextClient* client = context->client;
320
321  ListAlgoType::const_iterator itb = listAlgos_.begin(),
322                               ite = listAlgos_.end(), it;
323
324  CGenericAlgorithmTransformation* algo = 0;
325  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
326  for (it = itb; it != ite; ++it)
327  {
328    int elementPositionInGrid = it->first;
329    ETranformationType transType = (it->second).first;
330    int transformationOrder = (it->second).second;
331    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
332
333    // Create a temporary grid destination which contains transformed element of grid destination and
334    // non-transformed elements fo grid source
335    setUpGridDestination(elementPositionInGrid, transType, nbAgloTransformation);
336
337    // First of all, select an algorithm
338    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
339    {
340      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
341      algo = algoTransformation_.back();
342    }
343    else
344      algo = algoTransformation_[std::distance(itb, it)];
345
346    if (0 != algo) // Only registered transformation can be executed
347    {
348      algo->computeIndexSourceMapping(dataAuxInputs);
349
350      // ComputeTransformation of global index of each element
351      int elementPosition = it->first;
352      algo->computeGlobalSourceIndex(elementPosition,
353                                     gridSource_,
354                                     tmpGridDestination_,
355                                     globaIndexWeightFromSrcToDst);
356
357      // Compute transformation of global indexes among grids
358      computeTransformationMapping(globaIndexWeightFromSrcToDst);
359
360      if (1 < nbNormalAlgos_)
361      {
362        // Now grid destination becomes grid source in a new transformation
363        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid, transType, nbAgloTransformation);
364      }
365      ++nbAgloTransformation;
366    }
367  }
368}
369
370/*!
371  Compute exchange index between grid source and grid destination
372  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
373*/
374void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
375{
376  CContext* context = CContext::getCurrent();
377  CContextClient* client = context->client;
378  int nbClient = client->clientSize;
379  int clientRank = client->clientRank;
380
381  // Recalculate the distribution of grid destination
382  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
383  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
384
385  // Update number of local index on each transformation
386  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
387  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
388  localMaskOnGridDest_.push_back(std::vector<bool>());
389  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
390  tmpMask.resize(nbLocalIndex,false);
391
392  // Find out number of index sent from grid source and number of index received on grid destination
393  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
394                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
395  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
396  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
397  int connectedClient = globaIndexWeightFromSrcToDst.size();
398  int* recvCount=new int[nbClient];
399  int* displ=new int[nbClient];
400  int* sendRankBuff=new int[connectedClient];
401  int* sendSizeBuff=new int[connectedClient];
402  int n = 0;
403  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
404  {
405    sendRankBuff[n] = itIndex->first;
406    const SendIndexMap& sendIndexMap = itIndex->second;
407    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
408    int sendSize = 0;
409    for (itSend = itbSend; itSend != iteSend; ++itSend)
410    {
411      sendSize += itSend->second.size();
412    }
413    sendSizeBuff[n] = sendSize;
414    sendRankSizeMap[itIndex->first] = sendSize;
415  }
416  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
417
418  displ[0]=0 ;
419  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
420  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
421  int* recvRankBuff=new int[recvSize];
422  int* recvSizeBuff=new int[recvSize];
423  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
424  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
425  for (int i = 0; i < nbClient; ++i)
426  {
427    int currentPos = displ[i];
428    for (int j = 0; j < recvCount[i]; ++j)
429      if (recvRankBuff[currentPos+j] == clientRank)
430      {
431        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
432      }
433  }
434
435  // Sending global index of grid source to corresponding process as well as the corresponding mask
436  std::vector<MPI_Request> requests;
437  std::vector<MPI_Status> status;
438  boost::unordered_map<int, unsigned char* > recvMaskDst;
439  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
440  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
441  {
442    int recvRank = itRecv->first;
443    int recvSize = itRecv->second;
444    recvMaskDst[recvRank] = new unsigned char [recvSize];
445    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
446
447    requests.push_back(MPI_Request());
448    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
449    requests.push_back(MPI_Request());
450    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
451  }
452
453  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
454  boost::unordered_map<int, CArray<double,1> > weightDst;
455  boost::unordered_map<int, unsigned char* > sendMaskDst;
456  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
457  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
458  {
459    int sendRank = itIndex->first;
460    int sendSize = sendRankSizeMap[sendRank];
461    const SendIndexMap& sendIndexMap = itIndex->second;
462    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
463    globalIndexDst[sendRank].resize(sendSize);
464    weightDst[sendRank].resize(sendSize);
465    sendMaskDst[sendRank] = new unsigned char [sendSize];
466    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
467    int countIndex = 0;
468    for (itSend = itbSend; itSend != iteSend; ++itSend)
469    {
470      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
471      for (int idx = 0; idx < dstWeight.size(); ++idx)
472      {
473        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
474        weightDst[sendRank](countIndex) = dstWeight[idx].second;
475        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
476          sendMaskDst[sendRank][countIndex] = 1;
477        else
478          sendMaskDst[sendRank][countIndex] = 0;
479        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
480        ++countIndex;
481      }
482    }
483
484    // Send global index source and mask
485    requests.push_back(MPI_Request());
486    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
487    requests.push_back(MPI_Request());
488    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
489  }
490
491  status.resize(requests.size());
492  MPI_Waitall(requests.size(), &requests[0], &status[0]);
493
494  // 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
495  std::vector<MPI_Request>().swap(requests);
496  std::vector<MPI_Status>().swap(status);
497  // Okie, on destination side, we will wait for information of masked index of source
498  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
499  {
500    int recvRank = itSend->first;
501    int recvSize = itSend->second;
502
503    requests.push_back(MPI_Request());
504    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
505  }
506
507  // Ok, now we fill in local index of grid source (we even count for masked index)
508  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
509  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
510  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
511  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
512  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
513  {
514    int recvRank = itRecv->first;
515    int recvSize = itRecv->second;
516    unsigned char* recvMask = recvMaskDst[recvRank];
517    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
518    int realSendSize = 0;
519    for (int idx = 0; idx < recvSize; ++idx)
520    {
521      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
522        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
523         ++realSendSize;
524        else // inform the destination that this index is masked
525         *(recvMask+idx) = 0;
526    }
527
528    tmpSend[recvRank].resize(realSendSize);
529    realSendSize = 0;
530    for (int idx = 0; idx < recvSize; ++idx)
531    {
532      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
533      {
534        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
535         ++realSendSize;
536      }
537    }
538
539    // Okie, now inform the destination which source index are masked
540    requests.push_back(MPI_Request());
541    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
542  }
543  status.resize(requests.size());
544  MPI_Waitall(requests.size(), &requests[0], &status[0]);
545
546  // Cool, now we can fill in local index of grid destination (counted for masked index)
547  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
548  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
549  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
550  {
551    int recvRank = itSend->first;
552    int recvSize = itSend->second;
553    unsigned char* recvMask = sendMaskDst[recvRank];
554
555    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
556    CArray<double,1>& recvWeightDst = weightDst[recvRank];
557    int realRecvSize = 0;
558    for (int idx = 0; idx < recvSize; ++idx)
559    {
560      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
561         ++realRecvSize;
562    }
563
564    int localIndexDst;
565    recvTmp[recvRank].resize(realRecvSize);
566    realRecvSize = 0;
567    for (int idx = 0; idx < recvSize; ++idx)
568    {
569      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
570      {
571        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
572        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
573        tmpMask[globalLocalIndexGridDestSendToServer[recvIndexDst(idx)]] = true;
574         ++realRecvSize;
575      }
576    }
577  }
578
579  delete [] recvCount;
580  delete [] displ;
581  delete [] sendRankBuff;
582  delete [] recvRankBuff;
583  delete [] sendSizeBuff;
584  delete [] recvSizeBuff;
585
586  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
587  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
588    delete [] itChar->second;
589  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
590    delete [] itChar->second;
591  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
592  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
593    delete [] itLong->second;
594  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
595    delete [] itLong->second;
596
597}
598
599/*!
600  Local index of data which need sending from the grid source
601  \return local index of data
602*/
603const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
604{
605  return localIndexToSendFromGridSource_;
606}
607
608/*!
609  Local index of data which will be received on the grid destination
610  \return local index of data
611*/
612const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
613{
614  return localIndexToReceiveOnGridDest_;
615}
616
617/*!
618 Number of index will be received on the grid destination
619  \return number of index of data
620*/
621const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
622{
623  return nbLocalIndexOnGridDest_;
624}
625
626/*!
627  Local mask of data which will be received on the grid destination
628  \return local mask of data
629*/
630const std::list<std::vector<bool> >& CGridTransformation::getLocalMaskIndexOnGridDest() const
631{
632  return localMaskOnGridDest_;
633}
634
635}
Note: See TracBrowser for help on using the repository browser.