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

Last change on this file since 2397 was 2397, checked in by ymipsl, 22 months ago
  • Optimize remote connector computation in case of read (reverse way).
  • don't compute anymore clientFromServerConnector (and all intermediate computation) for non reading case.

YM

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