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

Last change on this file since 1545 was 1545, checked in by yushan, 6 years ago

branch_openmp merged with trunk r1544

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