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

Last change on this file since 1294 was 1284, checked in by ymipsl, 7 years ago

Fix 1 side effect introduced by previous commit...

YM

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