source: XIOS/dev/XIOS_DEV_CMIP6/src/node/domain.cpp @ 1429

Last change on this file since 1429 was 1429, checked in by ymipsl, 3 years ago

Bug fix in computeLocalMask domain function.
Bad interaction with data index and domain index. When using domain defined with i_index and j_index, spurious missing data may apear after a transformation.

YM

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