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

Last change on this file since 1427 was 1427, checked in by oabramkina, 3 years ago

Bugfix: correcting erroneous memory deallocations of r1419.

  • 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.4 KB
Line 
1#include "domain.hpp"
2
3#include "attribute_template.hpp"
4#include "object_template.hpp"
5#include "group_template.hpp"
6
7#include "xios_spl.hpp"
8#include "event_client.hpp"
9#include "event_server.hpp"
10#include "buffer_in.hpp"
11#include "message.hpp"
12#include "type.hpp"
13#include "context.hpp"
14#include "context_client.hpp"
15#include "context_server.hpp"
16#include "array_new.hpp"
17#include "distribution_client.hpp"
18#include "server_distribution_description.hpp"
19#include "client_server_mapping_distributed.hpp"
20
21#include <algorithm>
22
23namespace xios {
24
25   /// ////////////////////// Définitions ////////////////////// ///
26
27   CDomain::CDomain(void)
28      : CObjectTemplate<CDomain>(), CDomainAttributes()
29      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_()
30      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
31      , isClientAfterTransformationChecked(false), hasLonLat(false)
32      , isRedistributed_(false), hasPole(false), doZoomByIndex_(false)
33      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
34      , globalLocalIndexMap_(), computedWrittenIndex_(false)
35          , clients()
36   {
37   }
38
39   CDomain::CDomain(const StdString & id)
40      : CObjectTemplate<CDomain>(id), CDomainAttributes()
41      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_() 
42      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
43      , isClientAfterTransformationChecked(false), hasLonLat(false)
44      , isRedistributed_(false), hasPole(false), doZoomByIndex_(false)
45      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
46      , globalLocalIndexMap_(), computedWrittenIndex_(false)
47          , clients()
48   {
49    }
50
51   CDomain::~CDomain(void)
52   {
53   }
54
55   ///---------------------------------------------------------------
56
57   void CDomain::assignMesh(const StdString meshName, const int nvertex)
58   {
59     mesh = CMesh::getMesh(meshName, nvertex);
60   }
61
62   CDomain* CDomain::createDomain()
63   {
64     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
65     return domain;
66   }
67
68   std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
69   bool CDomain::_dummyTransformationMapList = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
70
71   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
72   {
73     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
74     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
75     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
76     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
77     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
78   }
79
80   const std::set<StdString> & CDomain::getRelFiles(void) const
81   {
82      return (this->relFiles);
83   }
84
85
86   /*!
87     Returns the number of indexes written by each server.
88     \return the number of indexes written by each server
89   */
90   int CDomain::getNumberWrittenIndexes(MPI_Comm writtenCom)
91   {
92     int writtenSize;
93     MPI_Comm_size(writtenCom, &writtenSize);
94     return numberWrittenIndexes_[writtenSize];
95   }
96
97   /*!
98     Returns the total number of indexes written by the servers.
99     \return the total number of indexes written by the servers
100   */
101   int CDomain::getTotalNumberWrittenIndexes(MPI_Comm writtenCom)
102   {
103     int writtenSize;
104     MPI_Comm_size(writtenCom, &writtenSize);
105     return totalNumberWrittenIndexes_[writtenSize];
106   }
107
108   /*!
109     Returns the offset of indexes written by each server.
110     \return the offset of indexes written by each server
111   */
112   int CDomain::getOffsetWrittenIndexes(MPI_Comm writtenCom)
113   {
114     int writtenSize;
115     MPI_Comm_size(writtenCom, &writtenSize);
116     return offsetWrittenIndexes_[writtenSize];
117   }
118
119   CArray<int, 1>& CDomain::getCompressedIndexToWriteOnServer(MPI_Comm writtenCom)
120   {
121     int writtenSize;
122     MPI_Comm_size(writtenCom, &writtenSize);
123     return compressedIndexToWriteOnServer[writtenSize];
124   }
125
126   //----------------------------------------------------------------
127
128   /*!
129    * Compute the minimum buffer size required to send the attributes to the server(s).
130    *
131    * \return A map associating the server rank with its minimum buffer size.
132    */
133   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
134   {
135
136     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
137
138     if (client->isServerLeader())
139     {
140       // size estimation for sendDistributionAttribut
141       size_t size = 11 * sizeof(size_t);
142
143       const std::list<int>& ranks = client->getRanksServerLeader();
144       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
145       {
146         if (size > attributesSizes[*itRank])
147           attributesSizes[*itRank] = size;
148       }
149     }
150
151     boost::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->serverSize].end();
152     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
153     for (size_t k = 0; k < connectedServerRank_[client->serverSize].size(); ++k)
154     {
155       int rank = connectedServerRank_[client->serverSize][k];
156       boost::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->serverSize].find(rank);
157       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
158
159       // size estimation for sendIndex (and sendArea which is always smaller or equal)
160       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
161       // if (isCompressible_)
162       // {
163       //   std::map<int, std::vector<int> >::const_iterator itWritten = indWrittenSrv_.find(rank);
164       //   size_t writtenIdxCount = (itWritten != itWrittenIndexEnd) ? itWritten->second.size() : 0;
165       //   sizeIndexEvent += CArray<int,1>::size(writtenIdxCount);
166       // }
167
168       // size estimation for sendLonLat
169       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
170       if (hasBounds)
171         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
172
173       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
174       if (size > attributesSizes[rank])
175         attributesSizes[rank] = size;
176     }
177
178     return attributesSizes;
179   }
180
181   //----------------------------------------------------------------
182
183   bool CDomain::isEmpty(void) const
184   {
185      return ((this->zoom_i_index.isEmpty()) || (0 == this->zoom_i_index.numElements()));
186
187   }
188
189   //----------------------------------------------------------------
190
191   bool CDomain::IsWritten(const StdString & filename) const
192   {
193      return (this->relFiles.find(filename) != this->relFiles.end());
194   }
195
196   bool CDomain::isWrittenCompressed(const StdString& filename) const
197   {
198      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
199   }
200
201   //----------------------------------------------------------------
202
203   bool CDomain::isDistributed(void) const
204   {
205      bool distributed =  !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
206              (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
207      distributed |= (1 == CContext::getCurrent()->client->clientSize);
208
209      return distributed;
210   }
211
212   //----------------------------------------------------------------
213
214   /*!
215    * Test whether the data defined on the domain can be outputted in a compressed way.
216    *
217    * \return true if and only if a mask was defined for this domain
218    */
219   bool CDomain::isCompressible(void) const
220   {
221      return isCompressible_;
222   }
223
224   void CDomain::addRelFile(const StdString & filename)
225   {
226      this->relFiles.insert(filename);
227   }
228
229   void CDomain::addRelFileCompressed(const StdString& filename)
230   {
231      this->relFilesCompressed.insert(filename);
232   }
233
234   StdString CDomain::GetName(void)   { return (StdString("domain")); }
235   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
236   ENodeType CDomain::GetType(void)   { return (eDomain); }
237
238   //----------------------------------------------------------------
239
240   /*!
241      Verify if all distribution information of a domain are available
242      This checking verifies the definition of distribution attributes (ni, nj, ibegin, jbegin)
243   */
244   bool CDomain::distributionAttributesHaveValue() const
245   {
246      bool hasValues = true;
247
248      if (ni.isEmpty() && ibegin.isEmpty() && i_index.isEmpty())
249      {
250        hasValues = false;
251        return hasValues;
252      }
253
254      return hasValues;
255   }
256
257   /*!
258     Redistribute RECTILINEAR or CURVILINEAR domain with a number of local domains.
259   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
260   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
261   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
262    \param [in] nbLocalDomain number of local domain on the domain destination
263   */
264   void CDomain::redistribute(int nbLocalDomain)
265   {
266     if (this->isRedistributed_) return;
267
268     this->isRedistributed_ = true;
269     CContext* context = CContext::getCurrent();
270     // For now the assumption is that secondary server pools consist of the same number of procs.
271     // CHANGE the line below if the assumption changes.
272     CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[0] : context->client;
273     int rankClient = client->clientRank;
274     int rankOnDomain = rankClient%nbLocalDomain;
275
276     if (ni_glo.isEmpty() || ni_glo <= 0 )
277     {
278        ERROR("CDomain::redistribute(int nbLocalDomain)",
279           << "[ Id = " << this->getId() << " ] "
280           << "The global domain is badly defined,"
281           << " check the \'ni_glo\'  value !")
282     }
283
284     if (nj_glo.isEmpty() || nj_glo <= 0 )
285     {
286        ERROR("CDomain::redistribute(int nbLocalDomain)",
287           << "[ Id = " << this->getId() << " ] "
288           << "The global domain is badly defined,"
289           << " check the \'nj_glo\'  value !")
290     }
291
292     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
293     {
294        int globalDomainSize = ni_glo * nj_glo;
295        if (globalDomainSize <= nbLocalDomain)
296        {
297          for (int idx = 0; idx < nbLocalDomain; ++idx)
298          {
299            if (rankOnDomain < globalDomainSize)
300            {
301              int iIdx = rankOnDomain % ni_glo;
302              int jIdx = rankOnDomain / ni_glo;
303              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
304              ni.setValue(1); nj.setValue(1);
305            }
306            else
307            {
308              ibegin.setValue(0); jbegin.setValue(0);
309              ni.setValue(0); nj.setValue(0);
310            }
311          }
312        }
313        else
314        {
315          float njGlo = nj_glo.getValue();
316          float niGlo = ni_glo.getValue();
317          int nbProcOnX, nbProcOnY, range;
318
319          // Compute (approximately) number of segment on x and y axis
320          float yOverXRatio = njGlo/niGlo;
321
322          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
323          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
324
325          // Simple distribution: Sweep from top to bottom, left to right
326          // Calculate local begin on x
327          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
328          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
329          for (int i = 1; i < nbProcOnX; ++i)
330          {
331            range = ni_glo / nbProcOnX;
332            if (i < (ni_glo%nbProcOnX)) ++range;
333            niVec[i-1] = range;
334            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
335          }
336          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
337
338          // Calculate local begin on y
339          for (int j = 1; j < nbProcOnY; ++j)
340          {
341            range = nj_glo / nbProcOnY;
342            if (j < (nj_glo%nbProcOnY)) ++range;
343            njVec[j-1] = range;
344            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
345          }
346          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
347
348          // Now assign value to ni, ibegin, nj, jbegin
349          int iIdx = rankOnDomain % nbProcOnX;
350          int jIdx = rankOnDomain / nbProcOnX;
351
352          if (rankOnDomain != (nbLocalDomain-1))
353          {
354            ibegin.setValue(ibeginVec[iIdx]);
355            jbegin.setValue(jbeginVec[jIdx]);
356            nj.setValue(njVec[jIdx]);
357            ni.setValue(niVec[iIdx]);
358          }
359          else // just merge all the remaining rectangle into the last one
360          {
361            ibegin.setValue(ibeginVec[iIdx]);
362            jbegin.setValue(jbeginVec[jIdx]);
363            nj.setValue(njVec[jIdx]);
364            ni.setValue(ni_glo - ibeginVec[iIdx]);
365          }
366        } 
367     }
368     else  // unstructured domain
369     {
370       if (this->i_index.isEmpty())
371       {
372          int globalDomainSize = ni_glo * nj_glo;
373          if (globalDomainSize <= nbLocalDomain)
374          {
375            for (int idx = 0; idx < nbLocalDomain; ++idx)
376            {
377              if (rankOnDomain < globalDomainSize)
378              {
379                int iIdx = rankOnDomain % ni_glo;
380                int jIdx = rankOnDomain / ni_glo;
381                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
382                ni.setValue(1); nj.setValue(1);
383              }
384              else
385              {
386                ibegin.setValue(0); jbegin.setValue(0);
387                ni.setValue(0); nj.setValue(0);
388              }
389            }
390          }
391          else
392          {
393            float njGlo = nj_glo.getValue();
394            float niGlo = ni_glo.getValue();
395            std::vector<int> ibeginVec(nbLocalDomain,0);
396            std::vector<int> niVec(nbLocalDomain);
397            for (int i = 1; i < nbLocalDomain; ++i)
398            {
399              int range = ni_glo / nbLocalDomain;
400              if (i < (ni_glo%nbLocalDomain)) ++range;
401              niVec[i-1] = range;
402              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
403            }
404            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
405
406            int iIdx = rankOnDomain % nbLocalDomain;
407            ibegin.setValue(ibeginVec[iIdx]);
408            jbegin.setValue(0);
409            ni.setValue(niVec[iIdx]);
410            nj.setValue(1);
411          }
412
413          i_index.resize(ni);         
414          for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
415        }
416        else
417        {
418          ibegin.setValue(this->i_index(0));
419          jbegin.setValue(0);
420          ni.setValue(this->i_index.numElements());
421          nj.setValue(1);
422        }
423     }
424
425     checkDomain();
426   }
427
428   /*!
429     Fill in longitude and latitude whose values are read from file
430   */
431   void CDomain::fillInLonLat()
432   {
433     switch (type)
434     {
435      case type_attr::rectilinear:
436        fillInRectilinearLonLat();
437        break;
438      case type_attr::curvilinear:
439        fillInCurvilinearLonLat();
440        break;
441      case type_attr::unstructured:
442        fillInUnstructuredLonLat();
443        break;
444
445      default:
446      break;
447     }
448     completeLonLatClient() ;
449
450   }
451
452   /*!
453     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
454     Range of longitude value from 0 - 360
455     Range of latitude value from -90 - +90
456   */
457   void CDomain::fillInRectilinearLonLat()
458   {
459     if (!lonvalue_rectilinear_read_from_file.isEmpty())
460     {
461       lonvalue_1d.resize(ni);
462       for (int idx = 0; idx < ni; ++idx)
463         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
464       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
465       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
466     }
467     else
468     {
469       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
470       lonvalue_1d.resize(ni);
471       double lonRange = lon_end - lon_start;
472       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
473
474        // Assign lon value
475       for (int i = 0; i < ni; ++i)
476       {
477         if (0 == (ibegin + i))
478         {
479           lonvalue_1d(i) = lon_start;
480         }
481         else if (ni_glo == (ibegin + i + 1))
482         {
483           lonvalue_1d(i) = lon_end;
484         }
485         else
486         {
487           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
488         }
489       }
490     }
491
492
493     if (!latvalue_rectilinear_read_from_file.isEmpty())
494     {
495       latvalue_1d.resize(nj);
496       for (int idx = 0; idx < nj; ++idx)
497         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
498       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
499       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
500     }
501     else
502     {
503       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
504       latvalue_1d.resize(nj);
505
506       double latRange = lat_end - lat_start;
507       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
508
509       for (int j = 0; j < nj; ++j)
510       {
511         if (0 == (jbegin + j))
512         {
513            latvalue_1d(j) = lat_start;
514         }
515         else if (nj_glo == (jbegin + j + 1))
516         {
517            latvalue_1d(j) = lat_end;
518         }
519         else
520         {
521           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
522         }
523       }
524     }
525   }
526
527    /*
528      Fill in longitude and latitude of curvilinear domain read from a file
529      If there are already longitude and latitude defined by model. We just igonore reading value.
530    */
531   void CDomain::fillInCurvilinearLonLat()
532   {
533     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, jdx);
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, jdx);
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, jdx);
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, jdx);
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(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(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, idx);
616
617        // We dont need these values anymore, so just delete them
618        bounds_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, idx);
628
629        // We dont need these values anymore, so just delete them
630        bounds_latvalue_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      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
1915
1916      nbWritten = 0;
1917      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1918      {
1919        indGlo = *itSrv;
1920        if (ite != globalLocalIndexMap_.find(indGlo))
1921        {
1922          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
1923          ++nbWritten;
1924        }
1925        else
1926        {
1927          localIndexToWriteOnServer(nbWritten) = 0;
1928          ++nbWritten;
1929        }
1930      }
1931     
1932      // if (isCompressible())
1933      // {
1934      //   nbWritten = 0;
1935      //   boost::unordered_map<size_t,size_t> localGlobalIndexMap;
1936      //   for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1937      //   {
1938      //     indGlo = *itSrv;
1939      //     if (ite != globalLocalIndexMap_.find(indGlo))
1940      //     {
1941      //       localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
1942      //       ++nbWritten;
1943      //     }                 
1944      //   }
1945
1946      //   nbWritten = 0;
1947      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1948      //   {
1949      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1950      //     {
1951      //       ++nbWritten;
1952      //     }
1953      //   }
1954
1955      //   compressedIndexToWriteOnServer.resize(nbWritten);
1956      //   nbWritten = 0;
1957      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1958      //   {
1959      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1960      //     {
1961      //       compressedIndexToWriteOnServer(nbWritten) = localGlobalIndexMap[data_i_index(idx)];
1962      //       ++nbWritten;
1963      //     }
1964      //   }
1965
1966      //   numberWrittenIndexes_ = nbWritten;
1967      //   if (isDistributed())
1968      //   {           
1969      //     MPI_Allreduce(&numberWrittenIndexes_, &totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1970      //     MPI_Scan(&numberWrittenIndexes_, &offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1971      //     offsetWrittenIndexes_ -= numberWrittenIndexes_;
1972      //   }
1973      //   else
1974      //     totalNumberWrittenIndexes_ = numberWrittenIndexes_;
1975      // }     
1976   }
1977
1978  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
1979  {
1980    int writtenCommSize;
1981    MPI_Comm_size(writtenComm, &writtenCommSize);
1982    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
1983      return;
1984
1985    if (isCompressible())
1986    {
1987      size_t nbWritten = 0, indGlo;
1988      CContext* context=CContext::getCurrent();     
1989      CContextServer* server = context->server; 
1990
1991      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1992      nBegin[0]       = zoom_ibegin;  nBegin[1] = zoom_jbegin;
1993      nSize[0]        = zoom_ni;      nSize[1]  = zoom_nj;
1994      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1995      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1996      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
1997      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1998
1999      boost::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2000                                                          ite = globalLocalIndexMap_.end(), it;   
2001      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2002                                       itSrve = writtenGlobalIndex.end(), itSrv;
2003      boost::unordered_map<size_t,size_t> localGlobalIndexMap;
2004      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
2005      {
2006        indGlo = *itSrv;
2007        if (ite != globalLocalIndexMap_.find(indGlo))
2008        {
2009          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
2010          ++nbWritten;
2011        }                 
2012      }
2013
2014      nbWritten = 0;
2015      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2016      {
2017        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2018        {
2019          ++nbWritten;
2020        }
2021      }
2022
2023      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2024      nbWritten = 0;
2025      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2026      {
2027        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2028        {
2029          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2030          ++nbWritten;
2031        }
2032      }
2033
2034      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2035      if (isDistributed())
2036      {
2037             
2038        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2039        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2040        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2041      }
2042      else
2043        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2044      }
2045  }
2046
2047  /*!
2048    Send all attributes from client to connected clients
2049    The attributes will be rebuilt on receiving side
2050  */
2051  void CDomain::sendAttributes()
2052  {
2053    sendDistributionAttributes();
2054    sendIndex();       
2055    sendMask();
2056    sendLonLat();
2057    sendArea();   
2058    sendDataIndex();
2059  }
2060
2061  /*!
2062    Send global index and zoom index from client to connected client(s)
2063    zoom index can be smaller than global index
2064  */
2065  void CDomain::sendIndex()
2066  {
2067    int ns, n, i, j, ind, nv, idx;
2068    std::list<CContextClient*>::iterator it;
2069    for (it=clients.begin(); it!=clients.end(); ++it)
2070    {
2071      CContextClient* client = *it;
2072
2073      int serverSize = client->serverSize;
2074      CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2075
2076      list<CMessage> list_msgsIndex;
2077      list<CArray<int,1> > list_indZoom, list_writtenInd, list_indGlob;
2078
2079      boost::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2080      iteIndex = indSrv_[serverSize].end();
2081      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2082      {
2083        int nbIndGlob = 0;
2084        int rank = connectedServerRank_[serverSize][k];
2085        itIndex = indSrv_[serverSize].find(rank);
2086        if (iteIndex != itIndex)
2087          nbIndGlob = itIndex->second.size();
2088
2089        list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2090
2091        CArray<int,1>& indGlob = list_indGlob.back();
2092        for (n = 0; n < nbIndGlob; ++n)
2093        {
2094          indGlob(n) = static_cast<int>(itIndex->second[n]);
2095        }
2096
2097        list_msgsIndex.push_back(CMessage());
2098        list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
2099        list_msgsIndex.back() << isCurvilinear;
2100        list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2101       
2102        eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2103      }
2104
2105      client->sendEvent(eventIndex);
2106    }
2107  }
2108
2109  /*!
2110    Send distribution from client to other clients
2111    Because a client in a level knows correctly the grid distribution of client on the next level
2112    it calculates this distribution then sends it to the corresponding clients on the next level
2113  */
2114  void CDomain::sendDistributionAttributes(void)
2115  {
2116    std::list<CContextClient*>::iterator it;
2117    for (it=clients.begin(); it!=clients.end(); ++it)
2118    {
2119      CContextClient* client = *it;
2120      int nbServer = client->serverSize;
2121      std::vector<int> nGlobDomain(2);
2122      nGlobDomain[0] = this->ni_glo;
2123      nGlobDomain[1] = this->nj_glo;
2124
2125      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2126      if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2127      else serverDescription.computeServerDistribution(false, 1);
2128
2129      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2130      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2131
2132      CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2133      if (client->isServerLeader())
2134      {
2135        std::list<CMessage> msgs;
2136
2137        const std::list<int>& ranks = client->getRanksServerLeader();
2138        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2139        {
2140          // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2141          const int ibegin_srv = serverIndexBegin[*itRank][0];
2142          const int jbegin_srv = serverIndexBegin[*itRank][1];
2143          const int ni_srv = serverDimensionSizes[*itRank][0];
2144          const int nj_srv = serverDimensionSizes[*itRank][1];
2145
2146          msgs.push_back(CMessage());
2147          CMessage& msg = msgs.back();
2148          msg << this->getId() ;
2149          msg << isUnstructed_;
2150          msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2151          msg << global_zoom_ni.getValue() << global_zoom_ibegin.getValue() << global_zoom_nj.getValue() << global_zoom_jbegin.getValue();       
2152          msg << isCompressible_;
2153
2154          event.push(*itRank,1,msg);
2155        }
2156        client->sendEvent(event);
2157      }
2158      else client->sendEvent(event);
2159    }
2160  }
2161
2162  /*!
2163    Send mask index from client to connected(s) clients   
2164  */
2165  void CDomain::sendMask()
2166  {
2167    int ns, n, i, j, ind, nv, idx;
2168    std::list<CContextClient*>::iterator it;
2169    for (it=clients.begin(); it!=clients.end(); ++it)
2170    {
2171      CContextClient* client = *it;
2172      int serverSize = client->serverSize;
2173
2174      // send area for each connected server
2175      CEventClient eventMask(getType(), EVENT_ID_MASK);
2176
2177      list<CMessage> list_msgsMask;
2178      list<CArray<bool,1> > list_mask;
2179
2180      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2181      iteMap = indSrv_[serverSize].end();
2182      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2183      {
2184        int nbData = 0;
2185        int rank = connectedServerRank_[serverSize][k];
2186        it = indSrv_[serverSize].find(rank);
2187        if (iteMap != it)
2188          nbData = it->second.size();
2189        list_mask.push_back(CArray<bool,1>(nbData));
2190
2191        const std::vector<size_t>& temp = it->second;
2192        for (n = 0; n < nbData; ++n)
2193        {
2194          idx = static_cast<int>(it->second[n]);
2195          list_mask.back()(n) = domainMask(globalLocalIndexMap_[idx]);
2196        }
2197
2198        list_msgsMask.push_back(CMessage());
2199        list_msgsMask.back() << this->getId() << list_mask.back();
2200        eventMask.push(rank, nbSenders[serverSize][rank], list_msgsMask.back());
2201      }
2202      client->sendEvent(eventMask);
2203    }
2204  }
2205
2206  /*!
2207    Send area from client to connected client(s)
2208  */
2209  void CDomain::sendArea()
2210  {
2211    if (!hasArea) return;
2212
2213    int ns, n, i, j, ind, nv, idx;
2214    std::list<CContextClient*>::iterator it;
2215
2216    for (it=clients.begin(); it!=clients.end(); ++it)
2217    {
2218      CContextClient* client = *it;
2219      int serverSize = client->serverSize;
2220
2221      // send area for each connected server
2222      CEventClient eventArea(getType(), EVENT_ID_AREA);
2223
2224      list<CMessage> list_msgsArea;
2225      list<CArray<double,1> > list_area;
2226
2227      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2228      iteMap = indSrv_[serverSize].end();
2229      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2230      {
2231        int nbData = 0;
2232        int rank = connectedServerRank_[serverSize][k];
2233        it = indSrv_[serverSize].find(rank);
2234        if (iteMap != it)
2235          nbData = it->second.size();
2236        list_area.push_back(CArray<double,1>(nbData));
2237
2238        const std::vector<size_t>& temp = it->second;
2239        for (n = 0; n < nbData; ++n)
2240        {
2241          idx = static_cast<int>(it->second[n]);
2242          list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2243        }
2244
2245        list_msgsArea.push_back(CMessage());
2246        list_msgsArea.back() << this->getId() << hasArea;
2247        list_msgsArea.back() << list_area.back();
2248        eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2249      }
2250      client->sendEvent(eventArea);
2251    }
2252  }
2253
2254  /*!
2255    Send longitude and latitude from client to servers
2256    Each client send long and lat information to corresponding connected clients(s).
2257    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2258  */
2259  void CDomain::sendLonLat()
2260  {
2261    if (!hasLonLat) return;
2262
2263    int ns, n, i, j, ind, nv, idx;
2264    std::list<CContextClient*>::iterator it;
2265    for (it=clients.begin(); it!=clients.end(); ++it)
2266    {
2267      CContextClient* client = *it;
2268      int serverSize = client->serverSize;
2269
2270      // send lon lat for each connected server
2271      CEventClient eventLon(getType(), EVENT_ID_LON);
2272      CEventClient eventLat(getType(), EVENT_ID_LAT);
2273
2274      list<CMessage> list_msgsLon, list_msgsLat;
2275      list<CArray<double,1> > list_lon, list_lat;
2276      list<CArray<double,2> > list_boundslon, list_boundslat;
2277
2278      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2279      iteMap = indSrv_[serverSize].end();
2280      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2281      {
2282        int nbData = 0;
2283        int rank = connectedServerRank_[serverSize][k];
2284        it = indSrv_[serverSize].find(rank);
2285        if (iteMap != it)
2286          nbData = it->second.size();
2287
2288        list_lon.push_back(CArray<double,1>(nbData));
2289        list_lat.push_back(CArray<double,1>(nbData));
2290
2291        if (hasBounds)
2292        {
2293          list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2294          list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2295        }
2296
2297        CArray<double,1>& lon = list_lon.back();
2298        CArray<double,1>& lat = list_lat.back();
2299        const std::vector<size_t>& temp = it->second;
2300        for (n = 0; n < nbData; ++n)
2301        {
2302          idx = static_cast<int>(it->second[n]);
2303          int localInd = globalLocalIndexMap_[idx];
2304          lon(n) = lonvalue(localInd);
2305          lat(n) = latvalue(localInd);
2306
2307          if (hasBounds)
2308          {
2309            CArray<double,2>& boundslon = list_boundslon.back();
2310            CArray<double,2>& boundslat = list_boundslat.back();
2311
2312            for (nv = 0; nv < nvertex; ++nv)
2313            {
2314              boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2315              boundslat(nv, n) = bounds_latvalue(nv, localInd);
2316            }
2317          }
2318        }
2319
2320        list_msgsLon.push_back(CMessage());
2321        list_msgsLat.push_back(CMessage());
2322
2323        list_msgsLon.back() << this->getId() << hasLonLat;
2324        if (hasLonLat) 
2325          list_msgsLon.back() << list_lon.back();
2326        list_msgsLon.back()  << hasBounds;
2327        if (hasBounds)
2328        {
2329          list_msgsLon.back() << list_boundslon.back();
2330        }
2331
2332        list_msgsLat.back() << this->getId() << hasLonLat;
2333        if (hasLonLat)
2334          list_msgsLat.back() << list_lat.back();
2335        list_msgsLat.back() << hasBounds;
2336        if (hasBounds)
2337        {         
2338          list_msgsLat.back() << list_boundslat.back();
2339        }
2340
2341        eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2342        eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2343      }
2344      client->sendEvent(eventLon);
2345      client->sendEvent(eventLat);
2346    }
2347  }
2348
2349  /*!
2350    Send data index to corresponding connected clients.
2351    Data index can be compressed however, we always send decompressed data index
2352    and they will be compressed on receiving.
2353    The compressed index are represented with 1 and others are represented with -1
2354  */
2355  void CDomain::sendDataIndex()
2356  {
2357    int ns, n, i, j, ind, nv, idx;
2358    std::list<CContextClient*>::iterator it;
2359    for (it=clients.begin(); it!=clients.end(); ++it)
2360    {
2361      CContextClient* client = *it;
2362
2363      int serverSize = client->serverSize;
2364
2365      // send area for each connected server
2366      CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2367
2368      list<CMessage> list_msgsDataIndex;
2369      list<CArray<int,1> > list_data_i_index, list_data_j_index;
2370
2371      int nbIndex = i_index.numElements();
2372      int niByIndex = max(i_index) - min(i_index) + 1;
2373      int njByIndex = max(j_index) - min(j_index) + 1; 
2374      int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2375      int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2376
2377     
2378      CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2379      dataIIndex = -1; 
2380      dataJIndex = -1;
2381      ind = 0;
2382
2383      for (idx = 0; idx < data_i_index.numElements(); ++idx)
2384      {
2385        int dataIidx = data_i_index(idx) + data_ibegin;
2386        int dataJidx = data_j_index(idx) + data_jbegin;
2387        if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2388            (0 <= dataJidx) && (dataJidx < dataJindexBound))
2389        {
2390          dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2391          dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2392        }
2393      }
2394
2395      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2396      iteMap = indSrv_[serverSize].end();
2397      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2398      {
2399        int nbData = 0;
2400        int rank = connectedServerRank_[serverSize][k];
2401        it = indSrv_[serverSize].find(rank);
2402        if (iteMap != it)
2403          nbData = it->second.size();
2404        list_data_i_index.push_back(CArray<int,1>(nbData));
2405        list_data_j_index.push_back(CArray<int,1>(nbData));
2406
2407        const std::vector<size_t>& temp = it->second;
2408        for (n = 0; n < nbData; ++n)
2409        {
2410          idx = static_cast<int>(it->second[n]);
2411          i = globalLocalIndexMap_[idx];
2412          list_data_i_index.back()(n) = dataIIndex(i);
2413          list_data_j_index.back()(n) = dataJIndex(i);
2414        }
2415
2416        list_msgsDataIndex.push_back(CMessage());
2417        list_msgsDataIndex.back() << this->getId();
2418        list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2419        eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2420      }
2421      client->sendEvent(eventDataIndex);
2422    }
2423  }
2424 
2425  bool CDomain::dispatchEvent(CEventServer& event)
2426  {
2427    if (SuperClass::dispatchEvent(event)) return true;
2428    else
2429    {
2430      switch(event.type)
2431      {
2432        case EVENT_ID_SERVER_ATTRIBUT:
2433          recvDistributionAttributes(event);
2434          return true;
2435          break;
2436        case EVENT_ID_INDEX:
2437          recvIndex(event);
2438          return true;
2439          break;
2440        case EVENT_ID_MASK:
2441          recvMask(event);
2442          return true;
2443          break;
2444        case EVENT_ID_LON:
2445          recvLon(event);
2446          return true;
2447          break;
2448        case EVENT_ID_LAT:
2449          recvLat(event);
2450          return true;
2451          break;
2452        case EVENT_ID_AREA:
2453          recvArea(event);
2454          return true;
2455          break; 
2456        case EVENT_ID_DATA_INDEX:
2457          recvDataIndex(event);
2458          return true;
2459          break;
2460        default:
2461          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2462                << "Unknown Event");
2463          return false;
2464       }
2465    }
2466  }
2467
2468  /*!
2469    Receive index event from clients(s)
2470    \param[in] event event contain info about rank and associated index
2471  */
2472  void CDomain::recvIndex(CEventServer& event)
2473  {
2474    string domainId;
2475    std::map<int, CBufferIn*> rankBuffers;
2476
2477    list<CEventServer::SSubEvent>::iterator it;
2478    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2479    {     
2480      CBufferIn* buffer = it->buffer;
2481      *buffer >> domainId;
2482      rankBuffers[it->rank] = buffer;       
2483    }
2484    get(domainId)->recvIndex(rankBuffers);
2485  }
2486
2487  /*!
2488    Receive index information from client(s). We use the global index for mapping index between
2489    sending clients and receiving clients.
2490    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2491  */
2492  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2493  {
2494    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2495    recvClientRanks_.resize(nbReceived);       
2496
2497    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2498    ind = 0;
2499    for (ind = 0; it != ite; ++it, ++ind)
2500    {       
2501       recvClientRanks_[ind] = it->first;
2502       CBufferIn& buffer = *(it->second);
2503       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2504       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2505    }
2506    int nbIndGlob = 0;
2507    for (i = 0; i < nbReceived; ++i)
2508    {
2509      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2510    }
2511   
2512    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2513    i_index.resize(nbIndGlob);
2514    j_index.resize(nbIndGlob);   
2515    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2516
2517    nbIndGlob = 0;
2518    for (i = 0; i < nbReceived; ++i)
2519    {
2520      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2521      for (ind = 0; ind < tmp.numElements(); ++ind)
2522      {
2523         index = tmp(ind);
2524         if (0 == globalLocalIndexMap_.count(index))
2525         {
2526           iIndex = (index%ni_glo)-ibegin;
2527           iIndex = (iIndex < 0) ? 0 : iIndex;
2528           jIndex = (index/ni_glo)-jbegin;
2529           jIndex = (jIndex < 0) ? 0 : jIndex;
2530           nbIndLoc = iIndex + ni * jIndex;
2531           if (nbIndLoc < nbIndexGlobMax)
2532           {
2533             i_index(nbIndLoc) = index % ni_glo;
2534             j_index(nbIndLoc) = index / ni_glo;
2535             globalLocalIndexMap_[index] = nbIndLoc; 
2536             ++nbIndGlob;
2537           }
2538           // i_index(nbIndGlob) = index % ni_glo;
2539           // j_index(nbIndGlob) = index / ni_glo;
2540           // globalLocalIndexMap_[index] = nbIndGlob; 
2541           // ++nbIndGlob;
2542         } 
2543      } 
2544    } 
2545
2546    if (nbIndGlob==0)
2547    {
2548      i_index.resize(nbIndGlob);
2549      j_index.resize(nbIndGlob);
2550    }
2551    else
2552    {
2553      i_index.resizeAndPreserve(nbIndGlob);
2554      j_index.resizeAndPreserve(nbIndGlob);
2555    }
2556  }
2557
2558  /*!
2559    Receive attributes event from clients(s)
2560    \param[in] event event contain info about rank and associated attributes
2561  */
2562  void CDomain::recvDistributionAttributes(CEventServer& event)
2563  {
2564    CBufferIn* buffer=event.subEvents.begin()->buffer;
2565    string domainId ;
2566    *buffer>>domainId ;
2567    get(domainId)->recvDistributionAttributes(*buffer);
2568  }
2569
2570  /*!
2571    Receive attributes from client(s): zoom info and begin and n of each server
2572    \param[in] rank rank of client source
2573    \param[in] buffer message containing attributes info
2574  */
2575  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2576  {
2577    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2578    int global_zoom_ni_tmp, global_zoom_ibegin_tmp, global_zoom_nj_tmp, global_zoom_jbegin_tmp;
2579    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2580           >> global_zoom_ni_tmp >> global_zoom_ibegin_tmp >> global_zoom_nj_tmp >> global_zoom_jbegin_tmp           
2581           >> isCompressible_;
2582    ni.setValue(ni_tmp);
2583    ibegin.setValue(ibegin_tmp);
2584    nj.setValue(nj_tmp);
2585    jbegin.setValue(jbegin_tmp);
2586
2587    global_zoom_ni.setValue(global_zoom_ni_tmp);
2588    global_zoom_ibegin.setValue(global_zoom_ibegin_tmp);
2589    global_zoom_nj.setValue(global_zoom_nj_tmp);
2590    global_zoom_jbegin.setValue(global_zoom_jbegin_tmp);
2591
2592    int iend = ibegin + ni  - 1;
2593    int jend = jbegin + nj  - 1;
2594    int zoom_iend_glob = global_zoom_ibegin + global_zoom_ni - 1;
2595    int zoom_jend_glob = global_zoom_jbegin + global_zoom_nj - 1;
2596
2597    zoom_ibegin.setValue(global_zoom_ibegin > ibegin ? global_zoom_ibegin : ibegin);
2598    int zoom_iend = zoom_iend_glob < iend ? zoom_iend_glob : iend ;
2599    zoom_ni.setValue(zoom_iend-zoom_ibegin+1);
2600
2601    zoom_jbegin.setValue(global_zoom_jbegin > jbegin ? global_zoom_jbegin : jbegin);
2602    int zoom_jend = zoom_jend_glob < jend ? zoom_jend_glob : jend ;
2603    zoom_nj.setValue(zoom_jend-zoom_jbegin+1);
2604
2605    if (zoom_ni<=0 || zoom_nj<=0)
2606    {
2607      zoom_ni=0 ; zoom_ibegin=global_zoom_ibegin ; //=0; zoom_iend=0 ;
2608      zoom_nj=0 ; zoom_jbegin=global_zoom_jbegin ; //=0; zoom_jend=0 ;
2609    }
2610
2611  }
2612
2613  /*!
2614    Receive area event from clients(s)
2615    \param[in] event event contain info about rank and associated area
2616  */
2617  void CDomain::recvMask(CEventServer& event)
2618  {
2619    string domainId;
2620    std::map<int, CBufferIn*> rankBuffers;
2621
2622    list<CEventServer::SSubEvent>::iterator it;
2623    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2624    {     
2625      CBufferIn* buffer = it->buffer;
2626      *buffer >> domainId;
2627      rankBuffers[it->rank] = buffer;     
2628    }
2629    get(domainId)->recvMask(rankBuffers);
2630  }
2631
2632
2633  /*!
2634    Receive mask information from client(s)
2635    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2636  */
2637  void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)
2638  {
2639    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2640    if (nbReceived != recvClientRanks_.size())
2641      ERROR("void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)",
2642           << "The number of sending clients is not correct."
2643           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2644
2645    vector<CArray<bool,1> > recvMaskValue(nbReceived);     
2646    for (i = 0; i < recvClientRanks_.size(); ++i)
2647    {
2648      int rank = recvClientRanks_[i];
2649      CBufferIn& buffer = *(rankBuffers[rank]);     
2650      buffer >> recvMaskValue[i];
2651    }
2652
2653    int nbMaskInd = 0;
2654    for (i = 0; i < nbReceived; ++i)
2655    {
2656      nbMaskInd += recvMaskValue[i].numElements();
2657    }
2658 
2659    if (nbMaskInd != globalLocalIndexMap_.size())
2660      info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2661               << "Something must be wrong with mask index "<< std::endl;
2662
2663    nbMaskInd = globalLocalIndexMap_.size();
2664    mask_1d.resize(nbMaskInd);
2665    domainMask.resize(nbMaskInd);
2666    mask_1d = false;
2667   
2668    for (i = 0; i < nbReceived; ++i)
2669    {
2670      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2671      CArray<bool,1>& tmp = recvMaskValue[i];
2672      for (ind = 0; ind < tmp.numElements(); ++ind)
2673      {
2674        lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2675        if (!mask_1d(lInd)) // Only rewrite mask_1d if it's not true
2676          mask_1d(lInd) = tmp(ind);
2677      }
2678    }
2679    domainMask=mask_1d ;
2680  }
2681
2682  /*!
2683    Receive longitude event from clients(s)
2684    \param[in] event event contain info about rank and associated longitude
2685  */
2686  void CDomain::recvLon(CEventServer& event)
2687  {
2688    string domainId;
2689    std::map<int, CBufferIn*> rankBuffers;
2690
2691    list<CEventServer::SSubEvent>::iterator it;
2692    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2693    {     
2694      CBufferIn* buffer = it->buffer;
2695      *buffer >> domainId;
2696      rankBuffers[it->rank] = buffer;       
2697    }
2698    get(domainId)->recvLon(rankBuffers);
2699  }
2700
2701  /*!
2702    Receive longitude information from client(s)
2703    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2704  */
2705  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2706  {
2707    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2708    if (nbReceived != recvClientRanks_.size())
2709      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2710           << "The number of sending clients is not correct."
2711           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2712
2713    vector<CArray<double,1> > recvLonValue(nbReceived);
2714    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2715    for (i = 0; i < recvClientRanks_.size(); ++i)
2716    {
2717      int rank = recvClientRanks_[i];
2718      CBufferIn& buffer = *(rankBuffers[rank]);
2719      buffer >> hasLonLat;
2720      if (hasLonLat)
2721        buffer >> recvLonValue[i];
2722      buffer >> hasBounds;
2723      if (hasBounds)
2724        buffer >> recvBoundsLonValue[i];
2725    }
2726
2727    if (hasLonLat)
2728    {
2729      int nbLonInd = 0;
2730      for (i = 0; i < nbReceived; ++i)
2731      {
2732        nbLonInd += recvLonValue[i].numElements();
2733      }
2734   
2735      if (nbLonInd != globalLocalIndexMap_.size())
2736        info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2737                 << "Something must be wrong with longitude index "<< std::endl;
2738
2739      nbLonInd = globalLocalIndexMap_.size();
2740      lonvalue.resize(nbLonInd);
2741      if (hasBounds)
2742      {
2743        bounds_lonvalue.resize(nvertex,nbLonInd);
2744        bounds_lonvalue = 0.;
2745      }
2746
2747      nbLonInd = 0;
2748      for (i = 0; i < nbReceived; ++i)
2749      {
2750        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2751        CArray<double,1>& tmp = recvLonValue[i];
2752        for (ind = 0; ind < tmp.numElements(); ++ind)
2753        {
2754          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2755          lonvalue(lInd) = tmp(ind); 
2756           if (hasBounds)
2757           {         
2758            for (int nv = 0; nv < nvertex; ++nv)
2759              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2760           }                 
2761        }
2762      }       
2763    }
2764  }
2765
2766  /*!
2767    Receive latitude event from clients(s)
2768    \param[in] event event contain info about rank and associated latitude
2769  */
2770  void CDomain::recvLat(CEventServer& event)
2771  {
2772    string domainId;
2773    std::map<int, CBufferIn*> rankBuffers;
2774
2775    list<CEventServer::SSubEvent>::iterator it;
2776    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2777    {     
2778      CBufferIn* buffer = it->buffer;
2779      *buffer >> domainId;
2780      rankBuffers[it->rank] = buffer;   
2781    }
2782    get(domainId)->recvLat(rankBuffers);
2783  }
2784
2785  /*!
2786    Receive latitude information from client(s)
2787    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2788  */
2789  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2790  {
2791    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2792    if (nbReceived != recvClientRanks_.size())
2793      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2794           << "The number of sending clients is not correct."
2795           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2796
2797    vector<CArray<double,1> > recvLatValue(nbReceived);
2798    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2799    for (i = 0; i < recvClientRanks_.size(); ++i)
2800    {
2801      int rank = recvClientRanks_[i];
2802      CBufferIn& buffer = *(rankBuffers[rank]);
2803      buffer >> hasLonLat;
2804      if (hasLonLat)
2805        buffer >> recvLatValue[i];
2806      buffer >> hasBounds;
2807      if (hasBounds)
2808        buffer >> recvBoundsLatValue[i];
2809    }
2810
2811    if (hasLonLat)
2812    {
2813      int nbLatInd = 0;
2814      for (i = 0; i < nbReceived; ++i)
2815      {
2816        nbLatInd += recvLatValue[i].numElements();
2817      }
2818   
2819      if (nbLatInd != globalLocalIndexMap_.size())
2820        info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2821                << "Something must be wrong with latitude index "<< std::endl;
2822
2823      nbLatInd = globalLocalIndexMap_.size();
2824      latvalue.resize(nbLatInd);
2825      if (hasBounds)
2826      {
2827        bounds_latvalue.resize(nvertex,nbLatInd);
2828        bounds_latvalue = 0. ;
2829      }
2830
2831      nbLatInd = 0;
2832      for (i = 0; i < nbReceived; ++i)
2833      {
2834        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2835        CArray<double,1>& tmp = recvLatValue[i];
2836        for (ind = 0; ind < tmp.numElements(); ++ind)
2837        {
2838          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2839          latvalue(lInd) = tmp(ind);   
2840           if (hasBounds)
2841           {
2842            CArray<double,2>& boundslat = recvBoundsLatValue[i];
2843            for (int nv = 0; nv < nvertex; ++nv)
2844              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2845           }   
2846          ++nbLatInd;
2847        }
2848      }       
2849    }
2850  }
2851
2852  /*!
2853    Receive area event from clients(s)
2854    \param[in] event event contain info about rank and associated area
2855  */
2856  void CDomain::recvArea(CEventServer& event)
2857  {
2858    string domainId;
2859    std::map<int, CBufferIn*> rankBuffers;
2860
2861    list<CEventServer::SSubEvent>::iterator it;
2862    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2863    {     
2864      CBufferIn* buffer = it->buffer;
2865      *buffer >> domainId;
2866      rankBuffers[it->rank] = buffer;     
2867    }
2868    get(domainId)->recvArea(rankBuffers);
2869  }
2870
2871  /*!
2872    Receive area information from client(s)
2873    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2874  */
2875  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
2876  {
2877    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2878    if (nbReceived != recvClientRanks_.size())
2879      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2880           << "The number of sending clients is not correct."
2881           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2882
2883    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2884    for (i = 0; i < recvClientRanks_.size(); ++i)
2885    {
2886      int rank = recvClientRanks_[i];
2887      CBufferIn& buffer = *(rankBuffers[rank]);     
2888      buffer >> hasArea;
2889      if (hasArea)
2890        buffer >> recvAreaValue[i];
2891    }
2892
2893    if (hasArea)
2894    {
2895      int nbAreaInd = 0;
2896      for (i = 0; i < nbReceived; ++i)
2897      {     
2898        nbAreaInd += recvAreaValue[i].numElements();
2899      }
2900
2901      if (nbAreaInd != globalLocalIndexMap_.size())
2902        info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2903                 << "Something must be wrong with area index "<< std::endl;
2904
2905      nbAreaInd = globalLocalIndexMap_.size();
2906      areavalue.resize(nbAreaInd);
2907      nbAreaInd = 0;     
2908      for (i = 0; i < nbReceived; ++i)
2909      {
2910        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2911        CArray<double,1>& tmp = recvAreaValue[i];
2912        for (ind = 0; ind < tmp.numElements(); ++ind)
2913        {
2914          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2915          areavalue(lInd) = tmp(ind);         
2916        }
2917      }
2918     
2919    }
2920  }
2921
2922  /*!
2923    Compare two domain objects.
2924    They are equal if only if they have identical attributes as well as their values.
2925    Moreover, they must have the same transformations.
2926  \param [in] domain Compared domain
2927  \return result of the comparison
2928  */
2929  bool CDomain::isEqual(CDomain* obj)
2930  {
2931    vector<StdString> excludedAttr;
2932    excludedAttr.push_back("domain_ref");
2933    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2934    if (!objEqual) return objEqual;
2935
2936    TransMapTypes thisTrans = this->getAllTransformations();
2937    TransMapTypes objTrans  = obj->getAllTransformations();
2938
2939    TransMapTypes::const_iterator it, itb, ite;
2940    std::vector<ETranformationType> thisTransType, objTransType;
2941    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2942      thisTransType.push_back(it->first);
2943    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2944      objTransType.push_back(it->first);
2945
2946    if (thisTransType.size() != objTransType.size()) return false;
2947    for (int idx = 0; idx < thisTransType.size(); ++idx)
2948      objEqual &= (thisTransType[idx] == objTransType[idx]);
2949
2950    return objEqual;
2951  }
2952
2953  /*!
2954    Receive data index event from clients(s)
2955    \param[in] event event contain info about rank and associated index
2956  */
2957  void CDomain::recvDataIndex(CEventServer& event)
2958  {
2959    string domainId;
2960    std::map<int, CBufferIn*> rankBuffers;
2961
2962    list<CEventServer::SSubEvent>::iterator it;
2963    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2964    {     
2965      CBufferIn* buffer = it->buffer;
2966      *buffer >> domainId;
2967      rankBuffers[it->rank] = buffer;       
2968    }
2969    get(domainId)->recvDataIndex(rankBuffers);
2970  }
2971
2972  /*!
2973    Receive data index information from client(s)
2974    A client receives data index from different clients to rebuild its own data index.
2975    Because we use global index + mask info to calculate the sending data to client(s),
2976    this data index must be updated with mask info (maybe it will change in the future)
2977    Because the data index is local, to rebuild data index of received client, we should use global index along with.
2978
2979    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2980  */
2981  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
2982  {
2983    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
2984    if (nbReceived != recvClientRanks_.size())
2985      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
2986           << "The number of sending clients is not correct."
2987           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2988
2989    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
2990    for (i = 0; i < recvClientRanks_.size(); ++i)
2991    {
2992      int rank = recvClientRanks_[i];
2993      CBufferIn& buffer = *(rankBuffers[rank]);
2994      buffer >> recvDataIIndex[i];
2995      buffer >> recvDataJIndex[i];
2996    }
2997   
2998    int nbIndex = i_index.numElements();
2999    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3000    dataIIndex = -1; dataJIndex = -1;
3001     
3002    nbIndex = 0;
3003    for (i = 0; i < nbReceived; ++i)
3004    {     
3005      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3006      CArray<int,1>& tmpI = recvDataIIndex[i];   
3007      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3008      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3009          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3010             << "The number of global received index is not coherent with the number of received data index."
3011             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3012
3013      for (ind = 0; ind < tmpI.numElements(); ++ind)
3014      {
3015         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3016         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3017         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3018
3019         if (!domainMask(lInd))   // Include mask info into data index on the RECEIVE getServerDimensionSizes   
3020         {
3021           dataIIndex(lInd) = dataJIndex(lInd) = -1;
3022         }
3023      } 
3024    }
3025
3026    int nbCompressedData = 0; 
3027    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3028    {
3029       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3030       if ((0 <= indexI) && (0 <= indexJ))
3031         ++nbCompressedData;
3032    }       
3033 
3034    data_i_index.resize(nbCompressedData);
3035    data_j_index.resize(nbCompressedData);
3036
3037    nbCompressedData = 0; 
3038    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3039    {
3040       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3041       if ((0 <= indexI) && (0 <= indexJ))
3042       {
3043          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3044          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3045         ++nbCompressedData;
3046       }
3047    }
3048
3049    // Reset data_ibegin, data_jbegin
3050    data_ibegin.setValue(0);
3051    data_jbegin.setValue(0);
3052  }
3053
3054  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3055  {
3056    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3057    return transformationMap_.back().second;
3058  }
3059
3060  /*!
3061    Check whether a domain has transformation
3062    \return true if domain has transformation
3063  */
3064  bool CDomain::hasTransformation()
3065  {
3066    return (!transformationMap_.empty());
3067  }
3068
3069  /*!
3070    Set transformation for current domain. It's the method to move transformation in hierarchy
3071    \param [in] domTrans transformation on domain
3072  */
3073  void CDomain::setTransformations(const TransMapTypes& domTrans)
3074  {
3075    transformationMap_ = domTrans;
3076  }
3077
3078  /*!
3079    Get all transformation current domain has
3080    \return all transformation
3081  */
3082  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3083  {
3084    return transformationMap_;
3085  }
3086
3087  void CDomain::duplicateTransformation(CDomain* src)
3088  {
3089    if (src->hasTransformation())
3090    {
3091      this->setTransformations(src->getAllTransformations());
3092    }
3093  }
3094
3095  /*!
3096   * Go through the hierarchy to find the domain from which the transformations must be inherited
3097   */
3098  void CDomain::solveInheritanceTransformation()
3099  {
3100    if (hasTransformation() || !hasDirectDomainReference())
3101      return;
3102
3103    CDomain* domain = this;
3104    std::vector<CDomain*> refDomains;
3105    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3106    {
3107      refDomains.push_back(domain);
3108      domain = domain->getDirectDomainReference();
3109    }
3110
3111    if (domain->hasTransformation())
3112      for (size_t i = 0; i < refDomains.size(); ++i)
3113        refDomains[i]->setTransformations(domain->getAllTransformations());
3114  }
3115
3116  void CDomain::setContextClient(CContextClient* contextClient)
3117  {
3118    if (clientsSet.find(contextClient)==clientsSet.end())
3119    {
3120      clients.push_back(contextClient) ;
3121      clientsSet.insert(contextClient);
3122    }
3123  }
3124
3125  /*!
3126    Parse children nodes of a domain in xml file.
3127    Whenver there is a new transformation, its type and name should be added into this function
3128    \param node child node to process
3129  */
3130  void CDomain::parse(xml::CXMLNode & node)
3131  {
3132    SuperClass::parse(node);
3133
3134    if (node.goToChildElement())
3135    {
3136      StdString nodeElementName;
3137      do
3138      {
3139        StdString nodeId("");
3140        if (node.getAttributes().end() != node.getAttributes().find("id"))
3141        { nodeId = node.getAttributes()["id"]; }
3142
3143        nodeElementName = node.getElementName();
3144        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3145        it = transformationMapList_.find(nodeElementName);
3146        if (ite != it)
3147        {
3148          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3149                                                                                                                nodeId,
3150                                                                                                                &node)));
3151        }
3152        else
3153        {
3154          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3155                << "The transformation " << nodeElementName << " has not been supported yet.");
3156        }
3157      } while (node.goToNextElement()) ;
3158      node.goToParentElement();
3159    }
3160  }
3161   //----------------------------------------------------------------
3162
3163   DEFINE_REF_FUNC(Domain,domain)
3164
3165   ///---------------------------------------------------------------
3166
3167} // namespace xios
Note: See TracBrowser for help on using the repository browser.