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

Last change on this file since 975 was 975, checked in by mhnguyen, 6 years ago

Checking more strictly i_index, j_index of a domain

+) Make sure that i_index, j_index have elements

Test
+) No

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