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

Last change on this file since 971 was 971, checked in by mhnguyen, 8 years ago

Correcting the previous commit

Test
+) On Curie
+) Tests pass

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