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

Last change on this file since 1078 was 1078, checked in by mhnguyen, 5 years ago

Adding rectilinear and curvilinear domain for expand_domain transformation
-) Rectilinear/curvilinear is expanded not only locally but also globally, its global size ni_glo, nj_glo become ni_glo+2 and nj_glo+2
-) Two attributes i_periodic, j_periodic are only used for rectilinear/curvilinear to process priodic condition

+) Do some minor modification

Test
+) Add test_connectivity_expand
+) On Curie
+) Work (but need more real tests)

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