source: XIOS/dev/dev_ym/XIOS_COUPLING/src/node/domain.cpp @ 2304

Last change on this file since 2304 was 2304, checked in by ymipsl, 2 years ago

Fix problem in remote connector for read variable : supress redondance optimisation on remote connector in read case.
YM

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