source: XIOS/dev/dev_olga/src/node/domain.cpp @ 1568

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

Grid mask is now applied in the source filter of clients: unmasked values are replaced by NaN. It is not reconstructed any more by servers.

This needs to be tested more rigorously before commiting to trunk.

  • 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: 109.1 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)
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)
45      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
46      , globalLocalIndexMap_(), computedWrittenIndex_(false)
47      , clients()
48   {
49    }
50
51   CDomain::~CDomain(void)
52   {
53   }
54
55   ///---------------------------------------------------------------
56
57   void CDomain::assignMesh(const StdString meshName, const int nvertex)
58   {
59     mesh = CMesh::getMesh(meshName, nvertex);
60   }
61
62   CDomain* CDomain::createDomain()
63   {
64     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
65     return domain;
66   }
67
68   std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
69   bool CDomain::_dummyTransformationMapList = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
70
71   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
72   {
73     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
74     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
75     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
76     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
77     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
78     m["reorder_domain"] = TRANS_REORDER_DOMAIN;
79     m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
80   }
81
82   const std::set<StdString> & CDomain::getRelFiles(void) const
83   {
84      return (this->relFiles);
85   }
86
87
88   /*!
89     Returns the number of indexes written by each server.
90     \return the number of indexes written by each server
91   */
92   int CDomain::getNumberWrittenIndexes(MPI_Comm writtenCom)
93   {
94     int writtenSize;
95     MPI_Comm_size(writtenCom, &writtenSize);
96     return numberWrittenIndexes_[writtenSize];
97   }
98
99   /*!
100     Returns the total number of indexes written by the servers.
101     \return the total number of indexes written by the servers
102   */
103   int CDomain::getTotalNumberWrittenIndexes(MPI_Comm writtenCom)
104   {
105     int writtenSize;
106     MPI_Comm_size(writtenCom, &writtenSize);
107     return totalNumberWrittenIndexes_[writtenSize];
108   }
109
110   /*!
111     Returns the offset of indexes written by each server.
112     \return the offset of indexes written by each server
113   */
114   int CDomain::getOffsetWrittenIndexes(MPI_Comm writtenCom)
115   {
116     int writtenSize;
117     MPI_Comm_size(writtenCom, &writtenSize);
118     return offsetWrittenIndexes_[writtenSize];
119   }
120
121   CArray<int, 1>& CDomain::getCompressedIndexToWriteOnServer(MPI_Comm writtenCom)
122   {
123     int writtenSize;
124     MPI_Comm_size(writtenCom, &writtenSize);
125     return compressedIndexToWriteOnServer[writtenSize];
126   }
127
128   //----------------------------------------------------------------
129
130   /*!
131    * Compute the minimum buffer size required to send the attributes to the server(s).
132    *
133    * \return A map associating the server rank with its minimum buffer size.
134    */
135   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
136   {
137
138     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
139
140     if (client->isServerLeader())
141     {
142       // size estimation for sendDistributionAttribut
143       size_t size = 11 * sizeof(size_t);
144
145       const std::list<int>& ranks = client->getRanksServerLeader();
146       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
147       {
148         if (size > attributesSizes[*itRank])
149           attributesSizes[*itRank] = size;
150       }
151     }
152
153     std::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->serverSize].end();
154     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
155     for (size_t k = 0; k < connectedServerRank_[client->serverSize].size(); ++k)
156     {
157       int rank = connectedServerRank_[client->serverSize][k];
158       std::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->serverSize].find(rank);
159       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
160
161       // size estimation for sendIndex (and sendArea which is always smaller or equal)
162       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
163       // if (isCompressible_)
164       // {
165       //   std::map<int, std::vector<int> >::const_iterator itWritten = indWrittenSrv_.find(rank);
166       //   size_t writtenIdxCount = (itWritten != itWrittenIndexEnd) ? itWritten->second.size() : 0;
167       //   sizeIndexEvent += CArray<int,1>::size(writtenIdxCount);
168       // }
169
170       // size estimation for sendLonLat
171       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
172       if (hasBounds)
173         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
174
175       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
176       if (size > attributesSizes[rank])
177         attributesSizes[rank] = size;
178     }
179
180     return attributesSizes;
181   }
182
183   //----------------------------------------------------------------
184
185   bool CDomain::isEmpty(void) const
186   {
187     return ((this->i_index.isEmpty()) || (0 == this->i_index.numElements()));
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() && lonvalue_2d.isEmpty() && lonvalue_1d.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() && latvalue_2d.isEmpty() && latvalue_1d.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 2D longitude and latitude of curvilinear domain read from a file.
530      If there are already longitude and latitude defined by model. We just ignore read value.
531    */
532   void CDomain::fillInCurvilinearLonLat()
533   {
534     if (!lonvalue_curvilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.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() && latvalue_1d.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() && bounds_lon_1d.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() && bounds_lat_1d.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
869   size_t CDomain::getGlobalWrittenSize(void)
870   {
871     return ni_glo*nj_glo ;
872   }
873   //----------------------------------------------------------------
874
875   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
876   void CDomain::checkLocalIDomain(void)
877   {
878      // If ibegin and ni are provided then we use them to check the validity of local domain
879      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
880      {
881        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
882        {
883          ERROR("CDomain::checkLocalIDomain(void)",
884                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
885                << "The local domain is wrongly defined,"
886                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
887        }
888      }
889
890      // i_index has higher priority than ibegin and ni
891      if (!i_index.isEmpty())
892      {
893        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
894        if (ni.isEmpty()) 
895        {         
896         // No information about ni
897          int minIndex = ni_glo - 1;
898          int maxIndex = 0;
899          for (int idx = 0; idx < i_index.numElements(); ++idx)
900          {
901            if (i_index(idx) < minIndex) minIndex = i_index(idx);
902            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
903          }
904          ni = maxIndex - minIndex + 1; 
905          minIIndex = minIIndex;         
906        }
907
908        // It's not so correct but if ibegin is not the first value of i_index
909        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
910        if (ibegin.isEmpty()) ibegin = minIIndex;
911      }
912      else if (ibegin.isEmpty() && ni.isEmpty())
913      {
914        ibegin = 0;
915        ni = ni_glo;
916      }
917      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
918      {
919        ERROR("CDomain::checkLocalIDomain(void)",
920              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
921              << "The local domain is wrongly defined," << endl
922              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
923              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
924      }
925       
926
927      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
928      {
929        ERROR("CDomain::checkLocalIDomain(void)",
930              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
931              << "The local domain is wrongly defined,"
932              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
933      }
934   }
935
936   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
937   void CDomain::checkLocalJDomain(void)
938   {
939    // If jbegin and nj are provided then we use them to check the validity of local domain
940     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
941     {
942       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
943       {
944         ERROR("CDomain::checkLocalJDomain(void)",
945                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
946                << "The local domain is wrongly defined,"
947                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
948       }
949     }
950
951     if (!j_index.isEmpty())
952     {
953        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
954        if (nj.isEmpty()) 
955        {
956          // No information about nj
957          int minIndex = nj_glo - 1;
958          int maxIndex = 0;
959          for (int idx = 0; idx < j_index.numElements(); ++idx)
960          {
961            if (j_index(idx) < minIndex) minIndex = j_index(idx);
962            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
963          }
964          nj = maxIndex - minIndex + 1;
965          minJIndex = minIndex; 
966        } 
967        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
968        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
969       if (jbegin.isEmpty()) jbegin = minJIndex;       
970     }
971     else if (jbegin.isEmpty() && nj.isEmpty())
972     {
973       jbegin = 0;
974       nj = nj_glo;
975     }     
976
977
978     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
979     {
980       ERROR("CDomain::checkLocalJDomain(void)",
981              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
982              << "The local domain is wrongly defined,"
983              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
984     }
985   }
986
987   //----------------------------------------------------------------
988
989   void CDomain::checkMask(void)
990   {
991      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
992        ERROR("CDomain::checkMask(void)",
993              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
994              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
995              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
996
997      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
998      {
999        if (mask_1d.numElements() != i_index.numElements())
1000          ERROR("CDomain::checkMask(void)",
1001                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1002                << "'mask_1d' does not have the same size as the local domain." << std::endl
1003                << "Local size is " << i_index.numElements() << "." << std::endl
1004                << "Mask size is " << mask_1d.numElements() << ".");
1005      }
1006
1007      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
1008      {
1009        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1010          ERROR("CDomain::checkMask(void)",
1011                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1012                << "The mask does not have the same size as the local domain." << std::endl
1013                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1014                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
1015      }
1016
1017      if (!mask_2d.isEmpty())
1018      {
1019        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
1020        for (int j = 0; j < nj; ++j)
1021          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1022//        mask_2d.reset();
1023      }
1024      else if (mask_1d.isEmpty())
1025      {
1026        domainMask.resize(i_index.numElements());
1027        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
1028      }
1029      else
1030      {
1031      domainMask.resize(mask_1d.numElements());
1032      domainMask=mask_1d ;
1033     }
1034   }
1035
1036   //----------------------------------------------------------------
1037
1038   void CDomain::checkDomainData(void)
1039   {
1040      if (data_dim.isEmpty())
1041      {
1042        data_dim.setValue(1);
1043      }
1044      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
1045      {
1046        ERROR("CDomain::checkDomainData(void)",
1047              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1048              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
1049      }
1050
1051      if (data_ibegin.isEmpty())
1052         data_ibegin.setValue(0);
1053      if (data_jbegin.isEmpty())
1054         data_jbegin.setValue(0);
1055
1056      if (data_ni.isEmpty())
1057      {
1058        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
1059      }
1060      else if (data_ni.getValue() < 0)
1061      {
1062        ERROR("CDomain::checkDomainData(void)",
1063              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1064              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
1065      }
1066
1067      if (data_nj.isEmpty())
1068      {
1069        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
1070      }
1071      else if (data_nj.getValue() < 0)
1072      {
1073        ERROR("CDomain::checkDomainData(void)",
1074              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1075              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1076      }
1077   }
1078
1079   //----------------------------------------------------------------
1080
1081   void CDomain::checkCompression(void)
1082   {
1083      if (!data_i_index.isEmpty())
1084      {
1085        if (!data_j_index.isEmpty() &&
1086            data_j_index.numElements() != data_i_index.numElements())
1087        {
1088           ERROR("CDomain::checkCompression(void)",
1089                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1090                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1091                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1092                 << "'data_j_index' size = " << data_j_index.numElements());
1093        }
1094
1095        if (2 == data_dim)
1096        {
1097          if (data_j_index.isEmpty())
1098          {
1099             ERROR("CDomain::checkCompression(void)",
1100                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1101                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
1102          }
1103        }
1104        else // (1 == data_dim)
1105        {
1106          if (data_j_index.isEmpty())
1107          {
1108            data_j_index.resize(data_ni);
1109            for (int j = 0; j < data_ni; ++j) data_j_index(j) = 0;
1110          }
1111        }
1112      }
1113      else
1114      {
1115        if (data_dim == 2 && !data_j_index.isEmpty())
1116          ERROR("CDomain::checkCompression(void)",
1117                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1118                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
1119
1120        if (1 == data_dim)
1121        {
1122          data_i_index.resize(data_ni);
1123          data_j_index.resize(data_ni);
1124
1125          for (int i = 0; i < data_ni; ++i)
1126          {
1127            data_i_index(i) = i;
1128            data_j_index(i) = 0;
1129          }
1130        }
1131        else // (data_dim == 2)
1132        {
1133          const int dsize = data_ni * data_nj;
1134          data_i_index.resize(dsize);
1135          data_j_index.resize(dsize);
1136
1137          for(int count = 0, j = 0; j < data_nj; ++j)
1138          {
1139            for(int i = 0; i < data_ni; ++i, ++count)
1140            {
1141              data_i_index(count) = i;
1142              data_j_index(count) = j;
1143            }
1144          }
1145        }
1146      }
1147   }
1148
1149   //----------------------------------------------------------------
1150   void CDomain::computeLocalMask(void)
1151   {
1152     localMask.resize(i_index.numElements()) ;
1153     localMask=false ;
1154
1155     size_t dn=data_i_index.numElements() ;
1156     int i,j ;
1157     size_t k,ind ;
1158
1159     for(k=0;k<dn;k++)
1160     {
1161       if (data_dim==2)
1162       {
1163          i=data_i_index(k)+data_ibegin ;
1164          j=data_j_index(k)+data_jbegin ;
1165          if (i>=0 && i<ni && j>=0 && j<nj)
1166          {
1167            ind=j*ni+i ;
1168            localMask(ind)=domainMask(ind) ;
1169          }
1170       }
1171       else
1172       {
1173          i=data_i_index(k)+data_ibegin ;
1174          if (i>=0 && i<i_index.numElements())
1175          {
1176            ind=i ;
1177            localMask(ind)=domainMask(ind) ;
1178          }
1179       }
1180     }
1181   }
1182
1183   void CDomain::checkEligibilityForCompressedOutput(void)
1184   {
1185     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1186     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1187   }
1188
1189   //----------------------------------------------------------------
1190
1191   /*
1192     Fill in longitude, latitude, bounds, and area into internal values (lonvalue, latvalue, bounds_lonvalue, bounds_latvalue, areavalue)
1193     which will be used by XIOS.
1194   */
1195   void CDomain::completeLonLatClient(void)
1196   {
1197     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1198     checkBounds() ;
1199     checkArea() ;
1200
1201     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
1202     {
1203       lonvalue.resize(ni * nj);
1204       latvalue.resize(ni * nj);
1205       if (hasBounds)
1206       {
1207         bounds_lonvalue.resize(nvertex, ni * nj);
1208         bounds_latvalue.resize(nvertex, ni * nj);
1209       }
1210
1211       for (int j = 0; j < nj; ++j)
1212       {
1213         for (int i = 0; i < ni; ++i)
1214         {
1215           int k = j * ni + i;
1216
1217           lonvalue(k) = lonvalue_2d(i,j);
1218           latvalue(k) = latvalue_2d(i,j);
1219
1220           if (hasBounds)
1221           {
1222             for (int n = 0; n < nvertex; ++n)
1223             {
1224               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1225               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
1226             }
1227           }
1228         }
1229       }
1230     }
1231     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
1232     {
1233       if (type_attr::rectilinear == type)
1234       {
1235         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1236         {
1237           lonvalue.resize(ni * nj);
1238           latvalue.resize(ni * nj);
1239           if (hasBounds)
1240           {
1241             bounds_lonvalue.resize(nvertex, ni * nj);
1242             bounds_latvalue.resize(nvertex, ni * nj);
1243           }
1244
1245           for (int j = 0; j < nj; ++j)
1246           {
1247             for (int i = 0; i < ni; ++i)
1248             {
1249               int k = j * ni + i;
1250
1251               lonvalue(k) = lonvalue_1d(i);
1252               latvalue(k) = latvalue_1d(j);
1253
1254               if (hasBounds)
1255               {
1256                 for (int n = 0; n < nvertex; ++n)
1257                 {
1258                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1259                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1260                 }
1261               }
1262             }
1263           }
1264         }
1265         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1266         {
1267           lonvalue.reference(lonvalue_1d);
1268           latvalue.reference(latvalue_1d);
1269            if (hasBounds)
1270           {
1271             bounds_lonvalue.reference(bounds_lon_1d);
1272             bounds_latvalue.reference(bounds_lat_1d);
1273           }
1274         }
1275         else
1276           ERROR("CDomain::completeLonClient(void)",
1277                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1278                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1279                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1280                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1281                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1282                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1283       }
1284       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1285       {
1286         lonvalue.reference(lonvalue_1d);
1287         latvalue.reference(latvalue_1d);
1288         if (hasBounds)
1289         {
1290           bounds_lonvalue.reference(bounds_lon_1d);
1291           bounds_latvalue.reference(bounds_lat_1d);
1292         }
1293       }
1294     }
1295
1296     if (!area.isEmpty() && areavalue.isEmpty())
1297     {
1298        areavalue.resize(ni*nj);
1299       for (int j = 0; j < nj; ++j)
1300       {
1301         for (int i = 0; i < ni; ++i)
1302         {
1303           int k = j * ni + i;
1304           areavalue(k) = area(i,j);
1305         }
1306       }
1307     }
1308   }
1309
1310   /*
1311     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1312   */
1313   void CDomain::convertLonLatValue(void)
1314   {
1315     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1316     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1317     {
1318       lonvalue_2d.resize(ni,nj);
1319       latvalue_2d.resize(ni,nj);
1320       if (hasBounds)
1321       {
1322         bounds_lon_2d.resize(nvertex, ni, nj);
1323         bounds_lat_2d.resize(nvertex, ni, nj);
1324       }
1325
1326       for (int j = 0; j < nj; ++j)
1327       {
1328         for (int i = 0; i < ni; ++i)
1329         {
1330           int k = j * ni + i;
1331
1332           lonvalue_2d(i,j) = lonvalue(k);
1333           latvalue_2d(i,j) = latvalue(k);
1334
1335           if (hasBounds)
1336           {
1337             for (int n = 0; n < nvertex; ++n)
1338             {
1339               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1340               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1341             }
1342           }
1343         }
1344       }
1345     }
1346     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1347     {
1348       if (type_attr::rectilinear == type)
1349       {
1350         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1351         {
1352           lonvalue.resize(ni * nj);
1353           latvalue.resize(ni * nj);
1354           if (hasBounds)
1355           {
1356             bounds_lonvalue.resize(nvertex, ni * nj);
1357             bounds_latvalue.resize(nvertex, ni * nj);
1358           }
1359
1360           for (int j = 0; j < nj; ++j)
1361           {
1362             for (int i = 0; i < ni; ++i)
1363             {
1364               int k = j * ni + i;
1365
1366               lonvalue(k) = lonvalue_1d(i);
1367               latvalue(k) = latvalue_1d(j);
1368
1369               if (hasBounds)
1370               {
1371                 for (int n = 0; n < nvertex; ++n)
1372                 {
1373                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1374                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1375                 }
1376               }
1377             }
1378           }
1379         }
1380         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1381         {
1382           lonvalue.reference(lonvalue_1d);
1383           latvalue.reference(latvalue_1d);
1384            if (hasBounds)
1385           {
1386             bounds_lonvalue.reference(bounds_lon_1d);
1387             bounds_latvalue.reference(bounds_lat_1d);
1388           }
1389         }
1390         else
1391           ERROR("CDomain::completeLonClient(void)",
1392                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1393                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1394                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1395                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1396                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1397                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1398       }
1399       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1400       {
1401         lonvalue.reference(lonvalue_1d);
1402         latvalue.reference(latvalue_1d);
1403         if (hasBounds)
1404         {
1405           bounds_lonvalue.reference(bounds_lon_1d);
1406           bounds_latvalue.reference(bounds_lat_1d);
1407         }
1408       }
1409     }
1410   }
1411
1412
1413   void CDomain::checkBounds(void)
1414   {
1415     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1416     if (!nvertex.isEmpty() && nvertex > 0 && !hasBoundValues)
1417     {
1418       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1419         ERROR("CDomain::checkBounds(void)",
1420               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1421               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1422               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1423
1424       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1425         ERROR("CDomain::checkBounds(void)",
1426               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1427               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1428               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1429
1430       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1431       {
1432         ERROR("CDomain::checkBounds(void)",
1433               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1434               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1435               << "Please define either both attributes or none.");
1436       }
1437
1438       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1439       {
1440         ERROR("CDomain::checkBounds(void)",
1441               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1442               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1443               << "Please define either both attributes or none.");
1444       }
1445
1446       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1447         ERROR("CDomain::checkBounds(void)",
1448               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1449               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1450               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(0)
1451               << " but nvertex is " << nvertex.getValue() << ".");
1452
1453       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1454         ERROR("CDomain::checkBounds(void)",
1455               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1456               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1457               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(0)
1458               << " but nvertex is " << nvertex.getValue() << ".");
1459
1460       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1461         ERROR("CDomain::checkBounds(void)",
1462               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1463               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1464
1465       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1466         ERROR("CDomain::checkBounds(void)",
1467               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1468               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1469
1470       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1471         ERROR("CDomain::checkBounds(void)",
1472               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1473               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1474               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(0)
1475               << " but nvertex is " << nvertex.getValue() << ".");
1476
1477       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1478         ERROR("CDomain::checkBounds(void)",
1479               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1480               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1481               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(0)
1482               << " but nvertex is " << nvertex.getValue() << ".");
1483
1484       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1485         ERROR("CDomain::checkBounds(void)",
1486               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1487               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1488
1489       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1490         ERROR("CDomain::checkBounds(void)",
1491               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1492
1493       // In case of reading UGRID bounds values are not required
1494       hasBounds = (!bounds_lat_1d.isEmpty() || !bounds_lat_2d.isEmpty() );
1495     }
1496     else if (hasBoundValues)
1497     {
1498       hasBounds = true;       
1499     }
1500     else
1501     {
1502       hasBounds = false;
1503     }
1504   }
1505
1506   void CDomain::checkArea(void)
1507   {
1508     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
1509     hasArea = !area.isEmpty();
1510     if (hasArea && !hasAreaValue)
1511     {
1512       if (area.extent(0) != ni || area.extent(1) != nj)
1513       {
1514         ERROR("CDomain::checkArea(void)",
1515               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1516               << "The area does not have the same size as the local domain." << std::endl
1517               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1518               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1519       }
1520//       if (areavalue.isEmpty())
1521//       {
1522//          areavalue.resize(ni*nj);
1523//         for (int j = 0; j < nj; ++j)
1524//         {
1525//           for (int i = 0; i < ni; ++i)
1526//           {
1527//             int k = j * ni + i;
1528//             areavalue(k) = area(i,j);
1529//           }
1530//         }
1531//       }
1532     }
1533   }
1534
1535   void CDomain::checkLonLat()
1536   {
1537     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1538                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1539     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1540     if (hasLonLat && !hasLonLatValue)
1541     {
1542       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1543         ERROR("CDomain::checkLonLat()",
1544               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1545               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1546               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1547
1548       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1549       {
1550         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1551           ERROR("CDomain::checkLonLat()",
1552                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1553                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1554                 << "Local size is " << i_index.numElements() << "." << std::endl
1555                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1556       }
1557
1558       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1559       {
1560         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1561           ERROR("CDomain::checkLonLat()",
1562                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1563                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1564                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1565                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1566       }
1567
1568       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1569         ERROR("CDomain::checkLonLat()",
1570               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1571               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1572               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1573
1574       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1575       {
1576         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1577           ERROR("CDomain::checkLonLat()",
1578                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1579                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1580                 << "Local size is " << i_index.numElements() << "." << std::endl
1581                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1582       }
1583
1584       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1585       {
1586         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1587           ERROR("CDomain::checkLonLat()",
1588                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1589                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1590                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1591                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1592       }
1593     }
1594   }
1595
1596   void CDomain::checkAttributesOnClientAfterTransformation()
1597   {
1598     CContext* context=CContext::getCurrent() ;
1599
1600     if (this->isClientAfterTransformationChecked) return;
1601     if (context->hasClient)
1602     {
1603      this->computeConnectedClients();
1604       if (hasLonLat)
1605         if (!context->hasServer)
1606           this->completeLonLatClient();
1607     }
1608
1609     this->isClientAfterTransformationChecked = true;
1610   }
1611
1612   //----------------------------------------------------------------
1613   // Divide function checkAttributes into 2 seperate ones
1614   // This function only checks all attributes of current domain
1615   void CDomain::checkAttributesOnClient()
1616   {
1617     if (this->isClientChecked) return;
1618     CContext* context=CContext::getCurrent();
1619
1620      if (context->hasClient && !context->hasServer)
1621      {
1622        this->checkDomain();
1623        this->checkBounds();
1624        this->checkArea();
1625        this->checkLonLat();
1626      }
1627
1628      if (context->hasClient && !context->hasServer)
1629      { // Ct client uniquement
1630         this->checkMask();
1631         this->checkDomainData();
1632         this->checkCompression();
1633         this->computeLocalMask() ;
1634      }
1635      else
1636      { // Ct serveur uniquement
1637      }
1638
1639      this->isClientChecked = true;
1640   }
1641
1642   // Send all checked attributes to server
1643   void CDomain::sendCheckedAttributes()
1644   {
1645     if (!this->isClientChecked) checkAttributesOnClient();
1646     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
1647     CContext* context=CContext::getCurrent() ;
1648
1649     if (this->isChecked) return;
1650     if (context->hasClient)
1651     {
1652       sendAttributes();
1653     }
1654     this->isChecked = true;
1655   }
1656
1657   void CDomain::checkAttributes(void)
1658   {
1659      if (this->isChecked) return;
1660      CContext* context=CContext::getCurrent() ;
1661
1662      this->checkDomain();
1663      this->checkLonLat();
1664      this->checkBounds();
1665      this->checkArea();
1666
1667      if (context->hasClient)
1668      { // Ct client uniquement
1669         this->checkMask();
1670         this->checkDomainData();
1671         this->checkCompression();
1672         this->computeLocalMask() ;
1673
1674      }
1675      else
1676      { // Ct serveur uniquement
1677      }
1678
1679      if (context->hasClient)
1680      {
1681        this->computeConnectedClients();
1682        this->completeLonLatClient();
1683      }
1684
1685      this->isChecked = true;
1686   }
1687
1688  /*!
1689     Compute the connection of a client to other clients to determine which clients to send attributes to.
1690     The sending clients are supposed to already know the distribution of receiving clients (In simple cases, it's band)
1691     The connection among clients is calculated by using global index.
1692     A client connects to other clients which holds the same global index as it.     
1693  */
1694  void CDomain::computeConnectedClients()
1695  {
1696    CContext* context=CContext::getCurrent() ;
1697   
1698    // This line should be changed soon.
1699    int nbSrvPools = (context->hasServer) ? (context->hasClient ? context->clientPrimServer.size() : 0) : 1;
1700
1701    nbSenders.clear();
1702    connectedServerRank_.clear();
1703
1704    for (int p = 0; p < nbSrvPools; ++p)
1705    {
1706      CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[p] : context->client;
1707      int nbServer = client->serverSize;
1708      int nbClient = client->clientSize;
1709      int rank     = client->clientRank;
1710      bool doComputeGlobalIndexServer = true;
1711
1712      if (connectedServerRank_.find(nbServer) == connectedServerRank_.end())
1713      {
1714
1715        if (indSrv_.find(nbServer) == indSrv_.end())
1716        {
1717          int i,j,i_ind,j_ind, nbIndex=i_index.numElements();
1718          int globalIndexCount = i_index.numElements();
1719          // Fill in index
1720          CArray<size_t,1> globalIndexDomain(nbIndex);
1721          size_t globalIndex;
1722
1723          for (i = 0; i < nbIndex; ++i)
1724          {
1725            i_ind=i_index(i);
1726            j_ind=j_index(i);
1727            globalIndex = i_ind + j_ind * ni_glo;
1728            globalIndexDomain(i) = globalIndex;
1729          }
1730
1731          if (globalLocalIndexMap_.empty())
1732          {
1733            for (i = 0; i < nbIndex; ++i)
1734              globalLocalIndexMap_[globalIndexDomain(i)] = i;
1735          }
1736
1737          size_t globalSizeIndex = 1, indexBegin, indexEnd;
1738          int range, clientSize = client->clientSize;
1739          std::vector<int> nGlobDomain(2);
1740          nGlobDomain[0] = this->ni_glo;
1741          nGlobDomain[1] = this->nj_glo;
1742          for (int i = 0; i < nGlobDomain.size(); ++i) globalSizeIndex *= nGlobDomain[i];
1743          indexBegin = 0;
1744          if (globalSizeIndex <= clientSize)
1745          {
1746            indexBegin = rank%globalSizeIndex;
1747            indexEnd = indexBegin;
1748          }
1749          else
1750          {
1751            for (int i = 0; i < clientSize; ++i)
1752            {
1753              range = globalSizeIndex / clientSize;
1754              if (i < (globalSizeIndex%clientSize)) ++range;
1755              if (i == client->clientRank) break;
1756              indexBegin += range;
1757            }
1758            indexEnd = indexBegin + range - 1;
1759          }
1760
1761          // Even if servers have no index, they must received something from client
1762          // We only use several client to send "empty" message to these servers
1763          CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
1764          std::vector<int> serverZeroIndex;
1765          if (isUnstructed_) serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 0);
1766          else serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t&,size_t&>(indexBegin, indexEnd), 1);
1767
1768          std::list<int> serverZeroIndexLeader;
1769          std::list<int> serverZeroIndexNotLeader;
1770          CContextClient::computeLeader(client->clientRank, client->clientSize, serverZeroIndex.size(), serverZeroIndexLeader, serverZeroIndexNotLeader);
1771          for (std::list<int>::iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
1772            *it = serverZeroIndex[*it];
1773
1774          CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(), client->intraComm);
1775          clientServerMap->computeServerIndexMapping(globalIndexDomain, nbServer);
1776          CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
1777
1778          CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
1779                 ite = globalIndexDomainOnServer.end();
1780          indSrv_[nbServer].swap(globalIndexDomainOnServer);
1781          connectedServerRank_[nbServer].clear();
1782          for (it = indSrv_[nbServer].begin(); it != ite; ++it)
1783            connectedServerRank_[nbServer].push_back(it->first);
1784
1785          for (std::list<int>::const_iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
1786            connectedServerRank_[nbServer].push_back(*it);
1787
1788          // Even if a client has no index, it must connect to at least one server and
1789          // send an "empty" data to this server
1790          if (connectedServerRank_[nbServer].empty())
1791            connectedServerRank_[nbServer].push_back(client->clientRank % client->serverSize);
1792
1793          nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
1794          delete clientServerMap;
1795        }
1796      }
1797    }
1798  }
1799
1800   /*!
1801     Compute index to write data. We only write data on the zoomed region, therefore, there should
1802     be a map between the complete grid and the reduced grid where we write data.
1803     By using global index we can easily create this kind of mapping.
1804   */
1805   void CDomain::computeWrittenIndex()
1806   { 
1807      if (computedWrittenIndex_) return;
1808      computedWrittenIndex_ = true;
1809
1810      CContext* context=CContext::getCurrent();     
1811      CContextServer* server = context->server; 
1812
1813      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1814      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
1815      nSize[0]        = ni;      nSize[1]  = nj;
1816      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1817      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1818      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
1819      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1820
1821      size_t nbWritten = 0, indGlo;     
1822      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
1823                                                          ite = globalLocalIndexMap_.end(), it;         
1824      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
1825                                       itSrve = writtenGlobalIndex.end(), itSrv;
1826
1827//      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1828//      {
1829//        indGlo = *itSrv;
1830//        if (ite != globalLocalIndexMap_.find(indGlo))
1831//        {
1832//          ++nbWritten;
1833//        }
1834//      }
1835
1836//      localIndexToWriteOnServer.resize(nbWritten);
1837      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
1838
1839      nbWritten = 0;
1840      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1841      {
1842        indGlo = *itSrv;
1843        if (ite != globalLocalIndexMap_.find(indGlo))
1844        {
1845          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
1846          ++nbWritten;
1847        }
1848        else
1849        {
1850          localIndexToWriteOnServer(nbWritten) = 0;
1851          ++nbWritten;
1852        }
1853      }
1854     
1855      // if (isCompressible())
1856      // {
1857      //   nbWritten = 0;
1858      //   std::unordered_map<size_t,size_t> localGlobalIndexMap;
1859      //   for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1860      //   {
1861      //     indGlo = *itSrv;
1862      //     if (ite != globalLocalIndexMap_.find(indGlo))
1863      //     {
1864      //       localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
1865      //       ++nbWritten;
1866      //     }                 
1867      //   }
1868
1869      //   nbWritten = 0;
1870      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1871      //   {
1872      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1873      //     {
1874      //       ++nbWritten;
1875      //     }
1876      //   }
1877
1878      //   compressedIndexToWriteOnServer.resize(nbWritten);
1879      //   nbWritten = 0;
1880      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1881      //   {
1882      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1883      //     {
1884      //       compressedIndexToWriteOnServer(nbWritten) = localGlobalIndexMap[data_i_index(idx)];
1885      //       ++nbWritten;
1886      //     }
1887      //   }
1888
1889      //   numberWrittenIndexes_ = nbWritten;
1890      //   if (isDistributed())
1891      //   {           
1892      //     MPI_Allreduce(&numberWrittenIndexes_, &totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1893      //     MPI_Scan(&numberWrittenIndexes_, &offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1894      //     offsetWrittenIndexes_ -= numberWrittenIndexes_;
1895      //   }
1896      //   else
1897      //     totalNumberWrittenIndexes_ = numberWrittenIndexes_;
1898      // }     
1899   }
1900
1901  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
1902  {
1903    int writtenCommSize;
1904    MPI_Comm_size(writtenComm, &writtenCommSize);
1905    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
1906      return;
1907
1908    if (isCompressible())
1909    {
1910      size_t nbWritten = 0, indGlo;
1911      CContext* context=CContext::getCurrent();     
1912      CContextServer* server = context->server; 
1913
1914      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1915      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
1916      nSize[0]        = ni;      nSize[1]  = nj;
1917      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1918      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1919      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
1920      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1921
1922      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
1923                                                          ite = globalLocalIndexMap_.end(), it;   
1924      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
1925                                       itSrve = writtenGlobalIndex.end(), itSrv;
1926      std::unordered_map<size_t,size_t> localGlobalIndexMap;
1927      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1928      {
1929        indGlo = *itSrv;
1930        if (ite != globalLocalIndexMap_.find(indGlo))
1931        {
1932          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
1933          ++nbWritten;
1934        }                 
1935      }
1936
1937      nbWritten = 0;
1938      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1939      {
1940        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1941        {
1942          ++nbWritten;
1943        }
1944      }
1945
1946      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
1947      nbWritten = 0;
1948      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1949      {
1950        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1951        {
1952          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
1953          ++nbWritten;
1954        }
1955      }
1956
1957      numberWrittenIndexes_[writtenCommSize] = nbWritten;
1958      if (isDistributed())
1959      {
1960             
1961        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
1962        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
1963        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
1964      }
1965      else
1966        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
1967      }
1968  }
1969
1970  /*!
1971    Send all attributes from client to connected clients
1972    The attributes will be rebuilt on receiving side
1973  */
1974  void CDomain::sendAttributes()
1975  {
1976    sendDistributionAttributes();
1977    sendIndex();       
1978    sendMask();
1979    sendLonLat();
1980    sendArea();   
1981    sendDataIndex();
1982  }
1983
1984  /*!
1985    Send global index from client to connected client(s)
1986  */
1987  void CDomain::sendIndex()
1988  {
1989    int ns, n, i, j, ind, nv, idx;
1990    std::list<CContextClient*>::iterator it;
1991    for (it=clients.begin(); it!=clients.end(); ++it)
1992    {
1993      CContextClient* client = *it;
1994
1995      int serverSize = client->serverSize;
1996      CEventClient eventIndex(getType(), EVENT_ID_INDEX);
1997
1998      list<CMessage> list_msgsIndex;
1999      list<CArray<int,1> > list_indGlob;
2000
2001      std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2002      iteIndex = indSrv_[serverSize].end();
2003      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2004      {
2005        int nbIndGlob = 0;
2006        int rank = connectedServerRank_[serverSize][k];
2007        itIndex = indSrv_[serverSize].find(rank);
2008        if (iteIndex != itIndex)
2009          nbIndGlob = itIndex->second.size();
2010
2011        list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2012
2013        CArray<int,1>& indGlob = list_indGlob.back();
2014        for (n = 0; n < nbIndGlob; ++n)
2015        {
2016          indGlob(n) = static_cast<int>(itIndex->second[n]);
2017        }
2018
2019        list_msgsIndex.push_back(CMessage());
2020        list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
2021        list_msgsIndex.back() << isCurvilinear;
2022        list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2023       
2024        eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2025      }
2026
2027      client->sendEvent(eventIndex);
2028    }
2029  }
2030
2031  /*!
2032    Send distribution from client to other clients
2033    Because a client in a level knows correctly the grid distribution of client on the next level
2034    it calculates this distribution then sends it to the corresponding clients on the next level
2035  */
2036  void CDomain::sendDistributionAttributes(void)
2037  {
2038    std::list<CContextClient*>::iterator it;
2039    for (it=clients.begin(); it!=clients.end(); ++it)
2040    {
2041      CContextClient* client = *it;
2042      int nbServer = client->serverSize;
2043      std::vector<int> nGlobDomain(2);
2044      nGlobDomain[0] = this->ni_glo;
2045      nGlobDomain[1] = this->nj_glo;
2046
2047      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2048      if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2049      else serverDescription.computeServerDistribution(false, 1);
2050
2051      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2052      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2053
2054      CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2055      if (client->isServerLeader())
2056      {
2057        std::list<CMessage> msgs;
2058
2059        const std::list<int>& ranks = client->getRanksServerLeader();
2060        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2061        {
2062          // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2063          const int ibegin_srv = serverIndexBegin[*itRank][0];
2064          const int jbegin_srv = serverIndexBegin[*itRank][1];
2065          const int ni_srv = serverDimensionSizes[*itRank][0];
2066          const int nj_srv = serverDimensionSizes[*itRank][1];
2067
2068          msgs.push_back(CMessage());
2069          CMessage& msg = msgs.back();
2070          msg << this->getId() ;
2071          msg << isUnstructed_;
2072          msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2073          msg << ni_glo.getValue() << nj_glo.getValue();
2074          msg << isCompressible_;
2075
2076          event.push(*itRank,1,msg);
2077        }
2078        client->sendEvent(event);
2079      }
2080      else client->sendEvent(event);
2081    }
2082  }
2083
2084  /*!
2085    Send mask index from client to connected(s) clients   
2086  */
2087  void CDomain::sendMask()
2088  {
2089    int ns, n, i, j, ind, nv, idx;
2090    std::list<CContextClient*>::iterator it;
2091    for (it=clients.begin(); it!=clients.end(); ++it)
2092    {
2093      CContextClient* client = *it;
2094      int serverSize = client->serverSize;
2095
2096      // send area for each connected server
2097      CEventClient eventMask(getType(), EVENT_ID_MASK);
2098
2099      list<CMessage> list_msgsMask;
2100      list<CArray<bool,1> > list_mask;
2101
2102      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2103      iteMap = indSrv_[serverSize].end();
2104      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2105      {
2106        int nbData = 0;
2107        int rank = connectedServerRank_[serverSize][k];
2108        it = indSrv_[serverSize].find(rank);
2109        if (iteMap != it)
2110          nbData = it->second.size();
2111        list_mask.push_back(CArray<bool,1>(nbData));
2112
2113        const std::vector<size_t>& temp = it->second;
2114        for (n = 0; n < nbData; ++n)
2115        {
2116          idx = static_cast<int>(it->second[n]);
2117          list_mask.back()(n) = domainMask(globalLocalIndexMap_[idx]);
2118        }
2119
2120        list_msgsMask.push_back(CMessage());
2121        list_msgsMask.back() << this->getId() << list_mask.back();
2122        eventMask.push(rank, nbSenders[serverSize][rank], list_msgsMask.back());
2123      }
2124      client->sendEvent(eventMask);
2125    }
2126  }
2127
2128  /*!
2129    Send area from client to connected client(s)
2130  */
2131  void CDomain::sendArea()
2132  {
2133    if (!hasArea) return;
2134
2135    int ns, n, i, j, ind, nv, idx;
2136    std::list<CContextClient*>::iterator it;
2137
2138    for (it=clients.begin(); it!=clients.end(); ++it)
2139    {
2140      CContextClient* client = *it;
2141      int serverSize = client->serverSize;
2142
2143      // send area for each connected server
2144      CEventClient eventArea(getType(), EVENT_ID_AREA);
2145
2146      list<CMessage> list_msgsArea;
2147      list<CArray<double,1> > list_area;
2148
2149      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2150      iteMap = indSrv_[serverSize].end();
2151      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2152      {
2153        int nbData = 0;
2154        int rank = connectedServerRank_[serverSize][k];
2155        it = indSrv_[serverSize].find(rank);
2156        if (iteMap != it)
2157          nbData = it->second.size();
2158        list_area.push_back(CArray<double,1>(nbData));
2159
2160        const std::vector<size_t>& temp = it->second;
2161        for (n = 0; n < nbData; ++n)
2162        {
2163          idx = static_cast<int>(it->second[n]);
2164          list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2165        }
2166
2167        list_msgsArea.push_back(CMessage());
2168        list_msgsArea.back() << this->getId() << hasArea;
2169        list_msgsArea.back() << list_area.back();
2170        eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2171      }
2172      client->sendEvent(eventArea);
2173    }
2174  }
2175
2176  /*!
2177    Send longitude and latitude from client to servers
2178    Each client send long and lat information to corresponding connected clients(s).
2179    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2180  */
2181  void CDomain::sendLonLat()
2182  {
2183    if (!hasLonLat) return;
2184
2185    int ns, n, i, j, ind, nv, idx;
2186    std::list<CContextClient*>::iterator it;
2187    for (it=clients.begin(); it!=clients.end(); ++it)
2188    {
2189      CContextClient* client = *it;
2190      int serverSize = client->serverSize;
2191
2192      // send lon lat for each connected server
2193      CEventClient eventLon(getType(), EVENT_ID_LON);
2194      CEventClient eventLat(getType(), EVENT_ID_LAT);
2195
2196      list<CMessage> list_msgsLon, list_msgsLat;
2197      list<CArray<double,1> > list_lon, list_lat;
2198      list<CArray<double,2> > list_boundslon, list_boundslat;
2199
2200      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2201      iteMap = indSrv_[serverSize].end();
2202      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2203      {
2204        int nbData = 0;
2205        int rank = connectedServerRank_[serverSize][k];
2206        it = indSrv_[serverSize].find(rank);
2207        if (iteMap != it)
2208          nbData = it->second.size();
2209
2210        list_lon.push_back(CArray<double,1>(nbData));
2211        list_lat.push_back(CArray<double,1>(nbData));
2212
2213        if (hasBounds)
2214        {
2215          list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2216          list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2217        }
2218
2219        CArray<double,1>& lon = list_lon.back();
2220        CArray<double,1>& lat = list_lat.back();
2221        const std::vector<size_t>& temp = it->second;
2222        for (n = 0; n < nbData; ++n)
2223        {
2224          idx = static_cast<int>(it->second[n]);
2225          int localInd = globalLocalIndexMap_[idx];
2226          lon(n) = lonvalue(localInd);
2227          lat(n) = latvalue(localInd);
2228
2229          if (hasBounds)
2230          {
2231            CArray<double,2>& boundslon = list_boundslon.back();
2232            CArray<double,2>& boundslat = list_boundslat.back();
2233
2234            for (nv = 0; nv < nvertex; ++nv)
2235            {
2236              boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2237              boundslat(nv, n) = bounds_latvalue(nv, localInd);
2238            }
2239          }
2240        }
2241
2242        list_msgsLon.push_back(CMessage());
2243        list_msgsLat.push_back(CMessage());
2244
2245        list_msgsLon.back() << this->getId() << hasLonLat;
2246        if (hasLonLat) 
2247          list_msgsLon.back() << list_lon.back();
2248        list_msgsLon.back()  << hasBounds;
2249        if (hasBounds)
2250        {
2251          list_msgsLon.back() << list_boundslon.back();
2252        }
2253
2254        list_msgsLat.back() << this->getId() << hasLonLat;
2255        if (hasLonLat)
2256          list_msgsLat.back() << list_lat.back();
2257        list_msgsLat.back() << hasBounds;
2258        if (hasBounds)
2259        {         
2260          list_msgsLat.back() << list_boundslat.back();
2261        }
2262
2263        eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2264        eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2265      }
2266      client->sendEvent(eventLon);
2267      client->sendEvent(eventLat);
2268    }
2269  }
2270
2271  /*!
2272    Send data index to corresponding connected clients.
2273    Data index can be compressed however, we always send decompressed data index
2274    and they will be compressed on receiving.
2275    The compressed index are represented with 1 and others are represented with -1
2276  */
2277  void CDomain::sendDataIndex()
2278  {
2279    int ns, n, i, j, ind, nv, idx;
2280    std::list<CContextClient*>::iterator it;
2281    for (it=clients.begin(); it!=clients.end(); ++it)
2282    {
2283      CContextClient* client = *it;
2284
2285      int serverSize = client->serverSize;
2286
2287      // send area for each connected server
2288      CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2289
2290      list<CMessage> list_msgsDataIndex;
2291      list<CArray<int,1> > list_data_i_index, list_data_j_index;
2292
2293      int nbIndex = i_index.numElements();
2294      int niByIndex = max(i_index) - min(i_index) + 1;
2295      int njByIndex = max(j_index) - min(j_index) + 1; 
2296      int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2297      int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2298
2299     
2300      CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2301      dataIIndex = -1; 
2302      dataJIndex = -1;
2303      ind = 0;
2304
2305      for (idx = 0; idx < data_i_index.numElements(); ++idx)
2306      {
2307        int dataIidx = data_i_index(idx) + data_ibegin;
2308        int dataJidx = data_j_index(idx) + data_jbegin;
2309        if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2310            (0 <= dataJidx) && (dataJidx < dataJindexBound))
2311        {
2312          dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2313          dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2314        }
2315      }
2316
2317      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2318      iteMap = indSrv_[serverSize].end();
2319      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2320      {
2321        int nbData = 0;
2322        int rank = connectedServerRank_[serverSize][k];
2323        it = indSrv_[serverSize].find(rank);
2324        if (iteMap != it)
2325          nbData = it->second.size();
2326        list_data_i_index.push_back(CArray<int,1>(nbData));
2327        list_data_j_index.push_back(CArray<int,1>(nbData));
2328
2329        const std::vector<size_t>& temp = it->second;
2330        for (n = 0; n < nbData; ++n)
2331        {
2332          idx = static_cast<int>(it->second[n]);
2333          i = globalLocalIndexMap_[idx];
2334          list_data_i_index.back()(n) = dataIIndex(i);
2335          list_data_j_index.back()(n) = dataJIndex(i);
2336        }
2337
2338        list_msgsDataIndex.push_back(CMessage());
2339        list_msgsDataIndex.back() << this->getId();
2340        list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2341        eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2342      }
2343      client->sendEvent(eventDataIndex);
2344    }
2345  }
2346 
2347  bool CDomain::dispatchEvent(CEventServer& event)
2348  {
2349    if (SuperClass::dispatchEvent(event)) return true;
2350    else
2351    {
2352      switch(event.type)
2353      {
2354        case EVENT_ID_SERVER_ATTRIBUT:
2355          recvDistributionAttributes(event);
2356          return true;
2357          break;
2358        case EVENT_ID_INDEX:
2359          recvIndex(event);
2360          return true;
2361          break;
2362        case EVENT_ID_MASK:
2363          recvMask(event);
2364          return true;
2365          break;
2366        case EVENT_ID_LON:
2367          recvLon(event);
2368          return true;
2369          break;
2370        case EVENT_ID_LAT:
2371          recvLat(event);
2372          return true;
2373          break;
2374        case EVENT_ID_AREA:
2375          recvArea(event);
2376          return true;
2377          break; 
2378        case EVENT_ID_DATA_INDEX:
2379          recvDataIndex(event);
2380          return true;
2381          break;
2382        default:
2383          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2384                << "Unknown Event");
2385          return false;
2386       }
2387    }
2388  }
2389
2390  /*!
2391    Receive index event from clients(s)
2392    \param[in] event event contain info about rank and associated index
2393  */
2394  void CDomain::recvIndex(CEventServer& event)
2395  {
2396    string domainId;
2397    std::map<int, CBufferIn*> rankBuffers;
2398
2399    list<CEventServer::SSubEvent>::iterator it;
2400    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2401    {     
2402      CBufferIn* buffer = it->buffer;
2403      *buffer >> domainId;
2404      rankBuffers[it->rank] = buffer;       
2405    }
2406    get(domainId)->recvIndex(rankBuffers);
2407  }
2408
2409  /*!
2410    Receive index information from client(s). We use the global index for mapping index between
2411    sending clients and receiving clients.
2412    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2413  */
2414  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2415  {
2416    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2417    recvClientRanks_.resize(nbReceived);       
2418
2419    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2420    ind = 0;
2421    for (ind = 0; it != ite; ++it, ++ind)
2422    {       
2423       recvClientRanks_[ind] = it->first;
2424       CBufferIn& buffer = *(it->second);
2425       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2426       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2427    }
2428    int nbIndGlob = 0;
2429    for (i = 0; i < nbReceived; ++i)
2430    {
2431      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2432    }
2433   
2434    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2435    i_index.resize(nbIndGlob);
2436    j_index.resize(nbIndGlob);   
2437    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2438
2439    nbIndGlob = 0;
2440    for (i = 0; i < nbReceived; ++i)
2441    {
2442      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2443      for (ind = 0; ind < tmp.numElements(); ++ind)
2444      {
2445         index = tmp(ind);
2446         if (0 == globalLocalIndexMap_.count(index))
2447         {
2448           iIndex = (index%ni_glo)-ibegin;
2449           iIndex = (iIndex < 0) ? 0 : iIndex;    // ??
2450           jIndex = (index/ni_glo)-jbegin;
2451           jIndex = (jIndex < 0) ? 0 : jIndex;    // ??
2452           nbIndLoc = iIndex + ni * jIndex;
2453           if (nbIndLoc < nbIndexGlobMax)
2454           {
2455             i_index(nbIndLoc) = index % ni_glo;
2456             j_index(nbIndLoc) = index / ni_glo;
2457             globalLocalIndexMap_[index] = nbIndLoc; 
2458             ++nbIndGlob;
2459           }
2460           // i_index(nbIndGlob) = index % ni_glo;
2461           // j_index(nbIndGlob) = index / ni_glo;
2462           // globalLocalIndexMap_[index] = nbIndGlob; 
2463           // ++nbIndGlob;
2464         } 
2465      } 
2466    } 
2467
2468    if (nbIndGlob==0)
2469    {
2470      i_index.resize(nbIndGlob);
2471      j_index.resize(nbIndGlob);
2472    }
2473    else
2474    {
2475      i_index.resizeAndPreserve(nbIndGlob);
2476      j_index.resizeAndPreserve(nbIndGlob);
2477    }
2478  }
2479
2480  /*!
2481    Receive attributes event from clients(s)
2482    \param[in] event event contain info about rank and associated attributes
2483  */
2484  void CDomain::recvDistributionAttributes(CEventServer& event)
2485  {
2486    CBufferIn* buffer=event.subEvents.begin()->buffer;
2487    string domainId ;
2488    *buffer>>domainId ;
2489    get(domainId)->recvDistributionAttributes(*buffer);
2490  }
2491
2492  /*!
2493    Receive attributes from client(s)
2494    \param[in] rank rank of client source
2495    \param[in] buffer message containing attributes info
2496  */
2497  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2498  {
2499    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2500    int ni_glo_tmp, nj_glo_tmp;
2501    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2502           >> ni_glo_tmp >> nj_glo_tmp
2503           >> isCompressible_;
2504
2505    ni.setValue(ni_tmp);
2506    ibegin.setValue(ibegin_tmp);
2507    nj.setValue(nj_tmp);
2508    jbegin.setValue(jbegin_tmp);
2509    ni_glo.setValue(ni_glo_tmp);
2510    nj_glo.setValue(nj_glo_tmp);
2511
2512  }
2513
2514  /*!
2515    Receive area event from clients(s)
2516    \param[in] event event contain info about rank and associated area
2517  */
2518  void CDomain::recvMask(CEventServer& event)
2519  {
2520    string domainId;
2521    std::map<int, CBufferIn*> rankBuffers;
2522
2523    list<CEventServer::SSubEvent>::iterator it;
2524    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2525    {     
2526      CBufferIn* buffer = it->buffer;
2527      *buffer >> domainId;
2528      rankBuffers[it->rank] = buffer;     
2529    }
2530    get(domainId)->recvMask(rankBuffers);
2531  }
2532
2533
2534  /*!
2535    Receive mask information from client(s)
2536    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2537  */
2538  void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)
2539  {
2540    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2541    if (nbReceived != recvClientRanks_.size())
2542      ERROR("void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)",
2543           << "The number of sending clients is not correct."
2544           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2545
2546    vector<CArray<bool,1> > recvMaskValue(nbReceived);     
2547    for (i = 0; i < recvClientRanks_.size(); ++i)
2548    {
2549      int rank = recvClientRanks_[i];
2550      CBufferIn& buffer = *(rankBuffers[rank]);     
2551      buffer >> recvMaskValue[i];
2552    }
2553
2554    int nbMaskInd = 0;
2555    for (i = 0; i < nbReceived; ++i)
2556    {
2557      nbMaskInd += recvMaskValue[i].numElements();
2558    }
2559 
2560    if (nbMaskInd != globalLocalIndexMap_.size())
2561      info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2562               << "something must be wrong with mask index "<< std::endl;
2563
2564    nbMaskInd = globalLocalIndexMap_.size();
2565    mask_1d.resize(nbMaskInd);
2566    domainMask.resize(nbMaskInd);
2567    mask_1d = false;
2568   
2569    for (i = 0; i < nbReceived; ++i)
2570    {
2571      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2572      CArray<bool,1>& tmp = recvMaskValue[i];
2573      for (ind = 0; ind < tmp.numElements(); ++ind)
2574      {
2575        lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2576        if (!mask_1d(lInd)) // Only rewrite mask_1d if it's not true
2577          mask_1d(lInd) = tmp(ind);
2578      }
2579    }
2580    domainMask=mask_1d ;
2581  }
2582
2583  /*!
2584    Receive longitude event from clients(s)
2585    \param[in] event event contain info about rank and associated longitude
2586  */
2587  void CDomain::recvLon(CEventServer& event)
2588  {
2589    string domainId;
2590    std::map<int, CBufferIn*> rankBuffers;
2591
2592    list<CEventServer::SSubEvent>::iterator it;
2593    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2594    {     
2595      CBufferIn* buffer = it->buffer;
2596      *buffer >> domainId;
2597      rankBuffers[it->rank] = buffer;       
2598    }
2599    get(domainId)->recvLon(rankBuffers);
2600  }
2601
2602  /*!
2603    Receive longitude information from client(s)
2604    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2605  */
2606  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2607  {
2608    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2609    if (nbReceived != recvClientRanks_.size())
2610      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2611           << "The number of sending clients is not correct."
2612           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2613
2614    vector<CArray<double,1> > recvLonValue(nbReceived);
2615    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2616    for (i = 0; i < recvClientRanks_.size(); ++i)
2617    {
2618      int rank = recvClientRanks_[i];
2619      CBufferIn& buffer = *(rankBuffers[rank]);
2620      buffer >> hasLonLat;
2621      if (hasLonLat)
2622        buffer >> recvLonValue[i];
2623      buffer >> hasBounds;
2624      if (hasBounds)
2625        buffer >> recvBoundsLonValue[i];
2626    }
2627
2628    if (hasLonLat)
2629    {
2630      int nbLonInd = 0;
2631      for (i = 0; i < nbReceived; ++i)
2632      {
2633        nbLonInd += recvLonValue[i].numElements();
2634      }
2635   
2636      if (nbLonInd != globalLocalIndexMap_.size())
2637        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2638                 << "something must be wrong with longitude index "<< std::endl;
2639
2640      nbLonInd = globalLocalIndexMap_.size();
2641      lonvalue.resize(nbLonInd);
2642      if (hasBounds)
2643      {
2644        bounds_lonvalue.resize(nvertex,nbLonInd);
2645        bounds_lonvalue = 0.;
2646      }
2647
2648      nbLonInd = 0;
2649      for (i = 0; i < nbReceived; ++i)
2650      {
2651        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2652        CArray<double,1>& tmp = recvLonValue[i];
2653        for (ind = 0; ind < tmp.numElements(); ++ind)
2654        {
2655          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2656          lonvalue(lInd) = tmp(ind); 
2657           if (hasBounds)
2658           {         
2659            for (int nv = 0; nv < nvertex; ++nv)
2660              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2661           }                 
2662        }
2663      }       
2664    }
2665  }
2666
2667  /*!
2668    Receive latitude event from clients(s)
2669    \param[in] event event contain info about rank and associated latitude
2670  */
2671  void CDomain::recvLat(CEventServer& event)
2672  {
2673    string domainId;
2674    std::map<int, CBufferIn*> rankBuffers;
2675
2676    list<CEventServer::SSubEvent>::iterator it;
2677    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2678    {     
2679      CBufferIn* buffer = it->buffer;
2680      *buffer >> domainId;
2681      rankBuffers[it->rank] = buffer;   
2682    }
2683    get(domainId)->recvLat(rankBuffers);
2684  }
2685
2686  /*!
2687    Receive latitude information from client(s)
2688    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2689  */
2690  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2691  {
2692    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2693    if (nbReceived != recvClientRanks_.size())
2694      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2695           << "The number of sending clients is not correct."
2696           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2697
2698    vector<CArray<double,1> > recvLatValue(nbReceived);
2699    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2700    for (i = 0; i < recvClientRanks_.size(); ++i)
2701    {
2702      int rank = recvClientRanks_[i];
2703      CBufferIn& buffer = *(rankBuffers[rank]);
2704      buffer >> hasLonLat;
2705      if (hasLonLat)
2706        buffer >> recvLatValue[i];
2707      buffer >> hasBounds;
2708      if (hasBounds)
2709        buffer >> recvBoundsLatValue[i];
2710    }
2711
2712    if (hasLonLat)
2713    {
2714      int nbLatInd = 0;
2715      for (i = 0; i < nbReceived; ++i)
2716      {
2717        nbLatInd += recvLatValue[i].numElements();
2718      }
2719   
2720      if (nbLatInd != globalLocalIndexMap_.size())
2721        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2722                << "something must be wrong with latitude index "<< std::endl;
2723
2724      nbLatInd = globalLocalIndexMap_.size();
2725      latvalue.resize(nbLatInd);
2726      if (hasBounds)
2727      {
2728        bounds_latvalue.resize(nvertex,nbLatInd);
2729        bounds_latvalue = 0. ;
2730      }
2731
2732      nbLatInd = 0;
2733      for (i = 0; i < nbReceived; ++i)
2734      {
2735        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2736        CArray<double,1>& tmp = recvLatValue[i];
2737        for (ind = 0; ind < tmp.numElements(); ++ind)
2738        {
2739          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2740          latvalue(lInd) = tmp(ind);   
2741           if (hasBounds)
2742           {
2743            CArray<double,2>& boundslat = recvBoundsLatValue[i];
2744            for (int nv = 0; nv < nvertex; ++nv)
2745              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2746           }   
2747          ++nbLatInd;
2748        }
2749      }       
2750    }
2751  }
2752
2753  /*!
2754    Receive area event from clients(s)
2755    \param[in] event event contain info about rank and associated area
2756  */
2757  void CDomain::recvArea(CEventServer& event)
2758  {
2759    string domainId;
2760    std::map<int, CBufferIn*> rankBuffers;
2761
2762    list<CEventServer::SSubEvent>::iterator it;
2763    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2764    {     
2765      CBufferIn* buffer = it->buffer;
2766      *buffer >> domainId;
2767      rankBuffers[it->rank] = buffer;     
2768    }
2769    get(domainId)->recvArea(rankBuffers);
2770  }
2771
2772  /*!
2773    Receive area information from client(s)
2774    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2775  */
2776  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
2777  {
2778    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2779    if (nbReceived != recvClientRanks_.size())
2780      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2781           << "The number of sending clients is not correct."
2782           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2783
2784    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2785    for (i = 0; i < recvClientRanks_.size(); ++i)
2786    {
2787      int rank = recvClientRanks_[i];
2788      CBufferIn& buffer = *(rankBuffers[rank]);     
2789      buffer >> hasArea;
2790      if (hasArea)
2791        buffer >> recvAreaValue[i];
2792    }
2793
2794    if (hasArea)
2795    {
2796      int nbAreaInd = 0;
2797      for (i = 0; i < nbReceived; ++i)
2798      {     
2799        nbAreaInd += recvAreaValue[i].numElements();
2800      }
2801
2802      if (nbAreaInd != globalLocalIndexMap_.size())
2803        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2804                 << "something must be wrong with area index "<< std::endl;
2805
2806      nbAreaInd = globalLocalIndexMap_.size();
2807      areavalue.resize(nbAreaInd);
2808      nbAreaInd = 0;     
2809      for (i = 0; i < nbReceived; ++i)
2810      {
2811        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2812        CArray<double,1>& tmp = recvAreaValue[i];
2813        for (ind = 0; ind < tmp.numElements(); ++ind)
2814        {
2815          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2816          areavalue(lInd) = tmp(ind);         
2817        }
2818      }
2819     
2820    }
2821  }
2822
2823  /*!
2824    Compare two domain objects.
2825    They are equal if only if they have identical attributes as well as their values.
2826    Moreover, they must have the same transformations.
2827  \param [in] domain Compared domain
2828  \return result of the comparison
2829  */
2830  bool CDomain::isEqual(CDomain* obj)
2831  {
2832    vector<StdString> excludedAttr;
2833    excludedAttr.push_back("domain_ref");
2834    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2835    if (!objEqual) return objEqual;
2836
2837    TransMapTypes thisTrans = this->getAllTransformations();
2838    TransMapTypes objTrans  = obj->getAllTransformations();
2839
2840    TransMapTypes::const_iterator it, itb, ite;
2841    std::vector<ETranformationType> thisTransType, objTransType;
2842    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2843      thisTransType.push_back(it->first);
2844    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2845      objTransType.push_back(it->first);
2846
2847    if (thisTransType.size() != objTransType.size()) return false;
2848    for (int idx = 0; idx < thisTransType.size(); ++idx)
2849      objEqual &= (thisTransType[idx] == objTransType[idx]);
2850
2851    return objEqual;
2852  }
2853
2854  /*!
2855    Receive data index event from clients(s)
2856    \param[in] event event contain info about rank and associated index
2857  */
2858  void CDomain::recvDataIndex(CEventServer& event)
2859  {
2860    string domainId;
2861    std::map<int, CBufferIn*> rankBuffers;
2862
2863    list<CEventServer::SSubEvent>::iterator it;
2864    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2865    {     
2866      CBufferIn* buffer = it->buffer;
2867      *buffer >> domainId;
2868      rankBuffers[it->rank] = buffer;       
2869    }
2870    get(domainId)->recvDataIndex(rankBuffers);
2871  }
2872
2873  /*!
2874    Receive data index information from client(s)
2875    A client receives data index from different clients to rebuild its own data index.
2876    Because we use global index + mask info to calculate the sending data to client(s),
2877    this data index must be updated with mask info (maybe it will change in the future)
2878    Because the data index is local, to rebuild data index of received client, we should use global index along with.
2879
2880    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2881  */
2882  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
2883  {
2884    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
2885    if (nbReceived != recvClientRanks_.size())
2886      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
2887           << "The number of sending clients is not correct."
2888           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2889
2890    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
2891    for (i = 0; i < recvClientRanks_.size(); ++i)
2892    {
2893      int rank = recvClientRanks_[i];
2894      CBufferIn& buffer = *(rankBuffers[rank]);
2895      buffer >> recvDataIIndex[i];
2896      buffer >> recvDataJIndex[i];
2897    }
2898
2899    int nbIndex = i_index.numElements();
2900    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2901    dataIIndex = -1; dataJIndex = -1;
2902
2903    nbIndex = 0;
2904    for (i = 0; i < nbReceived; ++i)
2905    {     
2906      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2907      CArray<int,1>& tmpI = recvDataIIndex[i];   
2908      CArray<int,1>& tmpJ = recvDataJIndex[i];     
2909      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
2910          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
2911             << "The number of global received index is not coherent with the number of received data index."
2912             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
2913
2914      for (ind = 0; ind < tmpI.numElements(); ++ind)
2915      {
2916         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2917         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
2918         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
2919
2920         if (!domainMask(lInd))   // Include mask info into data index on the RECEIVE getServerDimensionSizes
2921         {
2922           dataIIndex(lInd) = dataJIndex(lInd) = -1;
2923         }
2924      } 
2925    }
2926
2927    int nbCompressedData = 0; 
2928    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
2929    {
2930       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
2931       if ((0 <= indexI) && (0 <= indexJ))
2932         ++nbCompressedData;
2933    }       
2934 
2935    data_i_index.resize(nbCompressedData);
2936    data_j_index.resize(nbCompressedData);
2937
2938    nbCompressedData = 0; 
2939    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
2940    {
2941       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
2942       if ((0 <= indexI) && (0 <= indexJ))
2943       {
2944          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
2945          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
2946         ++nbCompressedData;
2947       }
2948    }
2949
2950    // Reset data_ibegin, data_jbegin
2951    data_ibegin.setValue(0);
2952    data_jbegin.setValue(0);
2953  }
2954
2955  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2956  {
2957    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2958    return transformationMap_.back().second;
2959  }
2960
2961  /*!
2962    Check whether a domain has transformation
2963    \return true if domain has transformation
2964  */
2965  bool CDomain::hasTransformation()
2966  {
2967    return (!transformationMap_.empty());
2968  }
2969
2970  /*!
2971    Set transformation for current domain. It's the method to move transformation in hierarchy
2972    \param [in] domTrans transformation on domain
2973  */
2974  void CDomain::setTransformations(const TransMapTypes& domTrans)
2975  {
2976    transformationMap_ = domTrans;
2977  }
2978
2979  /*!
2980    Get all transformation current domain has
2981    \return all transformation
2982  */
2983  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2984  {
2985    return transformationMap_;
2986  }
2987
2988  void CDomain::duplicateTransformation(CDomain* src)
2989  {
2990    if (src->hasTransformation())
2991    {
2992      this->setTransformations(src->getAllTransformations());
2993    }
2994  }
2995
2996  /*!
2997   * Go through the hierarchy to find the domain from which the transformations must be inherited
2998   */
2999  void CDomain::solveInheritanceTransformation()
3000  {
3001    if (hasTransformation() || !hasDirectDomainReference())
3002      return;
3003
3004    CDomain* domain = this;
3005    std::vector<CDomain*> refDomains;
3006    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3007    {
3008      refDomains.push_back(domain);
3009      domain = domain->getDirectDomainReference();
3010    }
3011
3012    if (domain->hasTransformation())
3013      for (size_t i = 0; i < refDomains.size(); ++i)
3014        refDomains[i]->setTransformations(domain->getAllTransformations());
3015  }
3016
3017  void CDomain::setContextClient(CContextClient* contextClient)
3018  {
3019    if (clientsSet.find(contextClient)==clientsSet.end())
3020    {
3021      clients.push_back(contextClient) ;
3022      clientsSet.insert(contextClient);
3023    }
3024  }
3025
3026  /*!
3027    Parse children nodes of a domain in xml file.
3028    Whenver there is a new transformation, its type and name should be added into this function
3029    \param node child node to process
3030  */
3031  void CDomain::parse(xml::CXMLNode & node)
3032  {
3033    SuperClass::parse(node);
3034
3035    if (node.goToChildElement())
3036    {
3037      StdString nodeElementName;
3038      do
3039      {
3040        StdString nodeId("");
3041        if (node.getAttributes().end() != node.getAttributes().find("id"))
3042        { nodeId = node.getAttributes()["id"]; }
3043
3044        nodeElementName = node.getElementName();
3045        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3046        it = transformationMapList_.find(nodeElementName);
3047        if (ite != it)
3048        {
3049          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3050                                                                                                                nodeId,
3051                                                                                                                &node)));
3052        }
3053        else
3054        {
3055          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3056                << "The transformation " << nodeElementName << " has not been supported yet.");
3057        }
3058      } while (node.goToNextElement()) ;
3059      node.goToParentElement();
3060    }
3061  }
3062   //----------------------------------------------------------------
3063
3064   DEFINE_REF_FUNC(Domain,domain)
3065
3066   ///---------------------------------------------------------------
3067
3068} // namespace xios
Note: See TracBrowser for help on using the repository browser.