source: XIOS3/branches/xios-3.0-beta/src/node/domain.cpp @ 2423

Last change on this file since 2423 was 2422, checked in by jderouillat, 20 months ago

Fix a division by 0 in BANDS domain distribution on servers

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