source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.cpp @ 1985

Last change on this file since 1985 was 1985, checked in by ymipsl, 4 years ago

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