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

Last change on this file since 1337 was 1337, checked in by oabramkina, 6 years ago

Avoiding recalculation of domain/axis distributions for server-pools of the same size. This commit is complimentary to r1263.

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