source: XIOS/trunk/src/transformation/domain_algorithm_expand.cpp @ 1982

Last change on this file since 1982 was 1972, checked in by yushan, 3 years ago

trunk : debug domain_expand
File size: 27.2 KB
Line 
1/*!
2   \file domain_algorithm_expand.cpp
3   \author Ha NGUYEN
4   \since 08 Aug 2016
5   \date 19 Sep 2016
6
7   \brief Algorithm for expanding an domain.
8 */
9#include "domain_algorithm_expand.hpp"
10#include "expand_domain.hpp"
11#include "mesh.hpp"
12#include "domain.hpp"
13#include "grid.hpp"
14#include "grid_transformation_factory_impl.hpp"
15#include "context.hpp"
16#include "context_client.hpp"
17
18namespace xios {
19CGenericAlgorithmTransformation* CDomainAlgorithmExpand::create(CGrid* gridDst, CGrid* gridSrc,
20                                                               CTransformation<CDomain>* transformation,
21                                                               int elementPositionInGrid,
22                                                               std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
23                                                               std::map<int, int>& elementPositionInGridSrc2AxisPosition,
24                                                               std::map<int, int>& elementPositionInGridSrc2DomainPosition,
25                                                               std::map<int, int>& elementPositionInGridDst2ScalarPosition,
26                                                               std::map<int, int>& elementPositionInGridDst2AxisPosition,
27                                                               std::map<int, int>& elementPositionInGridDst2DomainPosition)
28TRY
29{
30  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
31  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
32
33  CExpandDomain* expandDomain = dynamic_cast<CExpandDomain*> (transformation);
34  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
35  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
36
37  return (new CDomainAlgorithmExpand(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], expandDomain));
38}
39CATCH
40
41bool CDomainAlgorithmExpand::registerTrans()
42TRY
43{
44  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXPAND_DOMAIN, create);
45}
46CATCH
47
48CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,
49                                               CDomain* domainSource,
50                                               CExpandDomain* expandDomain)
51: CDomainAlgorithmTransformation(domainDestination, domainSource),
52  isXPeriodic_(false), isYPeriodic_(false)
53TRY
54{
55  if (domainDestination == domainSource)
56  {
57    ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",
58           << "Domain source and domain destination are the same. Please make sure domain destination refers to domain source" << std::endl
59           << "Domain source " <<domainSource->getId() << std::endl
60           << "Domain destination " <<domainDestination->getId() << std::endl);
61  }
62
63  this->type_ = (ELEMENT_MODIFICATION_WITH_DATA);
64  // Make sure domain source have all valid attributes
65  // domainSource->checkAllAttributes();
66  domainSource->initLonLatValue();
67  expandDomain->checkValid(domainDestination);
68  if (!expandDomain->i_periodic.isEmpty()) isXPeriodic_ = expandDomain->i_periodic;
69  if (!expandDomain->j_periodic.isEmpty()) isYPeriodic_ = expandDomain->j_periodic;
70 
71  switch (expandDomain->type)
72  {
73    case CExpandDomain::type_attr::node :
74      expandDomainNodeConnectivity(domainDestination,
75                                   domainSource);
76      break;
77    case CExpandDomain::type_attr::edge :
78      expandDomainEdgeConnectivity(domainDestination,
79                                   domainSource);
80      break;
81    default:
82      break;
83  }
84}
85CATCH
86
87/*!
88 *  Expand domain with edge-type neighbor
89 *  \param[in/out] domainDestination domain destination and will be modified
90 *  \param[in] domainSource domain source
91 */
92void CDomainAlgorithmExpand::expandDomainEdgeConnectivity(CDomain* domainDestination,
93                                                          CDomain* domainSource)
94TRY
95{
96  CContext* context = CContext::getCurrent();
97  CContextClient* client=context->client;
98
99  int type = 1; // For edge
100  CMesh mesh;
101  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
102  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
103  CArray<int,2> neighborsSrc;
104  switch (domainSource->type) {
105   case CDomain::type_attr::unstructured:     
106      mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
107      updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc);
108      break;
109   default:
110      updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc);
111      break;     
112  } 
113}
114CATCH
115
116/*!
117 *  Expand domain with node-type neighbor
118 *  \param[in/out] domainDestination domain destination and will be modified
119 *  \param[in] domainSource domain source
120 */
121void CDomainAlgorithmExpand::expandDomainNodeConnectivity(CDomain* domainDestination,
122                                                          CDomain* domainSource)
123TRY
124{
125  CContext* context = CContext::getCurrent();
126  CContextClient* client=context->client;
127
128  int type = 1; // For edge
129  CMesh mesh;
130  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
131  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
132  CArray<int,2> neighborsSrc;
133  switch (domainSource->type) {
134   case CDomain::type_attr::unstructured:     
135      mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
136      updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc);
137      break;
138   default:
139      updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc);
140      break;     
141  } 
142}
143CATCH
144
145/*!
146 *  Extend rectilinear or curvilinear domain destination and update its attributes
147 *  Suppose that domain destination and domain source have the same values for all attributes (by inheritance)
148 *  \param [in/out] domainDestination domain destination
149 *  \param [in] domainSource domain source
150 *  \param [in] neighborsDomainSrc neighbor of domain source. For now, we don't need it for rectilinear
151 */
152void CDomainAlgorithmExpand::updateRectilinearDomainAttributes(CDomain* domainDestination,
153                                                               CDomain* domainSource,
154                                                               CArray<int,2>& neighborsDomainSrc)
155TRY
156{
157  int index, globalIndex, idx;
158  int iindexDst, jindexDst, globIndexDst;
159  int iindexSrc, jindexSrc, globIndexSrc;
160  CContext* context = CContext::getCurrent();
161  CContextClient* client=context->client;
162
163  // First of all, "copy" all attributes of domain source to domain destination
164  StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue()
165                                                                      : "";
166  if (domainDstRef != domainSource->getId())
167  {
168    domainDestination->domain_ref.setValue(domainSource->getId());
169    domainDestination->solveRefInheritance(true);
170  }
171
172  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
173  else domainDestination->domain_ref.setValue(domainDstRef);
174
175 
176  // Here are attributes of source need tranfering
177  int niGloSrc = domainSource->ni_glo;
178  int njGloSrc = domainSource->nj_glo;
179  int niSrc = domainSource->ni, ibegin = domainSource->ibegin;
180  int njSrc = domainSource->nj, jbegin = domainSource->jbegin;
181  int dataDimSrc = domainSource->data_dim;
182  CArray<bool,1>& mask_1d_src = domainSource->domainMask;
183  CArray<int,1>& i_index_src = domainSource->i_index;
184  CArray<int,1>& j_index_src = domainSource->j_index;
185  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
186  CArray<int,1>& data_j_index_src = domainSource->data_j_index;
187  int data_i_begin_src = domainSource->data_ibegin;
188  int data_j_begin_src = domainSource->data_jbegin;
189  CArray<double,1>& lon_src = domainSource->lonvalue;
190  CArray<double,1>& lat_src = domainSource->latvalue;
191  CArray<double,1> lon_g;
192  CArray<double,1> lat_g;
193
194  // We need to generate boundary for longitude and latitude
195  if (domainSource->bounds_lon_1d.isEmpty() || domainSource->bounds_lat_1d.isEmpty())
196  {
197    CArray<double,1> lon = lon_src(Range(0,niSrc-1));
198    CArray<double,1> lat = lat_src(Range(0,lat_src.numElements()-niSrc,niSrc));
199    CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
200    CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
201    domainSource->AllgatherRectilinearLonLat(lon_src, lat_src, lon_g,lat_g);
202    domainSource->fillInRectilinearBoundLonLat(lon_g, lat_g, bounds_lon_src, bounds_lat_src);
203  }
204
205
206  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
207  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
208
209  int nVertex       = bounds_lon_src.shape()[0];
210  int oldNbLocal = i_index_src.numElements();
211  // Calculate ni, nj by using i_index and j_index
212  int niSrcByIndex = max(i_index_src) - min(i_index_src) + 1;
213  int njSrcByIndex = max(j_index_src) - min(j_index_src) + 1; 
214  int dataIindexBoundSrc = (1 == dataDimSrc) ? (niSrcByIndex * njSrcByIndex) : niSrcByIndex;
215  int dataJindexBoundSrc = (1 == dataDimSrc) ? (niSrcByIndex * njSrcByIndex) : njSrcByIndex;
216
217  // Uncompress data_i_index, data_j_index
218  CArray<int,1> data_i_index_src_full(oldNbLocal);
219  CArray<int,1> data_j_index_src_full(oldNbLocal);
220  int nbUnMaskedPointOnLocalDomain = 0;
221  data_i_index_src_full = -1; // Suppose all values are masked
222  data_j_index_src_full = -1; // Suppose all values are masked
223  for (idx = 0; idx < data_i_index_src.numElements(); ++idx)
224  {
225    int dataIidx = data_i_index_src(idx) + data_i_begin_src;
226    int dataJidx = data_j_index_src(idx) + data_j_begin_src;
227    if ((0 <= dataIidx) && (dataIidx < dataIindexBoundSrc) &&
228        (0 <= dataJidx) && (dataJidx < dataJindexBoundSrc))
229    {
230      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIidx;
231      data_j_index_src_full(nbUnMaskedPointOnLocalDomain) = dataJidx;
232      ++nbUnMaskedPointOnLocalDomain;
233    }
234  }
235
236  // Expand domain destination, not only local but also global
237  int niGloDst = niGloSrc + 2;
238  int njGloDst = njGloSrc + 2;
239  int niDst = niSrc + 2;
240  int njDst = njSrc + 2;
241  domainDestination->ni_glo.setValue(niGloDst);
242  domainDestination->nj_glo.setValue(njGloDst);
243  domainDestination->ni.setValue(niDst);
244  domainDestination->nj.setValue(njDst);
245
246  CArray<bool,1>& mask_1d_dst = domainDestination->domainMask;
247  CArray<int,1>& i_index_dst  = domainDestination->i_index;
248  CArray<int,1>& j_index_dst  = domainDestination->j_index; 
249  CArray<int,1>& data_i_index_dst  = domainDestination->data_i_index;
250  CArray<int,1>& data_j_index_dst  = domainDestination->data_j_index;
251 
252  // Make sure that we use only lonvalue_client, latvalue_client
253  if (!domainDestination->lonvalue_1d.isEmpty()) domainDestination->lonvalue_1d.reset();
254  if (!domainDestination->latvalue_1d.isEmpty()) domainDestination->latvalue_1d.reset();
255  if (!domainDestination->lonvalue_2d.isEmpty()) domainDestination->lonvalue_2d.reset();
256  if (!domainDestination->latvalue_2d.isEmpty()) domainDestination->latvalue_2d.reset();
257
258  // Recalculate i_index, j_index of extended domain
259  // Should be enough for common case, but if we have arbitrary distribution?
260  int newNbLocalDst = niDst * njDst;     
261
262  mask_1d_dst.resize(newNbLocalDst);
263  i_index_dst.resize(newNbLocalDst);
264  j_index_dst.resize(newNbLocalDst);
265  CArray<int,1> data_i_index_dst_full(newNbLocalDst);
266  CArray<int,1> data_j_index_dst_full(newNbLocalDst);
267
268  domainDestination->lonvalue.resizeAndPreserve(newNbLocalDst);
269  domainDestination->latvalue.resizeAndPreserve(newNbLocalDst);
270  domainDestination->bounds_lon_1d.resizeAndPreserve(nVertex, newNbLocalDst);
271  domainDestination->bounds_lat_1d.resizeAndPreserve(nVertex, newNbLocalDst);
272  CArray<double,1>& lon_dst   = domainDestination->lonvalue;
273  CArray<double,1>& lat_dst   = domainDestination->latvalue;
274  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
275  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
276
277  // Update i_index, j_index
278  for (int j = 0; j < njDst; ++j)
279    for (int i = 0; i < niDst; ++i)
280    {
281      idx = j * niDst + i; 
282      i_index_dst(idx) = i + ibegin;
283      j_index_dst(idx) = j + jbegin;
284    }
285 
286
287  // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...)
288  // Global index mapping between destination and source
289  this->transformationMapping_.resize(1);
290  this->transformationWeight_.resize(1);
291  TransformationIndexMap& transMap = this->transformationMapping_[0];
292  TransformationWeightMap& transWeight = this->transformationWeight_[0];
293
294  transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor()));
295  transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor()));
296 
297  // Index mapping for local domain
298  // Mapping global index of expanded domain into original one
299  // (Representing global index of expanded domain in form of global index of original one)
300  CArray<size_t,1> globalIndexSrcOnDstDomain(newNbLocalDst); 
301  for (idx = 0; idx < newNbLocalDst; ++idx)
302  {
303    iindexDst = i_index_dst(idx);
304    jindexDst = j_index_dst(idx);
305    globIndexDst = jindexDst * niGloDst + iindexDst;
306    globIndexSrc = (((jindexDst-1)+njGloSrc) % njGloSrc) * niGloSrc + (((iindexDst-1)+niGloSrc) % niGloSrc) ;
307    globalIndexSrcOnDstDomain(idx) = globIndexSrc;
308
309    transMap[globIndexDst].push_back(globIndexSrc);
310    transWeight[globIndexDst].push_back(1.0); 
311  }
312
313  // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...)
314  CClientClientDHTDouble::Index2VectorInfoTypeMap localData;
315  localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor()));
316
317  // Information exchanged among domains (attention to their order), number in parentheses presents size of data
318  // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1)
319  int dataPackageSize =  1 + 1 + // lon + lat
320                         nVertex + nVertex + //bounds_lon + bounds_lat
321                         1 + // mask_1d_dst;
322                         1 + 1; // data_i_index + data_j_index
323  // Initialize database
324  for (int idx = 0; idx < oldNbLocal; ++idx)
325  {
326    index = i_index_src(idx) + j_index_src(idx) * niGloSrc;
327    localData[index].resize(dataPackageSize);
328    std::vector<double>& data = localData[index];
329
330    //Pack data
331    int dataIdx = 0;
332    data[dataIdx] = lon_src(idx);++dataIdx;
333    data[dataIdx] = lat_src(idx);++dataIdx;
334    for (int i = 0; i < nVertex; ++i)
335    {
336      data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx;
337    }
338    for (int i = 0; i < nVertex; ++i)
339    {
340      data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
341    }
342    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx;
343    data[dataIdx] = data_i_index_src_full(idx);++dataIdx;
344    data[dataIdx] = data_j_index_src_full(idx);
345  }
346
347  CClientClientDHTDouble dhtData(localData,client->intraComm);
348  dhtData.computeIndexInfoMapping(globalIndexSrcOnDstDomain);
349  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
350  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
351
352  // Ok get all data for destination
353  // If domain is not periodic, then we mask all extended part.
354  int nbUnMaskedPointOnExtendedPart = 0, remainder = 0, dataIIndex, dataJIndex;
355  size_t nIdx;
356  double maskValue = 1.0;
357  for (index = 0; index < newNbLocalDst; ++index)
358  {
359     nIdx = globalIndexSrcOnDstDomain(index);
360     it = neighborData.find(nIdx);
361     if (ite != it)
362     {
363        std::vector<double>& data = it->second;
364        // Unpack data
365        int dataIdx = 0;
366        lon_dst(index) = data[dataIdx]; ++dataIdx;
367        lat_dst(index) = data[dataIdx]; ++dataIdx;
368        for (int i = 0; i < nVertex; ++i)
369        {
370          bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx;
371        }
372        for (int i = 0; i < nVertex; ++i)
373        {
374          bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx;
375        }
376       
377        // Check whether we have x periodic. If we don't, we should mask all point at 0 and niGloDst-1
378        maskValue = data[dataIdx];
379        if (!isXPeriodic_) 
380        {
381          remainder = i_index_dst(index) % (niGloDst-1);
382          if (0 == remainder) 
383          {
384            maskValue = -1.0;
385          }
386        }
387
388        if (!isYPeriodic_) 
389        {
390          remainder = j_index_dst(index) % (njGloDst-1);
391          if (0 == remainder) 
392          {
393            maskValue = -1.0;
394          }
395        }
396
397        mask_1d_dst(index) = (1.0 == maskValue) ? true : false; ++dataIdx;
398
399        dataIIndex = (int) data[dataIdx];
400        if (!isXPeriodic_) 
401        {
402          remainder = i_index_dst(index) % (niGloDst-1);
403          if (0 == remainder) 
404          {
405            dataIIndex = -1;
406          }
407        }
408        data_i_index_dst_full(index) = dataIIndex; ++dataIdx;
409       
410        dataJIndex = (int) data[dataIdx];
411        if (!isYPeriodic_) 
412        {
413          remainder = j_index_dst(index) % (njGloDst-1);
414          if (0 == remainder) 
415          {
416            dataJIndex = -1;
417          }
418        }       
419        data_j_index_dst_full(index) = dataJIndex;
420
421        if ((0 <= data_i_index_dst_full(index)) && (0 <= data_j_index_dst_full(index)))
422        {
423          ++nbUnMaskedPointOnExtendedPart;
424        }
425     }
426  }
427
428 
429  // Finally, update data_i_index, data_j_index
430  int dataDstDim = domainDestination->data_dim;
431  data_i_index_dst.resize(nbUnMaskedPointOnExtendedPart);
432  data_j_index_dst.resize(nbUnMaskedPointOnExtendedPart); 
433  int count = 0; 
434  for (idx = 0; idx < newNbLocalDst; ++idx)
435  {
436    dataIIndex = data_i_index_dst_full(idx);
437    dataJIndex = data_j_index_dst_full(idx);
438    if ((0 <= dataIIndex) && (0 <= dataJIndex))
439    {
440      data_i_index_dst(count) = (1 == dataDstDim) ? idx : i_index_dst(idx) - i_index_dst(0);
441      data_j_index_dst(count) = (1 == dataDstDim) ? 0   : j_index_dst(idx) - j_index_dst(0);
442      ++count;
443    }
444  }
445
446  // Update data_ni, data_nj
447 
448  domainDestination->data_ni.setValue((1==dataDstDim) ? niDst * njDst : niDst);
449  domainDestination->data_nj.setValue((1==dataDstDim) ? niDst * njDst : njDst); 
450  domainDestination->data_ibegin.setValue(0);
451  domainDestination->data_jbegin.setValue(0);
452
453  // Update longitude and latitude
454  if (niSrc == domainSource->lonvalue_1d.numElements() && njSrc == domainSource->latvalue_1d.numElements()) // Ok, we have rectilinear here
455  {
456     domainDestination->lonvalue_1d.resize(niDst);
457     domainDestination->lonvalue_1d = lon_dst(Range(0,niDst-1));
458     domainDestination->latvalue_1d.resize(njDst);
459     domainDestination->latvalue_1d = lat_dst(Range(0,lat_dst.numElements()-niDst,niDst));
460  }
461  else // It should be curvilinear
462  {
463     domainDestination->lonvalue_1d.resize(lon_dst.numElements());
464     domainDestination->lonvalue_1d = lon_dst;
465     domainDestination->latvalue_1d.resize(lat_dst.numElements());
466     domainDestination->latvalue_1d = (lat_dst);
467  }
468   domainDestination->mask_1d.resize(domainDestination->domainMask.numElements()) ;
469   domainDestination->mask_1d=domainDestination->domainMask ;
470   domainDestination->computeLocalMask() ;
471}
472CATCH
473
474/*!
475 *  Extend domain destination and update its attributes
476 *  Suppose that domain destination and domain source have the same values for all attributes (by inheritance)
477 *  \param [in/out] domainDestination domain destination
478 *  \param [in] domainSource domain source
479 *  \param [in] neighborsDomainSrc domain extended part
480 */
481void CDomainAlgorithmExpand::updateUnstructuredDomainAttributes(CDomain* domainDestination,
482                                                                CDomain* domainSource,
483                                                                CArray<int,2>& neighborsDomainSrc)
484TRY
485{
486
487  CContext* context = CContext::getCurrent();
488  CContextClient* client=context->client;
489
490  // First of all, "copy" all attributes of domain source to domain destination
491  StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue()
492                                                                      : "";
493  if (domainDstRef != domainSource->getId())
494  {
495    domainDestination->domain_ref.setValue(domainSource->getId());
496    domainDestination->solveRefInheritance(true);
497  }
498
499  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
500  else domainDestination->domain_ref.setValue(domainDstRef);
501
502  // Now extend domain destination
503  int niGlob = domainSource->ni_glo;
504  CArray<bool,1>& mask_1d_src = domainSource->domainMask;
505  CArray<int,1>& i_index_src = domainSource->i_index;
506  CArray<double,1>& lon_src = domainSource->lonvalue_1d;
507  CArray<double,1>& lat_src = domainSource->latvalue_1d;
508  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
509  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
510  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
511
512  int oldNbLocal = i_index_src.numElements(), index, globalIndex;
513  // Uncompress data_i_index
514  CArray<int,1> data_i_index_src_full(oldNbLocal);
515  int nbUnMaskedPointOnLocalDomain = 0;
516  data_i_index_src_full = -1; // Suppose all values are masked
517  for (int idx = 0; idx < data_i_index_src.numElements(); ++idx)
518  {
519    int dataIdx = data_i_index_src(idx);
520    if ((0 <= dataIdx) && (dataIdx < oldNbLocal))
521    {
522      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIdx;
523      ++nbUnMaskedPointOnLocalDomain;
524    }
525  }
526
527  CArray<bool,1>& mask_1d_dst = domainDestination->domainMask;
528  CArray<int,1>& i_index_dst = domainDestination->i_index;
529  CArray<int,1>& j_index_dst = domainDestination->j_index;
530  CArray<double,1>& lon_dst = domainDestination->lonvalue_1d;
531  CArray<double,1>& lat_dst = domainDestination->latvalue_1d;
532  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
533  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
534  CArray<int,1>& data_i_index_dst = domainDestination->data_i_index;
535  CArray<int,1>& data_j_index_dst = domainDestination->data_j_index;
536
537  // Resize all array-like attributes of domain destination
538  int nbNeighbor    = neighborsDomainSrc.shape()[1];
539  int newNbLocalDst = nbNeighbor + oldNbLocal;
540  int nVertex       = bounds_lon_dst.shape()[0];
541
542  mask_1d_dst.resizeAndPreserve(newNbLocalDst);
543  i_index_dst.resizeAndPreserve(newNbLocalDst);
544  j_index_dst.resizeAndPreserve(newNbLocalDst);
545  lon_dst.resizeAndPreserve(newNbLocalDst);
546  lat_dst.resizeAndPreserve(newNbLocalDst);
547  bounds_lon_dst.resizeAndPreserve(nVertex, newNbLocalDst);
548  bounds_lat_dst.resizeAndPreserve(nVertex, newNbLocalDst);
549  CArray<int,1> data_i_index_dst_full(newNbLocalDst);
550  data_i_index_dst_full(Range(0,oldNbLocal-1)) = data_i_index_src_full;
551  data_i_index_dst_full(Range(oldNbLocal,newNbLocalDst-1)) = -1;
552
553  // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...)
554  // Global index mapping between destination and source
555  this->transformationMapping_.resize(1);
556  this->transformationWeight_.resize(1);
557  TransformationIndexMap& transMap = this->transformationMapping_[0];
558  TransformationWeightMap& transWeight = this->transformationWeight_[0];
559
560  transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor()));
561  transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor()));
562  // First, index mapping for local domain
563  for (int idx = 0; idx < oldNbLocal; ++idx)
564  {
565    index = i_index_dst(idx);
566    transMap[index].push_back(index);
567    transWeight[index].push_back(1.0);
568  }
569  // Then, index mapping for extended part
570  for (int idx = 0; idx < nbNeighbor; ++idx)
571  {
572    index = idx + oldNbLocal;
573    globalIndex = neighborsDomainSrc(0,idx);
574    i_index_dst(index) = globalIndex;
575    j_index_dst(index) = 0;
576    transMap[globalIndex].push_back(globalIndex);
577    transWeight[globalIndex].push_back(1.0);
578  }
579
580  // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...)
581  CClientClientDHTDouble::Index2VectorInfoTypeMap localData;
582  localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor()));
583  // Information exchanged among domains (attention to their order), number in parentheses presents size of data
584  // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1)
585  int dataPackageSize =  1 + 1 + // lon + lat
586                         nVertex + nVertex + //bounds_lon + bounds_lat
587                         1 + // mask_1d_dst;
588                         1; // data_i_index
589  // Initialize database
590  for (int idx = 0; idx < oldNbLocal; ++idx)
591  {
592    index = i_index_src(idx);
593    localData[index].resize(dataPackageSize);
594    std::vector<double>& data = localData[index];
595
596    //Pack data
597    int dataIdx = 0;
598    data[dataIdx] = lon_src(idx);++dataIdx;
599    data[dataIdx] = lat_src(idx);++dataIdx;
600    for (int i = 0; i < nVertex; ++i)
601    {
602      data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx;
603    }
604    for (int i = 0; i < nVertex; ++i)
605    {
606      data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
607    }
608    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx;
609    data[dataIdx] = data_i_index_src_full(idx);
610  }
611
612  CClientClientDHTDouble dhtData(localData,client->intraComm);
613  CArray<size_t,1> neighborInd(nbNeighbor);
614  for (int idx = 0; idx < nbNeighbor; ++idx)
615    neighborInd(idx) = neighborsDomainSrc(0,idx);
616
617  // Compute local data on other domains
618  dhtData.computeIndexInfoMapping(neighborInd);
619  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
620  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
621  // Ok get neighbor data
622  size_t nIdx;
623  int nbUnMaskedPointOnExtendedPart = 0;
624  for (int idx = 0; idx < nbNeighbor; ++idx)
625  {
626    nIdx  = neighborInd(idx);
627    it = neighborData.find(nIdx);
628    if (ite != it)
629    {
630      index = idx + oldNbLocal;
631      std::vector<double>& data = it->second;
632      // Unpack data
633      int dataIdx = 0;
634      lon_dst(index) = data[dataIdx]; ++dataIdx;
635      lat_dst(index) = data[dataIdx]; ++dataIdx;
636      for (int i = 0; i < nVertex; ++i)
637      {
638        bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx;
639      }
640      for (int i = 0; i < nVertex; ++i)
641      {
642        bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx;
643      }
644      mask_1d_dst(index) = (1.0 == data[dataIdx]) ? true : false; ++dataIdx;
645      data_i_index_dst_full(index) = (int)(data[dataIdx]);
646      if (0 <= data_i_index_dst_full(index))
647      {
648        data_i_index_dst_full(index) = index;
649        ++nbUnMaskedPointOnExtendedPart;
650      }
651    }
652  }
653
654
655  // Finally, update data_i_index
656  int nbUnMaskedPointOnNewDstDomain = (nbUnMaskedPointOnExtendedPart + nbUnMaskedPointOnLocalDomain);
657  int count = 0, dataIdx;
658  for (int idx = 0; idx < newNbLocalDst; ++idx)
659  {
660    dataIdx = data_i_index_dst_full(idx);
661    if ((0 <= dataIdx))
662    {
663      ++count;
664    }
665  }
666
667  data_i_index_dst.resize(count);
668  data_j_index_dst.resize(count);
669  data_j_index_dst = 0;
670
671  count = 0;
672  for (int idx = 0; idx < newNbLocalDst; ++idx)
673  {
674    dataIdx = data_i_index_dst_full(idx);
675    if ((0 <= dataIdx))
676    {
677      data_i_index_dst(count) = dataIdx;
678      ++count;
679    }
680  }
681
682  // Update ni
683  domainDestination->ni.setValue(newNbLocalDst);
684  domainDestination->mask_1d.resize(domainDestination->domainMask.numElements()) ;
685  domainDestination->mask_1d=domainDestination->domainMask ;
686  domainDestination->computeLocalMask() ;
687}
688CATCH
689
690/*!
691  Compute the index mapping between domain on grid source and one on grid destination
692*/
693void CDomainAlgorithmExpand::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
694{
695
696}
697
698}
Note: See TracBrowser for help on using the repository browser.