source: XIOS/trunk/src/node/domain.cpp @ 2029

Last change on this file since 2029 was 1972, checked in by yushan, 4 years ago

trunk : debug domain_expand

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