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

Last change on this file was 2632, checked in by ymipsl, 4 weeks ago
  • Recheck grid/domain/axis/scalar after reading grid from file, otherwise some attributes are not taking into account.

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