source: XIOS/trunk/src/node/domain.cpp @ 1549

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

Adding a new domain transformation "domain_extract".

For now it has only the functionalities of zoom.

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