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

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

Big update on on going work related to data distribution and transfer between clients and servers.
Revisite of the source and store filter using "connectors".

-> inputs work again

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.0 KB
Line 
1#include "domain.hpp"
2
3#include "attribute_template.hpp"
4#include "object_template.hpp"
5#include "group_template.hpp"
6
7#include "xios_spl.hpp"
8#include "event_client.hpp"
9#include "event_server.hpp"
10#include "buffer_in.hpp"
11#include "message.hpp"
12#include "type.hpp"
13#include "context.hpp"
14#include "context_client.hpp"
15#include "context_server.hpp"
16#include "array_new.hpp"
17#include "distribution_client.hpp"
18#include "server_distribution_description.hpp"
19#include "client_server_mapping_distributed.hpp"
20#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    }
2279    else if (type==EDistributionType::NONE) // domain is not distributed ie all servers get the same local domain
2280    {
2281      int nbServer = client->serverSize;
2282      int nglo=ni_glo*nj_glo ;
2283      CArray<size_t,1> indGlo ;
2284      for(size_t i=0;i<nglo;i++) indGlo(i) = i ;
2285      for (auto& rankServer : client->getRanksServerLeader()) globalIndex[rankServer] = indGlo ; 
2286    }
2287    remoteElement_[client] = new CDistributedElement(ni_glo*nj_glo, globalIndex) ;
2288    remoteElement_[client]->addFullView() ;
2289  }
2290  CATCH
2291
2292 
2293
2294  void CDomain::distributeToServer(CContextClient* client, map<int, CArray<size_t,1>>& globalIndex, const string& domainId)
2295  TRY
2296  {
2297    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2298    CContext* context = CContext::getCurrent();
2299
2300    this->sendAllAttributesToServer(client, serverDomainId)  ;
2301
2302    CDistributedElement scatteredElement(ni_glo*nj_glo, globalIndex) ;
2303    scatteredElement.addFullView() ;
2304    CScattererConnector scattererConnector(localElement_->getView(CElementView::FULL), scatteredElement.getView(CElementView::FULL), context->getIntraComm()) ;
2305    scattererConnector.computeConnector() ;
2306
2307    // phase 0
2308    // send remote element to construct the full view on server, ie without hole
2309    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2310    CMessage message0 ;
2311    message0<<serverDomainId<<0 ; 
2312    remoteElement_[client]->sendToServer(client,event0,message0) ; 
2313   
2314    // phase 1
2315    // send the full view of element to construct the connector which connect distributed data coming from client to the full local view
2316    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2317    CMessage message1 ;
2318    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2319    scattererConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2320   
2321    sendDistributedAttributes(client, scattererConnector, domainId) ;
2322
2323 
2324    // phase 2 send the mask : data index + mask2D
2325    CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2326    CArray<bool,1> maskOut ;
2327    CLocalConnector workflowToFull(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2328    workflowToFull.computeConnector() ;
2329    maskIn=true ;
2330    workflowToFull.transfer(maskIn,maskOut,false) ;
2331
2332
2333    // phase 3 : prepare grid scatterer connector to send data from client to server
2334    map<int,CArray<size_t,1>> workflowGlobalIndex ;
2335    map<int,CArray<bool,1>> maskOut2 ; 
2336    scattererConnector.transfer(maskOut, maskOut2, false) ;
2337    scatteredElement.addView(CElementView::WORKFLOW, maskOut2) ;
2338    scatteredElement.getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2339    // create new workflow view for scattered element
2340    CDistributedElement clientToServerElement(scatteredElement.getGlobalSize(), workflowGlobalIndex) ;
2341    clientToServerElement.addFullView() ;
2342    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2343    CMessage message2 ;
2344    message2<<serverDomainId<<2 ; 
2345    clientToServerElement.sendToServer(client, event2, message2) ; 
2346    clientToServerConnector_[client] = new CScattererConnector(localElement_->getView(CElementView::WORKFLOW),
2347                                                              clientToServerElement.getView(CElementView::FULL), context->getIntraComm()) ;
2348    clientToServerConnector_[client]->computeConnector() ;
2349
2350
2351    CEventClient event3(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2352    CMessage message3 ;
2353    message3<<serverDomainId<<3 ; 
2354    clientToServerConnector_[client]->transfer(maskIn,client,event3,message3) ; 
2355   
2356    clientFromServerConnector_[client] = new CGathererConnector(clientToServerElement.getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW));
2357    clientFromServerConnector_[client]->computeConnector() ;
2358
2359  }
2360  CATCH
2361 
2362  void CDomain::recvDomainDistribution(CEventServer& event)
2363  TRY
2364  {
2365    string domainId;
2366    int phasis ;
2367    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
2368    get(domainId)->receivedDomainDistribution(event, phasis);
2369  }
2370  CATCH
2371
2372  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2373  TRY
2374  {
2375    CContext* context = CContext::getCurrent();
2376    if (phasis==0) // receive the remote element to construct the full view
2377    {
2378      localElement_ = new  CLocalElement(context->getIntraCommRank(),event) ;
2379      localElement_->addFullView() ;
2380      // construct the local dimension and indexes
2381      auto& globalIndex=localElement_->getGlobalIndex() ;
2382      int nij=globalIndex.numElements() ;
2383      int minI=ni_glo,maxI=-1,minJ=nj_glo,maxJ=-1 ;
2384      int i,j ;
2385      int niGlo=ni_glo, njGlo=njGlo ;
2386      for(int ij=0;ij<nij;ij++)
2387      {
2388        j=globalIndex(ij)/niGlo ;
2389        i=globalIndex(ij)%niGlo ;
2390        if (i<minI) minI=i ;
2391        if (i>maxI) maxI=i ;
2392        if (j<minJ) minJ=j ;
2393        if (j>maxJ) maxJ=j ;
2394      } 
2395      if (maxI>=minI) { ibegin=minI ; ni=maxI-minI+1 ; }
2396      else {ibegin=0; ni=0 ;}
2397      if (maxJ>=minJ) { jbegin=minJ ; nj=maxJ-minJ+1 ; }
2398      else {jbegin=0; nj=0 ;}
2399
2400    }
2401    else if (phasis==1) // receive the sent view from client to construct the full distributed full view on server
2402    {
2403      CContext* context = CContext::getCurrent();
2404      CDistributedElement* elementFrom = new  CDistributedElement(event) ;
2405      elementFrom->addFullView() ;
2406      gathererConnector_ = new CGathererConnector(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2407      gathererConnector_->computeConnector() ; 
2408    }
2409    else if (phasis==2)
2410    {
2411      delete gathererConnector_ ;
2412      elementFrom_ = new  CDistributedElement(event) ;
2413      elementFrom_->addFullView() ;
2414      gathererConnector_ =  new CGathererConnector(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2415      gathererConnector_ -> computeConnector() ;
2416    }
2417    else if (phasis==3)
2418    {
2419      CArray<bool,1> localMask ;
2420      gathererConnector_->transfer(event,localMask,false) ;
2421      localElement_->addView(CElementView::WORKFLOW, localMask) ;
2422      mask_1d.reference(localMask.copy()) ;
2423 
2424      serverFromClientConnector_ = new CGathererConnector(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW)) ;
2425      serverFromClientConnector_->computeConnector() ;
2426     
2427      serverToClientConnector_ = new CScattererConnector(localElement_->getView(CElementView::WORKFLOW), elementFrom_->getView(CElementView::FULL),
2428                                                         context->getIntraComm()) ;
2429      serverToClientConnector_->computeConnector() ;
2430 
2431    }
2432  }
2433  CATCH
2434
2435
2436  void CDomain::sendDistributedAttributes(CContextClient* client, CScattererConnector& scattererConnector,  const string& domainId)
2437  {
2438    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2439    CContext* context = CContext::getCurrent();
2440
2441    if (hasLonLat)
2442    {
2443      { // send longitude
2444        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2445        CMessage message ;
2446        message<<serverDomainId<<string("lon") ; 
2447        scattererConnector.transfer(lonvalue, client, event,message) ;
2448      }
2449     
2450      { // send latitude
2451        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2452        CMessage message ;
2453        message<<serverDomainId<<string("lat") ; 
2454        scattererConnector.transfer(latvalue, client, event, message) ;
2455      }
2456    }
2457
2458    if (hasBounds)
2459    { 
2460      { // send longitude boudaries
2461        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2462        CMessage message ;
2463        message<<serverDomainId<<string("boundslon") ; 
2464        scattererConnector.transfer(nvertex, bounds_lonvalue, client, event, message ) ;
2465      }
2466
2467      { // send latitude boudaries
2468        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2469        CMessage message ;
2470        message<<serverDomainId<<string("boundslat") ; 
2471        scattererConnector.transfer(nvertex, bounds_latvalue, client, event, message ) ;
2472      }
2473    }
2474
2475    if (hasArea)
2476    {  // send area
2477      CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2478      CMessage message ;
2479      message<<serverDomainId<<string("area") ; 
2480      scattererConnector.transfer(areavalue, client, event,message) ;
2481    }
2482  }
2483
2484  void CDomain::recvDistributedAttributes(CEventServer& event)
2485  TRY
2486  {
2487    string domainId;
2488    string type ;
2489    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> type ;
2490    get(domainId)->recvDistributedAttributes(event, type);
2491  }
2492  CATCH
2493
2494  void CDomain::recvDistributedAttributes(CEventServer& event, const string& type)
2495  TRY
2496  {
2497    if (type=="lon") 
2498    {
2499      CArray<double,1> value ;
2500      gathererConnector_->transfer(event, value, 0.); 
2501      lonvalue_2d.resize(ni,nj) ;
2502      lonvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2503    }
2504    else if (type=="lat")
2505    {
2506      CArray<double,1> value ;
2507      gathererConnector_->transfer(event, value, 0.); 
2508      latvalue_2d.resize(ni,nj) ;
2509      latvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2510    }
2511    else if (type=="boundslon")
2512    {
2513      CArray<double,1> value ;
2514      gathererConnector_->transfer(event, nvertex, value, 0.); 
2515      bounds_lon_2d.resize(nvertex,ni,nj) ;
2516      bounds_lon_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2517    }
2518    else if (type=="boundslat")
2519    {
2520      CArray<double,1> value ;
2521      gathererConnector_->transfer(event, nvertex, value, 0.); 
2522      bounds_lat_2d.resize(nvertex,ni,nj) ;
2523      bounds_lat_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2524    }
2525    else if (type=="area") 
2526    {
2527      CArray<double,1> value ;
2528      gathererConnector_->transfer(event, value, 0.); 
2529      area.resize(ni,nj) ;
2530      area=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2531    }
2532  }
2533  CATCH
2534
2535  void CDomain::sendDomainDistribution(CContextClient* client, const string& domainId)
2536  TRY
2537  {
2538    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2539    CContext* context = CContext::getCurrent();
2540    int nbServer = client->serverSize;
2541    std::vector<int> nGlobDomain(2);
2542    nGlobDomain[0] = this->ni_glo;
2543    nGlobDomain[1] = this->nj_glo;
2544
2545    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2546    int distributedPosition ;
2547    if (isUnstructed_) distributedPosition = 0 ;
2548    else distributedPosition = 1 ;
2549
2550    serverDescription.computeServerDistribution(false, distributedPosition);
2551   
2552    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2553    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2554 
2555    vector<unordered_map<size_t,vector<int>>> indexServerOnElement ;
2556    CArray<int,1> axisDomainOrder(1) ; axisDomainOrder(0)=2 ;
2557    auto zeroIndex=serverDescription.computeServerGlobalByElement(indexServerOnElement, context->getIntraCommRank(), context->getIntraCommSize(),
2558                                                                  axisDomainOrder,distributedPosition) ;
2559    // distribution is very bad => to redo
2560    // convert indexServerOnElement => map<int,CArray<size_t,1>> - need to be changed later
2561    map<int, vector<size_t>> vectGlobalIndex ;
2562    for(auto& indexRanks : indexServerOnElement[0])
2563    {
2564      size_t index=indexRanks.first ;
2565      auto& ranks=indexRanks.second ;
2566      for(int rank : ranks) vectGlobalIndex[rank].push_back(index) ;
2567    }
2568    map<int, CArray<size_t,1>> globalIndex ;
2569    for(auto& vect : vectGlobalIndex ) globalIndex.emplace(vect.first, CArray<size_t,1>(vect.second.data(), shape(vect.second.size()))) ; 
2570
2571    CDistributedElement remoteElement(ni_glo*nj_glo, globalIndex) ;
2572    remoteElement.addFullView() ;
2573
2574    CRemoteConnector remoteConnector(localElement_->getView(CElementView::FULL), remoteElement.getView(CElementView::FULL),context->getIntraComm()) ;
2575    remoteConnector.computeConnector() ;
2576    CDistributedElement scatteredElement(remoteElement.getGlobalSize(), remoteConnector.getDistributedGlobalIndex()) ;
2577    scatteredElement.addFullView() ;
2578    CScattererConnector scatterConnector(localElement_->getView(CElementView::FULL), scatteredElement.getView(CElementView::FULL), context->getIntraComm()) ;
2579    scatterConnector.computeConnector() ;
2580    CGridScattererConnector gridScatter({&scatterConnector}) ;
2581
2582    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2583    CMessage message0 ;
2584    message0<<serverDomainId<<0 ; 
2585    remoteElement.sendToServer(client,event0,message0) ;
2586
2587    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2588    CMessage message1 ;
2589    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2590    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2591   
2592    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2593    CMessage message2 ;
2594    message2<<serverDomainId<<2 ; 
2595//    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2596    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2597
2598/*
2599    localElement_->getView(CElementView::FULL)->sendRemoteElement(remoteConnector, client, event1, message1) ;
2600    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2601    CMessage message2 ;
2602    message2<<serverDomainId<<2 ;
2603    remoteConnector.transferToServer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2604*/
2605
2606  }
2607  CATCH
2608 
2609
2610 
2611 
2612
2613  /*!
2614    Send all attributes from client to connected clients
2615    The attributes will be rebuilt on receiving side
2616  */
2617  // ym obsolete to be removed
2618  void CDomain::sendAttributes()
2619  TRY
2620  {
2621    //sendDistributionAttributes();
2622    //sendIndex();       
2623    //sendLonLat();
2624    //sendArea();   
2625    //sendDataIndex();
2626  }
2627  CATCH
2628  /*!
2629    Send global index from client to connected client(s)
2630  */
2631  void CDomain::sendIndex(CContextClient* client, const string& domainId)
2632  TRY
2633  {
2634    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2635   
2636    int ns, n, i, j, ind, nv, idx;
2637    int serverSize = client->serverSize;
2638    CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2639
2640    list<CMessage> list_msgsIndex;
2641    list<CArray<int,1> > list_indGlob;
2642
2643    std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2644    iteIndex = indSrv_[serverSize].end();
2645    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2646    {
2647      int nbIndGlob = 0;
2648      int rank = connectedServerRank_[serverSize][k];
2649      itIndex = indSrv_[serverSize].find(rank);
2650      if (iteIndex != itIndex)
2651        nbIndGlob = itIndex->second.size();
2652
2653      list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2654
2655      CArray<int,1>& indGlob = list_indGlob.back();
2656      for (n = 0; n < nbIndGlob; ++n) indGlob(n) = static_cast<int>(itIndex->second[n]);
2657
2658      list_msgsIndex.push_back(CMessage());
2659      list_msgsIndex.back() << serverDomainId << (int)type; // enum ne fonctionne pour les message => ToFix
2660      list_msgsIndex.back() << isCurvilinear;
2661      list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2662       
2663      eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2664    }
2665    client->sendEvent(eventIndex);
2666  }
2667  CATCH_DUMP_ATTR
2668
2669  /*!
2670    Send distribution from client to other clients
2671    Because a client in a level knows correctly the grid distribution of client on the next level
2672    it calculates this distribution then sends it to the corresponding clients on the next level
2673  */
2674  void CDomain::sendDistributionAttributes(CContextClient* client, const string& domainId)
2675  TRY
2676  {
2677    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2678
2679    int nbServer = client->serverSize;
2680    std::vector<int> nGlobDomain(2);
2681    nGlobDomain[0] = this->ni_glo;
2682    nGlobDomain[1] = this->nj_glo;
2683
2684    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2685    if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2686    else serverDescription.computeServerDistribution(false, 1);
2687
2688    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2689    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2690
2691    CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2692    if (client->isServerLeader())
2693    {
2694      std::list<CMessage> msgs;
2695
2696      const std::list<int>& ranks = client->getRanksServerLeader();
2697      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2698      {
2699        // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2700        const int ibegin_srv = serverIndexBegin[*itRank][0];
2701        const int jbegin_srv = serverIndexBegin[*itRank][1];
2702        const int ni_srv = serverDimensionSizes[*itRank][0];
2703        const int nj_srv = serverDimensionSizes[*itRank][1];
2704
2705        msgs.push_back(CMessage());
2706        CMessage& msg = msgs.back();
2707        msg << serverDomainId ;
2708        msg << isUnstructed_;
2709        msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2710        msg << ni_glo.getValue() << nj_glo.getValue();
2711        msg << isCompressible_;
2712
2713        event.push(*itRank,1,msg);
2714      }
2715      client->sendEvent(event);
2716    }
2717    else client->sendEvent(event);
2718  }
2719  CATCH_DUMP_ATTR
2720
2721  /*!
2722    Send area from client to connected client(s)
2723  */
2724  void CDomain::sendArea(CContextClient* client, const string& domainId)
2725  TRY
2726  {
2727    if (!hasArea) return;
2728    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2729   
2730    int ns, n, i, j, ind, nv, idx;
2731    int serverSize = client->serverSize;
2732
2733    // send area for each connected server
2734    CEventClient eventArea(getType(), EVENT_ID_AREA);
2735
2736    list<CMessage> list_msgsArea;
2737    list<CArray<double,1> > list_area;
2738
2739    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2740    iteMap = indSrv_[serverSize].end();
2741    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2742    {
2743      int nbData = 0;
2744      int rank = connectedServerRank_[serverSize][k];
2745      it = indSrv_[serverSize].find(rank);
2746      if (iteMap != it)
2747        nbData = it->second.size();
2748      list_area.push_back(CArray<double,1>(nbData));
2749
2750      const std::vector<size_t>& temp = it->second;
2751      for (n = 0; n < nbData; ++n)
2752      {
2753        idx = static_cast<int>(it->second[n]);
2754        list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2755      }
2756
2757      list_msgsArea.push_back(CMessage());
2758      list_msgsArea.back() <<serverDomainId << hasArea;
2759      list_msgsArea.back() << list_area.back();
2760      eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2761    }
2762    client->sendEvent(eventArea);
2763  }
2764  CATCH_DUMP_ATTR
2765
2766  /*!
2767    Send longitude and latitude from client to servers
2768    Each client send long and lat information to corresponding connected clients(s).
2769    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2770  */
2771  void CDomain::sendLonLat(CContextClient* client, const string& domainId)
2772  TRY
2773  {
2774    if (!hasLonLat) return;
2775    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2776   
2777    int ns, n, i, j, ind, nv, idx;
2778    int serverSize = client->serverSize;
2779
2780    // send lon lat for each connected server
2781    CEventClient eventLon(getType(), EVENT_ID_LON);
2782    CEventClient eventLat(getType(), EVENT_ID_LAT);
2783
2784    list<CMessage> list_msgsLon, list_msgsLat;
2785    list<CArray<double,1> > list_lon, list_lat;
2786    list<CArray<double,2> > list_boundslon, list_boundslat;
2787
2788    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2789    iteMap = indSrv_[serverSize].end();
2790    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2791    {
2792      int nbData = 0;
2793      int rank = connectedServerRank_[serverSize][k];
2794      it = indSrv_[serverSize].find(rank);
2795      if (iteMap != it)
2796        nbData = it->second.size();
2797
2798      list_lon.push_back(CArray<double,1>(nbData));
2799      list_lat.push_back(CArray<double,1>(nbData));
2800
2801      if (hasBounds)
2802      {
2803        list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2804        list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2805      }
2806
2807      CArray<double,1>& lon = list_lon.back();
2808      CArray<double,1>& lat = list_lat.back();
2809      const std::vector<size_t>& temp = it->second;
2810      for (n = 0; n < nbData; ++n)
2811      {
2812        idx = static_cast<int>(it->second[n]);
2813        int localInd = globalLocalIndexMap_[idx];
2814        lon(n) = lonvalue(localInd);
2815        lat(n) = latvalue(localInd);
2816
2817        if (hasBounds)
2818        {
2819          CArray<double,2>& boundslon = list_boundslon.back();
2820          CArray<double,2>& boundslat = list_boundslat.back();
2821
2822          for (nv = 0; nv < nvertex; ++nv)
2823          {
2824            boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2825            boundslat(nv, n) = bounds_latvalue(nv, localInd);
2826          }
2827        }
2828      }
2829
2830      list_msgsLon.push_back(CMessage());
2831      list_msgsLat.push_back(CMessage());
2832
2833      list_msgsLon.back() << serverDomainId << hasLonLat;
2834      if (hasLonLat) 
2835        list_msgsLon.back() << list_lon.back();
2836      list_msgsLon.back()  << hasBounds;
2837      if (hasBounds)
2838      {
2839        list_msgsLon.back() << list_boundslon.back();
2840      }
2841
2842      list_msgsLat.back() << serverDomainId << hasLonLat;
2843      if (hasLonLat)
2844        list_msgsLat.back() << list_lat.back();
2845      list_msgsLat.back() << hasBounds;
2846      if (hasBounds)
2847      {         
2848        list_msgsLat.back() << list_boundslat.back();
2849      }
2850
2851      eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2852      eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2853    }
2854    client->sendEvent(eventLon);
2855    client->sendEvent(eventLat);
2856  }
2857  CATCH_DUMP_ATTR
2858
2859  /*!
2860    Send data index to corresponding connected clients.
2861    Data index can be compressed however, we always send decompressed data index
2862    and they will be compressed on receiving.
2863    The compressed index are represented with 1 and others are represented with -1
2864  */
2865  void CDomain::sendDataIndex(CContextClient* client, const string& domainId)
2866  TRY
2867  {
2868    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2869   
2870    int ns, n, i, j, ind, nv, idx;
2871    int serverSize = client->serverSize;
2872
2873    // send area for each connected server
2874    CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2875
2876    list<CMessage> list_msgsDataIndex;
2877    list<CArray<int,1> > list_data_i_index, list_data_j_index;
2878
2879    int nbIndex = i_index.numElements();
2880    int niByIndex = max(i_index) - min(i_index) + 1;
2881    int njByIndex = max(j_index) - min(j_index) + 1; 
2882    int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2883    int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2884
2885   
2886    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2887    dataIIndex = -1; 
2888    dataJIndex = -1;
2889    ind = 0;
2890
2891    for (idx = 0; idx < data_i_index.numElements(); ++idx)
2892    {
2893      int dataIidx = data_i_index(idx) + data_ibegin;
2894      int dataJidx = data_j_index(idx) + data_jbegin;
2895      if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2896          (0 <= dataJidx) && (dataJidx < dataJindexBound))
2897      {
2898        dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2899        dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2900      }
2901    }
2902
2903    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2904    iteMap = indSrv_[serverSize].end();
2905    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2906    {
2907      int nbData = 0;
2908      int rank = connectedServerRank_[serverSize][k];
2909      it = indSrv_[serverSize].find(rank);
2910      if (iteMap != it)
2911        nbData = it->second.size();
2912      list_data_i_index.push_back(CArray<int,1>(nbData));
2913      list_data_j_index.push_back(CArray<int,1>(nbData));
2914
2915      const std::vector<size_t>& temp = it->second;
2916      for (n = 0; n < nbData; ++n)
2917      {
2918        idx = static_cast<int>(it->second[n]);
2919        i = globalLocalIndexMap_[idx];
2920        list_data_i_index.back()(n) = dataIIndex(i);
2921        list_data_j_index.back()(n) = dataJIndex(i);
2922      }
2923
2924      list_msgsDataIndex.push_back(CMessage());
2925      list_msgsDataIndex.back() << serverDomainId ;
2926      list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2927      eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2928    }
2929    client->sendEvent(eventDataIndex);
2930  }
2931  CATCH
2932 
2933  bool CDomain::dispatchEvent(CEventServer& event)
2934  TRY
2935  {
2936    if (SuperClass::dispatchEvent(event)) return true;
2937    else
2938    {
2939      switch(event.type)
2940      {
2941        case EVENT_ID_SERVER_ATTRIBUT:
2942          recvDistributionAttributes(event);
2943          return true;
2944          break;
2945        case EVENT_ID_INDEX:
2946          recvIndex(event);
2947          return true;
2948          break;
2949        case EVENT_ID_LON:
2950          recvLon(event);
2951          return true;
2952          break;
2953        case EVENT_ID_LAT:
2954          recvLat(event);
2955          return true;
2956          break;
2957        case EVENT_ID_AREA:
2958          recvArea(event);
2959          return true;
2960          break; 
2961        case EVENT_ID_DATA_INDEX:
2962          recvDataIndex(event);
2963          return true;
2964          break;
2965        case EVENT_ID_DOMAIN_DISTRIBUTION:
2966          recvDomainDistribution(event);
2967          return true;
2968          break;
2969        case EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE:
2970          recvDistributedAttributes(event);
2971          return true;
2972          break; 
2973        default:
2974          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2975                << "Unknown Event");
2976          return false;
2977       }
2978    }
2979  }
2980  CATCH
2981
2982  /*!
2983    Receive index event from clients(s)
2984    \param[in] event event contain info about rank and associated index
2985  */
2986  void CDomain::recvIndex(CEventServer& event)
2987  TRY
2988  {
2989    string domainId;
2990    std::map<int, CBufferIn*> rankBuffers;
2991
2992    list<CEventServer::SSubEvent>::iterator it;
2993    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2994    {     
2995      CBufferIn* buffer = it->buffer;
2996      *buffer >> domainId;
2997      rankBuffers[it->rank] = buffer;       
2998    }
2999    get(domainId)->recvIndex(rankBuffers);
3000  }
3001  CATCH
3002
3003  /*!
3004    Receive index information from client(s). We use the global index for mapping index between
3005    sending clients and receiving clients.
3006    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
3007  */
3008  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
3009  TRY
3010  {
3011    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
3012    recvClientRanks_.resize(nbReceived);       
3013
3014    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
3015    ind = 0;
3016    for (ind = 0; it != ite; ++it, ++ind)
3017    {       
3018       recvClientRanks_[ind] = it->first;
3019       CBufferIn& buffer = *(it->second);
3020       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
3021       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
3022    }
3023    int nbIndGlob = 0;
3024    for (i = 0; i < nbReceived; ++i)
3025    {
3026      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
3027    }
3028   
3029    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
3030    i_index.resize(nbIndGlob);
3031    j_index.resize(nbIndGlob);   
3032    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
3033
3034    nbIndGlob = 0;
3035    for (i = 0; i < nbReceived; ++i)
3036    {
3037      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
3038      for (ind = 0; ind < tmp.numElements(); ++ind)
3039      {
3040         index = tmp(ind);
3041         if (0 == globalLocalIndexMap_.count(index))
3042         {
3043           iIndex = (index%ni_glo)-ibegin;
3044           iIndex = (iIndex < 0) ? 0 : iIndex;
3045           jIndex = (index/ni_glo)-jbegin;
3046           jIndex = (jIndex < 0) ? 0 : jIndex;
3047           nbIndLoc = iIndex + ni * jIndex;
3048           i_index(nbIndGlob) = index % ni_glo;
3049           j_index(nbIndGlob) = index / ni_glo;
3050           globalLocalIndexMap_[index] = nbIndGlob;
3051           ++nbIndGlob;
3052         } 
3053      } 
3054    } 
3055
3056    if (nbIndGlob==0)
3057    {
3058      i_index.resize(nbIndGlob);
3059      j_index.resize(nbIndGlob);
3060    }
3061    else
3062    {
3063      i_index.resizeAndPreserve(nbIndGlob);
3064      j_index.resizeAndPreserve(nbIndGlob);
3065    }
3066
3067    domainMask.resize(0); // Mask is not defined anymore on servers
3068  }
3069  CATCH
3070
3071  /*!
3072    Receive attributes event from clients(s)
3073    \param[in] event event contain info about rank and associated attributes
3074  */
3075  void CDomain::recvDistributionAttributes(CEventServer& event)
3076  TRY
3077  {
3078    CBufferIn* buffer=event.subEvents.begin()->buffer;
3079    string domainId ;
3080    *buffer>>domainId ;
3081    get(domainId)->recvDistributionAttributes(*buffer);
3082  }
3083  CATCH
3084
3085  /*!
3086    Receive attributes from client(s)
3087    \param[in] rank rank of client source
3088    \param[in] buffer message containing attributes info
3089  */
3090  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
3091  TRY
3092  {
3093    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
3094    int ni_glo_tmp, nj_glo_tmp;
3095    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
3096           >> ni_glo_tmp >> nj_glo_tmp
3097           >> isCompressible_;
3098
3099    ni.setValue(ni_tmp);
3100    ibegin.setValue(ibegin_tmp);
3101    nj.setValue(nj_tmp);
3102    jbegin.setValue(jbegin_tmp);
3103    ni_glo.setValue(ni_glo_tmp);
3104    nj_glo.setValue(nj_glo_tmp);
3105
3106  }
3107 CATCH_DUMP_ATTR
3108  /*!
3109    Receive longitude event from clients(s)
3110    \param[in] event event contain info about rank and associated longitude
3111  */
3112  void CDomain::recvLon(CEventServer& event)
3113  TRY
3114  {
3115    string domainId;
3116    std::map<int, CBufferIn*> rankBuffers;
3117
3118    list<CEventServer::SSubEvent>::iterator it;
3119    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3120    {     
3121      CBufferIn* buffer = it->buffer;
3122      *buffer >> domainId;
3123      rankBuffers[it->rank] = buffer;       
3124    }
3125    get(domainId)->recvLon(rankBuffers);
3126  }
3127  CATCH
3128
3129  /*!
3130    Receive longitude information from client(s)
3131    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
3132  */
3133  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
3134  TRY
3135  {
3136    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
3137    if (nbReceived != recvClientRanks_.size())
3138      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
3139           << "The number of sending clients is not correct."
3140           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3141   
3142    int nbLonInd = 0;
3143    vector<CArray<double,1> > recvLonValue(nbReceived);
3144    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
3145    for (i = 0; i < recvClientRanks_.size(); ++i)
3146    {
3147      int rank = recvClientRanks_[i];
3148      CBufferIn& buffer = *(rankBuffers[rank]);
3149      buffer >> hasLonLat;
3150      if (hasLonLat)
3151        buffer >> recvLonValue[i];
3152      buffer >> hasBounds;
3153      if (hasBounds)
3154        buffer >> recvBoundsLonValue[i];
3155    }
3156
3157    if (hasLonLat)
3158    {
3159      for (i = 0; i < nbReceived; ++i)
3160      {
3161        nbLonInd += recvLonValue[i].numElements();
3162      }
3163   
3164      if (nbLonInd != globalLocalIndexMap_.size())
3165        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3166                 << "something must be wrong with longitude index "<< std::endl;
3167
3168      nbLonInd = globalLocalIndexMap_.size();
3169      lonvalue.resize(nbLonInd);
3170      if (hasBounds)
3171      {
3172        bounds_lonvalue.resize(nvertex,nbLonInd);
3173        bounds_lonvalue = 0.;
3174      }
3175
3176      nbLonInd = 0;
3177      for (i = 0; i < nbReceived; ++i)
3178      {
3179        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3180        CArray<double,1>& tmp = recvLonValue[i];
3181        for (ind = 0; ind < tmp.numElements(); ++ind)
3182        {
3183          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3184          lonvalue(lInd) = tmp(ind); 
3185           if (hasBounds)
3186           {         
3187            for (int nv = 0; nv < nvertex; ++nv)
3188              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
3189           }                 
3190        }
3191      }       
3192    }
3193
3194   // setup attribute depending the type of domain
3195    if (hasLonLat)
3196    {
3197      nbLonInd = globalLocalIndexMap_.size();
3198      if (ni*nj != nbLonInd)
3199        ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
3200             << "The number of index received is not coherent with the given resolution"
3201             << "nbLonInd=" << nbLonInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
3202
3203      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
3204      {
3205        lonvalue_2d.resize(ni,nj);
3206          for(int ij=0, j=0 ; j<nj ; j++)
3207            for(int i=0 ; i<ni; i++, ij++) lonvalue_2d(i,j) = lonvalue(ij) ;
3208       
3209        if (hasBounds)
3210        {
3211          bounds_lon_2d.resize(nvertex, ni, nj) ;
3212          for(int ij=0, j=0 ; j<nj ; j++)
3213            for(int i=0 ; i<ni; i++, ij++) 
3214              for(int nv=0; nv<nvertex; nv++) bounds_lon_2d(nv,i,j) = bounds_lonvalue(nv,ij) ;
3215        }
3216      }
3217      else if (type == type_attr::unstructured || type == type_attr::gaussian)
3218      {
3219        lonvalue_1d.resize(nbLonInd);
3220        lonvalue_1d = lonvalue ;
3221        if (hasBounds)
3222        {
3223          bounds_lon_1d.resize(nvertex, nbLonInd) ;
3224          bounds_lon_1d = bounds_lonvalue ;
3225        }
3226      }
3227    }
3228  }
3229  CATCH_DUMP_ATTR
3230
3231  /*!
3232    Receive latitude event from clients(s)
3233    \param[in] event event contain info about rank and associated latitude
3234  */
3235  void CDomain::recvLat(CEventServer& event)
3236  TRY
3237  {
3238    string domainId;
3239    std::map<int, CBufferIn*> rankBuffers;
3240
3241    list<CEventServer::SSubEvent>::iterator it;
3242    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3243    {     
3244      CBufferIn* buffer = it->buffer;
3245      *buffer >> domainId;
3246      rankBuffers[it->rank] = buffer;   
3247    }
3248    get(domainId)->recvLat(rankBuffers);
3249  }
3250  CATCH
3251
3252  /*!
3253    Receive latitude information from client(s)
3254    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
3255  */
3256  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
3257  TRY
3258  {
3259    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
3260    if (nbReceived != recvClientRanks_.size())
3261      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
3262           << "The number of sending clients is not correct."
3263           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3264
3265    vector<CArray<double,1> > recvLatValue(nbReceived);
3266    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
3267    for (i = 0; i < recvClientRanks_.size(); ++i)
3268    {
3269      int rank = recvClientRanks_[i];
3270      CBufferIn& buffer = *(rankBuffers[rank]);
3271      buffer >> hasLonLat;
3272      if (hasLonLat)
3273        buffer >> recvLatValue[i];
3274      buffer >> hasBounds;
3275      if (hasBounds)
3276        buffer >> recvBoundsLatValue[i];
3277    }
3278
3279    int nbLatInd = 0;
3280    if (hasLonLat)
3281    {
3282      for (i = 0; i < nbReceived; ++i)
3283      {
3284        nbLatInd += recvLatValue[i].numElements();
3285      }
3286   
3287      if (nbLatInd != globalLocalIndexMap_.size())
3288        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3289                << "something must be wrong with latitude index "<< std::endl;
3290
3291      nbLatInd = globalLocalIndexMap_.size();
3292      latvalue.resize(nbLatInd);
3293      if (hasBounds)
3294      {
3295        bounds_latvalue.resize(nvertex,nbLatInd);
3296        bounds_latvalue = 0. ;
3297      }
3298
3299      nbLatInd = 0;
3300      for (i = 0; i < nbReceived; ++i)
3301      {
3302        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3303        CArray<double,1>& tmp = recvLatValue[i];
3304        for (ind = 0; ind < tmp.numElements(); ++ind)
3305        {
3306          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3307          latvalue(lInd) = tmp(ind);   
3308           if (hasBounds)
3309           {
3310            CArray<double,2>& boundslat = recvBoundsLatValue[i];
3311            for (int nv = 0; nv < nvertex; ++nv)
3312              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
3313           }   
3314          ++nbLatInd;
3315        }
3316      }       
3317    }
3318      // setup attribute depending the type of domain
3319    if (hasLonLat)
3320    {
3321      nbLatInd = globalLocalIndexMap_.size();
3322     
3323      if (ni*nj != nbLatInd)
3324        ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
3325             << "The number of index received is not coherent with the given resolution"
3326            << "nbLonInd=" << nbLatInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
3327
3328      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
3329      {
3330        {
3331          latvalue_2d.resize(ni,nj);
3332          for(int ij=0, j=0 ; j<nj ; j++)
3333            for(int i=0 ; i<ni; i++, ij++) latvalue_2d(i,j) = latvalue(ij) ;
3334       
3335          if (hasBounds)
3336          {
3337            bounds_lat_2d.resize(nvertex, ni, nj) ;
3338            for(int ij=0, j=0 ; j<nj ; j++)
3339              for(int i=0 ; i<ni; i++, ij++) 
3340                for(int nv=0; nv<nvertex; nv++) bounds_lat_2d(nv,i,j) = bounds_latvalue(nv,ij) ;
3341          }
3342        }
3343      }
3344      else if (type == type_attr::unstructured || type == type_attr::gaussian)
3345      {
3346        if (hasLonLat)
3347        {
3348          latvalue_1d.resize(nbLatInd);
3349          latvalue_1d = latvalue ;
3350          if (hasBounds)
3351          {
3352            bounds_lat_1d.resize(nvertex, nbLatInd) ;
3353            bounds_lat_1d = bounds_latvalue ;
3354          }
3355        }
3356      }
3357    } 
3358  }
3359  CATCH_DUMP_ATTR
3360
3361  /*!
3362    Receive area event from clients(s)
3363    \param[in] event event contain info about rank and associated area
3364  */
3365  void CDomain::recvArea(CEventServer& event)
3366  TRY
3367  {
3368    string domainId;
3369    std::map<int, CBufferIn*> rankBuffers;
3370
3371    list<CEventServer::SSubEvent>::iterator it;
3372    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3373    {     
3374      CBufferIn* buffer = it->buffer;
3375      *buffer >> domainId;
3376      rankBuffers[it->rank] = buffer;     
3377    }
3378    get(domainId)->recvArea(rankBuffers);
3379  }
3380  CATCH
3381
3382  /*!
3383    Receive area information from client(s)
3384    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3385  */
3386  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
3387  TRY
3388  {
3389    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
3390    if (nbReceived != recvClientRanks_.size())
3391      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
3392           << "The number of sending clients is not correct."
3393           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3394
3395    vector<CArray<double,1> > recvAreaValue(nbReceived);     
3396    for (i = 0; i < recvClientRanks_.size(); ++i)
3397    {
3398      int rank = recvClientRanks_[i];
3399      CBufferIn& buffer = *(rankBuffers[rank]);     
3400      buffer >> hasArea;
3401      if (hasArea)
3402        buffer >> recvAreaValue[i];
3403    }
3404
3405    if (hasArea)
3406    {
3407      int nbAreaInd = 0;
3408      for (i = 0; i < nbReceived; ++i)
3409      {     
3410        nbAreaInd += recvAreaValue[i].numElements();
3411      }
3412
3413      if (nbAreaInd != globalLocalIndexMap_.size())
3414        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3415                 << "something must be wrong with area index "<< std::endl;
3416
3417      nbAreaInd = globalLocalIndexMap_.size();
3418      areavalue.resize(nbAreaInd);
3419      nbAreaInd = 0;     
3420      for (i = 0; i < nbReceived; ++i)
3421      {
3422        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3423        CArray<double,1>& tmp = recvAreaValue[i];
3424        for (ind = 0; ind < tmp.numElements(); ++ind)
3425        {
3426          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3427          areavalue(lInd) = tmp(ind);         
3428        }
3429      }
3430     
3431      // fill the area attribute
3432      area.resize(ni,nj);
3433      for (int j = 0; j < nj; ++j)
3434      {
3435        for (int i = 0; i < ni; ++i)
3436        {
3437          int k = j * ni + i;
3438          area(i,j) = areavalue(k) ;
3439        }
3440      }
3441    }
3442  }
3443  CATCH_DUMP_ATTR
3444
3445  /*!
3446    Compare two domain objects.
3447    They are equal if only if they have identical attributes as well as their values.
3448    Moreover, they must have the same transformations.
3449  \param [in] domain Compared domain
3450  \return result of the comparison
3451  */
3452  bool CDomain::isEqual(CDomain* obj)
3453  TRY
3454  {
3455    vector<StdString> excludedAttr;
3456    excludedAttr.push_back("domain_ref");
3457    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
3458    if (!objEqual) return objEqual;
3459
3460    TransMapTypes thisTrans = this->getAllTransformations();
3461    TransMapTypes objTrans  = obj->getAllTransformations();
3462
3463    TransMapTypes::const_iterator it, itb, ite;
3464    std::vector<ETranformationType> thisTransType, objTransType;
3465    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
3466      thisTransType.push_back(it->first);
3467    for (it = objTrans.begin(); it != objTrans.end(); ++it)
3468      objTransType.push_back(it->first);
3469
3470    if (thisTransType.size() != objTransType.size()) return false;
3471    for (int idx = 0; idx < thisTransType.size(); ++idx)
3472      objEqual &= (thisTransType[idx] == objTransType[idx]);
3473
3474    return objEqual;
3475  }
3476  CATCH_DUMP_ATTR
3477
3478  /*!
3479    Receive data index event from clients(s)
3480    \param[in] event event contain info about rank and associated index
3481  */
3482  void CDomain::recvDataIndex(CEventServer& event)
3483  TRY
3484  {
3485    string domainId;
3486    std::map<int, CBufferIn*> rankBuffers;
3487
3488    list<CEventServer::SSubEvent>::iterator it;
3489    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3490    {     
3491      CBufferIn* buffer = it->buffer;
3492      *buffer >> domainId;
3493      rankBuffers[it->rank] = buffer;       
3494    }
3495    get(domainId)->recvDataIndex(rankBuffers);
3496  }
3497  CATCH
3498
3499  /*!
3500    Receive data index information from client(s)
3501    A client receives data index from different clients to rebuild its own data index.
3502    Because we use global index + mask info to calculate the sending data to client(s),
3503    this data index must be updated with mask info (maybe it will change in the future)
3504    Because the data index is local, to rebuild data index of received client, we should use global index along with.
3505
3506    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3507  */
3508  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
3509  TRY
3510  {
3511    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
3512    if (nbReceived != recvClientRanks_.size())
3513      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3514           << "The number of sending clients is not correct."
3515           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3516
3517    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3518    for (i = 0; i < recvClientRanks_.size(); ++i)
3519    {
3520      int rank = recvClientRanks_[i];
3521      CBufferIn& buffer = *(rankBuffers[rank]);
3522      buffer >> recvDataIIndex[i];
3523      buffer >> recvDataJIndex[i];
3524    }
3525   
3526    int nbIndex = i_index.numElements();
3527    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3528    dataIIndex = -1; dataJIndex = -1;
3529     
3530    nbIndex = 0;
3531    for (i = 0; i < nbReceived; ++i)
3532    {     
3533      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3534      CArray<int,1>& tmpI = recvDataIIndex[i];   
3535      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3536      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3537          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3538             << "The number of global received index is not coherent with the number of received data index."
3539             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3540
3541      for (ind = 0; ind < tmpI.numElements(); ++ind)
3542      {
3543         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3544         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3545         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3546      } 
3547    }
3548
3549    int nbCompressedData = 0; 
3550    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3551    {
3552       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3553       if ((0 <= indexI) && (0 <= indexJ))
3554         ++nbCompressedData;
3555    }       
3556 
3557    data_i_index.resize(nbCompressedData);
3558    data_j_index.resize(nbCompressedData);
3559
3560    nbCompressedData = 0; 
3561    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3562    {
3563       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3564       if ((0 <= indexI) && (0 <= indexJ))
3565       {
3566          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3567          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3568         ++nbCompressedData;
3569       }
3570    }
3571
3572    // Reset data_ibegin, data_jbegin
3573    data_ibegin.setValue(0);
3574    data_jbegin.setValue(0);
3575  }
3576  CATCH_DUMP_ATTR
3577
3578  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3579  TRY
3580  {
3581    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3582    return transformationMap_.back().second;
3583  }
3584  CATCH_DUMP_ATTR
3585
3586  /*!
3587    Check whether a domain has transformation
3588    \return true if domain has transformation
3589  */
3590  bool CDomain::hasTransformation()
3591  TRY
3592  {
3593    return (!transformationMap_.empty());
3594  }
3595  CATCH_DUMP_ATTR
3596
3597  /*!
3598    Set transformation for current domain. It's the method to move transformation in hierarchy
3599    \param [in] domTrans transformation on domain
3600  */
3601  void CDomain::setTransformations(const TransMapTypes& domTrans)
3602  TRY
3603  {
3604    transformationMap_ = domTrans;
3605  }
3606  CATCH_DUMP_ATTR
3607
3608  /*!
3609    Get all transformation current domain has
3610    \return all transformation
3611  */
3612  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3613  TRY
3614  {
3615    return transformationMap_;
3616  }
3617  CATCH_DUMP_ATTR
3618
3619  void CDomain::duplicateTransformation(CDomain* src)
3620  TRY
3621  {
3622    if (src->hasTransformation())
3623    {
3624      this->setTransformations(src->getAllTransformations());
3625    }
3626  }
3627  CATCH_DUMP_ATTR
3628   
3629  /*!
3630   * Go through the hierarchy to find the domain from which the transformations must be inherited
3631   */
3632  void CDomain::solveInheritanceTransformation()
3633  TRY
3634  {
3635    if (hasTransformation() || !hasDirectDomainReference())
3636      return;
3637
3638    CDomain* domain = this;
3639    std::vector<CDomain*> refDomains;
3640    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3641    {
3642      refDomains.push_back(domain);
3643      domain = domain->getDirectDomainReference();
3644    }
3645
3646    if (domain->hasTransformation())
3647      for (size_t i = 0; i < refDomains.size(); ++i)
3648        refDomains[i]->setTransformations(domain->getAllTransformations());
3649  }
3650  CATCH_DUMP_ATTR
3651
3652  void CDomain::setContextClient(CContextClient* contextClient)
3653  TRY
3654  {
3655    if (clientsSet.find(contextClient)==clientsSet.end())
3656    {
3657      clients.push_back(contextClient) ;
3658      clientsSet.insert(contextClient);
3659    }
3660  }
3661  CATCH_DUMP_ATTR
3662
3663  /*!
3664    Parse children nodes of a domain in xml file.
3665    Whenver there is a new transformation, its type and name should be added into this function
3666    \param node child node to process
3667  */
3668  void CDomain::parse(xml::CXMLNode & node)
3669  TRY
3670  {
3671    SuperClass::parse(node);
3672
3673    if (node.goToChildElement())
3674    {
3675      StdString nodeElementName;
3676      do
3677      {
3678        StdString nodeId("");
3679        if (node.getAttributes().end() != node.getAttributes().find("id"))
3680        { nodeId = node.getAttributes()["id"]; }
3681
3682        nodeElementName = node.getElementName();
3683        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3684        it = transformationMapList_.find(nodeElementName);
3685        if (ite != it)
3686        {
3687          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3688                                                                                                                nodeId,
3689                                                                                                                &node)));
3690        }
3691        else
3692        {
3693          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3694                << "The transformation " << nodeElementName << " has not been supported yet.");
3695        }
3696      } while (node.goToNextElement()) ;
3697      node.goToParentElement();
3698    }
3699  }
3700  CATCH_DUMP_ATTR
3701   //----------------------------------------------------------------
3702
3703   DEFINE_REF_FUNC(Domain,domain)
3704
3705   ///---------------------------------------------------------------
3706
3707} // namespace xios
Note: See TracBrowser for help on using the repository browser.