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

Last change on this file since 2545 was 2545, checked in by ymipsl, 10 months ago

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