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

Last change on this file since 1364 was 1364, checked in by oabramkina, 4 years ago

Fixing a bug in case of a varying number of processes per pool and unstructured grid.
Commit has no implications in case of one process per pool.

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