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

Last change on this file since 1900 was 1900, checked in by yushan, 4 years ago

xios-2.5 : minormodif for Jean-Zay after system maintenance to RH8 and complr update frr om 8 à to 8.3.1. To be checked

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