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

Last change on this file since 1202 was 1202, checked in by mhnguyen, 5 years ago

Porting non-continuous axis zoom to dev branch

+) Port axis zoom
+) Resolve some merge conflicts
+) Revert some codes

Test
+) On Curie
+) Ok

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