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

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

XIOS coupling branch
Some updates.

First coupling test is beginning to work...

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: 111.6 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::checkAttributes(void)
1724   TRY
1725   {
1726      if (this->checkAttributes_done_) return;
1727      this->checkDomain();
1728      this->checkLonLat();
1729      this->checkBounds();
1730      this->checkArea();
1731      this->checkMask();
1732      this->checkDomainData();
1733      this->checkCompression();
1734      this->computeLocalMask() ;
1735      this->completeLonLatClient();
1736      this->checkAttributes_done_ = true;
1737   }
1738   CATCH_DUMP_ATTR
1739
1740   //----------------------------------------------------------------
1741   // Divide function checkAttributes into 2 seperate ones
1742   // This function only checks all attributes of current domain
1743   void CDomain::checkAttributesOnClient()
1744   TRY
1745   {
1746     if (this->isClientChecked) return;
1747     CContext* context=CContext::getCurrent();
1748
1749      if (context->getServiceType()==CServicesManager::CLIENT)
1750      {
1751        this->checkDomain();
1752        this->checkBounds();
1753        this->checkArea();
1754        this->checkLonLat();
1755      }
1756
1757      if (context->getServiceType()==CServicesManager::CLIENT)
1758      { // Ct client uniquement
1759         this->checkMask();
1760         this->checkDomainData();
1761         this->checkCompression();
1762         this->computeLocalMask() ;
1763      }
1764      else
1765      { // Ct serveur uniquement
1766      }
1767
1768      this->isClientChecked = true;
1769   }
1770   CATCH_DUMP_ATTR
1771   
1772   // ym obselete, to be removed
1773   void CDomain::checkAttributesOnClientAfterTransformation()
1774   TRY
1775   {
1776     CContext* context=CContext::getCurrent() ;
1777
1778     if (this->isClientAfterTransformationChecked) return;
1779     if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1780     {
1781       // this->computeConnectedClients();
1782       if (hasLonLat)
1783         if (context->getServiceType()==CServicesManager::CLIENT)
1784           this->completeLonLatClient();
1785     }
1786
1787     this->isClientAfterTransformationChecked = true;
1788   }
1789   CATCH_DUMP_ATTR
1790
1791      // Send all checked attributes to server
1792   void CDomain::sendCheckedAttributes()
1793   TRY
1794   {
1795     if (!this->isClientChecked) checkAttributesOnClient();
1796     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1797     CContext* context=CContext::getCurrent() ;
1798
1799     if (this->isChecked) return;
1800     if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1801     {
1802       sendAttributes();
1803     }
1804     this->isChecked = true;
1805   }
1806   CATCH_DUMP_ATTR
1807
1808/* old version
1809   void CDomain::checkAttributes(void)
1810   TRY
1811   {
1812      if (this->isChecked) return;
1813      CContext* context=CContext::getCurrent() ;
1814
1815      this->checkDomain();
1816      this->checkLonLat();
1817      this->checkBounds();
1818      this->checkArea();
1819
1820      if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1821      { // Ct client uniquement
1822         this->checkMask();
1823         this->checkDomainData();
1824         this->checkCompression();
1825         this->computeLocalMask() ;
1826
1827      }
1828      else
1829      { // Ct serveur uniquement
1830      }
1831
1832      if (context->getServiceType()==CServicesManager::CLIENT || context->getServiceType()==CServicesManager::GATHERER)
1833      {
1834        this->computeConnectedClients();
1835        this->completeLonLatClient();
1836      }
1837
1838      this->isChecked = true;
1839   }
1840   CATCH_DUMP_ATTR
1841*/
1842   
1843  /*!
1844     Compute the connection of a client to other clients to determine which clients to send attributes to.
1845     The sending clients are supposed to already know the distribution of receiving clients (In simple cases, it's band)
1846     The connection among clients is calculated by using global index.
1847     A client connects to other clients which holds the same global index as it.     
1848  */
1849  void CDomain::computeConnectedClients(CContextClient* client)
1850  TRY
1851  {
1852    if (computeConnectedClients_done_.count(client)!=0) return ;
1853    else computeConnectedClients_done_.insert(client) ;
1854   
1855    CContext* context=CContext::getCurrent() ;
1856   
1857    int nbServer = client->serverSize;
1858    int nbClient = client->clientSize;
1859    int rank     = client->clientRank;
1860       
1861    if (listNbServer_.count(nbServer) == 0)
1862    {
1863      listNbServer_.insert(nbServer) ;
1864 
1865      if (connectedServerRank_.find(nbServer) != connectedServerRank_.end())
1866      {
1867        nbSenders.erase(nbServer);
1868        connectedServerRank_.erase(nbServer);
1869      }
1870
1871      if (indSrv_.find(nbServer) == indSrv_.end())
1872      {
1873        int i,j,i_ind,j_ind, nbIndex=i_index.numElements();
1874        int globalIndexCount = i_index.numElements();
1875        // Fill in index
1876        CArray<size_t,1> globalIndexDomain(nbIndex);
1877        size_t globalIndex;
1878
1879        for (i = 0; i < nbIndex; ++i)
1880        {
1881          i_ind=i_index(i);
1882          j_ind=j_index(i);
1883          globalIndex = i_ind + j_ind * ni_glo;
1884          globalIndexDomain(i) = globalIndex;
1885        }
1886
1887        if (globalLocalIndexMap_.empty()) 
1888          for (i = 0; i < nbIndex; ++i)  globalLocalIndexMap_[globalIndexDomain(i)] = i;
1889         
1890
1891        size_t globalSizeIndex = 1, indexBegin, indexEnd;
1892        int range, clientSize = client->clientSize;
1893        std::vector<int> nGlobDomain(2);
1894        nGlobDomain[0] = this->ni_glo;
1895        nGlobDomain[1] = this->nj_glo;
1896        for (int i = 0; i < nGlobDomain.size(); ++i) globalSizeIndex *= nGlobDomain[i];
1897        indexBegin = 0;
1898        if (globalSizeIndex <= clientSize)
1899        {
1900          indexBegin = rank%globalSizeIndex;
1901          indexEnd = indexBegin;
1902        }
1903        else
1904        {
1905          for (int i = 0; i < clientSize; ++i)
1906          {
1907            range = globalSizeIndex / clientSize;
1908            if (i < (globalSizeIndex%clientSize)) ++range;
1909            if (i == client->clientRank) break;
1910            indexBegin += range;
1911          }
1912          indexEnd = indexBegin + range - 1;
1913        }
1914
1915        // Even if servers have no index, they must received something from client
1916        // We only use several client to send "empty" message to these servers
1917        CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
1918        std::vector<int> serverZeroIndex;
1919        if (isUnstructed_) serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 0);
1920        else serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 1);
1921
1922        std::list<int> serverZeroIndexLeader;
1923        std::list<int> serverZeroIndexNotLeader;
1924        CContextClient::computeLeader(client->clientRank, client->clientSize, serverZeroIndex.size(), serverZeroIndexLeader, serverZeroIndexNotLeader);
1925        for (std::list<int>::iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
1926          *it = serverZeroIndex[*it];
1927
1928        CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(), client->intraComm);
1929        clientServerMap->computeServerIndexMapping(globalIndexDomain, nbServer);
1930        CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
1931
1932        CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(), ite = globalIndexDomainOnServer.end();
1933        indSrv_[nbServer].swap(globalIndexDomainOnServer);
1934        connectedServerRank_[nbServer].clear();
1935        for (it = indSrv_[nbServer].begin(); it != ite; ++it) connectedServerRank_[nbServer].push_back(it->first);
1936
1937        for (std::list<int>::const_iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
1938          connectedServerRank_[nbServer].push_back(*it);
1939
1940        // Even if a client has no index, it must connect to at least one server and
1941        // send an "empty" data to this server
1942        if (connectedServerRank_[nbServer].empty())
1943          connectedServerRank_[nbServer].push_back(client->clientRank % client->serverSize);
1944
1945        // Now check if all servers have data to receive. If not, master client will send empty data.
1946        // This ensures that all servers will participate in collective calls upon receiving even if they have no date to receive.
1947        std::vector<int> counts (clientSize);
1948        std::vector<int> displs (clientSize);
1949        displs[0] = 0;
1950        int localCount = connectedServerRank_[nbServer].size() ;
1951        MPI_Gather(&localCount, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, client->intraComm) ;
1952        for (int i = 0; i < clientSize-1; ++i) displs[i+1] = displs[i] + counts[i];
1953        std::vector<int> allConnectedServers(displs[clientSize-1]+counts[clientSize-1]);
1954        MPI_Gatherv(&(connectedServerRank_[nbServer])[0], localCount, MPI_INT, &allConnectedServers[0], &counts[0], &displs[0], MPI_INT, 0, client->intraComm);
1955
1956        if ((allConnectedServers.size() != nbServer) && (rank == 0))
1957        {
1958          std::vector<bool> isSrvConnected (nbServer, false);
1959          for (int i = 0; i < allConnectedServers.size(); ++i) isSrvConnected[allConnectedServers[i]] = true;
1960          for (int i = 0; i < nbServer; ++i) if (!isSrvConnected[i]) connectedServerRank_[nbServer].push_back(i);
1961        }
1962        nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
1963        delete clientServerMap;
1964      }
1965    }
1966  }
1967  CATCH_DUMP_ATTR
1968
1969   /*!
1970     Compute index to write data. We only write data on the zoomed region, therefore, there should
1971     be a map between the complete grid and the reduced grid where we write data.
1972     By using global index we can easily create this kind of mapping.
1973   */
1974   void CDomain::computeWrittenIndex()
1975   TRY
1976   { 
1977      if (computedWrittenIndex_) return;
1978      computedWrittenIndex_ = true;
1979
1980      CContext* context=CContext::getCurrent();     
1981
1982      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1983      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
1984      nSize[0]        = ni;      nSize[1]  = nj;
1985      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1986      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1987      CDistributionServer srvDist(context->intraCommSize_, nBegin, nSize, nBeginGlobal, nGlob); 
1988      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1989
1990      size_t nbWritten = 0, indGlo;     
1991      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
1992                                                          ite = globalLocalIndexMap_.end(), it;         
1993      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
1994                                       itSrve = writtenGlobalIndex.end(), itSrv;
1995
1996      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
1997      nbWritten = 0;
1998      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1999      {
2000        indGlo = *itSrv;
2001        if (ite != globalLocalIndexMap_.find(indGlo))
2002        {
2003          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
2004        }
2005        else
2006        {
2007          localIndexToWriteOnServer(nbWritten) = -1;
2008        }
2009        ++nbWritten;
2010      }
2011   }
2012   CATCH_DUMP_ATTR
2013
2014  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
2015  TRY
2016  {
2017    int writtenCommSize;
2018    MPI_Comm_size(writtenComm, &writtenCommSize);
2019    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
2020      return;
2021
2022    if (isCompressible())
2023    {
2024      size_t nbWritten = 0, indGlo;
2025      CContext* context=CContext::getCurrent();     
2026
2027      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2028      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2029      nSize[0]        = ni;      nSize[1]  = nj;
2030      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2031      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2032      CDistributionServer srvDist(context->intraCommSize_, nBegin, nSize, nBeginGlobal, nGlob); 
2033      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2034
2035      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2036                                                          ite = globalLocalIndexMap_.end(), it;   
2037      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2038                                       itSrve = writtenGlobalIndex.end(), itSrv;
2039      std::unordered_map<size_t,size_t> localGlobalIndexMap;
2040      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2041      {
2042        indGlo = *itSrv;
2043        if (ite != globalLocalIndexMap_.find(indGlo))
2044        {
2045          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
2046          ++nbWritten;
2047        }                 
2048      }
2049
2050      nbWritten = 0;
2051      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2052      {
2053        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2054        {
2055          ++nbWritten;
2056        }
2057      }
2058
2059      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2060      nbWritten = 0;
2061      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2062      {
2063        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2064        {
2065          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2066          ++nbWritten;
2067        }
2068      }
2069
2070      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2071      bool distributed_glo, distributed=isDistributed() ;
2072      MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, writtenComm) ;
2073     
2074      if (distributed_glo)
2075      {
2076             
2077        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2078        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2079        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2080      }
2081      else
2082        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2083      }
2084  }
2085  CATCH_DUMP_ATTR
2086
2087  void CDomain::sendDomainToFileServer(CContextClient* client)
2088  {
2089    if (sendDomainToFileServer_done_.count(client)!=0) return ;
2090    else sendDomainToFileServer_done_.insert(client) ;
2091
2092    StdString domDefRoot("domain_definition");
2093    CDomainGroup* domPtr = CDomainGroup::get(domDefRoot);
2094    domPtr->sendCreateChild(this->getId(), client);
2095    this->sendAllAttributesToServer(client)  ;
2096    this->sendDistributionAttributes(client);   
2097    this->sendIndex(client);       
2098    this->sendLonLat(client);
2099    this->sendArea(client);   
2100    this->sendDataIndex(client);
2101  }
2102
2103  void CDomain::sendDomainToCouplerOut(CContextClient* client, const string& fieldId, int posInGrid)
2104  {
2105    if (sendDomainToFileServer_done_.count(client)!=0) return ;
2106    else sendDomainToFileServer_done_.insert(client) ;
2107   
2108    const string domainId = "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
2109   
2110    if (!domain_ref.isEmpty())
2111    {
2112      auto domain_ref_tmp=domain_ref.getValue() ;
2113      domain_ref.reset() ; // remove the reference, find an other way to do that more cleanly
2114      this->sendAllAttributesToServer(client, domainId)  ;
2115      domain_ref = domain_ref_tmp ;
2116    }
2117    else this->sendAllAttributesToServer(client, domainId)  ;
2118
2119    this->sendDistributionAttributes(client, domainId);   
2120    this->sendIndex(client, domainId);       
2121    this->sendLonLat(client, domainId);
2122    this->sendArea(client, domainId);   
2123    this->sendDataIndex(client, domainId);
2124  }
2125
2126  void CDomain::makeAliasForCoupling(const string& fieldId, int posInGrid)
2127  {
2128    const string domainId = "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
2129    this->createAlias(domainId) ;
2130  }
2131
2132  /*!
2133    Send all attributes from client to connected clients
2134    The attributes will be rebuilt on receiving side
2135  */
2136  // ym obsolete to be removed
2137  void CDomain::sendAttributes()
2138  TRY
2139  {
2140    //sendDistributionAttributes();
2141    //sendIndex();       
2142    //sendLonLat();
2143    //sendArea();   
2144    //sendDataIndex();
2145  }
2146  CATCH
2147  /*!
2148    Send global index from client to connected client(s)
2149  */
2150  void CDomain::sendIndex(CContextClient* client, const string& domainId)
2151  TRY
2152  {
2153    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2154   
2155    int ns, n, i, j, ind, nv, idx;
2156    int serverSize = client->serverSize;
2157    CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2158
2159    list<CMessage> list_msgsIndex;
2160    list<CArray<int,1> > list_indGlob;
2161
2162    std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2163    iteIndex = indSrv_[serverSize].end();
2164    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2165    {
2166      int nbIndGlob = 0;
2167      int rank = connectedServerRank_[serverSize][k];
2168      itIndex = indSrv_[serverSize].find(rank);
2169      if (iteIndex != itIndex)
2170        nbIndGlob = itIndex->second.size();
2171
2172      list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2173
2174      CArray<int,1>& indGlob = list_indGlob.back();
2175      for (n = 0; n < nbIndGlob; ++n) indGlob(n) = static_cast<int>(itIndex->second[n]);
2176
2177      list_msgsIndex.push_back(CMessage());
2178      list_msgsIndex.back() << serverDomainId << (int)type; // enum ne fonctionne pour les message => ToFix
2179      list_msgsIndex.back() << isCurvilinear;
2180      list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2181       
2182      eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2183    }
2184    client->sendEvent(eventIndex);
2185  }
2186  CATCH_DUMP_ATTR
2187
2188  /*!
2189    Send distribution from client to other clients
2190    Because a client in a level knows correctly the grid distribution of client on the next level
2191    it calculates this distribution then sends it to the corresponding clients on the next level
2192  */
2193  void CDomain::sendDistributionAttributes(CContextClient* client, const string& domainId)
2194  TRY
2195  {
2196    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2197
2198    int nbServer = client->serverSize;
2199    std::vector<int> nGlobDomain(2);
2200    nGlobDomain[0] = this->ni_glo;
2201    nGlobDomain[1] = this->nj_glo;
2202
2203    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2204    if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2205    else serverDescription.computeServerDistribution(false, 1);
2206
2207    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2208    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2209
2210    CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2211    if (client->isServerLeader())
2212    {
2213      std::list<CMessage> msgs;
2214
2215      const std::list<int>& ranks = client->getRanksServerLeader();
2216      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2217      {
2218        // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2219        const int ibegin_srv = serverIndexBegin[*itRank][0];
2220        const int jbegin_srv = serverIndexBegin[*itRank][1];
2221        const int ni_srv = serverDimensionSizes[*itRank][0];
2222        const int nj_srv = serverDimensionSizes[*itRank][1];
2223
2224        msgs.push_back(CMessage());
2225        CMessage& msg = msgs.back();
2226        msg << serverDomainId ;
2227        msg << isUnstructed_;
2228        msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2229        msg << ni_glo.getValue() << nj_glo.getValue();
2230        msg << isCompressible_;
2231
2232        event.push(*itRank,1,msg);
2233      }
2234      client->sendEvent(event);
2235    }
2236    else client->sendEvent(event);
2237  }
2238  CATCH_DUMP_ATTR
2239
2240  /*!
2241    Send area from client to connected client(s)
2242  */
2243  void CDomain::sendArea(CContextClient* client, const string& domainId)
2244  TRY
2245  {
2246    if (!hasArea) return;
2247    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2248   
2249    int ns, n, i, j, ind, nv, idx;
2250    int serverSize = client->serverSize;
2251
2252    // send area for each connected server
2253    CEventClient eventArea(getType(), EVENT_ID_AREA);
2254
2255    list<CMessage> list_msgsArea;
2256    list<CArray<double,1> > list_area;
2257
2258    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2259    iteMap = indSrv_[serverSize].end();
2260    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2261    {
2262      int nbData = 0;
2263      int rank = connectedServerRank_[serverSize][k];
2264      it = indSrv_[serverSize].find(rank);
2265      if (iteMap != it)
2266        nbData = it->second.size();
2267      list_area.push_back(CArray<double,1>(nbData));
2268
2269      const std::vector<size_t>& temp = it->second;
2270      for (n = 0; n < nbData; ++n)
2271      {
2272        idx = static_cast<int>(it->second[n]);
2273        list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2274      }
2275
2276      list_msgsArea.push_back(CMessage());
2277      list_msgsArea.back() <<serverDomainId << hasArea;
2278      list_msgsArea.back() << list_area.back();
2279      eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2280    }
2281    client->sendEvent(eventArea);
2282  }
2283  CATCH_DUMP_ATTR
2284
2285  /*!
2286    Send longitude and latitude from client to servers
2287    Each client send long and lat information to corresponding connected clients(s).
2288    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2289  */
2290  void CDomain::sendLonLat(CContextClient* client, const string& domainId)
2291  TRY
2292  {
2293    if (!hasLonLat) return;
2294    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2295   
2296    int ns, n, i, j, ind, nv, idx;
2297    int serverSize = client->serverSize;
2298
2299    // send lon lat for each connected server
2300    CEventClient eventLon(getType(), EVENT_ID_LON);
2301    CEventClient eventLat(getType(), EVENT_ID_LAT);
2302
2303    list<CMessage> list_msgsLon, list_msgsLat;
2304    list<CArray<double,1> > list_lon, list_lat;
2305    list<CArray<double,2> > list_boundslon, list_boundslat;
2306
2307    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2308    iteMap = indSrv_[serverSize].end();
2309    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2310    {
2311      int nbData = 0;
2312      int rank = connectedServerRank_[serverSize][k];
2313      it = indSrv_[serverSize].find(rank);
2314      if (iteMap != it)
2315        nbData = it->second.size();
2316
2317      list_lon.push_back(CArray<double,1>(nbData));
2318      list_lat.push_back(CArray<double,1>(nbData));
2319
2320      if (hasBounds)
2321      {
2322        list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2323        list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2324      }
2325
2326      CArray<double,1>& lon = list_lon.back();
2327      CArray<double,1>& lat = list_lat.back();
2328      const std::vector<size_t>& temp = it->second;
2329      for (n = 0; n < nbData; ++n)
2330      {
2331        idx = static_cast<int>(it->second[n]);
2332        int localInd = globalLocalIndexMap_[idx];
2333        lon(n) = lonvalue(localInd);
2334        lat(n) = latvalue(localInd);
2335
2336        if (hasBounds)
2337        {
2338          CArray<double,2>& boundslon = list_boundslon.back();
2339          CArray<double,2>& boundslat = list_boundslat.back();
2340
2341          for (nv = 0; nv < nvertex; ++nv)
2342          {
2343            boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2344            boundslat(nv, n) = bounds_latvalue(nv, localInd);
2345          }
2346        }
2347      }
2348
2349      list_msgsLon.push_back(CMessage());
2350      list_msgsLat.push_back(CMessage());
2351
2352      list_msgsLon.back() << serverDomainId << hasLonLat;
2353      if (hasLonLat) 
2354        list_msgsLon.back() << list_lon.back();
2355      list_msgsLon.back()  << hasBounds;
2356      if (hasBounds)
2357      {
2358        list_msgsLon.back() << list_boundslon.back();
2359      }
2360
2361      list_msgsLat.back() << serverDomainId << hasLonLat;
2362      if (hasLonLat)
2363        list_msgsLat.back() << list_lat.back();
2364      list_msgsLat.back() << hasBounds;
2365      if (hasBounds)
2366      {         
2367        list_msgsLat.back() << list_boundslat.back();
2368      }
2369
2370      eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2371      eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2372    }
2373    client->sendEvent(eventLon);
2374    client->sendEvent(eventLat);
2375  }
2376  CATCH_DUMP_ATTR
2377
2378  /*!
2379    Send data index to corresponding connected clients.
2380    Data index can be compressed however, we always send decompressed data index
2381    and they will be compressed on receiving.
2382    The compressed index are represented with 1 and others are represented with -1
2383  */
2384  void CDomain::sendDataIndex(CContextClient* client, const string& domainId)
2385  TRY
2386  {
2387    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2388   
2389    int ns, n, i, j, ind, nv, idx;
2390    int serverSize = client->serverSize;
2391
2392    // send area for each connected server
2393    CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2394
2395    list<CMessage> list_msgsDataIndex;
2396    list<CArray<int,1> > list_data_i_index, list_data_j_index;
2397
2398    int nbIndex = i_index.numElements();
2399    int niByIndex = max(i_index) - min(i_index) + 1;
2400    int njByIndex = max(j_index) - min(j_index) + 1; 
2401    int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2402    int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2403
2404   
2405    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2406    dataIIndex = -1; 
2407    dataJIndex = -1;
2408    ind = 0;
2409
2410    for (idx = 0; idx < data_i_index.numElements(); ++idx)
2411    {
2412      int dataIidx = data_i_index(idx) + data_ibegin;
2413      int dataJidx = data_j_index(idx) + data_jbegin;
2414      if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2415          (0 <= dataJidx) && (dataJidx < dataJindexBound))
2416      {
2417        dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2418        dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2419      }
2420    }
2421
2422    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2423    iteMap = indSrv_[serverSize].end();
2424    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2425    {
2426      int nbData = 0;
2427      int rank = connectedServerRank_[serverSize][k];
2428      it = indSrv_[serverSize].find(rank);
2429      if (iteMap != it)
2430        nbData = it->second.size();
2431      list_data_i_index.push_back(CArray<int,1>(nbData));
2432      list_data_j_index.push_back(CArray<int,1>(nbData));
2433
2434      const std::vector<size_t>& temp = it->second;
2435      for (n = 0; n < nbData; ++n)
2436      {
2437        idx = static_cast<int>(it->second[n]);
2438        i = globalLocalIndexMap_[idx];
2439        list_data_i_index.back()(n) = dataIIndex(i);
2440        list_data_j_index.back()(n) = dataJIndex(i);
2441      }
2442
2443      list_msgsDataIndex.push_back(CMessage());
2444      list_msgsDataIndex.back() << serverDomainId ;
2445      list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2446      eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2447    }
2448    client->sendEvent(eventDataIndex);
2449  }
2450  CATCH
2451 
2452  bool CDomain::dispatchEvent(CEventServer& event)
2453  TRY
2454  {
2455    if (SuperClass::dispatchEvent(event)) return true;
2456    else
2457    {
2458      switch(event.type)
2459      {
2460        case EVENT_ID_SERVER_ATTRIBUT:
2461          recvDistributionAttributes(event);
2462          return true;
2463          break;
2464        case EVENT_ID_INDEX:
2465          recvIndex(event);
2466          return true;
2467          break;
2468        case EVENT_ID_LON:
2469          recvLon(event);
2470          return true;
2471          break;
2472        case EVENT_ID_LAT:
2473          recvLat(event);
2474          return true;
2475          break;
2476        case EVENT_ID_AREA:
2477          recvArea(event);
2478          return true;
2479          break; 
2480        case EVENT_ID_DATA_INDEX:
2481          recvDataIndex(event);
2482          return true;
2483          break;
2484        default:
2485          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2486                << "Unknown Event");
2487          return false;
2488       }
2489    }
2490  }
2491  CATCH
2492
2493  /*!
2494    Receive index event from clients(s)
2495    \param[in] event event contain info about rank and associated index
2496  */
2497  void CDomain::recvIndex(CEventServer& event)
2498  TRY
2499  {
2500    string domainId;
2501    std::map<int, CBufferIn*> rankBuffers;
2502
2503    list<CEventServer::SSubEvent>::iterator it;
2504    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2505    {     
2506      CBufferIn* buffer = it->buffer;
2507      *buffer >> domainId;
2508      rankBuffers[it->rank] = buffer;       
2509    }
2510    get(domainId)->recvIndex(rankBuffers);
2511  }
2512  CATCH
2513
2514  /*!
2515    Receive index information from client(s). We use the global index for mapping index between
2516    sending clients and receiving clients.
2517    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2518  */
2519  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2520  TRY
2521  {
2522    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2523    recvClientRanks_.resize(nbReceived);       
2524
2525    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2526    ind = 0;
2527    for (ind = 0; it != ite; ++it, ++ind)
2528    {       
2529       recvClientRanks_[ind] = it->first;
2530       CBufferIn& buffer = *(it->second);
2531       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2532       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2533    }
2534    int nbIndGlob = 0;
2535    for (i = 0; i < nbReceived; ++i)
2536    {
2537      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2538    }
2539   
2540    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2541    i_index.resize(nbIndGlob);
2542    j_index.resize(nbIndGlob);   
2543    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2544
2545    nbIndGlob = 0;
2546    for (i = 0; i < nbReceived; ++i)
2547    {
2548      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2549      for (ind = 0; ind < tmp.numElements(); ++ind)
2550      {
2551         index = tmp(ind);
2552         if (0 == globalLocalIndexMap_.count(index))
2553         {
2554           iIndex = (index%ni_glo)-ibegin;
2555           iIndex = (iIndex < 0) ? 0 : iIndex;
2556           jIndex = (index/ni_glo)-jbegin;
2557           jIndex = (jIndex < 0) ? 0 : jIndex;
2558           nbIndLoc = iIndex + ni * jIndex;
2559           i_index(nbIndGlob) = index % ni_glo;
2560           j_index(nbIndGlob) = index / ni_glo;
2561           globalLocalIndexMap_[index] = nbIndGlob;
2562           ++nbIndGlob;
2563         } 
2564      } 
2565    } 
2566
2567    if (nbIndGlob==0)
2568    {
2569      i_index.resize(nbIndGlob);
2570      j_index.resize(nbIndGlob);
2571    }
2572    else
2573    {
2574      i_index.resizeAndPreserve(nbIndGlob);
2575      j_index.resizeAndPreserve(nbIndGlob);
2576    }
2577
2578    domainMask.resize(0); // Mask is not defined anymore on servers
2579  }
2580  CATCH
2581
2582  /*!
2583    Receive attributes event from clients(s)
2584    \param[in] event event contain info about rank and associated attributes
2585  */
2586  void CDomain::recvDistributionAttributes(CEventServer& event)
2587  TRY
2588  {
2589    CBufferIn* buffer=event.subEvents.begin()->buffer;
2590    string domainId ;
2591    *buffer>>domainId ;
2592    get(domainId)->recvDistributionAttributes(*buffer);
2593  }
2594  CATCH
2595
2596  /*!
2597    Receive attributes from client(s)
2598    \param[in] rank rank of client source
2599    \param[in] buffer message containing attributes info
2600  */
2601  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2602  TRY
2603  {
2604    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2605    int ni_glo_tmp, nj_glo_tmp;
2606    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2607           >> ni_glo_tmp >> nj_glo_tmp
2608           >> isCompressible_;
2609
2610    ni.setValue(ni_tmp);
2611    ibegin.setValue(ibegin_tmp);
2612    nj.setValue(nj_tmp);
2613    jbegin.setValue(jbegin_tmp);
2614    ni_glo.setValue(ni_glo_tmp);
2615    nj_glo.setValue(nj_glo_tmp);
2616
2617  }
2618 CATCH_DUMP_ATTR
2619  /*!
2620    Receive longitude event from clients(s)
2621    \param[in] event event contain info about rank and associated longitude
2622  */
2623  void CDomain::recvLon(CEventServer& event)
2624  TRY
2625  {
2626    string domainId;
2627    std::map<int, CBufferIn*> rankBuffers;
2628
2629    list<CEventServer::SSubEvent>::iterator it;
2630    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2631    {     
2632      CBufferIn* buffer = it->buffer;
2633      *buffer >> domainId;
2634      rankBuffers[it->rank] = buffer;       
2635    }
2636    get(domainId)->recvLon(rankBuffers);
2637  }
2638  CATCH
2639
2640  /*!
2641    Receive longitude information from client(s)
2642    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2643  */
2644  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2645  TRY
2646  {
2647    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2648    if (nbReceived != recvClientRanks_.size())
2649      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2650           << "The number of sending clients is not correct."
2651           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2652   
2653    int nbLonInd = 0;
2654    vector<CArray<double,1> > recvLonValue(nbReceived);
2655    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2656    for (i = 0; i < recvClientRanks_.size(); ++i)
2657    {
2658      int rank = recvClientRanks_[i];
2659      CBufferIn& buffer = *(rankBuffers[rank]);
2660      buffer >> hasLonLat;
2661      if (hasLonLat)
2662        buffer >> recvLonValue[i];
2663      buffer >> hasBounds;
2664      if (hasBounds)
2665        buffer >> recvBoundsLonValue[i];
2666    }
2667
2668    if (hasLonLat)
2669    {
2670      for (i = 0; i < nbReceived; ++i)
2671      {
2672        nbLonInd += recvLonValue[i].numElements();
2673      }
2674   
2675      if (nbLonInd != globalLocalIndexMap_.size())
2676        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2677                 << "something must be wrong with longitude index "<< std::endl;
2678
2679      nbLonInd = globalLocalIndexMap_.size();
2680      lonvalue.resize(nbLonInd);
2681      if (hasBounds)
2682      {
2683        bounds_lonvalue.resize(nvertex,nbLonInd);
2684        bounds_lonvalue = 0.;
2685      }
2686
2687      nbLonInd = 0;
2688      for (i = 0; i < nbReceived; ++i)
2689      {
2690        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2691        CArray<double,1>& tmp = recvLonValue[i];
2692        for (ind = 0; ind < tmp.numElements(); ++ind)
2693        {
2694          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2695          lonvalue(lInd) = tmp(ind); 
2696           if (hasBounds)
2697           {         
2698            for (int nv = 0; nv < nvertex; ++nv)
2699              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2700           }                 
2701        }
2702      }       
2703    }
2704
2705   // setup attribute depending the type of domain
2706    if (hasLonLat)
2707    {
2708      nbLonInd = globalLocalIndexMap_.size();
2709      if (ni*nj != nbLonInd)
2710        ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2711             << "The number of index received is not coherent with the given resolution"
2712             << "nbLonInd=" << nbLonInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
2713
2714      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
2715      {
2716        lonvalue_2d.resize(ni,nj);
2717          for(int ij=0, j=0 ; j<nj ; j++)
2718            for(int i=0 ; i<ni; i++, ij++) lonvalue_2d(i,j) = lonvalue(ij) ;
2719       
2720        if (hasBounds)
2721        {
2722          bounds_lon_2d.resize(nvertex, ni, nj) ;
2723          for(int ij=0, j=0 ; j<nj ; j++)
2724            for(int i=0 ; i<ni; i++, ij++) 
2725              for(int nv=0; nv<nvertex; nv++) bounds_lon_2d(nv,i,j) = bounds_lonvalue(nv,ij) ;
2726        }
2727      }
2728    }
2729    else if (type == type_attr::unstructured || type == type_attr::gaussian)
2730    {
2731      lonvalue_1d.resize(nbLonInd);
2732      lonvalue_1d = lonvalue ;
2733      if (hasBounds)
2734      {
2735        bounds_lon_1d.resize(nvertex, nbLonInd) ;
2736        bounds_lon_1d = bounds_lonvalue ;
2737      }
2738    }
2739  }
2740  CATCH_DUMP_ATTR
2741
2742  /*!
2743    Receive latitude event from clients(s)
2744    \param[in] event event contain info about rank and associated latitude
2745  */
2746  void CDomain::recvLat(CEventServer& event)
2747  TRY
2748  {
2749    string domainId;
2750    std::map<int, CBufferIn*> rankBuffers;
2751
2752    list<CEventServer::SSubEvent>::iterator it;
2753    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2754    {     
2755      CBufferIn* buffer = it->buffer;
2756      *buffer >> domainId;
2757      rankBuffers[it->rank] = buffer;   
2758    }
2759    get(domainId)->recvLat(rankBuffers);
2760  }
2761  CATCH
2762
2763  /*!
2764    Receive latitude information from client(s)
2765    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2766  */
2767  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2768  TRY
2769  {
2770    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2771    if (nbReceived != recvClientRanks_.size())
2772      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2773           << "The number of sending clients is not correct."
2774           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2775
2776    vector<CArray<double,1> > recvLatValue(nbReceived);
2777    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2778    for (i = 0; i < recvClientRanks_.size(); ++i)
2779    {
2780      int rank = recvClientRanks_[i];
2781      CBufferIn& buffer = *(rankBuffers[rank]);
2782      buffer >> hasLonLat;
2783      if (hasLonLat)
2784        buffer >> recvLatValue[i];
2785      buffer >> hasBounds;
2786      if (hasBounds)
2787        buffer >> recvBoundsLatValue[i];
2788    }
2789
2790    int nbLatInd = 0;
2791    if (hasLonLat)
2792    {
2793      for (i = 0; i < nbReceived; ++i)
2794      {
2795        nbLatInd += recvLatValue[i].numElements();
2796      }
2797   
2798      if (nbLatInd != globalLocalIndexMap_.size())
2799        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2800                << "something must be wrong with latitude index "<< std::endl;
2801
2802      nbLatInd = globalLocalIndexMap_.size();
2803      latvalue.resize(nbLatInd);
2804      if (hasBounds)
2805      {
2806        bounds_latvalue.resize(nvertex,nbLatInd);
2807        bounds_latvalue = 0. ;
2808      }
2809
2810      nbLatInd = 0;
2811      for (i = 0; i < nbReceived; ++i)
2812      {
2813        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2814        CArray<double,1>& tmp = recvLatValue[i];
2815        for (ind = 0; ind < tmp.numElements(); ++ind)
2816        {
2817          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2818          latvalue(lInd) = tmp(ind);   
2819           if (hasBounds)
2820           {
2821            CArray<double,2>& boundslat = recvBoundsLatValue[i];
2822            for (int nv = 0; nv < nvertex; ++nv)
2823              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2824           }   
2825          ++nbLatInd;
2826        }
2827      }       
2828    }
2829      // setup attribute depending the type of domain
2830    if (hasLonLat)
2831    {
2832      nbLatInd = globalLocalIndexMap_.size();
2833     
2834      if (ni*nj != nbLatInd)
2835        ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2836             << "The number of index received is not coherent with the given resolution"
2837            << "nbLonInd=" << nbLatInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
2838
2839      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
2840      {
2841        {
2842          latvalue_2d.resize(ni,nj);
2843          for(int ij=0, j=0 ; j<nj ; j++)
2844            for(int i=0 ; i<ni; i++, ij++) latvalue_2d(i,j) = latvalue(ij) ;
2845       
2846          if (hasBounds)
2847          {
2848            bounds_lat_2d.resize(nvertex, ni, nj) ;
2849            for(int ij=0, j=0 ; j<nj ; j++)
2850              for(int i=0 ; i<ni; i++, ij++) 
2851                for(int nv=0; nv<nvertex; nv++) bounds_lat_2d(nv,i,j) = bounds_latvalue(nv,ij) ;
2852          }
2853        }
2854      }
2855      else if (type == type_attr::unstructured || type == type_attr::gaussian)
2856      {
2857        if (hasLonLat)
2858        {
2859          latvalue_1d.resize(nbLatInd);
2860          latvalue_1d = latvalue ;
2861          if (hasBounds)
2862          {
2863            bounds_lat_1d.resize(nvertex, nbLatInd) ;
2864            bounds_lon_1d = bounds_lonvalue ;
2865          }
2866        }
2867      }
2868    } 
2869  }
2870  CATCH_DUMP_ATTR
2871
2872  /*!
2873    Receive area event from clients(s)
2874    \param[in] event event contain info about rank and associated area
2875  */
2876  void CDomain::recvArea(CEventServer& event)
2877  TRY
2878  {
2879    string domainId;
2880    std::map<int, CBufferIn*> rankBuffers;
2881
2882    list<CEventServer::SSubEvent>::iterator it;
2883    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2884    {     
2885      CBufferIn* buffer = it->buffer;
2886      *buffer >> domainId;
2887      rankBuffers[it->rank] = buffer;     
2888    }
2889    get(domainId)->recvArea(rankBuffers);
2890  }
2891  CATCH
2892
2893  /*!
2894    Receive area information from client(s)
2895    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2896  */
2897  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
2898  TRY
2899  {
2900    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2901    if (nbReceived != recvClientRanks_.size())
2902      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2903           << "The number of sending clients is not correct."
2904           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2905
2906    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2907    for (i = 0; i < recvClientRanks_.size(); ++i)
2908    {
2909      int rank = recvClientRanks_[i];
2910      CBufferIn& buffer = *(rankBuffers[rank]);     
2911      buffer >> hasArea;
2912      if (hasArea)
2913        buffer >> recvAreaValue[i];
2914    }
2915
2916    if (hasArea)
2917    {
2918      int nbAreaInd = 0;
2919      for (i = 0; i < nbReceived; ++i)
2920      {     
2921        nbAreaInd += recvAreaValue[i].numElements();
2922      }
2923
2924      if (nbAreaInd != globalLocalIndexMap_.size())
2925        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2926                 << "something must be wrong with area index "<< std::endl;
2927
2928      nbAreaInd = globalLocalIndexMap_.size();
2929      areavalue.resize(nbAreaInd);
2930      nbAreaInd = 0;     
2931      for (i = 0; i < nbReceived; ++i)
2932      {
2933        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2934        CArray<double,1>& tmp = recvAreaValue[i];
2935        for (ind = 0; ind < tmp.numElements(); ++ind)
2936        {
2937          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2938          areavalue(lInd) = tmp(ind);         
2939        }
2940      }
2941     
2942    }
2943  }
2944  CATCH_DUMP_ATTR
2945
2946  /*!
2947    Compare two domain objects.
2948    They are equal if only if they have identical attributes as well as their values.
2949    Moreover, they must have the same transformations.
2950  \param [in] domain Compared domain
2951  \return result of the comparison
2952  */
2953  bool CDomain::isEqual(CDomain* obj)
2954  TRY
2955  {
2956    vector<StdString> excludedAttr;
2957    excludedAttr.push_back("domain_ref");
2958    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2959    if (!objEqual) return objEqual;
2960
2961    TransMapTypes thisTrans = this->getAllTransformations();
2962    TransMapTypes objTrans  = obj->getAllTransformations();
2963
2964    TransMapTypes::const_iterator it, itb, ite;
2965    std::vector<ETranformationType> thisTransType, objTransType;
2966    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2967      thisTransType.push_back(it->first);
2968    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2969      objTransType.push_back(it->first);
2970
2971    if (thisTransType.size() != objTransType.size()) return false;
2972    for (int idx = 0; idx < thisTransType.size(); ++idx)
2973      objEqual &= (thisTransType[idx] == objTransType[idx]);
2974
2975    return objEqual;
2976  }
2977  CATCH_DUMP_ATTR
2978
2979  /*!
2980    Receive data index event from clients(s)
2981    \param[in] event event contain info about rank and associated index
2982  */
2983  void CDomain::recvDataIndex(CEventServer& event)
2984  TRY
2985  {
2986    string domainId;
2987    std::map<int, CBufferIn*> rankBuffers;
2988
2989    list<CEventServer::SSubEvent>::iterator it;
2990    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2991    {     
2992      CBufferIn* buffer = it->buffer;
2993      *buffer >> domainId;
2994      rankBuffers[it->rank] = buffer;       
2995    }
2996    get(domainId)->recvDataIndex(rankBuffers);
2997  }
2998  CATCH
2999
3000  /*!
3001    Receive data index information from client(s)
3002    A client receives data index from different clients to rebuild its own data index.
3003    Because we use global index + mask info to calculate the sending data to client(s),
3004    this data index must be updated with mask info (maybe it will change in the future)
3005    Because the data index is local, to rebuild data index of received client, we should use global index along with.
3006
3007    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3008  */
3009  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
3010  TRY
3011  {
3012    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
3013    if (nbReceived != recvClientRanks_.size())
3014      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3015           << "The number of sending clients is not correct."
3016           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3017
3018    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3019    for (i = 0; i < recvClientRanks_.size(); ++i)
3020    {
3021      int rank = recvClientRanks_[i];
3022      CBufferIn& buffer = *(rankBuffers[rank]);
3023      buffer >> recvDataIIndex[i];
3024      buffer >> recvDataJIndex[i];
3025    }
3026   
3027    int nbIndex = i_index.numElements();
3028    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3029    dataIIndex = -1; dataJIndex = -1;
3030     
3031    nbIndex = 0;
3032    for (i = 0; i < nbReceived; ++i)
3033    {     
3034      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3035      CArray<int,1>& tmpI = recvDataIIndex[i];   
3036      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3037      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3038          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3039             << "The number of global received index is not coherent with the number of received data index."
3040             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3041
3042      for (ind = 0; ind < tmpI.numElements(); ++ind)
3043      {
3044         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3045         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3046         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3047      } 
3048    }
3049
3050    int nbCompressedData = 0; 
3051    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3052    {
3053       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3054       if ((0 <= indexI) && (0 <= indexJ))
3055         ++nbCompressedData;
3056    }       
3057 
3058    data_i_index.resize(nbCompressedData);
3059    data_j_index.resize(nbCompressedData);
3060
3061    nbCompressedData = 0; 
3062    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3063    {
3064       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3065       if ((0 <= indexI) && (0 <= indexJ))
3066       {
3067          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3068          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3069         ++nbCompressedData;
3070       }
3071    }
3072
3073    // Reset data_ibegin, data_jbegin
3074    data_ibegin.setValue(0);
3075    data_jbegin.setValue(0);
3076  }
3077  CATCH_DUMP_ATTR
3078
3079  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3080  TRY
3081  {
3082    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3083    return transformationMap_.back().second;
3084  }
3085  CATCH_DUMP_ATTR
3086
3087  /*!
3088    Check whether a domain has transformation
3089    \return true if domain has transformation
3090  */
3091  bool CDomain::hasTransformation()
3092  TRY
3093  {
3094    return (!transformationMap_.empty());
3095  }
3096  CATCH_DUMP_ATTR
3097
3098  /*!
3099    Set transformation for current domain. It's the method to move transformation in hierarchy
3100    \param [in] domTrans transformation on domain
3101  */
3102  void CDomain::setTransformations(const TransMapTypes& domTrans)
3103  TRY
3104  {
3105    transformationMap_ = domTrans;
3106  }
3107  CATCH_DUMP_ATTR
3108
3109  /*!
3110    Get all transformation current domain has
3111    \return all transformation
3112  */
3113  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3114  TRY
3115  {
3116    return transformationMap_;
3117  }
3118  CATCH_DUMP_ATTR
3119
3120  void CDomain::duplicateTransformation(CDomain* src)
3121  TRY
3122  {
3123    if (src->hasTransformation())
3124    {
3125      this->setTransformations(src->getAllTransformations());
3126    }
3127  }
3128  CATCH_DUMP_ATTR
3129   
3130  /*!
3131   * Go through the hierarchy to find the domain from which the transformations must be inherited
3132   */
3133  void CDomain::solveInheritanceTransformation()
3134  TRY
3135  {
3136    if (hasTransformation() || !hasDirectDomainReference())
3137      return;
3138
3139    CDomain* domain = this;
3140    std::vector<CDomain*> refDomains;
3141    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3142    {
3143      refDomains.push_back(domain);
3144      domain = domain->getDirectDomainReference();
3145    }
3146
3147    if (domain->hasTransformation())
3148      for (size_t i = 0; i < refDomains.size(); ++i)
3149        refDomains[i]->setTransformations(domain->getAllTransformations());
3150  }
3151  CATCH_DUMP_ATTR
3152
3153  void CDomain::setContextClient(CContextClient* contextClient)
3154  TRY
3155  {
3156    if (clientsSet.find(contextClient)==clientsSet.end())
3157    {
3158      clients.push_back(contextClient) ;
3159      clientsSet.insert(contextClient);
3160    }
3161  }
3162  CATCH_DUMP_ATTR
3163
3164  /*!
3165    Parse children nodes of a domain in xml file.
3166    Whenver there is a new transformation, its type and name should be added into this function
3167    \param node child node to process
3168  */
3169  void CDomain::parse(xml::CXMLNode & node)
3170  TRY
3171  {
3172    SuperClass::parse(node);
3173
3174    if (node.goToChildElement())
3175    {
3176      StdString nodeElementName;
3177      do
3178      {
3179        StdString nodeId("");
3180        if (node.getAttributes().end() != node.getAttributes().find("id"))
3181        { nodeId = node.getAttributes()["id"]; }
3182
3183        nodeElementName = node.getElementName();
3184        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3185        it = transformationMapList_.find(nodeElementName);
3186        if (ite != it)
3187        {
3188          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3189                                                                                                                nodeId,
3190                                                                                                                &node)));
3191        }
3192        else
3193        {
3194          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3195                << "The transformation " << nodeElementName << " has not been supported yet.");
3196        }
3197      } while (node.goToNextElement()) ;
3198      node.goToParentElement();
3199    }
3200  }
3201  CATCH_DUMP_ATTR
3202   //----------------------------------------------------------------
3203
3204   DEFINE_REF_FUNC(Domain,domain)
3205
3206   ///---------------------------------------------------------------
3207
3208} // namespace xios
Note: See TracBrowser for help on using the repository browser.