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

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

Further simplifications on sending data/grid indexes.

(1) Domain/axis mask is not sent anymore. It has been incorporated into data index.

(2) Creating a map that holds grid mask only in case if grid mask is defined.

Still to fix: a hole in a domain or axis.

  • 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: 110.9 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          nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
1854          delete clientServerMap;
1855        }
1856      }
1857    }
1858  }
1859
1860   /*!
1861     Compute index to write data. We only write data on the zoomed region, therefore, there should
1862     be a map between the complete grid and the reduced grid where we write data.
1863     By using global index we can easily create this kind of mapping.
1864   */
1865   void CDomain::computeWrittenIndex()
1866   { 
1867      if (computedWrittenIndex_) return;
1868      computedWrittenIndex_ = true;
1869
1870      CContext* context=CContext::getCurrent();     
1871      CContextServer* server = context->server; 
1872
1873      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1874      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
1875      nSize[0]        = ni;      nSize[1]  = nj;
1876      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1877      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1878      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
1879      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1880
1881      size_t nbWritten = 0, indGlo;     
1882      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
1883                                                          ite = globalLocalIndexMap_.end(), it;         
1884      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
1885                                       itSrve = writtenGlobalIndex.end(), itSrv;
1886
1887//      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1888//      {
1889//        indGlo = *itSrv;
1890//        if (ite != globalLocalIndexMap_.find(indGlo))
1891//        {
1892//          ++nbWritten;
1893//        }
1894//      }
1895
1896//      localIndexToWriteOnServer.resize(nbWritten);
1897      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
1898
1899      nbWritten = 0;
1900      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1901      {
1902        indGlo = *itSrv;
1903        if (ite != globalLocalIndexMap_.find(indGlo))
1904        {
1905          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
1906          ++nbWritten;
1907        }
1908        else
1909        {
1910          localIndexToWriteOnServer(nbWritten) = 0;
1911          ++nbWritten;
1912        }
1913      }
1914     
1915      // if (isCompressible())
1916      // {
1917      //   nbWritten = 0;
1918      //   std::unordered_map<size_t,size_t> localGlobalIndexMap;
1919      //   for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1920      //   {
1921      //     indGlo = *itSrv;
1922      //     if (ite != globalLocalIndexMap_.find(indGlo))
1923      //     {
1924      //       localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
1925      //       ++nbWritten;
1926      //     }                 
1927      //   }
1928
1929      //   nbWritten = 0;
1930      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1931      //   {
1932      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1933      //     {
1934      //       ++nbWritten;
1935      //     }
1936      //   }
1937
1938      //   compressedIndexToWriteOnServer.resize(nbWritten);
1939      //   nbWritten = 0;
1940      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1941      //   {
1942      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1943      //     {
1944      //       compressedIndexToWriteOnServer(nbWritten) = localGlobalIndexMap[data_i_index(idx)];
1945      //       ++nbWritten;
1946      //     }
1947      //   }
1948
1949      //   numberWrittenIndexes_ = nbWritten;
1950      //   if (isDistributed())
1951      //   {           
1952      //     MPI_Allreduce(&numberWrittenIndexes_, &totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1953      //     MPI_Scan(&numberWrittenIndexes_, &offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1954      //     offsetWrittenIndexes_ -= numberWrittenIndexes_;
1955      //   }
1956      //   else
1957      //     totalNumberWrittenIndexes_ = numberWrittenIndexes_;
1958      // }     
1959   }
1960
1961  void CDomain::computeWrittenCompressedIndex(MPI_Comm writtenComm)
1962  {
1963    int writtenCommSize;
1964    MPI_Comm_size(writtenComm, &writtenCommSize);
1965    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
1966      return;
1967
1968    if (isCompressible())
1969    {
1970      size_t nbWritten = 0, indGlo;
1971      CContext* context=CContext::getCurrent();     
1972      CContextServer* server = context->server; 
1973
1974      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1975      nBegin[0]       = ibegin;  nBegin[1] = jbegin;
1976      nSize[0]        = ni;      nSize[1]  = nj;
1977      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1978      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1979      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
1980      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1981
1982      std::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
1983                                                          ite = globalLocalIndexMap_.end(), it;   
1984      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
1985                                       itSrve = writtenGlobalIndex.end(), itSrv;
1986      std::unordered_map<size_t,size_t> localGlobalIndexMap;
1987      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1988      {
1989        indGlo = *itSrv;
1990        if (ite != globalLocalIndexMap_.find(indGlo))
1991        {
1992          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
1993          ++nbWritten;
1994        }                 
1995      }
1996
1997      nbWritten = 0;
1998      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1999      {
2000        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2001        {
2002          ++nbWritten;
2003        }
2004      }
2005
2006      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2007      nbWritten = 0;
2008      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2009      {
2010        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2011        {
2012          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2013          ++nbWritten;
2014        }
2015      }
2016
2017      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2018      if (isDistributed())
2019      {
2020             
2021        MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2022        MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2023        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2024      }
2025      else
2026        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2027      }
2028  }
2029
2030  /*!
2031    Send all attributes from client to connected clients
2032    The attributes will be rebuilt on receiving side
2033  */
2034  void CDomain::sendAttributes()
2035  {
2036    sendDistributionAttributes();
2037    sendIndex();       
2038//    sendMask();
2039    sendLonLat();
2040    sendArea();   
2041    sendDataIndex();
2042  }
2043
2044  /*!
2045    Send global index from client to connected client(s)
2046  */
2047  void CDomain::sendIndex()
2048  {
2049    int ns, n, i, j, ind, nv, idx;
2050    std::list<CContextClient*>::iterator it;
2051    for (it=clients.begin(); it!=clients.end(); ++it)
2052    {
2053      CContextClient* client = *it;
2054
2055      int serverSize = client->serverSize;
2056      CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2057
2058      list<CMessage> list_msgsIndex;
2059      list<CArray<int,1> > list_indGlob;
2060
2061      std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2062      iteIndex = indSrv_[serverSize].end();
2063      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2064      {
2065        int nbIndGlob = 0;
2066        int rank = connectedServerRank_[serverSize][k];
2067        itIndex = indSrv_[serverSize].find(rank);
2068        if (iteIndex != itIndex)
2069          nbIndGlob = itIndex->second.size();
2070
2071        list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2072
2073        CArray<int,1>& indGlob = list_indGlob.back();
2074        for (n = 0; n < nbIndGlob; ++n)
2075        {
2076          indGlob(n) = static_cast<int>(itIndex->second[n]);
2077        }
2078
2079        list_msgsIndex.push_back(CMessage());
2080        list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
2081        list_msgsIndex.back() << isCurvilinear;
2082        list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2083       
2084        eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2085      }
2086
2087      client->sendEvent(eventIndex);
2088    }
2089  }
2090
2091  /*!
2092    Send distribution from client to other clients
2093    Because a client in a level knows correctly the grid distribution of client on the next level
2094    it calculates this distribution then sends it to the corresponding clients on the next level
2095  */
2096  void CDomain::sendDistributionAttributes(void)
2097  {
2098    std::list<CContextClient*>::iterator it;
2099    for (it=clients.begin(); it!=clients.end(); ++it)
2100    {
2101      CContextClient* client = *it;
2102      int nbServer = client->serverSize;
2103      std::vector<int> nGlobDomain(2);
2104      nGlobDomain[0] = this->ni_glo;
2105      nGlobDomain[1] = this->nj_glo;
2106
2107      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2108      if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2109      else serverDescription.computeServerDistribution(false, 1);
2110
2111      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2112      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2113
2114      CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2115      if (client->isServerLeader())
2116      {
2117        std::list<CMessage> msgs;
2118
2119        const std::list<int>& ranks = client->getRanksServerLeader();
2120        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2121        {
2122          // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2123          const int ibegin_srv = serverIndexBegin[*itRank][0];
2124          const int jbegin_srv = serverIndexBegin[*itRank][1];
2125          const int ni_srv = serverDimensionSizes[*itRank][0];
2126          const int nj_srv = serverDimensionSizes[*itRank][1];
2127
2128          msgs.push_back(CMessage());
2129          CMessage& msg = msgs.back();
2130          msg << this->getId() ;
2131          msg << isUnstructed_;
2132          msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2133          msg << ni_glo.getValue() << nj_glo.getValue();
2134          msg << isCompressible_;
2135
2136          event.push(*itRank,1,msg);
2137        }
2138        client->sendEvent(event);
2139      }
2140      else client->sendEvent(event);
2141    }
2142  }
2143
2144  /*!
2145    Send mask index from client to connected(s) clients   
2146  */
2147  void CDomain::sendMask()
2148  {
2149    int ns, n, i, j, ind, nv, idx;
2150    std::list<CContextClient*>::iterator it;
2151    for (it=clients.begin(); it!=clients.end(); ++it)
2152    {
2153      CContextClient* client = *it;
2154      int serverSize = client->serverSize;
2155
2156      // send area for each connected server
2157      CEventClient eventMask(getType(), EVENT_ID_MASK);
2158
2159      list<CMessage> list_msgsMask;
2160      list<CArray<bool,1> > list_mask;
2161
2162      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2163      iteMap = indSrv_[serverSize].end();
2164      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2165      {
2166        int nbData = 0;
2167        int rank = connectedServerRank_[serverSize][k];
2168        it = indSrv_[serverSize].find(rank);
2169        if (iteMap != it)
2170          nbData = it->second.size();
2171        list_mask.push_back(CArray<bool,1>(nbData));
2172
2173        const std::vector<size_t>& temp = it->second;
2174        for (n = 0; n < nbData; ++n)
2175        {
2176          idx = static_cast<int>(it->second[n]);
2177          list_mask.back()(n) = domainMask(globalLocalIndexMap_[idx]);
2178        }
2179
2180        list_msgsMask.push_back(CMessage());
2181        list_msgsMask.back() << this->getId() << list_mask.back();
2182        eventMask.push(rank, nbSenders[serverSize][rank], list_msgsMask.back());
2183      }
2184      client->sendEvent(eventMask);
2185    }
2186  }
2187
2188  /*!
2189    Send area from client to connected client(s)
2190  */
2191  void CDomain::sendArea()
2192  {
2193    if (!hasArea) return;
2194
2195    int ns, n, i, j, ind, nv, idx;
2196    std::list<CContextClient*>::iterator it;
2197
2198    for (it=clients.begin(); it!=clients.end(); ++it)
2199    {
2200      CContextClient* client = *it;
2201      int serverSize = client->serverSize;
2202
2203      // send area for each connected server
2204      CEventClient eventArea(getType(), EVENT_ID_AREA);
2205
2206      list<CMessage> list_msgsArea;
2207      list<CArray<double,1> > list_area;
2208
2209      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2210      iteMap = indSrv_[serverSize].end();
2211      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2212      {
2213        int nbData = 0;
2214        int rank = connectedServerRank_[serverSize][k];
2215        it = indSrv_[serverSize].find(rank);
2216        if (iteMap != it)
2217          nbData = it->second.size();
2218        list_area.push_back(CArray<double,1>(nbData));
2219
2220        const std::vector<size_t>& temp = it->second;
2221        for (n = 0; n < nbData; ++n)
2222        {
2223          idx = static_cast<int>(it->second[n]);
2224          list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2225        }
2226
2227        list_msgsArea.push_back(CMessage());
2228        list_msgsArea.back() << this->getId() << hasArea;
2229        list_msgsArea.back() << list_area.back();
2230        eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2231      }
2232      client->sendEvent(eventArea);
2233    }
2234  }
2235
2236  /*!
2237    Send longitude and latitude from client to servers
2238    Each client send long and lat information to corresponding connected clients(s).
2239    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2240  */
2241  void CDomain::sendLonLat()
2242  {
2243    if (!hasLonLat) return;
2244
2245    int ns, n, i, j, ind, nv, idx;
2246    std::list<CContextClient*>::iterator it;
2247    for (it=clients.begin(); it!=clients.end(); ++it)
2248    {
2249      CContextClient* client = *it;
2250      int serverSize = client->serverSize;
2251
2252      // send lon lat for each connected server
2253      CEventClient eventLon(getType(), EVENT_ID_LON);
2254      CEventClient eventLat(getType(), EVENT_ID_LAT);
2255
2256      list<CMessage> list_msgsLon, list_msgsLat;
2257      list<CArray<double,1> > list_lon, list_lat;
2258      list<CArray<double,2> > list_boundslon, list_boundslat;
2259
2260      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2261      iteMap = indSrv_[serverSize].end();
2262      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2263      {
2264        int nbData = 0;
2265        int rank = connectedServerRank_[serverSize][k];
2266        it = indSrv_[serverSize].find(rank);
2267        if (iteMap != it)
2268          nbData = it->second.size();
2269
2270        list_lon.push_back(CArray<double,1>(nbData));
2271        list_lat.push_back(CArray<double,1>(nbData));
2272
2273        if (hasBounds)
2274        {
2275          list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2276          list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2277        }
2278
2279        CArray<double,1>& lon = list_lon.back();
2280        CArray<double,1>& lat = list_lat.back();
2281        const std::vector<size_t>& temp = it->second;
2282        for (n = 0; n < nbData; ++n)
2283        {
2284          idx = static_cast<int>(it->second[n]);
2285          int localInd = globalLocalIndexMap_[idx];
2286          lon(n) = lonvalue(localInd);
2287          lat(n) = latvalue(localInd);
2288
2289          if (hasBounds)
2290          {
2291            CArray<double,2>& boundslon = list_boundslon.back();
2292            CArray<double,2>& boundslat = list_boundslat.back();
2293
2294            for (nv = 0; nv < nvertex; ++nv)
2295            {
2296              boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2297              boundslat(nv, n) = bounds_latvalue(nv, localInd);
2298            }
2299          }
2300        }
2301
2302        list_msgsLon.push_back(CMessage());
2303        list_msgsLat.push_back(CMessage());
2304
2305        list_msgsLon.back() << this->getId() << hasLonLat;
2306        if (hasLonLat) 
2307          list_msgsLon.back() << list_lon.back();
2308        list_msgsLon.back()  << hasBounds;
2309        if (hasBounds)
2310        {
2311          list_msgsLon.back() << list_boundslon.back();
2312        }
2313
2314        list_msgsLat.back() << this->getId() << hasLonLat;
2315        if (hasLonLat)
2316          list_msgsLat.back() << list_lat.back();
2317        list_msgsLat.back() << hasBounds;
2318        if (hasBounds)
2319        {         
2320          list_msgsLat.back() << list_boundslat.back();
2321        }
2322
2323        eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2324        eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2325      }
2326      client->sendEvent(eventLon);
2327      client->sendEvent(eventLat);
2328    }
2329  }
2330
2331  /*!
2332    Send data index to corresponding connected clients.
2333    Data index can be compressed however, we always send decompressed data index
2334    and they will be compressed on receiving.
2335    The compressed index are represented with 1 and others are represented with -1
2336  */
2337  void CDomain::sendDataIndex()
2338  {
2339    int ns, n, i, j, ind, nv, idx;
2340    std::list<CContextClient*>::iterator it;
2341    for (it=clients.begin(); it!=clients.end(); ++it)
2342    {
2343      CContextClient* client = *it;
2344
2345      int serverSize = client->serverSize;
2346
2347      // send area for each connected server
2348      CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2349
2350      list<CMessage> list_msgsDataIndex;
2351      list<CArray<int,1> > list_data_i_index, list_data_j_index;
2352
2353      int nbIndex = i_index.numElements();
2354      int niByIndex = max(i_index) - min(i_index) + 1;
2355      int njByIndex = max(j_index) - min(j_index) + 1; 
2356      int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2357      int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2358
2359     
2360      CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2361      dataIIndex = -1; 
2362      dataJIndex = -1;
2363      ind = 0;
2364
2365      for (idx = 0; idx < data_i_index.numElements(); ++idx)
2366      {
2367        int dataIidx = data_i_index(idx) + data_ibegin;
2368        int dataJidx = data_j_index(idx) + data_jbegin;
2369        if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2370            (0 <= dataJidx) && (dataJidx < dataJindexBound))
2371        {
2372          dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2373          dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2374        }
2375      }
2376
2377      std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2378      iteMap = indSrv_[serverSize].end();
2379      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2380      {
2381        int nbData = 0;
2382        int rank = connectedServerRank_[serverSize][k];
2383        it = indSrv_[serverSize].find(rank);
2384        if (iteMap != it)
2385          nbData = it->second.size();
2386        list_data_i_index.push_back(CArray<int,1>(nbData));
2387        list_data_j_index.push_back(CArray<int,1>(nbData));
2388
2389        const std::vector<size_t>& temp = it->second;
2390        for (n = 0; n < nbData; ++n)
2391        {
2392          idx = static_cast<int>(it->second[n]);
2393          i = globalLocalIndexMap_[idx];
2394          list_data_i_index.back()(n) = dataIIndex(i);
2395          list_data_j_index.back()(n) = dataJIndex(i);
2396        }
2397
2398        list_msgsDataIndex.push_back(CMessage());
2399        list_msgsDataIndex.back() << this->getId();
2400        list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2401        eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2402      }
2403      client->sendEvent(eventDataIndex);
2404    }
2405  }
2406 
2407  bool CDomain::dispatchEvent(CEventServer& event)
2408  {
2409    if (SuperClass::dispatchEvent(event)) return true;
2410    else
2411    {
2412      switch(event.type)
2413      {
2414        case EVENT_ID_SERVER_ATTRIBUT:
2415          recvDistributionAttributes(event);
2416          return true;
2417          break;
2418        case EVENT_ID_INDEX:
2419          recvIndex(event);
2420          return true;
2421          break;
2422        case EVENT_ID_MASK:
2423          recvMask(event);
2424          return true;
2425          break;
2426        case EVENT_ID_LON:
2427          recvLon(event);
2428          return true;
2429          break;
2430        case EVENT_ID_LAT:
2431          recvLat(event);
2432          return true;
2433          break;
2434        case EVENT_ID_AREA:
2435          recvArea(event);
2436          return true;
2437          break; 
2438        case EVENT_ID_DATA_INDEX:
2439          recvDataIndex(event);
2440          return true;
2441          break;
2442        default:
2443          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2444                << "Unknown Event");
2445          return false;
2446       }
2447    }
2448  }
2449
2450  /*!
2451    Receive index event from clients(s)
2452    \param[in] event event contain info about rank and associated index
2453  */
2454  void CDomain::recvIndex(CEventServer& event)
2455  {
2456    string domainId;
2457    std::map<int, CBufferIn*> rankBuffers;
2458
2459    list<CEventServer::SSubEvent>::iterator it;
2460    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2461    {     
2462      CBufferIn* buffer = it->buffer;
2463      *buffer >> domainId;
2464      rankBuffers[it->rank] = buffer;       
2465    }
2466    get(domainId)->recvIndex(rankBuffers);
2467  }
2468
2469  /*!
2470    Receive index information from client(s). We use the global index for mapping index between
2471    sending clients and receiving clients.
2472    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2473  */
2474  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2475  {
2476    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2477    recvClientRanks_.resize(nbReceived);       
2478
2479    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2480    ind = 0;
2481    for (ind = 0; it != ite; ++it, ++ind)
2482    {       
2483       recvClientRanks_[ind] = it->first;
2484       CBufferIn& buffer = *(it->second);
2485       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2486       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2487    }
2488    int nbIndGlob = 0;
2489    for (i = 0; i < nbReceived; ++i)
2490    {
2491      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2492    }
2493   
2494    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2495    i_index.resize(nbIndGlob);
2496    j_index.resize(nbIndGlob);   
2497    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2498
2499    nbIndGlob = 0;
2500    for (i = 0; i < nbReceived; ++i)
2501    {
2502      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2503      for (ind = 0; ind < tmp.numElements(); ++ind)
2504      {
2505         index = tmp(ind);
2506         if (0 == globalLocalIndexMap_.count(index))
2507         {
2508           iIndex = (index%ni_glo)-ibegin;
2509           iIndex = (iIndex < 0) ? 0 : iIndex;
2510           jIndex = (index/ni_glo)-jbegin;
2511           jIndex = (jIndex < 0) ? 0 : jIndex;
2512           nbIndLoc = iIndex + ni * jIndex;
2513//           if (nbIndLoc < nbIndexGlobMax)  // THIS CONDITION IMPEDES THE CASE OF A HOLE
2514           {
2515             i_index(nbIndLoc) = index % ni_glo;
2516             j_index(nbIndLoc) = index / ni_glo;
2517             globalLocalIndexMap_[index] = nbIndLoc; 
2518             ++nbIndGlob;
2519           }
2520           // i_index(nbIndGlob) = index % ni_glo;
2521           // j_index(nbIndGlob) = index / ni_glo;
2522           // globalLocalIndexMap_[index] = nbIndGlob; 
2523           // ++nbIndGlob;
2524         } 
2525      } 
2526    } 
2527
2528    if (nbIndGlob==0)
2529    {
2530      i_index.resize(nbIndGlob);
2531      j_index.resize(nbIndGlob);
2532    }
2533    else
2534    {
2535      i_index.resizeAndPreserve(nbIndGlob);
2536      j_index.resizeAndPreserve(nbIndGlob);
2537    }
2538
2539    domainMask.resize(0); // Mask is not defined anymore on servers
2540  }
2541
2542  /*!
2543    Receive attributes event from clients(s)
2544    \param[in] event event contain info about rank and associated attributes
2545  */
2546  void CDomain::recvDistributionAttributes(CEventServer& event)
2547  {
2548    CBufferIn* buffer=event.subEvents.begin()->buffer;
2549    string domainId ;
2550    *buffer>>domainId ;
2551    get(domainId)->recvDistributionAttributes(*buffer);
2552  }
2553
2554  /*!
2555    Receive attributes from client(s)
2556    \param[in] rank rank of client source
2557    \param[in] buffer message containing attributes info
2558  */
2559  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
2560  {
2561    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
2562    int ni_glo_tmp, nj_glo_tmp;
2563    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2564           >> ni_glo_tmp >> nj_glo_tmp
2565           >> isCompressible_;
2566
2567    ni.setValue(ni_tmp);
2568    ibegin.setValue(ibegin_tmp);
2569    nj.setValue(nj_tmp);
2570    jbegin.setValue(jbegin_tmp);
2571    ni_glo.setValue(ni_glo_tmp);
2572    nj_glo.setValue(nj_glo_tmp);
2573
2574  }
2575
2576  /*!
2577    Receive area event from clients(s)
2578    \param[in] event event contain info about rank and associated area
2579  */
2580  void CDomain::recvMask(CEventServer& event)
2581  {
2582    string domainId;
2583    std::map<int, CBufferIn*> rankBuffers;
2584
2585    list<CEventServer::SSubEvent>::iterator it;
2586    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2587    {     
2588      CBufferIn* buffer = it->buffer;
2589      *buffer >> domainId;
2590      rankBuffers[it->rank] = buffer;     
2591    }
2592    get(domainId)->recvMask(rankBuffers);
2593  }
2594
2595
2596  /*!
2597    Receive mask information from client(s)
2598    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2599  */
2600  void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)
2601  {
2602    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2603    if (nbReceived != recvClientRanks_.size())
2604      ERROR("void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)",
2605           << "The number of sending clients is not correct."
2606           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2607
2608    vector<CArray<bool,1> > recvMaskValue(nbReceived);     
2609    for (i = 0; i < recvClientRanks_.size(); ++i)
2610    {
2611      int rank = recvClientRanks_[i];
2612      CBufferIn& buffer = *(rankBuffers[rank]);     
2613      buffer >> recvMaskValue[i];
2614    }
2615
2616    int nbMaskInd = 0;
2617    for (i = 0; i < nbReceived; ++i)
2618    {
2619      nbMaskInd += recvMaskValue[i].numElements();
2620    }
2621 
2622    if (nbMaskInd != globalLocalIndexMap_.size())
2623      info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2624               << "something must be wrong with mask index "<< std::endl;
2625
2626    nbMaskInd = globalLocalIndexMap_.size();
2627    mask_1d.resize(nbMaskInd);
2628    domainMask.resize(nbMaskInd);
2629    mask_1d = false;
2630   
2631    for (i = 0; i < nbReceived; ++i)
2632    {
2633      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2634      CArray<bool,1>& tmp = recvMaskValue[i];
2635      for (ind = 0; ind < tmp.numElements(); ++ind)
2636      {
2637        lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2638        if (!mask_1d(lInd)) // Only rewrite mask_1d if it's not true
2639          mask_1d(lInd) = tmp(ind);
2640      }
2641    }
2642    domainMask=mask_1d ;
2643  }
2644
2645  /*!
2646    Receive longitude event from clients(s)
2647    \param[in] event event contain info about rank and associated longitude
2648  */
2649  void CDomain::recvLon(CEventServer& event)
2650  {
2651    string domainId;
2652    std::map<int, CBufferIn*> rankBuffers;
2653
2654    list<CEventServer::SSubEvent>::iterator it;
2655    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2656    {     
2657      CBufferIn* buffer = it->buffer;
2658      *buffer >> domainId;
2659      rankBuffers[it->rank] = buffer;       
2660    }
2661    get(domainId)->recvLon(rankBuffers);
2662  }
2663
2664  /*!
2665    Receive longitude information from client(s)
2666    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2667  */
2668  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
2669  {
2670    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2671    if (nbReceived != recvClientRanks_.size())
2672      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2673           << "The number of sending clients is not correct."
2674           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2675
2676    vector<CArray<double,1> > recvLonValue(nbReceived);
2677    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2678    for (i = 0; i < recvClientRanks_.size(); ++i)
2679    {
2680      int rank = recvClientRanks_[i];
2681      CBufferIn& buffer = *(rankBuffers[rank]);
2682      buffer >> hasLonLat;
2683      if (hasLonLat)
2684        buffer >> recvLonValue[i];
2685      buffer >> hasBounds;
2686      if (hasBounds)
2687        buffer >> recvBoundsLonValue[i];
2688    }
2689
2690    if (hasLonLat)
2691    {
2692      int nbLonInd = 0;
2693      for (i = 0; i < nbReceived; ++i)
2694      {
2695        nbLonInd += recvLonValue[i].numElements();
2696      }
2697   
2698      if (nbLonInd != globalLocalIndexMap_.size())
2699        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2700                 << "something must be wrong with longitude index "<< std::endl;
2701
2702      nbLonInd = globalLocalIndexMap_.size();
2703      lonvalue.resize(nbLonInd);
2704      if (hasBounds)
2705      {
2706        bounds_lonvalue.resize(nvertex,nbLonInd);
2707        bounds_lonvalue = 0.;
2708      }
2709
2710      nbLonInd = 0;
2711      for (i = 0; i < nbReceived; ++i)
2712      {
2713        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2714        CArray<double,1>& tmp = recvLonValue[i];
2715        for (ind = 0; ind < tmp.numElements(); ++ind)
2716        {
2717          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2718          lonvalue(lInd) = tmp(ind); 
2719           if (hasBounds)
2720           {         
2721            for (int nv = 0; nv < nvertex; ++nv)
2722              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2723           }                 
2724        }
2725      }       
2726    }
2727  }
2728
2729  /*!
2730    Receive latitude event from clients(s)
2731    \param[in] event event contain info about rank and associated latitude
2732  */
2733  void CDomain::recvLat(CEventServer& event)
2734  {
2735    string domainId;
2736    std::map<int, CBufferIn*> rankBuffers;
2737
2738    list<CEventServer::SSubEvent>::iterator it;
2739    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2740    {     
2741      CBufferIn* buffer = it->buffer;
2742      *buffer >> domainId;
2743      rankBuffers[it->rank] = buffer;   
2744    }
2745    get(domainId)->recvLat(rankBuffers);
2746  }
2747
2748  /*!
2749    Receive latitude information from client(s)
2750    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2751  */
2752  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
2753  {
2754    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2755    if (nbReceived != recvClientRanks_.size())
2756      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2757           << "The number of sending clients is not correct."
2758           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2759
2760    vector<CArray<double,1> > recvLatValue(nbReceived);
2761    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2762    for (i = 0; i < recvClientRanks_.size(); ++i)
2763    {
2764      int rank = recvClientRanks_[i];
2765      CBufferIn& buffer = *(rankBuffers[rank]);
2766      buffer >> hasLonLat;
2767      if (hasLonLat)
2768        buffer >> recvLatValue[i];
2769      buffer >> hasBounds;
2770      if (hasBounds)
2771        buffer >> recvBoundsLatValue[i];
2772    }
2773
2774    if (hasLonLat)
2775    {
2776      int nbLatInd = 0;
2777      for (i = 0; i < nbReceived; ++i)
2778      {
2779        nbLatInd += recvLatValue[i].numElements();
2780      }
2781   
2782      if (nbLatInd != globalLocalIndexMap_.size())
2783        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2784                << "something must be wrong with latitude index "<< std::endl;
2785
2786      nbLatInd = globalLocalIndexMap_.size();
2787      latvalue.resize(nbLatInd);
2788      if (hasBounds)
2789      {
2790        bounds_latvalue.resize(nvertex,nbLatInd);
2791        bounds_latvalue = 0. ;
2792      }
2793
2794      nbLatInd = 0;
2795      for (i = 0; i < nbReceived; ++i)
2796      {
2797        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2798        CArray<double,1>& tmp = recvLatValue[i];
2799        for (ind = 0; ind < tmp.numElements(); ++ind)
2800        {
2801          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2802          latvalue(lInd) = tmp(ind);   
2803           if (hasBounds)
2804           {
2805            CArray<double,2>& boundslat = recvBoundsLatValue[i];
2806            for (int nv = 0; nv < nvertex; ++nv)
2807              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2808           }   
2809          ++nbLatInd;
2810        }
2811      }       
2812    }
2813  }
2814
2815  /*!
2816    Receive area event from clients(s)
2817    \param[in] event event contain info about rank and associated area
2818  */
2819  void CDomain::recvArea(CEventServer& event)
2820  {
2821    string domainId;
2822    std::map<int, CBufferIn*> rankBuffers;
2823
2824    list<CEventServer::SSubEvent>::iterator it;
2825    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2826    {     
2827      CBufferIn* buffer = it->buffer;
2828      *buffer >> domainId;
2829      rankBuffers[it->rank] = buffer;     
2830    }
2831    get(domainId)->recvArea(rankBuffers);
2832  }
2833
2834  /*!
2835    Receive area information from client(s)
2836    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2837  */
2838  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
2839  {
2840    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2841    if (nbReceived != recvClientRanks_.size())
2842      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2843           << "The number of sending clients is not correct."
2844           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2845
2846    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2847    for (i = 0; i < recvClientRanks_.size(); ++i)
2848    {
2849      int rank = recvClientRanks_[i];
2850      CBufferIn& buffer = *(rankBuffers[rank]);     
2851      buffer >> hasArea;
2852      if (hasArea)
2853        buffer >> recvAreaValue[i];
2854    }
2855
2856    if (hasArea)
2857    {
2858      int nbAreaInd = 0;
2859      for (i = 0; i < nbReceived; ++i)
2860      {     
2861        nbAreaInd += recvAreaValue[i].numElements();
2862      }
2863
2864      if (nbAreaInd != globalLocalIndexMap_.size())
2865        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
2866                 << "something must be wrong with area index "<< std::endl;
2867
2868      nbAreaInd = globalLocalIndexMap_.size();
2869      areavalue.resize(nbAreaInd);
2870      nbAreaInd = 0;     
2871      for (i = 0; i < nbReceived; ++i)
2872      {
2873        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2874        CArray<double,1>& tmp = recvAreaValue[i];
2875        for (ind = 0; ind < tmp.numElements(); ++ind)
2876        {
2877          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2878          areavalue(lInd) = tmp(ind);         
2879        }
2880      }
2881     
2882    }
2883  }
2884
2885  /*!
2886    Compare two domain objects.
2887    They are equal if only if they have identical attributes as well as their values.
2888    Moreover, they must have the same transformations.
2889  \param [in] domain Compared domain
2890  \return result of the comparison
2891  */
2892  bool CDomain::isEqual(CDomain* obj)
2893  {
2894    vector<StdString> excludedAttr;
2895    excludedAttr.push_back("domain_ref");
2896    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2897    if (!objEqual) return objEqual;
2898
2899    TransMapTypes thisTrans = this->getAllTransformations();
2900    TransMapTypes objTrans  = obj->getAllTransformations();
2901
2902    TransMapTypes::const_iterator it, itb, ite;
2903    std::vector<ETranformationType> thisTransType, objTransType;
2904    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2905      thisTransType.push_back(it->first);
2906    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2907      objTransType.push_back(it->first);
2908
2909    if (thisTransType.size() != objTransType.size()) return false;
2910    for (int idx = 0; idx < thisTransType.size(); ++idx)
2911      objEqual &= (thisTransType[idx] == objTransType[idx]);
2912
2913    return objEqual;
2914  }
2915
2916  /*!
2917    Receive data index event from clients(s)
2918    \param[in] event event contain info about rank and associated index
2919  */
2920  void CDomain::recvDataIndex(CEventServer& event)
2921  {
2922    string domainId;
2923    std::map<int, CBufferIn*> rankBuffers;
2924
2925    list<CEventServer::SSubEvent>::iterator it;
2926    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2927    {     
2928      CBufferIn* buffer = it->buffer;
2929      *buffer >> domainId;
2930      rankBuffers[it->rank] = buffer;       
2931    }
2932    get(domainId)->recvDataIndex(rankBuffers);
2933  }
2934
2935  /*!
2936    Receive data index information from client(s)
2937    A client receives data index from different clients to rebuild its own data index.
2938    Because we use global index + mask info to calculate the sending data to client(s),
2939    this data index must be updated with mask info (maybe it will change in the future)
2940    Because the data index is local, to rebuild data index of received client, we should use global index along with.
2941
2942    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2943  */
2944  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
2945  {
2946    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
2947    if (nbReceived != recvClientRanks_.size())
2948      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
2949           << "The number of sending clients is not correct."
2950           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2951
2952    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
2953    for (i = 0; i < recvClientRanks_.size(); ++i)
2954    {
2955      int rank = recvClientRanks_[i];
2956      CBufferIn& buffer = *(rankBuffers[rank]);
2957      buffer >> recvDataIIndex[i];
2958      buffer >> recvDataJIndex[i];
2959    }
2960   
2961    int nbIndex = i_index.numElements();
2962    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2963    dataIIndex = -1; dataJIndex = -1;
2964     
2965    nbIndex = 0;
2966    for (i = 0; i < nbReceived; ++i)
2967    {     
2968      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2969      CArray<int,1>& tmpI = recvDataIIndex[i];   
2970      CArray<int,1>& tmpJ = recvDataJIndex[i];     
2971      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
2972          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
2973             << "The number of global received index is not coherent with the number of received data index."
2974             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
2975
2976      for (ind = 0; ind < tmpI.numElements(); ++ind)
2977      {
2978         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2979         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
2980         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
2981
2982//         if (!domainMask(lInd))   // Include mask info into data index on the RECEIVE getServerDimensionSizes
2983//         {
2984//           dataIIndex(lInd) = dataJIndex(lInd) = -1;
2985//         }
2986      } 
2987    }
2988
2989    int nbCompressedData = 0; 
2990    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
2991    {
2992       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
2993       if ((0 <= indexI) && (0 <= indexJ))
2994         ++nbCompressedData;
2995    }       
2996 
2997    data_i_index.resize(nbCompressedData);
2998    data_j_index.resize(nbCompressedData);
2999
3000    nbCompressedData = 0; 
3001    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3002    {
3003       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3004       if ((0 <= indexI) && (0 <= indexJ))
3005       {
3006          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3007          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3008         ++nbCompressedData;
3009       }
3010    }
3011
3012    // Reset data_ibegin, data_jbegin
3013    data_ibegin.setValue(0);
3014    data_jbegin.setValue(0);
3015  }
3016
3017  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3018  {
3019    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3020    return transformationMap_.back().second;
3021  }
3022
3023  /*!
3024    Check whether a domain has transformation
3025    \return true if domain has transformation
3026  */
3027  bool CDomain::hasTransformation()
3028  {
3029    return (!transformationMap_.empty());
3030  }
3031
3032  /*!
3033    Set transformation for current domain. It's the method to move transformation in hierarchy
3034    \param [in] domTrans transformation on domain
3035  */
3036  void CDomain::setTransformations(const TransMapTypes& domTrans)
3037  {
3038    transformationMap_ = domTrans;
3039  }
3040
3041  /*!
3042    Get all transformation current domain has
3043    \return all transformation
3044  */
3045  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3046  {
3047    return transformationMap_;
3048  }
3049
3050  void CDomain::duplicateTransformation(CDomain* src)
3051  {
3052    if (src->hasTransformation())
3053    {
3054      this->setTransformations(src->getAllTransformations());
3055    }
3056  }
3057
3058  /*!
3059   * Go through the hierarchy to find the domain from which the transformations must be inherited
3060   */
3061  void CDomain::solveInheritanceTransformation()
3062  {
3063    if (hasTransformation() || !hasDirectDomainReference())
3064      return;
3065
3066    CDomain* domain = this;
3067    std::vector<CDomain*> refDomains;
3068    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3069    {
3070      refDomains.push_back(domain);
3071      domain = domain->getDirectDomainReference();
3072    }
3073
3074    if (domain->hasTransformation())
3075      for (size_t i = 0; i < refDomains.size(); ++i)
3076        refDomains[i]->setTransformations(domain->getAllTransformations());
3077  }
3078
3079  void CDomain::setContextClient(CContextClient* contextClient)
3080  {
3081    if (clientsSet.find(contextClient)==clientsSet.end())
3082    {
3083      clients.push_back(contextClient) ;
3084      clientsSet.insert(contextClient);
3085    }
3086  }
3087
3088  /*!
3089    Parse children nodes of a domain in xml file.
3090    Whenver there is a new transformation, its type and name should be added into this function
3091    \param node child node to process
3092  */
3093  void CDomain::parse(xml::CXMLNode & node)
3094  {
3095    SuperClass::parse(node);
3096
3097    if (node.goToChildElement())
3098    {
3099      StdString nodeElementName;
3100      do
3101      {
3102        StdString nodeId("");
3103        if (node.getAttributes().end() != node.getAttributes().find("id"))
3104        { nodeId = node.getAttributes()["id"]; }
3105
3106        nodeElementName = node.getElementName();
3107        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3108        it = transformationMapList_.find(nodeElementName);
3109        if (ite != it)
3110        {
3111          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3112                                                                                                                nodeId,
3113                                                                                                                &node)));
3114        }
3115        else
3116        {
3117          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3118                << "The transformation " << nodeElementName << " has not been supported yet.");
3119        }
3120      } while (node.goToNextElement()) ;
3121      node.goToParentElement();
3122    }
3123  }
3124   //----------------------------------------------------------------
3125
3126   DEFINE_REF_FUNC(Domain,domain)
3127
3128   ///---------------------------------------------------------------
3129
3130} // namespace xios
Note: See TracBrowser for help on using the repository browser.