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

Last change on this file since 987 was 987, checked in by mhnguyen, 7 years ago

First working version.
There are more things to do with it

  • 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: 77.7 KB
Line 
1#include "domain.hpp"
2
3#include "attribute_template.hpp"
4#include "object_template.hpp"
5#include "group_template.hpp"
6
7#include "xios_spl.hpp"
8#include "event_client.hpp"
9#include "event_server.hpp"
10#include "buffer_in.hpp"
11#include "message.hpp"
12#include "type.hpp"
13#include "context.hpp"
14#include "context_client.hpp"
15#include "context_server.hpp"
16#include "array_new.hpp"
17#include "distribution_client.hpp"
18#include "server_distribution_description.hpp"
19#include "client_server_mapping_distributed.hpp"
20#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     CContext* context=CContext::getCurrent() ;
1306
1307     if (this->isClientAfterTransformationChecked) return;
1308     if (context->hasClient)
1309     {
1310       this->checkMask();
1311       if (hasLonLat || hasArea || isCompressible_) this->computeConnectedServer();
1312       if (hasLonLat) this->completeLonLatClient();
1313     }
1314
1315     this->isClientAfterTransformationChecked = true;
1316   }
1317
1318   //----------------------------------------------------------------
1319   // Divide function checkAttributes into 2 seperate ones
1320   // This function only checks all attributes of current domain
1321   void CDomain::checkAttributesOnClient()
1322   {
1323     if (this->isClientChecked) return;
1324     CContext* context=CContext::getCurrent();
1325
1326      this->checkDomain();
1327      this->checkBounds();
1328      this->checkArea();
1329      this->checkLonLat();
1330
1331      if (context->hasClient)
1332      { // CÃŽté client uniquement
1333         this->checkMask();
1334         this->checkDomainData();
1335         this->checkCompression();
1336         this->computeLocalMask() ;
1337      }
1338      else
1339      { // CÃŽté serveur uniquement
1340      }
1341
1342      this->isClientChecked = true;
1343   }
1344
1345   // Send all checked attributes to server
1346   void CDomain::sendCheckedAttributes()
1347   {
1348     if (!this->isClientChecked) checkAttributesOnClient();
1349     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1350     CContext* context=CContext::getCurrent() ;
1351
1352     if (this->isChecked) return;
1353     if (context->hasClient)
1354     {
1355       sendServerAttribut();
1356       if (hasLonLat || hasArea || isCompressible_) sendLonLatArea();
1357     }
1358     this->isChecked = true;
1359   }
1360
1361   void CDomain::checkAttributes(void)
1362   {
1363      if (this->isChecked) return;
1364      CContext* context=CContext::getCurrent() ;
1365
1366      this->checkDomain();
1367      this->checkLonLat();
1368      this->checkBounds();
1369      this->checkArea();
1370
1371      if (context->hasClient)
1372      { // CÃŽté client uniquement
1373         this->checkMask();
1374         this->checkDomainData();
1375         this->checkCompression();
1376         this->computeLocalMask() ;
1377
1378      }
1379      else
1380      { // CÃŽté serveur uniquement
1381      }
1382
1383      if (context->hasClient)
1384      {
1385        this->computeConnectedServer();
1386        this->completeLonLatClient();
1387        this->sendServerAttribut();
1388        this->sendLonLatArea();
1389      }
1390
1391      this->isChecked = true;
1392   }
1393
1394  void CDomain::sendServerAttribut(void)
1395  {
1396    CContext* context = CContext::getCurrent();
1397     // Use correct context client to send message
1398     CContextClient* contextClientTmp = (0 != context->clientPrimServer) ? context->clientPrimServer
1399                                                                         : context->client;
1400
1401    // CContextClient* client = context->client;
1402    int nbServer = contextClientTmp->serverSize;
1403
1404    CServerDistributionDescription serverDescription(nGlobDomain_, nbServer);
1405    if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
1406    else serverDescription.computeServerDistribution(false, 1);
1407
1408    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
1409    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
1410
1411    CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
1412    if (contextClientTmp->isServerLeader())
1413    {
1414      std::list<CMessage> msgs;
1415
1416      const std::list<int>& ranks = contextClientTmp->getRanksServerLeader();
1417      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
1418      {
1419        // Use const int to ensure CMessage holds a copy of the value instead of just a reference
1420        const int ibegin_srv = serverIndexBegin[*itRank][0];
1421        const int jbegin_srv = serverIndexBegin[*itRank][1];
1422        const int ni_srv = serverDimensionSizes[*itRank][0];
1423        const int nj_srv = serverDimensionSizes[*itRank][1];
1424        const int iend_srv = ibegin_srv + ni_srv - 1;
1425        const int jend_srv = jbegin_srv + nj_srv - 1;
1426
1427        msgs.push_back(CMessage());
1428        CMessage& msg = msgs.back();
1429        msg << this->getId() ;
1430        msg << ni_srv << ibegin_srv << iend_srv << nj_srv << jbegin_srv << jend_srv;
1431        msg << global_zoom_ni.getValue() << global_zoom_ibegin.getValue() << global_zoom_nj.getValue() << global_zoom_jbegin.getValue();
1432        msg << isCompressible_;
1433
1434        event.push(*itRank,1,msg);
1435      }
1436      contextClientTmp->sendEvent(event);
1437    }
1438    else contextClientTmp->sendEvent(event);
1439  }
1440
1441  void CDomain::computeNGlobDomain()
1442  {
1443    nGlobDomain_.resize(2);
1444    nGlobDomain_[0] = ni_glo.getValue();
1445    nGlobDomain_[1] = nj_glo.getValue();
1446  }
1447
1448  void CDomain::computeConnectedServer(void)
1449  {
1450    CContext* context=CContext::getCurrent() ;
1451    CContextClient* client=context->client ;
1452    int nbServer=client->serverSize;
1453    int rank = client->clientRank;
1454    bool doComputeGlobalIndexServer = true;
1455
1456    int i,j,i_ind,j_ind, nbIndex;
1457    int global_zoom_iend=global_zoom_ibegin+global_zoom_ni-1 ;
1458    int global_zoom_jend=global_zoom_jbegin+global_zoom_nj-1 ;
1459
1460    // Precompute number of index
1461    int globalIndexCountZoom = 0;
1462    nbIndex = i_index.numElements();
1463    for (i = 0; i < nbIndex; ++i)
1464    {
1465      i_ind=i_index(i);
1466      j_ind=j_index(i);
1467
1468      if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1469      {
1470        ++globalIndexCountZoom;
1471      }
1472    }
1473
1474    int globalIndexWrittenCount = 0;
1475    if (isCompressible_)
1476    {
1477      for (i = 0; i < data_i_index.numElements(); ++i)
1478      {
1479        i_ind = CDistributionClient::getDomainIndex(data_i_index(i), data_j_index(i),
1480                                                    data_ibegin, data_jbegin, data_dim, ni,
1481                                                    j_ind);
1482        if (i_ind >= 0 && i_ind < ni && j_ind >= 0 && j_ind < nj && mask_1d(i_ind + j_ind * ni))
1483        {
1484          i_ind += ibegin;
1485          j_ind += jbegin;
1486          if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1487            ++globalIndexWrittenCount;
1488        }
1489      }
1490    }
1491
1492    // Fill in index
1493    CArray<size_t,1> globalIndexDomainZoom(globalIndexCountZoom);
1494    CArray<size_t,1> localIndexDomainZoom(globalIndexCountZoom);
1495    CArray<size_t,1> globalIndexDomain(nbIndex);
1496    size_t globalIndex;
1497    int globalIndexCount = 0;
1498    globalIndexCountZoom = 0;
1499
1500    for (i = 0; i < nbIndex; ++i)
1501    {
1502      i_ind=i_index(i);
1503      j_ind=j_index(i);
1504      globalIndex = i_ind + j_ind * ni_glo;
1505      globalIndexDomain(globalIndexCount) = globalIndex;
1506      ++globalIndexCount;
1507      if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1508      {
1509        globalIndexDomainZoom(globalIndexCountZoom) = globalIndex;
1510        localIndexDomainZoom(globalIndexCountZoom) = i;
1511        ++globalIndexCountZoom;
1512      }
1513    }
1514
1515    CArray<int,1> globalIndexWrittenDomain(globalIndexWrittenCount);
1516    if (isCompressible_)
1517    {
1518      globalIndexWrittenCount = 0;
1519      for (i = 0; i < data_i_index.numElements(); ++i)
1520      {
1521        i_ind = CDistributionClient::getDomainIndex(data_i_index(i), data_j_index(i),
1522                                                    data_ibegin, data_jbegin, data_dim, ni,
1523                                                    j_ind);
1524        if (i_ind >= 0 && i_ind < ni && j_ind >= 0 && j_ind < nj && mask_1d(i_ind + j_ind * ni))
1525        {
1526          i_ind += ibegin;
1527          j_ind += jbegin;
1528          if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1529          {
1530            globalIndexWrittenDomain(globalIndexWrittenCount) = i_ind + j_ind * ni_glo;
1531            ++globalIndexWrittenCount;
1532          }
1533        }
1534      }
1535    }
1536
1537    size_t globalSizeIndex = 1, indexBegin, indexEnd;
1538    int range, clientSize = client->clientSize;
1539    for (int i = 0; i < nGlobDomain_.size(); ++i) globalSizeIndex *= nGlobDomain_[i];
1540    indexBegin = 0;
1541    if (globalSizeIndex <= clientSize)
1542    {
1543      indexBegin = rank%globalSizeIndex;
1544      indexEnd = indexBegin;
1545    }
1546    else
1547    {
1548      for (int i = 0; i < clientSize; ++i)
1549      {
1550        range = globalSizeIndex / clientSize;
1551        if (i < (globalSizeIndex%clientSize)) ++range;
1552        if (i == client->clientRank) break;
1553        indexBegin += range;
1554      }
1555      indexEnd = indexBegin + range - 1;
1556    }
1557
1558    CServerDistributionDescription serverDescription(nGlobDomain_, nbServer);
1559    if (isUnstructed_) serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 0);
1560    else serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 1);
1561
1562    CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(),
1563                                                                                client->intraComm);
1564    clientServerMap->computeServerIndexMapping(globalIndexDomain);
1565    const CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
1566
1567    CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
1568                                                         ite = globalIndexDomainOnServer.end();
1569    typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
1570    std::vector<int>::iterator itVec;
1571
1572    indSrv_.clear();
1573    indWrittenSrv_.clear();
1574    for (; it != ite; ++it)
1575    {
1576      int rank = it->first;
1577      int indexSize = it->second.size();
1578      std::vector<int> permutIndex(indexSize);
1579      XIOSAlgorithms::fillInIndex(indexSize, permutIndex);
1580      XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(it->second, permutIndex);
1581      BinarySearch binSearch(it->second);
1582      int nb = globalIndexDomainZoom.numElements();
1583      for (int i = 0; i < nb; ++i)
1584      {
1585        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndexDomainZoom(i), itVec))
1586        {
1587          indSrv_[rank].push_back(localIndexDomainZoom(i));
1588        }
1589      }
1590      for (int i = 0; i < globalIndexWrittenDomain.numElements(); ++i)
1591      {
1592        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndexWrittenDomain(i), itVec))
1593        {
1594          indWrittenSrv_[rank].push_back(globalIndexWrittenDomain(i));
1595        }
1596      }
1597    }
1598
1599    connectedServerRank_.clear();
1600    for (it = globalIndexDomainOnServer.begin(); it != ite; ++it) {
1601      connectedServerRank_.push_back(it->first);
1602    }
1603
1604    nbConnectedClients_ = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_);
1605
1606    delete clientServerMap;
1607  }
1608
1609  const std::map<int, vector<size_t> >& CDomain::getIndexServer() const
1610  {
1611    return indSrv_;
1612  }
1613
1614  /*!
1615    Send index from client to server(s)
1616  */
1617  void CDomain::sendIndex()
1618  {
1619    int ns, n, i, j, ind, nv, idx;
1620    CContext* context = CContext::getCurrent();
1621    CContextClient* client=context->client;
1622
1623    CEventClient eventIndex(getType(), EVENT_ID_INDEX);
1624
1625    list<CMessage> list_msgsIndex;
1626    list<CArray<int,1> > list_indi, list_indj, list_writtenInd;
1627
1628    std::map<int, std::vector<size_t> >::const_iterator it, iteMap;
1629    iteMap = indSrv_.end();
1630    for (int k = 0; k < connectedServerRank_.size(); ++k)
1631    {
1632      int nbData = 0;
1633      int rank = connectedServerRank_[k];
1634      it = indSrv_.find(rank);
1635      if (iteMap != it)
1636        nbData = it->second.size();
1637
1638      list_indi.push_back(CArray<int,1>(nbData));
1639      list_indj.push_back(CArray<int,1>(nbData));
1640
1641      CArray<int,1>& indi = list_indi.back();
1642      CArray<int,1>& indj = list_indj.back();
1643      const std::vector<size_t>& temp = it->second;
1644      for (n = 0; n < nbData; ++n)
1645      {
1646        idx = static_cast<int>(it->second[n]);
1647        indi(n) = i_index(idx);
1648        indj(n) = j_index(idx);
1649      }
1650
1651      list_msgsIndex.push_back(CMessage());
1652
1653      list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
1654      list_msgsIndex.back() << isCurvilinear;
1655      list_msgsIndex.back() << list_indi.back() << list_indj.back();
1656
1657      if (isCompressible_)
1658      {
1659        std::vector<int>& writtenIndSrc = indWrittenSrv_[rank];
1660        list_writtenInd.push_back(CArray<int,1>(writtenIndSrc.size()));
1661        CArray<int,1>& writtenInd = list_writtenInd.back();
1662
1663        for (n = 0; n < writtenInd.numElements(); ++n)
1664          writtenInd(n) = writtenIndSrc[n];
1665
1666        list_msgsIndex.back() << writtenInd;
1667      }
1668
1669      eventIndex.push(rank, nbConnectedClients_[rank], list_msgsIndex.back());
1670    }
1671
1672    client->sendEvent(eventIndex);
1673  }
1674
1675  /*!
1676    Send area from client to server(s)
1677  */
1678  void CDomain::sendArea()
1679  {
1680    if (!hasArea) return;
1681
1682    int ns, n, i, j, ind, nv, idx;
1683    CContext* context = CContext::getCurrent();
1684    CContextClient* client=context->client;
1685
1686    // send area for each connected server
1687    CEventClient eventArea(getType(), EVENT_ID_AREA);
1688
1689    list<CMessage> list_msgsArea;
1690    list<CArray<double,1> > list_area;
1691
1692    std::map<int, std::vector<size_t> >::const_iterator it, iteMap;
1693    iteMap = indSrv_.end();
1694    for (int k = 0; k < connectedServerRank_.size(); ++k)
1695    {
1696      int nbData = 0;
1697      int rank = connectedServerRank_[k];
1698      it = indSrv_.find(rank);
1699      if (iteMap != it)
1700        nbData = it->second.size();
1701      list_area.push_back(CArray<double,1>(nbData));
1702
1703      const std::vector<size_t>& temp = it->second;
1704      for (n = 0; n < nbData; ++n)
1705      {
1706        idx = static_cast<int>(it->second[n]);
1707        i = i_index(idx);
1708        j = j_index(idx);
1709        if (hasArea)
1710          list_area.back()(n) = area(i - ibegin, j - jbegin);
1711      }
1712
1713      list_msgsArea.push_back(CMessage());
1714      list_msgsArea.back() << this->getId() << list_area.back();
1715      eventArea.push(rank, nbConnectedClients_[rank], list_msgsArea.back());
1716    }
1717    client->sendEvent(eventArea);
1718  }
1719
1720  /*!
1721    Send longitude and latitude from client to servers
1722    Each client send long and lat information to corresponding connected server(s).
1723    Because longitude and latitude are optional, this function only called if latitude and longitude exist
1724  */
1725  void CDomain::sendLonLat()
1726  {
1727    if (!hasLonLat) return;
1728
1729    int ns, n, i, j, ind, nv, idx;
1730    CContext* context = CContext::getCurrent();
1731    CContextClient* client=context->client;
1732
1733    // send lon lat for each connected server
1734    CEventClient eventLon(getType(), EVENT_ID_LON);
1735    CEventClient eventLat(getType(), EVENT_ID_LAT);
1736
1737    list<CMessage> list_msgsLon, list_msgsLat;
1738    list<CArray<double,1> > list_lon, list_lat;
1739    list<CArray<double,2> > list_boundslon, list_boundslat;
1740
1741    std::map<int, std::vector<size_t> >::const_iterator it, iteMap;
1742    iteMap = indSrv_.end();
1743    for (int k = 0; k < connectedServerRank_.size(); ++k)
1744    {
1745      int nbData = 0;
1746      int rank = connectedServerRank_[k];
1747      it = indSrv_.find(rank);
1748      if (iteMap != it)
1749        nbData = it->second.size();
1750
1751      list_lon.push_back(CArray<double,1>(nbData));
1752      list_lat.push_back(CArray<double,1>(nbData));
1753
1754      if (hasBounds)
1755      {
1756        list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
1757        list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
1758      }
1759
1760      CArray<double,1>& lon = list_lon.back();
1761      CArray<double,1>& lat = list_lat.back();
1762      const std::vector<size_t>& temp = it->second;
1763      for (n = 0; n < nbData; ++n)
1764      {
1765        idx = static_cast<int>(it->second[n]);
1766        lon(n) = lonvalue_client(idx);
1767        lat(n) = latvalue_client(idx);
1768
1769        if (hasBounds)
1770        {
1771          CArray<double,2>& boundslon = list_boundslon.back();
1772          CArray<double,2>& boundslat = list_boundslat.back();
1773
1774          for (nv = 0; nv < nvertex; ++nv)
1775          {
1776            boundslon(nv, n) = bounds_lon_client(nv, idx);
1777            boundslat(nv, n) = bounds_lat_client(nv, idx);
1778          }
1779        }
1780      }
1781
1782      list_msgsLon.push_back(CMessage());
1783      list_msgsLat.push_back(CMessage());
1784
1785      list_msgsLon.back() << this->getId() << list_lon.back();
1786      list_msgsLat.back() << this->getId() << list_lat.back();
1787
1788      if (hasBounds)
1789      {
1790        list_msgsLon.back() << list_boundslon.back();
1791        list_msgsLat.back() << list_boundslat.back();
1792      }
1793
1794      eventLon.push(rank, nbConnectedClients_[rank], list_msgsLon.back());
1795      eventLat.push(rank, nbConnectedClients_[rank], list_msgsLat.back());
1796    }
1797
1798    client->sendEvent(eventLon);
1799    client->sendEvent(eventLat);
1800  }
1801
1802  /*!
1803    Send some optional information to server(s)
1804    In the future, this function can be extended with more optional information to send
1805  */
1806  void CDomain::sendLonLatArea(void)
1807  {
1808    sendIndex();
1809    sendLonLat();
1810    sendArea();
1811  }
1812
1813  bool CDomain::dispatchEvent(CEventServer& event)
1814  {
1815    if (SuperClass::dispatchEvent(event)) return true;
1816    else
1817    {
1818      switch(event.type)
1819      {
1820        case EVENT_ID_SERVER_ATTRIBUT:
1821          recvServerAttribut(event);
1822          return true;
1823          break;
1824        case EVENT_ID_INDEX:
1825          recvIndex(event);
1826          return true;
1827          break;
1828        case EVENT_ID_LON:
1829          recvLon(event);
1830          return true;
1831          break;
1832        case EVENT_ID_LAT:
1833          recvLat(event);
1834          return true;
1835          break;
1836        case EVENT_ID_AREA:
1837          recvArea(event);
1838          return true;
1839          break;
1840        default:
1841          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
1842                << "Unknown Event");
1843          return false;
1844       }
1845    }
1846  }
1847
1848  /*!
1849    Receive attributes event from clients(s)
1850    \param[in] event event contain info about rank and associated attributes
1851  */
1852  void CDomain::recvServerAttribut(CEventServer& event)
1853  {
1854    CBufferIn* buffer=event.subEvents.begin()->buffer;
1855    string domainId ;
1856    *buffer>>domainId ;
1857    get(domainId)->recvServerAttribut(*buffer) ;
1858   
1859    CContext* context = CContext::getCurrent();
1860    if (context->hasClient && context->hasServer)
1861    {
1862      get(domainId)->sendServerAttribut();
1863    }
1864
1865  }
1866
1867  /*!
1868    Receive attributes from client(s): zoom info and begin and n of each server
1869    \param[in] rank rank of client source
1870    \param[in] buffer message containing attributes info
1871  */
1872  void CDomain::recvServerAttribut(CBufferIn& buffer)
1873  {
1874    int global_zoom_ni_tmp, global_zoom_ibegin_tmp, global_zoom_nj_tmp, global_zoom_jbegin_tmp;
1875    buffer >> ni_srv >> ibegin_srv >> iend_srv >> nj_srv >> jbegin_srv >> jend_srv
1876           >> global_zoom_ni_tmp >> global_zoom_ibegin_tmp >> global_zoom_nj_tmp >> global_zoom_jbegin_tmp
1877           >> isCompressible_;
1878
1879    global_zoom_ni.setValue(global_zoom_ni_tmp);
1880    global_zoom_ibegin.setValue(global_zoom_ibegin_tmp);
1881    global_zoom_nj.setValue(global_zoom_nj_tmp);
1882    global_zoom_jbegin.setValue(global_zoom_jbegin_tmp);
1883
1884    int zoom_iend = global_zoom_ibegin + global_zoom_ni - 1;
1885    int zoom_jend = global_zoom_jbegin + global_zoom_nj - 1;
1886
1887    zoom_ibegin_srv = global_zoom_ibegin > ibegin_srv ? global_zoom_ibegin : ibegin_srv ;
1888    zoom_iend_srv = zoom_iend < iend_srv ? zoom_iend : iend_srv ;
1889    zoom_ni_srv=zoom_iend_srv-zoom_ibegin_srv+1 ;
1890
1891    zoom_jbegin_srv = global_zoom_jbegin > jbegin_srv ? global_zoom_jbegin : jbegin_srv ;
1892    zoom_jend_srv = zoom_jend < jend_srv ? zoom_jend : jend_srv ;
1893    zoom_nj_srv=zoom_jend_srv-zoom_jbegin_srv+1 ;
1894
1895    if (zoom_ni_srv<=0 || zoom_nj_srv<=0)
1896    {
1897      zoom_ibegin_srv=0 ; zoom_iend_srv=0 ; zoom_ni_srv=0 ;
1898      zoom_jbegin_srv=0 ; zoom_jend_srv=0 ; zoom_nj_srv=0 ;
1899    }
1900    lonvalue_srv.resize(zoom_ni_srv*zoom_nj_srv) ;
1901    lonvalue_srv = 0. ;
1902    latvalue_srv.resize(zoom_ni_srv*zoom_nj_srv) ;
1903    latvalue_srv = 0. ;
1904    if (hasBounds)
1905    {
1906      bounds_lon_srv.resize(nvertex,zoom_ni_srv*zoom_nj_srv) ;
1907      bounds_lon_srv = 0. ;
1908      bounds_lat_srv.resize(nvertex,zoom_ni_srv*zoom_nj_srv) ;
1909      bounds_lat_srv = 0. ;
1910    }
1911
1912    if (hasArea)
1913    {
1914      area_srv.resize(zoom_ni_srv * zoom_nj_srv);
1915      area_srv = 0.;
1916    }
1917
1918  }
1919
1920  /*!
1921    Receive index event from clients(s)
1922    \param[in] event event contain info about rank and associated index
1923  */
1924  void CDomain::recvIndex(CEventServer& event)
1925  {
1926    CDomain* domain;
1927
1928    list<CEventServer::SSubEvent>::iterator it;
1929    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
1930    {
1931      CBufferIn* buffer = it->buffer;
1932      string domainId;
1933      *buffer >> domainId;
1934      domain = get(domainId);
1935      domain->recvIndex(it->rank, *buffer);
1936    }
1937
1938    if (domain->isCompressible_)
1939    {
1940      std::sort(domain->indexesToWrite.begin(), domain->indexesToWrite.end());
1941
1942      CContextServer* server = CContext::getCurrent()->server;
1943      domain->numberWrittenIndexes_ = domain->indexesToWrite.size();
1944      MPI_Allreduce(&domain->numberWrittenIndexes_, &domain->totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1945      MPI_Scan(&domain->numberWrittenIndexes_, &domain->offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1946      domain->offsetWrittenIndexes_ -= domain->numberWrittenIndexes_;
1947    }
1948  }
1949
1950  /*!
1951    Receive index information from client(s)
1952    \param[in] rank rank of client source
1953    \param[in] buffer message containing index info
1954  */
1955  void CDomain::recvIndex(int rank, CBufferIn& buffer)
1956  {
1957    int type_int;
1958    buffer >> type_int >> isCurvilinear >> indiSrv[rank] >> indjSrv[rank];
1959    type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
1960
1961    if (isCompressible_)
1962    {
1963      CArray<int, 1> writtenIndexes;
1964      buffer >> writtenIndexes;
1965      indexesToWrite.reserve(indexesToWrite.size() + writtenIndexes.numElements());
1966      for (int i = 0; i < writtenIndexes.numElements(); ++i)
1967        indexesToWrite.push_back(writtenIndexes(i));
1968    }
1969  }
1970
1971  /*!
1972    Receive longitude event from clients(s)
1973    \param[in] event event contain info about rank and associated longitude
1974  */
1975  void CDomain::recvLon(CEventServer& event)
1976  {
1977    list<CEventServer::SSubEvent>::iterator it;
1978    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
1979    {
1980      CBufferIn* buffer = it->buffer;
1981      string domainId;
1982      *buffer >> domainId;
1983      get(domainId)->recvLon(it->rank, *buffer);
1984    }
1985  }
1986
1987  /*!
1988    Receive longitude information from client(s)
1989    \param[in] rank rank of client source
1990    \param[in] buffer message containing longitude info
1991  */
1992  void CDomain::recvLon(int rank, CBufferIn& buffer)
1993  {
1994    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
1995    CArray<double,1> lon;
1996    CArray<double,2> boundslon;
1997
1998    buffer >> lon;
1999
2000    if (hasBounds) buffer >> boundslon;
2001
2002    int i, j, ind_srv;
2003    for (int ind = 0; ind < indi.numElements(); ind++)
2004    {
2005      i = indi(ind); j = indj(ind);
2006      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2007      lonvalue_srv(ind_srv) = lon(ind);
2008      if (hasBounds)
2009      {
2010        for (int nv = 0; nv < nvertex; ++nv)
2011          bounds_lon_srv(nv, ind_srv) = boundslon(nv, ind);
2012      }
2013    }
2014  }
2015
2016  /*!
2017    Receive latitude event from clients(s)
2018    \param[in] event event contain info about rank and associated latitude
2019  */
2020  void CDomain::recvLat(CEventServer& event)
2021  {
2022    list<CEventServer::SSubEvent>::iterator it;
2023    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2024    {
2025      CBufferIn* buffer = it->buffer;
2026      string domainId;
2027      *buffer >> domainId;
2028      get(domainId)->recvLat(it->rank, *buffer);
2029    }
2030  }
2031
2032  /*!
2033    Receive latitude information from client(s)
2034    \param[in] rank rank of client source
2035    \param[in] buffer message containing latitude info
2036  */
2037  void CDomain::recvLat(int rank, CBufferIn& buffer)
2038  {
2039    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2040    CArray<double,1> lat;
2041    CArray<double,2> boundslat;
2042
2043    buffer >> lat;
2044    if (hasBounds) buffer >> boundslat;
2045
2046    int i, j, ind_srv;
2047    for (int ind = 0; ind < indi.numElements(); ind++)
2048    {
2049      i = indi(ind); j = indj(ind);
2050      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2051      latvalue_srv(ind_srv) = lat(ind);
2052      if (hasBounds)
2053      {
2054        for (int nv = 0; nv < nvertex; nv++)
2055          bounds_lat_srv(nv, ind_srv) = boundslat(nv, ind);
2056      }
2057    }
2058  }
2059
2060  /*!
2061    Receive area event from clients(s)
2062    \param[in] event event contain info about rank and associated area
2063  */
2064  void CDomain::recvArea(CEventServer& event)
2065  {
2066    list<CEventServer::SSubEvent>::iterator it;
2067    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2068    {
2069      CBufferIn* buffer = it->buffer;
2070      string domainId;
2071      *buffer >> domainId;
2072      get(domainId)->recvArea(it->rank, *buffer);
2073    }
2074  }
2075
2076  /*!
2077    Receive area information from client(s)
2078    \param[in] rank rank of client source
2079    \param[in] buffer message containing area info
2080  */
2081  void CDomain::recvArea(int rank, CBufferIn& buffer)
2082  {
2083    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2084    CArray<double,1> clientArea;
2085
2086    buffer >> clientArea;
2087
2088    int i, j, ind_srv;
2089    for (int ind = 0; ind < indi.numElements(); ind++)
2090    {
2091      i = indi(ind); j = indj(ind);
2092      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2093      area_srv(ind_srv) = clientArea(ind);
2094    }
2095  }
2096
2097  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2098  {
2099    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2100    return transformationMap_.back().second;
2101  }
2102
2103  /*!
2104    Check whether a domain has transformation
2105    \return true if domain has transformation
2106  */
2107  bool CDomain::hasTransformation()
2108  {
2109    return (!transformationMap_.empty());
2110  }
2111
2112  /*!
2113    Set transformation for current domain. It's the method to move transformation in hierarchy
2114    \param [in] domTrans transformation on domain
2115  */
2116  void CDomain::setTransformations(const TransMapTypes& domTrans)
2117  {
2118    transformationMap_ = domTrans;
2119  }
2120
2121  /*!
2122    Get all transformation current domain has
2123    \return all transformation
2124  */
2125  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2126  {
2127    return transformationMap_;
2128  }
2129
2130  /*!
2131    Check the validity of all transformations applied on domain
2132  This functions is called AFTER all inherited attributes are solved
2133  */
2134  void CDomain::checkTransformations()
2135  {
2136    TransMapTypes::const_iterator itb = transformationMap_.begin(), it,
2137                                  ite = transformationMap_.end();
2138//    for (it = itb; it != ite; ++it)
2139//    {
2140//      (it->second)->checkValid(this);
2141//    }
2142  }
2143
2144  void CDomain::duplicateTransformation(CDomain* src)
2145  {
2146    if (src->hasTransformation())
2147    {
2148      this->setTransformations(src->getAllTransformations());
2149    }
2150  }
2151
2152  /*!
2153   * Go through the hierarchy to find the domain from which the transformations must be inherited
2154   */
2155  void CDomain::solveInheritanceTransformation()
2156  {
2157    if (hasTransformation() || !hasDirectDomainReference())
2158      return;
2159
2160    CDomain* domain = this;
2161    std::vector<CDomain*> refDomains;
2162    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2163    {
2164      refDomains.push_back(domain);
2165      domain = domain->getDirectDomainReference();
2166    }
2167
2168    if (domain->hasTransformation())
2169      for (size_t i = 0; i < refDomains.size(); ++i)
2170        refDomains[i]->setTransformations(domain->getAllTransformations());
2171  }
2172
2173  /*!
2174    Parse children nodes of a domain in xml file.
2175    Whenver there is a new transformation, its type and name should be added into this function
2176    \param node child node to process
2177  */
2178  void CDomain::parse(xml::CXMLNode & node)
2179  {
2180    SuperClass::parse(node);
2181
2182    if (node.goToChildElement())
2183    {
2184      StdString nodeElementName;
2185      do
2186      {
2187        StdString nodeId("");
2188        if (node.getAttributes().end() != node.getAttributes().find("id"))
2189        { nodeId = node.getAttributes()["id"]; }
2190
2191        nodeElementName = node.getElementName();
2192        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2193        it = transformationMapList_.find(nodeElementName);
2194        if (ite != it)
2195        {
2196          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2197                                                                                                                nodeId,
2198                                                                                                                &node)));
2199        }
2200        else
2201        {
2202          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2203                << "The transformation " << nodeElementName << " has not been supported yet.");
2204        }
2205      } while (node.goToNextElement()) ;
2206      node.goToParentElement();
2207    }
2208  }
2209   //----------------------------------------------------------------
2210
2211   DEFINE_REF_FUNC(Domain,domain)
2212
2213   ///---------------------------------------------------------------
2214
2215} // namespace xios
Note: See TracBrowser for help on using the repository browser.