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

Last change on this file since 1966 was 1966, checked in by oabramkina, 2 years ago

dev_oa: first working version for tiled domains

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