source: XIOS/dev/dev_trunk_omp/src/node/domain.cpp @ 1661

Last change on this file since 1661 was 1661, checked in by yushan, 5 years ago

MARK: branch merged with trunk @1660. Test (test_complete, test_remap) on ADA with IntelMPI and _usingEP/_usingMPI as switch.

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