source: XIOS/dev/dev_olga/src/node/domain.cpp @ 1597

Last change on this file since 1597 was 1597, checked in by oabramkina, 3 years ago

Bug(typo)fix for cases of 1D data with only data_i_index defined and not data_j_index.

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