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

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

Taking care of cases when a process doesn't possess a grid (its domain or axe size is zero).

Tested with CMIP6 toy models and IPSL model.

To do: remove grid mask from transformations.

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