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

Last change on this file since 1930 was 1930, 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".

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