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

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

Add a BANDS distribution in CDomain::computeRemoteElement regarding the number of servers and the number of clients.

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