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

Last change on this file since 2206 was 2206, checked in by ymipsl, 5 months ago

New feature : when can now use the syntax :
fieldId:domainId[n], in domain reference inside the workflow (XML). Same for axis and scalar.

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: 87.8 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 type)
1873  TRY
1874  {
1875    CContext* context = CContext::getCurrent();
1876    map<int, CArray<size_t,1>> globalIndex ;
1877
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    else if (type==EDistributionType::NONE) // domain is not distributed ie all servers get the same local domain
1912    {
1913      int nbServer = client->serverSize;
1914      int nglo=ni_glo*nj_glo ;
1915      CArray<size_t,1> indGlo ;
1916      for(size_t i=0;i<nglo;i++) indGlo(i) = i ;
1917      for (auto& rankServer : client->getRanksServerLeader()) globalIndex[rankServer] = indGlo ; 
1918    }
1919    remoteElement_[client] = new CDistributedElement(ni_glo*nj_glo, globalIndex) ;
1920    remoteElement_[client]->addFullView() ;
1921  }
1922  CATCH
1923
1924 
1925
1926  void CDomain::distributeToServer(CContextClient* client, map<int, CArray<size_t,1>>& globalIndex,
1927                                   CScattererConnector* &scattererConnector, const string& domainId)
1928  TRY
1929  {
1930    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
1931    CContext* context = CContext::getCurrent();
1932
1933    this->sendAllAttributesToServer(client, serverDomainId)  ;
1934
1935    CDistributedElement scatteredElement(ni_glo*nj_glo, globalIndex) ;
1936    scatteredElement.addFullView() ;
1937    scattererConnector = new CScattererConnector(localElement_->getView(CElementView::FULL), scatteredElement.getView(CElementView::FULL), 
1938                                           context->getIntraComm(), client->getRemoteSize()) ;
1939    scattererConnector->computeConnector() ;
1940
1941    // phase 0
1942    // send remote element to construct the full view on server, ie without hole
1943    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
1944    CMessage message0 ;
1945    message0<<serverDomainId<<0 ; 
1946    remoteElement_[client]->sendToServer(client,event0,message0) ; 
1947   
1948    // phase 1
1949    // send the full view of element to construct the connector which connect distributed data coming from client to the full local view
1950    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
1951    CMessage message1 ;
1952    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
1953    scattererConnector->transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
1954   
1955    sendDistributedAttributes(client, *scattererConnector, domainId) ;
1956
1957 
1958    // phase 2 send the mask : data index + mask2D
1959    CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
1960    CArray<bool,1> maskOut ;
1961    CLocalConnector workflowToFull(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
1962    workflowToFull.computeConnector() ;
1963    maskIn=true ;
1964    workflowToFull.transfer(maskIn,maskOut,false) ;
1965
1966
1967    // phase 3 : prepare grid scatterer connector to send data from client to server
1968    map<int,CArray<size_t,1>> workflowGlobalIndex ;
1969    map<int,CArray<bool,1>> maskOut2 ; 
1970    scattererConnector->transfer(maskOut, maskOut2, false) ;
1971    scatteredElement.addView(CElementView::WORKFLOW, maskOut2) ;
1972    scatteredElement.getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
1973    // create new workflow view for scattered element
1974    CDistributedElement clientToServerElement(scatteredElement.getGlobalSize(), workflowGlobalIndex) ;
1975    clientToServerElement.addFullView() ;
1976    CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
1977    CMessage message2 ;
1978    message2<<serverDomainId<<2 ; 
1979    clientToServerElement.sendToServer(client, event2, message2) ; 
1980    clientToServerConnector_[client] = new CScattererConnector(localElement_->getView(CElementView::WORKFLOW), clientToServerElement.getView(CElementView::FULL),
1981                                                               context->getIntraComm(), client->getRemoteSize()) ;
1982    clientToServerConnector_[client]->computeConnector() ;
1983
1984    clientFromServerConnector_[client] = new CGathererConnector(clientToServerElement.getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW));
1985    clientFromServerConnector_[client]->computeConnector() ;
1986
1987  }
1988  CATCH
1989 
1990  void CDomain::recvDomainDistribution(CEventServer& event)
1991  TRY
1992  {
1993    string domainId;
1994    int phasis ;
1995    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
1996    get(domainId)->receivedDomainDistribution(event, phasis);
1997  }
1998  CATCH
1999
2000  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2001  TRY
2002  {
2003    CContext* context = CContext::getCurrent();
2004    if (phasis==0) // receive the remote element to construct the full view
2005    {
2006      localElement_ = new  CLocalElement(context->getIntraCommRank(),event) ;
2007      localElement_->addFullView() ;
2008      // construct the local dimension and indexes
2009      auto& globalIndex=localElement_->getGlobalIndex() ;
2010      int nij=globalIndex.numElements() ;
2011      int minI=ni_glo,maxI=-1,minJ=nj_glo,maxJ=-1 ;
2012      int i,j ;
2013      int niGlo=ni_glo, njGlo=njGlo ;
2014      for(int ij=0;ij<nij;ij++)
2015      {
2016        j=globalIndex(ij)/niGlo ;
2017        i=globalIndex(ij)%niGlo ;
2018        if (i<minI) minI=i ;
2019        if (i>maxI) maxI=i ;
2020        if (j<minJ) minJ=j ;
2021        if (j>maxJ) maxJ=j ;
2022      } 
2023      if (maxI>=minI) { ibegin=minI ; ni=maxI-minI+1 ; }
2024      else {ibegin=0; ni=0 ;}
2025      if (maxJ>=minJ) { jbegin=minJ ; nj=maxJ-minJ+1 ; }
2026      else {jbegin=0; nj=0 ;}
2027
2028    }
2029    else if (phasis==1) // receive the sent view from client to construct the full distributed full view on server
2030    {
2031      CContext* context = CContext::getCurrent();
2032      CDistributedElement* elementFrom = new  CDistributedElement(event) ;
2033      elementFrom->addFullView() ;
2034      gathererConnector_ = new CGathererConnector(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2035      gathererConnector_->computeConnector() ; 
2036    }
2037    else if (phasis==2)
2038    {
2039//      delete gathererConnector_ ;
2040      elementFrom_ = new  CDistributedElement(event) ;
2041      elementFrom_->addFullView() ;
2042//      gathererConnector_ =  new CGathererConnector(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2043//      gathererConnector_ -> computeConnector() ;
2044    }
2045  }
2046  CATCH
2047
2048  void CDomain::setServerMask(CArray<bool,1>& serverMask, CContextClient* client)
2049  TRY
2050  {
2051    // nota : the client is needed to get the remote size for the scatterer connector. Maybe it is not the good place for this
2052    // Later, server to client connector can be computed on demand, with "client" as argument
2053    CContext* context = CContext::getCurrent();
2054    localElement_->addView(CElementView::WORKFLOW, serverMask) ;
2055    mask_1d.reference(serverMask.copy()) ;
2056 
2057    serverFromClientConnector_ = new CGathererConnector(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW)) ;
2058    serverFromClientConnector_->computeConnector() ;
2059     
2060    serverToClientConnector_ = new CScattererConnector(localElement_->getView(CElementView::WORKFLOW), elementFrom_->getView(CElementView::FULL),
2061                                                       context->getIntraComm(), client->getRemoteSize()) ;
2062    serverToClientConnector_->computeConnector() ;
2063  }
2064  CATCH_DUMP_ATTR
2065
2066
2067  void CDomain::sendDistributedAttributes(CContextClient* client, CScattererConnector& scattererConnector,  const string& domainId)
2068  {
2069    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2070    CContext* context = CContext::getCurrent();
2071
2072    if (hasLonLat)
2073    {
2074      { // send longitude
2075        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2076        CMessage message ;
2077        message<<serverDomainId<<string("lon") ; 
2078        scattererConnector.transfer(lonvalue, client, event,message) ;
2079      }
2080     
2081      { // send latitude
2082        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2083        CMessage message ;
2084        message<<serverDomainId<<string("lat") ; 
2085        scattererConnector.transfer(latvalue, client, event, message) ;
2086      }
2087    }
2088
2089    if (hasBounds)
2090    { 
2091      { // send longitude boudaries
2092        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2093        CMessage message ;
2094        message<<serverDomainId<<string("boundslon") ; 
2095        scattererConnector.transfer(nvertex, bounds_lonvalue, client, event, message ) ;
2096      }
2097
2098      { // send latitude boudaries
2099        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2100        CMessage message ;
2101        message<<serverDomainId<<string("boundslat") ; 
2102        scattererConnector.transfer(nvertex, bounds_latvalue, client, event, message ) ;
2103      }
2104    }
2105
2106    if (hasArea)
2107    {  // send area
2108      CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2109      CMessage message ;
2110      message<<serverDomainId<<string("area") ; 
2111      scattererConnector.transfer(areavalue, client, event,message) ;
2112    }
2113  }
2114
2115  void CDomain::recvDistributedAttributes(CEventServer& event)
2116  TRY
2117  {
2118    string domainId;
2119    string type ;
2120    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> type ;
2121    get(domainId)->recvDistributedAttributes(event, type);
2122  }
2123  CATCH
2124
2125  void CDomain::recvDistributedAttributes(CEventServer& event, const string& type)
2126  TRY
2127  {
2128    if (type=="lon") 
2129    {
2130      CArray<double,1> value ;
2131      gathererConnector_->transfer(event, value, 0.); 
2132      lonvalue_2d.resize(ni,nj) ;
2133      if (lonvalue_2d.numElements()>0) lonvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2134    }
2135    else if (type=="lat")
2136    {
2137      CArray<double,1> value ;
2138      gathererConnector_->transfer(event, value, 0.); 
2139      latvalue_2d.resize(ni,nj) ;
2140      if (latvalue_2d.numElements()>0) latvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2141    }
2142    else if (type=="boundslon")
2143    {
2144      CArray<double,1> value ;
2145      gathererConnector_->transfer(event, nvertex, value, 0.); 
2146      bounds_lon_2d.resize(nvertex,ni,nj) ;
2147      if (bounds_lon_2d.numElements()>0) bounds_lon_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2148    }
2149    else if (type=="boundslat")
2150    {
2151      CArray<double,1> value ;
2152      gathererConnector_->transfer(event, nvertex, value, 0.); 
2153      bounds_lat_2d.resize(nvertex,ni,nj) ;
2154      if (bounds_lat_2d.numElements()>0) bounds_lat_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2155    }
2156    else if (type=="area") 
2157    {
2158      CArray<double,1> value ;
2159      gathererConnector_->transfer(event, value, 0.); 
2160      area.resize(ni,nj) ;
2161      if (area.numElements()>0) area=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2162    }
2163  }
2164  CATCH
2165   
2166  bool CDomain::dispatchEvent(CEventServer& event)
2167  TRY
2168  {
2169    if (SuperClass::dispatchEvent(event)) return true;
2170    else
2171    {
2172      switch(event.type)
2173      {
2174        case EVENT_ID_DOMAIN_DISTRIBUTION:
2175          recvDomainDistribution(event);
2176          return true;
2177          break;
2178        case EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE:
2179          recvDistributedAttributes(event);
2180          return true;
2181          break; 
2182        default:
2183          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2184                << "Unknown Event");
2185          return false;
2186       }
2187    }
2188  }
2189  CATCH
2190
2191 
2192  /*!
2193    Compare two domain objects.
2194    They are equal if only if they have identical attributes as well as their values.
2195    Moreover, they must have the same transformations.
2196  \param [in] domain Compared domain
2197  \return result of the comparison
2198  */
2199  bool CDomain::isEqual(CDomain* obj)
2200  TRY
2201  {
2202    vector<StdString> excludedAttr;
2203    excludedAttr.push_back("domain_ref");
2204    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2205    if (!objEqual) return objEqual;
2206
2207    TransMapTypes thisTrans = this->getAllTransformations();
2208    TransMapTypes objTrans  = obj->getAllTransformations();
2209
2210    TransMapTypes::const_iterator it, itb, ite;
2211    std::vector<ETranformationType> thisTransType, objTransType;
2212    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2213      thisTransType.push_back(it->first);
2214    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2215      objTransType.push_back(it->first);
2216
2217    if (thisTransType.size() != objTransType.size()) return false;
2218    for (int idx = 0; idx < thisTransType.size(); ++idx)
2219      objEqual &= (thisTransType[idx] == objTransType[idx]);
2220
2221    return objEqual;
2222  }
2223  CATCH_DUMP_ATTR
2224
2225/////////////////////////////////////////////////////////////////////////
2226///////////////             TRANSFORMATIONS                    //////////
2227/////////////////////////////////////////////////////////////////////////
2228
2229  std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
2230  bool CDomain::dummyTransformationMapList_ = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
2231
2232  bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
2233  TRY
2234  {
2235    m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
2236    m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
2237    m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
2238    m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
2239    m["expand_domain"] = TRANS_EXPAND_DOMAIN;
2240    m["reorder_domain"] = TRANS_REORDER_DOMAIN;
2241    m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
2242    return true;
2243  }
2244  CATCH
2245
2246
2247  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2248  TRY
2249  {
2250    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2251    return transformationMap_.back().second;
2252  }
2253  CATCH_DUMP_ATTR
2254
2255  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, CTransformation<CDomain>* transformation)
2256  TRY
2257  {
2258    transformationMap_.push_back(std::make_pair(transType, transformation));
2259    return transformationMap_.back().second;
2260  }
2261  CATCH_DUMP_ATTR
2262  /*!
2263    Check whether a domain has transformation
2264    \return true if domain has transformation
2265  */
2266  bool CDomain::hasTransformation()
2267  TRY
2268  {
2269    return (!transformationMap_.empty());
2270  }
2271  CATCH_DUMP_ATTR
2272
2273  /*!
2274    Set transformation for current domain. It's the method to move transformation in hierarchy
2275    \param [in] domTrans transformation on domain
2276  */
2277  void CDomain::setTransformations(const TransMapTypes& domTrans)
2278  TRY
2279  {
2280    transformationMap_ = domTrans;
2281  }
2282  CATCH_DUMP_ATTR
2283
2284  /*!
2285    Get all transformation current domain has
2286    \return all transformation
2287  */
2288  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2289  TRY
2290  {
2291    return transformationMap_;
2292  }
2293  CATCH_DUMP_ATTR
2294
2295  void CDomain::duplicateTransformation(CDomain* src)
2296  TRY
2297  {
2298    if (src->hasTransformation())
2299    {
2300      this->setTransformations(src->getAllTransformations());
2301    }
2302  }
2303  CATCH_DUMP_ATTR
2304   
2305  /*!
2306   * Go through the hierarchy to find the domain from which the transformations must be inherited
2307   */
2308  void CDomain::solveInheritanceTransformation_old()
2309  TRY
2310  {
2311    if (hasTransformation() || !hasDirectDomainReference())
2312      return;
2313
2314    CDomain* domain = this;
2315    std::vector<CDomain*> refDomains;
2316    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2317    {
2318      refDomains.push_back(domain);
2319      domain = domain->getDirectDomainReference();
2320    }
2321
2322    if (domain->hasTransformation())
2323      for (size_t i = 0; i < refDomains.size(); ++i)
2324        refDomains[i]->setTransformations(domain->getAllTransformations());
2325  }
2326  CATCH_DUMP_ATTR
2327
2328
2329  void CDomain::solveInheritanceTransformation()
2330  TRY
2331  {
2332    if (solveInheritanceTransformation_done_) return;
2333    else solveInheritanceTransformation_done_=true ;
2334
2335    CDomain* domain = this;
2336    CDomain* Lastdomain ;
2337    std::list<CDomain*> refDomains;
2338    bool out=false ;
2339    vector<StdString> excludedAttr;
2340    excludedAttr.push_back("domain_ref");
2341   
2342    refDomains.push_front(domain) ;
2343    while (domain->hasDirectDomainReference() && !out)
2344    {
2345      CDomain* lastDomain=domain ;
2346      domain = domain->getDirectDomainReference();
2347      domain->solveRefInheritance() ;
2348      if (!domain->SuperClass::isEqual(lastDomain,excludedAttr)) out=true ;
2349      refDomains.push_front(domain) ;
2350    }
2351
2352    CTransformationPaths::TPath path ;
2353    auto& pathList = std::get<2>(path) ;
2354    std::get<0>(path) = EElement::DOMAIN ;
2355    std::get<1>(path) = refDomains.front()->getId() ;
2356    for (auto& domain : refDomains)
2357    {
2358      CDomain::TransMapTypes transformations = domain->getAllTransformations();
2359      for(auto& transformation : transformations) pathList.push_back({transformation.second->getTransformationType(), 
2360                                                                      transformation.second->getId()}) ;
2361    }
2362    transformationPaths_.addPath(path) ;
2363
2364  }
2365  CATCH_DUMP_ATTR
2366 
2367
2368  bool CDomain::activateFieldWorkflow(CGarbageCollector& gc)
2369  TRY
2370  {
2371    if (!domain_ref.isEmpty())
2372    {
2373      CField* field=getFieldFromId(domain_ref) ;
2374      if (field!=nullptr)
2375      {
2376        bool ret = field->buildWorkflowGraph(gc) ;
2377        if (!ret) return false ; // cannot build workflow graph at this state
2378      }
2379      else 
2380      {
2381        CDomain* domain = get(domain_ref) ;
2382        bool ret = domain->activateFieldWorkflow(gc) ;
2383        if (!ret) return false ; // cannot build workflow graph at this state
2384        domain_ref=domain->getId() ; // replace domain_ref by solved reference
2385      }
2386    }
2387    activateFieldWorkflow_done_=true ;
2388    return true ;
2389  }
2390  CATCH_DUMP_ATTR
2391/////////////////////////////////////////////////////////////////////////////////////////////
2392/////////////////////////////////////////////////////////////////////////////////////////////
2393
2394  void CDomain::setContextClient(CContextClient* contextClient)
2395  TRY
2396  {
2397    if (clientsSet.find(contextClient)==clientsSet.end())
2398    {
2399      clients.push_back(contextClient) ;
2400      clientsSet.insert(contextClient);
2401    }
2402  }
2403  CATCH_DUMP_ATTR
2404
2405  /*!
2406    Parse children nodes of a domain in xml file.
2407    Whenver there is a new transformation, its type and name should be added into this function
2408    \param node child node to process
2409  */
2410  void CDomain::parse(xml::CXMLNode & node)
2411  TRY
2412  {
2413    SuperClass::parse(node);
2414
2415    if (node.goToChildElement())
2416    {
2417      StdString nodeElementName;
2418      do
2419      {
2420        StdString nodeId("");
2421        if (node.getAttributes().end() != node.getAttributes().find("id"))
2422        { nodeId = node.getAttributes()["id"]; }
2423
2424        nodeElementName = node.getElementName();
2425        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2426        it = transformationMapList_.find(nodeElementName);
2427        if (ite != it)
2428        {
2429          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2430                                                                                                                nodeId,
2431                                                                                                                &node)));
2432        }
2433        else
2434        {
2435          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2436                << "The transformation " << nodeElementName << " has not been supported yet.");
2437        }
2438      } while (node.goToNextElement()) ;
2439      node.goToParentElement();
2440    }
2441  }
2442  CATCH_DUMP_ATTR
2443   //----------------------------------------------------------------
2444
2445   DEFINE_REF_FUNC(Domain,domain)
2446
2447   ///---------------------------------------------------------------
2448
2449} // namespace xios
Note: See TracBrowser for help on using the repository browser.