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

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

Coupling branch : replace hasServer and hasClient combination by the name of correct service : CLIENT, GATHERER or OUT_SERVER.

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