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

Last change on this file since 2247 was 2247, checked in by ymipsl, 9 months ago

Rewrite domain distribution for servers. Previously it was leading that each client was communicating with each servers to send them a small part of the domain.

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