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

Last change on this file since 1336 was 1336, checked in by oabramkina, 5 years ago

Changes in case of different number of processes per secondary-server pool.

No incorrect behaviour had been found in the previous versions; this commit is just to play safe.

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