source: XIOS/dev/dev_ym/XIOS_COUPLING/src/node/domain.cpp @ 1938

Last change on this file since 1938 was 1938, checked in by ymipsl, 4 years ago

XIOS Coupling branch : Solve spurious situation :

  • when client have no data on their local grid or local grid is 0 sized
  • when sever have no data on their local grid or local grid is 0 sized
  • holes in grid (missing global point) that cover a full client or a full server

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