source: XIOS/dev/branch_openmp/src/node/domain.cpp @ 1491

Last change on this file since 1491 was 1491, checked in by yushan, 3 years ago

Branch EP merged with Dev_cmip6 @r1490

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