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

Last change on this file since 2613 was 2613, checked in by jderouillat, 4 months ago

Fix the attached mode for scalar output, and some bugs revealed by the adastra porting in debug mode

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