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

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

XIOS coupling branch :
Blitz array resizeAndPreserve crash when size=0 => use resize instead

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