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

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

Masked interpolated points will have default_value (if defined)

+) Default value is used to initialize spatial transform filter
+) Add a new case for testing masked points

Test
+) On Curie
+) All interpolation tests work

  • Property svn:executable set to *
File size: 31.2 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      std::list<std::vector<bool> >().swap(localMaskOnGridDest_);
441    }
442    else
443      return;
444  }
445
446  CContext* context = CContext::getCurrent();
447  CContextClient* client = context->client;
448
449  ListAlgoType::const_iterator itb = listAlgos_.begin(),
450                               ite = listAlgos_.end(), it;
451
452  CGenericAlgorithmTransformation* algo = 0;
453  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
454  for (it = itb; it != ite; ++it)
455  {
456    int elementPositionInGrid = it->first;
457    ETranformationType transType = (it->second).first;
458    int transformationOrder = (it->second).second;
459    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
460
461    // Create a temporary grid destination which contains transformed element of grid destination and
462    // non-transformed elements fo grid source
463    setUpGridDestination(elementPositionInGrid, transType, nbAgloTransformation);
464
465    // First of all, select an algorithm
466    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
467    {
468      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
469      algo = algoTransformation_.back();
470    }
471    else
472      algo = algoTransformation_[std::distance(itb, it)];
473
474    if (0 != algo) // Only registered transformation can be executed
475    {
476      algo->computeIndexSourceMapping(dataAuxInputs);
477
478      // ComputeTransformation of global index of each element
479      int elementPosition = it->first;
480      algo->computeGlobalSourceIndex(elementPosition,
481                                     gridSource_,
482                                     tmpGridDestination_,
483                                     globaIndexWeightFromSrcToDst);
484
485      // Compute transformation of global indexes among grids
486      computeTransformationMapping(globaIndexWeightFromSrcToDst);
487
488      if (1 < nbAlgos_)
489      {
490        // Now grid destination becomes grid source in a new transformation
491        if (nbAgloTransformation != (nbAlgos_-1)) setUpGridSource(elementPositionInGrid, transType, nbAgloTransformation);
492      }
493      ++nbAgloTransformation;
494    }
495  }
496}
497
498/*!
499  Compute exchange index between grid source and grid destination
500  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
501*/
502void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
503{
504  CContext* context = CContext::getCurrent();
505  CContextClient* client = context->client;
506  int nbClient = client->clientSize;
507  int clientRank = client->clientRank;
508
509  // Recalculate the distribution of grid destination
510  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
511  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
512
513  // Update number of local index on each transformation
514  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
515  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
516  localMaskOnGridDest_.push_back(std::vector<bool>());
517  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
518  tmpMask.resize(nbLocalIndex,false);
519
520  // Find out number of index sent from grid source and number of index received on grid destination
521  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
522                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
523  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
524  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
525  int connectedClient = globaIndexWeightFromSrcToDst.size();
526  int* recvCount=new int[nbClient];
527  int* displ=new int[nbClient];
528  int* sendRankBuff=new int[connectedClient];
529  int* sendSizeBuff=new int[connectedClient];
530  int n = 0;
531  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
532  {
533    sendRankBuff[n] = itIndex->first;
534    const SendIndexMap& sendIndexMap = itIndex->second;
535    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
536    int sendSize = 0;
537    for (itSend = itbSend; itSend != iteSend; ++itSend)
538    {
539      sendSize += itSend->second.size();
540    }
541    sendSizeBuff[n] = sendSize;
542    sendRankSizeMap[itIndex->first] = sendSize;
543  }
544  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
545
546  displ[0]=0 ;
547  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
548  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
549  int* recvRankBuff=new int[recvSize];
550  int* recvSizeBuff=new int[recvSize];
551  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
552  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
553  for (int i = 0; i < nbClient; ++i)
554  {
555    int currentPos = displ[i];
556    for (int j = 0; j < recvCount[i]; ++j)
557      if (recvRankBuff[currentPos+j] == clientRank)
558      {
559        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
560      }
561  }
562
563  // Sending global index of grid source to corresponding process as well as the corresponding mask
564  std::vector<MPI_Request> requests;
565  std::vector<MPI_Status> status;
566  boost::unordered_map<int, unsigned char* > recvMaskDst;
567  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
568  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
569  {
570    int recvRank = itRecv->first;
571    int recvSize = itRecv->second;
572    recvMaskDst[recvRank] = new unsigned char [recvSize];
573    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
574
575    requests.push_back(MPI_Request());
576    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
577    requests.push_back(MPI_Request());
578    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
579  }
580
581  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
582  boost::unordered_map<int, CArray<double,1> > weightDst;
583  boost::unordered_map<int, unsigned char* > sendMaskDst;
584  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
585  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
586  {
587    int sendRank = itIndex->first;
588    int sendSize = sendRankSizeMap[sendRank];
589    const SendIndexMap& sendIndexMap = itIndex->second;
590    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
591    globalIndexDst[sendRank].resize(sendSize);
592    weightDst[sendRank].resize(sendSize);
593    sendMaskDst[sendRank] = new unsigned char [sendSize];
594    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
595    int countIndex = 0;
596    for (itSend = itbSend; itSend != iteSend; ++itSend)
597    {
598      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
599      for (int idx = 0; idx < dstWeight.size(); ++idx)
600      {
601        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
602        weightDst[sendRank](countIndex) = dstWeight[idx].second;
603        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
604          sendMaskDst[sendRank][countIndex] = 1;
605        else
606          sendMaskDst[sendRank][countIndex] = 0;
607        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
608        ++countIndex;
609      }
610    }
611
612    // Send global index source and mask
613    requests.push_back(MPI_Request());
614    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
615    requests.push_back(MPI_Request());
616    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
617  }
618
619  status.resize(requests.size());
620  MPI_Waitall(requests.size(), &requests[0], &status[0]);
621
622  // 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
623  std::vector<MPI_Request>().swap(requests);
624  std::vector<MPI_Status>().swap(status);
625  // Okie, on destination side, we will wait for information of masked index of source
626  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
627  {
628    int recvRank = itSend->first;
629    int recvSize = itSend->second;
630
631    requests.push_back(MPI_Request());
632    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
633  }
634
635  // Ok, now we fill in local index of grid source (we even count for masked index)
636  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
637  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
638  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
639  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
640  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
641  {
642    int recvRank = itRecv->first;
643    int recvSize = itRecv->second;
644    unsigned char* recvMask = recvMaskDst[recvRank];
645    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
646    int realSendSize = 0;
647    for (int idx = 0; idx < recvSize; ++idx)
648    {
649      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
650        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
651         ++realSendSize;
652        else // inform the destination that this index is masked
653         *(recvMask+idx) = 0;
654    }
655
656    tmpSend[recvRank].resize(realSendSize);
657    realSendSize = 0;
658    for (int idx = 0; idx < recvSize; ++idx)
659    {
660      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
661      {
662        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
663         ++realSendSize;
664      }
665    }
666
667    // Okie, now inform the destination which source index are masked
668    requests.push_back(MPI_Request());
669    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
670  }
671  status.resize(requests.size());
672  MPI_Waitall(requests.size(), &requests[0], &status[0]);
673
674  // Cool, now we can fill in local index of grid destination (counted for masked index)
675  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
676  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
677  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
678  {
679    int recvRank = itSend->first;
680    int recvSize = itSend->second;
681    unsigned char* recvMask = sendMaskDst[recvRank];
682
683    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
684    CArray<double,1>& recvWeightDst = weightDst[recvRank];
685    int realRecvSize = 0;
686    for (int idx = 0; idx < recvSize; ++idx)
687    {
688      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
689         ++realRecvSize;
690    }
691
692    int localIndexDst;
693    recvTmp[recvRank].resize(realRecvSize);
694    realRecvSize = 0;
695    for (int idx = 0; idx < recvSize; ++idx)
696    {
697      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
698      {
699        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
700        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
701        tmpMask[globalLocalIndexGridDestSendToServer[recvIndexDst(idx)]] = true;
702         ++realRecvSize;
703      }
704    }
705  }
706
707  delete [] recvCount;
708  delete [] displ;
709  delete [] sendRankBuff;
710  delete [] recvRankBuff;
711  delete [] sendSizeBuff;
712  delete [] recvSizeBuff;
713
714  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
715  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
716    delete [] itChar->second;
717  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
718    delete [] itChar->second;
719  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
720  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
721    delete [] itLong->second;
722  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
723    delete [] itLong->second;
724
725}
726
727bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
728{
729  bool res;
730  switch (transType)
731  {
732    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
733     res = true;
734     break;
735    default:
736     res = false;
737     break;
738  }
739
740  return res;
741}
742
743/*!
744  Local index of data which need sending from the grid source
745  \return local index of data
746*/
747const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
748{
749  return localIndexToSendFromGridSource_;
750}
751
752/*!
753  Local index of data which will be received on the grid destination
754  \return local index of data
755*/
756const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
757{
758  return localIndexToReceiveOnGridDest_;
759}
760
761/*!
762 Number of index will be received on the grid destination
763  \return number of index of data
764*/
765const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
766{
767  return nbLocalIndexOnGridDest_;
768}
769
770/*!
771  Local mask of data which will be received on the grid destination
772  \return local mask of data
773*/
774const std::list<std::vector<bool> >& CGridTransformation::getLocalMaskIndexOnGridDest() const
775{
776  return localMaskOnGridDest_;
777}
778
779}
Note: See TracBrowser for help on using the repository browser.