source: XIOS/dev/branch_openmp/src/node/domain.cpp @ 1642

Last change on this file since 1642 was 1642, checked in by yushan, 5 years ago

dev on ADA. add flag switch _usingEP/_usingMPI

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