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

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

Replacing Boost's unordered_map and shared_pointer by its STL counterparts.

Two notes for Curie:

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