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

Last change on this file since 1062 was 1062, checked in by ymipsl, 4 years ago

Bug fix in interpolation for cell boundaries generation.

YM

  • 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.5 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_lat_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     // A stupid condition to make sure that if there is only one client, domain
705     // should be considered to be distributed. This should be a temporary solution     
706     isDistributed_ |= (1 == CContext::getCurrent()->client->clientSize);
707   }
708
709   // Check global zoom of a domain
710   // If there is no zoom defined for the domain, zoom will have value of global doamin
711   void CDomain::checkZoom(void)
712   {
713     if (global_zoom_ibegin.isEmpty())
714      global_zoom_ibegin.setValue(0);
715     if (global_zoom_ni.isEmpty())
716      global_zoom_ni.setValue(ni_glo);
717     if (global_zoom_jbegin.isEmpty())
718      global_zoom_jbegin.setValue(0);
719     if (global_zoom_nj.isEmpty())
720      global_zoom_nj.setValue(nj_glo);
721   }
722
723   //----------------------------------------------------------------
724
725   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
726   void CDomain::checkLocalIDomain(void)
727   {
728      // If ibegin and ni are provided then we use them to check the validity of local domain
729      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
730      {
731        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
732        {
733          ERROR("CDomain::checkLocalIDomain(void)",
734                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
735                << "The local domain is wrongly defined,"
736                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
737        }
738      }
739
740      // i_index has higher priority than ibegin and ni
741      if (!i_index.isEmpty())
742      {
743        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
744        if (ni.isEmpty()) 
745        {         
746         // No information about ni
747          int minIndex = ni_glo - 1;
748          int maxIndex = 0;
749          for (int idx = 0; idx < i_index.numElements(); ++idx)
750          {
751            if (i_index(idx) < minIndex) minIndex = i_index(idx);
752            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
753          }
754          ni = maxIndex - minIndex + 1; 
755          minIIndex = minIIndex;         
756        }
757
758        // It's not so correct but if ibegin is not the first value of i_index
759        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
760        if (ibegin.isEmpty()) ibegin = minIIndex;
761      }
762      else if (ibegin.isEmpty() && ni.isEmpty())
763      {
764        ibegin = 0;
765        ni = ni_glo;
766      }
767      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
768      {
769        ERROR("CDomain::checkLocalIDomain(void)",
770              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
771              << "The local domain is wrongly defined," << endl
772              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
773              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
774      }
775       
776
777      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
778      {
779        ERROR("CDomain::checkLocalIDomain(void)",
780              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
781              << "The local domain is wrongly defined,"
782              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
783      }
784   }
785
786   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
787   void CDomain::checkLocalJDomain(void)
788   {
789    // If jbegin and nj are provided then we use them to check the validity of local domain
790     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
791     {
792       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
793       {
794         ERROR("CDomain::checkLocalJDomain(void)",
795                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
796                << "The local domain is wrongly defined,"
797                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
798       }
799     }
800
801     if (!j_index.isEmpty())
802     {
803        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
804        if (nj.isEmpty()) 
805        {
806          // No information about nj
807          int minIndex = nj_glo - 1;
808          int maxIndex = 0;
809          for (int idx = 0; idx < j_index.numElements(); ++idx)
810          {
811            if (j_index(idx) < minIndex) minIndex = j_index(idx);
812            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
813          }
814          nj = maxIndex - minIndex + 1;
815          minJIndex = minIndex; 
816        } 
817        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
818        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
819       if (jbegin.isEmpty()) jbegin = minJIndex;       
820     }
821     else if (jbegin.isEmpty() && nj.isEmpty())
822     {
823       jbegin = 0;
824       nj = nj_glo;
825     }     
826
827
828     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
829     {
830       ERROR("CDomain::checkLocalJDomain(void)",
831              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
832              << "The local domain is wrongly defined,"
833              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
834     }
835   }
836
837   //----------------------------------------------------------------
838
839   void CDomain::checkMask(void)
840   {
841      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
842        ERROR("CDomain::checkMask(void)",
843              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
844              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
845              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
846
847      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
848      {
849        if (mask_1d.numElements() != i_index.numElements())
850          ERROR("CDomain::checkMask(void)",
851                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
852                << "'mask_1d' does not have the same size as the local domain." << std::endl
853                << "Local size is " << i_index.numElements() << "." << std::endl
854                << "Mask size is " << mask_1d.numElements() << ".");
855      }
856
857      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
858      {
859        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
860          ERROR("CDomain::checkMask(void)",
861                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
862                << "The mask does not have the same size as the local domain." << std::endl
863                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
864                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
865      }
866
867      if (!mask_2d.isEmpty())
868      {
869        mask_1d.resize(mask_2d.extent(0) * mask_2d.extent(1));
870        for (int j = 0; j < nj; ++j)
871          for (int i = 0; i < ni; ++i) mask_1d(i+j*ni) = mask_2d(i,j);
872        mask_2d.reset();
873      }
874      else if (mask_1d.isEmpty())
875      {
876        mask_1d.resize(i_index.numElements());
877        for (int i = 0; i < i_index.numElements(); ++i) mask_1d(i) = true;
878      }
879   }
880
881   //----------------------------------------------------------------
882
883   void CDomain::checkDomainData(void)
884   {
885      if (data_dim.isEmpty())
886      {
887        data_dim.setValue(1);
888      }
889      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
890      {
891        ERROR("CDomain::checkDomainData(void)",
892              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
893              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
894      }
895
896      if (data_ibegin.isEmpty())
897         data_ibegin.setValue(0);
898      if (data_jbegin.isEmpty())
899         data_jbegin.setValue(0);
900
901      if (data_ni.isEmpty())
902      {
903        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
904      }
905      else if (data_ni.getValue() < 0)
906      {
907        ERROR("CDomain::checkDomainData(void)",
908              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
909              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
910      }
911
912      if (data_nj.isEmpty())
913      {
914        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
915      }
916      else if (data_nj.getValue() < 0)
917      {
918        ERROR("CDomain::checkDomainData(void)",
919              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
920              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
921      }
922   }
923
924   //----------------------------------------------------------------
925
926   void CDomain::checkCompression(void)
927   {
928      if (!data_i_index.isEmpty())
929      {
930        if (!data_j_index.isEmpty() &&
931            data_j_index.numElements() != data_i_index.numElements())
932        {
933           ERROR("CDomain::checkCompression(void)",
934                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
935                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
936                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
937                 << "'data_j_index' size = " << data_j_index.numElements());
938        }
939
940        if (2 == data_dim)
941        {
942          if (data_j_index.isEmpty())
943          {
944             ERROR("CDomain::checkCompression(void)",
945                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
946                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
947          }
948        }
949        else // (1 == data_dim)
950        {
951          if (data_j_index.isEmpty())
952          {
953            data_j_index.resize(data_ni);
954            for (int j = 0; j < data_ni; ++j) data_j_index(j) = 0;
955          }
956        }
957      }
958      else
959      {
960        if (data_dim == 2 && !data_j_index.isEmpty())
961          ERROR("CDomain::checkCompression(void)",
962                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
963                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
964
965        if (1 == data_dim)
966        {
967          data_i_index.resize(data_ni);
968          data_j_index.resize(data_ni);
969
970          for (int i = 0; i < data_ni; ++i)
971          {
972            data_i_index(i) = i;
973            data_j_index(i) = 0;
974          }
975        }
976        else // (data_dim == 2)
977        {
978          const int dsize = data_ni * data_nj;
979          data_i_index.resize(dsize);
980          data_j_index.resize(dsize);
981
982          for(int count = 0, j = 0; j < data_nj; ++j)
983          {
984            for(int i = 0; i < data_ni; ++i, ++count)
985            {
986              data_i_index(count) = i;
987              data_j_index(count) = j;
988            }
989          }
990        }
991      }
992   }
993
994   //----------------------------------------------------------------
995   void CDomain::computeLocalMask(void)
996   {
997     localMask.resize(ni*nj) ;
998     localMask=false ;
999     size_t zoom_ibegin=global_zoom_ibegin ;
1000     size_t zoom_iend=global_zoom_ibegin+global_zoom_ni-1 ;
1001     size_t zoom_jbegin=global_zoom_jbegin ;
1002     size_t zoom_jend=global_zoom_jbegin+global_zoom_nj-1 ;
1003
1004
1005     size_t dn=data_i_index.numElements() ;
1006     int i,j ;
1007     size_t k,ind ;
1008
1009     for(k=0;k<dn;k++)
1010     {
1011       if (data_dim==2)
1012       {
1013          i=data_i_index(k)+data_ibegin ;
1014          j=data_j_index(k)+data_jbegin ;
1015       }
1016       else
1017       {
1018          i=(data_i_index(k)+data_ibegin)%ni ;
1019          j=(data_i_index(k)+data_ibegin)/ni ;
1020       }
1021
1022       if (i>=0 && i<ni && j>=0 && j<nj)
1023         if (i+ibegin>=zoom_ibegin && i+ibegin<=zoom_iend && j+jbegin>=zoom_jbegin && j+jbegin<=zoom_jend)
1024         {
1025           ind=i+ni*j ;
1026           localMask(ind)=mask_1d(ind) ;
1027         }
1028     }
1029   }
1030
1031
1032
1033
1034
1035
1036
1037   void CDomain::checkEligibilityForCompressedOutput(void)
1038   {
1039     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1040     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1041   }
1042
1043   //----------------------------------------------------------------
1044
1045   void CDomain::completeLonLatClient(void)
1046   {
1047     if (!lonvalue_2d.isEmpty())
1048     {
1049       lonvalue_client.resize(ni * nj);
1050       latvalue_client.resize(ni * nj);
1051       if (hasBounds)
1052       {
1053         bounds_lon_client.resize(nvertex, ni * nj);
1054         bounds_lat_client.resize(nvertex, ni * nj);
1055       }
1056
1057       for (int j = 0; j < nj; ++j)
1058       {
1059         for (int i = 0; i < ni; ++i)
1060         {
1061           int k = j * ni + i;
1062
1063           lonvalue_client(k) = lonvalue_2d(i,j);
1064           latvalue_client(k) = latvalue_2d(i,j);
1065
1066           if (hasBounds)
1067           {
1068             for (int n = 0; n < nvertex; ++n)
1069             {
1070               bounds_lon_client(n,k) = bounds_lon_2d(n,i,j);
1071               bounds_lat_client(n,k) = bounds_lat_2d(n,i,j);
1072             }
1073           }
1074         }
1075       }
1076     }
1077     else if (!lonvalue_1d.isEmpty())
1078     {
1079       if (type_attr::rectilinear == type)
1080       {
1081         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1082         {
1083           lonvalue_client.resize(ni * nj);
1084           latvalue_client.resize(ni * nj);
1085           if (hasBounds)
1086           {
1087             bounds_lon_client.resize(nvertex, ni * nj);
1088             bounds_lat_client.resize(nvertex, ni * nj);
1089           }
1090
1091           for (int j = 0; j < nj; ++j)
1092           {
1093             for (int i = 0; i < ni; ++i)
1094             {
1095               int k = j * ni + i;
1096
1097               lonvalue_client(k) = lonvalue_1d(i);
1098               latvalue_client(k) = latvalue_1d(j);
1099
1100               if (hasBounds)
1101               {
1102                 for (int n = 0; n < nvertex; ++n)
1103                 {
1104                   bounds_lon_client(n,k) = bounds_lon_1d(n,i);
1105                   bounds_lat_client(n,k) = bounds_lat_1d(n,j);
1106                 }
1107               }
1108             }
1109           }
1110         }
1111         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements())
1112         {
1113           lonvalue_client.reference(lonvalue_1d);
1114           latvalue_client.reference(latvalue_1d);
1115            if (hasBounds)
1116           {
1117             bounds_lon_client.reference(bounds_lon_1d);
1118             bounds_lat_client.reference(bounds_lat_1d);
1119           }
1120         }
1121         else
1122           ERROR("CDomain::completeLonClient(void)",
1123                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1124                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1125                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1126                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1127                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1128                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1129       }
1130       else if (type == type_attr::curvilinear || type == type_attr::unstructured)
1131       {
1132         lonvalue_client.reference(lonvalue_1d);
1133         latvalue_client.reference(latvalue_1d);
1134         if (hasBounds)
1135         {
1136           bounds_lon_client.reference(bounds_lon_1d);
1137           bounds_lat_client.reference(bounds_lat_1d);
1138         }
1139       }
1140     }
1141   }
1142
1143   void CDomain::checkBounds(void)
1144   {
1145     if (!nvertex.isEmpty() && nvertex > 0)
1146     {
1147       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1148         ERROR("CDomain::checkBounds(void)",
1149               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1150               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1151               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1152
1153       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1154         ERROR("CDomain::checkBounds(void)",
1155               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1156               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1157               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1158
1159       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1160       {
1161         ERROR("CDomain::checkBounds(void)",
1162               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1163               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1164               << "Please define either both attributes or none.");
1165       }
1166
1167       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1168       {
1169         ERROR("CDomain::checkBounds(void)",
1170               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1171               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1172               << "Please define either both attributes or none.");
1173       }
1174
1175       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1176         ERROR("CDomain::checkBounds(void)",
1177               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1178               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1179               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(1)
1180               << " but nvertex is " << nvertex.getValue() << ".");
1181
1182       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1183         ERROR("CDomain::checkBounds(void)",
1184               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1185               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1186               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(2)
1187               << " but nvertex is " << nvertex.getValue() << ".");
1188
1189       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1190         ERROR("CDomain::checkBounds(void)",
1191               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1192               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1193
1194       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1195         ERROR("CDomain::checkBounds(void)",
1196               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1197               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1198
1199       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1200         ERROR("CDomain::checkBounds(void)",
1201               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1202               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1203               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(1)
1204               << " but nvertex is " << nvertex.getValue() << ".");
1205
1206       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1207         ERROR("CDomain::checkBounds(void)",
1208               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1209               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1210               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(2)
1211               << " but nvertex is " << nvertex.getValue() << ".");
1212
1213       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1214         ERROR("CDomain::checkBounds(void)",
1215               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1216               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1217
1218       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1219         ERROR("CDomain::checkBounds(void)",
1220               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1221
1222       hasBounds = true;
1223     }
1224     else
1225     {
1226       hasBounds = false;
1227       nvertex = 0;
1228     }
1229   }
1230
1231   void CDomain::checkArea(void)
1232   {
1233     hasArea = !area.isEmpty();
1234     if (hasArea)
1235     {
1236       if (area.extent(0) != ni || area.extent(1) != nj)
1237       {
1238         ERROR("CDomain::checkArea(void)",
1239               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1240               << "The area does not have the same size as the local domain." << std::endl
1241               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1242               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1243       }
1244     }
1245   }
1246
1247   void CDomain::checkLonLat()
1248   {
1249     hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1250                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1251     if (hasLonLat)
1252     {
1253       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1254         ERROR("CDomain::checkLonLat()",
1255               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1256               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1257               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1258
1259       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1260       {
1261         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1262           ERROR("CDomain::checkLonLat()",
1263                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1264                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1265                 << "Local size is " << i_index.numElements() << "." << std::endl
1266                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1267       }
1268
1269       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1270       {
1271         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1272           ERROR("CDomain::checkLonLat()",
1273                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1274                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1275                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1276                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1277       }
1278
1279       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1280         ERROR("CDomain::checkLonLat()",
1281               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1282               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1283               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1284
1285       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1286       {
1287         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1288           ERROR("CDomain::checkLonLat()",
1289                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1290                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1291                 << "Local size is " << i_index.numElements() << "." << std::endl
1292                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1293       }
1294
1295       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1296       {
1297         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1298           ERROR("CDomain::checkLonLat()",
1299                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1300                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1301                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1302                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1303       }
1304     }
1305   }
1306
1307   void CDomain::checkAttributesOnClientAfterTransformation()
1308   {
1309     CContext* context=CContext::getCurrent() ;
1310
1311     if (this->isClientAfterTransformationChecked) return;
1312     if (context->hasClient)
1313     {
1314       this->checkMask();
1315       if (hasLonLat || hasArea || isCompressible_) this->computeConnectedServer();
1316       if (hasLonLat) this->completeLonLatClient();
1317     }
1318
1319     this->isClientAfterTransformationChecked = true;
1320   }
1321
1322   //----------------------------------------------------------------
1323   // Divide function checkAttributes into 2 seperate ones
1324   // This function only checks all attributes of current domain
1325   void CDomain::checkAttributesOnClient()
1326   {
1327     if (this->isClientChecked) return;
1328     CContext* context=CContext::getCurrent();
1329
1330      this->checkDomain();
1331      this->checkBounds();
1332      this->checkArea();
1333      this->checkLonLat();
1334
1335      if (context->hasClient)
1336      { // CÃŽté client uniquement
1337         this->checkMask();
1338         this->checkDomainData();
1339         this->checkCompression();
1340         this->computeLocalMask() ;
1341      }
1342      else
1343      { // CÃŽté serveur uniquement
1344      }
1345
1346      this->isClientChecked = true;
1347   }
1348
1349   // Send all checked attributes to server
1350   void CDomain::sendCheckedAttributes()
1351   {
1352     if (!this->isClientChecked) checkAttributesOnClient();
1353     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1354     CContext* context=CContext::getCurrent() ;
1355
1356     if (this->isChecked) return;
1357     if (context->hasClient)
1358     {
1359       sendServerAttribut();
1360       if (hasLonLat || hasArea || isCompressible_) sendLonLatArea();
1361     }
1362     this->isChecked = true;
1363   }
1364
1365   void CDomain::checkAttributes(void)
1366   {
1367      if (this->isChecked) return;
1368      CContext* context=CContext::getCurrent() ;
1369
1370      this->checkDomain();
1371      this->checkLonLat();
1372      this->checkBounds();
1373      this->checkArea();
1374
1375      if (context->hasClient)
1376      { // CÃŽté client uniquement
1377         this->checkMask();
1378         this->checkDomainData();
1379         this->checkCompression();
1380         this->computeLocalMask() ;
1381
1382      }
1383      else
1384      { // CÃŽté serveur uniquement
1385      }
1386
1387      if (context->hasClient)
1388      {
1389        this->computeConnectedServer();
1390        this->completeLonLatClient();
1391        this->sendServerAttribut();
1392        this->sendLonLatArea();
1393      }
1394
1395      this->isChecked = true;
1396   }
1397
1398  void CDomain::sendServerAttribut(void)
1399  {
1400    CContext* context = CContext::getCurrent();
1401    CContextClient* client = context->client;
1402    int nbServer = client->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 (client->isServerLeader())
1413    {
1414      std::list<CMessage> msgs;
1415
1416      const std::list<int>& ranks = client->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      client->sendEvent(event);
1437    }
1438    else client->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
1860  /*!
1861    Receive attributes from client(s): zoom info and begin and n of each server
1862    \param[in] rank rank of client source
1863    \param[in] buffer message containing attributes info
1864  */
1865  void CDomain::recvServerAttribut(CBufferIn& buffer)
1866  {
1867    int global_zoom_ni_tmp, global_zoom_ibegin_tmp, global_zoom_nj_tmp, global_zoom_jbegin_tmp;
1868    buffer >> ni_srv >> ibegin_srv >> iend_srv >> nj_srv >> jbegin_srv >> jend_srv
1869           >> global_zoom_ni_tmp >> global_zoom_ibegin_tmp >> global_zoom_nj_tmp >> global_zoom_jbegin_tmp
1870           >> isCompressible_;
1871
1872    global_zoom_ni.setValue(global_zoom_ni_tmp);
1873    global_zoom_ibegin.setValue(global_zoom_ibegin_tmp);
1874    global_zoom_nj.setValue(global_zoom_nj_tmp);
1875    global_zoom_jbegin.setValue(global_zoom_jbegin_tmp);
1876
1877    int zoom_iend = global_zoom_ibegin + global_zoom_ni - 1;
1878    int zoom_jend = global_zoom_jbegin + global_zoom_nj - 1;
1879
1880    zoom_ibegin_srv = global_zoom_ibegin > ibegin_srv ? global_zoom_ibegin : ibegin_srv ;
1881    zoom_iend_srv = zoom_iend < iend_srv ? zoom_iend : iend_srv ;
1882    zoom_ni_srv=zoom_iend_srv-zoom_ibegin_srv+1 ;
1883
1884    zoom_jbegin_srv = global_zoom_jbegin > jbegin_srv ? global_zoom_jbegin : jbegin_srv ;
1885    zoom_jend_srv = zoom_jend < jend_srv ? zoom_jend : jend_srv ;
1886    zoom_nj_srv=zoom_jend_srv-zoom_jbegin_srv+1 ;
1887
1888    if (zoom_ni_srv<=0 || zoom_nj_srv<=0)
1889    {
1890      zoom_ibegin_srv=0 ; zoom_iend_srv=0 ; zoom_ni_srv=0 ;
1891      zoom_jbegin_srv=0 ; zoom_jend_srv=0 ; zoom_nj_srv=0 ;
1892    }
1893    lonvalue_srv.resize(zoom_ni_srv*zoom_nj_srv) ;
1894    lonvalue_srv = 0. ;
1895    latvalue_srv.resize(zoom_ni_srv*zoom_nj_srv) ;
1896    latvalue_srv = 0. ;
1897    if (hasBounds)
1898    {
1899      bounds_lon_srv.resize(nvertex,zoom_ni_srv*zoom_nj_srv) ;
1900      bounds_lon_srv = 0. ;
1901      bounds_lat_srv.resize(nvertex,zoom_ni_srv*zoom_nj_srv) ;
1902      bounds_lat_srv = 0. ;
1903    }
1904
1905    if (hasArea)
1906    {
1907      area_srv.resize(zoom_ni_srv * zoom_nj_srv);
1908      area_srv = 0.;
1909    }
1910
1911  }
1912
1913  /*!
1914    Receive index event from clients(s)
1915    \param[in] event event contain info about rank and associated index
1916  */
1917  void CDomain::recvIndex(CEventServer& event)
1918  {
1919    CDomain* domain;
1920
1921    list<CEventServer::SSubEvent>::iterator it;
1922    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
1923    {
1924      CBufferIn* buffer = it->buffer;
1925      string domainId;
1926      *buffer >> domainId;
1927      domain = get(domainId);
1928      domain->recvIndex(it->rank, *buffer);
1929    }
1930
1931    if (domain->isCompressible_)
1932    {
1933      std::sort(domain->indexesToWrite.begin(), domain->indexesToWrite.end());
1934
1935      CContextServer* server = CContext::getCurrent()->server;
1936      domain->numberWrittenIndexes_ = domain->indexesToWrite.size();
1937      MPI_Allreduce(&domain->numberWrittenIndexes_, &domain->totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1938      MPI_Scan(&domain->numberWrittenIndexes_, &domain->offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1939      domain->offsetWrittenIndexes_ -= domain->numberWrittenIndexes_;
1940    }
1941  }
1942
1943  /*!
1944    Receive index information from client(s)
1945    \param[in] rank rank of client source
1946    \param[in] buffer message containing index info
1947  */
1948  void CDomain::recvIndex(int rank, CBufferIn& buffer)
1949  {
1950    int type_int;
1951    buffer >> type_int >> isCurvilinear >> indiSrv[rank] >> indjSrv[rank];
1952    type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
1953
1954    if (isCompressible_)
1955    {
1956      CArray<int, 1> writtenIndexes;
1957      buffer >> writtenIndexes;
1958      indexesToWrite.reserve(indexesToWrite.size() + writtenIndexes.numElements());
1959      for (int i = 0; i < writtenIndexes.numElements(); ++i)
1960        indexesToWrite.push_back(writtenIndexes(i));
1961    }
1962  }
1963
1964  /*!
1965    Receive longitude event from clients(s)
1966    \param[in] event event contain info about rank and associated longitude
1967  */
1968  void CDomain::recvLon(CEventServer& event)
1969  {
1970    list<CEventServer::SSubEvent>::iterator it;
1971    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
1972    {
1973      CBufferIn* buffer = it->buffer;
1974      string domainId;
1975      *buffer >> domainId;
1976      get(domainId)->recvLon(it->rank, *buffer);
1977    }
1978  }
1979
1980  /*!
1981    Receive longitude information from client(s)
1982    \param[in] rank rank of client source
1983    \param[in] buffer message containing longitude info
1984  */
1985  void CDomain::recvLon(int rank, CBufferIn& buffer)
1986  {
1987    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
1988    CArray<double,1> lon;
1989    CArray<double,2> boundslon;
1990
1991    buffer >> lon;
1992
1993    if (hasBounds) buffer >> boundslon;
1994
1995    int i, j, ind_srv;
1996    for (int ind = 0; ind < indi.numElements(); ind++)
1997    {
1998      i = indi(ind); j = indj(ind);
1999      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2000      lonvalue_srv(ind_srv) = lon(ind);
2001      if (hasBounds)
2002      {
2003        for (int nv = 0; nv < nvertex; ++nv)
2004          bounds_lon_srv(nv, ind_srv) = boundslon(nv, ind);
2005      }
2006    }
2007  }
2008
2009  /*!
2010    Receive latitude event from clients(s)
2011    \param[in] event event contain info about rank and associated latitude
2012  */
2013  void CDomain::recvLat(CEventServer& event)
2014  {
2015    list<CEventServer::SSubEvent>::iterator it;
2016    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2017    {
2018      CBufferIn* buffer = it->buffer;
2019      string domainId;
2020      *buffer >> domainId;
2021      get(domainId)->recvLat(it->rank, *buffer);
2022    }
2023  }
2024
2025  /*!
2026    Receive latitude information from client(s)
2027    \param[in] rank rank of client source
2028    \param[in] buffer message containing latitude info
2029  */
2030  void CDomain::recvLat(int rank, CBufferIn& buffer)
2031  {
2032    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2033    CArray<double,1> lat;
2034    CArray<double,2> boundslat;
2035
2036    buffer >> lat;
2037    if (hasBounds) buffer >> boundslat;
2038
2039    int i, j, ind_srv;
2040    for (int ind = 0; ind < indi.numElements(); ind++)
2041    {
2042      i = indi(ind); j = indj(ind);
2043      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2044      latvalue_srv(ind_srv) = lat(ind);
2045      if (hasBounds)
2046      {
2047        for (int nv = 0; nv < nvertex; nv++)
2048          bounds_lat_srv(nv, ind_srv) = boundslat(nv, ind);
2049      }
2050    }
2051  }
2052
2053  /*!
2054    Receive area event from clients(s)
2055    \param[in] event event contain info about rank and associated area
2056  */
2057  void CDomain::recvArea(CEventServer& event)
2058  {
2059    list<CEventServer::SSubEvent>::iterator it;
2060    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2061    {
2062      CBufferIn* buffer = it->buffer;
2063      string domainId;
2064      *buffer >> domainId;
2065      get(domainId)->recvArea(it->rank, *buffer);
2066    }
2067  }
2068
2069  /*!
2070    Receive area information from client(s)
2071    \param[in] rank rank of client source
2072    \param[in] buffer message containing area info
2073  */
2074  void CDomain::recvArea(int rank, CBufferIn& buffer)
2075  {
2076    CArray<int,1> &indi = indiSrv[rank], &indj = indjSrv[rank];
2077    CArray<double,1> clientArea;
2078
2079    buffer >> clientArea;
2080
2081    int i, j, ind_srv;
2082    for (int ind = 0; ind < indi.numElements(); ind++)
2083    {
2084      i = indi(ind); j = indj(ind);
2085      ind_srv = (i - zoom_ibegin_srv) + (j - zoom_jbegin_srv) * zoom_ni_srv;
2086      area_srv(ind_srv) = clientArea(ind);
2087    }
2088  }
2089
2090  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2091  {
2092    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2093    return transformationMap_.back().second;
2094  }
2095
2096  /*!
2097    Check whether a domain has transformation
2098    \return true if domain has transformation
2099  */
2100  bool CDomain::hasTransformation()
2101  {
2102    return (!transformationMap_.empty());
2103  }
2104
2105  /*!
2106    Set transformation for current domain. It's the method to move transformation in hierarchy
2107    \param [in] domTrans transformation on domain
2108  */
2109  void CDomain::setTransformations(const TransMapTypes& domTrans)
2110  {
2111    transformationMap_ = domTrans;
2112  }
2113
2114  /*!
2115    Get all transformation current domain has
2116    \return all transformation
2117  */
2118  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2119  {
2120    return transformationMap_;
2121  }
2122
2123  /*!
2124    Check the validity of all transformations applied on domain
2125  This functions is called AFTER all inherited attributes are solved
2126  */
2127  void CDomain::checkTransformations()
2128  {
2129    TransMapTypes::const_iterator itb = transformationMap_.begin(), it,
2130                                  ite = transformationMap_.end();
2131//    for (it = itb; it != ite; ++it)
2132//    {
2133//      (it->second)->checkValid(this);
2134//    }
2135  }
2136
2137  void CDomain::duplicateTransformation(CDomain* src)
2138  {
2139    if (src->hasTransformation())
2140    {
2141      this->setTransformations(src->getAllTransformations());
2142    }
2143  }
2144
2145  /*!
2146   * Go through the hierarchy to find the domain from which the transformations must be inherited
2147   */
2148  void CDomain::solveInheritanceTransformation()
2149  {
2150    if (hasTransformation() || !hasDirectDomainReference())
2151      return;
2152
2153    CDomain* domain = this;
2154    std::vector<CDomain*> refDomains;
2155    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2156    {
2157      refDomains.push_back(domain);
2158      domain = domain->getDirectDomainReference();
2159    }
2160
2161    if (domain->hasTransformation())
2162      for (size_t i = 0; i < refDomains.size(); ++i)
2163        refDomains[i]->setTransformations(domain->getAllTransformations());
2164  }
2165
2166  /*!
2167    Parse children nodes of a domain in xml file.
2168    Whenver there is a new transformation, its type and name should be added into this function
2169    \param node child node to process
2170  */
2171  void CDomain::parse(xml::CXMLNode & node)
2172  {
2173    SuperClass::parse(node);
2174
2175    if (node.goToChildElement())
2176    {
2177      StdString nodeElementName;
2178      do
2179      {
2180        StdString nodeId("");
2181        if (node.getAttributes().end() != node.getAttributes().find("id"))
2182        { nodeId = node.getAttributes()["id"]; }
2183
2184        nodeElementName = node.getElementName();
2185        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2186        it = transformationMapList_.find(nodeElementName);
2187        if (ite != it)
2188        {
2189          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2190                                                                                                                nodeId,
2191                                                                                                                &node)));
2192        }
2193        else
2194        {
2195          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2196                << "The transformation " << nodeElementName << " has not been supported yet.");
2197        }
2198      } while (node.goToNextElement()) ;
2199      node.goToParentElement();
2200    }
2201  }
2202   //----------------------------------------------------------------
2203
2204   DEFINE_REF_FUNC(Domain,domain)
2205
2206   ///---------------------------------------------------------------
2207
2208} // namespace xios
Note: See TracBrowser for help on using the repository browser.