source: XIOS2/trunk/src/node/domain.cpp

Last change on this file was 2511, checked in by ymipsl, 13 months ago

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