source: XIOS/dev/branch_yushan/src/node/domain.cpp @ 1037

Last change on this file since 1037 was 1037, checked in by yushan, 7 years ago

initialize the branch

  • 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: 78.6 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), nGlobDomain_(), isCompressible_(false), isUnstructed_(false)
34      , isClientAfterTransformationChecked(false), hasLonLat(false)
35      , lonvalue_client(), latvalue_client(), bounds_lon_client(), bounds_lat_client()
36      , isRedistributed_(false), hasPole(false)
37   {
38   }
39
40   CDomain::CDomain(const StdString & id)
41      : CObjectTemplate<CDomain>(id), CDomainAttributes()
42      , isChecked(false), relFiles(), isClientChecked(false), nbConnectedClients_(), indSrv_(), connectedServerRank_()
43      , hasBounds(false), hasArea(false), isDistributed_(false), nGlobDomain_(), isCompressible_(false), isUnstructed_(false)
44      , isClientAfterTransformationChecked(false), hasLonLat(false)
45      , lonvalue_client(), latvalue_client(), bounds_lon_client(), bounds_lat_client()
46      , isRedistributed_(false), hasPole(false)
47   {
48         }
49
50   CDomain::~CDomain(void)
51   {
52   }
53
54   ///---------------------------------------------------------------
55
56   void CDomain::assignMesh(const StdString meshName, const int nvertex)
57   {
58     mesh = CMesh::getMesh(meshName, nvertex);
59   }
60
61   CDomain* CDomain::createDomain()
62   {
63     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
64     return domain;
65   }
66
67   std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
68   bool CDomain::_dummyTransformationMapList = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
69
70   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
71   {
72     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
73     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
74     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
75     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
76     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
77   }
78
79   const std::set<StdString> & CDomain::getRelFiles(void) const
80   {
81      return (this->relFiles);
82   }
83
84
85   const std::vector<int>& CDomain::getIndexesToWrite(void) const
86   {
87     return indexesToWrite;
88   }
89
90   /*!
91     Returns the number of indexes written by each server.
92     \return the number of indexes written by each server
93   */
94   int CDomain::getNumberWrittenIndexes() const
95   {
96     return numberWrittenIndexes_;
97   }
98
99   /*!
100     Returns the total number of indexes written by the servers.
101     \return the total number of indexes written by the servers
102   */
103   int CDomain::getTotalNumberWrittenIndexes() const
104   {
105     return totalNumberWrittenIndexes_;
106   }
107
108   /*!
109     Returns the offset of indexes written by each server.
110     \return the offset of indexes written by each server
111   */
112   int CDomain::getOffsetWrittenIndexes() const
113   {
114     return offsetWrittenIndexes_;
115   }
116
117   //----------------------------------------------------------------
118
119   /*!
120    * Compute the minimum buffer size required to send the attributes to the server(s).
121    *
122    * \return A map associating the server rank with its minimum buffer size.
123    */
124   std::map<int, StdSize> CDomain::getAttributesBufferSize()
125   {
126     CContextClient* client = CContext::getCurrent()->client;
127
128     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes();
129
130     if (client->isServerLeader())
131     {
132       // size estimation for sendServerAttribut
133       size_t size = 11 * sizeof(size_t);
134
135       const std::list<int>& ranks = client->getRanksServerLeader();
136       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
137       {
138         if (size > attributesSizes[*itRank])
139           attributesSizes[*itRank] = size;
140       }
141     }
142
143     std::map<int, std::vector<size_t> >::const_iterator itIndexEnd = indSrv_.end();
144     std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
145     for (size_t k = 0; k < connectedServerRank_.size(); ++k)
146     {
147       int rank = connectedServerRank_[k];
148       std::map<int, std::vector<size_t> >::const_iterator it = indSrv_.find(rank);
149       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
150
151       // size estimation for sendIndex (and sendArea which is always smaller or equal)
152       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
153       if (isCompressible_)
154       {
155         std::map<int, std::vector<int> >::const_iterator itWritten = indWrittenSrv_.find(rank);
156         size_t writtenIdxCount = (itWritten != itWrittenIndexEnd) ? itWritten->second.size() : 0;
157         sizeIndexEvent += CArray<int,1>::size(writtenIdxCount);
158       }
159
160       // size estimation for sendLonLat
161       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
162       if (hasBounds)
163         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
164
165       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
166       if (size > attributesSizes[rank])
167         attributesSizes[rank] = size;
168     }
169
170     return attributesSizes;
171   }
172
173   //----------------------------------------------------------------
174
175   bool CDomain::isEmpty(void) const
176   {
177      return ((this->zoom_ni_srv == 0) ||
178              (this->zoom_nj_srv == 0));
179   }
180
181   //----------------------------------------------------------------
182
183   bool CDomain::IsWritten(const StdString & filename) const
184   {
185      return (this->relFiles.find(filename) != this->relFiles.end());
186   }
187
188   bool CDomain::isWrittenCompressed(const StdString& filename) const
189   {
190      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
191   }
192
193   //----------------------------------------------------------------
194
195   bool CDomain::isDistributed(void) const
196   {
197      return isDistributed_;
198   }
199
200   //----------------------------------------------------------------
201
202   /*!
203    * Test whether the data defined on the domain can be outputted in a compressed way.
204    *
205    * \return true if and only if a mask was defined for this domain
206    */
207   bool CDomain::isCompressible(void) const
208   {
209      return isCompressible_;
210   }
211
212   void CDomain::addRelFile(const StdString & filename)
213   {
214      this->relFiles.insert(filename);
215   }
216
217   void CDomain::addRelFileCompressed(const StdString& filename)
218   {
219      this->relFilesCompressed.insert(filename);
220   }
221
222   StdString CDomain::GetName(void)   { return (StdString("domain")); }
223   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
224   ENodeType CDomain::GetType(void)   { return (eDomain); }
225
226   //----------------------------------------------------------------
227
228   /*!
229     Redistribute RECTILINEAR domain with a number of local domains.
230   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
231   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
232   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
233    \param [in] nbLocalDomain number of local domain on the domain destination
234   */
235   void CDomain::redistribute(int nbLocalDomain)
236   {
237     if (this->isRedistributed_) return;
238
239     this->isRedistributed_ = true;
240     CContext* context = CContext::getCurrent();
241     CContextClient* client = context->client;
242     int rankClient = client->clientRank;
243     int rankOnDomain = rankClient%nbLocalDomain;
244
245     if (ni_glo.isEmpty() || ni_glo <= 0 )
246     {
247        ERROR("CDomain::redistribute(int nbLocalDomain)",
248           << "[ Id = " << this->getId() << " ] "
249           << "The global domain is badly defined,"
250           << " check the \'ni_glo\'  value !")
251     }
252
253     if (nj_glo.isEmpty() || nj_glo <= 0 )
254     {
255        ERROR("CDomain::redistribute(int nbLocalDomain)",
256           << "[ Id = " << this->getId() << " ] "
257           << "The global domain is badly defined,"
258           << " check the \'nj_glo\'  value !")
259     }
260
261     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
262     {
263        int globalDomainSize = ni_glo * nj_glo;
264        if (globalDomainSize <= nbLocalDomain)
265        {
266          for (int idx = 0; idx < nbLocalDomain; ++idx)
267          {
268            if (rankOnDomain < globalDomainSize)
269            {
270              int iIdx = rankOnDomain % ni_glo;
271              int jIdx = rankOnDomain / ni_glo;
272              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
273              ni.setValue(1); nj.setValue(1);
274            }
275            else
276            {
277              ibegin.setValue(0); jbegin.setValue(0);
278              ni.setValue(0); nj.setValue(0);
279            }
280          }
281        }
282        else
283        {
284          float njGlo = nj_glo.getValue();
285          float niGlo = ni_glo.getValue();
286          int nbProcOnX, nbProcOnY, range;
287
288          // Compute (approximately) number of segment on x and y axis
289          float yOverXRatio = njGlo/niGlo;
290
291          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
292          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
293
294          // Simple distribution: Sweep from top to bottom, left to right
295          // Calculate local begin on x
296          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
297          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
298          for (int i = 1; i < nbProcOnX; ++i)
299          {
300            range = ni_glo / nbProcOnX;
301            if (i < (ni_glo%nbProcOnX)) ++range;
302            niVec[i-1] = range;
303            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
304          }
305          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
306
307          // Calculate local begin on y
308          for (int j = 1; j < nbProcOnY; ++j)
309          {
310            range = nj_glo / nbProcOnY;
311            if (j < (nj_glo%nbProcOnY)) ++range;
312            njVec[j-1] = range;
313            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
314          }
315          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
316
317          // Now assign value to ni, ibegin, nj, jbegin
318          int iIdx = rankOnDomain % nbProcOnX;
319          int jIdx = rankOnDomain / nbProcOnX;
320
321          if (rankOnDomain != (nbLocalDomain-1))
322          {
323            ibegin.setValue(ibeginVec[iIdx]);
324            jbegin.setValue(jbeginVec[jIdx]);
325            nj.setValue(njVec[jIdx]);
326            ni.setValue(niVec[iIdx]);
327          }
328          else // just merge all the remaining rectangle into the last one
329          {
330            ibegin.setValue(ibeginVec[iIdx]);
331            jbegin.setValue(jbeginVec[jIdx]);
332            nj.setValue(njVec[jIdx]);
333            ni.setValue(ni_glo - ibeginVec[iIdx]);
334          }
335        }
336
337        // Now fill other attributes
338        if (type_attr::rectilinear == type) fillInRectilinearLonLat();
339     }
340     else  // unstructured domain
341     {
342       if (this->i_index.isEmpty())
343       {
344          int globalDomainSize = ni_glo * nj_glo;
345          if (globalDomainSize <= nbLocalDomain)
346          {
347            for (int idx = 0; idx < nbLocalDomain; ++idx)
348            {
349              if (rankOnDomain < globalDomainSize)
350              {
351                int iIdx = rankOnDomain % ni_glo;
352                int jIdx = rankOnDomain / ni_glo;
353                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
354                ni.setValue(1); nj.setValue(1);
355              }
356              else
357              {
358                ibegin.setValue(0); jbegin.setValue(0);
359                ni.setValue(0); nj.setValue(0);
360              }
361            }
362          }
363          else
364          {
365            float njGlo = nj_glo.getValue();
366            float niGlo = ni_glo.getValue();
367            std::vector<int> ibeginVec(nbLocalDomain,0);
368            std::vector<int> niVec(nbLocalDomain);
369            for (int i = 1; i < nbLocalDomain; ++i)
370            {
371              int range = ni_glo / nbLocalDomain;
372              if (i < (ni_glo%nbLocalDomain)) ++range;
373              niVec[i-1] = range;
374              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
375            }
376            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
377
378            int iIdx = rankOnDomain % nbLocalDomain;
379            ibegin.setValue(ibeginVec[iIdx]);
380            jbegin.setValue(0);
381            ni.setValue(niVec[iIdx]);
382            nj.setValue(1);
383          }
384        }
385        else
386        {
387          ibegin.setValue(this->i_index(0));
388          jbegin.setValue(0);
389          ni.setValue(this->i_index.numElements());
390          nj.setValue(1);
391        }
392     }
393
394     checkDomain();
395   }
396
397   /*!
398     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
399     Range of longitude value from 0 - 360
400     Range of latitude value from -90 - +90
401   */
402   void CDomain::fillInRectilinearLonLat()
403   {
404     if (!lonvalue_rectilinear_read_from_file.isEmpty())
405     {
406       lonvalue_1d.resize(ni);
407       for (int idx = 0; idx < ni; ++idx)
408         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
409       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
410       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
411     }
412     else
413     {
414       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
415       lonvalue_1d.resize(ni);
416       double lonRange = lon_end - lon_start;
417       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
418
419        // Assign lon value
420       for (int i = 0; i < ni; ++i)
421       {
422         if (0 == (ibegin + i))
423         {
424           lonvalue_1d(i) = lon_start;
425         }
426         else if (ni_glo == (ibegin + i + 1))
427         {
428           lonvalue_1d(i) = lon_end;
429         }
430         else
431         {
432           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
433         }
434       }
435     }
436
437
438     if (!latvalue_rectilinear_read_from_file.isEmpty())
439     {
440       latvalue_1d.resize(nj);
441       for (int idx = 0; idx < nj; ++idx)
442         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
443       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
444       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
445     }
446     else
447     {
448       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
449       latvalue_1d.resize(nj);
450
451       double latRange = lat_end - lat_start;
452       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
453
454       for (int j = 0; j < nj; ++j)
455       {
456         if (0 == (jbegin + j))
457         {
458            latvalue_1d(j) = lat_start;
459         }
460         else if (nj_glo == (jbegin + j + 1))
461         {
462            latvalue_1d(j) = lat_end;
463         }
464         else
465         {
466           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
467         }
468       }
469     }
470   }
471
472
473
474   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g)
475   {
476          CContext* context = CContext::getCurrent();
477      CContextClient* client = context->client;
478          lon_g.resize(ni_glo) ;
479          lat_g.resize(nj_glo) ;
480
481
482          int* ibegin_g = new int[client->clientSize] ;
483          int* jbegin_g = new int[client->clientSize] ;
484          int* ni_g = new int[client->clientSize] ;
485          int* nj_g = new int[client->clientSize] ;
486          int v ;
487          v=ibegin ;
488          MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,client->intraComm) ;
489          v=jbegin ;
490          MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,client->intraComm) ;
491          v=ni ;
492          MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,client->intraComm) ;
493          v=nj ;
494          MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,client->intraComm) ;
495
496          MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,client->intraComm) ;
497          MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,client->intraComm) ;
498
499      delete[] ibegin_g ;
500      delete[] jbegin_g ;
501      delete[] ni_g ;
502      delete[] nj_g ;
503   }
504
505   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
506                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
507   {
508     int i,j,k;
509
510     const int nvertexValue = 4;
511     boundsLon.resize(nvertexValue,ni*nj);
512
513     if (ni_glo>1)
514     {
515       double lonStepStart = lon(1)-lon(0);
516       bounds_lon_start=lon(0) - lonStepStart/2;
517       double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2);
518       bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
519       double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
520
521       // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
522       if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
523       {
524         bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
525         bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
526       }
527     }
528     else
529     {
530       if (bounds_lon_start.isEmpty()) bounds_lon_start=-180. ;
531       if (bounds_lon_end.isEmpty()) bounds_lon_end=180.-1e-8 ;
532     }
533
534     for(j=0;j<nj;++j)
535       for(i=0;i<ni;++i)
536       {
537         k=j*ni+i;
538         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start
539                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2;
540         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end
541                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2;
542       }
543
544
545    boundsLat.resize(nvertexValue,nj*ni);
546    bool isNorthPole=false ;
547    bool isSouthPole=false ;
548    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
549    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
550
551    // lat boundaries beyond pole the assimilate it to pole
552    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
553    if (nj_glo>1)
554    {
555      double latStepStart = lat(1)-lat(0);
556      if (isNorthPole) bounds_lat_start=lat(0);
557      else
558      {
559        bounds_lat_start=lat(0)-latStepStart/2;
560        if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
561        else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
562        else if (bounds_lat_start <= 90 && bounds_lat_start >= lat(0))
563        {
564          if ( std::abs(90-bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
565        }
566        else if (bounds_lat_start >= -90 && bounds_lat_start <= lat(0))
567        {
568          if ( std::abs(-90 - bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
569        }
570      }
571
572      double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
573      if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
574      else
575      {
576        bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
577
578        if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
579        else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
580        else if (bounds_lat_end <= 90 && bounds_lat_end >= lat(nj_glo-1))
581        {
582          if ( std::abs(90-bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
583        }
584        else if (bounds_lat_end >= -90 && bounds_lat_end <= lat(nj_glo-1))
585        {
586          if ( std::abs(-90 - bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
587        }
588      }
589    }
590    else
591    {
592      if (bounds_lat_start.isEmpty()) bounds_lon_start=-90. ;
593      if (bounds_lat_end.isEmpty()) bounds_lat_end=90 ;
594    }
595
596    for(j=0;j<nj;++j)
597      for(i=0;i<ni;++i)
598      {
599        k=j*ni+i;
600        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start
601                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2;
602        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end
603                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2;
604      }
605   }
606
607   void CDomain::checkDomain(void)
608   {
609     if (type.isEmpty())
610     {
611       ERROR("CDomain::checkDomain(void)",
612             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
613             << "The domain type is mandatory, "
614             << "please define the 'type' attribute.")
615     }
616
617     if (type == type_attr::gaussian) 
618     {
619           hasPole=true ;
620             type.setValue(type_attr::unstructured) ;
621           }
622           else if (type == type_attr::rectilinear) hasPole=true ;
623         
624     if (type == type_attr::unstructured)
625     {
626        if (ni_glo.isEmpty())
627        {
628          ERROR("CDomain::checkDomain(void)",
629                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
630                << "The global domain is badly defined, "
631                << "the mandatory 'ni_glo' attribute is missing.")
632        }
633        else if (ni_glo <= 0)
634        {
635          ERROR("CDomain::checkDomain(void)",
636                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
637                << "The global domain is badly defined, "
638                << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
639        }
640        isUnstructed_ = true;
641        nj_glo = 1;
642        nj = 1;
643        jbegin = 0;
644        if (!i_index.isEmpty()) ni = i_index.numElements();
645        j_index.resize(ni);
646        for(int i=0;i<ni;++i) j_index(i)=0;
647
648        if (!area.isEmpty())
649          area.transposeSelf(1, 0);
650     }
651
652     if (ni_glo.isEmpty())
653     {
654       ERROR("CDomain::checkDomain(void)",
655             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
656             << "The global domain is badly defined, "
657             << "the mandatory 'ni_glo' attribute is missing.")
658     }
659     else if (ni_glo <= 0)
660     {
661       ERROR("CDomain::checkDomain(void)",
662             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
663             << "The global domain is badly defined, "
664             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
665     }
666
667     if (nj_glo.isEmpty())
668     {
669       ERROR("CDomain::checkDomain(void)",
670             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
671             << "The global domain is badly defined, "
672             << "the mandatory 'nj_glo' attribute is missing.")
673     }
674     else if (nj_glo <= 0)
675     {
676       ERROR("CDomain::checkDomain(void)",
677             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
678             << "The global domain is badly defined, "
679             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
680     }
681
682     checkLocalIDomain();
683     checkLocalJDomain();
684
685     if (i_index.isEmpty())
686     {
687       i_index.resize(ni*nj);
688       for (int j = 0; j < nj; ++j)
689         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
690     }
691
692     if (j_index.isEmpty())
693     {
694       j_index.resize(ni*nj);
695       for (int j = 0; j < nj; ++j)
696         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
697     }
698     computeNGlobDomain();
699     checkZoom();
700
701     isDistributed_ = !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
702                        (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
703   }
704
705   // Check global zoom of a domain
706   // If there is no zoom defined for the domain, zoom will have value of global doamin
707   void CDomain::checkZoom(void)
708   {
709     if (global_zoom_ibegin.isEmpty())
710      global_zoom_ibegin.setValue(0);
711     if (global_zoom_ni.isEmpty())
712      global_zoom_ni.setValue(ni_glo);
713     if (global_zoom_jbegin.isEmpty())
714      global_zoom_jbegin.setValue(0);
715     if (global_zoom_nj.isEmpty())
716      global_zoom_nj.setValue(nj_glo);
717   }
718
719   //----------------------------------------------------------------
720
721   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
722   void CDomain::checkLocalIDomain(void)
723   {
724      // If ibegin and ni are provided then we use them to check the validity of local domain
725      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
726      {
727        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
728        {
729          ERROR("CDomain::checkLocalIDomain(void)",
730                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
731                << "The local domain is wrongly defined,"
732                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
733        }
734      }
735
736      // i_index has higher priority than ibegin and ni
737      if (!i_index.isEmpty())
738      {
739        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
740        if (ni.isEmpty()) 
741        {         
742         // No information about ni
743          int minIndex = ni_glo - 1;
744          int maxIndex = 0;
745          for (int idx = 0; idx < i_index.numElements(); ++idx)
746          {
747            if (i_index(idx) < minIndex) minIndex = i_index(idx);
748            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
749          }
750          ni = maxIndex - minIndex + 1; 
751          minIIndex = minIIndex;         
752        }
753
754        // It's not so correct but if ibegin is not the first value of i_index
755        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
756        if (ibegin.isEmpty()) ibegin = minIIndex;
757      }
758      else if (ibegin.isEmpty() && ni.isEmpty())
759      {
760        ibegin = 0;
761        ni = ni_glo;
762      }
763      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
764      {
765        ERROR("CDomain::checkLocalIDomain(void)",
766              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
767              << "The local domain is wrongly defined," << endl
768              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
769              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
770      }
771       
772
773      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
774      {
775        ERROR("CDomain::checkLocalIDomain(void)",
776              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
777              << "The local domain is wrongly defined,"
778              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
779      }
780   }
781
782   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
783   void CDomain::checkLocalJDomain(void)
784   {
785    // If jbegin and nj are provided then we use them to check the validity of local domain
786     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
787     {
788       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
789       {
790         ERROR("CDomain::checkLocalJDomain(void)",
791                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
792                << "The local domain is wrongly defined,"
793                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
794       }
795     }
796
797     if (!j_index.isEmpty())
798     {
799        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
800        if (nj.isEmpty()) 
801        {
802          // No information about nj
803          int minIndex = nj_glo - 1;
804          int maxIndex = 0;
805          for (int idx = 0; idx < j_index.numElements(); ++idx)
806          {
807            if (j_index(idx) < minIndex) minIndex = j_index(idx);
808            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
809          }
810          nj = maxIndex - minIndex + 1;
811          minJIndex = minIndex; 
812        } 
813        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
814        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
815       if (jbegin.isEmpty()) jbegin = minJIndex;       
816     }
817     else if (jbegin.isEmpty() && nj.isEmpty())
818     {
819       jbegin = 0;
820       nj = nj_glo;
821     }     
822
823
824     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
825     {
826       ERROR("CDomain::checkLocalJDomain(void)",
827              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
828              << "The local domain is wrongly defined,"
829              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
830     }
831   }
832
833   //----------------------------------------------------------------
834
835   void CDomain::checkMask(void)
836   {
837      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
838        ERROR("CDomain::checkMask(void)",
839              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
840              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
841              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
842
843      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
844      {
845        if (mask_1d.numElements() != i_index.numElements())
846          ERROR("CDomain::checkMask(void)",
847                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
848                << "'mask_1d' does not have the same size as the local domain." << std::endl
849                << "Local size is " << i_index.numElements() << "." << std::endl
850                << "Mask size is " << mask_1d.numElements() << ".");
851      }
852
853      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
854      {
855        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
856          ERROR("CDomain::checkMask(void)",
857                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
858                << "The mask does not have the same size as the local domain." << std::endl
859                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
860                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
861      }
862
863      if (!mask_2d.isEmpty())
864      {
865        mask_1d.resize(mask_2d.extent(0) * mask_2d.extent(1));
866        for (int j = 0; j < nj; ++j)
867          for (int i = 0; i < ni; ++i) mask_1d(i+j*ni) = mask_2d(i,j);
868        mask_2d.reset();
869      }
870      else if (mask_1d.isEmpty())
871      {
872        mask_1d.resize(i_index.numElements());
873        for (int i = 0; i < i_index.numElements(); ++i) mask_1d(i) = true;
874      }
875   }
876
877   //----------------------------------------------------------------
878
879   void CDomain::checkDomainData(void)
880   {
881      if (data_dim.isEmpty())
882      {
883        data_dim.setValue(1);
884      }
885      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
886      {
887        ERROR("CDomain::checkDomainData(void)",
888              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
889              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
890      }
891
892      if (data_ibegin.isEmpty())
893         data_ibegin.setValue(0);
894      if (data_jbegin.isEmpty())
895         data_jbegin.setValue(0);
896
897      if (data_ni.isEmpty())
898      {
899        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
900      }
901      else if (data_ni.getValue() < 0)
902      {
903        ERROR("CDomain::checkDomainData(void)",
904              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
905              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
906      }
907
908      if (data_nj.isEmpty())
909      {
910        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
911      }
912      else if (data_nj.getValue() < 0)
913      {
914        ERROR("CDomain::checkDomainData(void)",
915              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
916              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
917      }
918   }
919
920   //----------------------------------------------------------------
921
922   void CDomain::checkCompression(void)
923   {
924      if (!data_i_index.isEmpty())
925      {
926        if (!data_j_index.isEmpty() &&
927            data_j_index.numElements() != data_i_index.numElements())
928        {
929           ERROR("CDomain::checkCompression(void)",
930                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
931                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
932                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
933                 << "'data_j_index' size = " << data_j_index.numElements());
934        }
935
936        if (2 == data_dim)
937        {
938          if (data_j_index.isEmpty())
939          {
940             ERROR("CDomain::checkCompression(void)",
941                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
942                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
943          }
944        }
945        else // (1 == data_dim)
946        {
947          if (data_j_index.isEmpty())
948          {
949            data_j_index.resize(data_ni);
950            for (int j = 0; j < data_ni; ++j) data_j_index(j) = 0;
951          }
952        }
953      }
954      else
955      {
956        if (data_dim == 2 && !data_j_index.isEmpty())
957          ERROR("CDomain::checkCompression(void)",
958                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
959                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
960
961        if (1 == data_dim)
962        {
963          data_i_index.resize(data_ni);
964          data_j_index.resize(data_ni);
965
966          for (int i = 0; i < data_ni; ++i)
967          {
968            data_i_index(i) = i;
969            data_j_index(i) = 0;
970          }
971        }
972        else // (data_dim == 2)
973        {
974          const int dsize = data_ni * data_nj;
975          data_i_index.resize(dsize);
976          data_j_index.resize(dsize);
977
978          for(int count = 0, j = 0; j < data_nj; ++j)
979          {
980            for(int i = 0; i < data_ni; ++i, ++count)
981            {
982              data_i_index(count) = i;
983              data_j_index(count) = j;
984            }
985          }
986        }
987      }
988   }
989
990   //----------------------------------------------------------------
991   void CDomain::computeLocalMask(void)
992   {
993     localMask.resize(ni*nj) ;
994     localMask=false ;
995     size_t zoom_ibegin=global_zoom_ibegin ;
996     size_t zoom_iend=global_zoom_ibegin+global_zoom_ni-1 ;
997     size_t zoom_jbegin=global_zoom_jbegin ;
998     size_t zoom_jend=global_zoom_jbegin+global_zoom_nj-1 ;
999
1000
1001     size_t dn=data_i_index.numElements() ;
1002     int i,j ;
1003     size_t k,ind ;
1004
1005     for(k=0;k<dn;k++)
1006     {
1007       if (data_dim==2)
1008       {
1009          i=data_i_index(k)+data_ibegin ;
1010          j=data_j_index(k)+data_jbegin ;
1011       }
1012       else
1013       {
1014          i=(data_i_index(k)+data_ibegin)%ni ;
1015          j=(data_i_index(k)+data_ibegin)/ni ;
1016       }
1017
1018       if (i>=0 && i<ni && j>=0 && j<nj)
1019         if (i+ibegin>=zoom_ibegin && i+ibegin<=zoom_iend && j+jbegin>=zoom_jbegin && j+jbegin<=zoom_jend)
1020         {
1021           ind=i+ni*j ;
1022           localMask(ind)=mask_1d(ind) ;
1023         }
1024     }
1025   }
1026
1027
1028
1029
1030
1031
1032
1033   void CDomain::checkEligibilityForCompressedOutput(void)
1034   {
1035     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1036     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1037   }
1038
1039   //----------------------------------------------------------------
1040
1041   void CDomain::completeLonLatClient(void)
1042   {
1043     if (!lonvalue_2d.isEmpty())
1044     {
1045       lonvalue_client.resize(ni * nj);
1046       latvalue_client.resize(ni * nj);
1047       if (hasBounds)
1048       {
1049         bounds_lon_client.resize(nvertex, ni * nj);
1050         bounds_lat_client.resize(nvertex, ni * nj);
1051       }
1052
1053       for (int j = 0; j < nj; ++j)
1054       {
1055         for (int i = 0; i < ni; ++i)
1056         {
1057           int k = j * ni + i;
1058
1059           lonvalue_client(k) = lonvalue_2d(i,j);
1060           latvalue_client(k) = latvalue_2d(i,j);
1061
1062           if (hasBounds)
1063           {
1064             for (int n = 0; n < nvertex; ++n)
1065             {
1066               bounds_lon_client(n,k) = bounds_lon_2d(n,i,j);
1067               bounds_lat_client(n,k) = bounds_lat_2d(n,i,j);
1068             }
1069           }
1070         }
1071       }
1072     }
1073     else if (!lonvalue_1d.isEmpty())
1074     {
1075       if (type_attr::rectilinear == type)
1076       {
1077         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1078         {
1079           lonvalue_client.resize(ni * nj);
1080           latvalue_client.resize(ni * nj);
1081           if (hasBounds)
1082           {
1083             bounds_lon_client.resize(nvertex, ni * nj);
1084             bounds_lat_client.resize(nvertex, ni * nj);
1085           }
1086
1087           for (int j = 0; j < nj; ++j)
1088           {
1089             for (int i = 0; i < ni; ++i)
1090             {
1091               int k = j * ni + i;
1092
1093               lonvalue_client(k) = lonvalue_1d(i);
1094               latvalue_client(k) = latvalue_1d(j);
1095
1096               if (hasBounds)
1097               {
1098                 for (int n = 0; n < nvertex; ++n)
1099                 {
1100                   bounds_lon_client(n,k) = bounds_lon_1d(n,i);
1101                   bounds_lat_client(n,k) = bounds_lat_1d(n,j);
1102                 }
1103               }
1104             }
1105           }
1106         }
1107         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements())
1108         {
1109           lonvalue_client.reference(lonvalue_1d);
1110           latvalue_client.reference(latvalue_1d);
1111            if (hasBounds)
1112           {
1113             bounds_lon_client.reference(bounds_lon_1d);
1114             bounds_lat_client.reference(bounds_lat_1d);
1115           }
1116         }
1117         else
1118           ERROR("CDomain::completeLonClient(void)",
1119                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1120                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1121                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1122                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1123                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1124                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1125       }
1126       else if (type == type_attr::curvilinear || type == type_attr::unstructured)
1127       {
1128         lonvalue_client.reference(lonvalue_1d);
1129         latvalue_client.reference(latvalue_1d);
1130         if (hasBounds)
1131         {
1132           bounds_lon_client.reference(bounds_lon_1d);
1133           bounds_lat_client.reference(bounds_lat_1d);
1134         }
1135       }
1136     }
1137   }
1138
1139   void CDomain::checkBounds(void)
1140   {
1141     if (!nvertex.isEmpty() && nvertex > 0)
1142     {
1143       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1144         ERROR("CDomain::checkBounds(void)",
1145               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1146               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1147               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1148
1149       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1150         ERROR("CDomain::checkBounds(void)",
1151               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1152               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1153               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1154
1155       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1156       {
1157         ERROR("CDomain::checkBounds(void)",
1158               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1159               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1160               << "Please define either both attributes or none.");
1161       }
1162
1163       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1164       {
1165         ERROR("CDomain::checkBounds(void)",
1166               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1167               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1168               << "Please define either both attributes or none.");
1169       }
1170
1171       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1172         ERROR("CDomain::checkBounds(void)",
1173               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1174               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1175               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(1)
1176               << " but nvertex is " << nvertex.getValue() << ".");
1177
1178       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1179         ERROR("CDomain::checkBounds(void)",
1180               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1181               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1182               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(2)
1183               << " but nvertex is " << nvertex.getValue() << ".");
1184
1185       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1186         ERROR("CDomain::checkBounds(void)",
1187               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1188               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1189
1190       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1191         ERROR("CDomain::checkBounds(void)",
1192               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1193               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1194
1195       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1196         ERROR("CDomain::checkBounds(void)",
1197               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1198               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1199               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(1)
1200               << " but nvertex is " << nvertex.getValue() << ".");
1201
1202       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1203         ERROR("CDomain::checkBounds(void)",
1204               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1205               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1206               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(2)
1207               << " but nvertex is " << nvertex.getValue() << ".");
1208
1209       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1210         ERROR("CDomain::checkBounds(void)",
1211               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1212               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1213
1214       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1215         ERROR("CDomain::checkBounds(void)",
1216               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1217
1218       hasBounds = true;
1219     }
1220     else
1221     {
1222       hasBounds = false;
1223       nvertex = 0;
1224     }
1225   }
1226
1227   void CDomain::checkArea(void)
1228   {
1229     hasArea = !area.isEmpty();
1230     if (hasArea)
1231     {
1232       if (area.extent(0) != ni || area.extent(1) != nj)
1233       {
1234         ERROR("CDomain::checkArea(void)",
1235               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1236               << "The area does not have the same size as the local domain." << std::endl
1237               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1238               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1239       }
1240     }
1241   }
1242
1243   void CDomain::checkLonLat()
1244   {
1245     hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1246                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1247     if (hasLonLat)
1248     {
1249       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1250         ERROR("CDomain::checkLonLat()",
1251               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1252               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1253               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1254
1255       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1256       {
1257         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1258           ERROR("CDomain::checkLonLat()",
1259                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1260                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1261                 << "Local size is " << i_index.numElements() << "." << std::endl
1262                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1263       }
1264
1265       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1266       {
1267         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1268           ERROR("CDomain::checkLonLat()",
1269                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1270                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1271                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1272                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1273       }
1274
1275       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1276         ERROR("CDomain::checkLonLat()",
1277               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1278               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1279               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1280
1281       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1282       {
1283         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1284           ERROR("CDomain::checkLonLat()",
1285                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1286                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1287                 << "Local size is " << i_index.numElements() << "." << std::endl
1288                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1289       }
1290
1291       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1292       {
1293         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1294           ERROR("CDomain::checkLonLat()",
1295                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1296                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1297                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1298                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1299       }
1300     }
1301   }
1302
1303   void CDomain::checkAttributesOnClientAfterTransformation()
1304   {
1305     int myRank;
1306     MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
1307
1308     CContext* context=CContext::getCurrent() ;  //printf("myRank = %d, CContext::getCurrent OK\n", myRank);
1309
1310     //printf("myRank = %d, this->isClientAfterTransformationChecked = %d\n", myRank, this->isClientAfterTransformationChecked);
1311     if (this->isClientAfterTransformationChecked) return;
1312     //printf("myRank = %d, context->hasClient = %d\n", myRank, context->hasClient);
1313     if (context->hasClient)
1314     {
1315       this->checkMask();  //printf("myRank = %d, this->checkMask OK\n", myRank);
1316       //printf("myRank = %d, hasLonLat = %d, hasArea = %d, isCompressible_ = %d\n", myRank, hasLonLat, hasArea, isCompressible_);
1317       if (hasLonLat || hasArea || isCompressible_) 
1318        {
1319          //printf("myRank = %d, start this->computeConnectedServer\n", myRank);
1320          this->computeConnectedServer();
1321          //printf("myRank = %d, this->computeConnectedServer OK\n", myRank);
1322        }
1323       //printf("myRank = %d, hasLonLat = %d\n", myRank, hasLonLat);
1324       if (hasLonLat) 
1325        {
1326          this->completeLonLatClient();
1327          //printf("myRank = %d, this->completeLonLatClient OK\n", myRank);
1328        }
1329     }
1330
1331     this->isClientAfterTransformationChecked = true;
1332   }
1333
1334   //----------------------------------------------------------------
1335   // Divide function checkAttributes into 2 seperate ones
1336   // This function only checks all attributes of current domain
1337   void CDomain::checkAttributesOnClient()
1338   {
1339     if (this->isClientChecked) return;
1340     CContext* context=CContext::getCurrent();
1341
1342      this->checkDomain();
1343      this->checkBounds();
1344      this->checkArea();
1345      this->checkLonLat();
1346
1347      if (context->hasClient)
1348      { // CÃŽté client uniquement
1349         this->checkMask();
1350         this->checkDomainData();
1351         this->checkCompression();
1352         this->computeLocalMask() ;
1353      }
1354      else
1355      { // CÃŽté serveur uniquement
1356      }
1357
1358      this->isClientChecked = true;
1359   }
1360
1361   // Send all checked attributes to server
1362   void CDomain::sendCheckedAttributes()
1363   {
1364     if (!this->isClientChecked) checkAttributesOnClient();
1365     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1366     CContext* context=CContext::getCurrent() ;
1367
1368     if (this->isChecked) return;
1369     if (context->hasClient)
1370     {
1371       sendServerAttribut();
1372       if (hasLonLat || hasArea || isCompressible_) sendLonLatArea();
1373     }
1374     this->isChecked = true;
1375   }
1376
1377   void CDomain::checkAttributes(void)
1378   {
1379      if (this->isChecked) return;
1380      CContext* context=CContext::getCurrent() ;
1381
1382      this->checkDomain();
1383      this->checkLonLat();
1384      this->checkBounds();
1385      this->checkArea();
1386
1387      if (context->hasClient)
1388      { // CÃŽté client uniquement
1389         this->checkMask();
1390         this->checkDomainData();
1391         this->checkCompression();
1392         this->computeLocalMask() ;
1393
1394      }
1395      else
1396      { // CÃŽté serveur uniquement
1397      }
1398
1399      if (context->hasClient)
1400      {
1401        this->computeConnectedServer();
1402        this->completeLonLatClient();
1403        this->sendServerAttribut();
1404        this->sendLonLatArea();
1405      }
1406
1407      this->isChecked = true;
1408   }
1409
1410  void CDomain::sendServerAttribut(void)
1411  {
1412    CContext* context = CContext::getCurrent();
1413    CContextClient* client = context->client;
1414    int nbServer = client->serverSize;
1415
1416    CServerDistributionDescription serverDescription(nGlobDomain_, nbServer);
1417    if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
1418    else serverDescription.computeServerDistribution(false, 1);
1419
1420    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
1421    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
1422
1423    CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
1424    if (client->isServerLeader())
1425    {
1426      std::list<CMessage> msgs;
1427
1428      const std::list<int>& ranks = client->getRanksServerLeader();
1429      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
1430      {
1431        // Use const int to ensure CMessage holds a copy of the value instead of just a reference
1432        const int ibegin_srv = serverIndexBegin[*itRank][0];
1433        const int jbegin_srv = serverIndexBegin[*itRank][1];
1434        const int ni_srv = serverDimensionSizes[*itRank][0];
1435        const int nj_srv = serverDimensionSizes[*itRank][1];
1436        const int iend_srv = ibegin_srv + ni_srv - 1;
1437        const int jend_srv = jbegin_srv + nj_srv - 1;
1438
1439        msgs.push_back(CMessage());
1440        CMessage& msg = msgs.back();
1441        msg << this->getId() ;
1442        msg << ni_srv << ibegin_srv << iend_srv << nj_srv << jbegin_srv << jend_srv;
1443        msg << global_zoom_ni.getValue() << global_zoom_ibegin.getValue() << global_zoom_nj.getValue() << global_zoom_jbegin.getValue();
1444        msg << isCompressible_;
1445
1446        event.push(*itRank,1,msg);
1447      }
1448      client->sendEvent(event);
1449    }
1450    else client->sendEvent(event);
1451  }
1452
1453  void CDomain::computeNGlobDomain()
1454  {
1455    nGlobDomain_.resize(2);
1456    nGlobDomain_[0] = ni_glo.getValue();
1457    nGlobDomain_[1] = nj_glo.getValue();
1458  }
1459
1460  void CDomain::computeConnectedServer(void)
1461  {
1462    int myRank;
1463    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
1464
1465    CContext* context=CContext::getCurrent() ;
1466    CContextClient* client=context->client ;
1467    int nbServer=client->serverSize;
1468    int rank = client->clientRank;
1469    bool doComputeGlobalIndexServer = true;
1470
1471    int i,j,i_ind,j_ind, nbIndex;
1472    int global_zoom_iend=global_zoom_ibegin+global_zoom_ni-1 ;
1473    int global_zoom_jend=global_zoom_jbegin+global_zoom_nj-1 ;
1474
1475    // Precompute number of index
1476    int globalIndexCountZoom = 0;
1477    nbIndex = i_index.numElements();
1478    for (i = 0; i < nbIndex; ++i)
1479    {
1480      i_ind=i_index(i);
1481      j_ind=j_index(i);
1482
1483      if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1484      {
1485        ++globalIndexCountZoom;
1486      }
1487    }
1488
1489    int globalIndexWrittenCount = 0;
1490    if (isCompressible_)
1491    {
1492      for (i = 0; i < data_i_index.numElements(); ++i)
1493      {
1494        i_ind = CDistributionClient::getDomainIndex(data_i_index(i), data_j_index(i),
1495                                                    data_ibegin, data_jbegin, data_dim, ni,
1496                                                    j_ind);
1497        if (i_ind >= 0 && i_ind < ni && j_ind >= 0 && j_ind < nj && mask_1d(i_ind + j_ind * ni))
1498        {
1499          i_ind += ibegin;
1500          j_ind += jbegin;
1501          if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1502            ++globalIndexWrittenCount;
1503        }
1504      }
1505    }
1506
1507    // Fill in index
1508    CArray<size_t,1> globalIndexDomainZoom(globalIndexCountZoom);
1509    CArray<size_t,1> localIndexDomainZoom(globalIndexCountZoom);
1510    CArray<size_t,1> globalIndexDomain(nbIndex);
1511    size_t globalIndex;
1512    int globalIndexCount = 0;
1513    globalIndexCountZoom = 0;
1514
1515    for (i = 0; i < nbIndex; ++i)
1516    {
1517      i_ind=i_index(i);
1518      j_ind=j_index(i);
1519      globalIndex = i_ind + j_ind * ni_glo;
1520      globalIndexDomain(globalIndexCount) = globalIndex;
1521      ++globalIndexCount;
1522      if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1523      {
1524        globalIndexDomainZoom(globalIndexCountZoom) = globalIndex;
1525        localIndexDomainZoom(globalIndexCountZoom) = i;
1526        ++globalIndexCountZoom;
1527      }
1528    }
1529
1530    CArray<int,1> globalIndexWrittenDomain(globalIndexWrittenCount);
1531    if (isCompressible_)
1532    {
1533      globalIndexWrittenCount = 0;
1534      for (i = 0; i < data_i_index.numElements(); ++i)
1535      {
1536        i_ind = CDistributionClient::getDomainIndex(data_i_index(i), data_j_index(i),
1537                                                    data_ibegin, data_jbegin, data_dim, ni,
1538                                                    j_ind);
1539        if (i_ind >= 0 && i_ind < ni && j_ind >= 0 && j_ind < nj && mask_1d(i_ind + j_ind * ni))
1540        {
1541          i_ind += ibegin;
1542          j_ind += jbegin;
1543          if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1544          {
1545            globalIndexWrittenDomain(globalIndexWrittenCount) = i_ind + j_ind * ni_glo;
1546            ++globalIndexWrittenCount;
1547          }
1548        }
1549      }
1550    }
1551   
1552
1553    size_t globalSizeIndex = 1, indexBegin, indexEnd;
1554    int range, clientSize = client->clientSize;
1555    for (int i = 0; i < nGlobDomain_.size(); ++i) globalSizeIndex *= nGlobDomain_[i];
1556    indexBegin = 0;
1557    if (globalSizeIndex <= clientSize)
1558    {
1559      indexBegin = rank%globalSizeIndex;
1560      indexEnd = indexBegin;
1561    }
1562    else
1563    {
1564      for (int i = 0; i < clientSize; ++i)
1565      {
1566        range = globalSizeIndex / clientSize;
1567        if (i < (globalSizeIndex%clientSize)) ++range;
1568        if (i == client->clientRank) break;
1569        indexBegin += range;
1570      }
1571      indexEnd = indexBegin + range - 1;
1572    }
1573   
1574    CServerDistributionDescription serverDescription(nGlobDomain_, nbServer);
1575    if (isUnstructed_) serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 0);
1576    else serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 1);
1577
1578    //printf("myRank = %d, check 7 OK\n", myRank);
1579
1580    CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(),
1581                                                                                client->intraComm);
1582    //printf("myRank = %d new OK\n", myRank);
1583
1584    clientServerMap->computeServerIndexMapping(globalIndexDomain); 
1585    //printf("myRank = %d, clientServerMap->computeServerIndexMapping(globalIndexDomain) OK\n", myRank);
1586   
1587    const CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
1588    //printf("myRank = %d, clientServerMap->getGlobalIndexOnServer OK\n", myRank);
1589   
1590    //printf("myRank = %d, check 8 OK\n", myRank);
1591
1592    CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
1593                                                         ite = globalIndexDomainOnServer.end();
1594    typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
1595    std::vector<int>::iterator itVec;
1596
1597    indSrv_.clear();
1598    indWrittenSrv_.clear();
1599    for (; it != ite; ++it)
1600    {
1601      int rank = it->first;
1602      int indexSize = it->second.size();
1603      std::vector<int> permutIndex(indexSize);
1604      XIOSAlgorithms::fillInIndex(indexSize, permutIndex);
1605      XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(it->second, permutIndex);
1606      BinarySearch binSearch(it->second);
1607      int nb = globalIndexDomainZoom.numElements();
1608      for (int i = 0; i < nb; ++i)
1609      {
1610        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndexDomainZoom(i), itVec))
1611        {
1612          indSrv_[rank].push_back(localIndexDomainZoom(i));
1613        }
1614      }
1615      for (int i = 0; i < globalIndexWrittenDomain.numElements(); ++i)
1616      {
1617        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndexWrittenDomain(i), itVec))
1618        {
1619          indWrittenSrv_[rank].push_back(globalIndexWrittenDomain(i));
1620        }
1621      }
1622    }
1623
1624    //printf("myRank = %d, check 9 OK\n", myRank);
1625
1626    connectedServerRank_.clear();
1627    for (it = globalIndexDomainOnServer.begin(); it != ite; ++it) {
1628      connectedServerRank_.push_back(it->first);
1629    }
1630
1631    nbConnectedClients_ = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_);
1632
1633    delete clientServerMap;
1634    //printf("myRank = %d, check 10 OK\n", myRank);
1635  }
1636
1637  const std::map<int, vector<size_t> >& CDomain::getIndexServer() const
1638  {
1639    return indSrv_;
1640  }
1641
1642  /*!
1643    Send index from client to server(s)
1644  */
1645  void CDomain::sendIndex()
1646  {
1647    int ns, n, i, j, ind, nv, idx;
1648    CContext* context = CContext::getCurrent();
1649    CContextClient* client=context->client;
1650
1651    CEventClient eventIndex(getType(), EVENT_ID_INDEX);
1652
1653    list<CMessage> list_msgsIndex;
1654    list<CArray<int,1> > list_indi, list_indj, list_writtenInd;
1655
1656    std::map<int, std::vector<size_t> >::const_iterator it, iteMap;
1657    iteMap = indSrv_.end();
1658    for (int k = 0; k < connectedServerRank_.size(); ++k)
1659    {
1660      int nbData = 0;
1661      int rank = connectedServerRank_[k];
1662      it = indSrv_.find(rank);
1663      if (iteMap != it)
1664        nbData = it->second.size();
1665
1666      list_indi.push_back(CArray<int,1>(nbData));
1667      list_indj.push_back(CArray<int,1>(nbData));
1668
1669      CArray<int,1>& indi = list_indi.back();
1670      CArray<int,1>& indj = list_indj.back();
1671      const std::vector<size_t>& temp = it->second;
1672      for (n = 0; n < nbData; ++n)
1673      {
1674        idx = static_cast<int>(it->second[n]);
1675        indi(n) = i_index(idx);
1676        indj(n) = j_index(idx);
1677      }
1678
1679      list_msgsIndex.push_back(CMessage());
1680
1681      list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
1682      list_msgsIndex.back() << isCurvilinear;
1683      list_msgsIndex.back() << list_indi.back() << list_indj.back();
1684
1685      if (isCompressible_)
1686      {
1687        std::vector<int>& writtenIndSrc = indWrittenSrv_[rank];
1688        list_writtenInd.push_back(CArray<int,1>(writtenIndSrc.size()));
1689        CArray<int,1>& writtenInd = list_writtenInd.back();
1690
1691        for (n = 0; n < writtenInd.numElements(); ++n)
1692          writtenInd(n) = writtenIndSrc[n];
1693
1694        list_msgsIndex.back() << writtenInd;
1695      }
1696
1697      eventIndex.push(rank, nbConnectedClients_[rank], list_msgsIndex.back());
1698    }
1699
1700    client->sendEvent(eventIndex);
1701  }
1702
1703  /*!
1704    Send area from client to server(s)
1705  */
1706  void CDomain::sendArea()
1707  {
1708    if (!hasArea) return;
1709
1710    int ns, n, i, j, ind, nv, idx;
1711    CContext* context = CContext::getCurrent();
1712    CContextClient* client=context->client;
1713
1714    // send area for each connected server
1715    CEventClient eventArea(getType(), EVENT_ID_AREA);
1716
1717    list<CMessage> list_msgsArea;
1718    list<CArray<double,1> > list_area;
1719
1720    std::map<int, std::vector<size_t> >::const_iterator it, iteMap;
1721    iteMap = indSrv_.end();
1722    for (int k = 0; k < connectedServerRank_.size(); ++k)
1723    {
1724      int nbData = 0;
1725      int rank = connectedServerRank_[k];
1726      it = indSrv_.find(rank);
1727      if (iteMap != it)
1728        nbData = it->second.size();
1729      list_area.push_back(CArray<double,1>(nbData));
1730
1731      const std::vector<size_t>& temp = it->second;
1732      for (n = 0; n < nbData; ++n)
1733      {
1734        idx = static_cast<int>(it->second[n]);
1735        i = i_index(idx);
1736        j = j_index(idx);
1737        if (hasArea)
1738          list_area.back()(n) = area(i - ibegin, j - jbegin);
1739      }
1740
1741      list_msgsArea.push_back(CMessage());
1742      list_msgsArea.back() << this->getId() << list_area.back();
1743      eventArea.push(rank, nbConnectedClients_[rank], list_msgsArea.back());
1744    }
1745    client->sendEvent(eventArea);
1746  }
1747
1748  /*!
1749    Send longitude and latitude from client to servers
1750    Each client send long and lat information to corresponding connected server(s).
1751    Because longitude and latitude are optional, this function only called if latitude and longitude exist
1752  */
1753  void CDomain::sendLonLat()
1754  {
1755    if (!hasLonLat) return;
1756
1757    int ns, n, i, j, ind, nv, idx;
1758    CContext* context = CContext::getCurrent();
1759    CContextClient* client=context->client;
1760
1761    // send lon lat for each connected server
1762    CEventClient eventLon(getType(), EVENT_ID_LON);
1763    CEventClient eventLat(getType(), EVENT_ID_LAT);
1764
1765    list<CMessage> list_msgsLon, list_msgsLat;
1766    list<CArray<double,1> > list_lon, list_lat;
1767    list<CArray<double,2> > list_boundslon, list_boundslat;
1768
1769    std::map<int, std::vector<size_t> >::const_iterator it, iteMap;
1770    iteMap = indSrv_.end();
1771    for (int k = 0; k < connectedServerRank_.size(); ++k)
1772    {
1773      int nbData = 0;
1774      int rank = connectedServerRank_[k];
1775      it = indSrv_.find(rank);
1776      if (iteMap != it)
1777        nbData = it->second.size();
1778
1779      list_lon.push_back(CArray<double,1>(nbData));
1780      list_lat.push_back(CArray<double,1>(nbData));
1781
1782      if (hasBounds)
1783      {
1784        list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
1785        list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
1786      }
1787
1788      CArray<double,1>& lon = list_lon.back();
1789      CArray<double,1>& lat = list_lat.back();
1790      const std::vector<size_t>& temp = it->second;
1791      for (n = 0; n < nbData; ++n)
1792      {
1793        idx = static_cast<int>(it->second[n]);
1794        lon(n) = lonvalue_client(idx);
1795        lat(n) = latvalue_client(idx);
1796
1797        if (hasBounds)
1798        {
1799          CArray<double,2>& boundslon = list_boundslon.back();
1800          CArray<double,2>& boundslat = list_boundslat.back();
1801
1802          for (nv = 0; nv < nvertex; ++nv)
1803          {
1804            boundslon(nv, n) = bounds_lon_client(nv, idx);
1805            boundslat(nv, n) = bounds_lat_client(nv, idx);
1806          }
1807        }
1808      }
1809
1810      list_msgsLon.push_back(CMessage());
1811      list_msgsLat.push_back(CMessage());
1812
1813      list_msgsLon.back() << this->getId() << list_lon.back();
1814      list_msgsLat.back() << this->getId() << list_lat.back();
1815
1816      if (hasBounds)
1817      {
1818        list_msgsLon.back() << list_boundslon.back();
1819        list_msgsLat.back() << list_boundslat.back();
1820      }
1821
1822      eventLon.push(rank, nbConnectedClients_[rank], list_msgsLon.back());
1823      eventLat.push(rank, nbConnectedClients_[rank], list_msgsLat.back());
1824    }
1825
1826    client->sendEvent(eventLon);
1827    client->sendEvent(eventLat);
1828  }
1829
1830  /*!
1831    Send some optional information to server(s)
1832    In the future, this function can be extended with more optional information to send
1833  */
1834  void CDomain::sendLonLatArea(void)
1835  {
1836    sendIndex();
1837    sendLonLat();
1838    sendArea();
1839  }
1840
1841  bool CDomain::dispatchEvent(CEventServer& event)
1842  {
1843    if (SuperClass::dispatchEvent(event)) return true;
1844    else
1845    {
1846      switch(event.type)
1847      {
1848        case EVENT_ID_SERVER_ATTRIBUT:
1849          recvServerAttribut(event);
1850          return true;
1851          break;
1852        case EVENT_ID_INDEX:
1853          recvIndex(event);
1854          return true;
1855          break;
1856        case EVENT_ID_LON:
1857          recvLon(event);
1858          return true;
1859          break;
1860        case EVENT_ID_LAT:
1861          recvLat(event);
1862          return true;
1863          break;
1864        case EVENT_ID_AREA:
1865          recvArea(event);
1866          return true;
1867          break;
1868        default:
1869          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
1870                << "Unknown Event");
1871          return false;
1872       }
1873    }
1874  }
1875
1876  /*!
1877    Receive attributes event from clients(s)
1878    \param[in] event event contain info about rank and associated attributes
1879  */
1880  void CDomain::recvServerAttribut(CEventServer& event)
1881  {
1882    CBufferIn* buffer=event.subEvents.begin()->buffer;
1883    string domainId ;
1884    *buffer>>domainId ;
1885    get(domainId)->recvServerAttribut(*buffer) ;
1886  }
1887
1888  /*!
1889    Receive attributes from client(s): zoom info and begin and n of each server
1890    \param[in] rank rank of client source
1891    \param[in] buffer message containing attributes info
1892  */
1893  void CDomain::recvServerAttribut(CBufferIn& buffer)
1894  {
1895    int global_zoom_ni_tmp, global_zoom_ibegin_tmp, global_zoom_nj_tmp, global_zoom_jbegin_tmp;
1896    buffer >> ni_srv >> ibegin_srv >> iend_srv >> nj_srv >> jbegin_srv >> jend_srv
1897           >> global_zoom_ni_tmp >> global_zoom_ibegin_tmp >> global_zoom_nj_tmp >> global_zoom_jbegin_tmp
1898           >> isCompressible_;
1899
1900    global_zoom_ni.setValue(global_zoom_ni_tmp);
1901    global_zoom_ibegin.setValue(global_zoom_ibegin_tmp);
1902    global_zoom_nj.setValue(global_zoom_nj_tmp);
1903    global_zoom_jbegin.setValue(global_zoom_jbegin_tmp);
1904
1905    int zoom_iend = global_zoom_ibegin + global_zoom_ni - 1;
1906    int zoom_jend = global_zoom_jbegin + global_zoom_nj - 1;
1907
1908    zoom_ibegin_srv = global_zoom_ibegin > ibegin_srv ? global_zoom_ibegin : ibegin_srv ;
1909    zoom_iend_srv = zoom_iend < iend_srv ? zoom_iend : iend_srv ;
1910    zoom_ni_srv=zoom_iend_srv-zoom_ibegin_srv+1 ;
1911
1912    zoom_jbegin_srv = global_zoom_jbegin > jbegin_srv ? global_zoom_jbegin : jbegin_srv ;
1913    zoom_jend_srv = zoom_jend < jend_srv ? zoom_jend : jend_srv ;
1914    zoom_nj_srv=zoom_jend_srv-zoom_jbegin_srv+1 ;
1915
1916    if (zoom_ni_srv<=0 || zoom_nj_srv<=0)
1917    {
1918      zoom_ibegin_srv=0 ; zoom_iend_srv=0 ; zoom_ni_srv=0 ;
1919      zoom_jbegin_srv=0 ; zoom_jend_srv=0 ; zoom_nj_srv=0 ;
1920    }
1921    lonvalue_srv.resize(zoom_ni_srv*zoom_nj_srv) ;
1922    lonvalue_srv = 0. ;
1923    latvalue_srv.resize(zoom_ni_srv*zoom_nj_srv) ;
1924    latvalue_srv = 0. ;
1925    if (hasBounds)
1926    {
1927      bounds_lon_srv.resize(nvertex,zoom_ni_srv*zoom_nj_srv) ;
1928      bounds_lon_srv = 0. ;
1929      bounds_lat_srv.resize(nvertex,zoom_ni_srv*zoom_nj_srv) ;
1930      bounds_lat_srv = 0. ;
1931    }
1932
1933    if (hasArea)
1934    {
1935      area_srv.resize(zoom_ni_srv * zoom_nj_srv);
1936      area_srv = 0.;
1937    }
1938
1939  }
1940
1941  /*!
1942    Receive index event from clients(s)
1943    \param[in] event event contain info about rank and associated index
1944  */
1945  void CDomain::recvIndex(CEventServer& event)
1946  {
1947    CDomain* domain;
1948
1949    list<CEventServer::SSubEvent>::iterator it;
1950    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
1951    {
1952      CBufferIn* buffer = it->buffer;
1953      string domainId;
1954      *buffer >> domainId;
1955      domain = get(domainId);
1956      domain->recvIndex(it->rank, *buffer);
1957    }
1958
1959    if (domain->isCompressible_)
1960    {
1961      std::sort(domain->indexesToWrite.begin(), domain->indexesToWrite.end());
1962
1963      CContextServer* server = CContext::getCurrent()->server;
1964      domain->numberWrittenIndexes_ = domain->indexesToWrite.size();
1965      MPI_Allreduce(&domain->numberWrittenIndexes_, &domain->totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1966      MPI_Scan(&domain->numberWrittenIndexes_, &domain->offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1967      domain->offsetWrittenIndexes_ -= domain->numberWrittenIndexes_;
1968    }
1969  }
1970
1971  /*!
1972    Receive index information from client(s)
1973    \param[in] rank rank of client source
1974    \param[in] buffer message containing index info
1975  */
1976  void CDomain::recvIndex(int rank, CBufferIn& buffer)
1977  {
1978    int type_int;
1979    buffer >> type_int >> isCurvilinear >> indiSrv[rank] >> indjSrv[rank];
1980    type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
1981
1982    if (isCompressible_)
1983    {
1984      CArray<int, 1> writtenIndexes;
1985      buffer >> writtenIndexes;
1986      indexesToWrite.reserve(indexesToWrite.size() + writtenIndexes.numElements());
1987      for (int i = 0; i < writtenIndexes.numElements(); ++i)
1988        indexesToWrite.push_back(writtenIndexes(i));
1989    }
1990  }
1991
1992  /*!
1993    Receive longitude event from clients(s)
1994    \param[in] event event contain info about rank and associated longitude
1995  */
1996  void CDomain::recvLon(CEventServer& event)
1997  {
1998    list<CEventServer::SSubEvent>::iterator it;
1999    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2000    {
2001      CBufferIn* buffer = it->buffer;
2002      string domainId;
2003      *buffer >> domainId;
2004      get(domainId)->recvLon(it->rank, *buffer);
2005    }
2006  }
2007
2008  /*!
2009    Receive longitude information from client(s)
2010    \param[in] rank rank of client source
2011    \param[in] buffer message containing longitude info
2012  */
2013  void CDomain::recvLon(int rank, CBufferIn& buffer)
2014  {
2015    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2016    CArray<double,1> lon;
2017    CArray<double,2> boundslon;
2018
2019    buffer >> lon;
2020
2021    if (hasBounds) buffer >> boundslon;
2022
2023    int i, j, ind_srv;
2024    for (int ind = 0; ind < indi.numElements(); ind++)
2025    {
2026      i = indi(ind); j = indj(ind);
2027      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2028      lonvalue_srv(ind_srv) = lon(ind);
2029      if (hasBounds)
2030      {
2031        for (int nv = 0; nv < nvertex; ++nv)
2032          bounds_lon_srv(nv, ind_srv) = boundslon(nv, ind);
2033      }
2034    }
2035  }
2036
2037  /*!
2038    Receive latitude event from clients(s)
2039    \param[in] event event contain info about rank and associated latitude
2040  */
2041  void CDomain::recvLat(CEventServer& event)
2042  {
2043    list<CEventServer::SSubEvent>::iterator it;
2044    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2045    {
2046      CBufferIn* buffer = it->buffer;
2047      string domainId;
2048      *buffer >> domainId;
2049      get(domainId)->recvLat(it->rank, *buffer);
2050    }
2051  }
2052
2053  /*!
2054    Receive latitude information from client(s)
2055    \param[in] rank rank of client source
2056    \param[in] buffer message containing latitude info
2057  */
2058  void CDomain::recvLat(int rank, CBufferIn& buffer)
2059  {
2060    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2061    CArray<double,1> lat;
2062    CArray<double,2> boundslat;
2063
2064    buffer >> lat;
2065    if (hasBounds) buffer >> boundslat;
2066
2067    int i, j, ind_srv;
2068    for (int ind = 0; ind < indi.numElements(); ind++)
2069    {
2070      i = indi(ind); j = indj(ind);
2071      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2072      latvalue_srv(ind_srv) = lat(ind);
2073      if (hasBounds)
2074      {
2075        for (int nv = 0; nv < nvertex; nv++)
2076          bounds_lat_srv(nv, ind_srv) = boundslat(nv, ind);
2077      }
2078    }
2079  }
2080
2081  /*!
2082    Receive area event from clients(s)
2083    \param[in] event event contain info about rank and associated area
2084  */
2085  void CDomain::recvArea(CEventServer& event)
2086  {
2087    list<CEventServer::SSubEvent>::iterator it;
2088    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2089    {
2090      CBufferIn* buffer = it->buffer;
2091      string domainId;
2092      *buffer >> domainId;
2093      get(domainId)->recvArea(it->rank, *buffer);
2094    }
2095  }
2096
2097  /*!
2098    Receive area information from client(s)
2099    \param[in] rank rank of client source
2100    \param[in] buffer message containing area info
2101  */
2102  void CDomain::recvArea(int rank, CBufferIn& buffer)
2103  {
2104    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2105    CArray<double,1> clientArea;
2106
2107    buffer >> clientArea;
2108
2109    int i, j, ind_srv;
2110    for (int ind = 0; ind < indi.numElements(); ind++)
2111    {
2112      i = indi(ind); j = indj(ind);
2113      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2114      area_srv(ind_srv) = clientArea(ind);
2115    }
2116  }
2117
2118  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2119  {
2120    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2121    return transformationMap_.back().second;
2122  }
2123
2124  /*!
2125    Check whether a domain has transformation
2126    \return true if domain has transformation
2127  */
2128  bool CDomain::hasTransformation()
2129  {
2130    return (!transformationMap_.empty());
2131  }
2132
2133  /*!
2134    Set transformation for current domain. It's the method to move transformation in hierarchy
2135    \param [in] domTrans transformation on domain
2136  */
2137  void CDomain::setTransformations(const TransMapTypes& domTrans)
2138  {
2139    transformationMap_ = domTrans;
2140  }
2141
2142  /*!
2143    Get all transformation current domain has
2144    \return all transformation
2145  */
2146  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2147  {
2148    return transformationMap_;
2149  }
2150
2151  /*!
2152    Check the validity of all transformations applied on domain
2153  This functions is called AFTER all inherited attributes are solved
2154  */
2155  void CDomain::checkTransformations()
2156  {
2157    TransMapTypes::const_iterator itb = transformationMap_.begin(), it,
2158                                  ite = transformationMap_.end();
2159//    for (it = itb; it != ite; ++it)
2160//    {
2161//      (it->second)->checkValid(this);
2162//    }
2163  }
2164
2165  void CDomain::duplicateTransformation(CDomain* src)
2166  {
2167    if (src->hasTransformation())
2168    {
2169      this->setTransformations(src->getAllTransformations());
2170    }
2171  }
2172
2173  /*!
2174   * Go through the hierarchy to find the domain from which the transformations must be inherited
2175   */
2176  void CDomain::solveInheritanceTransformation()
2177  {
2178    if (hasTransformation() || !hasDirectDomainReference())
2179      return;
2180
2181    CDomain* domain = this;
2182    std::vector<CDomain*> refDomains;
2183    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2184    {
2185      refDomains.push_back(domain);
2186      domain = domain->getDirectDomainReference();
2187    }
2188
2189    if (domain->hasTransformation())
2190      for (size_t i = 0; i < refDomains.size(); ++i)
2191        refDomains[i]->setTransformations(domain->getAllTransformations());
2192  }
2193
2194  /*!
2195    Parse children nodes of a domain in xml file.
2196    Whenver there is a new transformation, its type and name should be added into this function
2197    \param node child node to process
2198  */
2199  void CDomain::parse(xml::CXMLNode & node)
2200  {
2201    SuperClass::parse(node);
2202
2203    if (node.goToChildElement())
2204    {
2205      StdString nodeElementName;
2206      do
2207      {
2208        StdString nodeId("");
2209        if (node.getAttributes().end() != node.getAttributes().find("id"))
2210        { nodeId = node.getAttributes()["id"]; }
2211
2212        nodeElementName = node.getElementName();
2213        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2214        it = transformationMapList_.find(nodeElementName);
2215        if (ite != it)
2216        {
2217          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2218                                                                                                                nodeId,
2219                                                                                                                &node)));
2220        }
2221        else
2222        {
2223          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2224                << "The transformation " << nodeElementName << " has not been supported yet.");
2225        }
2226      } while (node.goToNextElement()) ;
2227      node.goToParentElement();
2228    }
2229  }
2230   //----------------------------------------------------------------
2231
2232   DEFINE_REF_FUNC(Domain,domain)
2233
2234   ///---------------------------------------------------------------
2235
2236} // namespace xios
Note: See TracBrowser for help on using the repository browser.