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

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

Set the code structure to compute the hash value of an element based on its attributs, use for now before writing an element in a file

  • 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: 95.5 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   int CDomain::computeAttributesHash( MPI_Comm comm )
1773   {
1774     int domain_hash = 0;
1775     
1776     // Compute the hash of distributed attributs (value ...)
1777     int globalSize = this->ni_glo.getValue()*this->nj_glo.getValue();
1778     CArray<size_t,1> globalIndex; // No redundancy globalIndex will be computed with the connector
1779     shared_ptr<CGridTransformConnector> gridTransformConnector;
1780     // Compute a without redundancy element FULL view to enable a consistent hash computation
1781     this->getLocalView(CElementView::FULL)->createWithoutRedundancyFullViewConnector( globalSize, comm, gridTransformConnector, globalIndex );
1782     int localSize = globalIndex.numElements();
1783           
1784     CArray<double,1> lon_distributedValue, lat_distributedValue ;
1785     gridTransformConnector->transfer(this->lonvalue, lon_distributedValue );
1786     gridTransformConnector->transfer(this->latvalue, lat_distributedValue );
1787
1788     // Compute the distributed hash (v0) of the element
1789     // it will be associated to the default element name (= map key), and to the name really written
1790     int localHash = 0;
1791     for (int iloc=0; iloc<localSize ; iloc++ ) localHash+=globalIndex(iloc)*lon_distributedValue(iloc)*lat_distributedValue(iloc);
1792     int distributedHash = 0;
1793     MPI_Allreduce( &localHash, &distributedHash, 1, MPI_INT, MPI_SUM, comm  );
1794     
1795     // Compute the hash of global attributs (unit, prec ...)
1796     vector<StdString> excludedAttr;
1797     //excludedAttr.push_back("name");
1798     // internal attributs
1799     excludedAttr.insert(excludedAttr.end(), { "ibegin", "jbegin", "ni", "nj", "i_index", "j_index" });
1800     excludedAttr.insert(excludedAttr.end(), { "data_ni", "data_nj", "data_ibegin", "data_jbegin" });
1801     excludedAttr.insert(excludedAttr.end(), { "data_i_index", "data_j_index", "domain_ref" });
1802     // in distributed through lonvalue and latvalue
1803     excludedAttr.insert(excludedAttr.end(), { "lonvalue_1d", "latvalue_1d", "lonvalue_2d", "latvalue_2d" });
1804     // should be considered in distributed
1805     excludedAttr.insert(excludedAttr.end(), { "mask_1d", "mask_2d" }); // ???
1806     excludedAttr.insert(excludedAttr.end(), { "bounds_lon_1d", "bounds_lat_1d", "bounds_lon_2d", "bounds_lat_2d" });
1807     excludedAttr.insert(excludedAttr.end(), { "area" });
1808     // private
1809     excludedAttr.insert(excludedAttr.end(), { "lon_start", "lon_end", "lat_start", "lat_end" });
1810     excludedAttr.insert(excludedAttr.end(), { "bounds_lon_start", "bounds_lon_end", "bounds_lat_start", "bounds_lat_end" });
1811     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" });
1812     excludedAttr.insert(excludedAttr.end(), { "lonvalue_rectilinear_read_from_file", "latvalue_rectilinear_read_from_file", "lonvalue_curvilinear_read_from_file", "latvalue_curvilinear_read_from_file" });
1813     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" });
1814     excludedAttr.insert(excludedAttr.end(), { "bounds_lonvalue_unstructured_read_from_file", "bounds_latvalue_unstructured_read_from_file" });
1815     
1816     int globalHash = this->computeGlobalAttributesHash( excludedAttr );
1817
1818     return distributedHash + globalHash;
1819   }
1820
1821   void CDomain::initializeLocalElement(void)
1822   {
1823      // after checkDomain i_index and j_index of size (ni*nj)
1824      int nij = ni*nj ;
1825      CArray<size_t, 1> ij_index(ni*nj) ;
1826      for(int ij=0; ij<nij ; ij++) ij_index(ij) = i_index(ij)+j_index(ij)*ni_glo ;
1827      int rank = CContext::getCurrent()->getIntraCommRank() ;
1828      localElement_ = make_shared<CLocalElement>(rank, ni_glo*nj_glo, ij_index) ;
1829   }
1830
1831   void CDomain::addFullView(void)
1832   {
1833      CArray<int,1> index(ni*nj) ;
1834      int nij=ni*nj ;
1835      for(int ij=0; ij<nij ; ij++) index(ij)=ij ;
1836      localElement_ -> addView(CElementView::FULL, index) ;
1837   }
1838
1839   void CDomain::addWorkflowView(void)
1840   {
1841     // information for workflow view is stored in localMask
1842     int nij=ni*nj ;
1843     int nMask=0 ;
1844     for(int ij=0; ij<nij ; ij++) if (localMask(ij)) nMask++ ;
1845     CArray<int,1> index(nMask) ;
1846
1847     nMask=0 ;
1848     for(int ij=0; ij<nij ; ij++) 
1849      if (localMask(ij))
1850      {
1851        index(nMask)=ij ;
1852        nMask++ ;
1853      }
1854      localElement_ -> addView(CElementView::WORKFLOW, index) ;
1855   }
1856
1857   void CDomain::addModelView(void)
1858   {
1859     // information for model view is stored in data_i_index/data_j_index
1860     // very weird, do not mix data_i_index and data_i_begin => in future only keep data_i_index
1861     int dataSize = data_i_index.numElements() ;
1862     CArray<int,1> index(dataSize) ;
1863     int i,j ;
1864     for(int k=0;k<dataSize;k++)
1865     {
1866        if (data_dim==2)
1867        {
1868          i=data_i_index(k)+data_ibegin ; // bad
1869          j=data_j_index(k)+data_jbegin ; // bad
1870          if (i>=0 && i<ni && j>=0 && j<nj) index(k)=i+j*ni ;
1871          else index(k)=-1 ;
1872        }
1873        else if (data_dim==1)
1874        {
1875          i=data_i_index(k)+data_ibegin ; // bad
1876          if (i>=0 && i<ni*nj) index(k)=i ;
1877          else index(k)=-1 ;
1878        }
1879     }
1880     localElement_->addView(CElementView::MODEL, index) ;
1881   }
1882       
1883   void CDomain::computeModelToWorkflowConnector(void)
1884   { 
1885     shared_ptr<CLocalView> srcView=getLocalView(CElementView::MODEL) ;
1886     shared_ptr<CLocalView> dstView=getLocalView(CElementView::WORKFLOW) ;
1887     modelToWorkflowConnector_ = make_shared<CLocalConnector>(srcView, dstView); 
1888     modelToWorkflowConnector_->computeConnector() ;
1889   }
1890
1891
1892  string CDomain::getCouplingAlias(const string& fieldId, int posInGrid)
1893  {
1894    return "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
1895  }
1896   
1897  /* to be removed later when coupling will be reimplemented, just to  not forget */
1898  void CDomain::sendDomainToCouplerOut(CContextClient* client, const string& fieldId, int posInGrid)
1899  {
1900    if (sendDomainToFileServer_done_.count(client)!=0) return ;
1901    else sendDomainToFileServer_done_.insert(client) ;
1902   
1903    const string domainId = getCouplingAlias(fieldId, posInGrid) ;
1904   
1905    if (!domain_ref.isEmpty())
1906    {
1907      auto domain_ref_tmp=domain_ref.getValue() ;
1908      domain_ref.reset() ; // remove the reference, find an other way to do that more cleanly
1909      this->sendAllAttributesToServer(client, domainId)  ;
1910      domain_ref = domain_ref_tmp ;
1911    }
1912    else this->sendAllAttributesToServer(client, domainId)  ;
1913  }
1914
1915
1916
1917
1918  void CDomain::makeAliasForCoupling(const string& fieldId, int posInGrid)
1919  {
1920    const string domainId = getCouplingAlias(fieldId, posInGrid);
1921    this->createAlias(domainId) ;
1922  }
1923
1924
1925  void CDomain::computeRemoteElement(CContextClient* client, EDistributionType distType)
1926  TRY
1927  {
1928    CContext* context = CContext::getCurrent();
1929    map<int, CArray<size_t,1>> globalIndex ;
1930/* old method
1931    if (type==EDistributionType::BANDS) // Bands distribution to send to file server
1932    {
1933      int nbServer = client->getRemoteSize();
1934      std::vector<int> nGlobDomain(2);
1935      nGlobDomain[0] = this->ni_glo;
1936      nGlobDomain[1] = this->nj_glo;
1937
1938      // to be changed in future, need to rewrite more simply domain distribution
1939      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
1940      int distributedPosition ;
1941      if (isUnstructed_) distributedPosition = 0 ;
1942      else distributedPosition = 1 ;
1943     
1944      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
1945      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
1946      vector<unordered_map<size_t,vector<int>>> indexServerOnElement ;
1947      CArray<int,1> axisDomainOrder(1) ; axisDomainOrder(0)=2 ;
1948      auto zeroIndex=serverDescription.computeServerGlobalByElement(indexServerOnElement, context->getIntraCommRank(), context->getIntraCommSize(),
1949                                                                  axisDomainOrder,distributedPosition) ;
1950      // distribution is very bad => to redo
1951      // convert indexServerOnElement => map<int,CArray<size_t,1>> - need to be changed later
1952      map<int, vector<size_t>> vectGlobalIndex ;
1953      for(auto& indexRanks : indexServerOnElement[0])
1954      {
1955        size_t index=indexRanks.first ;
1956        auto& ranks=indexRanks.second ;
1957        for(int rank : ranks) vectGlobalIndex[rank].push_back(index) ;
1958      }
1959      for(auto& vect : vectGlobalIndex ) globalIndex.emplace(vect.first, CArray<size_t,1>(vect.second.data(), shape(vect.second.size()),duplicateData)) ;
1960    // some servers receves no index (zeroIndex array) => root process take them into account.
1961      if (context->getIntraCommRank()==0)
1962        for(auto& rank : zeroIndex) globalIndex[rank] = CArray<size_t,1>() ;
1963    }
1964*/   
1965    if (distType==EDistributionType::BANDS && isUnstructed_) distType=EDistributionType::COLUMNS ;
1966
1967    if (distType==EDistributionType::BANDS) // Bands distribution to send to file server
1968    {
1969      int nbServer = client->getRemoteSize();
1970      int nbClient = client->getIntraCommSize() ;
1971      int rankClient = client->getIntraCommRank() ;
1972      int size = nbServer / nbClient ;
1973      int start ;
1974      if (nbServer%nbClient > rankClient)
1975      {
1976       start = (size+1) * rankClient ;
1977       size++ ;
1978      }
1979      else start = size*rankClient + nbServer%nbClient ;
1980     
1981      for(int i=0; i<size; i++)
1982      { 
1983        int rank=start+i ; 
1984        size_t indSize = nj_glo/nbServer ;
1985        size_t indStart ;
1986        if (nj_glo % nbServer > rank)
1987        {
1988          indStart = (indSize+1) * rank ;
1989          indSize++ ;
1990        }
1991        else indStart = indSize*rank + nj_glo%nbServer ;
1992       
1993        indStart=indStart*ni_glo ;
1994        indSize=indSize*ni_glo ;
1995        auto& globalInd =  globalIndex[rank] ;
1996        globalInd.resize(indSize) ;
1997        for(size_t n = 0 ; n<indSize; n++) globalInd(n)=indStart+n ;
1998      }
1999    }
2000    else if (distType==EDistributionType::COLUMNS) // Bands distribution to send to file server
2001    {
2002      int nbServer = client->getRemoteSize();
2003      int nbClient = client->getIntraCommSize() ;
2004      int rankClient = client->getIntraCommRank() ;
2005      int size = nbServer / nbClient ;
2006      int start ;
2007      if (nbServer%nbClient > rankClient)
2008      {
2009       start = (size+1) * rankClient ;
2010       size++ ;
2011      }
2012      else start = size*rankClient + nbServer%nbClient ;
2013     
2014      for(int i=0; i<size; i++)
2015      { 
2016        int rank=start+i ; 
2017        size_t indSize = ni_glo/nbServer ;
2018        size_t indStart ;
2019        if (ni_glo % nbServer > rank)
2020        {
2021          indStart = (indSize+1) * rank ;
2022          indSize++ ;
2023        }
2024        else indStart = indSize*rank + ni_glo%nbServer ;
2025       
2026        auto& globalInd =  globalIndex[rank] ;
2027        globalInd.resize(indSize*nj_glo) ;
2028        size_t n=0 ;
2029        for(int j=0; j<nj_glo;j++)
2030        {
2031          for(int i=0; i<indSize; i++, n++)  globalInd(n)=indStart+i ;
2032          indStart=indStart+ni_glo ; 
2033        }
2034      }
2035    }
2036    else if (distType==EDistributionType::NONE) // domain is not distributed ie all servers get the same local domain
2037    {
2038      int nbServer = client->getRemoteSize();
2039      int nglo=ni_glo*nj_glo ;
2040      CArray<size_t,1> indGlo ;
2041      for(size_t i=0;i<nglo;i++) indGlo(i) = i ;
2042      for (auto& rankServer : client->getRanksServerLeader()) globalIndex[rankServer] = indGlo ; 
2043    }
2044    remoteElement_[client] = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndex) ;
2045    remoteElement_[client]->addFullView() ;
2046  }
2047  CATCH
2048
2049 
2050
2051  void CDomain::distributeToServer(CContextClient* client, map<int, CArray<size_t,1>>& globalIndexOut, std::map<int, CArray<size_t,1>>& globalIndexIn,
2052                                   shared_ptr<CScattererConnector> &scattererConnector, const string& domainId)
2053  TRY
2054  {
2055    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2056    CContext* context = CContext::getCurrent();
2057
2058    this->sendAllAttributesToServer(client, serverDomainId)  ;
2059
2060    auto scatteredElement = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndexOut) ;
2061    scatteredElement->addFullView() ;
2062    scattererConnector = make_shared<CScattererConnector>(localElement_->getView(CElementView::FULL), scatteredElement->getView(CElementView::FULL), 
2063                                                          context->getIntraComm(), client->getRemoteSize()) ;
2064    scattererConnector->computeConnector() ;
2065
2066    // phase 0
2067    // send remote element to construct the full view on server, ie without hole
2068    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2069    CMessage message0 ;
2070    message0<<serverDomainId<<0 ; 
2071    remoteElement_[client]->sendToServer(client,event0,message0) ; 
2072   
2073    // phase 1
2074    // send the full view of element to construct the connector which connect distributed data coming from client to the full local view
2075    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2076    CMessage message1 ;
2077    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2078    scattererConnector->transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2079   
2080    sendDistributedAttributes(client, scattererConnector, domainId) ;
2081
2082 
2083    // phase 2 send the mask : data index + mask2D
2084    {
2085      CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2086      CArray<bool,1> maskOut ;
2087      auto workflowToFull = make_shared<CLocalConnector>(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2088      workflowToFull->computeConnector() ;
2089      maskIn=true ;
2090      workflowToFull->transfer(maskIn,maskOut,false) ;
2091
2092
2093      // prepare grid scatterer connector to send data from client to server
2094      map<int,CArray<size_t,1>> workflowGlobalIndex ;
2095      map<int,CArray<bool,1>> maskOut2 ; 
2096      scattererConnector->transfer(maskOut, maskOut2, false) ;
2097      scatteredElement->addView(CElementView::WORKFLOW, maskOut2) ;
2098      scatteredElement->getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2099      // create new workflow view for scattered element
2100      auto clientToServerElement = make_shared<CDistributedElement>(scatteredElement->getGlobalSize(), workflowGlobalIndex) ;
2101      clientToServerElement->addFullView() ;
2102      CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2103      CMessage message2 ;
2104      message2<<serverDomainId<<2 ; 
2105      clientToServerElement->sendToServer(client, event2, message2) ; 
2106      clientToServerConnector_[client] = make_shared<CScattererConnector>(localElement_->getView(CElementView::WORKFLOW), clientToServerElement->getView(CElementView::FULL),
2107                                                                        context->getIntraComm(), client->getRemoteSize()) ;
2108    clientToServerConnector_[client]->computeConnector() ;
2109    }
2110    ////////////
2111    // phase 3 : compute connector to receive from server
2112    ////////////
2113    {
2114      auto scatteredElement = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndexIn) ;
2115      scatteredElement->addFullView() ;
2116      auto scattererConnector = make_shared<CScattererConnector>(localElement_->getView(CElementView::FULL), scatteredElement->getView(CElementView::FULL), 
2117                                                                 context->getIntraComm(), client->getRemoteSize()) ;
2118      scattererConnector->computeConnector() ;
2119
2120      CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2121      CArray<bool,1> maskOut ;
2122      auto workflowToFull = make_shared<CLocalConnector>(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2123      workflowToFull->computeConnector() ;
2124      maskIn=true ;
2125      workflowToFull->transfer(maskIn,maskOut,false) ;
2126
2127      map<int,CArray<size_t,1>> workflowGlobalIndex ;
2128      map<int,CArray<bool,1>> maskOut2 ; 
2129      scattererConnector->transfer(maskOut, maskOut2, false) ;
2130      scatteredElement->addView(CElementView::WORKFLOW, maskOut2) ;
2131      scatteredElement->getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2132      auto clientToServerElement = make_shared<CDistributedElement>(scatteredElement->getGlobalSize(), workflowGlobalIndex) ;
2133      clientToServerElement->addFullView() ;
2134      CEventClient event3(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2135      CMessage message3 ;
2136      message3<<serverDomainId<<3 ; 
2137      clientToServerElement->sendToServer(client, event3, message3) ; 
2138
2139      clientFromServerConnector_[client] = make_shared<CGathererConnector>(clientToServerElement->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW));
2140      clientFromServerConnector_[client]->computeConnector() ;     
2141    }
2142  }
2143  CATCH
2144 
2145  void CDomain::recvDomainDistribution(CEventServer& event)
2146  TRY
2147  {
2148    string domainId;
2149    int phasis ;
2150    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
2151    get(domainId)->receivedDomainDistribution(event, phasis);
2152  }
2153  CATCH
2154
2155  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2156  TRY
2157  {
2158    CContext* context = CContext::getCurrent();
2159    if (phasis==0) // receive the remote element to construct the full view
2160    {
2161      localElement_ = make_shared<CLocalElement>(context->getIntraCommRank(),event) ;
2162      localElement_->addFullView() ;
2163      // construct the local dimension and indexes
2164      auto& globalIndex=localElement_->getGlobalIndex() ;
2165      int nij=globalIndex.numElements() ;
2166      int minI=ni_glo,maxI=-1,minJ=nj_glo,maxJ=-1 ;
2167      int i,j ;
2168      int niGlo=ni_glo, njGlo=njGlo ;
2169      for(int ij=0;ij<nij;ij++)
2170      {
2171        j=globalIndex(ij)/niGlo ;
2172        i=globalIndex(ij)%niGlo ;
2173        if (i<minI) minI=i ;
2174        if (i>maxI) maxI=i ;
2175        if (j<minJ) minJ=j ;
2176        if (j>maxJ) maxJ=j ;
2177      } 
2178      if (maxI>=minI) { ibegin=minI ; ni=maxI-minI+1 ; }
2179      else {ibegin=0; ni=0 ;}
2180      if (maxJ>=minJ) { jbegin=minJ ; nj=maxJ-minJ+1 ; }
2181      else {jbegin=0; nj=0 ;}
2182
2183    }
2184    else if (phasis==1) // receive the sent view from client to construct the full distributed full view on server
2185    {
2186      CContext* context = CContext::getCurrent();
2187      shared_ptr<CDistributedElement> elementFrom = make_shared<CDistributedElement>(event) ;
2188      elementFrom->addFullView() ;
2189      gathererConnector_ = make_shared<CGathererConnector>(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2190      gathererConnector_->computeConnector() ; 
2191    }
2192    else if (phasis==2)
2193    {
2194      elementFrom_ = make_shared<CDistributedElement>(event) ;
2195      elementFrom_->addFullView() ;
2196    }
2197    else if (phasis==3)
2198    {
2199      elementTo_ = make_shared<CDistributedElement>(event) ;
2200      elementTo_->addFullView() ;
2201    }
2202  }
2203  CATCH
2204
2205  void CDomain::setServerMask(CArray<bool,1>& serverMask, CContextClient* client)
2206  TRY
2207  {
2208    // nota : the client is needed to get the remote size for the scatterer connector. Maybe it is not the good place for this
2209    // Later, server to client connector can be computed on demand, with "client" as argument
2210    CContext* context = CContext::getCurrent();
2211    localElement_->addView(CElementView::WORKFLOW, serverMask) ;
2212    mask_1d.reference(serverMask.copy()) ;
2213 
2214    serverFromClientConnector_ = make_shared<CGathererConnector>(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW)) ;
2215    serverFromClientConnector_->computeConnector() ;
2216     
2217    serverToClientConnector_ = make_shared<CScattererConnector>(localElement_->getView(CElementView::WORKFLOW), elementTo_->getView(CElementView::FULL),
2218                                                                context->getIntraComm(), client->getRemoteSize()) ;
2219    serverToClientConnector_->computeConnector() ;
2220  }
2221  CATCH_DUMP_ATTR
2222
2223
2224  void CDomain::sendDistributedAttributes(CContextClient* client, shared_ptr<CScattererConnector> scattererConnector,  const string& domainId)
2225  {
2226    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2227    CContext* context = CContext::getCurrent();
2228
2229    if (hasLonLat)
2230    {
2231      { // send longitude
2232        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2233        CMessage message ;
2234        message<<serverDomainId<<string("lon") ; 
2235        scattererConnector->transfer(lonvalue, client, event,message) ;
2236      }
2237     
2238      { // send latitude
2239        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2240        CMessage message ;
2241        message<<serverDomainId<<string("lat") ; 
2242        scattererConnector->transfer(latvalue, client, event, message) ;
2243      }
2244    }
2245
2246    if (hasBounds)
2247    { 
2248      { // send longitude boudaries
2249        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2250        CMessage message ;
2251        message<<serverDomainId<<string("boundslon") ; 
2252        scattererConnector->transfer(nvertex, bounds_lonvalue, client, event, message ) ;
2253      }
2254
2255      { // send latitude boudaries
2256        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2257        CMessage message ;
2258        message<<serverDomainId<<string("boundslat") ; 
2259        scattererConnector->transfer(nvertex, bounds_latvalue, client, event, message ) ;
2260      }
2261    }
2262
2263    if (hasArea)
2264    {  // send area
2265      CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2266      CMessage message ;
2267      message<<serverDomainId<<string("area") ; 
2268      scattererConnector->transfer(areavalue, client, event,message) ;
2269    }
2270  }
2271
2272  void CDomain::recvDistributedAttributes(CEventServer& event)
2273  TRY
2274  {
2275    string domainId;
2276    string type ;
2277    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> type ;
2278    get(domainId)->recvDistributedAttributes(event, type);
2279  }
2280  CATCH
2281
2282  void CDomain::recvDistributedAttributes(CEventServer& event, const string& type)
2283  TRY
2284  {
2285    if (type=="lon") 
2286    {
2287      CArray<double,1> value ;
2288      gathererConnector_->transfer(event, value, 0.); 
2289      lonvalue_2d.resize(ni,nj) ;
2290      if (lonvalue_2d.numElements()>0) lonvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2291    }
2292    else if (type=="lat")
2293    {
2294      CArray<double,1> value ;
2295      gathererConnector_->transfer(event, value, 0.); 
2296      latvalue_2d.resize(ni,nj) ;
2297      if (latvalue_2d.numElements()>0) latvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2298    }
2299    else if (type=="boundslon")
2300    {
2301      CArray<double,1> value ;
2302      gathererConnector_->transfer(event, nvertex, value, 0.); 
2303      bounds_lon_2d.resize(nvertex,ni,nj) ;
2304      if (bounds_lon_2d.numElements()>0) bounds_lon_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2305    }
2306    else if (type=="boundslat")
2307    {
2308      CArray<double,1> value ;
2309      gathererConnector_->transfer(event, nvertex, value, 0.); 
2310      bounds_lat_2d.resize(nvertex,ni,nj) ;
2311      if (bounds_lat_2d.numElements()>0) bounds_lat_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2312    }
2313    else if (type=="area") 
2314    {
2315      CArray<double,1> value ;
2316      gathererConnector_->transfer(event, value, 0.); 
2317      area.resize(ni,nj) ;
2318      if (area.numElements()>0) area=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2319    }
2320  }
2321  CATCH
2322   
2323  bool CDomain::dispatchEvent(CEventServer& event)
2324  TRY
2325  {
2326    if (SuperClass::dispatchEvent(event)) return true;
2327    else
2328    {
2329      switch(event.type)
2330      {
2331        case EVENT_ID_DOMAIN_DISTRIBUTION:
2332          recvDomainDistribution(event);
2333          return true;
2334          break;
2335        case EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE:
2336          recvDistributedAttributes(event);
2337          return true;
2338          break; 
2339        default:
2340          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2341                << "Unknown Event");
2342          return false;
2343       }
2344    }
2345  }
2346  CATCH
2347
2348 
2349  /*!
2350    Compare two domain objects.
2351    They are equal if only if they have identical attributes as well as their values.
2352    Moreover, they must have the same transformations.
2353  \param [in] domain Compared domain
2354  \return result of the comparison
2355  */
2356  bool CDomain::isEqual(CDomain* obj)
2357  TRY
2358  {
2359    vector<StdString> excludedAttr;
2360    excludedAttr.push_back("domain_ref");
2361    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2362    if (!objEqual) return objEqual;
2363
2364    TransMapTypes thisTrans = this->getAllTransformations();
2365    TransMapTypes objTrans  = obj->getAllTransformations();
2366
2367    TransMapTypes::const_iterator it, itb, ite;
2368    std::vector<ETranformationType> thisTransType, objTransType;
2369    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2370      thisTransType.push_back(it->first);
2371    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2372      objTransType.push_back(it->first);
2373
2374    if (thisTransType.size() != objTransType.size()) return false;
2375    for (int idx = 0; idx < thisTransType.size(); ++idx)
2376      objEqual &= (thisTransType[idx] == objTransType[idx]);
2377
2378    return objEqual;
2379  }
2380  CATCH_DUMP_ATTR
2381
2382/////////////////////////////////////////////////////////////////////////
2383///////////////             TRANSFORMATIONS                    //////////
2384/////////////////////////////////////////////////////////////////////////
2385
2386  std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
2387  bool CDomain::dummyTransformationMapList_ = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
2388
2389  bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
2390  TRY
2391  {
2392    m["zoom_domain"] = TRANS_EXTRACT_DOMAIN;
2393    m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
2394    m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
2395    m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
2396    m["expand_domain"] = TRANS_EXPAND_DOMAIN;
2397    m["reorder_domain"] = TRANS_REORDER_DOMAIN;
2398    m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
2399    return true;
2400  }
2401  CATCH
2402
2403
2404  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2405  TRY
2406  {
2407    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2408    return transformationMap_.back().second;
2409  }
2410  CATCH_DUMP_ATTR
2411
2412  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, CTransformation<CDomain>* transformation)
2413  TRY
2414  {
2415    transformationMap_.push_back(std::make_pair(transType, transformation));
2416    return transformationMap_.back().second;
2417  }
2418  CATCH_DUMP_ATTR
2419  /*!
2420    Check whether a domain has transformation
2421    \return true if domain has transformation
2422  */
2423  bool CDomain::hasTransformation()
2424  TRY
2425  {
2426    return (!transformationMap_.empty());
2427  }
2428  CATCH_DUMP_ATTR
2429
2430  /*!
2431    Set transformation for current domain. It's the method to move transformation in hierarchy
2432    \param [in] domTrans transformation on domain
2433  */
2434  void CDomain::setTransformations(const TransMapTypes& domTrans)
2435  TRY
2436  {
2437    transformationMap_ = domTrans;
2438  }
2439  CATCH_DUMP_ATTR
2440
2441  /*!
2442    Get all transformation current domain has
2443    \return all transformation
2444  */
2445  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2446  TRY
2447  {
2448    return transformationMap_;
2449  }
2450  CATCH_DUMP_ATTR
2451
2452  void CDomain::duplicateTransformation(CDomain* src)
2453  TRY
2454  {
2455    if (src->hasTransformation())
2456    {
2457      this->setTransformations(src->getAllTransformations());
2458    }
2459  }
2460  CATCH_DUMP_ATTR
2461   
2462  /*!
2463   * Go through the hierarchy to find the domain from which the transformations must be inherited
2464   */
2465  void CDomain::solveInheritanceTransformation_old()
2466  TRY
2467  {
2468    if (hasTransformation() || !hasDirectDomainReference())
2469      return;
2470
2471    CDomain* domain = this;
2472    std::vector<CDomain*> refDomains;
2473    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2474    {
2475      refDomains.push_back(domain);
2476      domain = domain->getDirectDomainReference();
2477    }
2478
2479    if (domain->hasTransformation())
2480      for (size_t i = 0; i < refDomains.size(); ++i)
2481        refDomains[i]->setTransformations(domain->getAllTransformations());
2482  }
2483  CATCH_DUMP_ATTR
2484
2485
2486  void CDomain::solveInheritanceTransformation()
2487  TRY
2488  {
2489    if (solveInheritanceTransformation_done_) return;
2490    else solveInheritanceTransformation_done_=true ;
2491
2492    CDomain* domain = this;
2493    CDomain* Lastdomain ;
2494    std::list<CDomain*> refDomains;
2495    bool out=false ;
2496    vector<StdString> excludedAttr;
2497    excludedAttr.push_back("domain_ref");
2498   
2499    refDomains.push_front(domain) ;
2500    while (domain->hasDirectDomainReference() && !out)
2501    {
2502      CDomain* lastDomain=domain ;
2503      domain = domain->getDirectDomainReference();
2504      domain->solveRefInheritance() ;
2505      if (!domain->SuperClass::isEqual(lastDomain,excludedAttr)) out=true ;
2506      refDomains.push_front(domain) ;
2507    }
2508
2509    CTransformationPaths::TPath path ;
2510    auto& pathList = std::get<2>(path) ;
2511    std::get<0>(path) = EElement::DOMAIN ;
2512    std::get<1>(path) = refDomains.front()->getId() ;
2513    for (auto& domain : refDomains)
2514    {
2515      CDomain::TransMapTypes transformations = domain->getAllTransformations();
2516      for(auto& transformation : transformations) pathList.push_back({transformation.second->getTransformationType(), 
2517                                                                      transformation.second->getId()}) ;
2518    }
2519    transformationPaths_.addPath(path) ;
2520
2521  }
2522  CATCH_DUMP_ATTR
2523 
2524
2525  bool CDomain::activateFieldWorkflow(CGarbageCollector& gc)
2526  TRY
2527  {
2528    if (!domain_ref.isEmpty())
2529    {
2530      CField* field=getFieldFromId(domain_ref) ;
2531      if (field!=nullptr)
2532      {
2533        bool ret = field->buildWorkflowGraph(gc) ;
2534        if (!ret) return false ; // cannot build workflow graph at this state
2535      }
2536      else 
2537      {
2538        CDomain* domain = get(domain_ref) ;
2539        bool ret = domain->activateFieldWorkflow(gc) ;
2540        if (!ret) return false ; // cannot build workflow graph at this state
2541        domain_ref=domain->getId() ; // replace domain_ref by solved reference
2542      }
2543    }
2544    activateFieldWorkflow_done_=true ;
2545    return true ;
2546  }
2547  CATCH_DUMP_ATTR
2548/////////////////////////////////////////////////////////////////////////////////////////////
2549/////////////////////////////////////////////////////////////////////////////////////////////
2550
2551  void CDomain::setContextClient(CContextClient* contextClient)
2552  TRY
2553  {
2554    if (clientsSet.find(contextClient)==clientsSet.end())
2555    {
2556      clients.push_back(contextClient) ;
2557      clientsSet.insert(contextClient);
2558    }
2559  }
2560  CATCH_DUMP_ATTR
2561
2562  /*!
2563    Parse children nodes of a domain in xml file.
2564    Whenver there is a new transformation, its type and name should be added into this function
2565    \param node child node to process
2566  */
2567  void CDomain::parse(xml::CXMLNode & node)
2568  TRY
2569  {
2570    SuperClass::parse(node);
2571
2572    if (node.goToChildElement())
2573    {
2574      StdString nodeElementName;
2575      do
2576      {
2577        StdString nodeId("");
2578        if (node.getAttributes().end() != node.getAttributes().find("id"))
2579        { nodeId = node.getAttributes()["id"]; }
2580
2581        nodeElementName = node.getElementName();
2582        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2583        it = transformationMapList_.find(nodeElementName);
2584        if (ite != it)
2585        {
2586          if (it->first == "zoom_domain")
2587          {
2588            info(0) << "WARNING : " << it->first << " is deprecated, replaced by extract_domain." << endl;
2589          }
2590          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2591                                                                                                                nodeId,
2592                                                                                                                &node)));
2593        }
2594        else
2595        {
2596          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2597                << "The transformation " << nodeElementName << " has not been supported yet.");
2598        }
2599      } while (node.goToNextElement()) ;
2600      node.goToParentElement();
2601    }
2602  }
2603  CATCH_DUMP_ATTR
2604   //----------------------------------------------------------------
2605
2606   DEFINE_REF_FUNC(Domain,domain)
2607
2608   ///---------------------------------------------------------------
2609
2610} // namespace xios
Note: See TracBrowser for help on using the repository browser.