source: XIOS/branchs/xios-2.5/src/node/domain.cpp @ 1624

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

Bugfix for NEMO-like land processes elimination.

Values describing a domain (latitude, longitude, bounds,...) are set to -1 on land processes.

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