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

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

Finishing the implementation of expand domain transformation

+) Make use of updated new functions of class Mesh to compute neighboring cells
+) Make change to some minor stuffs

Test
+) On Curie
+) The transformation works correctly

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