source: XIOS/dev/dev_oa/src/node/domain.cpp @ 2034

Last change on this file since 2034 was 2034, checked in by oabramkina, 20 months ago

Adding a possibility of tiled and non-tiled sent on the same domain.

  • 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: 112.9 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          area.transposeSelf(1, 0);
896     }
897
898     if (ni_glo.isEmpty())
899     {
900       ERROR("CDomain::checkDomain(void)",
901             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
902             << "The global domain is badly defined, "
903             << "the mandatory 'ni_glo' attribute is missing.")
904     }
905     else if (ni_glo <= 0)
906     {
907       ERROR("CDomain::checkDomain(void)",
908             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
909             << "The global domain is badly defined, "
910             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
911     }
912
913     if (nj_glo.isEmpty())
914     {
915       ERROR("CDomain::checkDomain(void)",
916             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
917             << "The global domain is badly defined, "
918             << "the mandatory 'nj_glo' attribute is missing.")
919     }
920     else if (nj_glo <= 0)
921     {
922       ERROR("CDomain::checkDomain(void)",
923             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
924             << "The global domain is badly defined, "
925             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
926     }
927
928     checkLocalIDomain();
929     checkLocalJDomain();
930
931     if (i_index.isEmpty())
932     {
933       i_index.resize(ni*nj);
934       for (int j = 0; j < nj; ++j)
935         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
936     }
937
938     if (j_index.isEmpty())
939     {
940       j_index.resize(ni*nj);
941       for (int j = 0; j < nj; ++j)
942         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
943     }
944   }
945   CATCH_DUMP_ATTR
946
947   size_t CDomain::getGlobalWrittenSize(void)
948   {
949     return ni_glo*nj_glo ;
950   }
951   //----------------------------------------------------------------
952
953   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
954   void CDomain::checkLocalIDomain(void)
955   TRY
956   {
957      // If ibegin and ni are provided then we use them to check the validity of local domain
958      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
959      {
960        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
961        {
962          ERROR("CDomain::checkLocalIDomain(void)",
963                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
964                << "The local domain is wrongly defined,"
965                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
966        }
967      }
968
969      // i_index has higher priority than ibegin and ni
970      if (!i_index.isEmpty())
971      {
972        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
973        if (ni.isEmpty()) 
974        {         
975         // No information about ni
976          int minIndex = ni_glo - 1;
977          int maxIndex = 0;
978          for (int idx = 0; idx < i_index.numElements(); ++idx)
979          {
980            if (i_index(idx) < minIndex) minIndex = i_index(idx);
981            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
982          }
983          ni = maxIndex - minIndex + 1; 
984          minIIndex = minIIndex;         
985        }
986
987        // It's not so correct but if ibegin is not the first value of i_index
988        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
989        if (ibegin.isEmpty()) ibegin = minIIndex;
990      }
991      else if (ibegin.isEmpty() && ni.isEmpty())
992      {
993        ibegin = 0;
994        ni = ni_glo;
995      }
996      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
997      {
998        ERROR("CDomain::checkLocalIDomain(void)",
999              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1000              << "The local domain is wrongly defined," << endl
1001              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
1002              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
1003      }
1004       
1005
1006      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
1007      {
1008        ERROR("CDomain::checkLocalIDomain(void)",
1009              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1010              << "The local domain is wrongly defined,"
1011              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
1012      }
1013   }
1014   CATCH_DUMP_ATTR
1015
1016   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
1017   void CDomain::checkLocalJDomain(void)
1018   TRY
1019   {
1020    // If jbegin and nj are provided then we use them to check the validity of local domain
1021     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
1022     {
1023       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
1024       {
1025         ERROR("CDomain::checkLocalJDomain(void)",
1026                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1027                << "The local domain is wrongly defined,"
1028                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1029       }
1030     }
1031
1032     if (!j_index.isEmpty())
1033     {
1034        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
1035        if (nj.isEmpty()) 
1036        {
1037          // No information about nj
1038          int minIndex = nj_glo - 1;
1039          int maxIndex = 0;
1040          for (int idx = 0; idx < j_index.numElements(); ++idx)
1041          {
1042            if (j_index(idx) < minIndex) minIndex = j_index(idx);
1043            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
1044          }
1045          nj = maxIndex - minIndex + 1;
1046          minJIndex = minIndex; 
1047        } 
1048        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
1049        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
1050       if (jbegin.isEmpty()) jbegin = minJIndex;       
1051     }
1052     else if (jbegin.isEmpty() && nj.isEmpty())
1053     {
1054       jbegin = 0;
1055       nj = nj_glo;
1056     }     
1057
1058
1059     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
1060     {
1061       ERROR("CDomain::checkLocalJDomain(void)",
1062              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1063              << "The local domain is wrongly defined,"
1064              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1065     }
1066   }
1067   CATCH_DUMP_ATTR
1068
1069   //----------------------------------------------------------------
1070
1071   void CDomain::checkMask(void)
1072   TRY
1073   {
1074      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
1075        ERROR("CDomain::checkMask(void)",
1076              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1077              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
1078              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
1079
1080      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
1081      {
1082        if (mask_1d.numElements() != i_index.numElements())
1083          ERROR("CDomain::checkMask(void)",
1084                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1085                << "'mask_1d' does not have the same size as the local domain." << std::endl
1086                << "Local size is " << i_index.numElements() << "." << std::endl
1087                << "Mask size is " << mask_1d.numElements() << ".");
1088      }
1089
1090      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
1091      {
1092        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1093          ERROR("CDomain::checkMask(void)",
1094                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1095                << "The mask does not have the same size as the local domain." << std::endl
1096                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1097                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
1098      }
1099
1100      if (!mask_2d.isEmpty())
1101      {
1102        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
1103        for (int j = 0; j < nj; ++j)
1104          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1105//        mask_2d.reset();
1106      }
1107      else if (mask_1d.isEmpty())
1108      {
1109        domainMask.resize(i_index.numElements());
1110        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
1111      }
1112      else
1113      {
1114      domainMask.resize(mask_1d.numElements());
1115      domainMask=mask_1d ;
1116     }
1117   }
1118   CATCH_DUMP_ATTR
1119
1120   //----------------------------------------------------------------
1121
1122   void CDomain::checkDomainData(void)
1123   TRY
1124   {
1125      if (data_dim.isEmpty())
1126      {
1127        data_dim.setValue(1);
1128      }
1129      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
1130      {
1131        ERROR("CDomain::checkDomainData(void)",
1132              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1133              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
1134      }
1135
1136      if (data_ibegin.isEmpty())
1137         data_ibegin.setValue(0);
1138      if (data_jbegin.isEmpty())
1139         data_jbegin.setValue(0);
1140
1141      if (data_ni.isEmpty())
1142      {
1143        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
1144      }
1145      else if (data_ni.getValue() < 0)
1146      {
1147        ERROR("CDomain::checkDomainData(void)",
1148              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1149              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
1150      }
1151
1152      if (data_nj.isEmpty())
1153      {
1154        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
1155      }
1156      else if (data_nj.getValue() < 0)
1157      {
1158        ERROR("CDomain::checkDomainData(void)",
1159              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1160              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1161      }
1162   }
1163   CATCH_DUMP_ATTR
1164
1165   //----------------------------------------------------------------
1166
1167   void CDomain::checkCompression(void)
1168   TRY
1169   {
1170     int i,j,ind;
1171      if (!data_i_index.isEmpty())
1172      {
1173        if (!data_j_index.isEmpty() &&
1174            data_j_index.numElements() != data_i_index.numElements())
1175        {
1176           ERROR("CDomain::checkCompression(void)",
1177                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1178                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1179                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1180                 << "'data_j_index' size = " << data_j_index.numElements());
1181        }
1182
1183        if (2 == data_dim)
1184        {
1185          if (data_j_index.isEmpty())
1186          {
1187             ERROR("CDomain::checkCompression(void)",
1188                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1189                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
1190          }
1191          for (int k=0; k<data_i_index.numElements(); ++k)
1192          {
1193            i = data_i_index(k)+data_ibegin ;
1194            j = data_j_index(k)+data_jbegin ;
1195            if (i>=0 && i<ni && j>=0 && j<nj)
1196            {
1197              ind=j*ni+i ;
1198              if (!domainMask(ind))
1199              {
1200                data_i_index(k) = -1;
1201                data_j_index(k) = -1;
1202              }
1203            }
1204            else
1205            {
1206              data_i_index(k) = -1;
1207              data_j_index(k) = -1;
1208            }
1209          }
1210        }
1211        else // (1 == data_dim)
1212        {
1213          if (data_j_index.isEmpty())
1214          {
1215            data_j_index.resize(data_ni);
1216            data_j_index = 0;
1217          }
1218          for (int k=0; k<data_i_index.numElements(); ++k)
1219          {
1220            i=data_i_index(k)+data_ibegin ;
1221            if (i>=0 && i < domainMask.size())
1222            {
1223              if (!domainMask(i)) data_i_index(k) = -1;
1224            }
1225            else
1226              data_i_index(k) = -1;
1227
1228            if (!domainMask(i)) data_i_index(k) = -1;
1229          }
1230        }
1231      }
1232      else
1233      {
1234        if (data_dim == 2 && !data_j_index.isEmpty())
1235          ERROR("CDomain::checkCompression(void)",
1236                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1237                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
1238
1239        if (1 == data_dim)
1240        {
1241          data_i_index.resize(data_ni);
1242          data_j_index.resize(data_ni);
1243          data_j_index = 0;
1244
1245          for (int k = 0; k < data_ni; ++k)
1246          {
1247            i=k+data_ibegin ;
1248            if (i>=0 && i < domainMask.size())
1249            {
1250              if (domainMask(i))
1251                data_i_index(k) = k;
1252              else
1253                data_i_index(k) = -1;
1254            }
1255            else
1256              data_i_index(k) = -1;
1257          }
1258        }
1259        else // (data_dim == 2)
1260        {
1261          const int dsize = data_ni * data_nj;
1262          data_i_index.resize(dsize);
1263          data_j_index.resize(dsize);
1264
1265          for(int count = 0, kj = 0; kj < data_nj; ++kj)
1266          {
1267            for(int ki = 0; ki < data_ni; ++ki, ++count)
1268            {
1269              i = ki + data_ibegin;
1270              j = kj + data_jbegin;
1271              ind=j*ni+i ;
1272              if (i>=0 && i<ni && j>=0 && j<nj)
1273              {
1274                if (domainMask(ind))
1275                {
1276                  data_i_index(count) = ki;
1277                  data_j_index(count) = kj;
1278                }
1279                else
1280                {
1281                  data_i_index(count) = -1;
1282                  data_j_index(count) = -1;
1283                }
1284              }
1285              else
1286              {
1287                data_i_index(count) = -1;
1288                data_j_index(count) = -1;
1289              }
1290            }
1291          }
1292        }
1293      }
1294   }
1295   CATCH_DUMP_ATTR
1296
1297   //----------------------------------------------------------------
1298
1299   /*!
1300    * For tiled domains, data_i/j_index should not take into
1301    * account parameters defining data (data_ni/nj, data_i/jbegin...)
1302    * \param [out] dataIndexI
1303    * \param [out] dataIndexJ
1304    * \param [out] infoIndexI
1305    * \param [out] infoIndexJ
1306    */
1307
1308   void CDomain::computeCompressionTiled(CArray<int,1>& dataIndexI, CArray<int,1>& dataIndexJ,
1309                                         CArray<int,1>& infoIndexI, CArray<int,1>& infoIndexJ)
1310   TRY
1311   {
1312     const int dsize = ni * nj;
1313     dataIndexI.resize(dsize);
1314     dataIndexJ.resize(dsize);
1315
1316     dataIndexJ = 0;
1317     for (int k = 0; k < ni; ++k)
1318       dataIndexI(k) = k;
1319
1320     infoIndexI.resize(ni*nj);
1321     for (int j = 0; j < nj; ++j)
1322       for (int i = 0; i < ni; ++i) infoIndexI(i+j*ni) = i+ibegin;
1323
1324     infoIndexJ.resize(ni*nj);
1325     for (int j = 0; j < nj; ++j)
1326       for (int i = 0; i < ni; ++i) infoIndexJ(i+j*ni) = j+jbegin;
1327   }
1328   CATCH_DUMP_ATTR
1329
1330   //----------------------------------------------------------------
1331   void CDomain::computeLocalMask(void)
1332   TRY
1333   {
1334     localMask.resize(i_index.numElements()) ;
1335     localMask=false ;
1336
1337     size_t dn=data_i_index.numElements() ;
1338     int i,j ;
1339     size_t k,ind ;
1340
1341     for(k=0;k<dn;k++)
1342     {
1343       if (data_dim==2)
1344       {
1345          i=data_i_index(k)+data_ibegin ;
1346          j=data_j_index(k)+data_jbegin ;
1347          if (i>=0 && i<ni && j>=0 && j<nj)
1348          {
1349            ind=j*ni+i ;
1350            localMask(ind)=domainMask(ind) ;
1351          }
1352       }
1353       else
1354       {
1355          i=data_i_index(k)+data_ibegin ;
1356          if (i>=0 && i<i_index.numElements())
1357          {
1358            ind=i ;
1359            localMask(ind)=domainMask(ind) ;
1360          }
1361       }
1362     }
1363   }
1364   CATCH_DUMP_ATTR
1365
1366   void CDomain::checkEligibilityForCompressedOutput(void)
1367   TRY
1368   {
1369     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1370     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1371   }
1372   CATCH_DUMP_ATTR
1373
1374   //----------------------------------------------------------------
1375
1376   /*
1377     Fill in longitude, latitude, bounds, and area into internal values (lonvalue, latvalue, bounds_lonvalue, bounds_latvalue, areavalue)
1378     which will be used by XIOS.
1379   */
1380   void CDomain::completeLonLatClient(void)
1381   TRY
1382   {
1383     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1384     checkBounds() ;
1385     checkArea() ;
1386
1387     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
1388     {
1389       lonvalue.resize(ni * nj);
1390       latvalue.resize(ni * nj);
1391       if (hasBounds)
1392       {
1393         bounds_lonvalue.resize(nvertex, ni * nj);
1394         bounds_latvalue.resize(nvertex, ni * nj);
1395       }
1396
1397       for (int j = 0; j < nj; ++j)
1398       {
1399         for (int i = 0; i < ni; ++i)
1400         {
1401           int k = j * ni + i;
1402
1403           lonvalue(k) = lonvalue_2d(i,j);
1404           latvalue(k) = latvalue_2d(i,j);
1405
1406           if (hasBounds)
1407           {
1408             for (int n = 0; n < nvertex; ++n)
1409             {
1410               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1411               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
1412             }
1413           }
1414         }
1415       }
1416     }
1417     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
1418     {
1419       if (type_attr::rectilinear == type)
1420       {
1421         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1422         {
1423           lonvalue.resize(ni * nj);
1424           latvalue.resize(ni * nj);
1425           if (hasBounds)
1426           {
1427             bounds_lonvalue.resize(nvertex, ni * nj);
1428             bounds_latvalue.resize(nvertex, ni * nj);
1429           }
1430
1431           for (int j = 0; j < nj; ++j)
1432           {
1433             for (int i = 0; i < ni; ++i)
1434             {
1435               int k = j * ni + i;
1436
1437               lonvalue(k) = lonvalue_1d(i);
1438               latvalue(k) = latvalue_1d(j);
1439
1440               if (hasBounds)
1441               {
1442                 for (int n = 0; n < nvertex; ++n)
1443                 {
1444                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1445                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1446                 }
1447               }
1448             }
1449           }
1450         }
1451         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1452         {
1453           lonvalue.reference(lonvalue_1d);
1454           latvalue.reference(latvalue_1d);
1455            if (hasBounds)
1456           {
1457             bounds_lonvalue.reference(bounds_lon_1d);
1458             bounds_latvalue.reference(bounds_lat_1d);
1459           }
1460         }
1461         else
1462           ERROR("CDomain::completeLonClient(void)",
1463                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1464                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1465                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1466                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1467                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1468                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1469       }
1470       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1471       {
1472         lonvalue.reference(lonvalue_1d);
1473         latvalue.reference(latvalue_1d);
1474         if (hasBounds)
1475         {
1476           bounds_lonvalue.reference(bounds_lon_1d);
1477           bounds_latvalue.reference(bounds_lat_1d);
1478         }
1479       }
1480     }
1481
1482     if (!area.isEmpty() && areavalue.isEmpty())
1483     {
1484        areavalue.resize(ni*nj);
1485       for (int j = 0; j < nj; ++j)
1486       {
1487         for (int i = 0; i < ni; ++i)
1488         {
1489           int k = j * ni + i;
1490           areavalue(k) = area(i,j);
1491         }
1492       }
1493     }
1494   }
1495   CATCH_DUMP_ATTR
1496
1497   /*
1498     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1499   */
1500   void CDomain::convertLonLatValue(void)
1501   TRY
1502   {
1503     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1504     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1505     {
1506       lonvalue_2d.resize(ni,nj);
1507       latvalue_2d.resize(ni,nj);
1508       if (hasBounds)
1509       {
1510         bounds_lon_2d.resize(nvertex, ni, nj);
1511         bounds_lat_2d.resize(nvertex, ni, nj);
1512       }
1513
1514       for (int j = 0; j < nj; ++j)
1515       {
1516         for (int i = 0; i < ni; ++i)
1517         {
1518           int k = j * ni + i;
1519
1520           lonvalue_2d(i,j) = lonvalue(k);
1521           latvalue_2d(i,j) = latvalue(k);
1522
1523           if (hasBounds)
1524           {
1525             for (int n = 0; n < nvertex; ++n)
1526             {
1527               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1528               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1529             }
1530           }
1531         }
1532       }
1533     }
1534     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1535     {
1536       if (type_attr::rectilinear == type)
1537       {
1538         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1539         {
1540           lonvalue.resize(ni * nj);
1541           latvalue.resize(ni * nj);
1542           if (hasBounds)
1543           {
1544             bounds_lonvalue.resize(nvertex, ni * nj);
1545             bounds_latvalue.resize(nvertex, ni * nj);
1546           }
1547
1548           for (int j = 0; j < nj; ++j)
1549           {
1550             for (int i = 0; i < ni; ++i)
1551             {
1552               int k = j * ni + i;
1553
1554               lonvalue(k) = lonvalue_1d(i);
1555               latvalue(k) = latvalue_1d(j);
1556
1557               if (hasBounds)
1558               {
1559                 for (int n = 0; n < nvertex; ++n)
1560                 {
1561                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1562                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1563                 }
1564               }
1565             }
1566           }
1567         }
1568         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1569         {
1570           lonvalue.reference(lonvalue_1d);
1571           latvalue.reference(latvalue_1d);
1572            if (hasBounds)
1573           {
1574             bounds_lonvalue.reference(bounds_lon_1d);
1575             bounds_latvalue.reference(bounds_lat_1d);
1576           }
1577         }
1578         else
1579           ERROR("CDomain::completeLonClient(void)",
1580                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1581                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1582                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1583                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1584                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1585                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1586       }
1587       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1588       {
1589         lonvalue.reference(lonvalue_1d);
1590         latvalue.reference(latvalue_1d);
1591         if (hasBounds)
1592         {
1593           bounds_lonvalue.reference(bounds_lon_1d);
1594           bounds_latvalue.reference(bounds_lat_1d);
1595         }
1596       }
1597     }
1598   }
1599   CATCH_DUMP_ATTR
1600
1601   void CDomain::checkBounds(void)
1602   TRY
1603   {
1604     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1605     if (!nvertex.isEmpty() && nvertex > 0 && !hasBoundValues)
1606     {
1607       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1608         ERROR("CDomain::checkBounds(void)",
1609               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1610               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1611               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1612
1613       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1614         ERROR("CDomain::checkBounds(void)",
1615               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1616               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1617               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1618
1619       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1620       {
1621         ERROR("CDomain::checkBounds(void)",
1622               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1623               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1624               << "Please define either both attributes or none.");
1625       }
1626
1627       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1628       {
1629         ERROR("CDomain::checkBounds(void)",
1630               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1631               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1632               << "Please define either both attributes or none.");
1633       }
1634
1635       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1636         ERROR("CDomain::checkBounds(void)",
1637               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1638               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1639               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(0)
1640               << " but nvertex is " << nvertex.getValue() << ".");
1641
1642       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1643         ERROR("CDomain::checkBounds(void)",
1644               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1645               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1646               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(0)
1647               << " but nvertex is " << nvertex.getValue() << ".");
1648
1649       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1650         ERROR("CDomain::checkBounds(void)",
1651               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1652               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1653
1654       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1655         ERROR("CDomain::checkBounds(void)",
1656               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1657               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1658
1659       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1660         ERROR("CDomain::checkBounds(void)",
1661               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1662               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1663               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(0)
1664               << " but nvertex is " << nvertex.getValue() << ".");
1665
1666       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1667         ERROR("CDomain::checkBounds(void)",
1668               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1669               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1670               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(0)
1671               << " but nvertex is " << nvertex.getValue() << ".");
1672
1673       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1674         ERROR("CDomain::checkBounds(void)",
1675               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1676               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1677
1678       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1679         ERROR("CDomain::checkBounds(void)",
1680               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1681
1682       // In case of reading UGRID bounds values are not required
1683       hasBounds = (!bounds_lat_1d.isEmpty() || !bounds_lat_2d.isEmpty() );
1684     }
1685     else if (hasBoundValues)
1686     {
1687       hasBounds = true;       
1688     }
1689     else
1690     {
1691       hasBounds = false;
1692     }
1693   }
1694   CATCH_DUMP_ATTR
1695
1696   void CDomain::checkArea(void)
1697   TRY
1698   {
1699     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
1700     hasArea = !area.isEmpty();
1701     if (hasArea && !hasAreaValue)
1702     {
1703       if (area.extent(0) != ni || area.extent(1) != nj)
1704       {
1705         ERROR("CDomain::checkArea(void)",
1706               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1707               << "The area does not have the same size as the local domain." << std::endl
1708               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1709               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1710       }
1711//       if (areavalue.isEmpty())
1712//       {
1713//          areavalue.resize(ni*nj);
1714//         for (int j = 0; j < nj; ++j)
1715//         {
1716//           for (int i = 0; i < ni; ++i)
1717//           {
1718//             int k = j * ni + i;
1719//             areavalue(k) = area(i,j);
1720//           }
1721//         }
1722//       }
1723     }
1724   }
1725   CATCH_DUMP_ATTR
1726
1727   void CDomain::checkLonLat()
1728   TRY
1729   {
1730     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1731                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1732     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1733     if (hasLonLat && !hasLonLatValue)
1734     {
1735       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1736         ERROR("CDomain::checkLonLat()",
1737               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1738               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1739               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1740
1741       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1742       {
1743         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1744           ERROR("CDomain::checkLonLat()",
1745                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1746                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1747                 << "Local size is " << i_index.numElements() << "." << std::endl
1748                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1749       }
1750
1751       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1752       {
1753         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1754           ERROR("CDomain::checkLonLat()",
1755                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1756                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1757                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1758                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1759       }
1760
1761       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1762         ERROR("CDomain::checkLonLat()",
1763               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1764               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1765               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1766
1767       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1768       {
1769         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1770           ERROR("CDomain::checkLonLat()",
1771                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1772                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1773                 << "Local size is " << i_index.numElements() << "." << std::endl
1774                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1775       }
1776
1777       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1778       {
1779         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1780           ERROR("CDomain::checkLonLat()",
1781                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1782                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1783                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1784                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1785       }
1786     }
1787   }
1788   CATCH_DUMP_ATTR
1789
1790   //----------------------------------------------------------------
1791   void CDomain::checkTiles()
1792   TRY
1793   {
1794     if (!ntiles.isEmpty() && ntiles.getValue() >=1) isTiled_ = true;
1795     if (!tile_only.isEmpty() && tile_only.getValue() == true) {
1796       isTiled_ = true;
1797       isTiledOnly_ = true;
1798     }
1799
1800     if (isTiled_)
1801     {
1802
1803       // Attributes tile_ni/nj and tile_i/jbegin are mandatory for tiled domains
1804       if (tile_ni.numElements() != ntiles || tile_ibegin.numElements() != ntiles)
1805           ERROR("CDomain::checkTiles()",
1806                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1807                 << "'tile_ni' or 'tile_ibegin' are ill defined: these attributes must be specified for tiled domains and be of the size 'ntiles'." << std::endl
1808                 << "The number of tiles is  " << ntiles.getValue() << " while the size of 'tile_ni' is " << tile_ni.numElements() 
1809                 << " and the size of 'tile_ibegin' is " << tile_ibegin.numElements() << ".");
1810
1811       if (tile_nj.numElements() != ntiles || tile_jbegin.numElements() != ntiles)
1812           ERROR("CDomain::checkTiles()",
1813                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1814                 << "'tile_nj' or 'tile_jbegin' are ill defined: these attributes must be specified for tiled domains and be of the size 'ntiles'." << std::endl
1815                 << "The number of tiles is  " << ntiles.getValue() << " while the size of 'tile_nj' is " << tile_nj.numElements()
1816                 << " and the size of 'tile_jbegin' is " << tile_jbegin.numElements() << ".");
1817
1818       // Check on consistency of individual tile sizes and local domain size
1819       int sizeTiles = 0;
1820       for (int iTile = 0; iTile < ntiles.getValue(); ++iTile)
1821       {
1822         sizeTiles += tile_ni(iTile) * tile_nj(iTile);
1823       }
1824       if (sizeTiles != (ni*nj))
1825           ERROR("CDomain::checkTiles()",
1826                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1827                 << "tiles should cover the entire local domain and cannot overlap." << std::endl << ".");
1828
1829
1830       // Fill in tile_data_ni/nj and tile_data_i/jbegin if they are not specified
1831       if (tile_data_ni.numElements() == 0 || tile_data_nj.numElements() == 0)
1832       {
1833         tile_data_ni.resize(ntiles);
1834         tile_data_nj.resize(ntiles);
1835         for (int iTile = 0; iTile < ntiles.getValue(); ++iTile)
1836         {
1837           tile_data_ni(iTile) = tile_ni(iTile);
1838           tile_data_nj(iTile) = tile_nj(iTile);
1839         }
1840       }
1841
1842       if (tile_data_ibegin.numElements() == 0 || tile_data_jbegin.numElements() == 0)
1843       {
1844         tile_data_ibegin.resize(ntiles);
1845         tile_data_jbegin.resize(ntiles);
1846         tile_data_ibegin = 0;
1847         tile_data_jbegin = 0;
1848       }
1849
1850     } // isTiled_
1851   }
1852   CATCH_DUMP_ATTR
1853
1854   //----------------------------------------------------------------
1855   int CDomain::getTileDataISize(int tileId) const
1856   TRY
1857   {
1858     int retvalue = (tile_data_ni(tileId) > tile_ni(tileId)) ? tile_data_ni(tileId) : tile_ni(tileId);
1859     return retvalue;
1860   }
1861   CATCH_DUMP_ATTR
1862
1863   //----------------------------------------------------------------
1864   int CDomain::getTileDataJSize(int tileId) const
1865   TRY
1866   {
1867     int retvalue = (tile_data_nj(tileId) > tile_nj(tileId)) ? tile_data_nj(tileId) : tile_nj(tileId);
1868     return retvalue;
1869   }
1870   CATCH_DUMP_ATTR
1871
1872   void CDomain::checkAttributesOnClientAfterTransformation()
1873   TRY
1874   {
1875     CContext* context=CContext::getCurrent() ;
1876
1877     if (this->isClientAfterTransformationChecked) return;
1878     if (context->hasClient)
1879     {
1880      this->computeConnectedClients();
1881       if (hasLonLat)
1882         if (!context->hasServer)
1883           this->completeLonLatClient();
1884     }
1885
1886     this->isClientAfterTransformationChecked = true;
1887   }
1888   CATCH_DUMP_ATTR
1889
1890   //----------------------------------------------------------------
1891   // Divide function checkAttributes into 2 seperate ones
1892   // This function only checks all attributes of current domain
1893   void CDomain::checkAttributesOnClient()
1894   TRY
1895   {
1896     if (this->isClientChecked) return;
1897     CContext* context=CContext::getCurrent();
1898
1899      if (context->hasClient && !context->hasServer)
1900      {
1901        this->checkDomain();
1902        this->checkBounds();
1903        this->checkArea();
1904        this->checkLonLat();
1905        this->checkTiles();
1906      }
1907
1908      if (context->hasClient && !context->hasServer)
1909      { // Ct client uniquement
1910         this->checkMask();
1911         this->checkDomainData();
1912         this->checkCompression();
1913         this->computeLocalMask() ;
1914      }
1915      else
1916      { // Ct serveur uniquement
1917      }
1918
1919      this->isClientChecked = true;
1920   }
1921   CATCH_DUMP_ATTR
1922
1923   // Send all checked attributes to server
1924   void CDomain::sendCheckedAttributes()
1925   TRY
1926   {
1927     if (!this->isClientChecked) checkAttributesOnClient();
1928     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1929     CContext* context=CContext::getCurrent() ;
1930
1931     if (this->isChecked) return;
1932     if (context->hasClient)
1933     {
1934       sendAttributes();
1935     }
1936     this->isChecked = true;
1937   }
1938   CATCH_DUMP_ATTR
1939
1940   void CDomain::checkAttributes(void)
1941   TRY
1942   {
1943      if (this->isChecked) return;
1944      CContext* context=CContext::getCurrent() ;
1945
1946      this->checkDomain();
1947      this->checkLonLat();
1948      this->checkBounds();
1949      this->checkArea();
1950
1951      if (context->hasClient)
1952      { // Ct client uniquement
1953         this->checkMask();
1954         this->checkDomainData();
1955         this->checkCompression();
1956         this->computeLocalMask() ;
1957
1958      }
1959      else
1960      { // Ct serveur uniquement
1961      }
1962
1963      if (context->hasClient)
1964      {
1965        this->computeConnectedClients();
1966        this->completeLonLatClient();
1967      }
1968
1969      this->isChecked = true;
1970   }
1971   CATCH_DUMP_ATTR
1972
1973  /*!
1974     Compute the connection of a client to other clients to determine which clients to send attributes to.
1975     The sending clients are supposed to already know the distribution of receiving clients (In simple cases, it's band)
1976     The connection among clients is calculated by using global index.
1977     A client connects to other clients which holds the same global index as it.     
1978  */
1979  void CDomain::computeConnectedClients()
1980  TRY
1981  {
1982    CContext* context=CContext::getCurrent() ;
1983   
1984    // This line should be changed soon.
1985    int nbSrvPools = (context->hasServer) ? (context->hasClient ? context->clientPrimServer.size() : 0) : 1;
1986
1987    nbSenders.clear();
1988    connectedServerRank_.clear();
1989
1990    for (int p = 0; p < nbSrvPools; ++p)
1991    {
1992      CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[p] : context->client;
1993      int nbServer = client->serverSize;
1994      int nbClient = client->clientSize;
1995      int rank     = client->clientRank;
1996      bool doComputeGlobalIndexServer = true;
1997
1998      if (connectedServerRank_.find(nbServer) == connectedServerRank_.end())
1999      {
2000
2001        if (indSrv_.find(nbServer) == indSrv_.end())
2002        {
2003          int i,j,i_ind,j_ind, nbIndex=i_index.numElements();
2004          int globalIndexCount = i_index.numElements();
2005          // Fill in index
2006          CArray<size_t,1> globalIndexDomain(nbIndex);
2007          size_t globalIndex;
2008
2009          for (i = 0; i < nbIndex; ++i)
2010          {
2011            i_ind=i_index(i);
2012            j_ind=j_index(i);
2013            globalIndex = i_ind + j_ind * ni_glo;
2014            globalIndexDomain(i) = globalIndex;
2015          }
2016
2017          if (globalLocalIndexMap_.empty())
2018          {
2019            for (i = 0; i < nbIndex; ++i)
2020              globalLocalIndexMap_[globalIndexDomain(i)] = i;
2021          }
2022
2023          size_t globalSizeIndex = 1, indexBegin, indexEnd;
2024          int range, clientSize = client->clientSize;
2025          std::vector<int> nGlobDomain(2);
2026          nGlobDomain[0] = this->ni_glo;
2027          nGlobDomain[1] = this->nj_glo;
2028          for (int i = 0; i < nGlobDomain.size(); ++i) globalSizeIndex *= nGlobDomain[i];
2029          indexBegin = 0;
2030          if (globalSizeIndex <= clientSize)
2031          {
2032            indexBegin = rank%globalSizeIndex;
2033            indexEnd = indexBegin;
2034          }
2035          else
2036          {
2037            for (int i = 0; i < clientSize; ++i)
2038            {
2039              range = globalSizeIndex / clientSize;
2040              if (i < (globalSizeIndex%clientSize)) ++range;
2041              if (i == client->clientRank) break;
2042              indexBegin += range;
2043            }
2044            indexEnd = indexBegin + range - 1;
2045          }
2046
2047          // Even if servers have no index, they must received something from client
2048          // We only use several client to send "empty" message to these servers
2049          CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2050          std::vector<int> serverZeroIndex;
2051          if (isUnstructed_) serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 0);
2052          else serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 1);
2053
2054          std::list<int> serverZeroIndexLeader;
2055          std::list<int> serverZeroIndexNotLeader;
2056          CContextClient::computeLeader(client->clientRank, client->clientSize, serverZeroIndex.size(), serverZeroIndexLeader, serverZeroIndexNotLeader);
2057          for (std::list<int>::iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
2058            *it = serverZeroIndex[*it];
2059
2060          CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(), client->intraComm);
2061          clientServerMap->computeServerIndexMapping(globalIndexDomain, nbServer);
2062          CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
2063
2064          CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
2065                 ite = globalIndexDomainOnServer.end();
2066          indSrv_[nbServer].swap(globalIndexDomainOnServer);
2067          connectedServerRank_[nbServer].clear();
2068          for (it = indSrv_[nbServer].begin(); it != ite; ++it)
2069            connectedServerRank_[nbServer].push_back(it->first);
2070
2071          for (std::list<int>::const_iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
2072            connectedServerRank_[nbServer].push_back(*it);
2073
2074          // Even if a client has no index, it must connect to at least one server and
2075          // send an "empty" data to this server
2076          if (connectedServerRank_[nbServer].empty())
2077            connectedServerRank_[nbServer].push_back(client->clientRank % client->serverSize);
2078
2079          // Now check if all servers have data to receive. If not, master client will send empty data.
2080          // This ensures that all servers will participate in collective calls upon receiving even if they have no date to receive.
2081          std::vector<int> counts (clientSize);
2082          std::vector<int> displs (clientSize);
2083          displs[0] = 0;
2084          int localCount = connectedServerRank_[nbServer].size() ;
2085          MPI_Gather(&localCount, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, client->intraComm) ;
2086          for (int i = 0; i < clientSize-1; ++i)
2087          {
2088            displs[i+1] = displs[i] + counts[i];
2089          }
2090          std::vector<int> allConnectedServers(displs[clientSize-1]+counts[clientSize-1]);
2091          MPI_Gatherv(&(connectedServerRank_[nbServer])[0], localCount, MPI_INT, &allConnectedServers[0], &counts[0], &displs[0], MPI_INT, 0, client->intraComm);
2092
2093          if ((allConnectedServers.size() != nbServer) && (rank == 0))
2094          {
2095            std::vector<bool> isSrvConnected (nbServer, false);
2096            for (int i = 0; i < allConnectedServers.size(); ++i) isSrvConnected[allConnectedServers[i]] = true;
2097            for (int i = 0; i < nbServer; ++i)
2098            {
2099              if (!isSrvConnected[i]) connectedServerRank_[nbServer].push_back(i);
2100            }
2101          }
2102          nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
2103          delete clientServerMap;
2104        }
2105      }
2106    }
2107  }
2108  CATCH_DUMP_ATTR
2109
2110   /*!
2111     Compute index to write data. We only write data on the zoomed region, therefore, there should
2112     be a map between the complete grid and the reduced grid where we write data.
2113     By using global index we can easily create this kind of mapping.
2114   */
2115   void CDomain::computeWrittenIndex()
2116   TRY
2117   { 
2118      if (computedWrittenIndex_) return;
2119      computedWrittenIndex_ = true;
2120
2121      CContext* context=CContext::getCurrent();     
2122      CContextServer* server = context->server; 
2123
2124      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2125      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2126      nSize[0]        = ni;      nSize[1]  = nj;
2127      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2128      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2129      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
2130      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2131
2132      size_t nbWritten = 0, indGlo;     
2133      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2134                                                          ite = globalLocalIndexMap_.end(), it;         
2135      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2136                                       itSrve = writtenGlobalIndex.end(), itSrv;
2137
2138      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
2139      nbWritten = 0;
2140      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2141      {
2142        indGlo = *itSrv;
2143        if (ite != globalLocalIndexMap_.find(indGlo))
2144        {
2145          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
2146        }
2147        else
2148        {
2149          localIndexToWriteOnServer(nbWritten) = -1;
2150        }
2151        ++nbWritten;
2152      }
2153   }
2154   CATCH_DUMP_ATTR
2155
2156  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
2157  TRY
2158  {
2159    int writtenCommSize;
2160    MPI_Comm_size(writtenComm, &writtenCommSize);
2161    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
2162      return;
2163
2164    if (isCompressible())
2165    {
2166      size_t nbWritten = 0, indGlo;
2167      CContext* context=CContext::getCurrent();     
2168      CContextServer* server = context->server; 
2169
2170      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2171      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2172      nSize[0]        = ni;      nSize[1]  = nj;
2173      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2174      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2175      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
2176      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2177
2178      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2179                                                          ite = globalLocalIndexMap_.end(), it;   
2180      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2181                                       itSrve = writtenGlobalIndex.end(), itSrv;
2182      std::unordered_map<size_t,size_t> localGlobalIndexMap;
2183      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2184      {
2185        indGlo = *itSrv;
2186        if (ite != globalLocalIndexMap_.find(indGlo))
2187        {
2188          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
2189          ++nbWritten;
2190        }                 
2191      }
2192
2193      nbWritten = 0;
2194      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2195      {
2196        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2197        {
2198          ++nbWritten;
2199        }
2200      }
2201
2202      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2203      nbWritten = 0;
2204      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2205      {
2206        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2207        {
2208          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2209          ++nbWritten;
2210        }
2211      }
2212
2213      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2214      bool distributed_glo, distributed=isDistributed() ;
2215      MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, writtenComm) ;
2216     
2217      if (distributed_glo)
2218      {
2219             
2220        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2221        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2222        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2223      }
2224      else
2225        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2226      }
2227  }
2228  CATCH_DUMP_ATTR
2229
2230  /*!
2231    Send all attributes from client to connected clients
2232    The attributes will be rebuilt on receiving side
2233  */
2234  void CDomain::sendAttributes()
2235  TRY
2236  {
2237    sendDistributionAttributes();
2238    sendIndex();       
2239    sendLonLat();
2240    sendArea();   
2241    sendDataIndex();
2242  }
2243  CATCH
2244  /*!
2245    Send global index from client to connected client(s)
2246  */
2247  void CDomain::sendIndex()
2248  TRY
2249  {
2250    int ns, n, i, j, ind, nv, idx;
2251    std::list<CContextClient*>::iterator it;
2252    for (it=clients.begin(); it!=clients.end(); ++it)
2253    {
2254      CContextClient* client = *it;
2255
2256      int serverSize = client->serverSize;
2257      CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2258
2259      list<CMessage> list_msgsIndex;
2260      list<CArray<int,1> > list_indGlob;
2261
2262      std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2263      iteIndex = indSrv_[serverSize].end();
2264      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2265      {
2266        int nbIndGlob = 0;
2267        int rank = connectedServerRank_[serverSize][k];
2268        itIndex = indSrv_[serverSize].find(rank);
2269        if (iteIndex != itIndex)
2270          nbIndGlob = itIndex->second.size();
2271
2272        list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2273
2274        CArray<int,1>& indGlob = list_indGlob.back();
2275        for (n = 0; n < nbIndGlob; ++n)
2276        {
2277          indGlob(n) = static_cast<int>(itIndex->second[n]);
2278        }
2279
2280        list_msgsIndex.push_back(CMessage());
2281        list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
2282        list_msgsIndex.back() << isCurvilinear;
2283        list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2284       
2285        eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2286      }
2287
2288      client->sendEvent(eventIndex);
2289    }
2290  }
2291  CATCH_DUMP_ATTR
2292
2293  /*!
2294    Send distribution from client to other clients
2295    Because a client in a level knows correctly the grid distribution of client on the next level
2296    it calculates this distribution then sends it to the corresponding clients on the next level
2297  */
2298  void CDomain::sendDistributionAttributes(void)
2299  TRY
2300  {
2301    std::list<CContextClient*>::iterator it;
2302    for (it=clients.begin(); it!=clients.end(); ++it)
2303    {
2304      CContextClient* client = *it;
2305      int nbServer = client->serverSize;
2306      std::vector<int> nGlobDomain(2);
2307      nGlobDomain[0] = this->ni_glo;
2308      nGlobDomain[1] = this->nj_glo;
2309
2310      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2311      if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2312      else serverDescription.computeServerDistribution(false, 1);
2313
2314      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2315      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2316
2317      CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2318      if (client->isServerLeader())
2319      {
2320        std::list<CMessage> msgs;
2321
2322        const std::list<int>& ranks = client->getRanksServerLeader();
2323        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2324        {
2325          // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2326          const int ibegin_srv = serverIndexBegin[*itRank][0];
2327          const int jbegin_srv = serverIndexBegin[*itRank][1];
2328          const int ni_srv = serverDimensionSizes[*itRank][0];
2329          const int nj_srv = serverDimensionSizes[*itRank][1];
2330
2331          msgs.push_back(CMessage());
2332          CMessage& msg = msgs.back();
2333          msg << this->getId() ;
2334          msg << isUnstructed_;
2335          msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2336          msg << ni_glo.getValue() << nj_glo.getValue();
2337          msg << isCompressible_;
2338
2339          event.push(*itRank,1,msg);
2340        }
2341        client->sendEvent(event);
2342      }
2343      else client->sendEvent(event);
2344    }
2345  }
2346  CATCH_DUMP_ATTR
2347
2348  /*!
2349    Send area from client to connected client(s)
2350  */
2351  void CDomain::sendArea()
2352  TRY
2353  {
2354    if (!hasArea) return;
2355
2356    int ns, n, i, j, ind, nv, idx;
2357    std::list<CContextClient*>::iterator it;
2358
2359    for (it=clients.begin(); it!=clients.end(); ++it)
2360    {
2361      CContextClient* client = *it;
2362      int serverSize = client->serverSize;
2363
2364      // send area for each connected server
2365      CEventClient eventArea(getType(), EVENT_ID_AREA);
2366
2367      list<CMessage> list_msgsArea;
2368      list<CArray<double,1> > list_area;
2369
2370      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2371      iteMap = indSrv_[serverSize].end();
2372      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2373      {
2374        int nbData = 0;
2375        int rank = connectedServerRank_[serverSize][k];
2376        it = indSrv_[serverSize].find(rank);
2377        if (iteMap != it)
2378          nbData = it->second.size();
2379        list_area.push_back(CArray<double,1>(nbData));
2380
2381        const std::vector<size_t>& temp = it->second;
2382        for (n = 0; n < nbData; ++n)
2383        {
2384          idx = static_cast<int>(it->second[n]);
2385          list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2386        }
2387
2388        list_msgsArea.push_back(CMessage());
2389        list_msgsArea.back() << this->getId() << hasArea;
2390        list_msgsArea.back() << list_area.back();
2391        eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2392      }
2393      client->sendEvent(eventArea);
2394    }
2395  }
2396  CATCH_DUMP_ATTR
2397
2398  /*!
2399    Send longitude and latitude from client to servers
2400    Each client send long and lat information to corresponding connected clients(s).
2401    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2402  */
2403  void CDomain::sendLonLat()
2404  TRY
2405  {
2406    if (!hasLonLat) return;
2407
2408    int ns, n, i, j, ind, nv, idx;
2409    std::list<CContextClient*>::iterator it;
2410    for (it=clients.begin(); it!=clients.end(); ++it)
2411    {
2412      CContextClient* client = *it;
2413      int serverSize = client->serverSize;
2414
2415      // send lon lat for each connected server
2416      CEventClient eventLon(getType(), EVENT_ID_LON);
2417      CEventClient eventLat(getType(), EVENT_ID_LAT);
2418
2419      list<CMessage> list_msgsLon, list_msgsLat;
2420      list<CArray<double,1> > list_lon, list_lat;
2421      list<CArray<double,2> > list_boundslon, list_boundslat;
2422
2423      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2424      iteMap = indSrv_[serverSize].end();
2425      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2426      {
2427        int nbData = 0;
2428        int rank = connectedServerRank_[serverSize][k];
2429        it = indSrv_[serverSize].find(rank);
2430        if (iteMap != it)
2431          nbData = it->second.size();
2432
2433        list_lon.push_back(CArray<double,1>(nbData));
2434        list_lat.push_back(CArray<double,1>(nbData));
2435
2436        if (hasBounds)
2437        {
2438          list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2439          list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2440        }
2441
2442        CArray<double,1>& lon = list_lon.back();
2443        CArray<double,1>& lat = list_lat.back();
2444        const std::vector<size_t>& temp = it->second;
2445        for (n = 0; n < nbData; ++n)
2446        {
2447          idx = static_cast<int>(it->second[n]);
2448          int localInd = globalLocalIndexMap_[idx];
2449          lon(n) = lonvalue(localInd);
2450          lat(n) = latvalue(localInd);
2451
2452          if (hasBounds)
2453          {
2454            CArray<double,2>& boundslon = list_boundslon.back();
2455            CArray<double,2>& boundslat = list_boundslat.back();
2456
2457            for (nv = 0; nv < nvertex; ++nv)
2458            {
2459              boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2460              boundslat(nv, n) = bounds_latvalue(nv, localInd);
2461            }
2462          }
2463        }
2464
2465        list_msgsLon.push_back(CMessage());
2466        list_msgsLat.push_back(CMessage());
2467
2468        list_msgsLon.back() << this->getId() << hasLonLat;
2469        if (hasLonLat) 
2470          list_msgsLon.back() << list_lon.back();
2471        list_msgsLon.back()  << hasBounds;
2472        if (hasBounds)
2473        {
2474          list_msgsLon.back() << list_boundslon.back();
2475        }
2476
2477        list_msgsLat.back() << this->getId() << hasLonLat;
2478        if (hasLonLat)
2479          list_msgsLat.back() << list_lat.back();
2480        list_msgsLat.back() << hasBounds;
2481        if (hasBounds)
2482        {         
2483          list_msgsLat.back() << list_boundslat.back();
2484        }
2485
2486        eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2487        eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2488      }
2489      client->sendEvent(eventLon);
2490      client->sendEvent(eventLat);
2491    }
2492  }
2493  CATCH_DUMP_ATTR
2494
2495  /*!
2496    Send data index to corresponding connected clients.
2497    Data index can be compressed however, we always send decompressed data index
2498    and they will be compressed on receiving.
2499    The compressed index are represented with 1 and others are represented with -1
2500  */
2501  void CDomain::sendDataIndex()
2502  TRY
2503  {
2504    int ns, n, i, j, ind, nv, idx;
2505    std::list<CContextClient*>::iterator it;
2506    for (it=clients.begin(); it!=clients.end(); ++it)
2507    {
2508      CContextClient* client = *it;
2509
2510      int serverSize = client->serverSize;
2511
2512      // send area for each connected server
2513      CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2514
2515      list<CMessage> list_msgsDataIndex;
2516      list<CArray<int,1> > list_data_i_index, list_data_j_index;
2517
2518      int nbIndex = i_index.numElements();
2519      int niByIndex = max(i_index) - min(i_index) + 1;
2520      int njByIndex = max(j_index) - min(j_index) + 1; 
2521      int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2522      int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2523
2524     
2525      CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2526      dataIIndex = -1; 
2527      dataJIndex = -1;
2528      ind = 0;
2529
2530      for (idx = 0; idx < data_i_index.numElements(); ++idx)
2531      {
2532        int dataIidx = data_i_index(idx) + data_ibegin;
2533        int dataJidx = data_j_index(idx) + data_jbegin;
2534        if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2535            (0 <= dataJidx) && (dataJidx < dataJindexBound))
2536        {
2537          dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2538          dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2539        }
2540      }
2541
2542      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2543      iteMap = indSrv_[serverSize].end();
2544      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2545      {
2546        int nbData = 0;
2547        int rank = connectedServerRank_[serverSize][k];
2548        it = indSrv_[serverSize].find(rank);
2549        if (iteMap != it)
2550          nbData = it->second.size();
2551        list_data_i_index.push_back(CArray<int,1>(nbData));
2552        list_data_j_index.push_back(CArray<int,1>(nbData));
2553
2554        const std::vector<size_t>& temp = it->second;
2555        for (n = 0; n < nbData; ++n)
2556        {
2557          idx = static_cast<int>(it->second[n]);
2558          i = globalLocalIndexMap_[idx];
2559          list_data_i_index.back()(n) = dataIIndex(i);
2560          list_data_j_index.back()(n) = dataJIndex(i);
2561        }
2562
2563        list_msgsDataIndex.push_back(CMessage());
2564        list_msgsDataIndex.back() << this->getId();
2565        list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2566        eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2567      }
2568      client->sendEvent(eventDataIndex);
2569    }
2570  }
2571  CATCH
2572 
2573  bool CDomain::dispatchEvent(CEventServer& event)
2574  TRY
2575  {
2576    if (SuperClass::dispatchEvent(event)) return true;
2577    else
2578    {
2579      switch(event.type)
2580      {
2581        case EVENT_ID_SERVER_ATTRIBUT:
2582          recvDistributionAttributes(event);
2583          return true;
2584          break;
2585        case EVENT_ID_INDEX:
2586          recvIndex(event);
2587          return true;
2588          break;
2589        case EVENT_ID_LON:
2590          recvLon(event);
2591          return true;
2592          break;
2593        case EVENT_ID_LAT:
2594          recvLat(event);
2595          return true;
2596          break;
2597        case EVENT_ID_AREA:
2598          recvArea(event);
2599          return true;
2600          break; 
2601        case EVENT_ID_DATA_INDEX:
2602          recvDataIndex(event);
2603          return true;
2604          break;
2605        default:
2606          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2607                << "Unknown Event");
2608          return false;
2609       }
2610    }
2611  }
2612  CATCH
2613
2614  /*!
2615    Receive index event from clients(s)
2616    \param[in] event event contain info about rank and associated index
2617  */
2618  void CDomain::recvIndex(CEventServer& event)
2619  TRY
2620  {
2621    string domainId;
2622    std::map<int, CBufferIn*> rankBuffers;
2623
2624    list<CEventServer::SSubEvent>::iterator it;
2625    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2626    {     
2627      CBufferIn* buffer = it->buffer;
2628      *buffer >> domainId;
2629      rankBuffers[it->rank] = buffer;       
2630    }
2631    get(domainId)->recvIndex(rankBuffers);
2632  }
2633  CATCH
2634
2635  /*!
2636    Receive index information from client(s). We use the global index for mapping index between
2637    sending clients and receiving clients.
2638    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2639  */
2640  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2641  TRY
2642  {
2643    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2644    recvClientRanks_.resize(nbReceived);       
2645
2646    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2647    ind = 0;
2648    for (ind = 0; it != ite; ++it, ++ind)
2649    {       
2650       recvClientRanks_[ind] = it->first;
2651       CBufferIn& buffer = *(it->second);
2652       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2653       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2654    }
2655    int nbIndGlob = 0;
2656    for (i = 0; i < nbReceived; ++i)
2657    {
2658      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2659    }
2660   
2661    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2662    i_index.resize(nbIndGlob);
2663    j_index.resize(nbIndGlob);   
2664    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2665
2666    nbIndGlob = 0;
2667    for (i = 0; i < nbReceived; ++i)
2668    {
2669      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2670      for (ind = 0; ind < tmp.numElements(); ++ind)
2671      {
2672         index = tmp(ind);
2673         if (0 == globalLocalIndexMap_.count(index))
2674         {
2675           iIndex = (index%ni_glo)-ibegin;
2676           iIndex = (iIndex < 0) ? 0 : iIndex;
2677           jIndex = (index/ni_glo)-jbegin;
2678           jIndex = (jIndex < 0) ? 0 : jIndex;
2679           nbIndLoc = iIndex + ni * jIndex;
2680           i_index(nbIndGlob) = index % ni_glo;
2681           j_index(nbIndGlob) = index / ni_glo;
2682           globalLocalIndexMap_[index] = nbIndGlob;
2683           ++nbIndGlob;
2684         } 
2685      } 
2686    } 
2687
2688    if (nbIndGlob==0)
2689    {
2690      i_index.resize(nbIndGlob);
2691      j_index.resize(nbIndGlob);
2692    }
2693    else
2694    {
2695      i_index.resizeAndPreserve(nbIndGlob);
2696      j_index.resizeAndPreserve(nbIndGlob);
2697    }
2698
2699    domainMask.resize(0); // Mask is not defined anymore on servers
2700  }
2701  CATCH
2702
2703  /*!
2704    Receive attributes event from clients(s)
2705    \param[in] event event contain info about rank and associated attributes
2706  */
2707  void CDomain::recvDistributionAttributes(CEventServer& event)
2708  TRY
2709  {
2710    CBufferIn* buffer=event.subEvents.begin()->buffer;
2711    string domainId ;
2712    *buffer>>domainId ;
2713    get(domainId)->recvDistributionAttributes(*buffer);
2714  }
2715  CATCH
2716
2717  /*!
2718    Receive attributes from client(s)
2719    \param[in] rank rank of client source
2720    \param[in] buffer message containing attributes info
2721  */
2722  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2723  TRY
2724  {
2725    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2726    int ni_glo_tmp, nj_glo_tmp;
2727    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2728           >> ni_glo_tmp >> nj_glo_tmp
2729           >> isCompressible_;
2730
2731    ni.setValue(ni_tmp);
2732    ibegin.setValue(ibegin_tmp);
2733    nj.setValue(nj_tmp);
2734    jbegin.setValue(jbegin_tmp);
2735    ni_glo.setValue(ni_glo_tmp);
2736    nj_glo.setValue(nj_glo_tmp);
2737
2738  }
2739 CATCH_DUMP_ATTR
2740  /*!
2741    Receive longitude event from clients(s)
2742    \param[in] event event contain info about rank and associated longitude
2743  */
2744  void CDomain::recvLon(CEventServer& event)
2745  TRY
2746  {
2747    string domainId;
2748    std::map<int, CBufferIn*> rankBuffers;
2749
2750    list<CEventServer::SSubEvent>::iterator it;
2751    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2752    {     
2753      CBufferIn* buffer = it->buffer;
2754      *buffer >> domainId;
2755      rankBuffers[it->rank] = buffer;       
2756    }
2757    get(domainId)->recvLon(rankBuffers);
2758  }
2759  CATCH
2760
2761  /*!
2762    Receive longitude information from client(s)
2763    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2764  */
2765  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2766  TRY
2767  {
2768    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2769    if (nbReceived != recvClientRanks_.size())
2770      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2771           << "The number of sending clients is not correct."
2772           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2773
2774    vector<CArray<double,1> > recvLonValue(nbReceived);
2775    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2776    for (i = 0; i < recvClientRanks_.size(); ++i)
2777    {
2778      int rank = recvClientRanks_[i];
2779      CBufferIn& buffer = *(rankBuffers[rank]);
2780      buffer >> hasLonLat;
2781      if (hasLonLat)
2782        buffer >> recvLonValue[i];
2783      buffer >> hasBounds;
2784      if (hasBounds)
2785        buffer >> recvBoundsLonValue[i];
2786    }
2787
2788    if (hasLonLat)
2789    {
2790      int nbLonInd = 0;
2791      for (i = 0; i < nbReceived; ++i)
2792      {
2793        nbLonInd += recvLonValue[i].numElements();
2794      }
2795   
2796      if (nbLonInd != globalLocalIndexMap_.size())
2797        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2798                 << "something must be wrong with longitude index "<< std::endl;
2799
2800      nbLonInd = globalLocalIndexMap_.size();
2801      lonvalue.resize(nbLonInd);
2802      if (hasBounds)
2803      {
2804        bounds_lonvalue.resize(nvertex,nbLonInd);
2805        bounds_lonvalue = 0.;
2806      }
2807
2808      nbLonInd = 0;
2809      for (i = 0; i < nbReceived; ++i)
2810      {
2811        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2812        CArray<double,1>& tmp = recvLonValue[i];
2813        for (ind = 0; ind < tmp.numElements(); ++ind)
2814        {
2815          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2816          lonvalue(lInd) = tmp(ind); 
2817           if (hasBounds)
2818           {         
2819            for (int nv = 0; nv < nvertex; ++nv)
2820              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2821           }                 
2822        }
2823      }       
2824    }
2825  }
2826  CATCH_DUMP_ATTR
2827
2828  /*!
2829    Receive latitude event from clients(s)
2830    \param[in] event event contain info about rank and associated latitude
2831  */
2832  void CDomain::recvLat(CEventServer& event)
2833  TRY
2834  {
2835    string domainId;
2836    std::map<int, CBufferIn*> rankBuffers;
2837
2838    list<CEventServer::SSubEvent>::iterator it;
2839    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2840    {     
2841      CBufferIn* buffer = it->buffer;
2842      *buffer >> domainId;
2843      rankBuffers[it->rank] = buffer;   
2844    }
2845    get(domainId)->recvLat(rankBuffers);
2846  }
2847  CATCH
2848
2849  /*!
2850    Receive latitude information from client(s)
2851    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2852  */
2853  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2854  TRY
2855  {
2856    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2857    if (nbReceived != recvClientRanks_.size())
2858      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2859           << "The number of sending clients is not correct."
2860           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2861
2862    vector<CArray<double,1> > recvLatValue(nbReceived);
2863    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2864    for (i = 0; i < recvClientRanks_.size(); ++i)
2865    {
2866      int rank = recvClientRanks_[i];
2867      CBufferIn& buffer = *(rankBuffers[rank]);
2868      buffer >> hasLonLat;
2869      if (hasLonLat)
2870        buffer >> recvLatValue[i];
2871      buffer >> hasBounds;
2872      if (hasBounds)
2873        buffer >> recvBoundsLatValue[i];
2874    }
2875
2876    if (hasLonLat)
2877    {
2878      int nbLatInd = 0;
2879      for (i = 0; i < nbReceived; ++i)
2880      {
2881        nbLatInd += recvLatValue[i].numElements();
2882      }
2883   
2884      if (nbLatInd != globalLocalIndexMap_.size())
2885        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2886                << "something must be wrong with latitude index "<< std::endl;
2887
2888      nbLatInd = globalLocalIndexMap_.size();
2889      latvalue.resize(nbLatInd);
2890      if (hasBounds)
2891      {
2892        bounds_latvalue.resize(nvertex,nbLatInd);
2893        bounds_latvalue = 0. ;
2894      }
2895
2896      nbLatInd = 0;
2897      for (i = 0; i < nbReceived; ++i)
2898      {
2899        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2900        CArray<double,1>& tmp = recvLatValue[i];
2901        for (ind = 0; ind < tmp.numElements(); ++ind)
2902        {
2903          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2904          latvalue(lInd) = tmp(ind);   
2905           if (hasBounds)
2906           {
2907            CArray<double,2>& boundslat = recvBoundsLatValue[i];
2908            for (int nv = 0; nv < nvertex; ++nv)
2909              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2910           }   
2911          ++nbLatInd;
2912        }
2913      }       
2914    }
2915  }
2916  CATCH_DUMP_ATTR
2917
2918  /*!
2919    Receive area event from clients(s)
2920    \param[in] event event contain info about rank and associated area
2921  */
2922  void CDomain::recvArea(CEventServer& event)
2923  TRY
2924  {
2925    string domainId;
2926    std::map<int, CBufferIn*> rankBuffers;
2927
2928    list<CEventServer::SSubEvent>::iterator it;
2929    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2930    {     
2931      CBufferIn* buffer = it->buffer;
2932      *buffer >> domainId;
2933      rankBuffers[it->rank] = buffer;     
2934    }
2935    get(domainId)->recvArea(rankBuffers);
2936  }
2937  CATCH
2938
2939  /*!
2940    Receive area information from client(s)
2941    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2942  */
2943  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
2944  TRY
2945  {
2946    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2947    if (nbReceived != recvClientRanks_.size())
2948      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2949           << "The number of sending clients is not correct."
2950           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2951
2952    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2953    for (i = 0; i < recvClientRanks_.size(); ++i)
2954    {
2955      int rank = recvClientRanks_[i];
2956      CBufferIn& buffer = *(rankBuffers[rank]);     
2957      buffer >> hasArea;
2958      if (hasArea)
2959        buffer >> recvAreaValue[i];
2960    }
2961
2962    if (hasArea)
2963    {
2964      int nbAreaInd = 0;
2965      for (i = 0; i < nbReceived; ++i)
2966      {     
2967        nbAreaInd += recvAreaValue[i].numElements();
2968      }
2969
2970      if (nbAreaInd != globalLocalIndexMap_.size())
2971        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2972                 << "something must be wrong with area index "<< std::endl;
2973
2974      nbAreaInd = globalLocalIndexMap_.size();
2975      areavalue.resize(nbAreaInd);
2976      nbAreaInd = 0;     
2977      for (i = 0; i < nbReceived; ++i)
2978      {
2979        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2980        CArray<double,1>& tmp = recvAreaValue[i];
2981        for (ind = 0; ind < tmp.numElements(); ++ind)
2982        {
2983          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2984          areavalue(lInd) = tmp(ind);         
2985        }
2986      }
2987     
2988    }
2989  }
2990  CATCH_DUMP_ATTR
2991
2992  /*!
2993    Compare two domain objects.
2994    They are equal if only if they have identical attributes as well as their values.
2995    Moreover, they must have the same transformations.
2996  \param [in] domain Compared domain
2997  \return result of the comparison
2998  */
2999  bool CDomain::isEqual(CDomain* obj)
3000  TRY
3001  {
3002    vector<StdString> excludedAttr;
3003    excludedAttr.push_back("domain_ref");
3004    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
3005    if (!objEqual) return objEqual;
3006
3007    TransMapTypes thisTrans = this->getAllTransformations();
3008    TransMapTypes objTrans  = obj->getAllTransformations();
3009
3010    TransMapTypes::const_iterator it, itb, ite;
3011    std::vector<ETranformationType> thisTransType, objTransType;
3012    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
3013      thisTransType.push_back(it->first);
3014    for (it = objTrans.begin(); it != objTrans.end(); ++it)
3015      objTransType.push_back(it->first);
3016
3017    if (thisTransType.size() != objTransType.size()) return false;
3018    for (int idx = 0; idx < thisTransType.size(); ++idx)
3019      objEqual &= (thisTransType[idx] == objTransType[idx]);
3020
3021    return objEqual;
3022  }
3023  CATCH_DUMP_ATTR
3024
3025  /*!
3026    Receive data index event from clients(s)
3027    \param[in] event event contain info about rank and associated index
3028  */
3029  void CDomain::recvDataIndex(CEventServer& event)
3030  TRY
3031  {
3032    string domainId;
3033    std::map<int, CBufferIn*> rankBuffers;
3034
3035    list<CEventServer::SSubEvent>::iterator it;
3036    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3037    {     
3038      CBufferIn* buffer = it->buffer;
3039      *buffer >> domainId;
3040      rankBuffers[it->rank] = buffer;       
3041    }
3042    get(domainId)->recvDataIndex(rankBuffers);
3043  }
3044  CATCH
3045
3046  /*!
3047    Receive data index information from client(s)
3048    A client receives data index from different clients to rebuild its own data index.
3049    Because we use global index + mask info to calculate the sending data to client(s),
3050    this data index must be updated with mask info (maybe it will change in the future)
3051    Because the data index is local, to rebuild data index of received client, we should use global index along with.
3052
3053    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3054  */
3055  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
3056  TRY
3057  {
3058    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
3059    if (nbReceived != recvClientRanks_.size())
3060      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3061           << "The number of sending clients is not correct."
3062           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3063
3064    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3065    for (i = 0; i < recvClientRanks_.size(); ++i)
3066    {
3067      int rank = recvClientRanks_[i];
3068      CBufferIn& buffer = *(rankBuffers[rank]);
3069      buffer >> recvDataIIndex[i];
3070      buffer >> recvDataJIndex[i];
3071    }
3072   
3073    int nbIndex = i_index.numElements();
3074    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3075    dataIIndex = -1; dataJIndex = -1;
3076     
3077    nbIndex = 0;
3078    for (i = 0; i < nbReceived; ++i)
3079    {     
3080      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3081      CArray<int,1>& tmpI = recvDataIIndex[i];   
3082      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3083      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3084          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3085             << "The number of global received index is not coherent with the number of received data index."
3086             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3087
3088      for (ind = 0; ind < tmpI.numElements(); ++ind)
3089      {
3090         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3091         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3092         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3093      } 
3094    }
3095
3096    int nbCompressedData = 0; 
3097    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3098    {
3099       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3100       if ((0 <= indexI) && (0 <= indexJ))
3101         ++nbCompressedData;
3102    }       
3103 
3104    data_i_index.resize(nbCompressedData);
3105    data_j_index.resize(nbCompressedData);
3106
3107    nbCompressedData = 0; 
3108    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3109    {
3110       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3111       if ((0 <= indexI) && (0 <= indexJ))
3112       {
3113          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3114          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3115         ++nbCompressedData;
3116       }
3117    }
3118
3119    // Reset data_ibegin, data_jbegin
3120    data_ibegin.setValue(0);
3121    data_jbegin.setValue(0);
3122  }
3123  CATCH_DUMP_ATTR
3124
3125  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3126  TRY
3127  {
3128    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3129    return transformationMap_.back().second;
3130  }
3131  CATCH_DUMP_ATTR
3132
3133  /*!
3134    Check whether a domain has transformation
3135    \return true if domain has transformation
3136  */
3137  bool CDomain::hasTransformation()
3138  TRY
3139  {
3140    return (!transformationMap_.empty());
3141  }
3142  CATCH_DUMP_ATTR
3143
3144  /*!
3145    Set transformation for current domain. It's the method to move transformation in hierarchy
3146    \param [in] domTrans transformation on domain
3147  */
3148  void CDomain::setTransformations(const TransMapTypes& domTrans)
3149  TRY
3150  {
3151    transformationMap_ = domTrans;
3152  }
3153  CATCH_DUMP_ATTR
3154
3155  /*!
3156    Get all transformation current domain has
3157    \return all transformation
3158  */
3159  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3160  TRY
3161  {
3162    return transformationMap_;
3163  }
3164  CATCH_DUMP_ATTR
3165
3166  void CDomain::duplicateTransformation(CDomain* src)
3167  TRY
3168  {
3169    if (src->hasTransformation())
3170    {
3171      this->setTransformations(src->getAllTransformations());
3172    }
3173  }
3174  CATCH_DUMP_ATTR
3175
3176  /*!
3177   * Go through the hierarchy to find the domain from which the transformations must be inherited
3178   */
3179  void CDomain::solveInheritanceTransformation()
3180  TRY
3181  {
3182    if (hasTransformation() || !hasDirectDomainReference())
3183      return;
3184
3185    CDomain* domain = this;
3186    std::vector<CDomain*> refDomains;
3187    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3188    {
3189      refDomains.push_back(domain);
3190      domain = domain->getDirectDomainReference();
3191    }
3192
3193    if (domain->hasTransformation())
3194      for (size_t i = 0; i < refDomains.size(); ++i)
3195        refDomains[i]->setTransformations(domain->getAllTransformations());
3196  }
3197  CATCH_DUMP_ATTR
3198
3199  void CDomain::setContextClient(CContextClient* contextClient)
3200  TRY
3201  {
3202    if (clientsSet.find(contextClient)==clientsSet.end())
3203    {
3204      clients.push_back(contextClient) ;
3205      clientsSet.insert(contextClient);
3206    }
3207  }
3208  CATCH_DUMP_ATTR
3209
3210  /*!
3211    Parse children nodes of a domain in xml file.
3212    Whenver there is a new transformation, its type and name should be added into this function
3213    \param node child node to process
3214  */
3215  void CDomain::parse(xml::CXMLNode & node)
3216  TRY
3217  {
3218    SuperClass::parse(node);
3219
3220    if (node.goToChildElement())
3221    {
3222      StdString nodeElementName;
3223      do
3224      {
3225        StdString nodeId("");
3226        if (node.getAttributes().end() != node.getAttributes().find("id"))
3227        { nodeId = node.getAttributes()["id"]; }
3228
3229        nodeElementName = node.getElementName();
3230        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3231        it = transformationMapList_.find(nodeElementName);
3232        if (ite != it)
3233        {
3234          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3235                                                                                                                nodeId,
3236                                                                                                                &node)));
3237        }
3238        else
3239        {
3240          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3241                << "The transformation " << nodeElementName << " has not been supported yet.");
3242        }
3243      } while (node.goToNextElement()) ;
3244      node.goToParentElement();
3245    }
3246  }
3247  CATCH_DUMP_ATTR
3248   //----------------------------------------------------------------
3249
3250   DEFINE_REF_FUNC(Domain,domain)
3251
3252   ///---------------------------------------------------------------
3253
3254} // namespace xios
Note: See TracBrowser for help on using the repository browser.