source: XIOS2/trunk/src/node/domain.cpp @ 2440

Last change on this file since 2440 was 2440, checked in by ymipsl, 19 months ago

Add bounds management for rectilinear domain.
First try.
YM

  • Property copyright set to
    Software name : XIOS (Xml I/O Server)
    http://forge.ipsl.jussieu.fr/ioserver
    Creation date : January 2009
    Licence : CeCCIL version2
    see license file in root directory : Licence_CeCILL_V2-en.txt
    or http://www.cecill.info/licences/Licence_CeCILL_V2-en.html
    Holder : CEA/LSCE (Laboratoire des Sciences du CLimat et de l'Environnement)
    CNRS/IPSL (Institut Pierre Simon Laplace)
    Project Manager : Yann Meurdesoif
    yann.meurdesoif@cea.fr
  • Property svn:executable set to *
File size: 114.8 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
21#include <algorithm>
22
23namespace xios {
24
25   /// ////////////////////// Définitions ////////////////////// ///
26
27   CDomain::CDomain(void)
28      : CObjectTemplate<CDomain>(), CDomainAttributes()
29      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_()
30      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
31      , isClientAfterTransformationChecked(false), hasLonLat(false)
32      , isRedistributed_(false), hasPole(false)
33      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
34      , globalLocalIndexMap_(), computedWrittenIndex_(false)
35      , clients(), hasLatInReadFile_(false), hasBoundsLatInReadFile_(false)
36      , hasLonInReadFile_(false), hasBoundsLonInReadFile_(false)
37      , isTiled_(false), isTiledOnly_(false)
38   {
39   }
40
41   CDomain::CDomain(const StdString & id)
42      : CObjectTemplate<CDomain>(id), CDomainAttributes()
43      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_() 
44      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
45      , isClientAfterTransformationChecked(false), hasLonLat(false)
46      , isRedistributed_(false), hasPole(false)
47      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
48      , globalLocalIndexMap_(), computedWrittenIndex_(false)
49      , clients(), hasLatInReadFile_(false), hasBoundsLatInReadFile_(false)
50      , hasLonInReadFile_(false), hasBoundsLonInReadFile_(false)
51      , isTiled_(false), isTiledOnly_(false)
52   {
53    }
54
55   CDomain::~CDomain(void)
56   {
57   }
58
59   ///---------------------------------------------------------------
60
61   void CDomain::assignMesh(const StdString meshName, const int nvertex)
62   TRY
63   {
64     mesh = CMesh::getMesh(meshName, nvertex);
65   }
66   CATCH_DUMP_ATTR
67
68   CDomain* CDomain::createDomain()
69   TRY
70   {
71     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
72     return domain;
73   }
74   CATCH
75
76   std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
77   bool CDomain::_dummyTransformationMapList = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
78
79   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
80   TRY
81   {
82     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
83     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
84     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
85     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
86     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
87     m["reorder_domain"] = TRANS_REORDER_DOMAIN;
88     m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
89
90     return true ;
91   }
92   CATCH
93
94   const std::set<StdString> & CDomain::getRelFiles(void) const
95   TRY
96   {
97      return (this->relFiles);
98   }
99   CATCH
100
101   /*!
102     Returns the number of indexes written by each server.
103     \return the number of indexes written by each server
104   */
105   int CDomain::getNumberWrittenIndexes(MPI_Comm writtenCom)
106   TRY
107   {
108     int writtenSize;
109     MPI_Comm_size(writtenCom, &writtenSize);
110     return numberWrittenIndexes_[writtenSize];
111   }
112   CATCH_DUMP_ATTR
113
114   /*!
115     Returns the total number of indexes written by the servers.
116     \return the total number of indexes written by the servers
117   */
118   int CDomain::getTotalNumberWrittenIndexes(MPI_Comm writtenCom)
119   TRY
120   {
121     int writtenSize;
122     MPI_Comm_size(writtenCom, &writtenSize);
123     return totalNumberWrittenIndexes_[writtenSize];
124   }
125   CATCH_DUMP_ATTR
126
127   /*!
128     Returns the offset of indexes written by each server.
129     \return the offset of indexes written by each server
130   */
131   int CDomain::getOffsetWrittenIndexes(MPI_Comm writtenCom)
132   TRY
133   {
134     int writtenSize;
135     MPI_Comm_size(writtenCom, &writtenSize);
136     return offsetWrittenIndexes_[writtenSize];
137   }
138   CATCH_DUMP_ATTR
139
140   CArray<int, 1>& CDomain::getCompressedIndexToWriteOnServer(MPI_Comm writtenCom)
141   TRY
142   {
143     int writtenSize;
144     MPI_Comm_size(writtenCom, &writtenSize);
145     return compressedIndexToWriteOnServer[writtenSize];
146   }
147   CATCH_DUMP_ATTR
148
149   //----------------------------------------------------------------
150
151   /*!
152    * Compute the minimum buffer size required to send the attributes to the server(s).
153    *
154    * \return A map associating the server rank with its minimum buffer size.
155    */
156   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
157   TRY
158   {
159
160     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
161
162     if (client->isServerLeader())
163     {
164       // size estimation for sendDistributionAttribut
165       size_t size = 11 * sizeof(size_t);
166
167       const std::list<int>& ranks = client->getRanksServerLeader();
168       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
169       {
170         if (size > attributesSizes[*itRank])
171           attributesSizes[*itRank] = size;
172       }
173     }
174
175     std::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->serverSize].end();
176     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
177     for (size_t k = 0; k < connectedServerRank_[client->serverSize].size(); ++k)
178     {
179       int rank = connectedServerRank_[client->serverSize][k];
180       std::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->serverSize].find(rank);
181       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
182
183       // size estimation for sendIndex (and sendArea which is always smaller or equal)
184       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
185
186       // size estimation for sendLonLat
187       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
188       if (hasBounds)
189       {
190         int nvertexValue ;
191         if (type==type_attr::rectilinear || type==type_attr::curvilinear) nvertexValue=4 ;
192         else nvertexValue=nvertex ;
193         sizeLonLatEvent += CArray<double,2>::size(nvertexValue * idxCount);
194       }
195
196       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
197       if (size > attributesSizes[rank])
198         attributesSizes[rank] = size;
199     }
200
201     return attributesSizes;
202   }
203   CATCH_DUMP_ATTR
204
205   //----------------------------------------------------------------
206
207   bool CDomain::isEmpty(void) const
208   TRY
209   {
210     return ((this->i_index.isEmpty()) || (0 == this->i_index.numElements()));
211   }
212   CATCH
213
214   //----------------------------------------------------------------
215
216   bool CDomain::IsWritten(const StdString & filename) const
217   TRY
218   {
219      return (this->relFiles.find(filename) != this->relFiles.end());
220   }
221   CATCH
222
223   bool CDomain::isWrittenCompressed(const StdString& filename) const
224   TRY
225   {
226      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
227   }
228   CATCH
229
230   //----------------------------------------------------------------
231
232   bool CDomain::isDistributed(void) const
233   TRY
234   {
235      bool distributed =  !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
236              (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
237      bool distributed_glo ;
238      distributed |= (1 == CContext::getCurrent()->client->clientSize);
239
240      return distributed;
241   }
242   CATCH
243
244   //----------------------------------------------------------------
245
246   /*!
247    * Test whether the data defined on the domain can be outputted in a compressed way.
248    *
249    * \return true if and only if a mask was defined for this domain
250    */
251   bool CDomain::isCompressible(void) const
252   TRY
253   {
254      return isCompressible_;
255   }
256   CATCH
257
258   bool CDomain::isTiled(void) const
259   TRY
260   {
261      return isTiled_;
262   }
263   CATCH
264
265   bool CDomain::isTiledOnly(void) const
266   TRY
267   {
268      return isTiledOnly_;
269   }
270   CATCH
271
272   int CDomain::getTileId(int i, int j) const
273   TRY
274   {
275     int tileId = 0;
276     bool stop = false;
277     while (!stop)
278     {
279       if ((i >= tile_ibegin(tileId) ) && (i < tile_ni(tileId)+tile_ibegin(tileId) )
280           && (j >= tile_jbegin(tileId) ) && (j < tile_nj(tileId)+tile_jbegin(tileId) ) )
281         stop = true;
282       ++tileId;
283     }
284     return (tileId - 1);
285   }
286   CATCH
287
288   void CDomain::addRelFile(const StdString & filename)
289   TRY
290   {
291      this->relFiles.insert(filename);
292   }
293   CATCH_DUMP_ATTR
294
295   void CDomain::addRelFileCompressed(const StdString& filename)
296   TRY
297   {
298      this->relFilesCompressed.insert(filename);
299   }
300   CATCH_DUMP_ATTR
301
302   StdString CDomain::GetName(void)   { return (StdString("domain")); }
303   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
304   ENodeType CDomain::GetType(void)   { return (eDomain); }
305
306   //----------------------------------------------------------------
307
308   /*!
309      Verify if all distribution information of a domain are available
310      This checking verifies the definition of distribution attributes (ni, nj, ibegin, jbegin)
311   */
312   bool CDomain::distributionAttributesHaveValue() const
313   TRY
314   {
315      bool hasValues = true;
316
317      if (ni.isEmpty() && ibegin.isEmpty() && i_index.isEmpty())
318      {
319        hasValues = false;
320        return hasValues;
321      }
322
323      return hasValues;
324   }
325   CATCH
326
327   /*!
328     Redistribute RECTILINEAR or CURVILINEAR domain with a number of local domains.
329   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
330   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
331   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
332    \param [in] nbLocalDomain number of local domain on the domain destination
333   */
334   void CDomain::redistribute(int nbLocalDomain)
335   TRY
336   {
337     if (this->isRedistributed_) return;
338
339     this->isRedistributed_ = true;
340     CContext* context = CContext::getCurrent();
341     // For now the assumption is that secondary server pools consist of the same number of procs.
342     // CHANGE the line below if the assumption changes.
343     CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[0] : context->client;
344     int rankClient = client->clientRank;
345     int rankOnDomain = rankClient%nbLocalDomain;
346
347     if (ni_glo.isEmpty() || ni_glo <= 0 )
348     {
349        ERROR("CDomain::redistribute(int nbLocalDomain)",
350           << "[ Id = " << this->getId() << " ] "
351           << "The global domain is badly defined,"
352           << " check the \'ni_glo\'  value !")
353     }
354
355     if (nj_glo.isEmpty() || nj_glo <= 0 )
356     {
357        ERROR("CDomain::redistribute(int nbLocalDomain)",
358           << "[ Id = " << this->getId() << " ] "
359           << "The global domain is badly defined,"
360           << " check the \'nj_glo\'  value !")
361     }
362
363     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
364     {
365        int globalDomainSize = ni_glo * nj_glo;
366        if (globalDomainSize <= nbLocalDomain)
367        {
368          for (int idx = 0; idx < nbLocalDomain; ++idx)
369          {
370            if (rankOnDomain < globalDomainSize)
371            {
372              int iIdx = rankOnDomain % ni_glo;
373              int jIdx = rankOnDomain / ni_glo;
374              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
375              ni.setValue(1); nj.setValue(1);
376            }
377            else
378            {
379              ibegin.setValue(0); jbegin.setValue(0);
380              ni.setValue(0); nj.setValue(0);
381            }
382          }
383        }
384        else
385        {
386          float njGlo = nj_glo.getValue();
387          float niGlo = ni_glo.getValue();
388          int nbProcOnX, nbProcOnY, range;
389
390          // Compute (approximately) number of segment on x and y axis
391          float yOverXRatio = njGlo/niGlo;
392
393          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
394          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
395
396          // Simple distribution: Sweep from top to bottom, left to right
397          // Calculate local begin on x
398          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
399          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
400          for (int i = 1; i < nbProcOnX; ++i)
401          {
402            range = ni_glo / nbProcOnX;
403            if (i < (ni_glo%nbProcOnX)) ++range;
404            niVec[i-1] = range;
405            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
406          }
407          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
408
409          // Calculate local begin on y
410          for (int j = 1; j < nbProcOnY; ++j)
411          {
412            range = nj_glo / nbProcOnY;
413            if (j < (nj_glo%nbProcOnY)) ++range;
414            njVec[j-1] = range;
415            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
416          }
417          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
418
419          // Now assign value to ni, ibegin, nj, jbegin
420          int iIdx = rankOnDomain % nbProcOnX;
421          int jIdx = rankOnDomain / nbProcOnX;
422
423          if (rankOnDomain != (nbLocalDomain-1))
424          {
425            ibegin.setValue(ibeginVec[iIdx]);
426            jbegin.setValue(jbeginVec[jIdx]);
427            nj.setValue(njVec[jIdx]);
428            ni.setValue(niVec[iIdx]);
429          }
430          else // just merge all the remaining rectangle into the last one
431          {
432            ibegin.setValue(ibeginVec[iIdx]);
433            jbegin.setValue(jbeginVec[jIdx]);
434            nj.setValue(njVec[jIdx]);
435            ni.setValue(ni_glo - ibeginVec[iIdx]);
436          }
437        } 
438     }
439     else  // unstructured domain
440     {
441       if (this->i_index.isEmpty())
442       {
443          int globalDomainSize = ni_glo * nj_glo;
444          if (globalDomainSize <= nbLocalDomain)
445          {
446            for (int idx = 0; idx < nbLocalDomain; ++idx)
447            {
448              if (rankOnDomain < globalDomainSize)
449              {
450                int iIdx = rankOnDomain % ni_glo;
451                int jIdx = rankOnDomain / ni_glo;
452                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
453                ni.setValue(1); nj.setValue(1);
454              }
455              else
456              {
457                ibegin.setValue(0); jbegin.setValue(0);
458                ni.setValue(0); nj.setValue(0);
459              }
460            }
461          }
462          else
463          {
464            float njGlo = nj_glo.getValue();
465            float niGlo = ni_glo.getValue();
466            std::vector<int> ibeginVec(nbLocalDomain,0);
467            std::vector<int> niVec(nbLocalDomain);
468            for (int i = 1; i < nbLocalDomain; ++i)
469            {
470              int range = ni_glo / nbLocalDomain;
471              if (i < (ni_glo%nbLocalDomain)) ++range;
472              niVec[i-1] = range;
473              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
474            }
475            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
476
477            int iIdx = rankOnDomain % nbLocalDomain;
478            ibegin.setValue(ibeginVec[iIdx]);
479            jbegin.setValue(0);
480            ni.setValue(niVec[iIdx]);
481            nj.setValue(1);
482          }
483
484          i_index.resize(ni);         
485          for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
486        }
487        else
488        {
489          ibegin.setValue(this->i_index(0));
490          jbegin.setValue(0);
491          ni.setValue(this->i_index.numElements());
492          nj.setValue(1);
493        }
494     }
495
496     checkDomain();
497   }
498   CATCH_DUMP_ATTR
499
500   /*!
501     Fill in longitude and latitude whose values are read from file
502   */
503   void CDomain::fillInLonLat()
504   TRY
505   {
506     switch (type)
507     {
508      case type_attr::rectilinear:
509        fillInRectilinearLonLat();
510        break;
511      case type_attr::curvilinear:
512        fillInCurvilinearLonLat();
513        break;
514      case type_attr::unstructured:
515        fillInUnstructuredLonLat();
516        break;
517
518      default:
519      break;
520     }
521     completeLonLatClient() ;
522   }
523   CATCH_DUMP_ATTR
524
525   /*!
526     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
527     Range of longitude value from 0 - 360
528     Range of latitude value from -90 - +90
529   */
530   void CDomain::fillInRectilinearLonLat()
531   TRY
532   {
533     if (!lonvalue_rectilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
534     {
535       lonvalue_1d.resize(ni);
536       for (int idx = 0; idx < ni; ++idx)
537         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
538       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
539       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
540     }
541     else if (!hasLonInReadFile_)
542     {
543       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
544       lonvalue_1d.resize(ni);
545       double lonRange = lon_end - lon_start;
546       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
547
548        // Assign lon value
549       for (int i = 0; i < ni; ++i)
550       {
551         if (0 == (ibegin + i))
552         {
553           lonvalue_1d(i) = lon_start;
554         }
555         else if (ni_glo == (ibegin + i + 1))
556         {
557           lonvalue_1d(i) = lon_end;
558         }
559         else
560         {
561           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
562         }
563       }
564     }
565
566
567     if (!latvalue_rectilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
568     {
569       latvalue_1d.resize(nj);
570       for (int idx = 0; idx < nj; ++idx)
571         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
572       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
573       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
574     }
575     else if (!hasLatInReadFile_)
576     {
577       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
578       latvalue_1d.resize(nj);
579
580       double latRange = lat_end - lat_start;
581       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
582
583       for (int j = 0; j < nj; ++j)
584       {
585         if (0 == (jbegin + j))
586         {
587            latvalue_1d(j) = lat_start;
588         }
589         else if (nj_glo == (jbegin + j + 1))
590         {
591            latvalue_1d(j) = lat_end;
592         }
593         else
594         {
595           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
596         }
597       }
598     }
599   }
600   CATCH_DUMP_ATTR
601
602    /*
603      Fill in 2D longitude and latitude of curvilinear domain read from a file.
604      If there are already longitude and latitude defined by model. We just ignore read value.
605    */
606   void CDomain::fillInCurvilinearLonLat()
607   TRY
608   {
609     if (!lonvalue_curvilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
610     {
611       lonvalue_2d.resize(ni,nj);
612       for (int jdx = 0; jdx < nj; ++jdx)
613        for (int idx = 0; idx < ni; ++idx)
614         lonvalue_2d(idx,jdx) = lonvalue_curvilinear_read_from_file(idx, jdx);
615
616       lonvalue_curvilinear_read_from_file.free();
617     }
618
619     if (!latvalue_curvilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
620     {
621       latvalue_2d.resize(ni,nj);
622       for (int jdx = 0; jdx < nj; ++jdx)
623        for (int idx = 0; idx < ni; ++idx)
624           latvalue_2d(idx,jdx) = latvalue_curvilinear_read_from_file(idx, jdx);
625
626       latvalue_curvilinear_read_from_file.free();
627     }
628
629     if (!bounds_lonvalue_curvilinear_read_from_file.isEmpty() && bounds_lon_2d.isEmpty() && bounds_lon_1d.isEmpty())
630     {
631       bounds_lon_2d.resize(nvertex,ni,nj);
632       for (int jdx = 0; jdx < nj; ++jdx)
633        for (int idx = 0; idx < ni; ++idx)
634          for (int ndx = 0; ndx < nvertex; ++ndx)
635            bounds_lon_2d(ndx,idx,jdx) = bounds_lonvalue_curvilinear_read_from_file(ndx,idx, jdx);
636
637       bounds_lonvalue_curvilinear_read_from_file.free();
638     }
639
640     if (!bounds_latvalue_curvilinear_read_from_file.isEmpty() && bounds_lat_2d.isEmpty() && bounds_lat_1d.isEmpty())
641     {
642       bounds_lat_2d.resize(nvertex,ni,nj);
643       for (int jdx = 0; jdx < nj; ++jdx)
644        for (int idx = 0; idx < ni; ++idx)
645          for (int ndx = 0; ndx < nvertex; ++ndx)
646            bounds_lat_2d(ndx,idx,jdx) = bounds_latvalue_curvilinear_read_from_file(ndx,idx, jdx);
647
648       bounds_latvalue_curvilinear_read_from_file.free();
649     }
650   }
651   CATCH_DUMP_ATTR
652
653    /*
654      Fill in longitude and latitude of unstructured domain read from a file
655      If there are already longitude and latitude defined by model. We just igonore reading value.
656    */
657   void CDomain::fillInUnstructuredLonLat()
658   TRY
659   {
660     if (i_index.isEmpty())
661     {
662       i_index.resize(ni);
663       for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
664     }
665
666     if (!lonvalue_unstructured_read_from_file.isEmpty() && lonvalue_1d.isEmpty())
667     {
668        lonvalue_1d.resize(ni);
669        for (int idx = 0; idx < ni; ++idx)
670          lonvalue_1d(idx) = lonvalue_unstructured_read_from_file(idx);
671
672        // We dont need these values anymore, so just delete them
673        lonvalue_unstructured_read_from_file.free();
674     } 
675
676     if (!latvalue_unstructured_read_from_file.isEmpty() && latvalue_1d.isEmpty())
677     {
678        latvalue_1d.resize(ni);
679        for (int idx = 0; idx < ni; ++idx)
680          latvalue_1d(idx) =  latvalue_unstructured_read_from_file(idx);
681
682        // We dont need these values anymore, so just delete them
683        latvalue_unstructured_read_from_file.free();
684     }
685
686     if (!bounds_lonvalue_unstructured_read_from_file.isEmpty() && bounds_lon_1d.isEmpty())
687     {
688        int nbVertex = nvertex;
689        bounds_lon_1d.resize(nbVertex,ni);
690        for (int idx = 0; idx < ni; ++idx)
691          for (int jdx = 0; jdx < nbVertex; ++jdx)
692            bounds_lon_1d(jdx,idx) = bounds_lonvalue_unstructured_read_from_file(jdx, idx);
693
694        // We dont need these values anymore, so just delete them
695        bounds_lonvalue_unstructured_read_from_file.free();
696     }
697
698     if (!bounds_latvalue_unstructured_read_from_file.isEmpty() && bounds_lat_1d.isEmpty())
699     {
700        int nbVertex = nvertex;
701        bounds_lat_1d.resize(nbVertex,ni);
702        for (int idx = 0; idx < ni; ++idx)
703          for (int jdx = 0; jdx < nbVertex; ++jdx)
704            bounds_lat_1d(jdx,idx) = bounds_latvalue_unstructured_read_from_file(jdx, idx);
705
706        // We dont need these values anymore, so just delete them
707        bounds_latvalue_unstructured_read_from_file.free();
708     }
709   }
710   CATCH_DUMP_ATTR
711
712  /*
713    Get global longitude and latitude of rectilinear domain.
714  */
715   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g)
716   TRY
717   {
718     CContext* context = CContext::getCurrent();
719    // For now the assumption is that secondary server pools consist of the same number of procs.
720    // CHANGE the line below if the assumption changes.
721    CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[0] : context->client;
722     lon_g.resize(ni_glo) ;
723     lat_g.resize(nj_glo) ;
724
725
726     int* ibegin_g = new int[client->clientSize] ;
727     int* jbegin_g = new int[client->clientSize] ;
728     int* ni_g = new int[client->clientSize] ;
729     int* nj_g = new int[client->clientSize] ;
730     int v ;
731     v=ibegin ;
732     MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,client->intraComm) ;
733     v=jbegin ;
734     MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,client->intraComm) ;
735     v=ni ;
736     MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,client->intraComm) ;
737     v=nj ;
738     MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,client->intraComm) ;
739
740     MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,client->intraComm) ;
741     MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,client->intraComm) ;
742
743      delete[] ibegin_g ;
744      delete[] jbegin_g ;
745      delete[] ni_g ;
746      delete[] nj_g ;
747   }
748   CATCH_DUMP_ATTR
749
750   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
751                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
752   TRY
753  {
754     int i,j,k;
755
756     const int nvertexValue = 4;
757     boundsLon.resize(nvertexValue,ni*nj);
758
759     if (ni_glo>1)
760     {
761       double lonStepStart = lon(1)-lon(0);
762       bounds_lon_start=lon(0) - lonStepStart/2;
763       double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2);
764       bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
765       double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
766
767       // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
768       if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
769       {
770         bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
771         bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
772       }
773     }
774     else
775     {
776       if (bounds_lon_start.isEmpty()) bounds_lon_start=-180. ;
777       if (bounds_lon_end.isEmpty()) bounds_lon_end=180.-1e-8 ;
778     }
779
780     for(j=0;j<nj;++j)
781       for(i=0;i<ni;++i)
782       {
783         k=j*ni+i;
784         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start
785                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2;
786         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end
787                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2;
788       }
789
790
791    boundsLat.resize(nvertexValue,nj*ni);
792    bool isNorthPole=false ;
793    bool isSouthPole=false ;
794    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
795    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
796
797    // lat boundaries beyond pole the assimilate it to pole
798    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
799    if (nj_glo>1)
800    {
801      double latStepStart = lat(1)-lat(0);
802      if (isNorthPole) bounds_lat_start=lat(0);
803      else
804      {
805        bounds_lat_start=lat(0)-latStepStart/2;
806        if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
807        else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
808        else if (bounds_lat_start <= 90 && bounds_lat_start >= lat(0))
809        {
810          if ( std::abs(90-bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
811        }
812        else if (bounds_lat_start >= -90 && bounds_lat_start <= lat(0))
813        {
814          if ( std::abs(-90 - bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
815        }
816      }
817
818      double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
819      if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
820      else
821      {
822        bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
823
824        if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
825        else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
826        else if (bounds_lat_end <= 90 && bounds_lat_end >= lat(nj_glo-1))
827        {
828          if ( std::abs(90-bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
829        }
830        else if (bounds_lat_end >= -90 && bounds_lat_end <= lat(nj_glo-1))
831        {
832          if ( std::abs(-90 - bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
833        }
834      }
835    }
836    else
837    {
838      if (bounds_lat_start.isEmpty()) bounds_lat_start=-90. ;
839      if (bounds_lat_end.isEmpty()) bounds_lat_end=90 ;
840    }
841
842    for(j=0;j<nj;++j)
843      for(i=0;i<ni;++i)
844      {
845        k=j*ni+i;
846        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start
847                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2;
848        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end
849                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2;
850      }
851   }
852   CATCH_DUMP_ATTR
853
854   /*
855     General check of the domain to verify its mandatory attributes
856   */
857   void CDomain::checkDomain(void)
858   TRY
859   {
860     if (type.isEmpty())
861     {
862       ERROR("CDomain::checkDomain(void)",
863             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
864             << "The domain type is mandatory, "
865             << "please define the 'type' attribute.")
866     }
867
868     if (type == type_attr::gaussian) 
869     {
870        hasPole=true ;
871        type.setValue(type_attr::unstructured) ;
872      }
873      else if (type == type_attr::rectilinear) hasPole=true ;
874
875     if (type == type_attr::unstructured)
876     {
877        if (ni_glo.isEmpty())
878        {
879          ERROR("CDomain::checkDomain(void)",
880                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
881                << "The global domain is badly defined, "
882                << "the mandatory 'ni_glo' attribute is missing.")
883        }
884        else if (ni_glo <= 0)
885        {
886          ERROR("CDomain::checkDomain(void)",
887                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
888                << "The global domain is badly defined, "
889                << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
890        }
891        isUnstructed_ = true;
892        nj_glo = 1;
893        nj = 1;
894        jbegin = 0;
895        if (!i_index.isEmpty()) ni = i_index.numElements();
896        j_index.resize(ni);
897        for(int i=0;i<ni;++i) j_index(i)=0;
898
899        if (!area.isEmpty())
900        {
901          CArray<double,2> newArea(area.dataFirst(),shape(ni,nj),duplicateData) ;
902          area.reference(newArea) ;
903        }
904     }
905
906     if (ni_glo.isEmpty())
907     {
908       ERROR("CDomain::checkDomain(void)",
909             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
910             << "The global domain is badly defined, "
911             << "the mandatory 'ni_glo' attribute is missing.")
912     }
913     else if (ni_glo <= 0)
914     {
915       ERROR("CDomain::checkDomain(void)",
916             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
917             << "The global domain is badly defined, "
918             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
919     }
920
921     if (nj_glo.isEmpty())
922     {
923       ERROR("CDomain::checkDomain(void)",
924             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
925             << "The global domain is badly defined, "
926             << "the mandatory 'nj_glo' attribute is missing.")
927     }
928     else if (nj_glo <= 0)
929     {
930       ERROR("CDomain::checkDomain(void)",
931             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
932             << "The global domain is badly defined, "
933             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
934     }
935
936     checkLocalIDomain();
937     checkLocalJDomain();
938
939     if (i_index.isEmpty())
940     {
941       i_index.resize(ni*nj);
942       for (int j = 0; j < nj; ++j)
943         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
944     }
945
946     if (j_index.isEmpty())
947     {
948       j_index.resize(ni*nj);
949       for (int j = 0; j < nj; ++j)
950         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
951     }
952   }
953   CATCH_DUMP_ATTR
954
955   size_t CDomain::getGlobalWrittenSize(void)
956   {
957     return ni_glo*nj_glo ;
958   }
959   //----------------------------------------------------------------
960
961   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
962   void CDomain::checkLocalIDomain(void)
963   TRY
964   {
965      // If ibegin and ni are provided then we use them to check the validity of local domain
966      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
967      {
968        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
969        {
970          ERROR("CDomain::checkLocalIDomain(void)",
971                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
972                << "The local domain is wrongly defined,"
973                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
974        }
975      }
976
977      // i_index has higher priority than ibegin and ni
978      if (!i_index.isEmpty())
979      {
980        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
981        if (ni.isEmpty()) 
982        {         
983         // No information about ni
984          int minIndex = ni_glo - 1;
985          int maxIndex = 0;
986          for (int idx = 0; idx < i_index.numElements(); ++idx)
987          {
988            if (i_index(idx) < minIndex) minIndex = i_index(idx);
989            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
990          }
991          ni = maxIndex - minIndex + 1; 
992          minIIndex = minIIndex;         
993        }
994
995        // It's not so correct but if ibegin is not the first value of i_index
996        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
997        if (ibegin.isEmpty()) ibegin = minIIndex;
998      }
999      else if (ibegin.isEmpty() && ni.isEmpty())
1000      {
1001        ibegin = 0;
1002        ni = ni_glo;
1003      }
1004      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
1005      {
1006        ERROR("CDomain::checkLocalIDomain(void)",
1007              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1008              << "The local domain is wrongly defined," << endl
1009              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
1010              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
1011      }
1012       
1013
1014      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
1015      {
1016        ERROR("CDomain::checkLocalIDomain(void)",
1017              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1018              << "The local domain is wrongly defined,"
1019              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
1020      }
1021   }
1022   CATCH_DUMP_ATTR
1023
1024   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
1025   void CDomain::checkLocalJDomain(void)
1026   TRY
1027   {
1028    // If jbegin and nj are provided then we use them to check the validity of local domain
1029     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
1030     {
1031       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
1032       {
1033         ERROR("CDomain::checkLocalJDomain(void)",
1034                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1035                << "The local domain is wrongly defined,"
1036                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1037       }
1038     }
1039
1040     if (!j_index.isEmpty())
1041     {
1042        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
1043        if (nj.isEmpty()) 
1044        {
1045          // No information about nj
1046          int minIndex = nj_glo - 1;
1047          int maxIndex = 0;
1048          for (int idx = 0; idx < j_index.numElements(); ++idx)
1049          {
1050            if (j_index(idx) < minIndex) minIndex = j_index(idx);
1051            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
1052          }
1053          nj = maxIndex - minIndex + 1;
1054          minJIndex = minIndex; 
1055        } 
1056        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
1057        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
1058       if (jbegin.isEmpty()) jbegin = minJIndex;       
1059     }
1060     else if (jbegin.isEmpty() && nj.isEmpty())
1061     {
1062       jbegin = 0;
1063       nj = nj_glo;
1064     }     
1065
1066
1067     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
1068     {
1069       ERROR("CDomain::checkLocalJDomain(void)",
1070              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1071              << "The local domain is wrongly defined,"
1072              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1073     }
1074   }
1075   CATCH_DUMP_ATTR
1076
1077   //----------------------------------------------------------------
1078
1079   void CDomain::checkMask(void)
1080   TRY
1081   {
1082      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
1083        ERROR("CDomain::checkMask(void)",
1084              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1085              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
1086              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
1087
1088      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
1089      {
1090        if (mask_1d.numElements() != i_index.numElements())
1091          ERROR("CDomain::checkMask(void)",
1092                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1093                << "'mask_1d' does not have the same size as the local domain." << std::endl
1094                << "Local size is " << i_index.numElements() << "." << std::endl
1095                << "Mask size is " << mask_1d.numElements() << ".");
1096      }
1097
1098      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
1099      {
1100        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1101          ERROR("CDomain::checkMask(void)",
1102                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1103                << "The mask does not have the same size as the local domain." << std::endl
1104                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1105                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
1106      }
1107
1108      if (!mask_2d.isEmpty())
1109      {
1110        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
1111        for (int j = 0; j < nj; ++j)
1112          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1113//        mask_2d.reset();
1114      }
1115      else if (mask_1d.isEmpty())
1116      {
1117        domainMask.resize(i_index.numElements());
1118        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
1119      }
1120      else
1121      {
1122      domainMask.resize(mask_1d.numElements());
1123      domainMask=mask_1d ;
1124     }
1125   }
1126   CATCH_DUMP_ATTR
1127
1128   //----------------------------------------------------------------
1129
1130   void CDomain::checkDomainData(void)
1131   TRY
1132   {
1133      if (data_dim.isEmpty())
1134      {
1135        data_dim.setValue(1);
1136      }
1137      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
1138      {
1139        ERROR("CDomain::checkDomainData(void)",
1140              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1141              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
1142      }
1143
1144      if (data_ibegin.isEmpty())
1145         data_ibegin.setValue(0);
1146      if (data_jbegin.isEmpty())
1147         data_jbegin.setValue(0);
1148
1149      if (data_ni.isEmpty())
1150      {
1151        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
1152      }
1153      else if (data_ni.getValue() < 0)
1154      {
1155        ERROR("CDomain::checkDomainData(void)",
1156              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1157              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
1158      }
1159
1160      if (data_nj.isEmpty())
1161      {
1162        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
1163      }
1164      else if (data_nj.getValue() < 0)
1165      {
1166        ERROR("CDomain::checkDomainData(void)",
1167              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1168              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1169      }
1170   }
1171   CATCH_DUMP_ATTR
1172
1173   //----------------------------------------------------------------
1174
1175   void CDomain::checkCompression(void)
1176   TRY
1177   {
1178     int i,j,ind;
1179      if (!data_i_index.isEmpty())
1180      {
1181        if (!data_j_index.isEmpty() &&
1182            data_j_index.numElements() != data_i_index.numElements())
1183        {
1184           ERROR("CDomain::checkCompression(void)",
1185                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1186                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1187                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1188                 << "'data_j_index' size = " << data_j_index.numElements());
1189        }
1190
1191        if (2 == data_dim)
1192        {
1193          if (data_j_index.isEmpty())
1194          {
1195             ERROR("CDomain::checkCompression(void)",
1196                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1197                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
1198          }
1199          for (int k=0; k<data_i_index.numElements(); ++k)
1200          {
1201            i = data_i_index(k)+data_ibegin ;
1202            j = data_j_index(k)+data_jbegin ;
1203            if (i>=0 && i<ni && j>=0 && j<nj)
1204            {
1205              ind=j*ni+i ;
1206              if (!domainMask(ind))
1207              {
1208                data_i_index(k) = -1;
1209                data_j_index(k) = -1;
1210              }
1211            }
1212            else
1213            {
1214              data_i_index(k) = -1;
1215              data_j_index(k) = -1;
1216            }
1217          }
1218        }
1219        else // (1 == data_dim)
1220        {
1221          if (data_j_index.isEmpty())
1222          {
1223            data_j_index.resize(data_ni);
1224            data_j_index = 0;
1225          }
1226          for (int k=0; k<data_i_index.numElements(); ++k)
1227          {
1228            i=data_i_index(k)+data_ibegin ;
1229            if (i>=0 && i < domainMask.size())
1230            {
1231              if (!domainMask(i)) data_i_index(k) = -1;
1232            }
1233            else
1234              data_i_index(k) = -1;
1235
1236            if (!domainMask(i)) data_i_index(k) = -1;
1237          }
1238        }
1239      }
1240      else
1241      {
1242        if (data_dim == 2 && !data_j_index.isEmpty())
1243          ERROR("CDomain::checkCompression(void)",
1244                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1245                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
1246
1247        if (1 == data_dim)
1248        {
1249          data_i_index.resize(data_ni);
1250          data_j_index.resize(data_ni);
1251          data_j_index = 0;
1252
1253          for (int k = 0; k < data_ni; ++k)
1254          {
1255            i=k+data_ibegin ;
1256            if (i>=0 && i < domainMask.size())
1257            {
1258              if (domainMask(i))
1259                data_i_index(k) = k;
1260              else
1261                data_i_index(k) = -1;
1262            }
1263            else
1264              data_i_index(k) = -1;
1265          }
1266        }
1267        else // (data_dim == 2)
1268        {
1269          const int dsize = data_ni * data_nj;
1270          data_i_index.resize(dsize);
1271          data_j_index.resize(dsize);
1272
1273          for(int count = 0, kj = 0; kj < data_nj; ++kj)
1274          {
1275            for(int ki = 0; ki < data_ni; ++ki, ++count)
1276            {
1277              i = ki + data_ibegin;
1278              j = kj + data_jbegin;
1279              ind=j*ni+i ;
1280              if (i>=0 && i<ni && j>=0 && j<nj)
1281              {
1282                if (domainMask(ind))
1283                {
1284                  data_i_index(count) = ki;
1285                  data_j_index(count) = kj;
1286                }
1287                else
1288                {
1289                  data_i_index(count) = -1;
1290                  data_j_index(count) = -1;
1291                }
1292              }
1293              else
1294              {
1295                data_i_index(count) = -1;
1296                data_j_index(count) = -1;
1297              }
1298            }
1299          }
1300        }
1301      }
1302   }
1303   CATCH_DUMP_ATTR
1304
1305   //----------------------------------------------------------------
1306
1307   /*!
1308    * For tiled domains, data_i/j_index should not take into
1309    * account parameters defining data (data_ni/nj, data_i/jbegin...)
1310    * \param [out] dataIndexI
1311    * \param [out] dataIndexJ
1312    * \param [out] infoIndexI
1313    * \param [out] infoIndexJ
1314    */
1315
1316   void CDomain::computeCompressionTiled(CArray<int,1>& dataIndexI, CArray<int,1>& dataIndexJ,
1317                                         CArray<int,1>& infoIndexI, CArray<int,1>& infoIndexJ)
1318   TRY
1319   {
1320     const int dsize = ni * nj;
1321     dataIndexI.resize(dsize);
1322     dataIndexJ.resize(dsize);
1323
1324     dataIndexJ = 0;
1325     for (int k = 0; k < ni * nj; ++k)
1326       dataIndexI(k) = k;
1327
1328     infoIndexI.resize(ni*nj);
1329     for (int j = 0; j < nj; ++j)
1330       for (int i = 0; i < ni; ++i) infoIndexI(i+j*ni) = i+ibegin;
1331
1332     infoIndexJ.resize(ni*nj);
1333     for (int j = 0; j < nj; ++j)
1334       for (int i = 0; i < ni; ++i) infoIndexJ(i+j*ni) = j+jbegin;
1335   }
1336   CATCH_DUMP_ATTR
1337
1338   //----------------------------------------------------------------
1339   void CDomain::computeLocalMask(void)
1340   TRY
1341   {
1342     localMask.resize(i_index.numElements()) ;
1343     localMask=false ;
1344
1345     size_t dn=data_i_index.numElements() ;
1346     int i,j ;
1347     size_t k,ind ;
1348
1349     for(k=0;k<dn;k++)
1350     {
1351       if (data_dim==2)
1352       {
1353          i=data_i_index(k)+data_ibegin ;
1354          j=data_j_index(k)+data_jbegin ;
1355          if (i>=0 && i<ni && j>=0 && j<nj)
1356          {
1357            ind=j*ni+i ;
1358            localMask(ind)=domainMask(ind) ;
1359          }
1360       }
1361       else
1362       {
1363          i=data_i_index(k)+data_ibegin ;
1364          if (i>=0 && i<i_index.numElements())
1365          {
1366            ind=i ;
1367            localMask(ind)=domainMask(ind) ;
1368          }
1369       }
1370     }
1371   }
1372   CATCH_DUMP_ATTR
1373
1374   void CDomain::checkEligibilityForCompressedOutput(void)
1375   TRY
1376   {
1377     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1378     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1379   }
1380   CATCH_DUMP_ATTR
1381
1382   //----------------------------------------------------------------
1383
1384   /*
1385     Fill in longitude, latitude, bounds, and area into internal values (lonvalue, latvalue, bounds_lonvalue, bounds_latvalue, areavalue)
1386     which will be used by XIOS.
1387   */
1388   void CDomain::completeLonLatClient(void)
1389   TRY
1390   {
1391     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1392     checkBounds() ;
1393     checkArea() ;
1394
1395     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
1396     {
1397       lonvalue.resize(ni * nj);
1398       latvalue.resize(ni * nj);
1399       if (hasBounds)
1400       {
1401         bounds_lonvalue.resize(nvertex, ni * nj);
1402         bounds_latvalue.resize(nvertex, ni * nj);
1403       }
1404
1405       for (int j = 0; j < nj; ++j)
1406       {
1407         for (int i = 0; i < ni; ++i)
1408         {
1409           int k = j * ni + i;
1410
1411           lonvalue(k) = lonvalue_2d(i,j);
1412           latvalue(k) = latvalue_2d(i,j);
1413
1414           if (hasBounds)
1415           {
1416             for (int n = 0; n < nvertex; ++n)
1417             {
1418               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1419               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
1420             }
1421           }
1422         }
1423       }
1424     }
1425     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
1426     {
1427       if (type_attr::rectilinear == type)
1428       {
1429         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1430         {
1431           lonvalue.resize(ni * nj);
1432           latvalue.resize(ni * nj);
1433           if (hasBounds)
1434           {
1435             bounds_lonvalue.resize(4, ni * nj);
1436             bounds_latvalue.resize(4, ni * nj);
1437           }
1438
1439           for (int j = 0; j < nj; ++j)
1440           {
1441             for (int i = 0; i < ni; ++i)
1442             {
1443               int k = j * ni + i;
1444
1445               lonvalue(k) = lonvalue_1d(i);
1446               latvalue(k) = latvalue_1d(j);
1447
1448               if (hasBounds)
1449               {
1450                  bounds_lonvalue(0,k) = bounds_lon_1d(0,i);
1451                  bounds_latvalue(0,k) = bounds_lat_1d(0,j);
1452                  bounds_lonvalue(1,k) = bounds_lon_1d(1,i);
1453                  bounds_latvalue(1,k) = bounds_lat_1d(0,j);
1454                  bounds_lonvalue(2,k) = bounds_lon_1d(1,i);
1455                  bounds_latvalue(2,k) = bounds_lat_1d(1,j);
1456                  bounds_lonvalue(3,k) = bounds_lon_1d(0,i);
1457                  bounds_latvalue(3,k) = bounds_lat_1d(1,j);
1458               }
1459             }
1460           }
1461         }
1462         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1463         {
1464           lonvalue.reference(lonvalue_1d);
1465           latvalue.reference(latvalue_1d);
1466            if (hasBounds)
1467           {
1468             bounds_lonvalue.reference(bounds_lon_1d);
1469             bounds_latvalue.reference(bounds_lat_1d);
1470           }
1471         }
1472         else
1473           ERROR("CDomain::completeLonClient(void)",
1474                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1475                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1476                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1477                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1478                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1479                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1480       }
1481       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1482       {
1483         lonvalue.reference(lonvalue_1d);
1484         latvalue.reference(latvalue_1d);
1485         if (hasBounds)
1486         {
1487           bounds_lonvalue.reference(bounds_lon_1d);
1488           bounds_latvalue.reference(bounds_lat_1d);
1489         }
1490       }
1491     }
1492
1493     if (!area.isEmpty() && areavalue.isEmpty())
1494     {
1495        areavalue.resize(ni*nj);
1496       for (int j = 0; j < nj; ++j)
1497       {
1498         for (int i = 0; i < ni; ++i)
1499         {
1500           int k = j * ni + i;
1501           areavalue(k) = area(i,j);
1502         }
1503       }
1504     }
1505   }
1506   CATCH_DUMP_ATTR
1507
1508   /*
1509     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1510   */
1511   void CDomain::convertLonLatValue(void)
1512   TRY
1513   {
1514     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1515     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1516     {
1517       lonvalue_2d.resize(ni,nj);
1518       latvalue_2d.resize(ni,nj);
1519       if (hasBounds)
1520       {
1521         bounds_lon_2d.resize(nvertex, ni, nj);
1522         bounds_lat_2d.resize(nvertex, ni, nj);
1523       }
1524
1525       for (int j = 0; j < nj; ++j)
1526       {
1527         for (int i = 0; i < ni; ++i)
1528         {
1529           int k = j * ni + i;
1530
1531           lonvalue_2d(i,j) = lonvalue(k);
1532           latvalue_2d(i,j) = latvalue(k);
1533
1534           if (hasBounds)
1535           {
1536             for (int n = 0; n < nvertex; ++n)
1537             {
1538               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1539               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1540             }
1541           }
1542         }
1543       }
1544     }
1545     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1546     {
1547       if (type_attr::rectilinear == type)
1548       {
1549         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1550         {
1551           lonvalue.resize(ni * nj);
1552           latvalue.resize(ni * nj);
1553           if (hasBounds)
1554           {
1555             bounds_lonvalue.resize(4, ni * nj);
1556             bounds_latvalue.resize(4, ni * nj);
1557           }
1558
1559           for (int j = 0; j < nj; ++j)
1560           {
1561             for (int i = 0; i < ni; ++i)
1562             {
1563               int k = j * ni + i;
1564
1565               lonvalue(k) = lonvalue_1d(i);
1566               latvalue(k) = latvalue_1d(j);
1567
1568               if (hasBounds)
1569               {
1570                  bounds_lonvalue(0,k) = bounds_lon_1d(0,i);
1571                  bounds_latvalue(0,k) = bounds_lat_1d(0,j);
1572                  bounds_lonvalue(1,k) = bounds_lon_1d(1,i);
1573                  bounds_latvalue(1,k) = bounds_lat_1d(0,j);
1574                  bounds_lonvalue(2,k) = bounds_lon_1d(1,i);
1575                  bounds_latvalue(2,k) = bounds_lat_1d(1,j);
1576                  bounds_lonvalue(3,k) = bounds_lon_1d(0,i);
1577                  bounds_latvalue(3,k) = bounds_lat_1d(1,j);
1578               }
1579             }
1580           }
1581         }
1582         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1583         {
1584           lonvalue.reference(lonvalue_1d);
1585           latvalue.reference(latvalue_1d);
1586            if (hasBounds)
1587           {
1588             bounds_lonvalue.reference(bounds_lon_1d);
1589             bounds_latvalue.reference(bounds_lat_1d);
1590           }
1591         }
1592         else
1593           ERROR("CDomain::completeLonClient(void)",
1594                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1595                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1596                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1597                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1598                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1599                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1600       }
1601       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1602       {
1603         lonvalue.reference(lonvalue_1d);
1604         latvalue.reference(latvalue_1d);
1605         if (hasBounds)
1606         {
1607           bounds_lonvalue.reference(bounds_lon_1d);
1608           bounds_latvalue.reference(bounds_lat_1d);
1609         }
1610       }
1611     }
1612   }
1613   CATCH_DUMP_ATTR
1614
1615   void CDomain::checkBounds(void)
1616   TRY
1617   {
1618     int nvertexValue ; 
1619     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1620     if (hasBoundValues) 
1621     {
1622       hasBounds=true ;
1623       return ;
1624      }
1625
1626     if (type_attr::rectilinear == type) nvertexValue=2 ;
1627     else if (type_attr::curvilinear == type) nvertexValue=4 ;
1628     else if (!nvertex.isEmpty() && nvertex>0) nvertexValue = nvertex ;
1629     else 
1630     {
1631       hasBounds=false ;
1632       return ;
1633     }
1634   
1635      if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1636        ERROR("CDomain::checkBounds(void)",
1637              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1638              << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1639              << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1640
1641      if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1642        ERROR("CDomain::checkBounds(void)",
1643              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1644              << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1645              << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1646
1647      if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1648      {
1649        ERROR("CDomain::checkBounds(void)",
1650              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1651              << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1652              << "Please define either both attributes or none.");
1653      }
1654
1655      if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1656      {
1657        ERROR("CDomain::checkBounds(void)",
1658              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1659              << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1660              << "Please define either both attributes or none.");
1661      }
1662
1663      if (!bounds_lon_1d.isEmpty() && nvertexValue != bounds_lon_1d.extent(0))
1664        ERROR("CDomain::checkBounds(void)",
1665              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1666              << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1667              << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(0)
1668              << " but nvertex is " << nvertexValue << ".");
1669
1670      if (!bounds_lon_2d.isEmpty() && nvertexValue != bounds_lon_2d.extent(0))
1671        ERROR("CDomain::checkBounds(void)",
1672              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1673              << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1674              << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(0)
1675              << " but nvertex is " << nvertexValue << ".");
1676
1677      if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1678        ERROR("CDomain::checkBounds(void)",
1679              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1680              << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1681
1682      if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1683        ERROR("CDomain::checkBounds(void)",
1684              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1685              << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1686
1687      if (!bounds_lat_1d.isEmpty() && nvertexValue != bounds_lat_1d.extent(0))
1688        ERROR("CDomain::checkBounds(void)",
1689              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1690              << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1691              << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(0)
1692              << " but nvertex is " << nvertexValue << ".");
1693
1694      if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1695        ERROR("CDomain::checkBounds(void)",
1696              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1697              << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1698              << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(0)
1699              << " but nvertex is " << nvertexValue << ".");
1700
1701      if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1702        ERROR("CDomain::checkBounds(void)",
1703              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1704              << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1705
1706      if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1707        ERROR("CDomain::checkBounds(void)",
1708              << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1709
1710      // In case of reading UGRID bounds values are not required
1711      hasBounds = (!bounds_lat_1d.isEmpty() || !bounds_lat_2d.isEmpty() );
1712   }
1713   CATCH_DUMP_ATTR
1714
1715   void CDomain::checkArea(void)
1716   TRY
1717   {
1718     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
1719     hasArea = !area.isEmpty();
1720     if (hasArea && !hasAreaValue)
1721     {
1722       if (area.extent(0) != ni || area.extent(1) != nj)
1723       {
1724         ERROR("CDomain::checkArea(void)",
1725               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1726               << "The area does not have the same size as the local domain." << std::endl
1727               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1728               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1729       }
1730//       if (areavalue.isEmpty())
1731//       {
1732//          areavalue.resize(ni*nj);
1733//         for (int j = 0; j < nj; ++j)
1734//         {
1735//           for (int i = 0; i < ni; ++i)
1736//           {
1737//             int k = j * ni + i;
1738//             areavalue(k) = area(i,j);
1739//           }
1740//         }
1741//       }
1742     }
1743   }
1744   CATCH_DUMP_ATTR
1745
1746   void CDomain::checkLonLat()
1747   TRY
1748   {
1749     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1750                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1751     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1752     if (hasLonLat && !hasLonLatValue)
1753     {
1754       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1755         ERROR("CDomain::checkLonLat()",
1756               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1757               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1758               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1759
1760       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1761       {
1762         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1763           ERROR("CDomain::checkLonLat()",
1764                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1765                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1766                 << "Local size is " << i_index.numElements() << "." << std::endl
1767                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1768       }
1769
1770       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1771       {
1772         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1773           ERROR("CDomain::checkLonLat()",
1774                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1775                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1776                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1777                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1778       }
1779
1780       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1781         ERROR("CDomain::checkLonLat()",
1782               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1783               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1784               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1785
1786       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1787       {
1788         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1789           ERROR("CDomain::checkLonLat()",
1790                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1791                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1792                 << "Local size is " << i_index.numElements() << "." << std::endl
1793                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1794       }
1795
1796       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1797       {
1798         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1799           ERROR("CDomain::checkLonLat()",
1800                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1801                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1802                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1803                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1804       }
1805     }
1806   }
1807   CATCH_DUMP_ATTR
1808
1809   //----------------------------------------------------------------
1810   void CDomain::checkTiles()
1811   TRY
1812   {
1813     if (!ntiles.isEmpty() && ntiles.getValue() >=1) isTiled_ = true;
1814     if (!tile_only.isEmpty() && tile_only.getValue() == true) {
1815       isTiled_ = true;
1816       isTiledOnly_ = true;
1817     }
1818
1819     if (isTiled_)
1820     {
1821
1822       // Attributes tile_ni/nj and tile_i/jbegin are mandatory for tiled domains
1823       if (tile_ni.numElements() != ntiles || tile_ibegin.numElements() != ntiles)
1824           ERROR("CDomain::checkTiles()",
1825                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1826                 << "'tile_ni' or 'tile_ibegin' are ill defined: these attributes must be specified for tiled domains and be of the size 'ntiles'." << std::endl
1827                 << "The number of tiles is  " << ntiles.getValue() << " while the size of 'tile_ni' is " << tile_ni.numElements() 
1828                 << " and the size of 'tile_ibegin' is " << tile_ibegin.numElements() << ".");
1829
1830       if (tile_nj.numElements() != ntiles || tile_jbegin.numElements() != ntiles)
1831           ERROR("CDomain::checkTiles()",
1832                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1833                 << "'tile_nj' or 'tile_jbegin' are ill defined: these attributes must be specified for tiled domains and be of the size 'ntiles'." << std::endl
1834                 << "The number of tiles is  " << ntiles.getValue() << " while the size of 'tile_nj' is " << tile_nj.numElements()
1835                 << " and the size of 'tile_jbegin' is " << tile_jbegin.numElements() << ".");
1836
1837       // Check on consistency of individual tile sizes and local domain size
1838       int sizeTiles = 0;
1839       for (int iTile = 0; iTile < ntiles.getValue(); ++iTile)
1840       {
1841         sizeTiles += tile_ni(iTile) * tile_nj(iTile);
1842       }
1843       if (sizeTiles != (ni*nj))
1844           ERROR("CDomain::checkTiles()",
1845                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1846                 << "tiles should cover the entire local domain and cannot overlap." << std::endl << ".");
1847
1848
1849       // Fill in tile_data_ni/nj and tile_data_i/jbegin if they are not specified
1850       if (tile_data_ni.numElements() == 0 || tile_data_nj.numElements() == 0)
1851       {
1852         tile_data_ni.resize(ntiles);
1853         tile_data_nj.resize(ntiles);
1854         for (int iTile = 0; iTile < ntiles.getValue(); ++iTile)
1855         {
1856           tile_data_ni(iTile) = tile_ni(iTile);
1857           tile_data_nj(iTile) = tile_nj(iTile);
1858         }
1859       }
1860
1861       if (tile_data_ibegin.numElements() == 0 || tile_data_jbegin.numElements() == 0)
1862       {
1863         tile_data_ibegin.resize(ntiles);
1864         tile_data_jbegin.resize(ntiles);
1865         tile_data_ibegin = 0;
1866         tile_data_jbegin = 0;
1867       }
1868
1869     } // isTiled_
1870   }
1871   CATCH_DUMP_ATTR
1872
1873   //----------------------------------------------------------------
1874   int CDomain::getTileDataISize(int tileId) const
1875   TRY
1876   {
1877     int retvalue = (tile_data_ni(tileId) > tile_ni(tileId)) ? tile_data_ni(tileId) : tile_ni(tileId);
1878     return retvalue;
1879   }
1880   CATCH
1881
1882   //----------------------------------------------------------------
1883   int CDomain::getTileDataJSize(int tileId) const
1884   TRY
1885   {
1886     int retvalue = (tile_data_nj(tileId) > tile_nj(tileId)) ? tile_data_nj(tileId) : tile_nj(tileId);
1887     return retvalue;
1888   }
1889   CATCH
1890
1891   void CDomain::checkAttributesOnClientAfterTransformation()
1892   TRY
1893   {
1894     CContext* context=CContext::getCurrent() ;
1895
1896     if (this->isClientAfterTransformationChecked) return;
1897     if (context->hasClient)
1898     {
1899      this->computeConnectedClients();
1900       if (hasLonLat)
1901         if (!context->hasServer)
1902           this->completeLonLatClient();
1903     }
1904
1905     this->isClientAfterTransformationChecked = true;
1906   }
1907   CATCH_DUMP_ATTR
1908
1909   //----------------------------------------------------------------
1910   // Divide function checkAttributes into 2 seperate ones
1911   // This function only checks all attributes of current domain
1912   void CDomain::checkAttributesOnClient()
1913   TRY
1914   {
1915     if (this->isClientChecked) return;
1916     CContext* context=CContext::getCurrent();
1917
1918      if (context->hasClient && !context->hasServer)
1919      {
1920        this->checkDomain();
1921        this->checkBounds();
1922        this->checkArea();
1923        this->checkLonLat();
1924        this->checkTiles();
1925      }
1926
1927      if (context->hasClient && !context->hasServer)
1928      { // Ct client uniquement
1929         this->checkMask();
1930         this->checkDomainData();
1931         this->checkCompression();
1932         this->computeLocalMask() ;
1933      }
1934      else
1935      { // Ct serveur uniquement
1936      }
1937
1938      this->isClientChecked = true;
1939   }
1940   CATCH_DUMP_ATTR
1941
1942   // Send all checked attributes to server
1943   void CDomain::sendCheckedAttributes()
1944   TRY
1945   {
1946     if (!this->isClientChecked) checkAttributesOnClient();
1947     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1948     CContext* context=CContext::getCurrent() ;
1949
1950     if (this->isChecked) return;
1951     if (context->hasClient)
1952     {
1953       sendAttributes();
1954     }
1955     this->isChecked = true;
1956   }
1957   CATCH_DUMP_ATTR
1958
1959   void CDomain::checkAttributes(void)
1960   TRY
1961   {
1962      if (this->isChecked) return;
1963      CContext* context=CContext::getCurrent() ;
1964
1965      this->checkDomain();
1966      this->checkLonLat();
1967      this->checkBounds();
1968      this->checkArea();
1969
1970      if (context->hasClient)
1971      { // Ct client uniquement
1972         this->checkMask();
1973         this->checkDomainData();
1974         this->checkCompression();
1975         this->computeLocalMask() ;
1976
1977      }
1978      else
1979      { // Ct serveur uniquement
1980      }
1981
1982      if (context->hasClient)
1983      {
1984        this->computeConnectedClients();
1985        this->completeLonLatClient();
1986      }
1987
1988      this->isChecked = true;
1989   }
1990   CATCH_DUMP_ATTR
1991
1992   void CDomain::initLonLatValue(void)
1993   TRY
1994   {
1995      CContext* context=CContext::getCurrent() ;
1996
1997      if (context->hasClient)
1998      {
1999        this->completeLonLatClient();
2000      }
2001
2002   }
2003   CATCH_DUMP_ATTR
2004  /*!
2005     Compute the connection of a client to other clients to determine which clients to send attributes to.
2006     The sending clients are supposed to already know the distribution of receiving clients (In simple cases, it's band)
2007     The connection among clients is calculated by using global index.
2008     A client connects to other clients which holds the same global index as it.     
2009  */
2010  void CDomain::computeConnectedClients()
2011  TRY
2012  {
2013    CContext* context=CContext::getCurrent() ;
2014   
2015    // This line should be changed soon.
2016    int nbSrvPools = (context->hasServer) ? (context->hasClient ? context->clientPrimServer.size() : 0) : 1;
2017
2018    nbSenders.clear();
2019    connectedServerRank_.clear();
2020
2021    for (int p = 0; p < nbSrvPools; ++p)
2022    {
2023      CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[p] : context->client;
2024      int nbServer = client->serverSize;
2025      int nbClient = client->clientSize;
2026      int rank     = client->clientRank;
2027      bool doComputeGlobalIndexServer = true;
2028
2029      if (connectedServerRank_.find(nbServer) == connectedServerRank_.end())
2030      {
2031
2032        if (indSrv_.find(nbServer) == indSrv_.end())
2033        {
2034          int i,j,i_ind,j_ind, nbIndex=i_index.numElements();
2035          int globalIndexCount = i_index.numElements();
2036          // Fill in index
2037          CArray<size_t,1> globalIndexDomain(nbIndex);
2038          size_t globalIndex;
2039
2040          for (i = 0; i < nbIndex; ++i)
2041          {
2042            i_ind=i_index(i);
2043            j_ind=j_index(i);
2044            globalIndex = i_ind + j_ind * ni_glo;
2045            globalIndexDomain(i) = globalIndex;
2046          }
2047
2048          if (globalLocalIndexMap_.empty())
2049          {
2050            for (i = 0; i < nbIndex; ++i)
2051              globalLocalIndexMap_[globalIndexDomain(i)] = i;
2052          }
2053
2054          size_t globalSizeIndex = 1, indexBegin, indexEnd;
2055          int range, clientSize = client->clientSize;
2056          std::vector<int> nGlobDomain(2);
2057          nGlobDomain[0] = this->ni_glo;
2058          nGlobDomain[1] = this->nj_glo;
2059          for (int i = 0; i < nGlobDomain.size(); ++i) globalSizeIndex *= nGlobDomain[i];
2060          indexBegin = 0;
2061          if (globalSizeIndex <= clientSize)
2062          {
2063            indexBegin = rank%globalSizeIndex;
2064            indexEnd = indexBegin;
2065          }
2066          else
2067          {
2068            for (int i = 0; i < clientSize; ++i)
2069            {
2070              range = globalSizeIndex / clientSize;
2071              if (i < (globalSizeIndex%clientSize)) ++range;
2072              if (i == client->clientRank) break;
2073              indexBegin += range;
2074            }
2075            indexEnd = indexBegin + range - 1;
2076          }
2077
2078          // Even if servers have no index, they must received something from client
2079          // We only use several client to send "empty" message to these servers
2080          CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2081          std::vector<int> serverZeroIndex;
2082          if (isUnstructed_) serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 0);
2083          else serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 1);
2084
2085          std::list<int> serverZeroIndexLeader;
2086          std::list<int> serverZeroIndexNotLeader;
2087          CContextClient::computeLeader(client->clientRank, client->clientSize, serverZeroIndex.size(), serverZeroIndexLeader, serverZeroIndexNotLeader);
2088          for (std::list<int>::iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
2089            *it = serverZeroIndex[*it];
2090
2091          CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(), client->intraComm);
2092          clientServerMap->computeServerIndexMapping(globalIndexDomain, nbServer);
2093          CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
2094
2095          CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
2096                 ite = globalIndexDomainOnServer.end();
2097          indSrv_[nbServer].swap(globalIndexDomainOnServer);
2098          connectedServerRank_[nbServer].clear();
2099          for (it = indSrv_[nbServer].begin(); it != ite; ++it)
2100            connectedServerRank_[nbServer].push_back(it->first);
2101
2102          for (std::list<int>::const_iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
2103            connectedServerRank_[nbServer].push_back(*it);
2104
2105          // Even if a client has no index, it must connect to at least one server and
2106          // send an "empty" data to this server
2107          if (connectedServerRank_[nbServer].empty())
2108            connectedServerRank_[nbServer].push_back(client->clientRank % client->serverSize);
2109
2110          // Now check if all servers have data to receive. If not, master client will send empty data.
2111          // This ensures that all servers will participate in collective calls upon receiving even if they have no date to receive.
2112          std::vector<int> counts (clientSize);
2113          std::vector<int> displs (clientSize);
2114          displs[0] = 0;
2115          int localCount = connectedServerRank_[nbServer].size() ;
2116          MPI_Gather(&localCount, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, client->intraComm) ;
2117          for (int i = 0; i < clientSize-1; ++i)
2118          {
2119            displs[i+1] = displs[i] + counts[i];
2120          }
2121          std::vector<int> allConnectedServers(displs[clientSize-1]+counts[clientSize-1]);
2122          MPI_Gatherv(&(connectedServerRank_[nbServer])[0], localCount, MPI_INT, &allConnectedServers[0], &counts[0], &displs[0], MPI_INT, 0, client->intraComm);
2123
2124          if ((allConnectedServers.size() != nbServer) && (rank == 0))
2125          {
2126            std::vector<bool> isSrvConnected (nbServer, false);
2127            for (int i = 0; i < allConnectedServers.size(); ++i) isSrvConnected[allConnectedServers[i]] = true;
2128            for (int i = 0; i < nbServer; ++i)
2129            {
2130              if (!isSrvConnected[i]) connectedServerRank_[nbServer].push_back(i);
2131            }
2132          }
2133          nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
2134          delete clientServerMap;
2135        }
2136      }
2137    }
2138  }
2139  CATCH_DUMP_ATTR
2140
2141   /*!
2142     Compute index to write data. We only write data on the zoomed region, therefore, there should
2143     be a map between the complete grid and the reduced grid where we write data.
2144     By using global index we can easily create this kind of mapping.
2145   */
2146   void CDomain::computeWrittenIndex()
2147   TRY
2148   { 
2149      if (computedWrittenIndex_) return;
2150      computedWrittenIndex_ = true;
2151
2152      CContext* context=CContext::getCurrent();     
2153      CContextServer* server = context->server; 
2154
2155      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2156      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2157      nSize[0]        = ni;      nSize[1]  = nj;
2158      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2159      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2160      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
2161      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2162
2163      size_t nbWritten = 0, indGlo;     
2164      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2165                                                          ite = globalLocalIndexMap_.end(), it;         
2166      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2167                                       itSrve = writtenGlobalIndex.end(), itSrv;
2168
2169      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
2170      nbWritten = 0;
2171      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2172      {
2173        indGlo = *itSrv;
2174        if (ite != globalLocalIndexMap_.find(indGlo))
2175        {
2176          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
2177        }
2178        else
2179        {
2180          localIndexToWriteOnServer(nbWritten) = -1;
2181        }
2182        ++nbWritten;
2183      }
2184   }
2185   CATCH_DUMP_ATTR
2186
2187  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
2188  TRY
2189  {
2190    int writtenCommSize;
2191    MPI_Comm_size(writtenComm, &writtenCommSize);
2192    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
2193      return;
2194
2195    if (isCompressible())
2196    {
2197      size_t nbWritten = 0, indGlo;
2198      CContext* context=CContext::getCurrent();     
2199      CContextServer* server = context->server; 
2200
2201      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2202      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
2203      nSize[0]        = ni;      nSize[1]  = nj;
2204      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2205      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2206      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
2207      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2208
2209      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2210                                                          ite = globalLocalIndexMap_.end(), it;   
2211      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2212                                       itSrve = writtenGlobalIndex.end(), itSrv;
2213      std::unordered_map<size_t,size_t> localGlobalIndexMap;
2214      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2215      {
2216        indGlo = *itSrv;
2217        if (ite != globalLocalIndexMap_.find(indGlo))
2218        {
2219          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
2220          ++nbWritten;
2221        }                 
2222      }
2223
2224      nbWritten = 0;
2225      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2226      {
2227        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2228        {
2229          ++nbWritten;
2230        }
2231      }
2232
2233      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2234      nbWritten = 0;
2235      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2236      {
2237        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2238        {
2239          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2240          ++nbWritten;
2241        }
2242      }
2243
2244      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2245      bool distributed_glo, distributed=isDistributed() ;
2246      MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, writtenComm) ;
2247     
2248      if (distributed_glo)
2249      {
2250             
2251        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2252        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2253        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2254      }
2255      else
2256        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2257      }
2258  }
2259  CATCH_DUMP_ATTR
2260
2261  /*!
2262    Send all attributes from client to connected clients
2263    The attributes will be rebuilt on receiving side
2264  */
2265  void CDomain::sendAttributes()
2266  TRY
2267  {
2268    sendDistributionAttributes();
2269    sendIndex();       
2270    sendLonLat();
2271    sendArea();   
2272    sendDataIndex();
2273  }
2274  CATCH
2275  /*!
2276    Send global index from client to connected client(s)
2277  */
2278  void CDomain::sendIndex()
2279  TRY
2280  {
2281    int ns, n, i, j, ind, nv, idx;
2282    std::list<CContextClient*>::iterator it;
2283    for (it=clients.begin(); it!=clients.end(); ++it)
2284    {
2285      CContextClient* client = *it;
2286
2287      int serverSize = client->serverSize;
2288      CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2289
2290      list<CMessage> list_msgsIndex;
2291      list<CArray<int,1> > list_indGlob;
2292
2293      std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2294      iteIndex = indSrv_[serverSize].end();
2295      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2296      {
2297        int nbIndGlob = 0;
2298        int rank = connectedServerRank_[serverSize][k];
2299        itIndex = indSrv_[serverSize].find(rank);
2300        if (iteIndex != itIndex)
2301          nbIndGlob = itIndex->second.size();
2302
2303        list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2304
2305        CArray<int,1>& indGlob = list_indGlob.back();
2306        for (n = 0; n < nbIndGlob; ++n)
2307        {
2308          indGlob(n) = static_cast<int>(itIndex->second[n]);
2309        }
2310
2311        list_msgsIndex.push_back(CMessage());
2312        list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
2313        list_msgsIndex.back() << isCurvilinear;
2314        list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2315       
2316        eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2317      }
2318
2319      client->sendEvent(eventIndex);
2320    }
2321  }
2322  CATCH_DUMP_ATTR
2323
2324  /*!
2325    Send distribution from client to other clients
2326    Because a client in a level knows correctly the grid distribution of client on the next level
2327    it calculates this distribution then sends it to the corresponding clients on the next level
2328  */
2329  void CDomain::sendDistributionAttributes(void)
2330  TRY
2331  {
2332    std::list<CContextClient*>::iterator it;
2333    for (it=clients.begin(); it!=clients.end(); ++it)
2334    {
2335      CContextClient* client = *it;
2336      int nbServer = client->serverSize;
2337      std::vector<int> nGlobDomain(2);
2338      nGlobDomain[0] = this->ni_glo;
2339      nGlobDomain[1] = this->nj_glo;
2340
2341      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2342      if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2343      else serverDescription.computeServerDistribution(false, 1);
2344
2345      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2346      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2347
2348      CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2349      if (client->isServerLeader())
2350      {
2351        std::list<CMessage> msgs;
2352
2353        const std::list<int>& ranks = client->getRanksServerLeader();
2354        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2355        {
2356          // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2357          const int ibegin_srv = serverIndexBegin[*itRank][0];
2358          const int jbegin_srv = serverIndexBegin[*itRank][1];
2359          const int ni_srv = serverDimensionSizes[*itRank][0];
2360          const int nj_srv = serverDimensionSizes[*itRank][1];
2361
2362          msgs.push_back(CMessage());
2363          CMessage& msg = msgs.back();
2364          msg << this->getId() ;
2365          msg << isUnstructed_;
2366          msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2367          msg << ni_glo.getValue() << nj_glo.getValue();
2368          msg << isCompressible_;
2369
2370          event.push(*itRank,1,msg);
2371        }
2372        client->sendEvent(event);
2373      }
2374      else client->sendEvent(event);
2375    }
2376  }
2377  CATCH_DUMP_ATTR
2378
2379  /*!
2380    Send area from client to connected client(s)
2381  */
2382  void CDomain::sendArea()
2383  TRY
2384  {
2385    if (!hasArea) return;
2386
2387    int ns, n, i, j, ind, nv, idx;
2388    std::list<CContextClient*>::iterator it;
2389
2390    for (it=clients.begin(); it!=clients.end(); ++it)
2391    {
2392      CContextClient* client = *it;
2393      int serverSize = client->serverSize;
2394
2395      // send area for each connected server
2396      CEventClient eventArea(getType(), EVENT_ID_AREA);
2397
2398      list<CMessage> list_msgsArea;
2399      list<CArray<double,1> > list_area;
2400
2401      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2402      iteMap = indSrv_[serverSize].end();
2403      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2404      {
2405        int nbData = 0;
2406        int rank = connectedServerRank_[serverSize][k];
2407        it = indSrv_[serverSize].find(rank);
2408        if (iteMap != it)
2409          nbData = it->second.size();
2410        list_area.push_back(CArray<double,1>(nbData));
2411
2412        const std::vector<size_t>& temp = it->second;
2413        for (n = 0; n < nbData; ++n)
2414        {
2415          idx = static_cast<int>(it->second[n]);
2416          list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2417        }
2418
2419        list_msgsArea.push_back(CMessage());
2420        list_msgsArea.back() << this->getId() << hasArea;
2421        list_msgsArea.back() << list_area.back();
2422        eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2423      }
2424      client->sendEvent(eventArea);
2425    }
2426  }
2427  CATCH_DUMP_ATTR
2428
2429  /*!
2430    Send longitude and latitude from client to servers
2431    Each client send long and lat information to corresponding connected clients(s).
2432    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2433  */
2434  void CDomain::sendLonLat()
2435  TRY
2436  {
2437    if (!hasLonLat) return;
2438
2439    int ns, n, i, j, ind, nv, idx;
2440    std::list<CContextClient*>::iterator it;
2441    for (it=clients.begin(); it!=clients.end(); ++it)
2442    {
2443      CContextClient* client = *it;
2444      int serverSize = client->serverSize;
2445
2446      // send lon lat for each connected server
2447      CEventClient eventLon(getType(), EVENT_ID_LON);
2448      CEventClient eventLat(getType(), EVENT_ID_LAT);
2449
2450      list<CMessage> list_msgsLon, list_msgsLat;
2451      list<CArray<double,1> > list_lon, list_lat;
2452      list<CArray<double,2> > list_boundslon, list_boundslat;
2453
2454      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2455      iteMap = indSrv_[serverSize].end();
2456      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2457      {
2458        int nbData = 0;
2459        int rank = connectedServerRank_[serverSize][k];
2460        it = indSrv_[serverSize].find(rank);
2461        if (iteMap != it)
2462          nbData = it->second.size();
2463
2464        list_lon.push_back(CArray<double,1>(nbData));
2465        list_lat.push_back(CArray<double,1>(nbData));
2466
2467        if (hasBounds)
2468        {
2469          int nvertexValue ;
2470          if (type==type_attr::rectilinear || type==type_attr::curvilinear) nvertexValue=4 ;
2471          else nvertexValue=nvertex ;
2472          list_boundslon.push_back(CArray<double,2>(nvertexValue, nbData));
2473          list_boundslat.push_back(CArray<double,2>(nvertexValue, nbData));
2474        }
2475
2476        CArray<double,1>& lon = list_lon.back();
2477        CArray<double,1>& lat = list_lat.back();
2478        const std::vector<size_t>& temp = it->second;
2479        for (n = 0; n < nbData; ++n)
2480        {
2481          idx = static_cast<int>(it->second[n]);
2482          int localInd = globalLocalIndexMap_[idx];
2483          lon(n) = lonvalue(localInd);
2484          lat(n) = latvalue(localInd);
2485
2486          if (hasBounds)
2487          {
2488            CArray<double,2>& boundslon = list_boundslon.back();
2489            CArray<double,2>& boundslat = list_boundslat.back();
2490           
2491            int nvertexValue ;
2492            if (type==type_attr::rectilinear || type==type_attr::curvilinear) nvertexValue=4 ;
2493            else nvertexValue=nvertex ;
2494           
2495            for (nv = 0; nv < nvertexValue; ++nv)
2496            {
2497              boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2498              boundslat(nv, n) = bounds_latvalue(nv, localInd);
2499            }
2500          }
2501        }
2502
2503        list_msgsLon.push_back(CMessage());
2504        list_msgsLat.push_back(CMessage());
2505
2506        list_msgsLon.back() << this->getId() << hasLonLat;
2507        if (hasLonLat) 
2508          list_msgsLon.back() << list_lon.back();
2509        list_msgsLon.back()  << hasBounds;
2510        if (hasBounds)
2511        {
2512          list_msgsLon.back() << list_boundslon.back();
2513        }
2514
2515        list_msgsLat.back() << this->getId() << hasLonLat;
2516        if (hasLonLat)
2517          list_msgsLat.back() << list_lat.back();
2518        list_msgsLat.back() << hasBounds;
2519        if (hasBounds)
2520        {         
2521          list_msgsLat.back() << list_boundslat.back();
2522        }
2523
2524        eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2525        eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2526      }
2527      client->sendEvent(eventLon);
2528      client->sendEvent(eventLat);
2529    }
2530  }
2531  CATCH_DUMP_ATTR
2532
2533  /*!
2534    Send data index to corresponding connected clients.
2535    Data index can be compressed however, we always send decompressed data index
2536    and they will be compressed on receiving.
2537    The compressed index are represented with 1 and others are represented with -1
2538  */
2539  void CDomain::sendDataIndex()
2540  TRY
2541  {
2542    int ns, n, i, j, ind, nv, idx;
2543    std::list<CContextClient*>::iterator it;
2544    for (it=clients.begin(); it!=clients.end(); ++it)
2545    {
2546      CContextClient* client = *it;
2547
2548      int serverSize = client->serverSize;
2549
2550      // send area for each connected server
2551      CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2552
2553      list<CMessage> list_msgsDataIndex;
2554      list<CArray<int,1> > list_data_i_index, list_data_j_index;
2555
2556      int nbIndex = i_index.numElements();
2557      int niByIndex = max(i_index) - min(i_index) + 1;
2558      int njByIndex = max(j_index) - min(j_index) + 1; 
2559      int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2560      int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2561
2562     
2563      CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2564      dataIIndex = -1; 
2565      dataJIndex = -1;
2566      ind = 0;
2567
2568      for (idx = 0; idx < data_i_index.numElements(); ++idx)
2569      {
2570        int dataIidx = data_i_index(idx) + data_ibegin;
2571        int dataJidx = data_j_index(idx) + data_jbegin;
2572        if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2573            (0 <= dataJidx) && (dataJidx < dataJindexBound))
2574        {
2575          dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2576          dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2577        }
2578      }
2579
2580      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2581      iteMap = indSrv_[serverSize].end();
2582      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2583      {
2584        int nbData = 0;
2585        int rank = connectedServerRank_[serverSize][k];
2586        it = indSrv_[serverSize].find(rank);
2587        if (iteMap != it)
2588          nbData = it->second.size();
2589        list_data_i_index.push_back(CArray<int,1>(nbData));
2590        list_data_j_index.push_back(CArray<int,1>(nbData));
2591
2592        const std::vector<size_t>& temp = it->second;
2593        for (n = 0; n < nbData; ++n)
2594        {
2595          idx = static_cast<int>(it->second[n]);
2596          i = globalLocalIndexMap_[idx];
2597          list_data_i_index.back()(n) = dataIIndex(i);
2598          list_data_j_index.back()(n) = dataJIndex(i);
2599        }
2600
2601        list_msgsDataIndex.push_back(CMessage());
2602        list_msgsDataIndex.back() << this->getId();
2603        list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2604        eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2605      }
2606      client->sendEvent(eventDataIndex);
2607    }
2608  }
2609  CATCH
2610 
2611  bool CDomain::dispatchEvent(CEventServer& event)
2612  TRY
2613  {
2614    if (SuperClass::dispatchEvent(event)) return true;
2615    else
2616    {
2617      switch(event.type)
2618      {
2619        case EVENT_ID_SERVER_ATTRIBUT:
2620          recvDistributionAttributes(event);
2621          return true;
2622          break;
2623        case EVENT_ID_INDEX:
2624          recvIndex(event);
2625          return true;
2626          break;
2627        case EVENT_ID_LON:
2628          recvLon(event);
2629          return true;
2630          break;
2631        case EVENT_ID_LAT:
2632          recvLat(event);
2633          return true;
2634          break;
2635        case EVENT_ID_AREA:
2636          recvArea(event);
2637          return true;
2638          break; 
2639        case EVENT_ID_DATA_INDEX:
2640          recvDataIndex(event);
2641          return true;
2642          break;
2643        default:
2644          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2645                << "Unknown Event");
2646          return false;
2647       }
2648    }
2649  }
2650  CATCH
2651
2652  /*!
2653    Receive index event from clients(s)
2654    \param[in] event event contain info about rank and associated index
2655  */
2656  void CDomain::recvIndex(CEventServer& event)
2657  TRY
2658  {
2659    string domainId;
2660    std::map<int, CBufferIn*> rankBuffers;
2661
2662    list<CEventServer::SSubEvent>::iterator it;
2663    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2664    {     
2665      CBufferIn* buffer = it->buffer;
2666      *buffer >> domainId;
2667      rankBuffers[it->rank] = buffer;       
2668    }
2669    get(domainId)->recvIndex(rankBuffers);
2670  }
2671  CATCH
2672
2673  /*!
2674    Receive index information from client(s). We use the global index for mapping index between
2675    sending clients and receiving clients.
2676    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2677  */
2678  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2679  TRY
2680  {
2681    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2682    recvClientRanks_.resize(nbReceived);       
2683
2684    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2685    ind = 0;
2686    for (ind = 0; it != ite; ++it, ++ind)
2687    {       
2688       recvClientRanks_[ind] = it->first;
2689       CBufferIn& buffer = *(it->second);
2690       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2691       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2692    }
2693    int nbIndGlob = 0;
2694    for (i = 0; i < nbReceived; ++i)
2695    {
2696      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2697    }
2698   
2699    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2700    i_index.resize(nbIndGlob);
2701    j_index.resize(nbIndGlob);   
2702    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2703
2704    nbIndGlob = 0;
2705    for (i = 0; i < nbReceived; ++i)
2706    {
2707      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2708      for (ind = 0; ind < tmp.numElements(); ++ind)
2709      {
2710         index = tmp(ind);
2711         if (0 == globalLocalIndexMap_.count(index))
2712         {
2713           iIndex = (index%ni_glo)-ibegin;
2714           iIndex = (iIndex < 0) ? 0 : iIndex;
2715           jIndex = (index/ni_glo)-jbegin;
2716           jIndex = (jIndex < 0) ? 0 : jIndex;
2717           nbIndLoc = iIndex + ni * jIndex;
2718           i_index(nbIndGlob) = index % ni_glo;
2719           j_index(nbIndGlob) = index / ni_glo;
2720           globalLocalIndexMap_[index] = nbIndGlob;
2721           ++nbIndGlob;
2722         } 
2723      } 
2724    } 
2725
2726    if (nbIndGlob==0)
2727    {
2728      i_index.resize(nbIndGlob);
2729      j_index.resize(nbIndGlob);
2730    }
2731    else
2732    {
2733      i_index.resizeAndPreserve(nbIndGlob);
2734      j_index.resizeAndPreserve(nbIndGlob);
2735    }
2736
2737    domainMask.resize(0); // Mask is not defined anymore on servers
2738  }
2739  CATCH
2740
2741  /*!
2742    Receive attributes event from clients(s)
2743    \param[in] event event contain info about rank and associated attributes
2744  */
2745  void CDomain::recvDistributionAttributes(CEventServer& event)
2746  TRY
2747  {
2748    CBufferIn* buffer=event.subEvents.begin()->buffer;
2749    string domainId ;
2750    *buffer>>domainId ;
2751    get(domainId)->recvDistributionAttributes(*buffer);
2752  }
2753  CATCH
2754
2755  /*!
2756    Receive attributes from client(s)
2757    \param[in] rank rank of client source
2758    \param[in] buffer message containing attributes info
2759  */
2760  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2761  TRY
2762  {
2763    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2764    int ni_glo_tmp, nj_glo_tmp;
2765    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2766           >> ni_glo_tmp >> nj_glo_tmp
2767           >> isCompressible_;
2768
2769    ni.setValue(ni_tmp);
2770    ibegin.setValue(ibegin_tmp);
2771    nj.setValue(nj_tmp);
2772    jbegin.setValue(jbegin_tmp);
2773    ni_glo.setValue(ni_glo_tmp);
2774    nj_glo.setValue(nj_glo_tmp);
2775
2776  }
2777 CATCH_DUMP_ATTR
2778  /*!
2779    Receive longitude event from clients(s)
2780    \param[in] event event contain info about rank and associated longitude
2781  */
2782  void CDomain::recvLon(CEventServer& event)
2783  TRY
2784  {
2785    string domainId;
2786    std::map<int, CBufferIn*> rankBuffers;
2787
2788    list<CEventServer::SSubEvent>::iterator it;
2789    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2790    {     
2791      CBufferIn* buffer = it->buffer;
2792      *buffer >> domainId;
2793      rankBuffers[it->rank] = buffer;       
2794    }
2795    get(domainId)->recvLon(rankBuffers);
2796  }
2797  CATCH
2798
2799  /*!
2800    Receive longitude information from client(s)
2801    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2802  */
2803  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2804  TRY
2805  {
2806    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2807    if (nbReceived != recvClientRanks_.size())
2808      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2809           << "The number of sending clients is not correct."
2810           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2811
2812    vector<CArray<double,1> > recvLonValue(nbReceived);
2813    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2814    for (i = 0; i < recvClientRanks_.size(); ++i)
2815    {
2816      int rank = recvClientRanks_[i];
2817      CBufferIn& buffer = *(rankBuffers[rank]);
2818      buffer >> hasLonLat;
2819      if (hasLonLat)
2820        buffer >> recvLonValue[i];
2821      buffer >> hasBounds;
2822      if (hasBounds)
2823        buffer >> recvBoundsLonValue[i];
2824    }
2825
2826    if (hasLonLat)
2827    {
2828      int nbLonInd = 0;
2829      for (i = 0; i < nbReceived; ++i)
2830      {
2831        nbLonInd += recvLonValue[i].numElements();
2832      }
2833   
2834      if (nbLonInd != globalLocalIndexMap_.size())
2835        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2836                 << "something must be wrong with longitude index "<< std::endl;
2837
2838      nbLonInd = globalLocalIndexMap_.size();
2839      lonvalue.resize(nbLonInd);
2840      if (hasBounds)
2841      {
2842        int nvertexValue = (type==type_attr::rectilinear || type==type_attr::curvilinear) ? 4 : nvertex ;
2843        bounds_lonvalue.resize(nvertexValue,nbLonInd);
2844        bounds_lonvalue = 0.;
2845      }
2846
2847      nbLonInd = 0;
2848      for (i = 0; i < nbReceived; ++i)
2849      {
2850        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2851        CArray<double,1>& tmp = recvLonValue[i];
2852        for (ind = 0; ind < tmp.numElements(); ++ind)
2853        {
2854          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2855          lonvalue(lInd) = tmp(ind); 
2856           if (hasBounds)
2857           {         
2858              int nvertexValue = (type==type_attr::rectilinear || type==type_attr::curvilinear) ? 4 : nvertex ;
2859              for (int nv = 0; nv < nvertexValue; ++nv)
2860                bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2861           }                 
2862        }
2863      }       
2864    }
2865  }
2866  CATCH_DUMP_ATTR
2867
2868  /*!
2869    Receive latitude event from clients(s)
2870    \param[in] event event contain info about rank and associated latitude
2871  */
2872  void CDomain::recvLat(CEventServer& event)
2873  TRY
2874  {
2875    string domainId;
2876    std::map<int, CBufferIn*> rankBuffers;
2877
2878    list<CEventServer::SSubEvent>::iterator it;
2879    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2880    {     
2881      CBufferIn* buffer = it->buffer;
2882      *buffer >> domainId;
2883      rankBuffers[it->rank] = buffer;   
2884    }
2885    get(domainId)->recvLat(rankBuffers);
2886  }
2887  CATCH
2888
2889  /*!
2890    Receive latitude information from client(s)
2891    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2892  */
2893  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2894  TRY
2895  {
2896    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2897    if (nbReceived != recvClientRanks_.size())
2898      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2899           << "The number of sending clients is not correct."
2900           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2901
2902    vector<CArray<double,1> > recvLatValue(nbReceived);
2903    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2904    for (i = 0; i < recvClientRanks_.size(); ++i)
2905    {
2906      int rank = recvClientRanks_[i];
2907      CBufferIn& buffer = *(rankBuffers[rank]);
2908      buffer >> hasLonLat;
2909      if (hasLonLat)
2910        buffer >> recvLatValue[i];
2911      buffer >> hasBounds;
2912      if (hasBounds)
2913        buffer >> recvBoundsLatValue[i];
2914    }
2915
2916    if (hasLonLat)
2917    {
2918      int nbLatInd = 0;
2919      for (i = 0; i < nbReceived; ++i)
2920      {
2921        nbLatInd += recvLatValue[i].numElements();
2922      }
2923   
2924      if (nbLatInd != globalLocalIndexMap_.size())
2925        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2926                << "something must be wrong with latitude index "<< std::endl;
2927
2928      nbLatInd = globalLocalIndexMap_.size();
2929      latvalue.resize(nbLatInd);
2930      if (hasBounds)
2931      {
2932        int nvertexValue = (type==type_attr::rectilinear || type==type_attr::curvilinear) ? 4 : nvertex ;
2933       
2934        bounds_latvalue.resize(nvertexValue,nbLatInd);
2935        bounds_latvalue = 0. ;
2936      }
2937
2938      nbLatInd = 0;
2939      for (i = 0; i < nbReceived; ++i)
2940      {
2941        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2942        CArray<double,1>& tmp = recvLatValue[i];
2943        for (ind = 0; ind < tmp.numElements(); ++ind)
2944        {
2945          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2946          latvalue(lInd) = tmp(ind);   
2947           if (hasBounds)
2948           {
2949             int nvertexValue = (type==type_attr::rectilinear || type==type_attr::curvilinear) ? 4 : nvertex ;
2950
2951             CArray<double,2>& boundslat = recvBoundsLatValue[i];
2952             for (int nv = 0; nv < nvertexValue; ++nv)
2953             bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2954           }   
2955          ++nbLatInd;
2956        }
2957      }       
2958    }
2959  }
2960  CATCH_DUMP_ATTR
2961
2962  /*!
2963    Receive area event from clients(s)
2964    \param[in] event event contain info about rank and associated area
2965  */
2966  void CDomain::recvArea(CEventServer& event)
2967  TRY
2968  {
2969    string domainId;
2970    std::map<int, CBufferIn*> rankBuffers;
2971
2972    list<CEventServer::SSubEvent>::iterator it;
2973    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2974    {     
2975      CBufferIn* buffer = it->buffer;
2976      *buffer >> domainId;
2977      rankBuffers[it->rank] = buffer;     
2978    }
2979    get(domainId)->recvArea(rankBuffers);
2980  }
2981  CATCH
2982
2983  /*!
2984    Receive area information from client(s)
2985    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2986  */
2987  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
2988  TRY
2989  {
2990    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2991    if (nbReceived != recvClientRanks_.size())
2992      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2993           << "The number of sending clients is not correct."
2994           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2995
2996    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2997    for (i = 0; i < recvClientRanks_.size(); ++i)
2998    {
2999      int rank = recvClientRanks_[i];
3000      CBufferIn& buffer = *(rankBuffers[rank]);     
3001      buffer >> hasArea;
3002      if (hasArea)
3003        buffer >> recvAreaValue[i];
3004    }
3005
3006    if (hasArea)
3007    {
3008      int nbAreaInd = 0;
3009      for (i = 0; i < nbReceived; ++i)
3010      {     
3011        nbAreaInd += recvAreaValue[i].numElements();
3012      }
3013
3014      if (nbAreaInd != globalLocalIndexMap_.size())
3015        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3016                 << "something must be wrong with area index "<< std::endl;
3017
3018      nbAreaInd = globalLocalIndexMap_.size();
3019      areavalue.resize(nbAreaInd);
3020      nbAreaInd = 0;     
3021      for (i = 0; i < nbReceived; ++i)
3022      {
3023        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3024        CArray<double,1>& tmp = recvAreaValue[i];
3025        for (ind = 0; ind < tmp.numElements(); ++ind)
3026        {
3027          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3028          areavalue(lInd) = tmp(ind);         
3029        }
3030      }
3031     
3032    }
3033  }
3034  CATCH_DUMP_ATTR
3035
3036  /*!
3037    Compare two domain objects.
3038    They are equal if only if they have identical attributes as well as their values.
3039    Moreover, they must have the same transformations.
3040  \param [in] domain Compared domain
3041  \return result of the comparison
3042  */
3043  bool CDomain::isEqual(CDomain* obj)
3044  TRY
3045  {
3046    vector<StdString> excludedAttr;
3047    excludedAttr.push_back("domain_ref");
3048    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
3049    if (!objEqual) return objEqual;
3050
3051    TransMapTypes thisTrans = this->getAllTransformations();
3052    TransMapTypes objTrans  = obj->getAllTransformations();
3053
3054    TransMapTypes::const_iterator it, itb, ite;
3055    std::vector<ETranformationType> thisTransType, objTransType;
3056    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
3057      thisTransType.push_back(it->first);
3058    for (it = objTrans.begin(); it != objTrans.end(); ++it)
3059      objTransType.push_back(it->first);
3060
3061    if (thisTransType.size() != objTransType.size()) return false;
3062    for (int idx = 0; idx < thisTransType.size(); ++idx)
3063      objEqual &= (thisTransType[idx] == objTransType[idx]);
3064
3065    return objEqual;
3066  }
3067  CATCH_DUMP_ATTR
3068
3069  /*!
3070    Receive data index event from clients(s)
3071    \param[in] event event contain info about rank and associated index
3072  */
3073  void CDomain::recvDataIndex(CEventServer& event)
3074  TRY
3075  {
3076    string domainId;
3077    std::map<int, CBufferIn*> rankBuffers;
3078
3079    list<CEventServer::SSubEvent>::iterator it;
3080    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3081    {     
3082      CBufferIn* buffer = it->buffer;
3083      *buffer >> domainId;
3084      rankBuffers[it->rank] = buffer;       
3085    }
3086    get(domainId)->recvDataIndex(rankBuffers);
3087  }
3088  CATCH
3089
3090  /*!
3091    Receive data index information from client(s)
3092    A client receives data index from different clients to rebuild its own data index.
3093    Because we use global index + mask info to calculate the sending data to client(s),
3094    this data index must be updated with mask info (maybe it will change in the future)
3095    Because the data index is local, to rebuild data index of received client, we should use global index along with.
3096
3097    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3098  */
3099  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
3100  TRY
3101  {
3102    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
3103    if (nbReceived != recvClientRanks_.size())
3104      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3105           << "The number of sending clients is not correct."
3106           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3107
3108    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3109    for (i = 0; i < recvClientRanks_.size(); ++i)
3110    {
3111      int rank = recvClientRanks_[i];
3112      CBufferIn& buffer = *(rankBuffers[rank]);
3113      buffer >> recvDataIIndex[i];
3114      buffer >> recvDataJIndex[i];
3115    }
3116   
3117    int nbIndex = i_index.numElements();
3118    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3119    dataIIndex = -1; dataJIndex = -1;
3120     
3121    nbIndex = 0;
3122    for (i = 0; i < nbReceived; ++i)
3123    {     
3124      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3125      CArray<int,1>& tmpI = recvDataIIndex[i];   
3126      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3127      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3128          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3129             << "The number of global received index is not coherent with the number of received data index."
3130             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3131
3132      for (ind = 0; ind < tmpI.numElements(); ++ind)
3133      {
3134         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3135         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3136         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3137      } 
3138    }
3139
3140    int nbCompressedData = 0; 
3141    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3142    {
3143       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3144       if ((0 <= indexI) && (0 <= indexJ))
3145         ++nbCompressedData;
3146    }       
3147 
3148    data_i_index.resize(nbCompressedData);
3149    data_j_index.resize(nbCompressedData);
3150
3151    nbCompressedData = 0; 
3152    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3153    {
3154       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3155       if ((0 <= indexI) && (0 <= indexJ))
3156       {
3157          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3158          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3159         ++nbCompressedData;
3160       }
3161    }
3162
3163    // Reset data_ibegin, data_jbegin
3164    data_ibegin.setValue(0);
3165    data_jbegin.setValue(0);
3166  }
3167  CATCH_DUMP_ATTR
3168
3169  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3170  TRY
3171  {
3172    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3173    return transformationMap_.back().second;
3174  }
3175  CATCH_DUMP_ATTR
3176
3177  /*!
3178    Check whether a domain has transformation
3179    \return true if domain has transformation
3180  */
3181  bool CDomain::hasTransformation()
3182  TRY
3183  {
3184    return (!transformationMap_.empty());
3185  }
3186  CATCH_DUMP_ATTR
3187
3188  /*!
3189    Set transformation for current domain. It's the method to move transformation in hierarchy
3190    \param [in] domTrans transformation on domain
3191  */
3192  void CDomain::setTransformations(const TransMapTypes& domTrans)
3193  TRY
3194  {
3195    transformationMap_ = domTrans;
3196  }
3197  CATCH_DUMP_ATTR
3198
3199  /*!
3200    Get all transformation current domain has
3201    \return all transformation
3202  */
3203  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3204  TRY
3205  {
3206    return transformationMap_;
3207  }
3208  CATCH_DUMP_ATTR
3209
3210  void CDomain::duplicateTransformation(CDomain* src)
3211  TRY
3212  {
3213    if (src->hasTransformation())
3214    {
3215      this->setTransformations(src->getAllTransformations());
3216    }
3217  }
3218  CATCH_DUMP_ATTR
3219
3220  /*!
3221   * Go through the hierarchy to find the domain from which the transformations must be inherited
3222   */
3223  void CDomain::solveInheritanceTransformation()
3224  TRY
3225  {
3226    if (hasTransformation() || !hasDirectDomainReference())
3227      return;
3228
3229    CDomain* domain = this;
3230    std::vector<CDomain*> refDomains;
3231    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3232    {
3233      refDomains.push_back(domain);
3234      domain = domain->getDirectDomainReference();
3235    }
3236
3237    if (domain->hasTransformation())
3238      for (size_t i = 0; i < refDomains.size(); ++i)
3239        refDomains[i]->setTransformations(domain->getAllTransformations());
3240  }
3241  CATCH_DUMP_ATTR
3242
3243  void CDomain::setContextClient(CContextClient* contextClient)
3244  TRY
3245  {
3246    if (clientsSet.find(contextClient)==clientsSet.end())
3247    {
3248      clients.push_back(contextClient) ;
3249      clientsSet.insert(contextClient);
3250    }
3251  }
3252  CATCH_DUMP_ATTR
3253
3254  /*!
3255    Parse children nodes of a domain in xml file.
3256    Whenver there is a new transformation, its type and name should be added into this function
3257    \param node child node to process
3258  */
3259  void CDomain::parse(xml::CXMLNode & node)
3260  TRY
3261  {
3262    SuperClass::parse(node);
3263
3264    if (node.goToChildElement())
3265    {
3266      StdString nodeElementName;
3267      do
3268      {
3269        StdString nodeId("");
3270        if (node.getAttributes().end() != node.getAttributes().find("id"))
3271        { nodeId = node.getAttributes()["id"]; }
3272
3273        nodeElementName = node.getElementName();
3274        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3275        it = transformationMapList_.find(nodeElementName);
3276        if (ite != it)
3277        {
3278          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3279                                                                                                                nodeId,
3280                                                                                                                &node)));
3281        }
3282        else
3283        {
3284          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3285                << "The transformation " << nodeElementName << " has not been supported yet.");
3286        }
3287      } while (node.goToNextElement()) ;
3288      node.goToParentElement();
3289    }
3290  }
3291  CATCH_DUMP_ATTR
3292   //----------------------------------------------------------------
3293
3294   DEFINE_REF_FUNC(Domain,domain)
3295
3296   ///---------------------------------------------------------------
3297
3298} // namespace xios
Note: See TracBrowser for help on using the repository browser.