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

Last change on this file since 1325 was 1325, checked in by ymipsl, 7 years ago

Fix problem in domain mask introduced in rev. 1313

YM

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