source: XIOS3/trunk/src/node/domain.cpp @ 2389

Last change on this file since 2389 was 2389, checked in by jderouillat, 2 years ago

Move renames of the attributes name in the element classes

  • 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: 96.3 KB
Line 
1#include "domain.hpp"
2#include "attribute_template.hpp"
3#include "object_template.hpp"
4#include "group_template.hpp"
5
6#include "xios_spl.hpp"
7#include "event_client.hpp"
8#include "event_server.hpp"
9#include "buffer_in.hpp"
10#include "message.hpp"
11#include "type.hpp"
12#include "context.hpp"
13#include "context_client.hpp"
14#include "context_server.hpp"
15#include "array_new.hpp"
16#include "distribution_client.hpp"
17#include "server_distribution_description.hpp"
18#include "client_server_mapping_distributed.hpp"
19#include "local_connector.hpp"
20#include "grid_local_connector.hpp"
21#include "remote_connector.hpp"
22#include "gatherer_connector.hpp"
23#include "scatterer_connector.hpp"
24#include "grid_scatterer_connector.hpp"
25#include "grid_gatherer_connector.hpp"
26#include "transformation_path.hpp"
27#include "grid_transformation_factory_impl.hpp"
28
29#include <algorithm>
30#include <regex>
31
32
33namespace xios {
34
35   /// ////////////////////// Définitions ////////////////////// ///
36
37   CDomain::CDomain(void)
38      : CObjectTemplate<CDomain>(), CDomainAttributes()
39      , isChecked(false), relFiles(),  indSrv_(), connectedServerRank_()
40      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
41      , hasLonLat(false)
42      , isRedistributed_(false), hasPole(false)
43      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
44      , clients()
45   {
46   }
47
48   CDomain::CDomain(const StdString & id)
49      : CObjectTemplate<CDomain>(id), CDomainAttributes()
50      , isChecked(false), relFiles(), indSrv_(), connectedServerRank_() 
51      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
52      , hasLonLat(false)
53      , isRedistributed_(false), hasPole(false)
54      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
55      , clients()
56    {
57    }
58
59   CDomain::~CDomain(void)
60   {
61   }
62
63   void CDomain::releaseStaticAllocation(void)
64   {
65     transformationMapList_.clear() ;
66     CTransformation<CDomain>::unregisterAllTransformations() ;
67     CGridTransformationFactory<CDomain>::unregisterAllTransformations() ;
68   }
69   ///---------------------------------------------------------------
70
71   void CDomain::assignMesh(const StdString meshName, const int nvertex)
72   TRY
73   {
74     mesh = CMesh::getMesh(meshName, nvertex);
75   }
76   CATCH_DUMP_ATTR
77
78   CDomain* CDomain::createDomain()
79   TRY
80   {
81     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
82     return domain;
83   }
84   CATCH
85
86   CDomain* CDomain::get(const string& id, bool noError)
87   {
88     const regex r("::");
89     smatch m;
90     if (regex_search(id, m, r))
91     {
92        if (m.size()!=1) ERROR("CDomain* CDomain::get(string& id)", <<" id = "<<id<< "  -> bad format id, separator :: append more than one time");
93        string fieldId=m.prefix() ;
94        if (fieldId.empty()) ERROR("CDomain* CDomain::get(string& id)", <<" id = "<<id<< "  -> bad format id, field name is empty");
95        string suffix=m.suffix() ;
96        if (!CField::has(fieldId)) 
97          if (noError)  return nullptr ;
98          else ERROR("CDomain* CDomain::get(string& id)", <<" id = "<<id<< "  -> field Id : < "<<fieldId<<" > doesn't exist");
99        CField* field=CField::get(fieldId) ;
100        return field->getAssociatedDomain(suffix, noError) ;
101     }
102     else 
103     {
104       if (noError) if(!CObjectFactory::HasObject<CDomain>(id)) return nullptr ;
105       return CObjectFactory::GetObject<CDomain>(id).get();
106     }
107   }
108
109   bool CDomain::has(const string& id)
110   {
111     if (CDomain::get(id,true)==nullptr) return false ;
112     else return true ;
113   }
114   
115   CField* CDomain::getFieldFromId(const string& id)
116   {
117     const regex r("::");
118     smatch m;
119     if (regex_search(id, m, r))
120     {
121        if (m.size()!=1) ERROR("CField* CDomain::getFieldFromId(const string& id)", <<" id = "<<id<< "  -> bad format id, separator :: append more than one time");
122        string fieldId=m.prefix() ;
123        if (fieldId.empty()) ERROR("CField* CDomain::getFieldFromId(const string& id)", <<" id = "<<id<< "  -> bad format id, field name is empty");
124        string suffix=m.suffix() ;
125        CField* field=CField::get(fieldId) ;
126        return field ;
127     }
128     else return nullptr;
129   }
130
131   const std::set<StdString> & CDomain::getRelFiles(void) const
132   TRY
133   {
134      return (this->relFiles);
135   }
136   CATCH
137
138     //----------------------------------------------------------------
139
140   /*!
141    * Compute the minimum buffer size required to send the attributes to the server(s).
142    *
143    * \return A map associating the server rank with its minimum buffer size.
144    */
145   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
146   TRY
147   {
148
149     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
150
151     if (client->isServerLeader())
152     {
153       // size estimation for sendDistributionAttribut
154       size_t size = 11 * sizeof(size_t);
155
156       const std::list<int>& ranks = client->getRanksServerLeader();
157       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
158       {
159         if (size > attributesSizes[*itRank])
160           attributesSizes[*itRank] = size;
161       }
162     }
163
164     std::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->getRemoteSize()].end();
165     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
166     for (size_t k = 0; k < connectedServerRank_[client->getRemoteSize()].size(); ++k)
167     {
168       int rank = connectedServerRank_[client->getRemoteSize()][k];
169       std::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->getRemoteSize()].find(rank);
170       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
171
172       // size estimation for sendIndex (and sendArea which is always smaller or equal)
173       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
174
175       // size estimation for sendLonLat
176       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
177       if (hasBounds)
178         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
179
180       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
181       if (size > attributesSizes[rank])
182         attributesSizes[rank] = size;
183     }
184
185     return attributesSizes;
186   }
187   CATCH_DUMP_ATTR
188
189   //----------------------------------------------------------------
190
191   bool CDomain::isEmpty(void) const
192   TRY
193   {
194     return ((this->i_index.isEmpty()) || (0 == this->i_index.numElements()));
195   }
196   CATCH
197
198   //----------------------------------------------------------------
199
200   bool CDomain::IsWritten(const StdString & filename) const
201   TRY
202   {
203      return (this->relFiles.find(filename) != this->relFiles.end());
204   }
205   CATCH
206
207   bool CDomain::isWrittenCompressed(const StdString& filename) const
208   TRY
209   {
210      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
211   }
212   CATCH
213
214   //----------------------------------------------------------------
215
216   bool CDomain::isDistributed(void) const
217   TRY
218   {
219      bool distributed =  !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
220              (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
221      bool distributed_glo ;
222      distributed |= (1 == CContext::getCurrent()->intraCommSize_);
223
224      return distributed;
225   }
226   CATCH
227
228   //----------------------------------------------------------------
229
230   /*!
231    * Compute if the domain can be ouput in a compressed way.
232    * In this case the workflow view on server side must be the same
233    * than the full view for all context rank. The result is stored on
234    * internal isCompressible_ attribute.
235    */
236   void CDomain::computeIsCompressible(void)
237   TRY
238   {
239     // mesh is compressible contains some masked or indexed value, ie if full view is different of workflow view.
240     // But now assume that the size of the 2 view must be equal for everybody. True on server side
241     int isSameView = getLocalView(CElementView::FULL)->getSize() ==  getLocalView(CElementView::WORKFLOW)->getSize();
242     MPI_Allreduce(MPI_IN_PLACE, &isSameView, 1, MPI_INT, MPI_LAND, CContext::getCurrent()->getIntraComm()) ;
243     if (isSameView) isCompressible_ = false ;
244     else isCompressible_ = true ;
245     isCompressibleComputed_=true ;
246   }
247   CATCH
248
249   void CDomain::addRelFile(const StdString & filename)
250   TRY
251   {
252      this->relFiles.insert(filename);
253   }
254   CATCH_DUMP_ATTR
255
256   void CDomain::addRelFileCompressed(const StdString& filename)
257   TRY
258   {
259      this->relFilesCompressed.insert(filename);
260   }
261   CATCH_DUMP_ATTR
262
263   StdString CDomain::GetName(void)   { return (StdString("domain")); }
264   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
265   ENodeType CDomain::GetType(void)   { return (eDomain); }
266
267   //----------------------------------------------------------------
268
269   /*!
270      Verify if all distribution information of a domain are available
271      This checking verifies the definition of distribution attributes (ni, nj, ibegin, jbegin)
272   */
273   bool CDomain::distributionAttributesHaveValue() const
274   TRY
275   {
276      bool hasValues = true;
277
278      if (ni.isEmpty() && ibegin.isEmpty() && i_index.isEmpty())
279      {
280        hasValues = false;
281        return hasValues;
282      }
283
284      return hasValues;
285   }
286   CATCH
287
288   /*!
289     Redistribute RECTILINEAR or CURVILINEAR domain with a number of local domains.
290   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
291   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
292   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
293    \param [in] nbLocalDomain number of local domain on the domain destination
294   */
295   void CDomain::redistribute(int nbLocalDomain)
296   TRY
297   {
298     if (this->isRedistributed_) return;
299
300     this->isRedistributed_ = true;
301     CContext* context = CContext::getCurrent();
302     // For now the assumption is that secondary server pools consist of the same number of procs.
303     // CHANGE the line below if the assumption changes.
304
305     int rankClient = context->intraCommRank_;
306     int rankOnDomain = rankClient%nbLocalDomain;
307
308     if (ni_glo.isEmpty() || ni_glo <= 0 )
309     {
310        ERROR("CDomain::redistribute(int nbLocalDomain)",
311           << "[ Id = " << this->getId() << " ] "
312           << "The global domain is badly defined,"
313           << " check the \'ni_glo\'  value !")
314     }
315
316     if (nj_glo.isEmpty() || nj_glo <= 0 )
317     {
318        ERROR("CDomain::redistribute(int nbLocalDomain)",
319           << "[ Id = " << this->getId() << " ] "
320           << "The global domain is badly defined,"
321           << " check the \'nj_glo\'  value !")
322     }
323
324     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
325     {
326        int globalDomainSize = ni_glo * nj_glo;
327        if (globalDomainSize <= nbLocalDomain)
328        {
329          for (int idx = 0; idx < nbLocalDomain; ++idx)
330          {
331            if (rankOnDomain < globalDomainSize)
332            {
333              int iIdx = rankOnDomain % ni_glo;
334              int jIdx = rankOnDomain / ni_glo;
335              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
336              ni.setValue(1); nj.setValue(1);
337            }
338            else
339            {
340              ibegin.setValue(0); jbegin.setValue(0);
341              ni.setValue(0); nj.setValue(0);
342            }
343          }
344        }
345        else
346        {
347          float njGlo = nj_glo.getValue();
348          float niGlo = ni_glo.getValue();
349          int nbProcOnX, nbProcOnY, range;
350
351          // Compute (approximately) number of segment on x and y axis
352          float yOverXRatio = njGlo/niGlo;
353
354          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
355          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
356
357          // Simple distribution: Sweep from top to bottom, left to right
358          // Calculate local begin on x
359          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
360          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
361          for (int i = 1; i < nbProcOnX; ++i)
362          {
363            range = ni_glo / nbProcOnX;
364            if (i < (ni_glo%nbProcOnX)) ++range;
365            niVec[i-1] = range;
366            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
367          }
368          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
369
370          // Calculate local begin on y
371          for (int j = 1; j < nbProcOnY; ++j)
372          {
373            range = nj_glo / nbProcOnY;
374            if (j < (nj_glo%nbProcOnY)) ++range;
375            njVec[j-1] = range;
376            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
377          }
378          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
379
380          // Now assign value to ni, ibegin, nj, jbegin
381          int iIdx = rankOnDomain % nbProcOnX;
382          int jIdx = rankOnDomain / nbProcOnX;
383
384          if (rankOnDomain != (nbLocalDomain-1))
385          {
386            ibegin.setValue(ibeginVec[iIdx]);
387            jbegin.setValue(jbeginVec[jIdx]);
388            nj.setValue(njVec[jIdx]);
389            ni.setValue(niVec[iIdx]);
390          }
391          else // just merge all the remaining rectangle into the last one
392          {
393            ibegin.setValue(ibeginVec[iIdx]);
394            jbegin.setValue(jbeginVec[jIdx]);
395            nj.setValue(njVec[jIdx]);
396            ni.setValue(ni_glo - ibeginVec[iIdx]);
397          }
398        } 
399     }
400     else  // unstructured domain
401     {
402       if (this->i_index.isEmpty())
403       {
404          int globalDomainSize = ni_glo * nj_glo;
405          if (globalDomainSize <= nbLocalDomain)
406          {
407            for (int idx = 0; idx < nbLocalDomain; ++idx)
408            {
409              if (rankOnDomain < globalDomainSize)
410              {
411                int iIdx = rankOnDomain % ni_glo;
412                int jIdx = rankOnDomain / ni_glo;
413                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
414                ni.setValue(1); nj.setValue(1);
415              }
416              else
417              {
418                ibegin.setValue(0); jbegin.setValue(0);
419                ni.setValue(0); nj.setValue(0);
420              }
421            }
422          }
423          else
424          {
425            float njGlo = nj_glo.getValue();
426            float niGlo = ni_glo.getValue();
427            std::vector<int> ibeginVec(nbLocalDomain,0);
428            std::vector<int> niVec(nbLocalDomain);
429            for (int i = 1; i < nbLocalDomain; ++i)
430            {
431              int range = ni_glo / nbLocalDomain;
432              if (i < (ni_glo%nbLocalDomain)) ++range;
433              niVec[i-1] = range;
434              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
435            }
436            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
437
438            int iIdx = rankOnDomain % nbLocalDomain;
439            ibegin.setValue(ibeginVec[iIdx]);
440            jbegin.setValue(0);
441            ni.setValue(niVec[iIdx]);
442            nj.setValue(1);
443          }
444
445          i_index.resize(ni);         
446          for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
447        }
448        else
449        {
450          ibegin.setValue(this->i_index(0));
451          jbegin.setValue(0);
452          ni.setValue(this->i_index.numElements());
453          nj.setValue(1);
454        }
455     }
456
457     checkDomain();
458   }
459   CATCH_DUMP_ATTR
460
461   /*!
462     Fill in longitude and latitude whose values are read from file
463   */
464   void CDomain::fillInLonLat()
465   TRY
466   {
467     switch (type)
468     {
469      case type_attr::rectilinear:
470        fillInRectilinearLonLat();
471        break;
472      case type_attr::curvilinear:
473        fillInCurvilinearLonLat();
474        break;
475      case type_attr::unstructured:
476        fillInUnstructuredLonLat();
477        break;
478
479      default:
480      break;
481     }
482     completeLonLatClient() ;
483   }
484   CATCH_DUMP_ATTR
485
486   /*!
487     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
488     Range of longitude value from 0 - 360
489     Range of latitude value from -90 - +90
490   */
491   void CDomain::fillInRectilinearLonLat()
492   TRY
493   {
494     if (!lonvalue_rectilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty() )
495     {
496       lonvalue_1d.resize(ni);
497       for (int idx = 0; idx < ni; ++idx)
498         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
499       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
500       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
501     }
502     else if (has_lon_in_read_file.isEmpty() || !has_lon_in_read_file)
503     {
504       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
505       lonvalue_1d.resize(ni);
506       double lonRange = lon_end - lon_start;
507       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
508
509        // Assign lon value
510       for (int i = 0; i < ni; ++i)
511       {
512         if (0 == (ibegin + i))
513         {
514           lonvalue_1d(i) = lon_start;
515         }
516         else if (ni_glo == (ibegin + i + 1))
517         {
518           lonvalue_1d(i) = lon_end;
519         }
520         else
521         {
522           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
523         }
524       }
525     }
526
527
528     if (!latvalue_rectilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
529     {
530       latvalue_1d.resize(nj);
531       for (int idx = 0; idx < nj; ++idx)
532         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
533       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
534       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
535     }
536     else if (has_lat_in_read_file.isEmpty() || !has_lat_in_read_file)
537     {
538       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
539       latvalue_1d.resize(nj);
540
541       double latRange = lat_end - lat_start;
542       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
543
544       for (int j = 0; j < nj; ++j)
545       {
546         if (0 == (jbegin + j))
547         {
548            latvalue_1d(j) = lat_start;
549         }
550         else if (nj_glo == (jbegin + j + 1))
551         {
552            latvalue_1d(j) = lat_end;
553         }
554         else
555         {
556           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
557         }
558       }
559     }
560   }
561   CATCH_DUMP_ATTR
562
563    /*
564      Fill in 2D longitude and latitude of curvilinear domain read from a file.
565      If there are already longitude and latitude defined by model. We just ignore read value.
566    */
567   void CDomain::fillInCurvilinearLonLat()
568   TRY
569   {
570     if (!lonvalue_curvilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty() && lonvalue_1d.isEmpty())
571     {
572       lonvalue_2d.resize(ni,nj);
573       for (int jdx = 0; jdx < nj; ++jdx)
574        for (int idx = 0; idx < ni; ++idx)
575         lonvalue_2d(idx,jdx) = lonvalue_curvilinear_read_from_file(idx, jdx);
576
577       lonvalue_curvilinear_read_from_file.free();
578     }
579
580     if (!latvalue_curvilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty() && latvalue_1d.isEmpty())
581     {
582       latvalue_2d.resize(ni,nj);
583       for (int jdx = 0; jdx < nj; ++jdx)
584        for (int idx = 0; idx < ni; ++idx)
585           latvalue_2d(idx,jdx) = latvalue_curvilinear_read_from_file(idx, jdx);
586
587       latvalue_curvilinear_read_from_file.free();
588     }
589
590     if (!bounds_lonvalue_curvilinear_read_from_file.isEmpty() && bounds_lon_2d.isEmpty() && bounds_lon_1d.isEmpty())
591     {
592       bounds_lon_2d.resize(nvertex,ni,nj);
593       for (int jdx = 0; jdx < nj; ++jdx)
594        for (int idx = 0; idx < ni; ++idx)
595          for (int ndx = 0; ndx < nvertex; ++ndx)
596            bounds_lon_2d(ndx,idx,jdx) = bounds_lonvalue_curvilinear_read_from_file(ndx,idx, jdx);
597
598       bounds_lonvalue_curvilinear_read_from_file.free();
599     }
600
601     if (!bounds_latvalue_curvilinear_read_from_file.isEmpty() && bounds_lat_2d.isEmpty() && bounds_lat_1d.isEmpty())
602     {
603       bounds_lat_2d.resize(nvertex,ni,nj);
604       for (int jdx = 0; jdx < nj; ++jdx)
605        for (int idx = 0; idx < ni; ++idx)
606          for (int ndx = 0; ndx < nvertex; ++ndx)
607            bounds_lat_2d(ndx,idx,jdx) = bounds_latvalue_curvilinear_read_from_file(ndx,idx, jdx);
608
609       bounds_latvalue_curvilinear_read_from_file.free();
610     }
611   }
612   CATCH_DUMP_ATTR
613
614    /*
615      Fill in longitude and latitude of unstructured domain read from a file
616      If there are already longitude and latitude defined by model. We just igonore reading value.
617    */
618   void CDomain::fillInUnstructuredLonLat()
619   TRY
620   {
621     if (i_index.isEmpty())
622     {
623       i_index.resize(ni);
624       for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
625     }
626
627     if (!lonvalue_unstructured_read_from_file.isEmpty() && lonvalue_1d.isEmpty())
628     {
629        lonvalue_1d.resize(ni);
630        for (int idx = 0; idx < ni; ++idx)
631          lonvalue_1d(idx) = lonvalue_unstructured_read_from_file(idx);
632
633        // We dont need these values anymore, so just delete them
634        lonvalue_unstructured_read_from_file.free();
635     } 
636
637     if (!latvalue_unstructured_read_from_file.isEmpty() && latvalue_1d.isEmpty())
638     {
639        latvalue_1d.resize(ni);
640        for (int idx = 0; idx < ni; ++idx)
641          latvalue_1d(idx) =  latvalue_unstructured_read_from_file(idx);
642
643        // We dont need these values anymore, so just delete them
644        latvalue_unstructured_read_from_file.free();
645     }
646
647     if (!bounds_lonvalue_unstructured_read_from_file.isEmpty() && bounds_lon_1d.isEmpty())
648     {
649        int nbVertex = nvertex;
650        bounds_lon_1d.resize(nbVertex,ni);
651        for (int idx = 0; idx < ni; ++idx)
652          for (int jdx = 0; jdx < nbVertex; ++jdx)
653            bounds_lon_1d(jdx,idx) = bounds_lonvalue_unstructured_read_from_file(jdx, idx);
654
655        // We dont need these values anymore, so just delete them
656        bounds_lonvalue_unstructured_read_from_file.free();
657     }
658
659     if (!bounds_latvalue_unstructured_read_from_file.isEmpty() && bounds_lat_1d.isEmpty())
660     {
661        int nbVertex = nvertex;
662        bounds_lat_1d.resize(nbVertex,ni);
663        for (int idx = 0; idx < ni; ++idx)
664          for (int jdx = 0; jdx < nbVertex; ++jdx)
665            bounds_lat_1d(jdx,idx) = bounds_latvalue_unstructured_read_from_file(jdx, idx);
666
667        // We dont need these values anymore, so just delete them
668        bounds_latvalue_unstructured_read_from_file.free();
669     }
670   }
671   CATCH_DUMP_ATTR
672
673  /*
674    Get global longitude and latitude of rectilinear domain.
675  */
676   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g)
677   TRY
678   {
679     CContext* context = CContext::getCurrent();
680    // For now the assumption is that secondary server pools consist of the same number of procs.
681    // CHANGE the line below if the assumption changes.
682     int clientSize = context->intraCommSize_ ;
683     lon_g.resize(ni_glo) ;
684     lat_g.resize(nj_glo) ;
685
686
687     int* ibegin_g = new int[clientSize] ;
688     int* jbegin_g = new int[clientSize] ;
689     int* ni_g = new int[clientSize] ;
690     int* nj_g = new int[clientSize] ;
691     int v ;
692     v=ibegin ;
693     MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,context->intraComm_) ;
694     v=jbegin ;
695     MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,context->intraComm_) ;
696     v=ni ;
697     MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,context->intraComm_) ;
698     v=nj ;
699     MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,context->intraComm_) ;
700
701     MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,context->intraComm_) ;
702     MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,context->intraComm_) ;
703
704      delete[] ibegin_g ;
705      delete[] jbegin_g ;
706      delete[] ni_g ;
707      delete[] nj_g ;
708   }
709   CATCH_DUMP_ATTR
710
711   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
712                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
713   TRY
714  {
715     int i,j,k;
716
717     const int nvertexValue = 4;
718     boundsLon.resize(nvertexValue,ni*nj);
719
720     if (ni_glo>1)
721     {
722       double lonStepStart = lon(1)-lon(0);
723       bounds_lon_start=lon(0) - lonStepStart/2;
724       double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2);
725       bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
726       double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
727
728       // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
729       if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
730       {
731         bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
732         bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
733       }
734     }
735     else
736     {
737       if (bounds_lon_start.isEmpty()) bounds_lon_start=-180. ;
738       if (bounds_lon_end.isEmpty()) bounds_lon_end=180.-1e-8 ;
739     }
740
741     for(j=0;j<nj;++j)
742       for(i=0;i<ni;++i)
743       {
744         k=j*ni+i;
745         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start
746                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2;
747         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end
748                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2;
749       }
750
751
752    boundsLat.resize(nvertexValue,nj*ni);
753    bool isNorthPole=false ;
754    bool isSouthPole=false ;
755    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
756    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
757
758    // lat boundaries beyond pole the assimilate it to pole
759    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
760    if (nj_glo>1)
761    {
762      double latStepStart = lat(1)-lat(0);
763      if (isNorthPole) bounds_lat_start=lat(0);
764      else
765      {
766        bounds_lat_start=lat(0)-latStepStart/2;
767        if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
768        else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
769        else if (bounds_lat_start <= 90 && bounds_lat_start >= lat(0))
770        {
771          if ( std::abs(90-bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
772        }
773        else if (bounds_lat_start >= -90 && bounds_lat_start <= lat(0))
774        {
775          if ( std::abs(-90 - bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
776        }
777      }
778
779      double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
780      if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
781      else
782      {
783        bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
784
785        if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
786        else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
787        else if (bounds_lat_end <= 90 && bounds_lat_end >= lat(nj_glo-1))
788        {
789          if ( std::abs(90-bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
790        }
791        else if (bounds_lat_end >= -90 && bounds_lat_end <= lat(nj_glo-1))
792        {
793          if ( std::abs(-90 - bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
794        }
795      }
796    }
797    else
798    {
799      if (bounds_lat_start.isEmpty()) bounds_lat_start=-90. ;
800      if (bounds_lat_end.isEmpty()) bounds_lat_end=90 ;
801    }
802
803    for(j=0;j<nj;++j)
804      for(i=0;i<ni;++i)
805      {
806        k=j*ni+i;
807        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start
808                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2;
809        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end
810                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2;
811      }
812   }
813   CATCH_DUMP_ATTR
814
815   /*
816     General check of the domain to verify its mandatory attributes
817   */
818   void CDomain::checkDomain(void)
819   TRY
820   {
821     if (type.isEmpty())
822     {
823       ERROR("CDomain::checkDomain(void)",
824             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
825             << "The domain type is mandatory, "
826             << "please define the 'type' attribute.")
827     }
828
829     if (type == type_attr::gaussian) 
830     {
831        hasPole=true ;
832        type.setValue(type_attr::unstructured) ;
833      }
834      else if (type == type_attr::rectilinear) hasPole=true ;
835
836     if (type == type_attr::unstructured)
837     {
838        if (ni_glo.isEmpty())
839        {
840          ERROR("CDomain::checkDomain(void)",
841                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
842                << "The global domain is badly defined, "
843                << "the mandatory 'ni_glo' attribute is missing.")
844        }
845        else if (ni_glo <= 0)
846        {
847          ERROR("CDomain::checkDomain(void)",
848                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
849                << "The global domain is badly defined, "
850                << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
851        }
852        isUnstructed_ = true;
853        nj_glo = 1;
854        nj = 1;
855        jbegin = 0;
856        if (!i_index.isEmpty()) 
857        {
858          ni = i_index.numElements();
859          j_index.resize(ni);
860          for(int i=0;i<ni;++i) j_index(i)=0;
861        }
862       
863
864        if (!area.isEmpty()) area.transposeSelf(1, 0); // => to be checked why is it transposed
865     }
866
867     if (ni_glo.isEmpty())
868     {
869       ERROR("CDomain::checkDomain(void)",
870             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
871             << "The global domain is badly defined, "
872             << "the mandatory 'ni_glo' attribute is missing.")
873     }
874     else if (ni_glo <= 0)
875     {
876       ERROR("CDomain::checkDomain(void)",
877             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
878             << "The global domain is badly defined, "
879             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
880     }
881
882     if (nj_glo.isEmpty())
883     {
884       ERROR("CDomain::checkDomain(void)",
885             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
886             << "The global domain is badly defined, "
887             << "the mandatory 'nj_glo' attribute is missing.")
888     }
889     else if (nj_glo <= 0)
890     {
891       ERROR("CDomain::checkDomain(void)",
892             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
893             << "The global domain is badly defined, "
894             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
895     }
896
897     checkLocalIDomain();
898     checkLocalJDomain();
899
900     if (i_index.isEmpty())
901     {
902       i_index.resize(ni*nj);
903       for (int j = 0; j < nj; ++j)
904         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
905     }
906
907     if (j_index.isEmpty())
908     {
909       j_index.resize(ni*nj);
910       for (int j = 0; j < nj; ++j)
911         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
912     }
913   }
914   CATCH_DUMP_ATTR
915
916   size_t CDomain::getGlobalWrittenSize(void)
917   {
918     return ni_glo*nj_glo ;
919   }
920   //----------------------------------------------------------------
921
922   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
923   void CDomain::checkLocalIDomain(void)
924   TRY
925   {
926      // If ibegin and ni are provided then we use them to check the validity of local domain
927      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
928      {
929        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
930        {
931          ERROR("CDomain::checkLocalIDomain(void)",
932                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
933                << "The local domain is wrongly defined,"
934                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
935        }
936      }
937
938      // i_index has higher priority than ibegin and ni
939      if (!i_index.isEmpty())
940      {
941        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
942        if (ni.isEmpty()) 
943        {         
944         // No information about ni
945          int minIndex = ni_glo - 1;
946          int maxIndex = 0;
947          for (int idx = 0; idx < i_index.numElements(); ++idx)
948          {
949            if (i_index(idx) < minIndex) minIndex = i_index(idx);
950            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
951          }
952          if (i_index.numElements()) {
953            ni = maxIndex - minIndex + 1; 
954            minIIndex = minIndex;
955          }         
956          else {
957            ni = 0;
958          }
959        }
960
961        // It's not so correct but if ibegin is not the first value of i_index
962        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
963        if (ibegin.isEmpty()) ibegin = minIIndex;
964      }
965      else if (ibegin.isEmpty() && ni.isEmpty())
966      {
967        ibegin = 0;
968        ni = ni_glo;
969      }
970      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
971      {
972        ERROR("CDomain::checkLocalIDomain(void)",
973              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
974              << "The local domain is wrongly defined," << endl
975              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
976              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
977      }
978       
979
980      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
981      {
982        ERROR("CDomain::checkLocalIDomain(void)",
983              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
984              << "The local domain is wrongly defined,"
985              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
986      }
987   }
988   CATCH_DUMP_ATTR
989
990   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
991   void CDomain::checkLocalJDomain(void)
992   TRY
993   {
994    // If jbegin and nj are provided then we use them to check the validity of local domain
995     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
996     {
997       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
998       {
999         ERROR("CDomain::checkLocalJDomain(void)",
1000                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1001                << "The local domain is wrongly defined,"
1002                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1003       }
1004     }
1005
1006     if (!j_index.isEmpty())
1007     {
1008        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
1009        if (nj.isEmpty()) 
1010        {
1011          // No information about nj
1012          int minIndex = nj_glo - 1;
1013          int maxIndex = 0;
1014          for (int idx = 0; idx < j_index.numElements(); ++idx)
1015          {
1016            if (j_index(idx) < minIndex) minIndex = j_index(idx);
1017            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
1018          }
1019          if (j_index.numElements()) {
1020            nj = maxIndex - minIndex + 1;
1021            minJIndex = minIndex; 
1022          }
1023          else
1024            nj = 0;
1025        } 
1026        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
1027        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
1028       if (jbegin.isEmpty()) jbegin = minJIndex;       
1029     }
1030     else if (jbegin.isEmpty() && nj.isEmpty())
1031     {
1032       jbegin = 0;
1033       nj = nj_glo;
1034     }     
1035
1036
1037     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
1038     {
1039       ERROR("CDomain::checkLocalJDomain(void)",
1040              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1041              << "The local domain is wrongly defined,"
1042              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
1043     }
1044   }
1045   CATCH_DUMP_ATTR
1046
1047   //----------------------------------------------------------------
1048
1049   void CDomain::checkMask(void)
1050   TRY
1051   {
1052      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
1053        ERROR("CDomain::checkMask(void)",
1054              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1055              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
1056              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
1057
1058      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
1059      {
1060        if (mask_1d.numElements() != i_index.numElements())
1061          ERROR("CDomain::checkMask(void)",
1062                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1063                << "'mask_1d' does not have the same size as the local domain." << std::endl
1064                << "Local size is " << i_index.numElements() << "." << std::endl
1065                << "Mask size is " << mask_1d.numElements() << ".");
1066      }
1067
1068      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
1069      {
1070        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1071          ERROR("CDomain::checkMask(void)",
1072                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1073                << "The mask does not have the same size as the local domain." << std::endl
1074                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1075                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
1076      }
1077
1078      if (!mask_2d.isEmpty())
1079      {
1080        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
1081        for (int j = 0; j < nj; ++j)
1082          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1083//        mask_2d.reset();
1084      }
1085      else if (mask_1d.isEmpty())
1086      {
1087        domainMask.resize(i_index.numElements());
1088        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
1089      }
1090      else
1091      {
1092      domainMask.resize(mask_1d.numElements());
1093      domainMask=mask_1d ;
1094     }
1095   }
1096   CATCH_DUMP_ATTR
1097
1098   //----------------------------------------------------------------
1099
1100   void CDomain::checkDomainData(void)
1101   TRY
1102   {
1103      if (data_dim.isEmpty())
1104      {
1105        data_dim.setValue(1);
1106      }
1107      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
1108      {
1109        ERROR("CDomain::checkDomainData(void)",
1110              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1111              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
1112      }
1113
1114      if (data_ibegin.isEmpty())
1115         data_ibegin.setValue(0);
1116      if (data_jbegin.isEmpty())
1117         data_jbegin.setValue(0);
1118
1119      if (data_ni.isEmpty())
1120      {
1121        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
1122      }
1123      else if (data_ni.getValue() < 0)
1124      {
1125        ERROR("CDomain::checkDomainData(void)",
1126              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1127              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
1128      }
1129
1130      if (data_nj.isEmpty())
1131      {
1132        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
1133      }
1134      else if (data_nj.getValue() < 0)
1135      {
1136        ERROR("CDomain::checkDomainData(void)",
1137              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1138              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1139      }
1140   }
1141   CATCH_DUMP_ATTR
1142
1143   //----------------------------------------------------------------
1144
1145   void CDomain::checkCompression(void)
1146   TRY
1147   {
1148     int i,j,ind;
1149      if (!data_i_index.isEmpty())
1150      {
1151        if (!data_j_index.isEmpty() &&
1152            data_j_index.numElements() != data_i_index.numElements())
1153        {
1154           ERROR("CDomain::checkCompression(void)",
1155                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1156                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1157                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1158                 << "'data_j_index' size = " << data_j_index.numElements());
1159        }
1160
1161        if (2 == data_dim)
1162        {
1163          if (data_j_index.isEmpty())
1164          {
1165             ERROR("CDomain::checkCompression(void)",
1166                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1167                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
1168          }
1169          for (int k=0; k<data_i_index.numElements(); ++k)
1170          {
1171            i = data_i_index(k)+data_ibegin ;
1172            j = data_j_index(k)+data_jbegin ;
1173            if (i>=0 && i<ni && j>=0 && j<nj)
1174            {
1175              ind=j*ni+i ;
1176              if ( (ind<0)||(!domainMask(ind)) )
1177              {
1178                data_i_index(k) = -1;
1179                data_j_index(k) = -1;
1180              }
1181            }
1182            else
1183            {
1184              data_i_index(k) = -1;
1185              data_j_index(k) = -1;
1186            }
1187          }
1188        }
1189        else // (1 == data_dim)
1190        {
1191          if (data_j_index.isEmpty())
1192          {
1193            data_j_index.resize(data_ni);
1194            data_j_index = 0;
1195          }
1196          for (int k=0; k<data_i_index.numElements(); ++k)
1197          {
1198            i=data_i_index(k)+data_ibegin ;
1199            if (i>=0 && i < domainMask.size())
1200            {
1201              if (!domainMask(i)) data_i_index(k) = -1;
1202            }
1203            else
1204              data_i_index(k) = -1;
1205
1206            if ( (i<0)||(!domainMask(i)) ) data_i_index(k) = -1;
1207          }
1208        }
1209      }
1210      else
1211      {
1212        if (data_dim == 2 && !data_j_index.isEmpty())
1213          ERROR("CDomain::checkCompression(void)",
1214                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1215                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
1216
1217        if (1 == data_dim)
1218        {
1219          data_i_index.resize(data_ni);
1220          data_j_index.resize(data_ni);
1221          data_j_index = 0;
1222
1223          for (int k = 0; k < data_ni; ++k)
1224          {
1225            i=k+data_ibegin ;
1226            if (i>=0 && i < domainMask.size())
1227            {
1228              if ( (i<0)||(!domainMask(i)) )
1229                data_i_index(k) = -1;
1230              else
1231                data_i_index(k) = k;
1232            }
1233            else
1234              data_i_index(k) = -1;
1235          }
1236        }
1237        else // (data_dim == 2)
1238        {
1239          const int dsize = data_ni * data_nj;
1240          data_i_index.resize(dsize);
1241          data_j_index.resize(dsize);
1242
1243          for(int count = 0, kj = 0; kj < data_nj; ++kj)
1244          {
1245            for(int ki = 0; ki < data_ni; ++ki, ++count)
1246            {
1247              i = ki + data_ibegin;
1248              j = kj + data_jbegin;
1249              ind=j*ni+i ;
1250              if (i>=0 && i<ni && j>=0 && j<nj)
1251              {
1252                if ( (ind<0)||(!domainMask(ind)) )
1253                {
1254                  data_i_index(count) = -1;
1255                  data_j_index(count) = -1;
1256                }
1257                else
1258                {
1259                  data_i_index(count) = ki;
1260                  data_j_index(count) = kj;
1261                }
1262              }
1263              else
1264              {
1265                data_i_index(count) = -1;
1266                data_j_index(count) = -1;
1267              }
1268            }
1269          }
1270        }
1271      }
1272   }
1273   CATCH_DUMP_ATTR
1274
1275   //----------------------------------------------------------------
1276   void CDomain::computeLocalMask(void)
1277   TRY
1278   {
1279     localMask.resize(i_index.numElements()) ;
1280     localMask=false ;
1281
1282     size_t dn=data_i_index.numElements() ;
1283     int i,j ;
1284     size_t k,ind ;
1285
1286     for(k=0;k<dn;k++)
1287     {
1288       if (data_dim==2)
1289       {
1290          i=data_i_index(k)+data_ibegin ;
1291          j=data_j_index(k)+data_jbegin ;
1292          if (i>=0 && i<ni && j>=0 && j<nj)
1293          {
1294            ind=j*ni+i ;
1295            localMask(ind)=domainMask(ind) ;
1296          }
1297       }
1298       else
1299       {
1300          i=data_i_index(k)+data_ibegin ;
1301          if (i>=0 && i<i_index.numElements())
1302          {
1303            ind=i ;
1304            localMask(ind)=domainMask(ind) ;
1305          }
1306       }
1307     }
1308   }
1309   CATCH_DUMP_ATTR
1310
1311
1312   //----------------------------------------------------------------
1313
1314   /*
1315     Fill in longitude, latitude, bounds, and area into internal values (lonvalue, latvalue, bounds_lonvalue, bounds_latvalue, areavalue)
1316     which will be used by XIOS.
1317   */
1318   void CDomain::completeLonLatClient(void)
1319   TRY
1320   {
1321     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1322     checkBounds() ;
1323     checkArea() ;
1324
1325     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
1326     {
1327       lonvalue.resize(ni * nj);
1328       latvalue.resize(ni * nj);
1329       if (hasBounds)
1330       {
1331         bounds_lonvalue.resize(nvertex, ni * nj);
1332         bounds_latvalue.resize(nvertex, ni * nj);
1333       }
1334
1335       for (int j = 0; j < nj; ++j)
1336       {
1337         for (int i = 0; i < ni; ++i)
1338         {
1339           int k = j * ni + i;
1340
1341           lonvalue(k) = lonvalue_2d(i,j);
1342           latvalue(k) = latvalue_2d(i,j);
1343
1344           if (hasBounds)
1345           {
1346             for (int n = 0; n < nvertex; ++n)
1347             {
1348               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1349               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
1350             }
1351           }
1352         }
1353       }
1354     }
1355     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
1356     {
1357       if (type_attr::rectilinear == type)
1358       {
1359         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1360         {
1361           lonvalue.resize(ni * nj);
1362           latvalue.resize(ni * nj);
1363           if (hasBounds)
1364           {
1365             bounds_lonvalue.resize(nvertex, ni * nj);
1366             bounds_latvalue.resize(nvertex, ni * nj);
1367           }
1368
1369           for (int j = 0; j < nj; ++j)
1370           {
1371             for (int i = 0; i < ni; ++i)
1372             {
1373               int k = j * ni + i;
1374
1375               lonvalue(k) = lonvalue_1d(i);
1376               latvalue(k) = latvalue_1d(j);
1377
1378               if (hasBounds)
1379               {
1380                 for (int n = 0; n < nvertex; ++n)
1381                 {
1382                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1383                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1384                 }
1385               }
1386             }
1387           }
1388         }
1389         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1390         {
1391           lonvalue.reference(lonvalue_1d);
1392           latvalue.reference(latvalue_1d);
1393            if (hasBounds)
1394           {
1395             bounds_lonvalue.reference(bounds_lon_1d);
1396             bounds_latvalue.reference(bounds_lat_1d);
1397           }
1398         }
1399         else
1400           ERROR("CDomain::completeLonClient(void)",
1401                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1402                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1403                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1404                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1405                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1406                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1407       }
1408       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1409       {
1410         lonvalue.reference(lonvalue_1d);
1411         latvalue.reference(latvalue_1d);
1412         if (hasBounds)
1413         {
1414           bounds_lonvalue.reference(bounds_lon_1d);
1415           bounds_latvalue.reference(bounds_lat_1d);
1416         }
1417       }
1418     }
1419
1420     if (!area.isEmpty() && areavalue.isEmpty())
1421     {
1422        areavalue.resize(ni*nj);
1423       for (int j = 0; j < nj; ++j)
1424       {
1425         for (int i = 0; i < ni; ++i)
1426         {
1427           int k = j * ni + i;
1428           areavalue(k) = area(i,j);
1429         }
1430       }
1431     }
1432   }
1433   CATCH_DUMP_ATTR
1434
1435   /*
1436     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1437   */
1438   void CDomain::convertLonLatValue(void)
1439   TRY
1440   {
1441     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1442     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1443     {
1444       lonvalue_2d.resize(ni,nj);
1445       latvalue_2d.resize(ni,nj);
1446       if (hasBounds)
1447       {
1448         bounds_lon_2d.resize(nvertex, ni, nj);
1449         bounds_lat_2d.resize(nvertex, ni, nj);
1450       }
1451
1452       for (int j = 0; j < nj; ++j)
1453       {
1454         for (int i = 0; i < ni; ++i)
1455         {
1456           int k = j * ni + i;
1457
1458           lonvalue_2d(i,j) = lonvalue(k);
1459           latvalue_2d(i,j) = latvalue(k);
1460
1461           if (hasBounds)
1462           {
1463             for (int n = 0; n < nvertex; ++n)
1464             {
1465               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1466               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1467             }
1468           }
1469         }
1470       }
1471     }
1472     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1473     {
1474       if (type_attr::rectilinear == type)
1475       {
1476         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1477         {
1478           lonvalue.resize(ni * nj);
1479           latvalue.resize(ni * nj);
1480           if (hasBounds)
1481           {
1482             bounds_lonvalue.resize(nvertex, ni * nj);
1483             bounds_latvalue.resize(nvertex, ni * nj);
1484           }
1485
1486           for (int j = 0; j < nj; ++j)
1487           {
1488             for (int i = 0; i < ni; ++i)
1489             {
1490               int k = j * ni + i;
1491
1492               lonvalue(k) = lonvalue_1d(i);
1493               latvalue(k) = latvalue_1d(j);
1494
1495               if (hasBounds)
1496               {
1497                 for (int n = 0; n < nvertex; ++n)
1498                 {
1499                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1500                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1501                 }
1502               }
1503             }
1504           }
1505         }
1506         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1507         {
1508           lonvalue.reference(lonvalue_1d);
1509           latvalue.reference(latvalue_1d);
1510            if (hasBounds)
1511           {
1512             bounds_lonvalue.reference(bounds_lon_1d);
1513             bounds_latvalue.reference(bounds_lat_1d);
1514           }
1515         }
1516         else
1517           ERROR("CDomain::completeLonClient(void)",
1518                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1519                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1520                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1521                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1522                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1523                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1524       }
1525       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1526       {
1527         lonvalue.reference(lonvalue_1d);
1528         latvalue.reference(latvalue_1d);
1529         if (hasBounds)
1530         {
1531           bounds_lonvalue.reference(bounds_lon_1d);
1532           bounds_latvalue.reference(bounds_lat_1d);
1533         }
1534       }
1535     }
1536   }
1537   CATCH_DUMP_ATTR
1538
1539   void CDomain::checkBounds(void)
1540   TRY
1541   {
1542     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1543     if (!nvertex.isEmpty() && nvertex > 0 && !hasBoundValues)
1544     {
1545       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1546         ERROR("CDomain::checkBounds(void)",
1547               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1548               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1549               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1550
1551       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1552         ERROR("CDomain::checkBounds(void)",
1553               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1554               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1555               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1556
1557       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1558       {
1559         ERROR("CDomain::checkBounds(void)",
1560               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1561               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1562               << "Please define either both attributes or none.");
1563       }
1564
1565       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1566       {
1567         ERROR("CDomain::checkBounds(void)",
1568               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1569               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1570               << "Please define either both attributes or none.");
1571       }
1572
1573       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1574         ERROR("CDomain::checkBounds(void)",
1575               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1576               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1577               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(0)
1578               << " but nvertex is " << nvertex.getValue() << ".");
1579
1580       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1581         ERROR("CDomain::checkBounds(void)",
1582               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1583               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1584               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(0)
1585               << " but nvertex is " << nvertex.getValue() << ".");
1586
1587       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1588         ERROR("CDomain::checkBounds(void)",
1589               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1590               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1591
1592       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1593         ERROR("CDomain::checkBounds(void)",
1594               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1595               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1596
1597       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1598         ERROR("CDomain::checkBounds(void)",
1599               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1600               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1601               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(0)
1602               << " but nvertex is " << nvertex.getValue() << ".");
1603
1604       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1605         ERROR("CDomain::checkBounds(void)",
1606               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1607               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1608               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(0)
1609               << " but nvertex is " << nvertex.getValue() << ".");
1610
1611       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1612         ERROR("CDomain::checkBounds(void)",
1613               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1614               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1615
1616       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1617         ERROR("CDomain::checkBounds(void)",
1618               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1619
1620       // In case of reading UGRID bounds values are not required
1621       hasBounds = (!bounds_lat_1d.isEmpty() || !bounds_lat_2d.isEmpty() );
1622     }
1623     else if (hasBoundValues)
1624     {
1625       hasBounds = true;       
1626     }
1627     else
1628     {
1629       hasBounds = false;
1630     }
1631   }
1632   CATCH_DUMP_ATTR
1633
1634   void CDomain::checkArea(void)
1635   TRY
1636   {
1637     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
1638     hasArea = !area.isEmpty();
1639     if (hasArea && !hasAreaValue)
1640     {
1641       if (area.extent(0) != ni || area.extent(1) != nj)
1642       {
1643         ERROR("CDomain::checkArea(void)",
1644               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1645               << "The area does not have the same size as the local domain." << std::endl
1646               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1647               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
1648       }
1649//       if (areavalue.isEmpty())
1650//       {
1651//          areavalue.resize(ni*nj);
1652//         for (int j = 0; j < nj; ++j)
1653//         {
1654//           for (int i = 0; i < ni; ++i)
1655//           {
1656//             int k = j * ni + i;
1657//             areavalue(k) = area(i,j);
1658//           }
1659//         }
1660//       }
1661     }
1662   }
1663   CATCH_DUMP_ATTR
1664
1665   void CDomain::checkLonLat()
1666   TRY
1667   {
1668     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1669                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1670     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1671     if (hasLonLat && !hasLonLatValue)
1672     {
1673       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1674         ERROR("CDomain::checkLonLat()",
1675               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1676               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1677               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1678
1679       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1680       {
1681         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1682           ERROR("CDomain::checkLonLat()",
1683                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1684                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1685                 << "Local size is " << i_index.numElements() << "." << std::endl
1686                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1687       }
1688
1689       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1690       {
1691         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1692           ERROR("CDomain::checkLonLat()",
1693                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1694                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1695                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1696                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1697       }
1698
1699       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1700         ERROR("CDomain::checkLonLat()",
1701               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1702               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1703               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1704
1705       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1706       {
1707         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1708           ERROR("CDomain::checkLonLat()",
1709                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1710                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1711                 << "Local size is " << i_index.numElements() << "." << std::endl
1712                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1713       }
1714
1715       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1716       {
1717         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1718           ERROR("CDomain::checkLonLat()",
1719                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1720                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1721                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1722                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1723       }
1724     }
1725   }
1726   CATCH_DUMP_ATTR
1727
1728   void CDomain::checkAttributes(void)
1729   TRY
1730   {
1731      if (this->checkAttributes_done_) return;
1732      this->checkDomain();
1733      this->checkLonLat();
1734      this->checkBounds();
1735      this->checkArea();
1736      this->checkMask();
1737      this->checkDomainData();
1738      this->checkCompression();
1739      this->computeLocalMask() ;
1740      this->completeLonLatClient();
1741      this->initializeLocalElement() ;
1742      this->addFullView() ; // probably do not automatically add View, but only if requested
1743      this->addWorkflowView() ; // probably do not automatically add View, but only if requested
1744      this->addModelView() ; // probably do not automatically add View, but only if requested
1745      // testing ?
1746     /*
1747      shared_ptr<CLocalView> local = localElement_->getView(CElementView::WORKFLOW) ;
1748      shared_ptr<CLocalView> model = localElement_->getView(CElementView::MODEL) ;
1749
1750      CLocalConnector test1(model, local) ;
1751      test1.computeConnector() ;
1752      CLocalConnector test2(local, model) ;
1753      test2.computeConnector() ;
1754      CGridLocalConnector gridTest1(vector<CLocalConnector*>{&test1}) ;
1755      CGridLocalConnector gridTest2(vector<CLocalConnector*>{&test2}) ;
1756     
1757     
1758      CArray<int,1> out1 ;
1759      CArray<int,1> out2 ;
1760      test1.transfer(data_i_index,out1,-111) ;
1761      test2.transfer(out1,out2,-111) ;
1762     
1763      out1 = 0 ;
1764      out2 = 0 ;
1765      gridTest1.transfer(data_i_index,out1,-111) ;
1766      gridTest2.transfer(out1, out2,-111) ;
1767    */ 
1768      this->checkAttributes_done_ = true;
1769   }
1770   CATCH_DUMP_ATTR
1771
1772   size_t CDomain::computeAttributesHash( MPI_Comm comm )
1773   {
1774     // Compute the hash of distributed attributs (value ...)
1775     int globalSize = this->ni_glo.getValue()*this->nj_glo.getValue();
1776     CArray<size_t,1> globalIndex; // No redundancy globalIndex will be computed with the connector
1777     shared_ptr<CGridTransformConnector> gridTransformConnector;
1778     // Compute a without redundancy element FULL view to enable a consistent hash computation
1779     this->getLocalView(CElementView::FULL)->createWithoutRedundancyFullViewConnector( globalSize, comm, gridTransformConnector, globalIndex );
1780     int localSize = globalIndex.numElements();
1781           
1782     CArray<double,1> lon_distributedValue, lat_distributedValue ;
1783     gridTransformConnector->transfer(this->lonvalue, lon_distributedValue );
1784     gridTransformConnector->transfer(this->latvalue, lat_distributedValue );
1785
1786     // Compute the distributed hash (v0) of the element
1787     // it will be associated to the default element name (= map key), and to the name really written
1788     size_t localHash = 0;
1789     for (int iloc=0; iloc<localSize ; iloc++ ) localHash+=globalIndex(iloc)*lon_distributedValue(iloc)*lat_distributedValue(iloc);
1790     size_t distributedHash = 0;
1791     MPI_Allreduce( &localHash, &distributedHash, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm  );
1792     
1793     // Compute the hash of global attributs (unit, prec ...)
1794     vector<StdString> excludedAttr;
1795     //excludedAttr.push_back("name");
1796     // internal attributs
1797     excludedAttr.insert(excludedAttr.end(), { "ibegin", "jbegin", "ni", "nj", "i_index", "j_index" });
1798     excludedAttr.insert(excludedAttr.end(), { "data_ni", "data_nj", "data_ibegin", "data_jbegin" });
1799     excludedAttr.insert(excludedAttr.end(), { "data_i_index", "data_j_index", "domain_ref" });
1800     // in distributed through lonvalue and latvalue
1801     excludedAttr.insert(excludedAttr.end(), { "lonvalue_1d", "latvalue_1d", "lonvalue_2d", "latvalue_2d" });
1802     // should be considered in distributed
1803     excludedAttr.insert(excludedAttr.end(), { "mask_1d", "mask_2d" }); // ???
1804     excludedAttr.insert(excludedAttr.end(), { "bounds_lon_1d", "bounds_lat_1d", "bounds_lon_2d", "bounds_lat_2d" });
1805     excludedAttr.insert(excludedAttr.end(), { "area" });
1806     // private
1807     excludedAttr.insert(excludedAttr.end(), { "lon_start", "lon_end", "lat_start", "lat_end" });
1808     excludedAttr.insert(excludedAttr.end(), { "bounds_lon_start", "bounds_lon_end", "bounds_lat_start", "bounds_lat_end" });
1809     excludedAttr.insert(excludedAttr.end(), { "has_lat_in_read_file", "has_lon_in_read_file", "has_bounds_lat_in_read_file", "has_bounds_lon_in_read_file" });
1810     excludedAttr.insert(excludedAttr.end(), { "lonvalue_rectilinear_read_from_file", "latvalue_rectilinear_read_from_file", "lonvalue_curvilinear_read_from_file", "latvalue_curvilinear_read_from_file" });
1811     excludedAttr.insert(excludedAttr.end(), { "bounds_lonvalue_curvilinear_read_from_file", "bounds_latvalue_curvilinear_read_from_file", "lonvalue_unstructured_read_from_file", "latvalue_unstructured_read_from_file" });
1812     excludedAttr.insert(excludedAttr.end(), { "bounds_lonvalue_unstructured_read_from_file", "bounds_latvalue_unstructured_read_from_file" });
1813     
1814     size_t globalHash = this->computeGlobalAttributesHash( excludedAttr );
1815
1816     return distributedHash + globalHash;
1817   }
1818 
1819   void CDomain::renameAttributesBeforeWriting(CDomain* writtenDomain)
1820   {
1821     if (writtenDomain!=NULL)
1822     {
1823       this->name = writtenDomain->getDomainOutputName();
1824       this->lon_name = writtenDomain->lon_name;
1825       this->lat_name = writtenDomain->lat_name;
1826       if (this->type == CDomain::type_attr::unstructured)
1827       {
1828         this->dim_i_name = writtenDomain->dim_i_name;
1829         this->dim_j_name = writtenDomain->dim_j_name;
1830       }
1831     }
1832     else
1833     {
1834       this->name =  this->getId();
1835       this->lon_name = "lon_"+this->getId();
1836       this->lat_name = "lat_"+this->getId();
1837       if (this->type == CDomain::type_attr::unstructured)
1838       {
1839         this->dim_i_name = "cell_"+this->getId();
1840         this->dim_j_name = "nvertex_"+this->getId();
1841       }
1842     }
1843   }
1844
1845   void CDomain::initializeLocalElement(void)
1846   {
1847      // after checkDomain i_index and j_index of size (ni*nj)
1848      int nij = ni*nj ;
1849      CArray<size_t, 1> ij_index(ni*nj) ;
1850      for(int ij=0; ij<nij ; ij++) ij_index(ij) = i_index(ij)+j_index(ij)*ni_glo ;
1851      int rank = CContext::getCurrent()->getIntraCommRank() ;
1852      localElement_ = make_shared<CLocalElement>(rank, ni_glo*nj_glo, ij_index) ;
1853   }
1854
1855   void CDomain::addFullView(void)
1856   {
1857      CArray<int,1> index(ni*nj) ;
1858      int nij=ni*nj ;
1859      for(int ij=0; ij<nij ; ij++) index(ij)=ij ;
1860      localElement_ -> addView(CElementView::FULL, index) ;
1861   }
1862
1863   void CDomain::addWorkflowView(void)
1864   {
1865     // information for workflow view is stored in localMask
1866     int nij=ni*nj ;
1867     int nMask=0 ;
1868     for(int ij=0; ij<nij ; ij++) if (localMask(ij)) nMask++ ;
1869     CArray<int,1> index(nMask) ;
1870
1871     nMask=0 ;
1872     for(int ij=0; ij<nij ; ij++) 
1873      if (localMask(ij))
1874      {
1875        index(nMask)=ij ;
1876        nMask++ ;
1877      }
1878      localElement_ -> addView(CElementView::WORKFLOW, index) ;
1879   }
1880
1881   void CDomain::addModelView(void)
1882   {
1883     // information for model view is stored in data_i_index/data_j_index
1884     // very weird, do not mix data_i_index and data_i_begin => in future only keep data_i_index
1885     int dataSize = data_i_index.numElements() ;
1886     CArray<int,1> index(dataSize) ;
1887     int i,j ;
1888     for(int k=0;k<dataSize;k++)
1889     {
1890        if (data_dim==2)
1891        {
1892          i=data_i_index(k)+data_ibegin ; // bad
1893          j=data_j_index(k)+data_jbegin ; // bad
1894          if (i>=0 && i<ni && j>=0 && j<nj) index(k)=i+j*ni ;
1895          else index(k)=-1 ;
1896        }
1897        else if (data_dim==1)
1898        {
1899          i=data_i_index(k)+data_ibegin ; // bad
1900          if (i>=0 && i<ni*nj) index(k)=i ;
1901          else index(k)=-1 ;
1902        }
1903     }
1904     localElement_->addView(CElementView::MODEL, index) ;
1905   }
1906       
1907   void CDomain::computeModelToWorkflowConnector(void)
1908   { 
1909     shared_ptr<CLocalView> srcView=getLocalView(CElementView::MODEL) ;
1910     shared_ptr<CLocalView> dstView=getLocalView(CElementView::WORKFLOW) ;
1911     modelToWorkflowConnector_ = make_shared<CLocalConnector>(srcView, dstView); 
1912     modelToWorkflowConnector_->computeConnector() ;
1913   }
1914
1915
1916  string CDomain::getCouplingAlias(const string& fieldId, int posInGrid)
1917  {
1918    return "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
1919  }
1920   
1921  /* to be removed later when coupling will be reimplemented, just to  not forget */
1922  void CDomain::sendDomainToCouplerOut(CContextClient* client, const string& fieldId, int posInGrid)
1923  {
1924    if (sendDomainToFileServer_done_.count(client)!=0) return ;
1925    else sendDomainToFileServer_done_.insert(client) ;
1926   
1927    const string domainId = getCouplingAlias(fieldId, posInGrid) ;
1928   
1929    if (!domain_ref.isEmpty())
1930    {
1931      auto domain_ref_tmp=domain_ref.getValue() ;
1932      domain_ref.reset() ; // remove the reference, find an other way to do that more cleanly
1933      this->sendAllAttributesToServer(client, domainId)  ;
1934      domain_ref = domain_ref_tmp ;
1935    }
1936    else this->sendAllAttributesToServer(client, domainId)  ;
1937  }
1938
1939
1940
1941
1942  void CDomain::makeAliasForCoupling(const string& fieldId, int posInGrid)
1943  {
1944    const string domainId = getCouplingAlias(fieldId, posInGrid);
1945    this->createAlias(domainId) ;
1946  }
1947
1948
1949  void CDomain::computeRemoteElement(CContextClient* client, EDistributionType distType)
1950  TRY
1951  {
1952    CContext* context = CContext::getCurrent();
1953    map<int, CArray<size_t,1>> globalIndex ;
1954/* old method
1955    if (type==EDistributionType::BANDS) // Bands distribution to send to file server
1956    {
1957      int nbServer = client->getRemoteSize();
1958      std::vector<int> nGlobDomain(2);
1959      nGlobDomain[0] = this->ni_glo;
1960      nGlobDomain[1] = this->nj_glo;
1961
1962      // to be changed in future, need to rewrite more simply domain distribution
1963      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
1964      int distributedPosition ;
1965      if (isUnstructed_) distributedPosition = 0 ;
1966      else distributedPosition = 1 ;
1967     
1968      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
1969      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
1970      vector<unordered_map<size_t,vector<int>>> indexServerOnElement ;
1971      CArray<int,1> axisDomainOrder(1) ; axisDomainOrder(0)=2 ;
1972      auto zeroIndex=serverDescription.computeServerGlobalByElement(indexServerOnElement, context->getIntraCommRank(), context->getIntraCommSize(),
1973                                                                  axisDomainOrder,distributedPosition) ;
1974      // distribution is very bad => to redo
1975      // convert indexServerOnElement => map<int,CArray<size_t,1>> - need to be changed later
1976      map<int, vector<size_t>> vectGlobalIndex ;
1977      for(auto& indexRanks : indexServerOnElement[0])
1978      {
1979        size_t index=indexRanks.first ;
1980        auto& ranks=indexRanks.second ;
1981        for(int rank : ranks) vectGlobalIndex[rank].push_back(index) ;
1982      }
1983      for(auto& vect : vectGlobalIndex ) globalIndex.emplace(vect.first, CArray<size_t,1>(vect.second.data(), shape(vect.second.size()),duplicateData)) ;
1984    // some servers receves no index (zeroIndex array) => root process take them into account.
1985      if (context->getIntraCommRank()==0)
1986        for(auto& rank : zeroIndex) globalIndex[rank] = CArray<size_t,1>() ;
1987    }
1988*/   
1989    if (distType==EDistributionType::BANDS && isUnstructed_) distType=EDistributionType::COLUMNS ;
1990
1991    if (distType==EDistributionType::BANDS) // Bands distribution to send to file server
1992    {
1993      int nbServer = client->getRemoteSize();
1994      int nbClient = client->getIntraCommSize() ;
1995      int rankClient = client->getIntraCommRank() ;
1996      int size = nbServer / nbClient ;
1997      int start ;
1998      if (nbServer%nbClient > rankClient)
1999      {
2000       start = (size+1) * rankClient ;
2001       size++ ;
2002      }
2003      else start = size*rankClient + nbServer%nbClient ;
2004     
2005      for(int i=0; i<size; i++)
2006      { 
2007        int rank=start+i ; 
2008        size_t indSize = nj_glo/nbServer ;
2009        size_t indStart ;
2010        if (nj_glo % nbServer > rank)
2011        {
2012          indStart = (indSize+1) * rank ;
2013          indSize++ ;
2014        }
2015        else indStart = indSize*rank + nj_glo%nbServer ;
2016       
2017        indStart=indStart*ni_glo ;
2018        indSize=indSize*ni_glo ;
2019        auto& globalInd =  globalIndex[rank] ;
2020        globalInd.resize(indSize) ;
2021        for(size_t n = 0 ; n<indSize; n++) globalInd(n)=indStart+n ;
2022      }
2023    }
2024    else if (distType==EDistributionType::COLUMNS) // Bands distribution to send to file server
2025    {
2026      int nbServer = client->getRemoteSize();
2027      int nbClient = client->getIntraCommSize() ;
2028      int rankClient = client->getIntraCommRank() ;
2029      int size = nbServer / nbClient ;
2030      int start ;
2031      if (nbServer%nbClient > rankClient)
2032      {
2033       start = (size+1) * rankClient ;
2034       size++ ;
2035      }
2036      else start = size*rankClient + nbServer%nbClient ;
2037     
2038      for(int i=0; i<size; i++)
2039      { 
2040        int rank=start+i ; 
2041        size_t indSize = ni_glo/nbServer ;
2042        size_t indStart ;
2043        if (ni_glo % nbServer > rank)
2044        {
2045          indStart = (indSize+1) * rank ;
2046          indSize++ ;
2047        }
2048        else indStart = indSize*rank + ni_glo%nbServer ;
2049       
2050        auto& globalInd =  globalIndex[rank] ;
2051        globalInd.resize(indSize*nj_glo) ;
2052        size_t n=0 ;
2053        for(int j=0; j<nj_glo;j++)
2054        {
2055          for(int i=0; i<indSize; i++, n++)  globalInd(n)=indStart+i ;
2056          indStart=indStart+ni_glo ; 
2057        }
2058      }
2059    }
2060    else if (distType==EDistributionType::NONE) // domain is not distributed ie all servers get the same local domain
2061    {
2062      int nbServer = client->getRemoteSize();
2063      int nglo=ni_glo*nj_glo ;
2064      CArray<size_t,1> indGlo ;
2065      for(size_t i=0;i<nglo;i++) indGlo(i) = i ;
2066      for (auto& rankServer : client->getRanksServerLeader()) globalIndex[rankServer] = indGlo ; 
2067    }
2068    remoteElement_[client] = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndex) ;
2069    remoteElement_[client]->addFullView() ;
2070  }
2071  CATCH
2072
2073 
2074
2075  void CDomain::distributeToServer(CContextClient* client, map<int, CArray<size_t,1>>& globalIndexOut, std::map<int, CArray<size_t,1>>& globalIndexIn,
2076                                   shared_ptr<CScattererConnector> &scattererConnector, const string& domainId)
2077  TRY
2078  {
2079    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2080    CContext* context = CContext::getCurrent();
2081
2082    this->sendAllAttributesToServer(client, serverDomainId)  ;
2083
2084    auto scatteredElement = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndexOut) ;
2085    scatteredElement->addFullView() ;
2086    scattererConnector = make_shared<CScattererConnector>(localElement_->getView(CElementView::FULL), scatteredElement->getView(CElementView::FULL), 
2087                                                          context->getIntraComm(), client->getRemoteSize()) ;
2088    scattererConnector->computeConnector() ;
2089
2090    // phase 0
2091    // send remote element to construct the full view on server, ie without hole
2092    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2093    CMessage message0 ;
2094    message0<<serverDomainId<<0 ; 
2095    remoteElement_[client]->sendToServer(client,event0,message0) ; 
2096   
2097    // phase 1
2098    // send the full view of element to construct the connector which connect distributed data coming from client to the full local view
2099    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2100    CMessage message1 ;
2101    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2102    scattererConnector->transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2103   
2104    sendDistributedAttributes(client, scattererConnector, domainId) ;
2105
2106 
2107    // phase 2 send the mask : data index + mask2D
2108    {
2109      CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2110      CArray<bool,1> maskOut ;
2111      auto workflowToFull = make_shared<CLocalConnector>(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2112      workflowToFull->computeConnector() ;
2113      maskIn=true ;
2114      workflowToFull->transfer(maskIn,maskOut,false) ;
2115
2116
2117      // prepare grid scatterer connector to send data from client to server
2118      map<int,CArray<size_t,1>> workflowGlobalIndex ;
2119      map<int,CArray<bool,1>> maskOut2 ; 
2120      scattererConnector->transfer(maskOut, maskOut2, false) ;
2121      scatteredElement->addView(CElementView::WORKFLOW, maskOut2) ;
2122      scatteredElement->getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2123      // create new workflow view for scattered element
2124      auto clientToServerElement = make_shared<CDistributedElement>(scatteredElement->getGlobalSize(), workflowGlobalIndex) ;
2125      clientToServerElement->addFullView() ;
2126      CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2127      CMessage message2 ;
2128      message2<<serverDomainId<<2 ; 
2129      clientToServerElement->sendToServer(client, event2, message2) ; 
2130      clientToServerConnector_[client] = make_shared<CScattererConnector>(localElement_->getView(CElementView::WORKFLOW), clientToServerElement->getView(CElementView::FULL),
2131                                                                        context->getIntraComm(), client->getRemoteSize()) ;
2132    clientToServerConnector_[client]->computeConnector() ;
2133    }
2134    ////////////
2135    // phase 3 : compute connector to receive from server
2136    ////////////
2137    {
2138      auto scatteredElement = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndexIn) ;
2139      scatteredElement->addFullView() ;
2140      auto scattererConnector = make_shared<CScattererConnector>(localElement_->getView(CElementView::FULL), scatteredElement->getView(CElementView::FULL), 
2141                                                                 context->getIntraComm(), client->getRemoteSize()) ;
2142      scattererConnector->computeConnector() ;
2143
2144      CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2145      CArray<bool,1> maskOut ;
2146      auto workflowToFull = make_shared<CLocalConnector>(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2147      workflowToFull->computeConnector() ;
2148      maskIn=true ;
2149      workflowToFull->transfer(maskIn,maskOut,false) ;
2150
2151      map<int,CArray<size_t,1>> workflowGlobalIndex ;
2152      map<int,CArray<bool,1>> maskOut2 ; 
2153      scattererConnector->transfer(maskOut, maskOut2, false) ;
2154      scatteredElement->addView(CElementView::WORKFLOW, maskOut2) ;
2155      scatteredElement->getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2156      auto clientToServerElement = make_shared<CDistributedElement>(scatteredElement->getGlobalSize(), workflowGlobalIndex) ;
2157      clientToServerElement->addFullView() ;
2158      CEventClient event3(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2159      CMessage message3 ;
2160      message3<<serverDomainId<<3 ; 
2161      clientToServerElement->sendToServer(client, event3, message3) ; 
2162
2163      clientFromServerConnector_[client] = make_shared<CGathererConnector>(clientToServerElement->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW));
2164      clientFromServerConnector_[client]->computeConnector() ;     
2165    }
2166  }
2167  CATCH
2168 
2169  void CDomain::recvDomainDistribution(CEventServer& event)
2170  TRY
2171  {
2172    string domainId;
2173    int phasis ;
2174    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
2175    get(domainId)->receivedDomainDistribution(event, phasis);
2176  }
2177  CATCH
2178
2179  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2180  TRY
2181  {
2182    CContext* context = CContext::getCurrent();
2183    if (phasis==0) // receive the remote element to construct the full view
2184    {
2185      localElement_ = make_shared<CLocalElement>(context->getIntraCommRank(),event) ;
2186      localElement_->addFullView() ;
2187      // construct the local dimension and indexes
2188      auto& globalIndex=localElement_->getGlobalIndex() ;
2189      int nij=globalIndex.numElements() ;
2190      int minI=ni_glo,maxI=-1,minJ=nj_glo,maxJ=-1 ;
2191      int i,j ;
2192      int niGlo=ni_glo, njGlo=njGlo ;
2193      for(int ij=0;ij<nij;ij++)
2194      {
2195        j=globalIndex(ij)/niGlo ;
2196        i=globalIndex(ij)%niGlo ;
2197        if (i<minI) minI=i ;
2198        if (i>maxI) maxI=i ;
2199        if (j<minJ) minJ=j ;
2200        if (j>maxJ) maxJ=j ;
2201      } 
2202      if (maxI>=minI) { ibegin=minI ; ni=maxI-minI+1 ; }
2203      else {ibegin=0; ni=0 ;}
2204      if (maxJ>=minJ) { jbegin=minJ ; nj=maxJ-minJ+1 ; }
2205      else {jbegin=0; nj=0 ;}
2206
2207    }
2208    else if (phasis==1) // receive the sent view from client to construct the full distributed full view on server
2209    {
2210      CContext* context = CContext::getCurrent();
2211      shared_ptr<CDistributedElement> elementFrom = make_shared<CDistributedElement>(event) ;
2212      elementFrom->addFullView() ;
2213      gathererConnector_ = make_shared<CGathererConnector>(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2214      gathererConnector_->computeConnector() ; 
2215    }
2216    else if (phasis==2)
2217    {
2218      elementFrom_ = make_shared<CDistributedElement>(event) ;
2219      elementFrom_->addFullView() ;
2220    }
2221    else if (phasis==3)
2222    {
2223      elementTo_ = make_shared<CDistributedElement>(event) ;
2224      elementTo_->addFullView() ;
2225    }
2226  }
2227  CATCH
2228
2229  void CDomain::setServerMask(CArray<bool,1>& serverMask, CContextClient* client)
2230  TRY
2231  {
2232    // nota : the client is needed to get the remote size for the scatterer connector. Maybe it is not the good place for this
2233    // Later, server to client connector can be computed on demand, with "client" as argument
2234    CContext* context = CContext::getCurrent();
2235    localElement_->addView(CElementView::WORKFLOW, serverMask) ;
2236    mask_1d.reference(serverMask.copy()) ;
2237 
2238    serverFromClientConnector_ = make_shared<CGathererConnector>(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW)) ;
2239    serverFromClientConnector_->computeConnector() ;
2240     
2241    serverToClientConnector_ = make_shared<CScattererConnector>(localElement_->getView(CElementView::WORKFLOW), elementTo_->getView(CElementView::FULL),
2242                                                                context->getIntraComm(), client->getRemoteSize()) ;
2243    serverToClientConnector_->computeConnector() ;
2244  }
2245  CATCH_DUMP_ATTR
2246
2247
2248  void CDomain::sendDistributedAttributes(CContextClient* client, shared_ptr<CScattererConnector> scattererConnector,  const string& domainId)
2249  {
2250    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2251    CContext* context = CContext::getCurrent();
2252
2253    if (hasLonLat)
2254    {
2255      { // send longitude
2256        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2257        CMessage message ;
2258        message<<serverDomainId<<string("lon") ; 
2259        scattererConnector->transfer(lonvalue, client, event,message) ;
2260      }
2261     
2262      { // send latitude
2263        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2264        CMessage message ;
2265        message<<serverDomainId<<string("lat") ; 
2266        scattererConnector->transfer(latvalue, client, event, message) ;
2267      }
2268    }
2269
2270    if (hasBounds)
2271    { 
2272      { // send longitude boudaries
2273        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2274        CMessage message ;
2275        message<<serverDomainId<<string("boundslon") ; 
2276        scattererConnector->transfer(nvertex, bounds_lonvalue, client, event, message ) ;
2277      }
2278
2279      { // send latitude boudaries
2280        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2281        CMessage message ;
2282        message<<serverDomainId<<string("boundslat") ; 
2283        scattererConnector->transfer(nvertex, bounds_latvalue, client, event, message ) ;
2284      }
2285    }
2286
2287    if (hasArea)
2288    {  // send area
2289      CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2290      CMessage message ;
2291      message<<serverDomainId<<string("area") ; 
2292      scattererConnector->transfer(areavalue, client, event,message) ;
2293    }
2294  }
2295
2296  void CDomain::recvDistributedAttributes(CEventServer& event)
2297  TRY
2298  {
2299    string domainId;
2300    string type ;
2301    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> type ;
2302    get(domainId)->recvDistributedAttributes(event, type);
2303  }
2304  CATCH
2305
2306  void CDomain::recvDistributedAttributes(CEventServer& event, const string& type)
2307  TRY
2308  {
2309    if (type=="lon") 
2310    {
2311      CArray<double,1> value ;
2312      gathererConnector_->transfer(event, value, 0.); 
2313      lonvalue_2d.resize(ni,nj) ;
2314      if (lonvalue_2d.numElements()>0) lonvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2315    }
2316    else if (type=="lat")
2317    {
2318      CArray<double,1> value ;
2319      gathererConnector_->transfer(event, value, 0.); 
2320      latvalue_2d.resize(ni,nj) ;
2321      if (latvalue_2d.numElements()>0) latvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2322    }
2323    else if (type=="boundslon")
2324    {
2325      CArray<double,1> value ;
2326      gathererConnector_->transfer(event, nvertex, value, 0.); 
2327      bounds_lon_2d.resize(nvertex,ni,nj) ;
2328      if (bounds_lon_2d.numElements()>0) bounds_lon_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2329    }
2330    else if (type=="boundslat")
2331    {
2332      CArray<double,1> value ;
2333      gathererConnector_->transfer(event, nvertex, value, 0.); 
2334      bounds_lat_2d.resize(nvertex,ni,nj) ;
2335      if (bounds_lat_2d.numElements()>0) bounds_lat_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2336    }
2337    else if (type=="area") 
2338    {
2339      CArray<double,1> value ;
2340      gathererConnector_->transfer(event, value, 0.); 
2341      area.resize(ni,nj) ;
2342      if (area.numElements()>0) area=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2343    }
2344  }
2345  CATCH
2346   
2347  bool CDomain::dispatchEvent(CEventServer& event)
2348  TRY
2349  {
2350    if (SuperClass::dispatchEvent(event)) return true;
2351    else
2352    {
2353      switch(event.type)
2354      {
2355        case EVENT_ID_DOMAIN_DISTRIBUTION:
2356          recvDomainDistribution(event);
2357          return true;
2358          break;
2359        case EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE:
2360          recvDistributedAttributes(event);
2361          return true;
2362          break; 
2363        default:
2364          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2365                << "Unknown Event");
2366          return false;
2367       }
2368    }
2369  }
2370  CATCH
2371
2372 
2373  /*!
2374    Compare two domain objects.
2375    They are equal if only if they have identical attributes as well as their values.
2376    Moreover, they must have the same transformations.
2377  \param [in] domain Compared domain
2378  \return result of the comparison
2379  */
2380  bool CDomain::isEqual(CDomain* obj)
2381  TRY
2382  {
2383    vector<StdString> excludedAttr;
2384    excludedAttr.push_back("domain_ref");
2385    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2386    if (!objEqual) return objEqual;
2387
2388    TransMapTypes thisTrans = this->getAllTransformations();
2389    TransMapTypes objTrans  = obj->getAllTransformations();
2390
2391    TransMapTypes::const_iterator it, itb, ite;
2392    std::vector<ETranformationType> thisTransType, objTransType;
2393    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2394      thisTransType.push_back(it->first);
2395    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2396      objTransType.push_back(it->first);
2397
2398    if (thisTransType.size() != objTransType.size()) return false;
2399    for (int idx = 0; idx < thisTransType.size(); ++idx)
2400      objEqual &= (thisTransType[idx] == objTransType[idx]);
2401
2402    return objEqual;
2403  }
2404  CATCH_DUMP_ATTR
2405
2406/////////////////////////////////////////////////////////////////////////
2407///////////////             TRANSFORMATIONS                    //////////
2408/////////////////////////////////////////////////////////////////////////
2409
2410  std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
2411  bool CDomain::dummyTransformationMapList_ = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
2412
2413  bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
2414  TRY
2415  {
2416    m["zoom_domain"] = TRANS_EXTRACT_DOMAIN;
2417    m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
2418    m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
2419    m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
2420    m["expand_domain"] = TRANS_EXPAND_DOMAIN;
2421    m["reorder_domain"] = TRANS_REORDER_DOMAIN;
2422    m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
2423    return true;
2424  }
2425  CATCH
2426
2427
2428  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2429  TRY
2430  {
2431    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2432    return transformationMap_.back().second;
2433  }
2434  CATCH_DUMP_ATTR
2435
2436  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, CTransformation<CDomain>* transformation)
2437  TRY
2438  {
2439    transformationMap_.push_back(std::make_pair(transType, transformation));
2440    return transformationMap_.back().second;
2441  }
2442  CATCH_DUMP_ATTR
2443  /*!
2444    Check whether a domain has transformation
2445    \return true if domain has transformation
2446  */
2447  bool CDomain::hasTransformation()
2448  TRY
2449  {
2450    return (!transformationMap_.empty());
2451  }
2452  CATCH_DUMP_ATTR
2453
2454  /*!
2455    Set transformation for current domain. It's the method to move transformation in hierarchy
2456    \param [in] domTrans transformation on domain
2457  */
2458  void CDomain::setTransformations(const TransMapTypes& domTrans)
2459  TRY
2460  {
2461    transformationMap_ = domTrans;
2462  }
2463  CATCH_DUMP_ATTR
2464
2465  /*!
2466    Get all transformation current domain has
2467    \return all transformation
2468  */
2469  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2470  TRY
2471  {
2472    return transformationMap_;
2473  }
2474  CATCH_DUMP_ATTR
2475
2476  void CDomain::duplicateTransformation(CDomain* src)
2477  TRY
2478  {
2479    if (src->hasTransformation())
2480    {
2481      this->setTransformations(src->getAllTransformations());
2482    }
2483  }
2484  CATCH_DUMP_ATTR
2485   
2486  /*!
2487   * Go through the hierarchy to find the domain from which the transformations must be inherited
2488   */
2489  void CDomain::solveInheritanceTransformation_old()
2490  TRY
2491  {
2492    if (hasTransformation() || !hasDirectDomainReference())
2493      return;
2494
2495    CDomain* domain = this;
2496    std::vector<CDomain*> refDomains;
2497    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2498    {
2499      refDomains.push_back(domain);
2500      domain = domain->getDirectDomainReference();
2501    }
2502
2503    if (domain->hasTransformation())
2504      for (size_t i = 0; i < refDomains.size(); ++i)
2505        refDomains[i]->setTransformations(domain->getAllTransformations());
2506  }
2507  CATCH_DUMP_ATTR
2508
2509
2510  void CDomain::solveInheritanceTransformation()
2511  TRY
2512  {
2513    if (solveInheritanceTransformation_done_) return;
2514    else solveInheritanceTransformation_done_=true ;
2515
2516    CDomain* domain = this;
2517    CDomain* Lastdomain ;
2518    std::list<CDomain*> refDomains;
2519    bool out=false ;
2520    vector<StdString> excludedAttr;
2521    excludedAttr.push_back("domain_ref");
2522   
2523    refDomains.push_front(domain) ;
2524    while (domain->hasDirectDomainReference() && !out)
2525    {
2526      CDomain* lastDomain=domain ;
2527      domain = domain->getDirectDomainReference();
2528      domain->solveRefInheritance() ;
2529      if (!domain->SuperClass::isEqual(lastDomain,excludedAttr)) out=true ;
2530      refDomains.push_front(domain) ;
2531    }
2532
2533    CTransformationPaths::TPath path ;
2534    auto& pathList = std::get<2>(path) ;
2535    std::get<0>(path) = EElement::DOMAIN ;
2536    std::get<1>(path) = refDomains.front()->getId() ;
2537    for (auto& domain : refDomains)
2538    {
2539      CDomain::TransMapTypes transformations = domain->getAllTransformations();
2540      for(auto& transformation : transformations) pathList.push_back({transformation.second->getTransformationType(), 
2541                                                                      transformation.second->getId()}) ;
2542    }
2543    transformationPaths_.addPath(path) ;
2544
2545  }
2546  CATCH_DUMP_ATTR
2547 
2548
2549  bool CDomain::activateFieldWorkflow(CGarbageCollector& gc)
2550  TRY
2551  {
2552    if (!domain_ref.isEmpty())
2553    {
2554      CField* field=getFieldFromId(domain_ref) ;
2555      if (field!=nullptr)
2556      {
2557        bool ret = field->buildWorkflowGraph(gc) ;
2558        if (!ret) return false ; // cannot build workflow graph at this state
2559      }
2560      else 
2561      {
2562        CDomain* domain = get(domain_ref) ;
2563        bool ret = domain->activateFieldWorkflow(gc) ;
2564        if (!ret) return false ; // cannot build workflow graph at this state
2565        domain_ref=domain->getId() ; // replace domain_ref by solved reference
2566      }
2567    }
2568    activateFieldWorkflow_done_=true ;
2569    return true ;
2570  }
2571  CATCH_DUMP_ATTR
2572/////////////////////////////////////////////////////////////////////////////////////////////
2573/////////////////////////////////////////////////////////////////////////////////////////////
2574
2575  void CDomain::setContextClient(CContextClient* contextClient)
2576  TRY
2577  {
2578    if (clientsSet.find(contextClient)==clientsSet.end())
2579    {
2580      clients.push_back(contextClient) ;
2581      clientsSet.insert(contextClient);
2582    }
2583  }
2584  CATCH_DUMP_ATTR
2585
2586  /*!
2587    Parse children nodes of a domain in xml file.
2588    Whenver there is a new transformation, its type and name should be added into this function
2589    \param node child node to process
2590  */
2591  void CDomain::parse(xml::CXMLNode & node)
2592  TRY
2593  {
2594    SuperClass::parse(node);
2595
2596    if (node.goToChildElement())
2597    {
2598      StdString nodeElementName;
2599      do
2600      {
2601        StdString nodeId("");
2602        if (node.getAttributes().end() != node.getAttributes().find("id"))
2603        { nodeId = node.getAttributes()["id"]; }
2604
2605        nodeElementName = node.getElementName();
2606        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2607        it = transformationMapList_.find(nodeElementName);
2608        if (ite != it)
2609        {
2610          if (it->first == "zoom_domain")
2611          {
2612            info(0) << "WARNING : " << it->first << " is deprecated, replaced by extract_domain." << endl;
2613          }
2614          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2615                                                                                                                nodeId,
2616                                                                                                                &node)));
2617        }
2618        else
2619        {
2620          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2621                << "The transformation " << nodeElementName << " has not been supported yet.");
2622        }
2623      } while (node.goToNextElement()) ;
2624      node.goToParentElement();
2625    }
2626  }
2627  CATCH_DUMP_ATTR
2628   //----------------------------------------------------------------
2629
2630   DEFINE_REF_FUNC(Domain,domain)
2631
2632   ///---------------------------------------------------------------
2633
2634} // namespace xios
Note: See TracBrowser for help on using the repository browser.