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

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

Solve issues for grid mask on server side.

YM

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