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

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

Xios coupling

  • fix problem when sending grid mask from client to server
  • remove methods about grid compression, which will be managed in other ways

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: 132.0 KB
Line 
1#include "domain.hpp"
2
3#include "attribute_template.hpp"
4#include "object_template.hpp"
5#include "group_template.hpp"
6
7#include "xios_spl.hpp"
8#include "event_client.hpp"
9#include "event_server.hpp"
10#include "buffer_in.hpp"
11#include "message.hpp"
12#include "type.hpp"
13#include "context.hpp"
14#include "context_client.hpp"
15#include "context_server.hpp"
16#include "array_new.hpp"
17#include "distribution_client.hpp"
18#include "server_distribution_description.hpp"
19#include "client_server_mapping_distributed.hpp"
20#include "local_connector.hpp"
21#include "grid_local_connector.hpp"
22#include "remote_connector.hpp"
23#include "gatherer_connector.hpp"
24#include "scatterer_connector.hpp"
25#include "grid_scatterer_connector.hpp"
26#include "grid_gatherer_connector.hpp"
27
28
29
30
31#include <algorithm>
32
33namespace xios {
34
35   /// ////////////////////// Définitions ////////////////////// ///
36
37   CDomain::CDomain(void)
38      : CObjectTemplate<CDomain>(), CDomainAttributes()
39      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_()
40      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
41      , isClientAfterTransformationChecked(false), hasLonLat(false)
42      , isRedistributed_(false), hasPole(false)
43      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
44      , globalLocalIndexMap_(), computedWrittenIndex_(false)
45      , clients(), hasLatInReadFile_(false), hasBoundsLatInReadFile_(false)
46      , hasLonInReadFile_(false), hasBoundsLonInReadFile_(false)
47   {
48   }
49
50   CDomain::CDomain(const StdString & id)
51      : CObjectTemplate<CDomain>(id), CDomainAttributes()
52      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_() 
53      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
54      , isClientAfterTransformationChecked(false), hasLonLat(false)
55      , isRedistributed_(false), hasPole(false)
56      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
57      , globalLocalIndexMap_(), computedWrittenIndex_(false)
58      , clients(), hasLatInReadFile_(false), hasBoundsLatInReadFile_(false)
59      , hasLonInReadFile_(false), hasBoundsLonInReadFile_(false)
60   {
61    }
62
63   CDomain::~CDomain(void)
64   {
65   }
66
67   ///---------------------------------------------------------------
68
69   void CDomain::assignMesh(const StdString meshName, const int nvertex)
70   TRY
71   {
72     mesh = CMesh::getMesh(meshName, nvertex);
73   }
74   CATCH_DUMP_ATTR
75
76   CDomain* CDomain::createDomain()
77   TRY
78   {
79     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
80     return domain;
81   }
82   CATCH
83
84   std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
85   bool CDomain::_dummyTransformationMapList = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
86
87   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
88   TRY
89   {
90     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
91     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
92     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
93     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
94     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
95     m["reorder_domain"] = TRANS_REORDER_DOMAIN;
96     m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
97   }
98   CATCH
99
100   const std::set<StdString> & CDomain::getRelFiles(void) const
101   TRY
102   {
103      return (this->relFiles);
104   }
105   CATCH
106
107   /*!
108     Returns the number of indexes written by each server.
109     \return the number of indexes written by each server
110   */
111   int CDomain::getNumberWrittenIndexes(MPI_Comm writtenCom)
112   TRY
113   {
114     int writtenSize;
115     MPI_Comm_size(writtenCom, &writtenSize);
116     return numberWrittenIndexes_[writtenSize];
117   }
118   CATCH_DUMP_ATTR
119
120   /*!
121     Returns the total number of indexes written by the servers.
122     \return the total number of indexes written by the servers
123   */
124   int CDomain::getTotalNumberWrittenIndexes(MPI_Comm writtenCom)
125   TRY
126   {
127     int writtenSize;
128     MPI_Comm_size(writtenCom, &writtenSize);
129     return totalNumberWrittenIndexes_[writtenSize];
130   }
131   CATCH_DUMP_ATTR
132
133   /*!
134     Returns the offset of indexes written by each server.
135     \return the offset of indexes written by each server
136   */
137   int CDomain::getOffsetWrittenIndexes(MPI_Comm writtenCom)
138   TRY
139   {
140     int writtenSize;
141     MPI_Comm_size(writtenCom, &writtenSize);
142     return offsetWrittenIndexes_[writtenSize];
143   }
144   CATCH_DUMP_ATTR
145
146   CArray<int, 1>& CDomain::getCompressedIndexToWriteOnServer(MPI_Comm writtenCom)
147   TRY
148   {
149     int writtenSize;
150     MPI_Comm_size(writtenCom, &writtenSize);
151     return compressedIndexToWriteOnServer[writtenSize];
152   }
153   CATCH_DUMP_ATTR
154
155   //----------------------------------------------------------------
156
157   /*!
158    * Compute the minimum buffer size required to send the attributes to the server(s).
159    *
160    * \return A map associating the server rank with its minimum buffer size.
161    */
162   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
163   TRY
164   {
165
166     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
167
168     if (client->isServerLeader())
169     {
170       // size estimation for sendDistributionAttribut
171       size_t size = 11 * sizeof(size_t);
172
173       const std::list<int>& ranks = client->getRanksServerLeader();
174       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
175       {
176         if (size > attributesSizes[*itRank])
177           attributesSizes[*itRank] = size;
178       }
179     }
180
181     std::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->serverSize].end();
182     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
183     for (size_t k = 0; k < connectedServerRank_[client->serverSize].size(); ++k)
184     {
185       int rank = connectedServerRank_[client->serverSize][k];
186       std::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->serverSize].find(rank);
187       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
188
189       // size estimation for sendIndex (and sendArea which is always smaller or equal)
190       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
191
192       // size estimation for sendLonLat
193       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
194       if (hasBounds)
195         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
196
197       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
198       if (size > attributesSizes[rank])
199         attributesSizes[rank] = size;
200     }
201
202     return attributesSizes;
203   }
204   CATCH_DUMP_ATTR
205
206   //----------------------------------------------------------------
207
208   bool CDomain::isEmpty(void) const
209   TRY
210   {
211     return ((this->i_index.isEmpty()) || (0 == this->i_index.numElements()));
212   }
213   CATCH
214
215   //----------------------------------------------------------------
216
217   bool CDomain::IsWritten(const StdString & filename) const
218   TRY
219   {
220      return (this->relFiles.find(filename) != this->relFiles.end());
221   }
222   CATCH
223
224   bool CDomain::isWrittenCompressed(const StdString& filename) const
225   TRY
226   {
227      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
228   }
229   CATCH
230
231   //----------------------------------------------------------------
232
233   bool CDomain::isDistributed(void) const
234   TRY
235   {
236      bool distributed =  !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
237              (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
238      bool distributed_glo ;
239      distributed |= (1 == CContext::getCurrent()->intraCommSize_);
240
241      return distributed;
242   }
243   CATCH
244
245   //----------------------------------------------------------------
246
247   /*!
248    * Compute if the domain can be ouput in a compressed way.
249    * In this case the workflow view on server side must be the same
250    * than the full view for all context rank. The result is stored on
251    * internal isCompressible_ attribute.
252    */
253   void CDomain::computeIsCompressible(void)
254   TRY
255   {
256     // mesh is compressible contains some masked or indexed value, ie if full view is different of workflow view.
257     // But now assume that the size of the 2 view must be equal for everybody. True on server side
258     int isSameView = getLocalView(CElementView::FULL)->getSize() ==  getLocalView(CElementView::WORKFLOW)->getSize();
259     MPI_Allreduce(MPI_IN_PLACE, &isSameView, 1, MPI_INT, MPI_LAND, CContext::getCurrent()->getIntraComm()) ;
260     if (isSameView) isCompressible_ = false ;
261     else isCompressible_ = true ;
262     isCompressibleComputed_=true ;
263   }
264   CATCH
265
266   void CDomain::addRelFile(const StdString & filename)
267   TRY
268   {
269      this->relFiles.insert(filename);
270   }
271   CATCH_DUMP_ATTR
272
273   void CDomain::addRelFileCompressed(const StdString& filename)
274   TRY
275   {
276      this->relFilesCompressed.insert(filename);
277   }
278   CATCH_DUMP_ATTR
279
280   StdString CDomain::GetName(void)   { return (StdString("domain")); }
281   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
282   ENodeType CDomain::GetType(void)   { return (eDomain); }
283
284   //----------------------------------------------------------------
285
286   /*!
287      Verify if all distribution information of a domain are available
288      This checking verifies the definition of distribution attributes (ni, nj, ibegin, jbegin)
289   */
290   bool CDomain::distributionAttributesHaveValue() const
291   TRY
292   {
293      bool hasValues = true;
294
295      if (ni.isEmpty() && ibegin.isEmpty() && i_index.isEmpty())
296      {
297        hasValues = false;
298        return hasValues;
299      }
300
301      return hasValues;
302   }
303   CATCH
304
305   /*!
306     Redistribute RECTILINEAR or CURVILINEAR domain with a number of local domains.
307   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
308   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
309   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
310    \param [in] nbLocalDomain number of local domain on the domain destination
311   */
312   void CDomain::redistribute(int nbLocalDomain)
313   TRY
314   {
315     if (this->isRedistributed_) return;
316
317     this->isRedistributed_ = true;
318     CContext* context = CContext::getCurrent();
319     // For now the assumption is that secondary server pools consist of the same number of procs.
320     // CHANGE the line below if the assumption changes.
321
322     int rankClient = context->intraCommRank_;
323     int rankOnDomain = rankClient%nbLocalDomain;
324
325     if (ni_glo.isEmpty() || ni_glo <= 0 )
326     {
327        ERROR("CDomain::redistribute(int nbLocalDomain)",
328           << "[ Id = " << this->getId() << " ] "
329           << "The global domain is badly defined,"
330           << " check the \'ni_glo\'  value !")
331     }
332
333     if (nj_glo.isEmpty() || nj_glo <= 0 )
334     {
335        ERROR("CDomain::redistribute(int nbLocalDomain)",
336           << "[ Id = " << this->getId() << " ] "
337           << "The global domain is badly defined,"
338           << " check the \'nj_glo\'  value !")
339     }
340
341     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
342     {
343        int globalDomainSize = ni_glo * nj_glo;
344        if (globalDomainSize <= nbLocalDomain)
345        {
346          for (int idx = 0; idx < nbLocalDomain; ++idx)
347          {
348            if (rankOnDomain < globalDomainSize)
349            {
350              int iIdx = rankOnDomain % ni_glo;
351              int jIdx = rankOnDomain / ni_glo;
352              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
353              ni.setValue(1); nj.setValue(1);
354            }
355            else
356            {
357              ibegin.setValue(0); jbegin.setValue(0);
358              ni.setValue(0); nj.setValue(0);
359            }
360          }
361        }
362        else
363        {
364          float njGlo = nj_glo.getValue();
365          float niGlo = ni_glo.getValue();
366          int nbProcOnX, nbProcOnY, range;
367
368          // Compute (approximately) number of segment on x and y axis
369          float yOverXRatio = njGlo/niGlo;
370
371          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
372          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
373
374          // Simple distribution: Sweep from top to bottom, left to right
375          // Calculate local begin on x
376          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
377          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
378          for (int i = 1; i < nbProcOnX; ++i)
379          {
380            range = ni_glo / nbProcOnX;
381            if (i < (ni_glo%nbProcOnX)) ++range;
382            niVec[i-1] = range;
383            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
384          }
385          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
386
387          // Calculate local begin on y
388          for (int j = 1; j < nbProcOnY; ++j)
389          {
390            range = nj_glo / nbProcOnY;
391            if (j < (nj_glo%nbProcOnY)) ++range;
392            njVec[j-1] = range;
393            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
394          }
395          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
396
397          // Now assign value to ni, ibegin, nj, jbegin
398          int iIdx = rankOnDomain % nbProcOnX;
399          int jIdx = rankOnDomain / nbProcOnX;
400
401          if (rankOnDomain != (nbLocalDomain-1))
402          {
403            ibegin.setValue(ibeginVec[iIdx]);
404            jbegin.setValue(jbeginVec[jIdx]);
405            nj.setValue(njVec[jIdx]);
406            ni.setValue(niVec[iIdx]);
407          }
408          else // just merge all the remaining rectangle into the last one
409          {
410            ibegin.setValue(ibeginVec[iIdx]);
411            jbegin.setValue(jbeginVec[jIdx]);
412            nj.setValue(njVec[jIdx]);
413            ni.setValue(ni_glo - ibeginVec[iIdx]);
414          }
415        } 
416     }
417     else  // unstructured domain
418     {
419       if (this->i_index.isEmpty())
420       {
421          int globalDomainSize = ni_glo * nj_glo;
422          if (globalDomainSize <= nbLocalDomain)
423          {
424            for (int idx = 0; idx < nbLocalDomain; ++idx)
425            {
426              if (rankOnDomain < globalDomainSize)
427              {
428                int iIdx = rankOnDomain % ni_glo;
429                int jIdx = rankOnDomain / ni_glo;
430                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
431                ni.setValue(1); nj.setValue(1);
432              }
433              else
434              {
435                ibegin.setValue(0); jbegin.setValue(0);
436                ni.setValue(0); nj.setValue(0);
437              }
438            }
439          }
440          else
441          {
442            float njGlo = nj_glo.getValue();
443            float niGlo = ni_glo.getValue();
444            std::vector<int> ibeginVec(nbLocalDomain,0);
445            std::vector<int> niVec(nbLocalDomain);
446            for (int i = 1; i < nbLocalDomain; ++i)
447            {
448              int range = ni_glo / nbLocalDomain;
449              if (i < (ni_glo%nbLocalDomain)) ++range;
450              niVec[i-1] = range;
451              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
452            }
453            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
454
455            int iIdx = rankOnDomain % nbLocalDomain;
456            ibegin.setValue(ibeginVec[iIdx]);
457            jbegin.setValue(0);
458            ni.setValue(niVec[iIdx]);
459            nj.setValue(1);
460          }
461
462          i_index.resize(ni);         
463          for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
464        }
465        else
466        {
467          ibegin.setValue(this->i_index(0));
468          jbegin.setValue(0);
469          ni.setValue(this->i_index.numElements());
470          nj.setValue(1);
471        }
472     }
473
474     checkDomain();
475   }
476   CATCH_DUMP_ATTR
477
478   /*!
479     Fill in longitude and latitude whose values are read from file
480   */
481   void CDomain::fillInLonLat()
482   TRY
483   {
484     switch (type)
485     {
486      case type_attr::rectilinear:
487        fillInRectilinearLonLat();
488        break;
489      case type_attr::curvilinear:
490        fillInCurvilinearLonLat();
491        break;
492      case type_attr::unstructured:
493        fillInUnstructuredLonLat();
494        break;
495
496      default:
497      break;
498     }
499     completeLonLatClient() ;
500   }
501   CATCH_DUMP_ATTR
502
503   /*!
504     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
505     Range of longitude value from 0 - 360
506     Range of latitude value from -90 - +90
507   */
508   void CDomain::fillInRectilinearLonLat()
509   TRY
510   {
511     if (!lonvalue_rectilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
512     {
513       lonvalue_1d.resize(ni);
514       for (int idx = 0; idx < ni; ++idx)
515         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
516       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
517       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
518     }
519     else if (!hasLonInReadFile_)
520     {
521       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
522       lonvalue_1d.resize(ni);
523       double lonRange = lon_end - lon_start;
524       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
525
526        // Assign lon value
527       for (int i = 0; i < ni; ++i)
528       {
529         if (0 == (ibegin + i))
530         {
531           lonvalue_1d(i) = lon_start;
532         }
533         else if (ni_glo == (ibegin + i + 1))
534         {
535           lonvalue_1d(i) = lon_end;
536         }
537         else
538         {
539           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
540         }
541       }
542     }
543
544
545     if (!latvalue_rectilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
546     {
547       latvalue_1d.resize(nj);
548       for (int idx = 0; idx < nj; ++idx)
549         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
550       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
551       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
552     }
553     else if (!hasLatInReadFile_)
554     {
555       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
556       latvalue_1d.resize(nj);
557
558       double latRange = lat_end - lat_start;
559       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
560
561       for (int j = 0; j < nj; ++j)
562       {
563         if (0 == (jbegin + j))
564         {
565            latvalue_1d(j) = lat_start;
566         }
567         else if (nj_glo == (jbegin + j + 1))
568         {
569            latvalue_1d(j) = lat_end;
570         }
571         else
572         {
573           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
574         }
575       }
576     }
577   }
578   CATCH_DUMP_ATTR
579
580    /*
581      Fill in 2D longitude and latitude of curvilinear domain read from a file.
582      If there are already longitude and latitude defined by model. We just ignore read value.
583    */
584   void CDomain::fillInCurvilinearLonLat()
585   TRY
586   {
587     if (!lonvalue_curvilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
588     {
589       lonvalue_2d.resize(ni,nj);
590       for (int jdx = 0; jdx < nj; ++jdx)
591        for (int idx = 0; idx < ni; ++idx)
592         lonvalue_2d(idx,jdx) = lonvalue_curvilinear_read_from_file(idx, jdx);
593
594       lonvalue_curvilinear_read_from_file.free();
595     }
596
597     if (!latvalue_curvilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
598     {
599       latvalue_2d.resize(ni,nj);
600       for (int jdx = 0; jdx < nj; ++jdx)
601        for (int idx = 0; idx < ni; ++idx)
602           latvalue_2d(idx,jdx) = latvalue_curvilinear_read_from_file(idx, jdx);
603
604       latvalue_curvilinear_read_from_file.free();
605     }
606
607     if (!bounds_lonvalue_curvilinear_read_from_file.isEmpty() && bounds_lon_2d.isEmpty() && bounds_lon_1d.isEmpty())
608     {
609       bounds_lon_2d.resize(nvertex,ni,nj);
610       for (int jdx = 0; jdx < nj; ++jdx)
611        for (int idx = 0; idx < ni; ++idx)
612          for (int ndx = 0; ndx < nvertex; ++ndx)
613            bounds_lon_2d(ndx,idx,jdx) = bounds_lonvalue_curvilinear_read_from_file(ndx,idx, jdx);
614
615       bounds_lonvalue_curvilinear_read_from_file.free();
616     }
617
618     if (!bounds_latvalue_curvilinear_read_from_file.isEmpty() && bounds_lat_2d.isEmpty() && bounds_lat_1d.isEmpty())
619     {
620       bounds_lat_2d.resize(nvertex,ni,nj);
621       for (int jdx = 0; jdx < nj; ++jdx)
622        for (int idx = 0; idx < ni; ++idx)
623          for (int ndx = 0; ndx < nvertex; ++ndx)
624            bounds_lat_2d(ndx,idx,jdx) = bounds_latvalue_curvilinear_read_from_file(ndx,idx, jdx);
625
626       bounds_latvalue_curvilinear_read_from_file.free();
627     }
628   }
629   CATCH_DUMP_ATTR
630
631    /*
632      Fill in longitude and latitude of unstructured domain read from a file
633      If there are already longitude and latitude defined by model. We just igonore reading value.
634    */
635   void CDomain::fillInUnstructuredLonLat()
636   TRY
637   {
638     if (i_index.isEmpty())
639     {
640       i_index.resize(ni);
641       for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
642     }
643
644     if (!lonvalue_unstructured_read_from_file.isEmpty() && lonvalue_1d.isEmpty())
645     {
646        lonvalue_1d.resize(ni);
647        for (int idx = 0; idx < ni; ++idx)
648          lonvalue_1d(idx) = lonvalue_unstructured_read_from_file(idx);
649
650        // We dont need these values anymore, so just delete them
651        lonvalue_unstructured_read_from_file.free();
652     } 
653
654     if (!latvalue_unstructured_read_from_file.isEmpty() && latvalue_1d.isEmpty())
655     {
656        latvalue_1d.resize(ni);
657        for (int idx = 0; idx < ni; ++idx)
658          latvalue_1d(idx) =  latvalue_unstructured_read_from_file(idx);
659
660        // We dont need these values anymore, so just delete them
661        latvalue_unstructured_read_from_file.free();
662     }
663
664     if (!bounds_lonvalue_unstructured_read_from_file.isEmpty() && bounds_lon_1d.isEmpty())
665     {
666        int nbVertex = nvertex;
667        bounds_lon_1d.resize(nbVertex,ni);
668        for (int idx = 0; idx < ni; ++idx)
669          for (int jdx = 0; jdx < nbVertex; ++jdx)
670            bounds_lon_1d(jdx,idx) = bounds_lonvalue_unstructured_read_from_file(jdx, idx);
671
672        // We dont need these values anymore, so just delete them
673        bounds_lonvalue_unstructured_read_from_file.free();
674     }
675
676     if (!bounds_latvalue_unstructured_read_from_file.isEmpty() && bounds_lat_1d.isEmpty())
677     {
678        int nbVertex = nvertex;
679        bounds_lat_1d.resize(nbVertex,ni);
680        for (int idx = 0; idx < ni; ++idx)
681          for (int jdx = 0; jdx < nbVertex; ++jdx)
682            bounds_lat_1d(jdx,idx) = bounds_latvalue_unstructured_read_from_file(jdx, idx);
683
684        // We dont need these values anymore, so just delete them
685        bounds_latvalue_unstructured_read_from_file.free();
686     }
687   }
688   CATCH_DUMP_ATTR
689
690  /*
691    Get global longitude and latitude of rectilinear domain.
692  */
693   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g)
694   TRY
695   {
696     CContext* context = CContext::getCurrent();
697    // For now the assumption is that secondary server pools consist of the same number of procs.
698    // CHANGE the line below if the assumption changes.
699     int clientSize = context->intraCommSize_ ;
700     lon_g.resize(ni_glo) ;
701     lat_g.resize(nj_glo) ;
702
703
704     int* ibegin_g = new int[clientSize] ;
705     int* jbegin_g = new int[clientSize] ;
706     int* ni_g = new int[clientSize] ;
707     int* nj_g = new int[clientSize] ;
708     int v ;
709     v=ibegin ;
710     MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,context->intraComm_) ;
711     v=jbegin ;
712     MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,context->intraComm_) ;
713     v=ni ;
714     MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,context->intraComm_) ;
715     v=nj ;
716     MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,context->intraComm_) ;
717
718     MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,context->intraComm_) ;
719     MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,context->intraComm_) ;
720
721      delete[] ibegin_g ;
722      delete[] jbegin_g ;
723      delete[] ni_g ;
724      delete[] nj_g ;
725   }
726   CATCH_DUMP_ATTR
727
728   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
729                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
730   TRY
731  {
732     int i,j,k;
733
734     const int nvertexValue = 4;
735     boundsLon.resize(nvertexValue,ni*nj);
736
737     if (ni_glo>1)
738     {
739       double lonStepStart = lon(1)-lon(0);
740       bounds_lon_start=lon(0) - lonStepStart/2;
741       double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2);
742       bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
743       double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
744
745       // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
746       if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
747       {
748         bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
749         bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
750       }
751     }
752     else
753     {
754       if (bounds_lon_start.isEmpty()) bounds_lon_start=-180. ;
755       if (bounds_lon_end.isEmpty()) bounds_lon_end=180.-1e-8 ;
756     }
757
758     for(j=0;j<nj;++j)
759       for(i=0;i<ni;++i)
760       {
761         k=j*ni+i;
762         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start
763                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2;
764         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end
765                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2;
766       }
767
768
769    boundsLat.resize(nvertexValue,nj*ni);
770    bool isNorthPole=false ;
771    bool isSouthPole=false ;
772    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
773    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
774
775    // lat boundaries beyond pole the assimilate it to pole
776    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
777    if (nj_glo>1)
778    {
779      double latStepStart = lat(1)-lat(0);
780      if (isNorthPole) bounds_lat_start=lat(0);
781      else
782      {
783        bounds_lat_start=lat(0)-latStepStart/2;
784        if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
785        else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
786        else if (bounds_lat_start <= 90 && bounds_lat_start >= lat(0))
787        {
788          if ( std::abs(90-bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
789        }
790        else if (bounds_lat_start >= -90 && bounds_lat_start <= lat(0))
791        {
792          if ( std::abs(-90 - bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
793        }
794      }
795
796      double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
797      if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
798      else
799      {
800        bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
801
802        if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
803        else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
804        else if (bounds_lat_end <= 90 && bounds_lat_end >= lat(nj_glo-1))
805        {
806          if ( std::abs(90-bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
807        }
808        else if (bounds_lat_end >= -90 && bounds_lat_end <= lat(nj_glo-1))
809        {
810          if ( std::abs(-90 - bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
811        }
812      }
813    }
814    else
815    {
816      if (bounds_lat_start.isEmpty()) bounds_lat_start=-90. ;
817      if (bounds_lat_end.isEmpty()) bounds_lat_end=90 ;
818    }
819
820    for(j=0;j<nj;++j)
821      for(i=0;i<ni;++i)
822      {
823        k=j*ni+i;
824        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start
825                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2;
826        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end
827                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2;
828      }
829   }
830   CATCH_DUMP_ATTR
831
832   /*
833     General check of the domain to verify its mandatory attributes
834   */
835   void CDomain::checkDomain(void)
836   TRY
837   {
838     if (type.isEmpty())
839     {
840       ERROR("CDomain::checkDomain(void)",
841             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
842             << "The domain type is mandatory, "
843             << "please define the 'type' attribute.")
844     }
845
846     if (type == type_attr::gaussian) 
847     {
848        hasPole=true ;
849        type.setValue(type_attr::unstructured) ;
850      }
851      else if (type == type_attr::rectilinear) hasPole=true ;
852
853     if (type == type_attr::unstructured)
854     {
855        if (ni_glo.isEmpty())
856        {
857          ERROR("CDomain::checkDomain(void)",
858                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
859                << "The global domain is badly defined, "
860                << "the mandatory 'ni_glo' attribute is missing.")
861        }
862        else if (ni_glo <= 0)
863        {
864          ERROR("CDomain::checkDomain(void)",
865                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
866                << "The global domain is badly defined, "
867                << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
868        }
869        isUnstructed_ = true;
870        nj_glo = 1;
871        nj = 1;
872        jbegin = 0;
873        if (!i_index.isEmpty()) ni = i_index.numElements();
874        j_index.resize(ni);
875        for(int i=0;i<ni;++i) j_index(i)=0;
876
877        if (!area.isEmpty())
878          area.transposeSelf(1, 0);
879     }
880
881     if (ni_glo.isEmpty())
882     {
883       ERROR("CDomain::checkDomain(void)",
884             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
885             << "The global domain is badly defined, "
886             << "the mandatory 'ni_glo' attribute is missing.")
887     }
888     else if (ni_glo <= 0)
889     {
890       ERROR("CDomain::checkDomain(void)",
891             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
892             << "The global domain is badly defined, "
893             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
894     }
895
896     if (nj_glo.isEmpty())
897     {
898       ERROR("CDomain::checkDomain(void)",
899             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
900             << "The global domain is badly defined, "
901             << "the mandatory 'nj_glo' attribute is missing.")
902     }
903     else if (nj_glo <= 0)
904     {
905       ERROR("CDomain::checkDomain(void)",
906             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
907             << "The global domain is badly defined, "
908             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
909     }
910
911     checkLocalIDomain();
912     checkLocalJDomain();
913
914     if (i_index.isEmpty())
915     {
916       i_index.resize(ni*nj);
917       for (int j = 0; j < nj; ++j)
918         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
919     }
920
921     if (j_index.isEmpty())
922     {
923       j_index.resize(ni*nj);
924       for (int j = 0; j < nj; ++j)
925         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
926     }
927   }
928   CATCH_DUMP_ATTR
929
930   size_t CDomain::getGlobalWrittenSize(void)
931   {
932     return ni_glo*nj_glo ;
933   }
934   //----------------------------------------------------------------
935
936   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
937   void CDomain::checkLocalIDomain(void)
938   TRY
939   {
940      // If ibegin and ni are provided then we use them to check the validity of local domain
941      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
942      {
943        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
944        {
945          ERROR("CDomain::checkLocalIDomain(void)",
946                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
947                << "The local domain is wrongly defined,"
948                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
949        }
950      }
951
952      // i_index has higher priority than ibegin and ni
953      if (!i_index.isEmpty())
954      {
955        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
956        if (ni.isEmpty()) 
957        {         
958         // No information about ni
959          int minIndex = ni_glo - 1;
960          int maxIndex = 0;
961          for (int idx = 0; idx < i_index.numElements(); ++idx)
962          {
963            if (i_index(idx) < minIndex) minIndex = i_index(idx);
964            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
965          }
966          ni = maxIndex - minIndex + 1; 
967          minIIndex = minIIndex;         
968        }
969
970        // It's not so correct but if ibegin is not the first value of i_index
971        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
972        if (ibegin.isEmpty()) ibegin = minIIndex;
973      }
974      else if (ibegin.isEmpty() && ni.isEmpty())
975      {
976        ibegin = 0;
977        ni = ni_glo;
978      }
979      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
980      {
981        ERROR("CDomain::checkLocalIDomain(void)",
982              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
983              << "The local domain is wrongly defined," << endl
984              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
985              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
986      }
987       
988
989      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
990      {
991        ERROR("CDomain::checkLocalIDomain(void)",
992              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
993              << "The local domain is wrongly defined,"
994              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
995      }
996   }
997   CATCH_DUMP_ATTR
998
999   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
1000   void CDomain::checkLocalJDomain(void)
1001   TRY
1002   {
1003    // If jbegin and nj are provided then we use them to check the validity of local domain
1004     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
1005     {
1006       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
1007       {
1008         ERROR("CDomain::checkLocalJDomain(void)",
1009                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1010                << "The local domain is wrongly defined,"
1011                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1012       }
1013     }
1014
1015     if (!j_index.isEmpty())
1016     {
1017        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
1018        if (nj.isEmpty()) 
1019        {
1020          // No information about nj
1021          int minIndex = nj_glo - 1;
1022          int maxIndex = 0;
1023          for (int idx = 0; idx < j_index.numElements(); ++idx)
1024          {
1025            if (j_index(idx) < minIndex) minIndex = j_index(idx);
1026            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
1027          }
1028          nj = maxIndex - minIndex + 1;
1029          minJIndex = minIndex; 
1030        } 
1031        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
1032        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
1033       if (jbegin.isEmpty()) jbegin = minJIndex;       
1034     }
1035     else if (jbegin.isEmpty() && nj.isEmpty())
1036     {
1037       jbegin = 0;
1038       nj = nj_glo;
1039     }     
1040
1041
1042     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
1043     {
1044       ERROR("CDomain::checkLocalJDomain(void)",
1045              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1046              << "The local domain is wrongly defined,"
1047              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1048     }
1049   }
1050   CATCH_DUMP_ATTR
1051
1052   //----------------------------------------------------------------
1053
1054   void CDomain::checkMask(void)
1055   TRY
1056   {
1057      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
1058        ERROR("CDomain::checkMask(void)",
1059              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1060              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
1061              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
1062
1063      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
1064      {
1065        if (mask_1d.numElements() != i_index.numElements())
1066          ERROR("CDomain::checkMask(void)",
1067                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1068                << "'mask_1d' does not have the same size as the local domain." << std::endl
1069                << "Local size is " << i_index.numElements() << "." << std::endl
1070                << "Mask size is " << mask_1d.numElements() << ".");
1071      }
1072
1073      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
1074      {
1075        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1076          ERROR("CDomain::checkMask(void)",
1077                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1078                << "The mask does not have the same size as the local domain." << std::endl
1079                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1080                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
1081      }
1082
1083      if (!mask_2d.isEmpty())
1084      {
1085        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
1086        for (int j = 0; j < nj; ++j)
1087          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1088//        mask_2d.reset();
1089      }
1090      else if (mask_1d.isEmpty())
1091      {
1092        domainMask.resize(i_index.numElements());
1093        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
1094      }
1095      else
1096      {
1097      domainMask.resize(mask_1d.numElements());
1098      domainMask=mask_1d ;
1099     }
1100   }
1101   CATCH_DUMP_ATTR
1102
1103   //----------------------------------------------------------------
1104
1105   void CDomain::checkDomainData(void)
1106   TRY
1107   {
1108      if (data_dim.isEmpty())
1109      {
1110        data_dim.setValue(1);
1111      }
1112      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
1113      {
1114        ERROR("CDomain::checkDomainData(void)",
1115              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1116              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
1117      }
1118
1119      if (data_ibegin.isEmpty())
1120         data_ibegin.setValue(0);
1121      if (data_jbegin.isEmpty())
1122         data_jbegin.setValue(0);
1123
1124      if (data_ni.isEmpty())
1125      {
1126        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
1127      }
1128      else if (data_ni.getValue() < 0)
1129      {
1130        ERROR("CDomain::checkDomainData(void)",
1131              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1132              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
1133      }
1134
1135      if (data_nj.isEmpty())
1136      {
1137        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
1138      }
1139      else if (data_nj.getValue() < 0)
1140      {
1141        ERROR("CDomain::checkDomainData(void)",
1142              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1143              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1144      }
1145   }
1146   CATCH_DUMP_ATTR
1147
1148   //----------------------------------------------------------------
1149
1150   void CDomain::checkCompression(void)
1151   TRY
1152   {
1153     int i,j,ind;
1154      if (!data_i_index.isEmpty())
1155      {
1156        if (!data_j_index.isEmpty() &&
1157            data_j_index.numElements() != data_i_index.numElements())
1158        {
1159           ERROR("CDomain::checkCompression(void)",
1160                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1161                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1162                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1163                 << "'data_j_index' size = " << data_j_index.numElements());
1164        }
1165
1166        if (2 == data_dim)
1167        {
1168          if (data_j_index.isEmpty())
1169          {
1170             ERROR("CDomain::checkCompression(void)",
1171                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1172                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
1173          }
1174          for (int k=0; k<data_i_index.numElements(); ++k)
1175          {
1176            i = data_i_index(k)+data_ibegin ;
1177            j = data_j_index(k)+data_jbegin ;
1178            if (i>=0 && i<ni && j>=0 && j<nj)
1179            {
1180              ind=j*ni+i ;
1181              if (!domainMask(ind))
1182              {
1183                data_i_index(k) = -1;
1184                data_j_index(k) = -1;
1185              }
1186            }
1187            else
1188            {
1189              data_i_index(k) = -1;
1190              data_j_index(k) = -1;
1191            }
1192          }
1193        }
1194        else // (1 == data_dim)
1195        {
1196          if (data_j_index.isEmpty())
1197          {
1198            data_j_index.resize(data_ni);
1199            data_j_index = 0;
1200          }
1201          for (int k=0; k<data_i_index.numElements(); ++k)
1202          {
1203            i=data_i_index(k)+data_ibegin ;
1204            if (i>=0 && i < domainMask.size())
1205            {
1206              if (!domainMask(i)) data_i_index(k) = -1;
1207            }
1208            else
1209              data_i_index(k) = -1;
1210
1211            if (!domainMask(i)) data_i_index(k) = -1;
1212          }
1213        }
1214      }
1215      else
1216      {
1217        if (data_dim == 2 && !data_j_index.isEmpty())
1218          ERROR("CDomain::checkCompression(void)",
1219                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1220                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
1221
1222        if (1 == data_dim)
1223        {
1224          data_i_index.resize(data_ni);
1225          data_j_index.resize(data_ni);
1226          data_j_index = 0;
1227
1228          for (int k = 0; k < data_ni; ++k)
1229          {
1230            i=k+data_ibegin ;
1231            if (i>=0 && i < domainMask.size())
1232            {
1233              if (domainMask(i))
1234                data_i_index(k) = k;
1235              else
1236                data_i_index(k) = -1;
1237            }
1238            else
1239              data_i_index(k) = -1;
1240          }
1241        }
1242        else // (data_dim == 2)
1243        {
1244          const int dsize = data_ni * data_nj;
1245          data_i_index.resize(dsize);
1246          data_j_index.resize(dsize);
1247
1248          for(int count = 0, kj = 0; kj < data_nj; ++kj)
1249          {
1250            for(int ki = 0; ki < data_ni; ++ki, ++count)
1251            {
1252              i = ki + data_ibegin;
1253              j = kj + data_jbegin;
1254              ind=j*ni+i ;
1255              if (i>=0 && i<ni && j>=0 && j<nj)
1256              {
1257                if (domainMask(ind))
1258                {
1259                  data_i_index(count) = ki;
1260                  data_j_index(count) = kj;
1261                }
1262                else
1263                {
1264                  data_i_index(count) = -1;
1265                  data_j_index(count) = -1;
1266                }
1267              }
1268              else
1269              {
1270                data_i_index(count) = -1;
1271                data_j_index(count) = -1;
1272              }
1273            }
1274          }
1275        }
1276      }
1277   }
1278   CATCH_DUMP_ATTR
1279
1280   //----------------------------------------------------------------
1281   void CDomain::computeLocalMask(void)
1282   TRY
1283   {
1284     localMask.resize(i_index.numElements()) ;
1285     localMask=false ;
1286
1287     size_t dn=data_i_index.numElements() ;
1288     int i,j ;
1289     size_t k,ind ;
1290
1291     for(k=0;k<dn;k++)
1292     {
1293       if (data_dim==2)
1294       {
1295          i=data_i_index(k)+data_ibegin ;
1296          j=data_j_index(k)+data_jbegin ;
1297          if (i>=0 && i<ni && j>=0 && j<nj)
1298          {
1299            ind=j*ni+i ;
1300            localMask(ind)=domainMask(ind) ;
1301          }
1302       }
1303       else
1304       {
1305          i=data_i_index(k)+data_ibegin ;
1306          if (i>=0 && i<i_index.numElements())
1307          {
1308            ind=i ;
1309            localMask(ind)=domainMask(ind) ;
1310          }
1311       }
1312     }
1313   }
1314   CATCH_DUMP_ATTR
1315
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,
2298                                   CScattererConnector* &scattererConnector, const string& domainId)
2299  TRY
2300  {
2301    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2302    CContext* context = CContext::getCurrent();
2303
2304    this->sendAllAttributesToServer(client, serverDomainId)  ;
2305
2306    CDistributedElement scatteredElement(ni_glo*nj_glo, globalIndex) ;
2307    scatteredElement.addFullView() ;
2308    scattererConnector = new CScattererConnector(localElement_->getView(CElementView::FULL), scatteredElement.getView(CElementView::FULL), 
2309                                           context->getIntraComm(), client->getRemoteSize()) ;
2310    scattererConnector->computeConnector() ;
2311
2312    // phase 0
2313    // send remote element to construct the full view on server, ie without hole
2314    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2315    CMessage message0 ;
2316    message0<<serverDomainId<<0 ; 
2317    remoteElement_[client]->sendToServer(client,event0,message0) ; 
2318   
2319    // phase 1
2320    // send the full view of element to construct the connector which connect distributed data coming from client to the full local view
2321    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2322    CMessage message1 ;
2323    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2324    scattererConnector->transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2325   
2326    sendDistributedAttributes(client, *scattererConnector, domainId) ;
2327
2328 
2329    // phase 2 send the mask : data index + mask2D
2330    CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2331    CArray<bool,1> maskOut ;
2332    CLocalConnector workflowToFull(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2333    workflowToFull.computeConnector() ;
2334    maskIn=true ;
2335    workflowToFull.transfer(maskIn,maskOut,false) ;
2336
2337
2338    // phase 3 : prepare grid scatterer connector to send data from client to server
2339    map<int,CArray<size_t,1>> workflowGlobalIndex ;
2340    map<int,CArray<bool,1>> maskOut2 ; 
2341    scattererConnector->transfer(maskOut, maskOut2, false) ;
2342    scatteredElement.addView(CElementView::WORKFLOW, maskOut2) ;
2343    scatteredElement.getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2344    // create new workflow view for scattered element
2345    CDistributedElement clientToServerElement(scatteredElement.getGlobalSize(), workflowGlobalIndex) ;
2346    clientToServerElement.addFullView() ;
2347    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2348    CMessage message2 ;
2349    message2<<serverDomainId<<2 ; 
2350    clientToServerElement.sendToServer(client, event2, message2) ; 
2351    clientToServerConnector_[client] = new CScattererConnector(localElement_->getView(CElementView::WORKFLOW), clientToServerElement.getView(CElementView::FULL),
2352                                                               context->getIntraComm(), client->getRemoteSize()) ;
2353    clientToServerConnector_[client]->computeConnector() ;
2354
2355    clientFromServerConnector_[client] = new CGathererConnector(clientToServerElement.getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW));
2356    clientFromServerConnector_[client]->computeConnector() ;
2357
2358  }
2359  CATCH
2360 
2361  void CDomain::recvDomainDistribution(CEventServer& event)
2362  TRY
2363  {
2364    string domainId;
2365    int phasis ;
2366    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
2367    get(domainId)->receivedDomainDistribution(event, phasis);
2368  }
2369  CATCH
2370
2371  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2372  TRY
2373  {
2374    CContext* context = CContext::getCurrent();
2375    if (phasis==0) // receive the remote element to construct the full view
2376    {
2377      localElement_ = new  CLocalElement(context->getIntraCommRank(),event) ;
2378      localElement_->addFullView() ;
2379      // construct the local dimension and indexes
2380      auto& globalIndex=localElement_->getGlobalIndex() ;
2381      int nij=globalIndex.numElements() ;
2382      int minI=ni_glo,maxI=-1,minJ=nj_glo,maxJ=-1 ;
2383      int i,j ;
2384      int niGlo=ni_glo, njGlo=njGlo ;
2385      for(int ij=0;ij<nij;ij++)
2386      {
2387        j=globalIndex(ij)/niGlo ;
2388        i=globalIndex(ij)%niGlo ;
2389        if (i<minI) minI=i ;
2390        if (i>maxI) maxI=i ;
2391        if (j<minJ) minJ=j ;
2392        if (j>maxJ) maxJ=j ;
2393      } 
2394      if (maxI>=minI) { ibegin=minI ; ni=maxI-minI+1 ; }
2395      else {ibegin=0; ni=0 ;}
2396      if (maxJ>=minJ) { jbegin=minJ ; nj=maxJ-minJ+1 ; }
2397      else {jbegin=0; nj=0 ;}
2398
2399    }
2400    else if (phasis==1) // receive the sent view from client to construct the full distributed full view on server
2401    {
2402      CContext* context = CContext::getCurrent();
2403      CDistributedElement* elementFrom = new  CDistributedElement(event) ;
2404      elementFrom->addFullView() ;
2405      gathererConnector_ = new CGathererConnector(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2406      gathererConnector_->computeConnector() ; 
2407    }
2408    else if (phasis==2)
2409    {
2410//      delete gathererConnector_ ;
2411      elementFrom_ = new  CDistributedElement(event) ;
2412      elementFrom_->addFullView() ;
2413//      gathererConnector_ =  new CGathererConnector(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2414//      gathererConnector_ -> computeConnector() ;
2415    }
2416  }
2417  CATCH
2418
2419  void CDomain::setServerMask(CArray<bool,1>& serverMask, CContextClient* client)
2420  TRY
2421  {
2422    // nota : the client is needed to get the remote size for the scatterer connector. Maybe it is not the good place for this
2423    // Later, server to client connector can be computed on demand, with "client" as argument
2424    CContext* context = CContext::getCurrent();
2425    localElement_->addView(CElementView::WORKFLOW, serverMask) ;
2426    mask_1d.reference(serverMask.copy()) ;
2427 
2428    serverFromClientConnector_ = new CGathererConnector(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW)) ;
2429    serverFromClientConnector_->computeConnector() ;
2430     
2431    serverToClientConnector_ = new CScattererConnector(localElement_->getView(CElementView::WORKFLOW), elementFrom_->getView(CElementView::FULL),
2432                                                       context->getIntraComm(), client->getRemoteSize()) ;
2433    serverToClientConnector_->computeConnector() ;
2434  }
2435  CATCH_DUMP_ATTR
2436
2437
2438  void CDomain::sendDistributedAttributes(CContextClient* client, CScattererConnector& scattererConnector,  const string& domainId)
2439  {
2440    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2441    CContext* context = CContext::getCurrent();
2442
2443    if (hasLonLat)
2444    {
2445      { // send longitude
2446        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2447        CMessage message ;
2448        message<<serverDomainId<<string("lon") ; 
2449        scattererConnector.transfer(lonvalue, client, event,message) ;
2450      }
2451     
2452      { // send latitude
2453        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2454        CMessage message ;
2455        message<<serverDomainId<<string("lat") ; 
2456        scattererConnector.transfer(latvalue, client, event, message) ;
2457      }
2458    }
2459
2460    if (hasBounds)
2461    { 
2462      { // send longitude boudaries
2463        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2464        CMessage message ;
2465        message<<serverDomainId<<string("boundslon") ; 
2466        scattererConnector.transfer(nvertex, bounds_lonvalue, client, event, message ) ;
2467      }
2468
2469      { // send latitude boudaries
2470        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2471        CMessage message ;
2472        message<<serverDomainId<<string("boundslat") ; 
2473        scattererConnector.transfer(nvertex, bounds_latvalue, client, event, message ) ;
2474      }
2475    }
2476
2477    if (hasArea)
2478    {  // send area
2479      CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2480      CMessage message ;
2481      message<<serverDomainId<<string("area") ; 
2482      scattererConnector.transfer(areavalue, client, event,message) ;
2483    }
2484  }
2485
2486  void CDomain::recvDistributedAttributes(CEventServer& event)
2487  TRY
2488  {
2489    string domainId;
2490    string type ;
2491    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> type ;
2492    get(domainId)->recvDistributedAttributes(event, type);
2493  }
2494  CATCH
2495
2496  void CDomain::recvDistributedAttributes(CEventServer& event, const string& type)
2497  TRY
2498  {
2499    if (type=="lon") 
2500    {
2501      CArray<double,1> value ;
2502      gathererConnector_->transfer(event, value, 0.); 
2503      lonvalue_2d.resize(ni,nj) ;
2504      if (lonvalue_2d.numElements()>0) lonvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2505    }
2506    else if (type=="lat")
2507    {
2508      CArray<double,1> value ;
2509      gathererConnector_->transfer(event, value, 0.); 
2510      latvalue_2d.resize(ni,nj) ;
2511      if (latvalue_2d.numElements()>0) latvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2512    }
2513    else if (type=="boundslon")
2514    {
2515      CArray<double,1> value ;
2516      gathererConnector_->transfer(event, nvertex, value, 0.); 
2517      bounds_lon_2d.resize(nvertex,ni,nj) ;
2518      if (bounds_lon_2d.numElements()>0) bounds_lon_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2519    }
2520    else if (type=="boundslat")
2521    {
2522      CArray<double,1> value ;
2523      gathererConnector_->transfer(event, nvertex, value, 0.); 
2524      bounds_lat_2d.resize(nvertex,ni,nj) ;
2525      if (bounds_lat_2d.numElements()>0) bounds_lat_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2526    }
2527    else if (type=="area") 
2528    {
2529      CArray<double,1> value ;
2530      gathererConnector_->transfer(event, value, 0.); 
2531      area.resize(ni,nj) ;
2532      if (area.numElements()>0) area=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2533    }
2534  }
2535  CATCH
2536
2537  void CDomain::sendDomainDistribution(CContextClient* client, const string& domainId)
2538  TRY
2539  {
2540    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2541    CContext* context = CContext::getCurrent();
2542    int nbServer = client->serverSize;
2543    std::vector<int> nGlobDomain(2);
2544    nGlobDomain[0] = this->ni_glo;
2545    nGlobDomain[1] = this->nj_glo;
2546
2547    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2548    int distributedPosition ;
2549    if (isUnstructed_) distributedPosition = 0 ;
2550    else distributedPosition = 1 ;
2551
2552    serverDescription.computeServerDistribution(false, distributedPosition);
2553   
2554    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2555    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2556 
2557    vector<unordered_map<size_t,vector<int>>> indexServerOnElement ;
2558    CArray<int,1> axisDomainOrder(1) ; axisDomainOrder(0)=2 ;
2559    auto zeroIndex=serverDescription.computeServerGlobalByElement(indexServerOnElement, context->getIntraCommRank(), context->getIntraCommSize(),
2560                                                                  axisDomainOrder,distributedPosition) ;
2561    // distribution is very bad => to redo
2562    // convert indexServerOnElement => map<int,CArray<size_t,1>> - need to be changed later
2563    map<int, vector<size_t>> vectGlobalIndex ;
2564    for(auto& indexRanks : indexServerOnElement[0])
2565    {
2566      size_t index=indexRanks.first ;
2567      auto& ranks=indexRanks.second ;
2568      for(int rank : ranks) vectGlobalIndex[rank].push_back(index) ;
2569    }
2570    map<int, CArray<size_t,1>> globalIndex ;
2571    for(auto& vect : vectGlobalIndex ) globalIndex.emplace(vect.first, CArray<size_t,1>(vect.second.data(), shape(vect.second.size()))) ; 
2572
2573    CDistributedElement remoteElement(ni_glo*nj_glo, globalIndex) ;
2574    remoteElement.addFullView() ;
2575
2576    CRemoteConnector remoteConnector(localElement_->getView(CElementView::FULL), remoteElement.getView(CElementView::FULL),context->getIntraComm()) ;
2577    remoteConnector.computeConnector() ;
2578    CDistributedElement scatteredElement(remoteElement.getGlobalSize(), remoteConnector.getDistributedGlobalIndex()) ;
2579    scatteredElement.addFullView() ;
2580    CScattererConnector scatterConnector(localElement_->getView(CElementView::FULL), scatteredElement.getView(CElementView::FULL), 
2581                                         context->getIntraComm(), client->getRemoteSize()) ;
2582    scatterConnector.computeConnector() ;
2583    CGridScattererConnector gridScatter({&scatterConnector}) ;
2584
2585    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2586    CMessage message0 ;
2587    message0<<serverDomainId<<0 ; 
2588    remoteElement.sendToServer(client,event0,message0) ;
2589
2590    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2591    CMessage message1 ;
2592    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2593    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2594   
2595    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2596    CMessage message2 ;
2597    message2<<serverDomainId<<2 ; 
2598//    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2599    scatterConnector.transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2600
2601/*
2602    localElement_->getView(CElementView::FULL)->sendRemoteElement(remoteConnector, client, event1, message1) ;
2603    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2604    CMessage message2 ;
2605    message2<<serverDomainId<<2 ;
2606    remoteConnector.transferToServer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event2,message2) ;
2607*/
2608
2609  }
2610  CATCH
2611 
2612
2613 
2614 
2615
2616  /*!
2617    Send all attributes from client to connected clients
2618    The attributes will be rebuilt on receiving side
2619  */
2620  // ym obsolete to be removed
2621  void CDomain::sendAttributes()
2622  TRY
2623  {
2624    //sendDistributionAttributes();
2625    //sendIndex();       
2626    //sendLonLat();
2627    //sendArea();   
2628    //sendDataIndex();
2629  }
2630  CATCH
2631  /*!
2632    Send global index from client to connected client(s)
2633  */
2634  void CDomain::sendIndex(CContextClient* client, const string& domainId)
2635  TRY
2636  {
2637    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2638   
2639    int ns, n, i, j, ind, nv, idx;
2640    int serverSize = client->serverSize;
2641    CEventClient eventIndex(getType(), EVENT_ID_INDEX);
2642
2643    list<CMessage> list_msgsIndex;
2644    list<CArray<int,1> > list_indGlob;
2645
2646    std::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2647    iteIndex = indSrv_[serverSize].end();
2648    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2649    {
2650      int nbIndGlob = 0;
2651      int rank = connectedServerRank_[serverSize][k];
2652      itIndex = indSrv_[serverSize].find(rank);
2653      if (iteIndex != itIndex)
2654        nbIndGlob = itIndex->second.size();
2655
2656      list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2657
2658      CArray<int,1>& indGlob = list_indGlob.back();
2659      for (n = 0; n < nbIndGlob; ++n) indGlob(n) = static_cast<int>(itIndex->second[n]);
2660
2661      list_msgsIndex.push_back(CMessage());
2662      list_msgsIndex.back() << serverDomainId << (int)type; // enum ne fonctionne pour les message => ToFix
2663      list_msgsIndex.back() << isCurvilinear;
2664      list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2665       
2666      eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2667    }
2668    client->sendEvent(eventIndex);
2669  }
2670  CATCH_DUMP_ATTR
2671
2672  /*!
2673    Send distribution from client to other clients
2674    Because a client in a level knows correctly the grid distribution of client on the next level
2675    it calculates this distribution then sends it to the corresponding clients on the next level
2676  */
2677  void CDomain::sendDistributionAttributes(CContextClient* client, const string& domainId)
2678  TRY
2679  {
2680    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2681
2682    int nbServer = client->serverSize;
2683    std::vector<int> nGlobDomain(2);
2684    nGlobDomain[0] = this->ni_glo;
2685    nGlobDomain[1] = this->nj_glo;
2686
2687    CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2688    if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2689    else serverDescription.computeServerDistribution(false, 1);
2690
2691    std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2692    std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2693
2694    CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2695    if (client->isServerLeader())
2696    {
2697      std::list<CMessage> msgs;
2698
2699      const std::list<int>& ranks = client->getRanksServerLeader();
2700      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2701      {
2702        // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2703        const int ibegin_srv = serverIndexBegin[*itRank][0];
2704        const int jbegin_srv = serverIndexBegin[*itRank][1];
2705        const int ni_srv = serverDimensionSizes[*itRank][0];
2706        const int nj_srv = serverDimensionSizes[*itRank][1];
2707
2708        msgs.push_back(CMessage());
2709        CMessage& msg = msgs.back();
2710        msg << serverDomainId ;
2711        msg << isUnstructed_;
2712        msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2713        msg << ni_glo.getValue() << nj_glo.getValue();
2714        msg << isCompressible_;
2715
2716        event.push(*itRank,1,msg);
2717      }
2718      client->sendEvent(event);
2719    }
2720    else client->sendEvent(event);
2721  }
2722  CATCH_DUMP_ATTR
2723
2724  /*!
2725    Send area from client to connected client(s)
2726  */
2727  void CDomain::sendArea(CContextClient* client, const string& domainId)
2728  TRY
2729  {
2730    if (!hasArea) return;
2731    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2732   
2733    int ns, n, i, j, ind, nv, idx;
2734    int serverSize = client->serverSize;
2735
2736    // send area for each connected server
2737    CEventClient eventArea(getType(), EVENT_ID_AREA);
2738
2739    list<CMessage> list_msgsArea;
2740    list<CArray<double,1> > list_area;
2741
2742    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2743    iteMap = indSrv_[serverSize].end();
2744    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2745    {
2746      int nbData = 0;
2747      int rank = connectedServerRank_[serverSize][k];
2748      it = indSrv_[serverSize].find(rank);
2749      if (iteMap != it)
2750        nbData = it->second.size();
2751      list_area.push_back(CArray<double,1>(nbData));
2752
2753      const std::vector<size_t>& temp = it->second;
2754      for (n = 0; n < nbData; ++n)
2755      {
2756        idx = static_cast<int>(it->second[n]);
2757        list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2758      }
2759
2760      list_msgsArea.push_back(CMessage());
2761      list_msgsArea.back() <<serverDomainId << hasArea;
2762      list_msgsArea.back() << list_area.back();
2763      eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
2764    }
2765    client->sendEvent(eventArea);
2766  }
2767  CATCH_DUMP_ATTR
2768
2769  /*!
2770    Send longitude and latitude from client to servers
2771    Each client send long and lat information to corresponding connected clients(s).
2772    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2773  */
2774  void CDomain::sendLonLat(CContextClient* client, const string& domainId)
2775  TRY
2776  {
2777    if (!hasLonLat) return;
2778    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2779   
2780    int ns, n, i, j, ind, nv, idx;
2781    int serverSize = client->serverSize;
2782
2783    // send lon lat for each connected server
2784    CEventClient eventLon(getType(), EVENT_ID_LON);
2785    CEventClient eventLat(getType(), EVENT_ID_LAT);
2786
2787    list<CMessage> list_msgsLon, list_msgsLat;
2788    list<CArray<double,1> > list_lon, list_lat;
2789    list<CArray<double,2> > list_boundslon, list_boundslat;
2790
2791    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2792    iteMap = indSrv_[serverSize].end();
2793    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2794    {
2795      int nbData = 0;
2796      int rank = connectedServerRank_[serverSize][k];
2797      it = indSrv_[serverSize].find(rank);
2798      if (iteMap != it)
2799        nbData = it->second.size();
2800
2801      list_lon.push_back(CArray<double,1>(nbData));
2802      list_lat.push_back(CArray<double,1>(nbData));
2803
2804      if (hasBounds)
2805      {
2806        list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2807        list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2808      }
2809
2810      CArray<double,1>& lon = list_lon.back();
2811      CArray<double,1>& lat = list_lat.back();
2812      const std::vector<size_t>& temp = it->second;
2813      for (n = 0; n < nbData; ++n)
2814      {
2815        idx = static_cast<int>(it->second[n]);
2816        int localInd = globalLocalIndexMap_[idx];
2817        lon(n) = lonvalue(localInd);
2818        lat(n) = latvalue(localInd);
2819
2820        if (hasBounds)
2821        {
2822          CArray<double,2>& boundslon = list_boundslon.back();
2823          CArray<double,2>& boundslat = list_boundslat.back();
2824
2825          for (nv = 0; nv < nvertex; ++nv)
2826          {
2827            boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2828            boundslat(nv, n) = bounds_latvalue(nv, localInd);
2829          }
2830        }
2831      }
2832
2833      list_msgsLon.push_back(CMessage());
2834      list_msgsLat.push_back(CMessage());
2835
2836      list_msgsLon.back() << serverDomainId << hasLonLat;
2837      if (hasLonLat) 
2838        list_msgsLon.back() << list_lon.back();
2839      list_msgsLon.back()  << hasBounds;
2840      if (hasBounds)
2841      {
2842        list_msgsLon.back() << list_boundslon.back();
2843      }
2844
2845      list_msgsLat.back() << serverDomainId << hasLonLat;
2846      if (hasLonLat)
2847        list_msgsLat.back() << list_lat.back();
2848      list_msgsLat.back() << hasBounds;
2849      if (hasBounds)
2850      {         
2851        list_msgsLat.back() << list_boundslat.back();
2852      }
2853
2854      eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2855      eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
2856    }
2857    client->sendEvent(eventLon);
2858    client->sendEvent(eventLat);
2859  }
2860  CATCH_DUMP_ATTR
2861
2862  /*!
2863    Send data index to corresponding connected clients.
2864    Data index can be compressed however, we always send decompressed data index
2865    and they will be compressed on receiving.
2866    The compressed index are represented with 1 and others are represented with -1
2867  */
2868  void CDomain::sendDataIndex(CContextClient* client, const string& domainId)
2869  TRY
2870  {
2871    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2872   
2873    int ns, n, i, j, ind, nv, idx;
2874    int serverSize = client->serverSize;
2875
2876    // send area for each connected server
2877    CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2878
2879    list<CMessage> list_msgsDataIndex;
2880    list<CArray<int,1> > list_data_i_index, list_data_j_index;
2881
2882    int nbIndex = i_index.numElements();
2883    int niByIndex = max(i_index) - min(i_index) + 1;
2884    int njByIndex = max(j_index) - min(j_index) + 1; 
2885    int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2886    int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2887
2888   
2889    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2890    dataIIndex = -1; 
2891    dataJIndex = -1;
2892    ind = 0;
2893
2894    for (idx = 0; idx < data_i_index.numElements(); ++idx)
2895    {
2896      int dataIidx = data_i_index(idx) + data_ibegin;
2897      int dataJidx = data_j_index(idx) + data_jbegin;
2898      if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2899          (0 <= dataJidx) && (dataJidx < dataJindexBound))
2900      {
2901        dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2902        dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2903      }
2904    }
2905
2906    std::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2907    iteMap = indSrv_[serverSize].end();
2908    for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2909    {
2910      int nbData = 0;
2911      int rank = connectedServerRank_[serverSize][k];
2912      it = indSrv_[serverSize].find(rank);
2913      if (iteMap != it)
2914        nbData = it->second.size();
2915      list_data_i_index.push_back(CArray<int,1>(nbData));
2916      list_data_j_index.push_back(CArray<int,1>(nbData));
2917
2918      const std::vector<size_t>& temp = it->second;
2919      for (n = 0; n < nbData; ++n)
2920      {
2921        idx = static_cast<int>(it->second[n]);
2922        i = globalLocalIndexMap_[idx];
2923        list_data_i_index.back()(n) = dataIIndex(i);
2924        list_data_j_index.back()(n) = dataJIndex(i);
2925      }
2926
2927      list_msgsDataIndex.push_back(CMessage());
2928      list_msgsDataIndex.back() << serverDomainId ;
2929      list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2930      eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2931    }
2932    client->sendEvent(eventDataIndex);
2933  }
2934  CATCH
2935 
2936  bool CDomain::dispatchEvent(CEventServer& event)
2937  TRY
2938  {
2939    if (SuperClass::dispatchEvent(event)) return true;
2940    else
2941    {
2942      switch(event.type)
2943      {
2944        case EVENT_ID_SERVER_ATTRIBUT:
2945          recvDistributionAttributes(event);
2946          return true;
2947          break;
2948        case EVENT_ID_INDEX:
2949          recvIndex(event);
2950          return true;
2951          break;
2952        case EVENT_ID_LON:
2953          recvLon(event);
2954          return true;
2955          break;
2956        case EVENT_ID_LAT:
2957          recvLat(event);
2958          return true;
2959          break;
2960        case EVENT_ID_AREA:
2961          recvArea(event);
2962          return true;
2963          break; 
2964        case EVENT_ID_DATA_INDEX:
2965          recvDataIndex(event);
2966          return true;
2967          break;
2968        case EVENT_ID_DOMAIN_DISTRIBUTION:
2969          recvDomainDistribution(event);
2970          return true;
2971          break;
2972        case EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE:
2973          recvDistributedAttributes(event);
2974          return true;
2975          break; 
2976        default:
2977          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2978                << "Unknown Event");
2979          return false;
2980       }
2981    }
2982  }
2983  CATCH
2984
2985  /*!
2986    Receive index event from clients(s)
2987    \param[in] event event contain info about rank and associated index
2988  */
2989  void CDomain::recvIndex(CEventServer& event)
2990  TRY
2991  {
2992    string domainId;
2993    std::map<int, CBufferIn*> rankBuffers;
2994
2995    list<CEventServer::SSubEvent>::iterator it;
2996    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2997    {     
2998      CBufferIn* buffer = it->buffer;
2999      *buffer >> domainId;
3000      rankBuffers[it->rank] = buffer;       
3001    }
3002    get(domainId)->recvIndex(rankBuffers);
3003  }
3004  CATCH
3005
3006  /*!
3007    Receive index information from client(s). We use the global index for mapping index between
3008    sending clients and receiving clients.
3009    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
3010  */
3011  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
3012  TRY
3013  {
3014    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
3015    recvClientRanks_.resize(nbReceived);       
3016
3017    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
3018    ind = 0;
3019    for (ind = 0; it != ite; ++it, ++ind)
3020    {       
3021       recvClientRanks_[ind] = it->first;
3022       CBufferIn& buffer = *(it->second);
3023       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
3024       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
3025    }
3026    int nbIndGlob = 0;
3027    for (i = 0; i < nbReceived; ++i)
3028    {
3029      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
3030    }
3031   
3032    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
3033    i_index.resize(nbIndGlob);
3034    j_index.resize(nbIndGlob);   
3035    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
3036
3037    nbIndGlob = 0;
3038    for (i = 0; i < nbReceived; ++i)
3039    {
3040      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
3041      for (ind = 0; ind < tmp.numElements(); ++ind)
3042      {
3043         index = tmp(ind);
3044         if (0 == globalLocalIndexMap_.count(index))
3045         {
3046           iIndex = (index%ni_glo)-ibegin;
3047           iIndex = (iIndex < 0) ? 0 : iIndex;
3048           jIndex = (index/ni_glo)-jbegin;
3049           jIndex = (jIndex < 0) ? 0 : jIndex;
3050           nbIndLoc = iIndex + ni * jIndex;
3051           i_index(nbIndGlob) = index % ni_glo;
3052           j_index(nbIndGlob) = index / ni_glo;
3053           globalLocalIndexMap_[index] = nbIndGlob;
3054           ++nbIndGlob;
3055         } 
3056      } 
3057    } 
3058
3059    if (nbIndGlob==0)
3060    {
3061      i_index.resize(nbIndGlob);
3062      j_index.resize(nbIndGlob);
3063    }
3064    else
3065    {
3066      i_index.resizeAndPreserve(nbIndGlob);
3067      j_index.resizeAndPreserve(nbIndGlob);
3068    }
3069
3070    domainMask.resize(0); // Mask is not defined anymore on servers
3071  }
3072  CATCH
3073
3074  /*!
3075    Receive attributes event from clients(s)
3076    \param[in] event event contain info about rank and associated attributes
3077  */
3078  void CDomain::recvDistributionAttributes(CEventServer& event)
3079  TRY
3080  {
3081    CBufferIn* buffer=event.subEvents.begin()->buffer;
3082    string domainId ;
3083    *buffer>>domainId ;
3084    get(domainId)->recvDistributionAttributes(*buffer);
3085  }
3086  CATCH
3087
3088  /*!
3089    Receive attributes from client(s)
3090    \param[in] rank rank of client source
3091    \param[in] buffer message containing attributes info
3092  */
3093  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
3094  TRY
3095  {
3096    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
3097    int ni_glo_tmp, nj_glo_tmp;
3098    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
3099           >> ni_glo_tmp >> nj_glo_tmp
3100           >> isCompressible_;
3101
3102    ni.setValue(ni_tmp);
3103    ibegin.setValue(ibegin_tmp);
3104    nj.setValue(nj_tmp);
3105    jbegin.setValue(jbegin_tmp);
3106    ni_glo.setValue(ni_glo_tmp);
3107    nj_glo.setValue(nj_glo_tmp);
3108
3109  }
3110 CATCH_DUMP_ATTR
3111  /*!
3112    Receive longitude event from clients(s)
3113    \param[in] event event contain info about rank and associated longitude
3114  */
3115  void CDomain::recvLon(CEventServer& event)
3116  TRY
3117  {
3118    string domainId;
3119    std::map<int, CBufferIn*> rankBuffers;
3120
3121    list<CEventServer::SSubEvent>::iterator it;
3122    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3123    {     
3124      CBufferIn* buffer = it->buffer;
3125      *buffer >> domainId;
3126      rankBuffers[it->rank] = buffer;       
3127    }
3128    get(domainId)->recvLon(rankBuffers);
3129  }
3130  CATCH
3131
3132  /*!
3133    Receive longitude information from client(s)
3134    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
3135  */
3136  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
3137  TRY
3138  {
3139    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
3140    if (nbReceived != recvClientRanks_.size())
3141      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
3142           << "The number of sending clients is not correct."
3143           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3144   
3145    int nbLonInd = 0;
3146    vector<CArray<double,1> > recvLonValue(nbReceived);
3147    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
3148    for (i = 0; i < recvClientRanks_.size(); ++i)
3149    {
3150      int rank = recvClientRanks_[i];
3151      CBufferIn& buffer = *(rankBuffers[rank]);
3152      buffer >> hasLonLat;
3153      if (hasLonLat)
3154        buffer >> recvLonValue[i];
3155      buffer >> hasBounds;
3156      if (hasBounds)
3157        buffer >> recvBoundsLonValue[i];
3158    }
3159
3160    if (hasLonLat)
3161    {
3162      for (i = 0; i < nbReceived; ++i)
3163      {
3164        nbLonInd += recvLonValue[i].numElements();
3165      }
3166   
3167      if (nbLonInd != globalLocalIndexMap_.size())
3168        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3169                 << "something must be wrong with longitude index "<< std::endl;
3170
3171      nbLonInd = globalLocalIndexMap_.size();
3172      lonvalue.resize(nbLonInd);
3173      if (hasBounds)
3174      {
3175        bounds_lonvalue.resize(nvertex,nbLonInd);
3176        bounds_lonvalue = 0.;
3177      }
3178
3179      nbLonInd = 0;
3180      for (i = 0; i < nbReceived; ++i)
3181      {
3182        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3183        CArray<double,1>& tmp = recvLonValue[i];
3184        for (ind = 0; ind < tmp.numElements(); ++ind)
3185        {
3186          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3187          lonvalue(lInd) = tmp(ind); 
3188           if (hasBounds)
3189           {         
3190            for (int nv = 0; nv < nvertex; ++nv)
3191              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
3192           }                 
3193        }
3194      }       
3195    }
3196
3197   // setup attribute depending the type of domain
3198    if (hasLonLat)
3199    {
3200      nbLonInd = globalLocalIndexMap_.size();
3201      if (ni*nj != nbLonInd)
3202        ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
3203             << "The number of index received is not coherent with the given resolution"
3204             << "nbLonInd=" << nbLonInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
3205
3206      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
3207      {
3208        lonvalue_2d.resize(ni,nj);
3209          for(int ij=0, j=0 ; j<nj ; j++)
3210            for(int i=0 ; i<ni; i++, ij++) lonvalue_2d(i,j) = lonvalue(ij) ;
3211       
3212        if (hasBounds)
3213        {
3214          bounds_lon_2d.resize(nvertex, ni, nj) ;
3215          for(int ij=0, j=0 ; j<nj ; j++)
3216            for(int i=0 ; i<ni; i++, ij++) 
3217              for(int nv=0; nv<nvertex; nv++) bounds_lon_2d(nv,i,j) = bounds_lonvalue(nv,ij) ;
3218        }
3219      }
3220      else if (type == type_attr::unstructured || type == type_attr::gaussian)
3221      {
3222        lonvalue_1d.resize(nbLonInd);
3223        lonvalue_1d = lonvalue ;
3224        if (hasBounds)
3225        {
3226          bounds_lon_1d.resize(nvertex, nbLonInd) ;
3227          bounds_lon_1d = bounds_lonvalue ;
3228        }
3229      }
3230    }
3231  }
3232  CATCH_DUMP_ATTR
3233
3234  /*!
3235    Receive latitude event from clients(s)
3236    \param[in] event event contain info about rank and associated latitude
3237  */
3238  void CDomain::recvLat(CEventServer& event)
3239  TRY
3240  {
3241    string domainId;
3242    std::map<int, CBufferIn*> rankBuffers;
3243
3244    list<CEventServer::SSubEvent>::iterator it;
3245    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3246    {     
3247      CBufferIn* buffer = it->buffer;
3248      *buffer >> domainId;
3249      rankBuffers[it->rank] = buffer;   
3250    }
3251    get(domainId)->recvLat(rankBuffers);
3252  }
3253  CATCH
3254
3255  /*!
3256    Receive latitude information from client(s)
3257    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
3258  */
3259  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
3260  TRY
3261  {
3262    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
3263    if (nbReceived != recvClientRanks_.size())
3264      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
3265           << "The number of sending clients is not correct."
3266           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3267
3268    vector<CArray<double,1> > recvLatValue(nbReceived);
3269    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
3270    for (i = 0; i < recvClientRanks_.size(); ++i)
3271    {
3272      int rank = recvClientRanks_[i];
3273      CBufferIn& buffer = *(rankBuffers[rank]);
3274      buffer >> hasLonLat;
3275      if (hasLonLat)
3276        buffer >> recvLatValue[i];
3277      buffer >> hasBounds;
3278      if (hasBounds)
3279        buffer >> recvBoundsLatValue[i];
3280    }
3281
3282    int nbLatInd = 0;
3283    if (hasLonLat)
3284    {
3285      for (i = 0; i < nbReceived; ++i)
3286      {
3287        nbLatInd += recvLatValue[i].numElements();
3288      }
3289   
3290      if (nbLatInd != globalLocalIndexMap_.size())
3291        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3292                << "something must be wrong with latitude index "<< std::endl;
3293
3294      nbLatInd = globalLocalIndexMap_.size();
3295      latvalue.resize(nbLatInd);
3296      if (hasBounds)
3297      {
3298        bounds_latvalue.resize(nvertex,nbLatInd);
3299        bounds_latvalue = 0. ;
3300      }
3301
3302      nbLatInd = 0;
3303      for (i = 0; i < nbReceived; ++i)
3304      {
3305        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3306        CArray<double,1>& tmp = recvLatValue[i];
3307        for (ind = 0; ind < tmp.numElements(); ++ind)
3308        {
3309          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3310          latvalue(lInd) = tmp(ind);   
3311           if (hasBounds)
3312           {
3313            CArray<double,2>& boundslat = recvBoundsLatValue[i];
3314            for (int nv = 0; nv < nvertex; ++nv)
3315              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
3316           }   
3317          ++nbLatInd;
3318        }
3319      }       
3320    }
3321      // setup attribute depending the type of domain
3322    if (hasLonLat)
3323    {
3324      nbLatInd = globalLocalIndexMap_.size();
3325     
3326      if (ni*nj != nbLatInd)
3327        ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
3328             << "The number of index received is not coherent with the given resolution"
3329            << "nbLonInd=" << nbLatInd << ", ni=" << ni <<", nj"<<ni<<" ni*nj="<<ni*nj);
3330
3331      if (type == type_attr::rectilinear || type == type_attr::curvilinear)
3332      {
3333        {
3334          latvalue_2d.resize(ni,nj);
3335          for(int ij=0, j=0 ; j<nj ; j++)
3336            for(int i=0 ; i<ni; i++, ij++) latvalue_2d(i,j) = latvalue(ij) ;
3337       
3338          if (hasBounds)
3339          {
3340            bounds_lat_2d.resize(nvertex, ni, nj) ;
3341            for(int ij=0, j=0 ; j<nj ; j++)
3342              for(int i=0 ; i<ni; i++, ij++) 
3343                for(int nv=0; nv<nvertex; nv++) bounds_lat_2d(nv,i,j) = bounds_latvalue(nv,ij) ;
3344          }
3345        }
3346      }
3347      else if (type == type_attr::unstructured || type == type_attr::gaussian)
3348      {
3349        if (hasLonLat)
3350        {
3351          latvalue_1d.resize(nbLatInd);
3352          latvalue_1d = latvalue ;
3353          if (hasBounds)
3354          {
3355            bounds_lat_1d.resize(nvertex, nbLatInd) ;
3356            bounds_lat_1d = bounds_latvalue ;
3357          }
3358        }
3359      }
3360    } 
3361  }
3362  CATCH_DUMP_ATTR
3363
3364  /*!
3365    Receive area event from clients(s)
3366    \param[in] event event contain info about rank and associated area
3367  */
3368  void CDomain::recvArea(CEventServer& event)
3369  TRY
3370  {
3371    string domainId;
3372    std::map<int, CBufferIn*> rankBuffers;
3373
3374    list<CEventServer::SSubEvent>::iterator it;
3375    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3376    {     
3377      CBufferIn* buffer = it->buffer;
3378      *buffer >> domainId;
3379      rankBuffers[it->rank] = buffer;     
3380    }
3381    get(domainId)->recvArea(rankBuffers);
3382  }
3383  CATCH
3384
3385  /*!
3386    Receive area information from client(s)
3387    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3388  */
3389  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
3390  TRY
3391  {
3392    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
3393    if (nbReceived != recvClientRanks_.size())
3394      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
3395           << "The number of sending clients is not correct."
3396           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3397
3398    vector<CArray<double,1> > recvAreaValue(nbReceived);     
3399    for (i = 0; i < recvClientRanks_.size(); ++i)
3400    {
3401      int rank = recvClientRanks_[i];
3402      CBufferIn& buffer = *(rankBuffers[rank]);     
3403      buffer >> hasArea;
3404      if (hasArea)
3405        buffer >> recvAreaValue[i];
3406    }
3407
3408    if (hasArea)
3409    {
3410      int nbAreaInd = 0;
3411      for (i = 0; i < nbReceived; ++i)
3412      {     
3413        nbAreaInd += recvAreaValue[i].numElements();
3414      }
3415
3416      if (nbAreaInd != globalLocalIndexMap_.size())
3417        info (0) << "If domain " << this->getDomainOutputName() <<" does not have overlapped regions between processes "
3418                 << "something must be wrong with area index "<< std::endl;
3419
3420      nbAreaInd = globalLocalIndexMap_.size();
3421      areavalue.resize(nbAreaInd);
3422      nbAreaInd = 0;     
3423      for (i = 0; i < nbReceived; ++i)
3424      {
3425        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3426        CArray<double,1>& tmp = recvAreaValue[i];
3427        for (ind = 0; ind < tmp.numElements(); ++ind)
3428        {
3429          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3430          areavalue(lInd) = tmp(ind);         
3431        }
3432      }
3433     
3434      // fill the area attribute
3435      area.resize(ni,nj);
3436      for (int j = 0; j < nj; ++j)
3437      {
3438        for (int i = 0; i < ni; ++i)
3439        {
3440          int k = j * ni + i;
3441          area(i,j) = areavalue(k) ;
3442        }
3443      }
3444    }
3445  }
3446  CATCH_DUMP_ATTR
3447
3448  /*!
3449    Compare two domain objects.
3450    They are equal if only if they have identical attributes as well as their values.
3451    Moreover, they must have the same transformations.
3452  \param [in] domain Compared domain
3453  \return result of the comparison
3454  */
3455  bool CDomain::isEqual(CDomain* obj)
3456  TRY
3457  {
3458    vector<StdString> excludedAttr;
3459    excludedAttr.push_back("domain_ref");
3460    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
3461    if (!objEqual) return objEqual;
3462
3463    TransMapTypes thisTrans = this->getAllTransformations();
3464    TransMapTypes objTrans  = obj->getAllTransformations();
3465
3466    TransMapTypes::const_iterator it, itb, ite;
3467    std::vector<ETranformationType> thisTransType, objTransType;
3468    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
3469      thisTransType.push_back(it->first);
3470    for (it = objTrans.begin(); it != objTrans.end(); ++it)
3471      objTransType.push_back(it->first);
3472
3473    if (thisTransType.size() != objTransType.size()) return false;
3474    for (int idx = 0; idx < thisTransType.size(); ++idx)
3475      objEqual &= (thisTransType[idx] == objTransType[idx]);
3476
3477    return objEqual;
3478  }
3479  CATCH_DUMP_ATTR
3480
3481  /*!
3482    Receive data index event from clients(s)
3483    \param[in] event event contain info about rank and associated index
3484  */
3485  void CDomain::recvDataIndex(CEventServer& event)
3486  TRY
3487  {
3488    string domainId;
3489    std::map<int, CBufferIn*> rankBuffers;
3490
3491    list<CEventServer::SSubEvent>::iterator it;
3492    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
3493    {     
3494      CBufferIn* buffer = it->buffer;
3495      *buffer >> domainId;
3496      rankBuffers[it->rank] = buffer;       
3497    }
3498    get(domainId)->recvDataIndex(rankBuffers);
3499  }
3500  CATCH
3501
3502  /*!
3503    Receive data index information from client(s)
3504    A client receives data index from different clients to rebuild its own data index.
3505    Because we use global index + mask info to calculate the sending data to client(s),
3506    this data index must be updated with mask info (maybe it will change in the future)
3507    Because the data index is local, to rebuild data index of received client, we should use global index along with.
3508
3509    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
3510  */
3511  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
3512  TRY
3513  {
3514    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
3515    if (nbReceived != recvClientRanks_.size())
3516      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3517           << "The number of sending clients is not correct."
3518           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
3519
3520    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3521    for (i = 0; i < recvClientRanks_.size(); ++i)
3522    {
3523      int rank = recvClientRanks_[i];
3524      CBufferIn& buffer = *(rankBuffers[rank]);
3525      buffer >> recvDataIIndex[i];
3526      buffer >> recvDataJIndex[i];
3527    }
3528   
3529    int nbIndex = i_index.numElements();
3530    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3531    dataIIndex = -1; dataJIndex = -1;
3532     
3533    nbIndex = 0;
3534    for (i = 0; i < nbReceived; ++i)
3535    {     
3536      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3537      CArray<int,1>& tmpI = recvDataIIndex[i];   
3538      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3539      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3540          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3541             << "The number of global received index is not coherent with the number of received data index."
3542             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3543
3544      for (ind = 0; ind < tmpI.numElements(); ++ind)
3545      {
3546         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3547         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3548         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3549      } 
3550    }
3551
3552    int nbCompressedData = 0; 
3553    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3554    {
3555       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3556       if ((0 <= indexI) && (0 <= indexJ))
3557         ++nbCompressedData;
3558    }       
3559 
3560    data_i_index.resize(nbCompressedData);
3561    data_j_index.resize(nbCompressedData);
3562
3563    nbCompressedData = 0; 
3564    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3565    {
3566       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3567       if ((0 <= indexI) && (0 <= indexJ))
3568       {
3569          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3570          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3571         ++nbCompressedData;
3572       }
3573    }
3574
3575    // Reset data_ibegin, data_jbegin
3576    data_ibegin.setValue(0);
3577    data_jbegin.setValue(0);
3578  }
3579  CATCH_DUMP_ATTR
3580
3581  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3582  TRY
3583  {
3584    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3585    return transformationMap_.back().second;
3586  }
3587  CATCH_DUMP_ATTR
3588
3589  /*!
3590    Check whether a domain has transformation
3591    \return true if domain has transformation
3592  */
3593  bool CDomain::hasTransformation()
3594  TRY
3595  {
3596    return (!transformationMap_.empty());
3597  }
3598  CATCH_DUMP_ATTR
3599
3600  /*!
3601    Set transformation for current domain. It's the method to move transformation in hierarchy
3602    \param [in] domTrans transformation on domain
3603  */
3604  void CDomain::setTransformations(const TransMapTypes& domTrans)
3605  TRY
3606  {
3607    transformationMap_ = domTrans;
3608  }
3609  CATCH_DUMP_ATTR
3610
3611  /*!
3612    Get all transformation current domain has
3613    \return all transformation
3614  */
3615  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3616  TRY
3617  {
3618    return transformationMap_;
3619  }
3620  CATCH_DUMP_ATTR
3621
3622  void CDomain::duplicateTransformation(CDomain* src)
3623  TRY
3624  {
3625    if (src->hasTransformation())
3626    {
3627      this->setTransformations(src->getAllTransformations());
3628    }
3629  }
3630  CATCH_DUMP_ATTR
3631   
3632  /*!
3633   * Go through the hierarchy to find the domain from which the transformations must be inherited
3634   */
3635  void CDomain::solveInheritanceTransformation()
3636  TRY
3637  {
3638    if (hasTransformation() || !hasDirectDomainReference())
3639      return;
3640
3641    CDomain* domain = this;
3642    std::vector<CDomain*> refDomains;
3643    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
3644    {
3645      refDomains.push_back(domain);
3646      domain = domain->getDirectDomainReference();
3647    }
3648
3649    if (domain->hasTransformation())
3650      for (size_t i = 0; i < refDomains.size(); ++i)
3651        refDomains[i]->setTransformations(domain->getAllTransformations());
3652  }
3653  CATCH_DUMP_ATTR
3654
3655  void CDomain::setContextClient(CContextClient* contextClient)
3656  TRY
3657  {
3658    if (clientsSet.find(contextClient)==clientsSet.end())
3659    {
3660      clients.push_back(contextClient) ;
3661      clientsSet.insert(contextClient);
3662    }
3663  }
3664  CATCH_DUMP_ATTR
3665
3666  /*!
3667    Parse children nodes of a domain in xml file.
3668    Whenver there is a new transformation, its type and name should be added into this function
3669    \param node child node to process
3670  */
3671  void CDomain::parse(xml::CXMLNode & node)
3672  TRY
3673  {
3674    SuperClass::parse(node);
3675
3676    if (node.goToChildElement())
3677    {
3678      StdString nodeElementName;
3679      do
3680      {
3681        StdString nodeId("");
3682        if (node.getAttributes().end() != node.getAttributes().find("id"))
3683        { nodeId = node.getAttributes()["id"]; }
3684
3685        nodeElementName = node.getElementName();
3686        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
3687        it = transformationMapList_.find(nodeElementName);
3688        if (ite != it)
3689        {
3690          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3691                                                                                                                nodeId,
3692                                                                                                                &node)));
3693        }
3694        else
3695        {
3696          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3697                << "The transformation " << nodeElementName << " has not been supported yet.");
3698        }
3699      } while (node.goToNextElement()) ;
3700      node.goToParentElement();
3701    }
3702  }
3703  CATCH_DUMP_ATTR
3704   //----------------------------------------------------------------
3705
3706   DEFINE_REF_FUNC(Domain,domain)
3707
3708   ///---------------------------------------------------------------
3709
3710} // namespace xios
Note: See TracBrowser for help on using the repository browser.