source: XIOS/dev/XIOS_DEV_CMIP6/src/node/domain.cpp @ 1410

Last change on this file since 1410 was 1410, checked in by oabramkina, 5 years ago

Improving protocol for reading : grid attributes such longitude, latitude, etc are read by clients locally.
It concerns only curvilinear and unstructured domains. Local attributes such as ni/nj, ibegin/jbegin are mandatory as before.

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