source: XIOS/dev/dev_ym/XIOS_COUPLING/src/node/domain.cpp @ 1918

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

Big update on on going work related to data distribution and transfer between clients and servers.

  • move all related file into distribution directorie
  • implement the concept of data "View"
  • implement the concept of "connector" which make the data transfer between 2 differents "Views"

YM

  • Property copyright set to
    Software name : XIOS (Xml I/O Server)
    http://forge.ipsl.jussieu.fr/ioserver
    Creation date : January 2009
    Licence : CeCCIL version2
    see license file in root directory : Licence_CeCILL_V2-en.txt
    or http://www.cecill.info/licences/Licence_CeCILL_V2-en.html
    Holder : CEA/LSCE (Laboratoire des Sciences du CLimat et de l'Environnement)
    CNRS/IPSL (Institut Pierre Simon Laplace)
    Project Manager : Yann Meurdesoif
    yann.meurdesoif@cea.fr
  • Property svn:executable set to *
File size: 120.0 KB
Line 
1#include "domain.hpp"
2
3#include "attribute_template.hpp"
4#include "object_template.hpp"
5#include "group_template.hpp"
6
7#include "xios_spl.hpp"
8#include "event_client.hpp"
9#include "event_server.hpp"
10#include "buffer_in.hpp"
11#include "message.hpp"
12#include "type.hpp"
13#include "context.hpp"
14#include "context_client.hpp"
15#include "context_server.hpp"
16#include "array_new.hpp"
17#include "distribution_client.hpp"
18#include "server_distribution_description.hpp"
19#include "client_server_mapping_distributed.hpp"
20#include "local_connector.hpp"
21#include "grid_local_connector.hpp"
22#include "remote_connector.hpp"
23#include "gatherer_connector.hpp"
24#include "scatterer_connector.hpp"
25#include "grid_scatterer_connector.hpp"
26#include "grid_gatherer_connector.hpp"
27
28
29
30
31#include <algorithm>
32
33namespace xios {
34
35   /// ////////////////////// Définitions ////////////////////// ///
36
37   CDomain::CDomain(void)
38      : CObjectTemplate<CDomain>(), CDomainAttributes()
39      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_()
40      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
41      , isClientAfterTransformationChecked(false), hasLonLat(false)
42      , isRedistributed_(false), hasPole(false)
43      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
44      , globalLocalIndexMap_(), computedWrittenIndex_(false)
45      , clients(), hasLatInReadFile_(false), hasBoundsLatInReadFile_(false)
46      , hasLonInReadFile_(false), hasBoundsLonInReadFile_(false)
47   {
48   }
49
50   CDomain::CDomain(const StdString & id)
51      : CObjectTemplate<CDomain>(id), CDomainAttributes()
52      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_() 
53      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
54      , isClientAfterTransformationChecked(false), hasLonLat(false)
55      , isRedistributed_(false), hasPole(false)
56      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
57      , globalLocalIndexMap_(), computedWrittenIndex_(false)
58      , clients(), hasLatInReadFile_(false), hasBoundsLatInReadFile_(false)
59      , hasLonInReadFile_(false), hasBoundsLonInReadFile_(false)
60   {
61    }
62
63   CDomain::~CDomain(void)
64   {
65   }
66
67   ///---------------------------------------------------------------
68
69   void CDomain::assignMesh(const StdString meshName, const int nvertex)
70   TRY
71   {
72     mesh = CMesh::getMesh(meshName, nvertex);
73   }
74   CATCH_DUMP_ATTR
75
76   CDomain* CDomain::createDomain()
77   TRY
78   {
79     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
80     return domain;
81   }
82   CATCH
83
84   std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
85   bool CDomain::_dummyTransformationMapList = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
86
87   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
88   TRY
89   {
90     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
91     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
92     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
93     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
94     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
95     m["reorder_domain"] = TRANS_REORDER_DOMAIN;
96     m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
97   }
98   CATCH
99
100   const std::set<StdString> & CDomain::getRelFiles(void) const
101   TRY
102   {
103      return (this->relFiles);
104   }
105   CATCH
106
107   /*!
108     Returns the number of indexes written by each server.
109     \return the number of indexes written by each server
110   */
111   int CDomain::getNumberWrittenIndexes(MPI_Comm writtenCom)
112   TRY
113   {
114     int writtenSize;
115     MPI_Comm_size(writtenCom, &writtenSize);
116     return numberWrittenIndexes_[writtenSize];
117   }
118   CATCH_DUMP_ATTR
119
120   /*!
121     Returns the total number of indexes written by the servers.
122     \return the total number of indexes written by the servers
123   */
124   int CDomain::getTotalNumberWrittenIndexes(MPI_Comm writtenCom)
125   TRY
126   {
127     int writtenSize;
128     MPI_Comm_size(writtenCom, &writtenSize);
129     return totalNumberWrittenIndexes_[writtenSize];
130   }
131   CATCH_DUMP_ATTR
132
133   /*!
134     Returns the offset of indexes written by each server.
135     \return the offset of indexes written by each server
136   */
137   int CDomain::getOffsetWrittenIndexes(MPI_Comm writtenCom)
138   TRY
139   {
140     int writtenSize;
141     MPI_Comm_size(writtenCom, &writtenSize);
142     return offsetWrittenIndexes_[writtenSize];
143   }
144   CATCH_DUMP_ATTR
145
146   CArray<int, 1>& CDomain::getCompressedIndexToWriteOnServer(MPI_Comm writtenCom)
147   TRY
148   {
149     int writtenSize;
150     MPI_Comm_size(writtenCom, &writtenSize);
151     return compressedIndexToWriteOnServer[writtenSize];
152   }
153   CATCH_DUMP_ATTR
154
155   //----------------------------------------------------------------
156
157   /*!
158    * Compute the minimum buffer size required to send the attributes to the server(s).
159    *
160    * \return A map associating the server rank with its minimum buffer size.
161    */
162   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
163   TRY
164   {
165
166     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
167
168     if (client->isServerLeader())
169     {
170       // size estimation for sendDistributionAttribut
171       size_t size = 11 * sizeof(size_t);
172
173       const std::list<int>& ranks = client->getRanksServerLeader();
174       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
175       {
176         if (size > attributesSizes[*itRank])
177           attributesSizes[*itRank] = size;
178       }
179     }
180
181     std::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->serverSize].end();
182     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
183     for (size_t k = 0; k < connectedServerRank_[client->serverSize].size(); ++k)
184     {
185       int rank = connectedServerRank_[client->serverSize][k];
186       std::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->serverSize].find(rank);
187       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
188
189       // size estimation for sendIndex (and sendArea which is always smaller or equal)
190       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
191
192       // size estimation for sendLonLat
193       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
194       if (hasBounds)
195         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
196
197       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
198       if (size > attributesSizes[rank])
199         attributesSizes[rank] = size;
200     }
201
202     return attributesSizes;
203   }
204   CATCH_DUMP_ATTR
205
206   //----------------------------------------------------------------
207
208   bool CDomain::isEmpty(void) const
209   TRY
210   {
211     return ((this->i_index.isEmpty()) || (0 == this->i_index.numElements()));
212   }
213   CATCH
214
215   //----------------------------------------------------------------
216
217   bool CDomain::IsWritten(const StdString & filename) const
218   TRY
219   {
220      return (this->relFiles.find(filename) != this->relFiles.end());
221   }
222   CATCH
223
224   bool CDomain::isWrittenCompressed(const StdString& filename) const
225   TRY
226   {
227      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
228   }
229   CATCH
230
231   //----------------------------------------------------------------
232
233   bool CDomain::isDistributed(void) const
234   TRY
235   {
236      bool distributed =  !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
237              (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
238      bool distributed_glo ;
239      distributed |= (1 == CContext::getCurrent()->intraCommSize_);
240
241      return distributed;
242   }
243   CATCH
244
245   //----------------------------------------------------------------
246
247   /*!
248    * Test whether the data defined on the domain can be outputted in a compressed way.
249    *
250    * \return true if and only if a mask was defined for this domain
251    */
252   bool CDomain::isCompressible(void) const
253   TRY
254   {
255      return isCompressible_;
256   }
257   CATCH
258
259   void CDomain::addRelFile(const StdString & filename)
260   TRY
261   {
262      this->relFiles.insert(filename);
263   }
264   CATCH_DUMP_ATTR
265
266   void CDomain::addRelFileCompressed(const StdString& filename)
267   TRY
268   {
269      this->relFilesCompressed.insert(filename);
270   }
271   CATCH_DUMP_ATTR
272
273   StdString CDomain::GetName(void)   { return (StdString("domain")); }
274   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
275   ENodeType CDomain::GetType(void)   { return (eDomain); }
276
277   //----------------------------------------------------------------
278
279   /*!
280      Verify if all distribution information of a domain are available
281      This checking verifies the definition of distribution attributes (ni, nj, ibegin, jbegin)
282   */
283   bool CDomain::distributionAttributesHaveValue() const
284   TRY
285   {
286      bool hasValues = true;
287
288      if (ni.isEmpty() && ibegin.isEmpty() && i_index.isEmpty())
289      {
290        hasValues = false;
291        return hasValues;
292      }
293
294      return hasValues;
295   }
296   CATCH
297
298   /*!
299     Redistribute RECTILINEAR or CURVILINEAR domain with a number of local domains.
300   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
301   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
302   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
303    \param [in] nbLocalDomain number of local domain on the domain destination
304   */
305   void CDomain::redistribute(int nbLocalDomain)
306   TRY
307   {
308     if (this->isRedistributed_) return;
309
310     this->isRedistributed_ = true;
311     CContext* context = CContext::getCurrent();
312     // For now the assumption is that secondary server pools consist of the same number of procs.
313     // CHANGE the line below if the assumption changes.
314
315     int rankClient = context->intraCommRank_;
316     int rankOnDomain = rankClient%nbLocalDomain;
317
318     if (ni_glo.isEmpty() || ni_glo <= 0 )
319     {
320        ERROR("CDomain::redistribute(int nbLocalDomain)",
321           << "[ Id = " << this->getId() << " ] "
322           << "The global domain is badly defined,"
323           << " check the \'ni_glo\'  value !")
324     }
325
326     if (nj_glo.isEmpty() || nj_glo <= 0 )
327     {
328        ERROR("CDomain::redistribute(int nbLocalDomain)",
329           << "[ Id = " << this->getId() << " ] "
330           << "The global domain is badly defined,"
331           << " check the \'nj_glo\'  value !")
332     }
333
334     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
335     {
336        int globalDomainSize = ni_glo * nj_glo;
337        if (globalDomainSize <= nbLocalDomain)
338        {
339          for (int idx = 0; idx < nbLocalDomain; ++idx)
340          {
341            if (rankOnDomain < globalDomainSize)
342            {
343              int iIdx = rankOnDomain % ni_glo;
344              int jIdx = rankOnDomain / ni_glo;
345              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
346              ni.setValue(1); nj.setValue(1);
347            }
348            else
349            {
350              ibegin.setValue(0); jbegin.setValue(0);
351              ni.setValue(0); nj.setValue(0);
352            }
353          }
354        }
355        else
356        {
357          float njGlo = nj_glo.getValue();
358          float niGlo = ni_glo.getValue();
359          int nbProcOnX, nbProcOnY, range;
360
361          // Compute (approximately) number of segment on x and y axis
362          float yOverXRatio = njGlo/niGlo;
363
364          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
365          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
366
367          // Simple distribution: Sweep from top to bottom, left to right
368          // Calculate local begin on x
369          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
370          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
371          for (int i = 1; i < nbProcOnX; ++i)
372          {
373            range = ni_glo / nbProcOnX;
374            if (i < (ni_glo%nbProcOnX)) ++range;
375            niVec[i-1] = range;
376            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
377          }
378          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
379
380          // Calculate local begin on y
381          for (int j = 1; j < nbProcOnY; ++j)
382          {
383            range = nj_glo / nbProcOnY;
384            if (j < (nj_glo%nbProcOnY)) ++range;
385            njVec[j-1] = range;
386            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
387          }
388          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
389
390          // Now assign value to ni, ibegin, nj, jbegin
391          int iIdx = rankOnDomain % nbProcOnX;
392          int jIdx = rankOnDomain / nbProcOnX;
393
394          if (rankOnDomain != (nbLocalDomain-1))
395          {
396            ibegin.setValue(ibeginVec[iIdx]);
397            jbegin.setValue(jbeginVec[jIdx]);
398            nj.setValue(njVec[jIdx]);
399            ni.setValue(niVec[iIdx]);
400          }
401          else // just merge all the remaining rectangle into the last one
402          {
403            ibegin.setValue(ibeginVec[iIdx]);
404            jbegin.setValue(jbeginVec[jIdx]);
405            nj.setValue(njVec[jIdx]);
406            ni.setValue(ni_glo - ibeginVec[iIdx]);
407          }
408        } 
409     }
410     else  // unstructured domain
411     {
412       if (this->i_index.isEmpty())
413       {
414          int globalDomainSize = ni_glo * nj_glo;
415          if (globalDomainSize <= nbLocalDomain)
416          {
417            for (int idx = 0; idx < nbLocalDomain; ++idx)
418            {
419              if (rankOnDomain < globalDomainSize)
420              {
421                int iIdx = rankOnDomain % ni_glo;
422                int jIdx = rankOnDomain / ni_glo;
423                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
424                ni.setValue(1); nj.setValue(1);
425              }
426              else
427              {
428                ibegin.setValue(0); jbegin.setValue(0);
429                ni.setValue(0); nj.setValue(0);
430              }
431            }
432          }
433          else
434          {
435            float njGlo = nj_glo.getValue();
436            float niGlo = ni_glo.getValue();
437            std::vector<int> ibeginVec(nbLocalDomain,0);
438            std::vector<int> niVec(nbLocalDomain);
439            for (int i = 1; i < nbLocalDomain; ++i)
440            {
441              int range = ni_glo / nbLocalDomain;
442              if (i < (ni_glo%nbLocalDomain)) ++range;
443              niVec[i-1] = range;
444              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
445            }
446            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
447
448            int iIdx = rankOnDomain % nbLocalDomain;
449            ibegin.setValue(ibeginVec[iIdx]);
450            jbegin.setValue(0);
451            ni.setValue(niVec[iIdx]);
452            nj.setValue(1);
453          }
454
455          i_index.resize(ni);         
456          for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
457        }
458        else
459        {
460          ibegin.setValue(this->i_index(0));
461          jbegin.setValue(0);
462          ni.setValue(this->i_index.numElements());
463          nj.setValue(1);
464        }
465     }
466
467     checkDomain();
468   }
469   CATCH_DUMP_ATTR
470
471   /*!
472     Fill in longitude and latitude whose values are read from file
473   */
474   void CDomain::fillInLonLat()
475   TRY
476   {
477     switch (type)
478     {
479      case type_attr::rectilinear:
480        fillInRectilinearLonLat();
481        break;
482      case type_attr::curvilinear:
483        fillInCurvilinearLonLat();
484        break;
485      case type_attr::unstructured:
486        fillInUnstructuredLonLat();
487        break;
488
489      default:
490      break;
491     }
492     completeLonLatClient() ;
493   }
494   CATCH_DUMP_ATTR
495
496   /*!
497     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
498     Range of longitude value from 0 - 360
499     Range of latitude value from -90 - +90
500   */
501   void CDomain::fillInRectilinearLonLat()
502   TRY
503   {
504     if (!lonvalue_rectilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
505     {
506       lonvalue_1d.resize(ni);
507       for (int idx = 0; idx < ni; ++idx)
508         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
509       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
510       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
511     }
512     else if (!hasLonInReadFile_)
513     {
514       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
515       lonvalue_1d.resize(ni);
516       double lonRange = lon_end - lon_start;
517       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
518
519        // Assign lon value
520       for (int i = 0; i < ni; ++i)
521       {
522         if (0 == (ibegin + i))
523         {
524           lonvalue_1d(i) = lon_start;
525         }
526         else if (ni_glo == (ibegin + i + 1))
527         {
528           lonvalue_1d(i) = lon_end;
529         }
530         else
531         {
532           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
533         }
534       }
535     }
536
537
538     if (!latvalue_rectilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
539     {
540       latvalue_1d.resize(nj);
541       for (int idx = 0; idx < nj; ++idx)
542         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
543       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
544       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
545     }
546     else if (!hasLatInReadFile_)
547     {
548       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
549       latvalue_1d.resize(nj);
550
551       double latRange = lat_end - lat_start;
552       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
553
554       for (int j = 0; j < nj; ++j)
555       {
556         if (0 == (jbegin + j))
557         {
558            latvalue_1d(j) = lat_start;
559         }
560         else if (nj_glo == (jbegin + j + 1))
561         {
562            latvalue_1d(j) = lat_end;
563         }
564         else
565         {
566           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
567         }
568       }
569     }
570   }
571   CATCH_DUMP_ATTR
572
573    /*
574      Fill in 2D longitude and latitude of curvilinear domain read from a file.
575      If there are already longitude and latitude defined by model. We just ignore read value.
576    */
577   void CDomain::fillInCurvilinearLonLat()
578   TRY
579   {
580     if (!lonvalue_curvilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
581     {
582       lonvalue_2d.resize(ni,nj);
583       for (int jdx = 0; jdx < nj; ++jdx)
584        for (int idx = 0; idx < ni; ++idx)
585         lonvalue_2d(idx,jdx) = lonvalue_curvilinear_read_from_file(idx, jdx);
586
587       lonvalue_curvilinear_read_from_file.free();
588     }
589
590     if (!latvalue_curvilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
591     {
592       latvalue_2d.resize(ni,nj);
593       for (int jdx = 0; jdx < nj; ++jdx)
594        for (int idx = 0; idx < ni; ++idx)
595           latvalue_2d(idx,jdx) = latvalue_curvilinear_read_from_file(idx, jdx);
596
597       latvalue_curvilinear_read_from_file.free();
598     }
599
600     if (!bounds_lonvalue_curvilinear_read_from_file.isEmpty() && bounds_lon_2d.isEmpty() && bounds_lon_1d.isEmpty())
601     {
602       bounds_lon_2d.resize(nvertex,ni,nj);
603       for (int jdx = 0; jdx < nj; ++jdx)
604        for (int idx = 0; idx < ni; ++idx)
605          for (int ndx = 0; ndx < nvertex; ++ndx)
606            bounds_lon_2d(ndx,idx,jdx) = bounds_lonvalue_curvilinear_read_from_file(ndx,idx, jdx);
607
608       bounds_lonvalue_curvilinear_read_from_file.free();
609     }
610
611     if (!bounds_latvalue_curvilinear_read_from_file.isEmpty() && bounds_lat_2d.isEmpty() && bounds_lat_1d.isEmpty())
612     {
613       bounds_lat_2d.resize(nvertex,ni,nj);
614       for (int jdx = 0; jdx < nj; ++jdx)
615        for (int idx = 0; idx < ni; ++idx)
616          for (int ndx = 0; ndx < nvertex; ++ndx)
617            bounds_lat_2d(ndx,idx,jdx) = bounds_latvalue_curvilinear_read_from_file(ndx,idx, jdx);
618
619       bounds_latvalue_curvilinear_read_from_file.free();
620     }
621   }
622   CATCH_DUMP_ATTR
623
624    /*
625      Fill in longitude and latitude of unstructured domain read from a file
626      If there are already longitude and latitude defined by model. We just igonore reading value.
627    */
628   void CDomain::fillInUnstructuredLonLat()
629   TRY
630   {
631     if (i_index.isEmpty())
632     {
633       i_index.resize(ni);
634       for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
635     }
636
637     if (!lonvalue_unstructured_read_from_file.isEmpty() && lonvalue_1d.isEmpty())
638     {
639        lonvalue_1d.resize(ni);
640        for (int idx = 0; idx < ni; ++idx)
641          lonvalue_1d(idx) = lonvalue_unstructured_read_from_file(idx);
642
643        // We dont need these values anymore, so just delete them
644        lonvalue_unstructured_read_from_file.free();
645     } 
646
647     if (!latvalue_unstructured_read_from_file.isEmpty() && latvalue_1d.isEmpty())
648     {
649        latvalue_1d.resize(ni);
650        for (int idx = 0; idx < ni; ++idx)
651          latvalue_1d(idx) =  latvalue_unstructured_read_from_file(idx);
652
653        // We dont need these values anymore, so just delete them
654        latvalue_unstructured_read_from_file.free();
655     }
656
657     if (!bounds_lonvalue_unstructured_read_from_file.isEmpty() && bounds_lon_1d.isEmpty())
658     {
659        int nbVertex = nvertex;
660        bounds_lon_1d.resize(nbVertex,ni);
661        for (int idx = 0; idx < ni; ++idx)
662          for (int jdx = 0; jdx < nbVertex; ++jdx)
663            bounds_lon_1d(jdx,idx) = bounds_lonvalue_unstructured_read_from_file(jdx, idx);
664
665        // We dont need these values anymore, so just delete them
666        bounds_lonvalue_unstructured_read_from_file.free();
667     }
668
669     if (!bounds_latvalue_unstructured_read_from_file.isEmpty() && bounds_lat_1d.isEmpty())
670     {
671        int nbVertex = nvertex;
672        bounds_lat_1d.resize(nbVertex,ni);
673        for (int idx = 0; idx < ni; ++idx)
674          for (int jdx = 0; jdx < nbVertex; ++jdx)
675            bounds_lat_1d(jdx,idx) = bounds_latvalue_unstructured_read_from_file(jdx, idx);
676
677        // We dont need these values anymore, so just delete them
678        bounds_latvalue_unstructured_read_from_file.free();
679     }
680   }
681   CATCH_DUMP_ATTR
682
683  /*
684    Get global longitude and latitude of rectilinear domain.
685  */
686   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g)
687   TRY
688   {
689     CContext* context = CContext::getCurrent();
690    // For now the assumption is that secondary server pools consist of the same number of procs.
691    // CHANGE the line below if the assumption changes.
692     int clientSize = context->intraCommSize_ ;
693     lon_g.resize(ni_glo) ;
694     lat_g.resize(nj_glo) ;
695
696
697     int* ibegin_g = new int[clientSize] ;
698     int* jbegin_g = new int[clientSize] ;
699     int* ni_g = new int[clientSize] ;
700     int* nj_g = new int[clientSize] ;
701     int v ;
702     v=ibegin ;
703     MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,context->intraComm_) ;
704     v=jbegin ;
705     MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,context->intraComm_) ;
706     v=ni ;
707     MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,context->intraComm_) ;
708     v=nj ;
709     MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,context->intraComm_) ;
710
711     MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,context->intraComm_) ;
712     MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,context->intraComm_) ;
713
714      delete[] ibegin_g ;
715      delete[] jbegin_g ;
716      delete[] ni_g ;
717      delete[] nj_g ;
718   }
719   CATCH_DUMP_ATTR
720
721   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
722                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
723   TRY
724  {
725     int i,j,k;
726
727     const int nvertexValue = 4;
728     boundsLon.resize(nvertexValue,ni*nj);
729
730     if (ni_glo>1)
731     {
732       double lonStepStart = lon(1)-lon(0);
733       bounds_lon_start=lon(0) - lonStepStart/2;
734       double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2);
735       bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
736       double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
737
738       // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
739       if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
740       {
741         bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
742         bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
743       }
744     }
745     else
746     {
747       if (bounds_lon_start.isEmpty()) bounds_lon_start=-180. ;
748       if (bounds_lon_end.isEmpty()) bounds_lon_end=180.-1e-8 ;
749     }
750
751     for(j=0;j<nj;++j)
752       for(i=0;i<ni;++i)
753       {
754         k=j*ni+i;
755         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start
756                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2;
757         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end
758                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2;
759       }
760
761
762    boundsLat.resize(nvertexValue,nj*ni);
763    bool isNorthPole=false ;
764    bool isSouthPole=false ;
765    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
766    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
767
768    // lat boundaries beyond pole the assimilate it to pole
769    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
770    if (nj_glo>1)
771    {
772      double latStepStart = lat(1)-lat(0);
773      if (isNorthPole) bounds_lat_start=lat(0);
774      else
775      {
776        bounds_lat_start=lat(0)-latStepStart/2;
777        if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
778        else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
779        else if (bounds_lat_start <= 90 && bounds_lat_start >= lat(0))
780        {
781          if ( std::abs(90-bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
782        }
783        else if (bounds_lat_start >= -90 && bounds_lat_start <= lat(0))
784        {
785          if ( std::abs(-90 - bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
786        }
787      }
788
789      double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
790      if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
791      else
792      {
793        bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
794
795        if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
796        else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
797        else if (bounds_lat_end <= 90 && bounds_lat_end >= lat(nj_glo-1))
798        {
799          if ( std::abs(90-bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
800        }
801        else if (bounds_lat_end >= -90 && bounds_lat_end <= lat(nj_glo-1))
802        {
803          if ( std::abs(-90 - bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
804        }
805      }
806    }
807    else
808    {
809      if (bounds_lat_start.isEmpty()) bounds_lat_start=-90. ;
810      if (bounds_lat_end.isEmpty()) bounds_lat_end=90 ;
811    }
812
813    for(j=0;j<nj;++j)
814      for(i=0;i<ni;++i)
815      {
816        k=j*ni+i;
817        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start
818                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2;
819        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end
820                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2;
821      }
822   }
823   CATCH_DUMP_ATTR
824
825   /*
826     General check of the domain to verify its mandatory attributes
827   */
828   void CDomain::checkDomain(void)
829   TRY
830   {
831     if (type.isEmpty())
832     {
833       ERROR("CDomain::checkDomain(void)",
834             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
835             << "The domain type is mandatory, "
836             << "please define the 'type' attribute.")
837     }
838
839     if (type == type_attr::gaussian) 
840     {
841        hasPole=true ;
842        type.setValue(type_attr::unstructured) ;
843      }
844      else if (type == type_attr::rectilinear) hasPole=true ;
845
846     if (type == type_attr::unstructured)
847     {
848        if (ni_glo.isEmpty())
849        {
850          ERROR("CDomain::checkDomain(void)",
851                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
852                << "The global domain is badly defined, "
853                << "the mandatory 'ni_glo' attribute is missing.")
854        }
855        else if (ni_glo <= 0)
856        {
857          ERROR("CDomain::checkDomain(void)",
858                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
859                << "The global domain is badly defined, "
860                << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
861        }
862        isUnstructed_ = true;
863        nj_glo = 1;
864        nj = 1;
865        jbegin = 0;
866        if (!i_index.isEmpty()) ni = i_index.numElements();
867        j_index.resize(ni);
868        for(int i=0;i<ni;++i) j_index(i)=0;
869
870        if (!area.isEmpty())
871          area.transposeSelf(1, 0);
872     }
873
874     if (ni_glo.isEmpty())
875     {
876       ERROR("CDomain::checkDomain(void)",
877             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
878             << "The global domain is badly defined, "
879             << "the mandatory 'ni_glo' attribute is missing.")
880     }
881     else if (ni_glo <= 0)
882     {
883       ERROR("CDomain::checkDomain(void)",
884             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
885             << "The global domain is badly defined, "
886             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
887     }
888
889     if (nj_glo.isEmpty())
890     {
891       ERROR("CDomain::checkDomain(void)",
892             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
893             << "The global domain is badly defined, "
894             << "the mandatory 'nj_glo' attribute is missing.")
895     }
896     else if (nj_glo <= 0)
897     {
898       ERROR("CDomain::checkDomain(void)",
899             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
900             << "The global domain is badly defined, "
901             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
902     }
903
904     checkLocalIDomain();
905     checkLocalJDomain();
906
907     if (i_index.isEmpty())
908     {
909       i_index.resize(ni*nj);
910       for (int j = 0; j < nj; ++j)
911         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
912     }
913
914     if (j_index.isEmpty())
915     {
916       j_index.resize(ni*nj);
917       for (int j = 0; j < nj; ++j)
918         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
919     }
920   }
921   CATCH_DUMP_ATTR
922
923   size_t CDomain::getGlobalWrittenSize(void)
924   {
925     return ni_glo*nj_glo ;
926   }
927   //----------------------------------------------------------------
928
929   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
930   void CDomain::checkLocalIDomain(void)
931   TRY
932   {
933      // If ibegin and ni are provided then we use them to check the validity of local domain
934      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
935      {
936        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
937        {
938          ERROR("CDomain::checkLocalIDomain(void)",
939                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
940                << "The local domain is wrongly defined,"
941                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
942        }
943      }
944
945      // i_index has higher priority than ibegin and ni
946      if (!i_index.isEmpty())
947      {
948        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
949        if (ni.isEmpty()) 
950        {         
951         // No information about ni
952          int minIndex = ni_glo - 1;
953          int maxIndex = 0;
954          for (int idx = 0; idx < i_index.numElements(); ++idx)
955          {
956            if (i_index(idx) < minIndex) minIndex = i_index(idx);
957            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
958          }
959          ni = maxIndex - minIndex + 1; 
960          minIIndex = minIIndex;         
961        }
962
963        // It's not so correct but if ibegin is not the first value of i_index
964        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
965        if (ibegin.isEmpty()) ibegin = minIIndex;
966      }
967      else if (ibegin.isEmpty() && ni.isEmpty())
968      {
969        ibegin = 0;
970        ni = ni_glo;
971      }
972      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
973      {
974        ERROR("CDomain::checkLocalIDomain(void)",
975              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
976              << "The local domain is wrongly defined," << endl
977              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
978              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
979      }
980       
981
982      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
983      {
984        ERROR("CDomain::checkLocalIDomain(void)",
985              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
986              << "The local domain is wrongly defined,"
987              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
988      }
989   }
990   CATCH_DUMP_ATTR
991
992   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
993   void CDomain::checkLocalJDomain(void)
994   TRY
995   {
996    // If jbegin and nj are provided then we use them to check the validity of local domain
997     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
998     {
999       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
1000       {
1001         ERROR("CDomain::checkLocalJDomain(void)",
1002                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1003                << "The local domain is wrongly defined,"
1004                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1005       }
1006     }
1007
1008     if (!j_index.isEmpty())
1009     {
1010        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
1011        if (nj.isEmpty()) 
1012        {
1013          // No information about nj
1014          int minIndex = nj_glo - 1;
1015          int maxIndex = 0;
1016          for (int idx = 0; idx < j_index.numElements(); ++idx)
1017          {
1018            if (j_index(idx) < minIndex) minIndex = j_index(idx);
1019            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
1020          }
1021          nj = maxIndex - minIndex + 1;
1022          minJIndex = minIndex; 
1023        } 
1024        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
1025        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
1026       if (jbegin.isEmpty()) jbegin = minJIndex;       
1027     }
1028     else if (jbegin.isEmpty() && nj.isEmpty())
1029     {
1030       jbegin = 0;
1031       nj = nj_glo;
1032     }     
1033
1034
1035     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
1036     {
1037       ERROR("CDomain::checkLocalJDomain(void)",
1038              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1039              << "The local domain is wrongly defined,"
1040              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1041     }
1042   }
1043   CATCH_DUMP_ATTR
1044
1045   //----------------------------------------------------------------
1046
1047   void CDomain::checkMask(void)
1048   TRY
1049   {
1050      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
1051        ERROR("CDomain::checkMask(void)",
1052              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1053              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
1054              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
1055
1056      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
1057      {
1058        if (mask_1d.numElements() != i_index.numElements())
1059          ERROR("CDomain::checkMask(void)",
1060                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1061                << "'mask_1d' does not have the same size as the local domain." << std::endl
1062                << "Local size is " << i_index.numElements() << "." << std::endl
1063                << "Mask size is " << mask_1d.numElements() << ".");
1064      }
1065
1066      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
1067      {
1068        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1069          ERROR("CDomain::checkMask(void)",
1070                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1071                << "The mask does not have the same size as the local domain." << std::endl
1072                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1073                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
1074      }
1075
1076      if (!mask_2d.isEmpty())
1077      {
1078        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
1079        for (int j = 0; j < nj; ++j)
1080          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1081//        mask_2d.reset();
1082      }
1083      else if (mask_1d.isEmpty())
1084      {
1085        domainMask.resize(i_index.numElements());
1086        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
1087      }
1088      else
1089      {
1090      domainMask.resize(mask_1d.numElements());
1091      domainMask=mask_1d ;
1092     }
1093   }
1094   CATCH_DUMP_ATTR
1095
1096   //----------------------------------------------------------------
1097
1098   void CDomain::checkDomainData(void)
1099   TRY
1100   {
1101      if (data_dim.isEmpty())
1102      {
1103        data_dim.setValue(1);
1104      }
1105      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
1106      {
1107        ERROR("CDomain::checkDomainData(void)",
1108              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1109              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
1110      }
1111
1112      if (data_ibegin.isEmpty())
1113         data_ibegin.setValue(0);
1114      if (data_jbegin.isEmpty())
1115         data_jbegin.setValue(0);
1116
1117      if (data_ni.isEmpty())
1118      {
1119        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
1120      }
1121      else if (data_ni.getValue() < 0)
1122      {
1123        ERROR("CDomain::checkDomainData(void)",
1124              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1125              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
1126      }
1127
1128      if (data_nj.isEmpty())
1129      {
1130        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
1131      }
1132      else if (data_nj.getValue() < 0)
1133      {
1134        ERROR("CDomain::checkDomainData(void)",
1135              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1136              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1137      }
1138   }
1139   CATCH_DUMP_ATTR
1140
1141   //----------------------------------------------------------------
1142
1143   void CDomain::checkCompression(void)
1144   TRY
1145   {
1146     int i,j,ind;
1147      if (!data_i_index.isEmpty())
1148      {
1149        if (!data_j_index.isEmpty() &&
1150            data_j_index.numElements() != data_i_index.numElements())
1151        {
1152           ERROR("CDomain::checkCompression(void)",
1153                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1154                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1155                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1156                 << "'data_j_index' size = " << data_j_index.numElements());
1157        }
1158
1159        if (2 == data_dim)
1160        {
1161          if (data_j_index.isEmpty())
1162          {
1163             ERROR("CDomain::checkCompression(void)",
1164                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1165                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
1166          }
1167          for (int k=0; k<data_i_index.numElements(); ++k)
1168          {
1169            i = data_i_index(k)+data_ibegin ;
1170            j = data_j_index(k)+data_jbegin ;
1171            if (i>=0 && i<ni && j>=0 && j<nj)
1172            {
1173              ind=j*ni+i ;
1174              if (!domainMask(ind))
1175              {
1176                data_i_index(k) = -1;
1177                data_j_index(k) = -1;
1178              }
1179            }
1180            else
1181            {
1182              data_i_index(k) = -1;
1183              data_j_index(k) = -1;
1184            }
1185          }
1186        }
1187        else // (1 == data_dim)
1188        {
1189          if (data_j_index.isEmpty())
1190          {
1191            data_j_index.resize(data_ni);
1192            data_j_index = 0;
1193          }
1194          for (int k=0; k<data_i_index.numElements(); ++k)
1195          {
1196            i=data_i_index(k)+data_ibegin ;
1197            if (i>=0 && i < domainMask.size())
1198            {
1199              if (!domainMask(i)) data_i_index(k) = -1;
1200            }
1201            else
1202              data_i_index(k) = -1;
1203
1204            if (!domainMask(i)) data_i_index(k) = -1;
1205          }
1206        }
1207      }
1208      else
1209      {
1210        if (data_dim == 2 && !data_j_index.isEmpty())
1211          ERROR("CDomain::checkCompression(void)",
1212                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1213                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
1214
1215        if (1 == data_dim)
1216        {
1217          data_i_index.resize(data_ni);
1218          data_j_index.resize(data_ni);
1219          data_j_index = 0;
1220
1221          for (int k = 0; k < data_ni; ++k)
1222          {
1223            i=k+data_ibegin ;
1224            if (i>=0 && i < domainMask.size())
1225            {
1226              if (domainMask(i))
1227                data_i_index(k) = k;
1228              else
1229                data_i_index(k) = -1;
1230            }
1231            else
1232              data_i_index(k) = -1;
1233          }
1234        }
1235        else // (data_dim == 2)
1236        {
1237          const int dsize = data_ni * data_nj;
1238          data_i_index.resize(dsize);
1239          data_j_index.resize(dsize);
1240
1241          for(int count = 0, kj = 0; kj < data_nj; ++kj)
1242          {
1243            for(int ki = 0; ki < data_ni; ++ki, ++count)
1244            {
1245              i = ki + data_ibegin;
1246              j = kj + data_jbegin;
1247              ind=j*ni+i ;
1248              if (i>=0 && i<ni && j>=0 && j<nj)
1249              {
1250                if (domainMask(ind))
1251                {
1252                  data_i_index(count) = ki;
1253                  data_j_index(count) = kj;
1254                }
1255                else
1256                {
1257                  data_i_index(count) = -1;
1258                  data_j_index(count) = -1;
1259                }
1260              }
1261              else
1262              {
1263                data_i_index(count) = -1;
1264                data_j_index(count) = -1;
1265              }
1266            }
1267          }
1268        }
1269      }
1270   }
1271   CATCH_DUMP_ATTR
1272
1273   //----------------------------------------------------------------
1274   void CDomain::computeLocalMask(void)
1275   TRY
1276   {
1277     localMask.resize(i_index.numElements()) ;
1278     localMask=false ;
1279
1280     size_t dn=data_i_index.numElements() ;
1281     int i,j ;
1282     size_t k,ind ;
1283
1284     for(k=0;k<dn;k++)
1285     {
1286       if (data_dim==2)
1287       {
1288          i=data_i_index(k)+data_ibegin ;
1289          j=data_j_index(k)+data_jbegin ;
1290          if (i>=0 && i<ni && j>=0 && j<nj)
1291          {
1292            ind=j*ni+i ;
1293            localMask(ind)=domainMask(ind) ;
1294          }
1295       }
1296       else
1297       {
1298          i=data_i_index(k)+data_ibegin ;
1299          if (i>=0 && i<i_index.numElements())
1300          {
1301            ind=i ;
1302            localMask(ind)=domainMask(ind) ;
1303          }
1304       }
1305     }
1306   }
1307   CATCH_DUMP_ATTR
1308
1309   void CDomain::checkEligibilityForCompressedOutput(void)
1310   TRY
1311   {
1312     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1313     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1314   }
1315   CATCH_DUMP_ATTR
1316
1317   //----------------------------------------------------------------
1318
1319   /*
1320     Fill in longitude, latitude, bounds, and area into internal values (lonvalue, latvalue, bounds_lonvalue, bounds_latvalue, areavalue)
1321     which will be used by XIOS.
1322   */
1323   void CDomain::completeLonLatClient(void)
1324   TRY
1325   {
1326     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1327     checkBounds() ;
1328     checkArea() ;
1329
1330     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
1331     {
1332       lonvalue.resize(ni * nj);
1333       latvalue.resize(ni * nj);
1334       if (hasBounds)
1335       {
1336         bounds_lonvalue.resize(nvertex, ni * nj);
1337         bounds_latvalue.resize(nvertex, ni * nj);
1338       }
1339
1340       for (int j = 0; j < nj; ++j)
1341       {
1342         for (int i = 0; i < ni; ++i)
1343         {
1344           int k = j * ni + i;
1345
1346           lonvalue(k) = lonvalue_2d(i,j);
1347           latvalue(k) = latvalue_2d(i,j);
1348
1349           if (hasBounds)
1350           {
1351             for (int n = 0; n < nvertex; ++n)
1352             {
1353               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1354               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
1355             }
1356           }
1357         }
1358       }
1359     }
1360     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
1361     {
1362       if (type_attr::rectilinear == type)
1363       {
1364         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1365         {
1366           lonvalue.resize(ni * nj);
1367           latvalue.resize(ni * nj);
1368           if (hasBounds)
1369           {
1370             bounds_lonvalue.resize(nvertex, ni * nj);
1371             bounds_latvalue.resize(nvertex, ni * nj);
1372           }
1373
1374           for (int j = 0; j < nj; ++j)
1375           {
1376             for (int i = 0; i < ni; ++i)
1377             {
1378               int k = j * ni + i;
1379
1380               lonvalue(k) = lonvalue_1d(i);
1381               latvalue(k) = latvalue_1d(j);
1382
1383               if (hasBounds)
1384               {
1385                 for (int n = 0; n < nvertex; ++n)
1386                 {
1387                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1388                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1389                 }
1390               }
1391             }
1392           }
1393         }
1394         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1395         {
1396           lonvalue.reference(lonvalue_1d);
1397           latvalue.reference(latvalue_1d);
1398            if (hasBounds)
1399           {
1400             bounds_lonvalue.reference(bounds_lon_1d);
1401             bounds_latvalue.reference(bounds_lat_1d);
1402           }
1403         }
1404         else
1405           ERROR("CDomain::completeLonClient(void)",
1406                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1407                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1408                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1409                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1410                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1411                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1412       }
1413       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1414       {
1415         lonvalue.reference(lonvalue_1d);
1416         latvalue.reference(latvalue_1d);
1417         if (hasBounds)
1418         {
1419           bounds_lonvalue.reference(bounds_lon_1d);
1420           bounds_latvalue.reference(bounds_lat_1d);
1421         }
1422       }
1423     }
1424
1425     if (!area.isEmpty() && areavalue.isEmpty())
1426     {
1427        areavalue.resize(ni*nj);
1428       for (int j = 0; j < nj; ++j)
1429       {
1430         for (int i = 0; i < ni; ++i)
1431         {
1432           int k = j * ni + i;
1433           areavalue(k) = area(i,j);
1434         }
1435       }
1436     }
1437   }
1438   CATCH_DUMP_ATTR
1439
1440   /*
1441     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1442   */
1443   void CDomain::convertLonLatValue(void)
1444   TRY
1445   {
1446     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1447     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1448     {
1449       lonvalue_2d.resize(ni,nj);
1450       latvalue_2d.resize(ni,nj);
1451       if (hasBounds)
1452       {
1453         bounds_lon_2d.resize(nvertex, ni, nj);
1454         bounds_lat_2d.resize(nvertex, ni, nj);
1455       }
1456
1457       for (int j = 0; j < nj; ++j)
1458       {
1459         for (int i = 0; i < ni; ++i)
1460         {
1461           int k = j * ni + i;
1462
1463           lonvalue_2d(i,j) = lonvalue(k);
1464           latvalue_2d(i,j) = latvalue(k);
1465
1466           if (hasBounds)
1467           {
1468             for (int n = 0; n < nvertex; ++n)
1469             {
1470               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1471               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1472             }
1473           }
1474         }
1475       }
1476     }
1477     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1478     {
1479       if (type_attr::rectilinear == type)
1480       {
1481         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1482         {
1483           lonvalue.resize(ni * nj);
1484           latvalue.resize(ni * nj);
1485           if (hasBounds)
1486           {
1487             bounds_lonvalue.resize(nvertex, ni * nj);
1488             bounds_latvalue.resize(nvertex, ni * nj);
1489           }
1490
1491           for (int j = 0; j < nj; ++j)
1492           {
1493             for (int i = 0; i < ni; ++i)
1494             {
1495               int k = j * ni + i;
1496
1497               lonvalue(k) = lonvalue_1d(i);
1498               latvalue(k) = latvalue_1d(j);
1499
1500               if (hasBounds)
1501               {
1502                 for (int n = 0; n < nvertex; ++n)
1503                 {
1504                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1505                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1506                 }
1507               }
1508             }
1509           }
1510         }
1511         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1512         {
1513           lonvalue.reference(lonvalue_1d);
1514           latvalue.reference(latvalue_1d);
1515            if (hasBounds)
1516           {
1517             bounds_lonvalue.reference(bounds_lon_1d);
1518             bounds_latvalue.reference(bounds_lat_1d);
1519           }
1520         }
1521         else
1522           ERROR("CDomain::completeLonClient(void)",
1523                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1524                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1525                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1526                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1527                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1528                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1529       }
1530       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1531       {
1532         lonvalue.reference(lonvalue_1d);
1533         latvalue.reference(latvalue_1d);
1534         if (hasBounds)
1535         {
1536           bounds_lonvalue.reference(bounds_lon_1d);
1537           bounds_latvalue.reference(bounds_lat_1d);
1538         }
1539       }
1540     }
1541   }
1542   CATCH_DUMP_ATTR
1543
1544   void CDomain::checkBounds(void)
1545   TRY
1546   {
1547     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1548     if (!nvertex.isEmpty() && nvertex > 0 && !hasBoundValues)
1549     {
1550       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1551         ERROR("CDomain::checkBounds(void)",
1552               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1553               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1554               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1555
1556       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1557         ERROR("CDomain::checkBounds(void)",
1558               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1559               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1560               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1561
1562       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1563       {
1564         ERROR("CDomain::checkBounds(void)",
1565               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1566               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1567               << "Please define either both attributes or none.");
1568       }
1569
1570       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1571       {
1572         ERROR("CDomain::checkBounds(void)",
1573               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1574               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1575               << "Please define either both attributes or none.");
1576       }
1577
1578       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1579         ERROR("CDomain::checkBounds(void)",
1580               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1581               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1582               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(0)
1583               << " but nvertex is " << nvertex.getValue() << ".");
1584
1585       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1586         ERROR("CDomain::checkBounds(void)",
1587               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1588               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1589               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(0)
1590               << " but nvertex is " << nvertex.getValue() << ".");
1591
1592       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1593         ERROR("CDomain::checkBounds(void)",
1594               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1595               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1596
1597       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1598         ERROR("CDomain::checkBounds(void)",
1599               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1600               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1601
1602       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1603         ERROR("CDomain::checkBounds(void)",
1604               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1605               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1606               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(0)
1607               << " but nvertex is " << nvertex.getValue() << ".");
1608
1609       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1610         ERROR("CDomain::checkBounds(void)",
1611               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1612               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1613               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(0)
1614               << " but nvertex is " << nvertex.getValue() << ".");
1615
1616       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1617         ERROR("CDomain::checkBounds(void)",
1618               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1619               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1620
1621       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1622         ERROR("CDomain::checkBounds(void)",
1623               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1624
1625       // In case of reading UGRID bounds values are not required
1626       hasBounds = (!bounds_lat_1d.isEmpty() || !bounds_lat_2d.isEmpty() );
1627     }
1628     else if (hasBoundValues)
1629     {
1630       hasBounds = true;       
1631     }
1632     else
1633     {
1634       hasBounds = false;
1635     }
1636   }
1637   CATCH_DUMP_ATTR
1638
1639   void CDomain::checkArea(void)
1640   TRY
1641   {
1642     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
1643     hasArea = !area.isEmpty();
1644     if (hasArea && !hasAreaValue)
1645     {
1646       if (area.extent(0) != ni || area.extent(1) != nj)
1647       {
1648         ERROR("CDomain::checkArea(void)",
1649               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1650               << "The area does not have the same size as the local domain." << std::endl
1651               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1652               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1653       }
1654//       if (areavalue.isEmpty())
1655//       {
1656//          areavalue.resize(ni*nj);
1657//         for (int j = 0; j < nj; ++j)
1658//         {
1659//           for (int i = 0; i < ni; ++i)
1660//           {
1661//             int k = j * ni + i;
1662//             areavalue(k) = area(i,j);
1663//           }
1664//         }
1665//       }
1666     }
1667   }
1668   CATCH_DUMP_ATTR
1669
1670   void CDomain::checkLonLat()
1671   TRY
1672   {
1673     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1674                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1675     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1676     if (hasLonLat && !hasLonLatValue)
1677     {
1678       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1679         ERROR("CDomain::checkLonLat()",
1680               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1681               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1682               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1683
1684       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1685       {
1686         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1687           ERROR("CDomain::checkLonLat()",
1688                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1689                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1690                 << "Local size is " << i_index.numElements() << "." << std::endl
1691                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1692       }
1693
1694       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1695       {
1696         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1697           ERROR("CDomain::checkLonLat()",
1698                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1699                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1700                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1701                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1702       }
1703
1704       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1705         ERROR("CDomain::checkLonLat()",
1706               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1707               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1708               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1709
1710       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1711       {
1712         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1713           ERROR("CDomain::checkLonLat()",
1714                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1715                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1716                 << "Local size is " << i_index.numElements() << "." << std::endl
1717                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1718       }
1719
1720       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1721       {
1722         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1723           ERROR("CDomain::checkLonLat()",
1724                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1725                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1726                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1727                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1728       }
1729     }
1730   }
1731   CATCH_DUMP_ATTR
1732
1733   void CDomain::checkAttributes(void)
1734   TRY
1735   {
1736      if (this->checkAttributes_done_) return;
1737      this->checkDomain();
1738      this->checkLonLat();
1739      this->checkBounds();
1740      this->checkArea();
1741      this->checkMask();
1742      this->checkDomainData();
1743      this->checkCompression();
1744      this->computeLocalMask() ;
1745      this->completeLonLatClient();
1746      this->initializeLocalElement() ;
1747      this->addFullView() ;
1748      this->addWorkflowView() ;
1749      this->addModelView() ;
1750      // testing ?
1751      CLocalView* local = localElement_->getView(CElementView::WORKFLOW) ;
1752      CLocalView* model = localElement_->getView(CElementView::MODEL) ;
1753
1754      CLocalConnector test1(model, local) ;
1755      test1.computeConnector() ;
1756      CLocalConnector test2(local, model) ;
1757      test2.computeConnector() ;
1758      CGridLocalConnector gridTest1(vector<CLocalConnector*>{&test1}) ;
1759      CGridLocalConnector gridTest2(vector<CLocalConnector*>{&test2}) ;
1760     
1761     
1762      CArray<int,1> out1 ;
1763      CArray<int,1> out2 ;
1764      test1.transfer(data_i_index,out1,-111) ;
1765      test2.transfer(out1,out2,-111) ;
1766     
1767      out1 = 0 ;
1768      out2 = 0 ;
1769      gridTest1.transfer(data_i_index,out1,-111) ;
1770      gridTest2.transfer(out1, out2,-111) ;
1771     
1772      this->checkAttributes_done_ = true;
1773   }
1774   CATCH_DUMP_ATTR
1775
1776
1777   void CDomain::initializeLocalElement(void)
1778   {
1779      // after checkDomain i_index and j_index of size (ni*nj)
1780      int nij = ni*nj ;
1781      CArray<size_t, 1> ij_index(ni*nj) ;
1782      for(int ij=0; ij<nij ; ij++) ij_index(ij) = i_index(ij)+j_index(ij)*ni_glo ;
1783      int rank = CContext::getCurrent()->getIntraCommRank() ;
1784      localElement_ = new CLocalElement(rank, ni_glo*nj_glo, ij_index) ;
1785   }
1786
1787   void CDomain::addFullView(void)
1788   {
1789      CArray<int,1> index(ni*nj) ;
1790      int nij=ni*nj ;
1791      for(int ij=0; ij<nij ; ij++) index(ij)=ij ;
1792      localElement_ -> addView(CElementView::FULL, index) ;
1793   }
1794
1795   void CDomain::addWorkflowView(void)
1796   {
1797     // information for workflow view is stored in localMask
1798     int nij=ni*nj ;
1799     int nMask=0 ;
1800     for(int ij=0; ij<nij ; ij++) if (localMask(ij)) nMask++ ;
1801     CArray<int,1> index(nMask) ;
1802
1803     nMask=0 ;
1804     for(int ij=0; ij<nij ; ij++) 
1805      if (localMask(ij))
1806      {
1807        index(nMask)=ij ;
1808        nMask++ ;
1809      }
1810      localElement_ -> addView(CElementView::WORKFLOW, index) ;
1811   }
1812
1813   void CDomain::addModelView(void)
1814   {
1815     // information for model view is stored in data_i_index/data_j_index
1816     // very weird, do not mix data_i_index and data_i_begin => in future only keep data_i_index
1817     int dataSize = data_i_index.numElements() ;
1818     CArray<int,1> index(dataSize) ;
1819     int i,j ;
1820     for(int k=0;k<dataSize;k++)
1821     {
1822        i=data_i_index(k)+data_ibegin ; // bad
1823        j=data_j_index(k)+data_jbegin ; // bad
1824        if (i>=0 && i<ni && j>=0 && j<nj) index(k)=i+j*ni ;
1825        else index(k)=-1 ;
1826     }
1827     localElement_->addView(CElementView::MODEL, index) ;
1828   }
1829       
1830   void CDomain::computeModelToWorkflowConnector(void)
1831   { 
1832     CLocalView* srcView=getLocalView(CElementView::MODEL) ;
1833     CLocalView* dstView=getLocalView(CElementView::WORKFLOW) ;
1834     modelToWorkflowConnector_ = new CLocalConnector(srcView, dstView); 
1835     modelToWorkflowConnector_->computeConnector() ;
1836   }
1837
1838
1839   //----------------------------------------------------------------
1840   // Divide function checkAttributes into 2 seperate ones
1841   // This function only checks all attributes of current domain
1842   void CDomain::checkAttributesOnClient()
1843   TRY
1844   {
1845     if (this->isClientChecked) return;
1846     CContext* context=CContext::getCurrent();
1847
1848      if (context->getServiceType()==CServicesManager::CLIENT)
1849      {
1850        this->checkDomain();
1851        this->checkBounds();
1852        this->checkArea();
1853        this->checkLonLat();
1854      }
1855
1856      if (context->getServiceType()==CServicesManager::CLIENT)
1857      { // Ct client uniquement
1858         this->checkMask();
1859         this->checkDomainData();
1860         this->checkCompression();
1861         this->computeLocalMask() ;
1862      }
1863      else
1864      { // Ct serveur uniquement
1865      }
1866
1867      this->isClientChecked = true;
1868   }
1869   CATCH_DUMP_ATTR
1870   
1871   // ym obselete, to be removed
1872   void CDomain::checkAttributesOnClientAfterTransformation()
1873   TRY
1874   {
1875     CContext* context=CContext::getCurrent() ;
1876
1877     if (this->isClientAfterTransformationChecked) return;
1878     if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1879     {
1880       // this->computeConnectedClients();
1881       if (hasLonLat)
1882         if (context->getServiceType()==CServicesManager::CLIENT)
1883           this->completeLonLatClient();
1884     }
1885
1886     this->isClientAfterTransformationChecked = true;
1887   }
1888   CATCH_DUMP_ATTR
1889
1890      // Send all checked attributes to server
1891   void CDomain::sendCheckedAttributes()
1892   TRY
1893   {
1894     if (!this->isClientChecked) checkAttributesOnClient();
1895     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1896     CContext* context=CContext::getCurrent() ;
1897
1898     if (this->isChecked) return;
1899     if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1900     {
1901       sendAttributes();
1902     }
1903     this->isChecked = true;
1904   }
1905   CATCH_DUMP_ATTR
1906
1907/* old version
1908   void CDomain::checkAttributes(void)
1909   TRY
1910   {
1911      if (this->isChecked) return;
1912      CContext* context=CContext::getCurrent() ;
1913
1914      this->checkDomain();
1915      this->checkLonLat();
1916      this->checkBounds();
1917      this->checkArea();
1918
1919      if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1920      { // Ct client uniquement
1921         this->checkMask();
1922         this->checkDomainData();
1923         this->checkCompression();
1924         this->computeLocalMask() ;
1925
1926      }
1927      else
1928      { // Ct serveur uniquement
1929      }
1930
1931      if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1932      {
1933        this->computeConnectedClients();
1934        this->completeLonLatClient();
1935      }
1936
1937      this->isChecked = true;
1938   }
1939   CATCH_DUMP_ATTR
1940*/
1941   
1942  /*!
1943     Compute the connection of a client to other clients to determine which clients to send attributes to.
1944     The sending clients are supposed to already know the distribution of receiving clients (In simple cases, it's band)
1945     The connection among clients is calculated by using global index.
1946     A client connects to other clients which holds the same global index as it.     
1947  */
1948  void CDomain::computeConnectedClients(CContextClient* client)
1949  TRY
1950  {
1951    if (computeConnectedClients_done_.count(client)!=0) return ;
1952    else computeConnectedClients_done_.insert(client) ;
1953   
1954    CContext* context=CContext::getCurrent() ;
1955   
1956    int nbServer = client->serverSize;
1957    int nbClient = client->clientSize;
1958    int rank     = client->clientRank;
1959       
1960    if (listNbServer_.count(nbServer) == 0)
1961    {
1962      listNbServer_.insert(nbServer) ;
1963 
1964      if (connectedServerRank_.find(nbServer) != connectedServerRank_.end())
1965      {
1966        nbSenders.erase(nbServer);
1967        connectedServerRank_.erase(nbServer);
1968      }
1969
1970      if (indSrv_.find(nbServer) == indSrv_.end())
1971      {
1972        int i,j,i_ind,j_ind, nbIndex=i_index.numElements();
1973        int globalIndexCount = i_index.numElements();
1974        // Fill in index
1975        CArray<size_t,1> globalIndexDomain(nbIndex);
1976        size_t globalIndex;
1977
1978        for (i = 0; i < nbIndex; ++i)
1979        {
1980          i_ind=i_index(i);
1981          j_ind=j_index(i);
1982          globalIndex = i_ind + j_ind * ni_glo;
1983          globalIndexDomain(i) = globalIndex;
1984        }
1985
1986        if (globalLocalIndexMap_.empty()) 
1987          for (i = 0; i < nbIndex; ++i)  globalLocalIndexMap_[globalIndexDomain(i)] = i;
1988         
1989
1990        size_t globalSizeIndex = 1, indexBegin, indexEnd;
1991        int range, clientSize = client->clientSize;
1992        std::vector<int> nGlobDomain(2);
1993        nGlobDomain[0] = this->ni_glo;
1994        nGlobDomain[1] = this->nj_glo;
1995        for (int i = 0; i < nGlobDomain.size(); ++i) globalSizeIndex *= nGlobDomain[i];
1996        indexBegin = 0;
1997        if (globalSizeIndex <= clientSize)
1998        {
1999          indexBegin = rank%globalSizeIndex;
2000          indexEnd = indexBegin;
2001        }
2002        else
2003        {
2004          for (int i = 0; i < clientSize; ++i)
2005          {
2006            range = globalSizeIndex / clientSize;
2007            if (i < (globalSizeIndex%clientSize)) ++range;
2008            if (i == client->clientRank) break;
2009            indexBegin += range;
2010          }
2011          indexEnd = indexBegin + range - 1;
2012        }
2013
2014        // Even if servers have no index, they must received something from client
2015        // We only use several client to send "empty" message to these servers
2016        CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2017        std::vector<int> serverZeroIndex;
2018        if (isUnstructed_) serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 0);
2019        else serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 1);
2020
2021        std::list<int> serverZeroIndexLeader;
2022        std::list<int> serverZeroIndexNotLeader;
2023        CContextClient::computeLeader(client->clientRank, client->clientSize, serverZeroIndex.size(), serverZeroIndexLeader, serverZeroIndexNotLeader);
2024        for (std::list<int>::iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
2025          *it = serverZeroIndex[*it];
2026
2027        CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(), client->intraComm);
2028        clientServerMap->computeServerIndexMapping(globalIndexDomain, nbServer);
2029        CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
2030
2031        CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(), ite = globalIndexDomainOnServer.end();
2032        indSrv_[nbServer].swap(globalIndexDomainOnServer);
2033        connectedServerRank_[nbServer].clear();
2034        for (it = indSrv_[nbServer].begin(); it != ite; ++it) connectedServerRank_[nbServer].push_back(it->first);
2035
2036        for (std::list<int>::const_iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
2037          connectedServerRank_[nbServer].push_back(*it);
2038
2039        // Even if a client has no index, it must connect to at least one server and
2040        // send an "empty" data to this server
2041        if (connectedServerRank_[nbServer].empty())
2042          connectedServerRank_[nbServer].push_back(client->clientRank % client->serverSize);
2043
2044        // Now check if all servers have data to receive. If not, master client will send empty data.
2045        // This ensures that all servers will participate in collective calls upon receiving even if they have no date to receive.
2046        std::vector<int> counts (clientSize);
2047        std::vector<int> displs (clientSize);
2048        displs[0] = 0;
2049        int localCount = connectedServerRank_[nbServer].size() ;
2050        MPI_Gather(&localCount, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, client->intraComm) ;
2051        for (int i = 0; i < clientSize-1; ++i) displs[i+1] = displs[i] + counts[i];
2052        std::vector<int> allConnectedServers(displs[clientSize-1]+counts[clientSize-1]);
2053        MPI_Gatherv(&(connectedServerRank_[nbServer])[0], localCount, MPI_INT, &allConnectedServers[0], &counts[0], &displs[0], MPI_INT, 0, client->intraComm);
2054
2055        if ((allConnectedServers.size() != nbServer) && (rank == 0))
2056        {
2057          std::vector<bool> isSrvConnected (nbServer, false);
2058          for (int i = 0; i < allConnectedServers.size(); ++i) isSrvConnected[allConnectedServers[i]] = true;
2059          for (int i = 0; i < nbServer; ++i) if (!isSrvConnected[i]) connectedServerRank_[nbServer].push_back(i);
2060        }
2061        nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
2062        delete clientServerMap;
2063      }
2064    }
2065  }
2066  CATCH_DUMP_ATTR
2067
2068   /*!
2069     Compute index to write data. We only write data on the zoomed region, therefore, there should
2070     be a map between the complete grid and the reduced grid where we write data.
2071     By using global index we can easily create this kind of mapping.
2072   */
2073   void CDomain::computeWrittenIndex()
2074   TRY
2075   { 
2076      if (computedWrittenIndex_) return;
2077      computedWrittenIndex_ = true;
2078
2079      CContext* context=CContext::getCurrent();     
2080
2081      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2082      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2083      nSize[0]        = ni;      nSize[1]  = nj;
2084      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2085      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2086      CDistributionServer srvDist(context->intraCommSize_, nBegin, nSize, nBeginGlobal, nGlob); 
2087      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2088
2089      size_t nbWritten = 0, indGlo;     
2090      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2091                                                          ite = globalLocalIndexMap_.end(), it;         
2092      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2093                                       itSrve = writtenGlobalIndex.end(), itSrv;
2094
2095      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
2096      nbWritten = 0;
2097      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2098      {
2099        indGlo = *itSrv;
2100        if (ite != globalLocalIndexMap_.find(indGlo))
2101        {
2102          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
2103        }
2104        else
2105        {
2106          localIndexToWriteOnServer(nbWritten) = -1;
2107        }
2108        ++nbWritten;
2109      }
2110   }
2111   CATCH_DUMP_ATTR
2112
2113  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
2114  TRY
2115  {
2116    int writtenCommSize;
2117    MPI_Comm_size(writtenComm, &writtenCommSize);
2118    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
2119      return;
2120
2121    if (isCompressible())
2122    {
2123      size_t nbWritten = 0, indGlo;
2124      CContext* context=CContext::getCurrent();     
2125
2126      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2127      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2128      nSize[0]        = ni;      nSize[1]  = nj;
2129      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2130      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2131      CDistributionServer srvDist(context->intraCommSize_, nBegin, nSize, nBeginGlobal, nGlob); 
2132      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2133
2134      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2135                                                          ite = globalLocalIndexMap_.end(), it;   
2136      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2137                                       itSrve = writtenGlobalIndex.end(), itSrv;
2138      std::unordered_map<size_t,size_t> localGlobalIndexMap;
2139      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2140      {
2141        indGlo = *itSrv;
2142        if (ite != globalLocalIndexMap_.find(indGlo))
2143        {
2144          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
2145          ++nbWritten;
2146        }                 
2147      }
2148
2149      nbWritten = 0;
2150      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2151      {
2152        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2153        {
2154          ++nbWritten;
2155        }
2156      }
2157
2158      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2159      nbWritten = 0;
2160      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2161      {
2162        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2163        {
2164          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2165          ++nbWritten;
2166        }
2167      }
2168
2169      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2170      bool distributed_glo, distributed=isDistributed() ;
2171      MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, writtenComm) ;
2172     
2173      if (distributed_glo)
2174      {
2175             
2176        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2177        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2178        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2179      }
2180      else
2181        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2182      }
2183  }
2184  CATCH_DUMP_ATTR
2185
2186  void CDomain::sendDomainToFileServer(CContextClient* client)
2187  {
2188    if (sendDomainToFileServer_done_.count(client)!=0) return ;
2189    else sendDomainToFileServer_done_.insert(client) ;
2190
2191    StdString domDefRoot("domain_definition");
2192    CDomainGroup* domPtr = CDomainGroup::get(domDefRoot);
2193    domPtr->sendCreateChild(this->getId(), client);
2194    this->sendAllAttributesToServer(client)  ;
2195    this->sendDistributionAttributes(client);   
2196    this->sendIndex(client);       
2197    this->sendLonLat(client);
2198    this->sendArea(client);   
2199    this->sendDataIndex(client);
2200
2201    // test new connector functionnality
2202    this->sendDomainDistribution(client) ;
2203  }
2204
2205  void CDomain::sendDomainToCouplerOut(CContextClient* client, const string& fieldId, int posInGrid)
2206  {
2207    if (sendDomainToFileServer_done_.count(client)!=0) return ;
2208    else sendDomainToFileServer_done_.insert(client) ;
2209   
2210    const string domainId = "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
2211   
2212    if (!domain_ref.isEmpty())
2213    {
2214      auto domain_ref_tmp=domain_ref.getValue() ;
2215      domain_ref.reset() ; // remove the reference, find an other way to do that more cleanly
2216      this->sendAllAttributesToServer(client, domainId)  ;
2217      domain_ref = domain_ref_tmp ;
2218    }
2219    else this->sendAllAttributesToServer(client, domainId)  ;
2220
2221    this->sendDistributionAttributes(client, domainId);   
2222    this->sendIndex(client, domainId);       
2223    this->sendLonLat(client, domainId);
2224    this->sendArea(client, domainId);   
2225    this->sendDataIndex(client, domainId);
2226  }
2227
2228  void CDomain::makeAliasForCoupling(const string& fieldId, int posInGrid)
2229  {
2230    const string domainId = "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
2231    this->createAlias(domainId) ;
2232  }
2233
2234
2235  void CDomain::sendDomainDistribution(CContextClient* client, const string& domainId)
2236  TRY
2237  {
2238    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2239    CContext* context = CContext::getCurrent();
2240    int nbServer = client->serverSize;
2241    std::vector<int> nGlobDomain(2);
2242    nGlobDomain[0] = this->ni_glo;
2243    nGlobDomain[1] = this->nj_glo;
2244
2245    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2246    int distributedPosition ;
2247    if (isUnstructed_) distributedPosition = 0 ;
2248    else distributedPosition = 1 ;
2249
2250    serverDescription.computeServerDistribution(false, distributedPosition);
2251   
2252    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2253    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2254 
2255    vector<unordered_map<size_t,vector<int>>> indexServerOnElement ;
2256    CArray<int,1> axisDomainOrder(1) ; axisDomainOrder(0)=2 ;
2257    auto zeroIndex=serverDescription.computeServerGlobalByElement(indexServerOnElement, context->getIntraCommRank(), context->getIntraCommSize(),
2258                                                                  axisDomainOrder,distributedPosition) ;
2259    // distribution is very bad => to redo
2260    // convert indexServerOnElement => map<int,CArray<size_t,1>> - need to be changed later
2261    map<int, vector<size_t>> vectGlobalIndex ;
2262    for(auto& indexRanks : indexServerOnElement[0])
2263    {
2264      size_t index=indexRanks.first ;
2265      auto& ranks=indexRanks.second ;
2266      for(int rank : ranks) vectGlobalIndex[rank].push_back(index) ;
2267    }
2268    map<int, CArray<size_t,1>> globalIndex ;
2269    for(auto& vect : vectGlobalIndex ) globalIndex.emplace(vect.first, CArray<size_t,1>(vect.second.data(), shape(vect.second.size()))) ; 
2270
2271    CDistributedElement remoteElement(ni_glo*nj_glo, globalIndex) ;
2272    remoteElement.addFullView() ;
2273
2274    CRemoteConnector remoteConnector(localElement_->getView(CElementView::FULL), remoteElement.getView(CElementView::FULL),context->getIntraComm()) ;
2275    remoteConnector.computeConnector() ;
2276    CDistributedElement scatteredElement(remoteElement.getGlobalSize(), remoteConnector.getDistributedGlobalIndex()) ;
2277    scatteredElement.addFullView() ;
2278    CScattererConnector scatterConnector(localElement_->getView(CElementView::FULL), scatteredElement.getView(CElementView::FULL), context->getIntraComm()) ;
2279    scatterConnector.computeConnector() ;
2280    CGridScattererConnector gridScatter({&scatterConnector}) ;
2281
2282    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2283    CMessage message0 ;
2284    message0<<serverDomainId<<0 ; 
2285    remoteElement.sendToServer(client,event0,message0) ;
2286
2287    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2288    CMessage message1 ;
2289    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2290    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2291   
2292    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2293    CMessage message2 ;
2294    message2<<serverDomainId<<2 ; 
2295//    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2296    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2297
2298/*
2299    localElement_->getView(CElementView::FULL)->sendRemoteElement(remoteConnector, client, event1, message1) ;
2300    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2301    CMessage message2 ;
2302    message2<<serverDomainId<<2 ;
2303    remoteConnector.transferToServer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2304*/
2305
2306  }
2307  CATCH
2308 
2309
2310  void CDomain::recvDomainDistribution(CEventServer& event)
2311  TRY
2312  {
2313    string domainId;
2314    int phasis ;
2315    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
2316    get(domainId)->receivedDomainDistribution(event, phasis);
2317  }
2318  CATCH
2319
2320  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2321  TRY
2322  {
2323    CContext* context = CContext::getCurrent();
2324    if (phasis==0)
2325    {
2326      localElement_ = new  CLocalElement(context->getIntraCommRank(),event) ;
2327      localElement_->addFullView() ;
2328    }
2329    else if (phasis==1)
2330    {
2331      CContext* context = CContext::getCurrent();
2332      CDistributedElement* elementFrom = new  CDistributedElement(event) ;
2333      elementFrom->addFullView() ;
2334      gathererConnector_ = new CGathererConnector(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2335      gathererConnector_->computeConnector() ; 
2336    }
2337    else if (phasis==2)
2338    {
2339      CArray<size_t,1> globalIndex ;
2340      //gathererConnector_->transfer(event,globalIndex) ;
2341      CGridGathererConnector gridGathererConnector({gathererConnector_}) ;
2342      gridGathererConnector.transfer(event, globalIndex) ;
2343    }
2344    else if (phasis==3)
2345    {
2346
2347    }
2348  }
2349  CATCH
2350
2351 
2352
2353  /*!
2354    Send all attributes from client to connected clients
2355    The attributes will be rebuilt on receiving side
2356  */
2357  // ym obsolete to be removed
2358  void CDomain::sendAttributes()
2359  TRY
2360  {
2361    //sendDistributionAttributes();
2362    //sendIndex();       
2363    //sendLonLat();
2364    //sendArea();   
2365    //sendDataIndex();
2366  }
2367  CATCH
2368  /*!
2369    Send global index from client to connected client(s)
2370  */
2371  void CDomain::sendIndex(CContextClient* client, const string& domainId)
2372  TRY
2373  {
2374    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2375   
2376    int ns, n, i, j, ind, nv, idx;
2377    int serverSize = client->serverSize;
2378    CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2379
2380    list<CMessage> list_msgsIndex;
2381    list<CArray<int,1> > list_indGlob;
2382
2383    std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2384    iteIndex = indSrv_[serverSize].end();
2385    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2386    {
2387      int nbIndGlob = 0;
2388      int rank = connectedServerRank_[serverSize][k];
2389      itIndex = indSrv_[serverSize].find(rank);
2390      if (iteIndex != itIndex)
2391        nbIndGlob = itIndex->second.size();
2392
2393      list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2394
2395      CArray<int,1>& indGlob = list_indGlob.back();
2396      for (n = 0; n < nbIndGlob; ++n) indGlob(n) = static_cast<int>(itIndex->second[n]);
2397
2398      list_msgsIndex.push_back(CMessage());
2399      list_msgsIndex.back() << serverDomainId << (int)type; // enum ne fonctionne pour les message => ToFix
2400      list_msgsIndex.back() << isCurvilinear;
2401      list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2402       
2403      eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2404    }
2405    client->sendEvent(eventIndex);
2406  }
2407  CATCH_DUMP_ATTR
2408
2409  /*!
2410    Send distribution from client to other clients
2411    Because a client in a level knows correctly the grid distribution of client on the next level
2412    it calculates this distribution then sends it to the corresponding clients on the next level
2413  */
2414  void CDomain::sendDistributionAttributes(CContextClient* client, const string& domainId)
2415  TRY
2416  {
2417    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2418
2419    int nbServer = client->serverSize;
2420    std::vector<int> nGlobDomain(2);
2421    nGlobDomain[0] = this->ni_glo;
2422    nGlobDomain[1] = this->nj_glo;
2423
2424    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2425    if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2426    else serverDescription.computeServerDistribution(false, 1);
2427
2428    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2429    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2430
2431    CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2432    if (client->isServerLeader())
2433    {
2434      std::list<CMessage> msgs;
2435
2436      const std::list<int>& ranks = client->getRanksServerLeader();
2437      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2438      {
2439        // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2440        const int ibegin_srv = serverIndexBegin[*itRank][0];
2441        const int jbegin_srv = serverIndexBegin[*itRank][1];
2442        const int ni_srv = serverDimensionSizes[*itRank][0];
2443        const int nj_srv = serverDimensionSizes[*itRank][1];
2444
2445        msgs.push_back(CMessage());
2446        CMessage& msg = msgs.back();
2447        msg << serverDomainId ;
2448        msg << isUnstructed_;
2449        msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2450        msg << ni_glo.getValue() << nj_glo.getValue();
2451        msg << isCompressible_;
2452
2453        event.push(*itRank,1,msg);
2454      }
2455      client->sendEvent(event);
2456    }
2457    else client->sendEvent(event);
2458  }
2459  CATCH_DUMP_ATTR
2460
2461  /*!
2462    Send area from client to connected client(s)
2463  */
2464  void CDomain::sendArea(CContextClient* client, const string& domainId)
2465  TRY
2466  {
2467    if (!hasArea) return;
2468    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2469   
2470    int ns, n, i, j, ind, nv, idx;
2471    int serverSize = client->serverSize;
2472
2473    // send area for each connected server
2474    CEventClient eventArea(getType(), EVENT_ID_AREA);
2475
2476    list<CMessage> list_msgsArea;
2477    list<CArray<double,1> > list_area;
2478
2479    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2480    iteMap = indSrv_[serverSize].end();
2481    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2482    {
2483      int nbData = 0;
2484      int rank = connectedServerRank_[serverSize][k];
2485      it = indSrv_[serverSize].find(rank);
2486      if (iteMap != it)
2487        nbData = it->second.size();
2488      list_area.push_back(CArray<double,1>(nbData));
2489
2490      const std::vector<size_t>& temp = it->second;
2491      for (n = 0; n < nbData; ++n)
2492      {
2493        idx = static_cast<int>(it->second[n]);
2494        list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2495      }
2496
2497      list_msgsArea.push_back(CMessage());
2498      list_msgsArea.back() <<serverDomainId << hasArea;
2499      list_msgsArea.back() << list_area.back();
2500      eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2501    }
2502    client->sendEvent(eventArea);
2503  }
2504  CATCH_DUMP_ATTR
2505
2506  /*!
2507    Send longitude and latitude from client to servers
2508    Each client send long and lat information to corresponding connected clients(s).
2509    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2510  */
2511  void CDomain::sendLonLat(CContextClient* client, const string& domainId)
2512  TRY
2513  {
2514    if (!hasLonLat) return;
2515    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2516   
2517    int ns, n, i, j, ind, nv, idx;
2518    int serverSize = client->serverSize;
2519
2520    // send lon lat for each connected server
2521    CEventClient eventLon(getType(), EVENT_ID_LON);
2522    CEventClient eventLat(getType(), EVENT_ID_LAT);
2523
2524    list<CMessage> list_msgsLon, list_msgsLat;
2525    list<CArray<double,1> > list_lon, list_lat;
2526    list<CArray<double,2> > list_boundslon, list_boundslat;
2527
2528    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2529    iteMap = indSrv_[serverSize].end();
2530    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2531    {
2532      int nbData = 0;
2533      int rank = connectedServerRank_[serverSize][k];
2534      it = indSrv_[serverSize].find(rank);
2535      if (iteMap != it)
2536        nbData = it->second.size();
2537
2538      list_lon.push_back(CArray<double,1>(nbData));
2539      list_lat.push_back(CArray<double,1>(nbData));
2540
2541      if (hasBounds)
2542      {
2543        list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2544        list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2545      }
2546
2547      CArray<double,1>& lon = list_lon.back();
2548      CArray<double,1>& lat = list_lat.back();
2549      const std::vector<size_t>& temp = it->second;
2550      for (n = 0; n < nbData; ++n)
2551      {
2552        idx = static_cast<int>(it->second[n]);
2553        int localInd = globalLocalIndexMap_[idx];
2554        lon(n) = lonvalue(localInd);
2555        lat(n) = latvalue(localInd);
2556
2557        if (hasBounds)
2558        {
2559          CArray<double,2>& boundslon = list_boundslon.back();
2560          CArray<double,2>& boundslat = list_boundslat.back();
2561
2562          for (nv = 0; nv < nvertex; ++nv)
2563          {
2564            boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2565            boundslat(nv, n) = bounds_latvalue(nv, localInd);
2566          }
2567        }
2568      }
2569
2570      list_msgsLon.push_back(CMessage());
2571      list_msgsLat.push_back(CMessage());
2572
2573      list_msgsLon.back() << serverDomainId << hasLonLat;
2574      if (hasLonLat) 
2575        list_msgsLon.back() << list_lon.back();
2576      list_msgsLon.back()  << hasBounds;
2577      if (hasBounds)
2578      {
2579        list_msgsLon.back() << list_boundslon.back();
2580      }
2581
2582      list_msgsLat.back() << serverDomainId << hasLonLat;
2583      if (hasLonLat)
2584        list_msgsLat.back() << list_lat.back();
2585      list_msgsLat.back() << hasBounds;
2586      if (hasBounds)
2587      {         
2588        list_msgsLat.back() << list_boundslat.back();
2589      }
2590
2591      eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2592      eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2593    }
2594    client->sendEvent(eventLon);
2595    client->sendEvent(eventLat);
2596  }
2597  CATCH_DUMP_ATTR
2598
2599  /*!
2600    Send data index to corresponding connected clients.
2601    Data index can be compressed however, we always send decompressed data index
2602    and they will be compressed on receiving.
2603    The compressed index are represented with 1 and others are represented with -1
2604  */
2605  void CDomain::sendDataIndex(CContextClient* client, const string& domainId)
2606  TRY
2607  {
2608    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2609   
2610    int ns, n, i, j, ind, nv, idx;
2611    int serverSize = client->serverSize;
2612
2613    // send area for each connected server
2614    CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2615
2616    list<CMessage> list_msgsDataIndex;
2617    list<CArray<int,1> > list_data_i_index, list_data_j_index;
2618
2619    int nbIndex = i_index.numElements();
2620    int niByIndex = max(i_index) - min(i_index) + 1;
2621    int njByIndex = max(j_index) - min(j_index) + 1; 
2622    int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2623    int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2624
2625   
2626    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2627    dataIIndex = -1; 
2628    dataJIndex = -1;
2629    ind = 0;
2630
2631    for (idx = 0; idx < data_i_index.numElements(); ++idx)
2632    {
2633      int dataIidx = data_i_index(idx) + data_ibegin;
2634      int dataJidx = data_j_index(idx) + data_jbegin;
2635      if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2636          (0 <= dataJidx) && (dataJidx < dataJindexBound))
2637      {
2638        dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2639        dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2640      }
2641    }
2642
2643    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2644    iteMap = indSrv_[serverSize].end();
2645    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2646    {
2647      int nbData = 0;
2648      int rank = connectedServerRank_[serverSize][k];
2649      it = indSrv_[serverSize].find(rank);
2650      if (iteMap != it)
2651        nbData = it->second.size();
2652      list_data_i_index.push_back(CArray<int,1>(nbData));
2653      list_data_j_index.push_back(CArray<int,1>(nbData));
2654
2655      const std::vector<size_t>& temp = it->second;
2656      for (n = 0; n < nbData; ++n)
2657      {
2658        idx = static_cast<int>(it->second[n]);
2659        i = globalLocalIndexMap_[idx];
2660        list_data_i_index.back()(n) = dataIIndex(i);
2661        list_data_j_index.back()(n) = dataJIndex(i);
2662      }
2663
2664      list_msgsDataIndex.push_back(CMessage());
2665      list_msgsDataIndex.back() << serverDomainId ;
2666      list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2667      eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2668    }
2669    client->sendEvent(eventDataIndex);
2670  }
2671  CATCH
2672 
2673  bool CDomain::dispatchEvent(CEventServer& event)
2674  TRY
2675  {
2676    if (SuperClass::dispatchEvent(event)) return true;
2677    else
2678    {
2679      switch(event.type)
2680      {
2681        case EVENT_ID_SERVER_ATTRIBUT:
2682          recvDistributionAttributes(event);
2683          return true;
2684          break;
2685        case EVENT_ID_INDEX:
2686          recvIndex(event);
2687          return true;
2688          break;
2689        case EVENT_ID_LON:
2690          recvLon(event);
2691          return true;
2692          break;
2693        case EVENT_ID_LAT:
2694          recvLat(event);
2695          return true;
2696          break;
2697        case EVENT_ID_AREA:
2698          recvArea(event);
2699          return true;
2700          break; 
2701        case EVENT_ID_DATA_INDEX:
2702          recvDataIndex(event);
2703          return true;
2704          break;
2705        case EVENT_ID_DOMAIN_DISTRIBUTION:
2706          recvDomainDistribution(event);
2707          return true;
2708          break;
2709        default:
2710          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2711                << "Unknown Event");
2712          return false;
2713       }
2714    }
2715  }
2716  CATCH
2717
2718  /*!
2719    Receive index event from clients(s)
2720    \param[in] event event contain info about rank and associated index
2721  */
2722  void CDomain::recvIndex(CEventServer& event)
2723  TRY
2724  {
2725    string domainId;
2726    std::map<int, CBufferIn*> rankBuffers;
2727
2728    list<CEventServer::SSubEvent>::iterator it;
2729    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2730    {     
2731      CBufferIn* buffer = it->buffer;
2732      *buffer >> domainId;
2733      rankBuffers[it->rank] = buffer;       
2734    }
2735    get(domainId)->recvIndex(rankBuffers);
2736  }
2737  CATCH
2738
2739  /*!
2740    Receive index information from client(s). We use the global index for mapping index between
2741    sending clients and receiving clients.
2742    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2743  */
2744  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2745  TRY
2746  {
2747    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2748    recvClientRanks_.resize(nbReceived);       
2749
2750    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2751    ind = 0;
2752    for (ind = 0; it != ite; ++it, ++ind)
2753    {       
2754       recvClientRanks_[ind] = it->first;
2755       CBufferIn& buffer = *(it->second);
2756       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2757       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2758    }
2759    int nbIndGlob = 0;
2760    for (i = 0; i < nbReceived; ++i)
2761    {
2762      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2763    }
2764   
2765    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2766    i_index.resize(nbIndGlob);
2767    j_index.resize(nbIndGlob);   
2768    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2769
2770    nbIndGlob = 0;
2771    for (i = 0; i < nbReceived; ++i)
2772    {
2773      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2774      for (ind = 0; ind < tmp.numElements(); ++ind)
2775      {
2776         index = tmp(ind);
2777         if (0 == globalLocalIndexMap_.count(index))
2778         {
2779           iIndex = (index%ni_glo)-ibegin;
2780           iIndex = (iIndex < 0) ? 0 : iIndex;
2781           jIndex = (index/ni_glo)-jbegin;
2782           jIndex = (jIndex < 0) ? 0 : jIndex;
2783           nbIndLoc = iIndex + ni * jIndex;
2784           i_index(nbIndGlob) = index % ni_glo;
2785           j_index(nbIndGlob) = index / ni_glo;
2786           globalLocalIndexMap_[index] = nbIndGlob;
2787           ++nbIndGlob;
2788         } 
2789      } 
2790    } 
2791
2792    if (nbIndGlob==0)
2793    {
2794      i_index.resize(nbIndGlob);
2795      j_index.resize(nbIndGlob);
2796    }
2797    else
2798    {
2799      i_index.resizeAndPreserve(nbIndGlob);
2800      j_index.resizeAndPreserve(nbIndGlob);
2801    }
2802
2803    domainMask.resize(0); // Mask is not defined anymore on servers
2804  }
2805  CATCH
2806
2807  /*!
2808    Receive attributes event from clients(s)
2809    \param[in] event event contain info about rank and associated attributes
2810  */
2811  void CDomain::recvDistributionAttributes(CEventServer& event)
2812  TRY
2813  {
2814    CBufferIn* buffer=event.subEvents.begin()->buffer;
2815    string domainId ;
2816    *buffer>>domainId ;
2817    get(domainId)->recvDistributionAttributes(*buffer);
2818  }
2819  CATCH
2820
2821  /*!
2822    Receive attributes from client(s)
2823    \param[in] rank rank of client source
2824    \param[in] buffer message containing attributes info
2825  */
2826  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2827  TRY
2828  {
2829    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2830    int ni_glo_tmp, nj_glo_tmp;
2831    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2832           >> ni_glo_tmp >> nj_glo_tmp
2833           >> isCompressible_;
2834
2835    ni.setValue(ni_tmp);
2836    ibegin.setValue(ibegin_tmp);
2837    nj.setValue(nj_tmp);
2838    jbegin.setValue(jbegin_tmp);
2839    ni_glo.setValue(ni_glo_tmp);
2840    nj_glo.setValue(nj_glo_tmp);
2841
2842  }
2843 CATCH_DUMP_ATTR
2844  /*!
2845    Receive longitude event from clients(s)
2846    \param[in] event event contain info about rank and associated longitude
2847  */
2848  void CDomain::recvLon(CEventServer& event)
2849  TRY
2850  {
2851    string domainId;
2852    std::map<int, CBufferIn*> rankBuffers;
2853
2854    list<CEventServer::SSubEvent>::iterator it;
2855    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2856    {     
2857      CBufferIn* buffer = it->buffer;
2858      *buffer >> domainId;
2859      rankBuffers[it->rank] = buffer;       
2860    }
2861    get(domainId)->recvLon(rankBuffers);
2862  }
2863  CATCH
2864
2865  /*!
2866    Receive longitude information from client(s)
2867    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2868  */
2869  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2870  TRY
2871  {
2872    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2873    if (nbReceived != recvClientRanks_.size())
2874      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2875           << "The number of sending clients is not correct."
2876           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2877   
2878    int nbLonInd = 0;
2879    vector<CArray<double,1> > recvLonValue(nbReceived);
2880    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2881    for (i = 0; i < recvClientRanks_.size(); ++i)
2882    {
2883      int rank = recvClientRanks_[i];
2884      CBufferIn& buffer = *(rankBuffers[rank]);
2885      buffer >> hasLonLat;
2886      if (hasLonLat)
2887        buffer >> recvLonValue[i];
2888      buffer >> hasBounds;
2889      if (hasBounds)
2890        buffer >> recvBoundsLonValue[i];
2891    }
2892
2893    if (hasLonLat)
2894    {
2895      for (i = 0; i < nbReceived; ++i)
2896      {
2897        nbLonInd += recvLonValue[i].numElements();
2898      }
2899   
2900      if (nbLonInd != globalLocalIndexMap_.size())
2901        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2902                 << "something must be wrong with longitude index "<< std::endl;
2903
2904      nbLonInd = globalLocalIndexMap_.size();
2905      lonvalue.resize(nbLonInd);
2906      if (hasBounds)
2907      {
2908        bounds_lonvalue.resize(nvertex,nbLonInd);
2909        bounds_lonvalue = 0.;
2910      }
2911
2912      nbLonInd = 0;
2913      for (i = 0; i < nbReceived; ++i)
2914      {
2915        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2916        CArray<double,1>& tmp = recvLonValue[i];
2917        for (ind = 0; ind < tmp.numElements(); ++ind)
2918        {
2919          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2920          lonvalue(lInd) = tmp(ind); 
2921           if (hasBounds)
2922           {         
2923            for (int nv = 0; nv < nvertex; ++nv)
2924              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2925           }                 
2926        }
2927      }       
2928    }
2929
2930   // setup attribute depending the type of domain
2931    if (hasLonLat)
2932    {
2933      nbLonInd = globalLocalIndexMap_.size();
2934      if (ni*nj != nbLonInd)
2935        ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2936             << "The number of index received is not coherent with the given resolution"
2937             << "nbLonInd=" << nbLonInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
2938
2939      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
2940      {
2941        lonvalue_2d.resize(ni,nj);
2942          for(int ij=0, j=0 ; j<nj ; j++)
2943            for(int i=0 ; i<ni; i++, ij++) lonvalue_2d(i,j) = lonvalue(ij) ;
2944       
2945        if (hasBounds)
2946        {
2947          bounds_lon_2d.resize(nvertex, ni, nj) ;
2948          for(int ij=0, j=0 ; j<nj ; j++)
2949            for(int i=0 ; i<ni; i++, ij++) 
2950              for(int nv=0; nv<nvertex; nv++) bounds_lon_2d(nv,i,j) = bounds_lonvalue(nv,ij) ;
2951        }
2952      }
2953      else if (type == type_attr::unstructured || type == type_attr::gaussian)
2954      {
2955        lonvalue_1d.resize(nbLonInd);
2956        lonvalue_1d = lonvalue ;
2957        if (hasBounds)
2958        {
2959          bounds_lon_1d.resize(nvertex, nbLonInd) ;
2960          bounds_lon_1d = bounds_lonvalue ;
2961        }
2962      }
2963    }
2964  }
2965  CATCH_DUMP_ATTR
2966
2967  /*!
2968    Receive latitude event from clients(s)
2969    \param[in] event event contain info about rank and associated latitude
2970  */
2971  void CDomain::recvLat(CEventServer& event)
2972  TRY
2973  {
2974    string domainId;
2975    std::map<int, CBufferIn*> rankBuffers;
2976
2977    list<CEventServer::SSubEvent>::iterator it;
2978    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2979    {     
2980      CBufferIn* buffer = it->buffer;
2981      *buffer >> domainId;
2982      rankBuffers[it->rank] = buffer;   
2983    }
2984    get(domainId)->recvLat(rankBuffers);
2985  }
2986  CATCH
2987
2988  /*!
2989    Receive latitude information from client(s)
2990    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2991  */
2992  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2993  TRY
2994  {
2995    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2996    if (nbReceived != recvClientRanks_.size())
2997      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2998           << "The number of sending clients is not correct."
2999           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3000
3001    vector<CArray<double,1> > recvLatValue(nbReceived);
3002    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
3003    for (i = 0; i < recvClientRanks_.size(); ++i)
3004    {
3005      int rank = recvClientRanks_[i];
3006      CBufferIn& buffer = *(rankBuffers[rank]);
3007      buffer >> hasLonLat;
3008      if (hasLonLat)
3009        buffer >> recvLatValue[i];
3010      buffer >> hasBounds;
3011      if (hasBounds)
3012        buffer >> recvBoundsLatValue[i];
3013    }
3014
3015    int nbLatInd = 0;
3016    if (hasLonLat)
3017    {
3018      for (i = 0; i < nbReceived; ++i)
3019      {
3020        nbLatInd += recvLatValue[i].numElements();
3021      }
3022   
3023      if (nbLatInd != globalLocalIndexMap_.size())
3024        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3025                << "something must be wrong with latitude index "<< std::endl;
3026
3027      nbLatInd = globalLocalIndexMap_.size();
3028      latvalue.resize(nbLatInd);
3029      if (hasBounds)
3030      {
3031        bounds_latvalue.resize(nvertex,nbLatInd);
3032        bounds_latvalue = 0. ;
3033      }
3034
3035      nbLatInd = 0;
3036      for (i = 0; i < nbReceived; ++i)
3037      {
3038        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3039        CArray<double,1>& tmp = recvLatValue[i];
3040        for (ind = 0; ind < tmp.numElements(); ++ind)
3041        {
3042          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3043          latvalue(lInd) = tmp(ind);   
3044           if (hasBounds)
3045           {
3046            CArray<double,2>& boundslat = recvBoundsLatValue[i];
3047            for (int nv = 0; nv < nvertex; ++nv)
3048              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
3049           }   
3050          ++nbLatInd;
3051        }
3052      }       
3053    }
3054      // setup attribute depending the type of domain
3055    if (hasLonLat)
3056    {
3057      nbLatInd = globalLocalIndexMap_.size();
3058     
3059      if (ni*nj != nbLatInd)
3060        ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
3061             << "The number of index received is not coherent with the given resolution"
3062            << "nbLonInd=" << nbLatInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
3063
3064      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
3065      {
3066        {
3067          latvalue_2d.resize(ni,nj);
3068          for(int ij=0, j=0 ; j<nj ; j++)
3069            for(int i=0 ; i<ni; i++, ij++) latvalue_2d(i,j) = latvalue(ij) ;
3070       
3071          if (hasBounds)
3072          {
3073            bounds_lat_2d.resize(nvertex, ni, nj) ;
3074            for(int ij=0, j=0 ; j<nj ; j++)
3075              for(int i=0 ; i<ni; i++, ij++) 
3076                for(int nv=0; nv<nvertex; nv++) bounds_lat_2d(nv,i,j) = bounds_latvalue(nv,ij) ;
3077          }
3078        }
3079      }
3080      else if (type == type_attr::unstructured || type == type_attr::gaussian)
3081      {
3082        if (hasLonLat)
3083        {
3084          latvalue_1d.resize(nbLatInd);
3085          latvalue_1d = latvalue ;
3086          if (hasBounds)
3087          {
3088            bounds_lat_1d.resize(nvertex, nbLatInd) ;
3089            bounds_lat_1d = bounds_latvalue ;
3090          }
3091        }
3092      }
3093    } 
3094  }
3095  CATCH_DUMP_ATTR
3096
3097  /*!
3098    Receive area event from clients(s)
3099    \param[in] event event contain info about rank and associated area
3100  */
3101  void CDomain::recvArea(CEventServer& event)
3102  TRY
3103  {
3104    string domainId;
3105    std::map<int, CBufferIn*> rankBuffers;
3106
3107    list<CEventServer::SSubEvent>::iterator it;
3108    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3109    {     
3110      CBufferIn* buffer = it->buffer;
3111      *buffer >> domainId;
3112      rankBuffers[it->rank] = buffer;     
3113    }
3114    get(domainId)->recvArea(rankBuffers);
3115  }
3116  CATCH
3117
3118  /*!
3119    Receive area information from client(s)
3120    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3121  */
3122  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
3123  TRY
3124  {
3125    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
3126    if (nbReceived != recvClientRanks_.size())
3127      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
3128           << "The number of sending clients is not correct."
3129           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3130
3131    vector<CArray<double,1> > recvAreaValue(nbReceived);     
3132    for (i = 0; i < recvClientRanks_.size(); ++i)
3133    {
3134      int rank = recvClientRanks_[i];
3135      CBufferIn& buffer = *(rankBuffers[rank]);     
3136      buffer >> hasArea;
3137      if (hasArea)
3138        buffer >> recvAreaValue[i];
3139    }
3140
3141    if (hasArea)
3142    {
3143      int nbAreaInd = 0;
3144      for (i = 0; i < nbReceived; ++i)
3145      {     
3146        nbAreaInd += recvAreaValue[i].numElements();
3147      }
3148
3149      if (nbAreaInd != globalLocalIndexMap_.size())
3150        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3151                 << "something must be wrong with area index "<< std::endl;
3152
3153      nbAreaInd = globalLocalIndexMap_.size();
3154      areavalue.resize(nbAreaInd);
3155      nbAreaInd = 0;     
3156      for (i = 0; i < nbReceived; ++i)
3157      {
3158        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3159        CArray<double,1>& tmp = recvAreaValue[i];
3160        for (ind = 0; ind < tmp.numElements(); ++ind)
3161        {
3162          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3163          areavalue(lInd) = tmp(ind);         
3164        }
3165      }
3166     
3167      // fill the area attribute
3168      area.resize(ni,nj);
3169      for (int j = 0; j < nj; ++j)
3170      {
3171        for (int i = 0; i < ni; ++i)
3172        {
3173          int k = j * ni + i;
3174          area(i,j) = areavalue(k) ;
3175        }
3176      }
3177    }
3178  }
3179  CATCH_DUMP_ATTR
3180
3181  /*!
3182    Compare two domain objects.
3183    They are equal if only if they have identical attributes as well as their values.
3184    Moreover, they must have the same transformations.
3185  \param [in] domain Compared domain
3186  \return result of the comparison
3187  */
3188  bool CDomain::isEqual(CDomain* obj)
3189  TRY
3190  {
3191    vector<StdString> excludedAttr;
3192    excludedAttr.push_back("domain_ref");
3193    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
3194    if (!objEqual) return objEqual;
3195
3196    TransMapTypes thisTrans = this->getAllTransformations();
3197    TransMapTypes objTrans  = obj->getAllTransformations();
3198
3199    TransMapTypes::const_iterator it, itb, ite;
3200    std::vector<ETranformationType> thisTransType, objTransType;
3201    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
3202      thisTransType.push_back(it->first);
3203    for (it = objTrans.begin(); it != objTrans.end(); ++it)
3204      objTransType.push_back(it->first);
3205
3206    if (thisTransType.size() != objTransType.size()) return false;
3207    for (int idx = 0; idx < thisTransType.size(); ++idx)
3208      objEqual &= (thisTransType[idx] == objTransType[idx]);
3209
3210    return objEqual;
3211  }
3212  CATCH_DUMP_ATTR
3213
3214  /*!
3215    Receive data index event from clients(s)
3216    \param[in] event event contain info about rank and associated index
3217  */
3218  void CDomain::recvDataIndex(CEventServer& event)
3219  TRY
3220  {
3221    string domainId;
3222    std::map<int, CBufferIn*> rankBuffers;
3223
3224    list<CEventServer::SSubEvent>::iterator it;
3225    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3226    {     
3227      CBufferIn* buffer = it->buffer;
3228      *buffer >> domainId;
3229      rankBuffers[it->rank] = buffer;       
3230    }
3231    get(domainId)->recvDataIndex(rankBuffers);
3232  }
3233  CATCH
3234
3235  /*!
3236    Receive data index information from client(s)
3237    A client receives data index from different clients to rebuild its own data index.
3238    Because we use global index + mask info to calculate the sending data to client(s),
3239    this data index must be updated with mask info (maybe it will change in the future)
3240    Because the data index is local, to rebuild data index of received client, we should use global index along with.
3241
3242    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3243  */
3244  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
3245  TRY
3246  {
3247    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
3248    if (nbReceived != recvClientRanks_.size())
3249      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3250           << "The number of sending clients is not correct."
3251           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3252
3253    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3254    for (i = 0; i < recvClientRanks_.size(); ++i)
3255    {
3256      int rank = recvClientRanks_[i];
3257      CBufferIn& buffer = *(rankBuffers[rank]);
3258      buffer >> recvDataIIndex[i];
3259      buffer >> recvDataJIndex[i];
3260    }
3261   
3262    int nbIndex = i_index.numElements();
3263    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3264    dataIIndex = -1; dataJIndex = -1;
3265     
3266    nbIndex = 0;
3267    for (i = 0; i < nbReceived; ++i)
3268    {     
3269      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3270      CArray<int,1>& tmpI = recvDataIIndex[i];   
3271      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3272      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3273          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3274             << "The number of global received index is not coherent with the number of received data index."
3275             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3276
3277      for (ind = 0; ind < tmpI.numElements(); ++ind)
3278      {
3279         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3280         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3281         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3282      } 
3283    }
3284
3285    int nbCompressedData = 0; 
3286    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3287    {
3288       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3289       if ((0 <= indexI) && (0 <= indexJ))
3290         ++nbCompressedData;
3291    }       
3292 
3293    data_i_index.resize(nbCompressedData);
3294    data_j_index.resize(nbCompressedData);
3295
3296    nbCompressedData = 0; 
3297    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3298    {
3299       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3300       if ((0 <= indexI) && (0 <= indexJ))
3301       {
3302          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3303          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3304         ++nbCompressedData;
3305       }
3306    }
3307
3308    // Reset data_ibegin, data_jbegin
3309    data_ibegin.setValue(0);
3310    data_jbegin.setValue(0);
3311  }
3312  CATCH_DUMP_ATTR
3313
3314  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3315  TRY
3316  {
3317    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3318    return transformationMap_.back().second;
3319  }
3320  CATCH_DUMP_ATTR
3321
3322  /*!
3323    Check whether a domain has transformation
3324    \return true if domain has transformation
3325  */
3326  bool CDomain::hasTransformation()
3327  TRY
3328  {
3329    return (!transformationMap_.empty());
3330  }
3331  CATCH_DUMP_ATTR
3332
3333  /*!
3334    Set transformation for current domain. It's the method to move transformation in hierarchy
3335    \param [in] domTrans transformation on domain
3336  */
3337  void CDomain::setTransformations(const TransMapTypes& domTrans)
3338  TRY
3339  {
3340    transformationMap_ = domTrans;
3341  }
3342  CATCH_DUMP_ATTR
3343
3344  /*!
3345    Get all transformation current domain has
3346    \return all transformation
3347  */
3348  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3349  TRY
3350  {
3351    return transformationMap_;
3352  }
3353  CATCH_DUMP_ATTR
3354
3355  void CDomain::duplicateTransformation(CDomain* src)
3356  TRY
3357  {
3358    if (src->hasTransformation())
3359    {
3360      this->setTransformations(src->getAllTransformations());
3361    }
3362  }
3363  CATCH_DUMP_ATTR
3364   
3365  /*!
3366   * Go through the hierarchy to find the domain from which the transformations must be inherited
3367   */
3368  void CDomain::solveInheritanceTransformation()
3369  TRY
3370  {
3371    if (hasTransformation() || !hasDirectDomainReference())
3372      return;
3373
3374    CDomain* domain = this;
3375    std::vector<CDomain*> refDomains;
3376    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3377    {
3378      refDomains.push_back(domain);
3379      domain = domain->getDirectDomainReference();
3380    }
3381
3382    if (domain->hasTransformation())
3383      for (size_t i = 0; i < refDomains.size(); ++i)
3384        refDomains[i]->setTransformations(domain->getAllTransformations());
3385  }
3386  CATCH_DUMP_ATTR
3387
3388  void CDomain::setContextClient(CContextClient* contextClient)
3389  TRY
3390  {
3391    if (clientsSet.find(contextClient)==clientsSet.end())
3392    {
3393      clients.push_back(contextClient) ;
3394      clientsSet.insert(contextClient);
3395    }
3396  }
3397  CATCH_DUMP_ATTR
3398
3399  /*!
3400    Parse children nodes of a domain in xml file.
3401    Whenver there is a new transformation, its type and name should be added into this function
3402    \param node child node to process
3403  */
3404  void CDomain::parse(xml::CXMLNode & node)
3405  TRY
3406  {
3407    SuperClass::parse(node);
3408
3409    if (node.goToChildElement())
3410    {
3411      StdString nodeElementName;
3412      do
3413      {
3414        StdString nodeId("");
3415        if (node.getAttributes().end() != node.getAttributes().find("id"))
3416        { nodeId = node.getAttributes()["id"]; }
3417
3418        nodeElementName = node.getElementName();
3419        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3420        it = transformationMapList_.find(nodeElementName);
3421        if (ite != it)
3422        {
3423          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3424                                                                                                                nodeId,
3425                                                                                                                &node)));
3426        }
3427        else
3428        {
3429          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3430                << "The transformation " << nodeElementName << " has not been supported yet.");
3431        }
3432      } while (node.goToNextElement()) ;
3433      node.goToParentElement();
3434    }
3435  }
3436  CATCH_DUMP_ATTR
3437   //----------------------------------------------------------------
3438
3439   DEFINE_REF_FUNC(Domain,domain)
3440
3441   ///---------------------------------------------------------------
3442
3443} // namespace xios
Note: See TracBrowser for help on using the repository browser.