source: XIOS/dev/dev_olga/src/transformation/domain_algorithm_expand.cpp @ 1620

Last change on this file since 1620 was 1612, checked in by oabramkina, 22 months ago

Dev: adding exception handling.

To activate it, compilation flag -DXIOS_EXCEPTION should be added.

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