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

Last change on this file since 1054 was 1054, checked in by oabramkina, 7 years ago

dev: intermediate commit.

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