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

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

modif for Curie, CurrContext?->CurrContext_ptr in object_factory and group_factory

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