source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/domain_algorithm_expand.cpp @ 1260

Last change on this file since 1260 was 1260, checked in by ymipsl, 3 years ago

Buf fix in reduction. Missing value update was not set correctly

YM

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