source: XIOS3/dev/XIOS_ATTACHED/src/node/domain.cpp

Last change on this file was 2482, checked in by ymipsl, 15 months ago

First guess in supression of attached mode replaced by online reader and write filters

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     checkBounds() ;
1356     checkArea() ;
1357
1358     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
1359     {
1360       lonvalue.resize(ni * nj);
1361       latvalue.resize(ni * nj);
1362       if (hasBounds)
1363       {
1364         bounds_lonvalue.resize(nvertex, ni * nj);
1365         bounds_latvalue.resize(nvertex, ni * nj);
1366       }
1367
1368       for (int j = 0; j < nj; ++j)
1369       {
1370         for (int i = 0; i < ni; ++i)
1371         {
1372           int k = j * ni + i;
1373
1374           lonvalue(k) = lonvalue_2d(i,j);
1375           latvalue(k) = latvalue_2d(i,j);
1376
1377           if (hasBounds)
1378           {
1379             for (int n = 0; n < nvertex; ++n)
1380             {
1381               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1382               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
1383             }
1384           }
1385         }
1386       }
1387     }
1388     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
1389     {
1390       if (type_attr::rectilinear == type)
1391       {
1392         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1393         {
1394           lonvalue.resize(ni * nj);
1395           latvalue.resize(ni * nj);
1396           if (hasBounds)
1397           {
1398             bounds_lonvalue.resize(nvertex, ni * nj);
1399             bounds_latvalue.resize(nvertex, ni * nj);
1400           }
1401
1402           for (int j = 0; j < nj; ++j)
1403           {
1404             for (int i = 0; i < ni; ++i)
1405             {
1406               int k = j * ni + i;
1407
1408               lonvalue(k) = lonvalue_1d(i);
1409               latvalue(k) = latvalue_1d(j);
1410
1411               if (hasBounds)
1412               {
1413                 for (int n = 0; n < nvertex; ++n)
1414                 {
1415                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1416                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1417                 }
1418               }
1419             }
1420           }
1421         }
1422         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1423         {
1424           lonvalue.reference(lonvalue_1d.copy());
1425           latvalue.reference(latvalue_1d.copy());
1426           if (hasBounds)
1427           {
1428             bounds_lonvalue.reference(bounds_lon_1d.copy());
1429             bounds_latvalue.reference(bounds_lat_1d.copy());
1430           }
1431         }
1432         else
1433           ERROR("CDomain::completeLonClient(void)",
1434                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1435                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1436                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1437                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1438                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1439                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1440       }
1441       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1442       {
1443         lonvalue.reference(lonvalue_1d.copy());
1444         latvalue.reference(latvalue_1d.copy());
1445         if (hasBounds)
1446         {
1447           bounds_lonvalue.reference(bounds_lon_1d.copy());
1448           bounds_latvalue.reference(bounds_lat_1d.copy());
1449         }
1450       }
1451     }
1452
1453     if (!area_2d.isEmpty() && areavalue.isEmpty())
1454     {
1455       areavalue.resize(ni*nj);
1456       for (int j = 0; j < nj; ++j)
1457       {
1458         for (int i = 0; i < ni; ++i)
1459         {
1460           int k = j * ni + i;
1461           areavalue(k) = area_2d(i,j);
1462         }
1463       }
1464     }
1465     else if (!area_1d.isEmpty() && areavalue.isEmpty()) areavalue.reference(area_1d.copy());
1466
1467   }
1468   CATCH_DUMP_ATTR
1469
1470   /*
1471     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1472   */
1473   void CDomain::convertLonLatValue(void)
1474   TRY
1475   {
1476     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1477     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1478     {
1479       lonvalue_2d.resize(ni,nj);
1480       latvalue_2d.resize(ni,nj);
1481       if (hasBounds)
1482       {
1483         bounds_lon_2d.resize(nvertex, ni, nj);
1484         bounds_lat_2d.resize(nvertex, ni, nj);
1485       }
1486
1487       for (int j = 0; j < nj; ++j)
1488       {
1489         for (int i = 0; i < ni; ++i)
1490         {
1491           int k = j * ni + i;
1492
1493           lonvalue_2d(i,j) = lonvalue(k);
1494           latvalue_2d(i,j) = latvalue(k);
1495
1496           if (hasBounds)
1497           {
1498             for (int n = 0; n < nvertex; ++n)
1499             {
1500               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1501               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1502             }
1503           }
1504         }
1505       }
1506     }
1507     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1508     {
1509       if (type_attr::rectilinear == type)
1510       {
1511         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1512         {
1513           lonvalue.resize(ni * nj);
1514           latvalue.resize(ni * nj);
1515           if (hasBounds)
1516           {
1517             bounds_lonvalue.resize(nvertex, ni * nj);
1518             bounds_latvalue.resize(nvertex, ni * nj);
1519           }
1520
1521           for (int j = 0; j < nj; ++j)
1522           {
1523             for (int i = 0; i < ni; ++i)
1524             {
1525               int k = j * ni + i;
1526
1527               lonvalue(k) = lonvalue_1d(i);
1528               latvalue(k) = latvalue_1d(j);
1529
1530               if (hasBounds)
1531               {
1532                 for (int n = 0; n < nvertex; ++n)
1533                 {
1534                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1535                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1536                 }
1537               }
1538             }
1539           }
1540         }
1541         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1542         {
1543           lonvalue.reference(lonvalue_1d);
1544           latvalue.reference(latvalue_1d);
1545            if (hasBounds)
1546           {
1547             bounds_lonvalue.reference(bounds_lon_1d);
1548             bounds_latvalue.reference(bounds_lat_1d);
1549           }
1550         }
1551         else
1552           ERROR("CDomain::completeLonClient(void)",
1553                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1554                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1555                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1556                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1557                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1558                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1559       }
1560       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1561       {
1562         lonvalue.reference(lonvalue_1d);
1563         latvalue.reference(latvalue_1d);
1564         if (hasBounds)
1565         {
1566           bounds_lonvalue.reference(bounds_lon_1d);
1567           bounds_latvalue.reference(bounds_lat_1d);
1568         }
1569       }
1570     }
1571   }
1572   CATCH_DUMP_ATTR
1573
1574   void CDomain::checkBounds(void)
1575   TRY
1576   {
1577     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1578     if (!nvertex.isEmpty() && nvertex > 0 && !hasBoundValues)
1579     {
1580       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1581         ERROR("CDomain::checkBounds(void)",
1582               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1583               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1584               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
1585
1586       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1587         ERROR("CDomain::checkBounds(void)",
1588               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1589               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1590               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
1591
1592       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1593       {
1594         ERROR("CDomain::checkBounds(void)",
1595               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1596               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1597               << "Please define either both attributes or none.");
1598       }
1599
1600       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1601       {
1602         ERROR("CDomain::checkBounds(void)",
1603               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1604               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1605               << "Please define either both attributes or none.");
1606       }
1607
1608       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
1609         ERROR("CDomain::checkBounds(void)",
1610               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1611               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1612               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(0)
1613               << " but nvertex is " << nvertex.getValue() << ".");
1614
1615       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
1616         ERROR("CDomain::checkBounds(void)",
1617               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1618               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1619               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(0)
1620               << " but nvertex is " << nvertex.getValue() << ".");
1621
1622       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
1623         ERROR("CDomain::checkBounds(void)",
1624               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1625               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
1626
1627       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
1628         ERROR("CDomain::checkBounds(void)",
1629               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1630               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
1631
1632       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
1633         ERROR("CDomain::checkBounds(void)",
1634               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1635               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1636               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(0)
1637               << " but nvertex is " << nvertex.getValue() << ".");
1638
1639       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
1640         ERROR("CDomain::checkBounds(void)",
1641               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1642               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1643               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(0)
1644               << " but nvertex is " << nvertex.getValue() << ".");
1645
1646       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
1647         ERROR("CDomain::checkBounds(void)",
1648               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1649               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
1650
1651       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
1652         ERROR("CDomain::checkBounds(void)",
1653               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1654
1655       // In case of reading UGRID bounds values are not required
1656       hasBounds = (!bounds_lat_1d.isEmpty() || !bounds_lat_2d.isEmpty() );
1657     }
1658     else if (hasBoundValues)
1659     {
1660       hasBounds = true;       
1661     }
1662     else
1663     {
1664       hasBounds = false;
1665     }
1666   }
1667   CATCH_DUMP_ATTR
1668
1669   void CDomain::checkArea(void)
1670   TRY
1671   {
1672     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
1673     hasArea = !area_1d.isEmpty() || !area_2d.isEmpty();
1674     if (hasArea && !hasAreaValue)
1675     {
1676       if (!area_2d.isEmpty() && (area_2d.extent(0) != ni || area_2d.extent(1) != nj))
1677       {
1678         ERROR("CDomain::checkArea(void)",
1679               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1680               << "The area does not have the same size as the local domain." << std::endl
1681               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1682               << "Area size is " << area_2d.extent(0) << " x " << area_2d.extent(1) << ".");
1683       }
1684       
1685       if (!area_1d.isEmpty() && area_1d.extent(0) != ni*nj)
1686       {
1687         ERROR("CDomain::checkArea(void)",
1688               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1689               << "The area does not have the same size as the local domain." << std::endl
1690               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1691               << "Area size is " << area_1d.extent(0) << " but must be ni*nj=" << ni*nj << " .");
1692       }
1693
1694     }
1695   }
1696   CATCH_DUMP_ATTR
1697
1698   void CDomain::checkLonLat()
1699   TRY
1700   {
1701     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1702                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1703     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1704     if (hasLonLat && !hasLonLatValue)
1705     {
1706       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1707         ERROR("CDomain::checkLonLat()",
1708               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1709               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1710               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
1711
1712       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1713       {
1714         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
1715           ERROR("CDomain::checkLonLat()",
1716                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1717                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1718                 << "Local size is " << i_index.numElements() << "." << std::endl
1719                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
1720       }
1721
1722       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1723       {
1724         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
1725           ERROR("CDomain::checkLonLat()",
1726                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1727                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1728                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1729                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
1730       }
1731
1732       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1733         ERROR("CDomain::checkLonLat()",
1734               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1735               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1736               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
1737
1738       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1739       {
1740         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
1741           ERROR("CDomain::checkLonLat()",
1742                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1743                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1744                 << "Local size is " << i_index.numElements() << "." << std::endl
1745                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
1746       }
1747
1748       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1749       {
1750         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
1751           ERROR("CDomain::checkLonLat()",
1752                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1753                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1754                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1755                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
1756       }
1757     }
1758   }
1759   CATCH_DUMP_ATTR
1760
1761   void CDomain::checkAttributes(void)
1762   TRY
1763   {
1764      if (this->checkAttributes_done_) return;
1765      this->checkDomain();
1766      this->compute2dBox() ;
1767      this->checkLonLat();
1768      this->checkBounds();
1769      this->checkArea();
1770      this->checkMask();
1771      this->checkDomainData();
1772      this->checkCompression();
1773      this->computeLocalMask() ;
1774      this->completeLonLatClient();
1775      this->initializeLocalElement() ;
1776      this->addFullView() ; // probably do not automatically add View, but only if requested
1777      this->addWorkflowView() ; // probably do not automatically add View, but only if requested
1778      this->addModelView() ; // probably do not automatically add View, but only if requested
1779      // testing ?
1780     /*
1781      shared_ptr<CLocalView> local = localElement_->getView(CElementView::WORKFLOW) ;
1782      shared_ptr<CLocalView> model = localElement_->getView(CElementView::MODEL) ;
1783
1784      CLocalConnector test1(model, local) ;
1785      test1.computeConnector() ;
1786      CLocalConnector test2(local, model) ;
1787      test2.computeConnector() ;
1788      CGridLocalConnector gridTest1(vector<CLocalConnector*>{&test1}) ;
1789      CGridLocalConnector gridTest2(vector<CLocalConnector*>{&test2}) ;
1790     
1791     
1792      CArray<int,1> out1 ;
1793      CArray<int,1> out2 ;
1794      test1.transfer(data_i_index,out1,-111) ;
1795      test2.transfer(out1,out2,-111) ;
1796     
1797      out1 = 0 ;
1798      out2 = 0 ;
1799      gridTest1.transfer(data_i_index,out1,-111) ;
1800      gridTest2.transfer(out1, out2,-111) ;
1801    */ 
1802      this->checkAttributes_done_ = true;
1803   }
1804   CATCH_DUMP_ATTR
1805
1806   size_t CDomain::computeAttributesHash( MPI_Comm comm )
1807   {
1808     // Compute the hash of distributed attributs (value ...)
1809     int globalSize = this->ni_glo.getValue()*this->nj_glo.getValue();
1810     CArray<size_t,1> globalIndex; // No redundancy globalIndex will be computed with the connector
1811     shared_ptr<CGridTransformConnector> gridTransformConnector;
1812     // Compute a without redundancy element FULL view to enable a consistent hash computation
1813     this->getLocalView(CElementView::FULL)->createWithoutRedundancyFullViewConnector( globalSize, comm, gridTransformConnector, globalIndex );
1814     int localSize = globalIndex.numElements();
1815           
1816     CArray<double,1> lon_distributedValue, lat_distributedValue ;
1817     gridTransformConnector->transfer(this->lonvalue, lon_distributedValue );
1818     gridTransformConnector->transfer(this->latvalue, lat_distributedValue );
1819
1820     // Compute the distributed hash (v0) of the element
1821     // it will be associated to the default element name (= map key), and to the name really written
1822     size_t localHash = 0;
1823     for (int iloc=0; iloc<localSize ; iloc++ ) localHash+=globalIndex(iloc)*lon_distributedValue(iloc)*lat_distributedValue(iloc);
1824     size_t distributedHash = 0;
1825     MPI_Allreduce( &localHash, &distributedHash, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm  );
1826     
1827     // Compute the hash of global attributs (unit, prec ...)
1828     vector<StdString> excludedAttr;
1829     //excludedAttr.push_back("name");
1830     // internal attributs
1831     excludedAttr.insert(excludedAttr.end(), { "ibegin", "jbegin", "ni", "nj", "i_index", "j_index" });
1832     excludedAttr.insert(excludedAttr.end(), { "data_ni", "data_nj", "data_ibegin", "data_jbegin" });
1833     excludedAttr.insert(excludedAttr.end(), { "data_i_index", "data_j_index", "domain_ref" });
1834     // in distributed through lonvalue and latvalue
1835     excludedAttr.insert(excludedAttr.end(), { "lonvalue_1d", "latvalue_1d", "lonvalue_2d", "latvalue_2d" });
1836     // should be considered in distributed
1837     excludedAttr.insert(excludedAttr.end(), { "mask_1d", "mask_2d" }); // ???
1838     excludedAttr.insert(excludedAttr.end(), { "bounds_lon_1d", "bounds_lat_1d", "bounds_lon_2d", "bounds_lat_2d" });
1839     excludedAttr.insert(excludedAttr.end(), { "area_1d", "area_2d" });
1840     // private
1841     excludedAttr.insert(excludedAttr.end(), { "lon_start", "lon_end", "lat_start", "lat_end" });
1842     excludedAttr.insert(excludedAttr.end(), { "bounds_lon_start", "bounds_lon_end", "bounds_lat_start", "bounds_lat_end" });
1843     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" });
1844     excludedAttr.insert(excludedAttr.end(), { "lonvalue_rectilinear_read_from_file", "latvalue_rectilinear_read_from_file", "lonvalue_curvilinear_read_from_file", "latvalue_curvilinear_read_from_file" });
1845     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" });
1846     excludedAttr.insert(excludedAttr.end(), { "bounds_lonvalue_unstructured_read_from_file", "bounds_latvalue_unstructured_read_from_file" });
1847     
1848     size_t globalHash = this->computeGlobalAttributesHash( excludedAttr );
1849
1850     return distributedHash + globalHash;
1851   }
1852 
1853   void CDomain::renameAttributesBeforeWriting(CDomain* writtenDomain)
1854   {
1855     if (writtenDomain!=NULL)
1856     {
1857       this->name = writtenDomain->getDomainOutputName();
1858       this->lon_name = writtenDomain->lon_name;
1859       this->lat_name = writtenDomain->lat_name;
1860       if (this->type == CDomain::type_attr::unstructured)
1861       {
1862         this->dim_i_name = writtenDomain->dim_i_name;
1863         this->nvertex_name = writtenDomain->nvertex_name;
1864       }
1865     }
1866     else
1867     {
1868       this->name =  this->getId();
1869       this->lon_name = "lon_"+this->getId();
1870       this->lat_name = "lat_"+this->getId();
1871       if (this->type == CDomain::type_attr::unstructured)
1872       {
1873         this->dim_i_name = "cell_"+this->getId();
1874         this->nvertex_name = "nvertex_"+this->getId();
1875       }
1876     }
1877   }
1878
1879   void CDomain::initializeLocalElement(void)
1880   {
1881      // after checkDomain i_index and j_index of size (ni*nj)
1882      int nij = ni*nj ;
1883      CArray<size_t, 1> ij_index(ni*nj) ;
1884      for(int ij=0; ij<nij ; ij++) ij_index(ij) = i_index(ij)+j_index(ij)*ni_glo ;
1885      int rank = CContext::getCurrent()->getIntraCommRank() ;
1886      localElement_ = make_shared<CLocalElement>(rank, ni_glo*nj_glo, ij_index) ;
1887   }
1888
1889   void CDomain::addFullView(void)
1890   {
1891      CArray<int,1> index(ni*nj) ;
1892      int nij=ni*nj ;
1893      for(int ij=0; ij<nij ; ij++) index(ij)=ij ;
1894      localElement_ -> addView(CElementView::FULL, index) ;
1895   }
1896
1897   void CDomain::addWorkflowView(void)
1898   {
1899     // information for workflow view is stored in localMask
1900     int nij=ni*nj ;
1901     int nMask=0 ;
1902     for(int ij=0; ij<nij ; ij++) if (localMask(ij)) nMask++ ;
1903     CArray<int,1> index(nMask) ;
1904
1905     nMask=0 ;
1906     for(int ij=0; ij<nij ; ij++) 
1907      if (localMask(ij))
1908      {
1909        index(nMask)=ij ;
1910        nMask++ ;
1911      }
1912      localElement_ -> addView(CElementView::WORKFLOW, index) ;
1913   }
1914
1915   void CDomain::addModelView(void)
1916   {
1917     // information for model view is stored in data_i_index/data_j_index
1918     // very weird, do not mix data_i_index and data_i_begin => in future only keep data_i_index
1919     int dataSize = data_i_index.numElements() ;
1920     CArray<int,1> index(dataSize) ;
1921     int i,j ;
1922     for(int k=0;k<dataSize;k++)
1923     {
1924        if (data_dim==2)
1925        {
1926          i=data_i_index(k)+data_ibegin ; // bad
1927          j=data_j_index(k)+data_jbegin ; // bad
1928          if (i>=0 && i<ni && j>=0 && j<nj) index(k)=i+j*ni ;
1929          else index(k)=-1 ;
1930        }
1931        else if (data_dim==1)
1932        {
1933          i=data_i_index(k)+data_ibegin ; // bad
1934          if (i>=0 && i<ni*nj) index(k)=i ;
1935          else index(k)=-1 ;
1936        }
1937     }
1938     localElement_->addView(CElementView::MODEL, index) ;
1939   }
1940       
1941   void CDomain::computeModelToWorkflowConnector(void)
1942   { 
1943     shared_ptr<CLocalView> srcView=getLocalView(CElementView::MODEL) ;
1944     shared_ptr<CLocalView> dstView=getLocalView(CElementView::WORKFLOW) ;
1945     modelToWorkflowConnector_ = make_shared<CLocalConnector>(srcView, dstView); 
1946     modelToWorkflowConnector_->computeConnector() ;
1947   }
1948
1949
1950  string CDomain::getCouplingAlias(const string& fieldId, int posInGrid)
1951  {
1952    return "_domain["+std::to_string(posInGrid)+"]_of_"+fieldId ;
1953  }
1954   
1955  /* to be removed later when coupling will be reimplemented, just to  not forget */
1956  void CDomain::sendDomainToCouplerOut(CContextClient* client, const string& fieldId, int posInGrid)
1957  {
1958    if (sendDomainToFileServer_done_.count(client)!=0) return ;
1959    else sendDomainToFileServer_done_.insert(client) ;
1960   
1961    const string domainId = getCouplingAlias(fieldId, posInGrid) ;
1962   
1963    if (!domain_ref.isEmpty())
1964    {
1965      auto domain_ref_tmp=domain_ref.getValue() ;
1966      domain_ref.reset() ; // remove the reference, find an other way to do that more cleanly
1967      this->sendAllAttributesToServer(client, domainId)  ;
1968      domain_ref = domain_ref_tmp ;
1969    }
1970    else this->sendAllAttributesToServer(client, domainId)  ;
1971  }
1972
1973
1974
1975
1976  void CDomain::makeAliasForCoupling(const string& fieldId, int posInGrid)
1977  {
1978    const string domainId = getCouplingAlias(fieldId, posInGrid);
1979    this->createAlias(domainId) ;
1980  }
1981
1982
1983  void CDomain::computeRemoteElement(CContextClient* client, EDistributionType distType)
1984  TRY
1985  {
1986    CContext* context = CContext::getCurrent();
1987    map<int, CArray<size_t,1>> globalIndex ;
1988
1989    if (distType==EDistributionType::BANDS && isUnstructed_) distType=EDistributionType::COLUMNS ;
1990
1991    if (distType==EDistributionType::BANDS) // Bands distribution to send to file server
1992    {
1993      int nbServer = client->getRemoteSize();
1994      int nbClient = client->getIntraCommSize() ;
1995      int rankClient = client->getIntraCommRank() ;
1996
1997      int nbChunk = max( nbServer, nbClient);
1998      int size ;//= nbChunk / nbClient ;
1999      int start ;
2000     
2001      if (nbServer<=nbClient) // nbChunk = nbClient
2002      {
2003        // distribution regarding client only : 1 chunk per client, one server can recv from many clients (see serverRank below)
2004        size = 1;
2005        start = rankClient;
2006      }
2007      else // nbChunk = nbServer
2008      {
2009        // distribution regarding servers : 1 client with size(+1) chunk will send to size(+1)  servers
2010        size = nbChunk/nbClient;
2011        int nClientWithAdditionalChunk = nbChunk - size*nbClient;
2012        if (rankClient<nClientWithAdditionalChunk) // distribute size+1 chunks per client
2013        {
2014          size++;
2015          start =  size*rankClient;
2016        }
2017        else // distribute the rest with size chunks
2018        {
2019          start =  (size+1)*nClientWithAdditionalChunk+size*(rankClient-nClientWithAdditionalChunk);
2020        }
2021      }
2022     
2023      for(int i=0; i<size; i++)
2024      { 
2025        int rank=start+i ; 
2026        size_t indSize = nj_glo/nbChunk ;
2027        size_t indStart ;
2028        if (nj_glo % nbChunk > rank)
2029        {
2030          indStart = (indSize+1) * rank ;
2031          indSize++ ;
2032        }
2033        else indStart = indSize*rank + nj_glo%nbChunk ;
2034       
2035        indStart=indStart*ni_glo ;
2036        indSize=indSize*ni_glo ;
2037        // compute the server rank concerns by the chunk
2038        int serverRank;
2039        if (nbServer<=nbClient)
2040        {
2041          // nChunkPerServer(+1) client will send to 1 server
2042          if (nj_glo<nbChunk)
2043              nbChunk=nj_glo;
2044          int nChunkPerServer = nbChunk/nbServer;
2045          int nServerWithAdditionalChunk = nbChunk - nChunkPerServer*nbServer;
2046          if (rankClient<nServerWithAdditionalChunk*(nChunkPerServer+1))
2047          {
2048            serverRank = rank/(nChunkPerServer+1);
2049          }
2050          else
2051          {
2052            // nServerWithAdditionalChunk servers have been affected above with (nChunkPerServer+1) chunks
2053            // the rest will recv nChunkPerServer
2054            if (nChunkPerServer>0)
2055            {
2056              serverRank = nServerWithAdditionalChunk+(rank-nServerWithAdditionalChunk*(nChunkPerServer+1))/nChunkPerServer;
2057            }
2058            else
2059            {
2060              // no chunk for all servers, current rank will not manage informations for this domain
2061              serverRank = client->getRemoteSize();
2062            }
2063          }
2064        }
2065        else // (nbServer > nbClient)
2066        {
2067          serverRank = rank;
2068        }
2069
2070        if (nbChunk<nbServer)
2071        {
2072            if ( (serverRank==client->getRemoteSize()) && (rankClient<nbServer) )
2073            {
2074                indSize = 0;
2075                serverRank = rank;
2076            }
2077        }
2078
2079        if (serverRank<client->getRemoteSize())
2080        {
2081          auto& globalInd =  globalIndex[serverRank] ;
2082          globalInd.resize(indSize) ;
2083          for(size_t n = 0 ; n<indSize; n++) globalInd(n)=indStart+n ;
2084        }
2085      }
2086    }
2087    else if (distType==EDistributionType::COLUMNS) // Bands distribution to send to file server
2088    {
2089      int nbServer = client->getRemoteSize();
2090      int nbClient = client->getIntraCommSize() ;
2091      int rankClient = client->getIntraCommRank() ;
2092      int size = nbServer / nbClient ;
2093      int start ;
2094      if (nbServer%nbClient > rankClient)
2095      {
2096       start = (size+1) * rankClient ;
2097       size++ ;
2098      }
2099      else start = size*rankClient + nbServer%nbClient ;
2100     
2101      for(int i=0; i<size; i++)
2102      { 
2103        int rank=start+i ; 
2104        size_t indSize = ni_glo/nbServer ;
2105        size_t indStart ;
2106        if (ni_glo % nbServer > rank)
2107        {
2108          indStart = (indSize+1) * rank ;
2109          indSize++ ;
2110        }
2111        else indStart = indSize*rank + ni_glo%nbServer ;
2112       
2113        auto& globalInd =  globalIndex[rank] ;
2114        globalInd.resize(indSize*nj_glo) ;
2115        size_t n=0 ;
2116        for(int j=0; j<nj_glo;j++)
2117        {
2118          for(int i=0; i<indSize; i++, n++)  globalInd(n)=indStart+i ;
2119          indStart=indStart+ni_glo ; 
2120        }
2121      }
2122    }
2123    else if (distType==EDistributionType::NONE) // domain is not distributed ie all servers get the same local domain
2124    {
2125      int nbServer = client->getRemoteSize();
2126      int nglo=ni_glo*nj_glo ;
2127      CArray<size_t,1> indGlo ;
2128      for(size_t i=0;i<nglo;i++) indGlo(i) = i ;
2129      for (auto& rankServer : client->getRanksServerLeader()) globalIndex[rankServer] = indGlo ; 
2130    }
2131    remoteElement_[client] = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndex) ;
2132    remoteElement_[client]->addFullView() ;
2133  }
2134  CATCH
2135
2136  void CDomain::distributeToServer(CContextClient* client, bool inOut, map<int, CArray<size_t,1>>& globalIndexOut, std::map<int, CArray<size_t,1>>& globalIndexIn,
2137                                   shared_ptr<CScattererConnector> &scattererConnector, const string& domainId)
2138  TRY
2139  {
2140    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2141    CContext* context = CContext::getCurrent();
2142
2143    this->sendAllAttributesToServer(client, serverDomainId)  ;
2144
2145    auto scatteredElement = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndexOut) ;
2146    scatteredElement->addFullView() ;
2147    scattererConnector = make_shared<CScattererConnector>(localElement_->getView(CElementView::FULL), scatteredElement->getView(CElementView::FULL), 
2148                                                          context->getIntraComm(), client->getRemoteSize()) ;
2149    scattererConnector->computeConnector() ;
2150
2151    // phase 0
2152    // send remote element to construct the full view on server, ie without hole
2153    CEventClient event0(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2154    CMessage message0 ;
2155    message0<<serverDomainId<<0 ; 
2156    remoteElement_[client]->sendToServer(client,event0,message0) ; 
2157   
2158    // phase 1
2159    // send the full view of element to construct the connector which connect distributed data coming from client to the full local view
2160    CEventClient event1(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2161    CMessage message1 ;
2162    message1<<serverDomainId<<1<<localElement_->getView(CElementView::FULL)->getGlobalSize() ; 
2163    scattererConnector->transfer(localElement_->getView(CElementView::FULL)->getGlobalIndex(),client,event1,message1) ;
2164   
2165    sendDistributedAttributes(client, scattererConnector, domainId) ;
2166
2167 
2168    // phase 2 send the mask : data index + mask2D
2169    {
2170      CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2171      CArray<bool,1> maskOut ;
2172      auto workflowToFull = make_shared<CLocalConnector>(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2173      workflowToFull->computeConnector() ;
2174      maskIn=true ;
2175      workflowToFull->transfer(maskIn,maskOut,false) ;
2176
2177
2178      // prepare grid scatterer connector to send data from client to server
2179      map<int,CArray<size_t,1>> workflowGlobalIndex ;
2180      map<int,CArray<bool,1>> maskOut2 ; 
2181      scattererConnector->transfer(maskOut, maskOut2, false) ;
2182      scatteredElement->addView(CElementView::WORKFLOW, maskOut2) ;
2183      scatteredElement->getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2184      // create new workflow view for scattered element
2185      auto clientToServerElement = make_shared<CDistributedElement>(scatteredElement->getGlobalSize(), workflowGlobalIndex) ;
2186      clientToServerElement->addFullView() ;
2187      CEventClient event2(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2188      CMessage message2 ;
2189      message2<<serverDomainId<<2 ; 
2190      clientToServerElement->sendToServer(client, event2, message2) ; 
2191      clientToServerConnector_[client] = make_shared<CScattererConnector>(localElement_->getView(CElementView::WORKFLOW), clientToServerElement->getView(CElementView::FULL),
2192                                                                        context->getIntraComm(), client->getRemoteSize()) ;
2193    clientToServerConnector_[client]->computeConnector() ;
2194    }
2195    ////////////
2196    // phase 3 : compute connector to receive from server
2197    ////////////
2198    if (inOut)
2199    {
2200      auto scatteredElement = make_shared<CDistributedElement>(ni_glo*nj_glo, globalIndexIn) ;
2201      scatteredElement->addFullView() ;
2202      auto scattererConnector = make_shared<CScattererConnector>(localElement_->getView(CElementView::FULL), scatteredElement->getView(CElementView::FULL), 
2203                                                                 context->getIntraComm(), client->getRemoteSize()) ;
2204      scattererConnector->computeConnector() ;
2205
2206      CArray<bool,1> maskIn(localElement_->getView(CElementView::WORKFLOW)->getSize());
2207      CArray<bool,1> maskOut ;
2208      auto workflowToFull = make_shared<CLocalConnector>(localElement_->getView(CElementView::WORKFLOW), localElement_->getView(CElementView::FULL)) ;
2209      workflowToFull->computeConnector() ;
2210      maskIn=true ;
2211      workflowToFull->transfer(maskIn,maskOut,false) ;
2212
2213      map<int,CArray<size_t,1>> workflowGlobalIndex ;
2214      map<int,CArray<bool,1>> maskOut2 ; 
2215      scattererConnector->transfer(maskOut, maskOut2, false) ;
2216      scatteredElement->addView(CElementView::WORKFLOW, maskOut2) ;
2217      scatteredElement->getView(CElementView::WORKFLOW)->getGlobalIndexView(workflowGlobalIndex) ;
2218      auto clientToServerElement = make_shared<CDistributedElement>(scatteredElement->getGlobalSize(), workflowGlobalIndex) ;
2219      clientToServerElement->addFullView() ;
2220      CEventClient event3(getType(), EVENT_ID_DOMAIN_DISTRIBUTION);
2221      CMessage message3 ;
2222      message3<<serverDomainId<<3 ; 
2223      clientToServerElement->sendToServer(client, event3, message3) ; 
2224
2225      clientFromServerConnector_[client] = make_shared<CGathererConnector>(clientToServerElement->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW));
2226      clientFromServerConnector_[client]->computeConnector() ;     
2227    }
2228  }
2229  CATCH
2230 
2231  void CDomain::recvDomainDistribution(CEventServer& event)
2232  TRY
2233  {
2234    string domainId;
2235    int phasis ;
2236    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> phasis ;
2237    get(domainId)->receivedDomainDistribution(event, phasis);
2238  }
2239  CATCH
2240
2241  void CDomain::receivedDomainDistribution(CEventServer& event, int phasis)
2242  TRY
2243  {
2244    CContext* context = CContext::getCurrent();
2245    if (phasis==0) // receive the remote element to construct the full view
2246    {
2247      localElement_ = make_shared<CLocalElement>(context->getIntraCommRank(),event) ;
2248      localElement_->addFullView() ;
2249      // construct the local dimension and indexes
2250      auto& globalIndex=localElement_->getGlobalIndex() ;
2251      int nij=globalIndex.numElements() ;
2252      int minI=ni_glo,maxI=-1,minJ=nj_glo,maxJ=-1 ;
2253      int i,j ;
2254      int niGlo=ni_glo, njGlo=njGlo ;
2255      for(int ij=0;ij<nij;ij++)
2256      {
2257        j=globalIndex(ij)/niGlo ;
2258        i=globalIndex(ij)%niGlo ;
2259        if (i<minI) minI=i ;
2260        if (i>maxI) maxI=i ;
2261        if (j<minJ) minJ=j ;
2262        if (j>maxJ) maxJ=j ;
2263      } 
2264      if (maxI>=minI) { ibegin=minI ; ni=maxI-minI+1 ; }
2265      else {ibegin=0; ni=0 ;}
2266      if (maxJ>=minJ) { jbegin=minJ ; nj=maxJ-minJ+1 ; }
2267      else {jbegin=0; nj=0 ;}
2268
2269    }
2270    else if (phasis==1) // receive the sent view from client to construct the full distributed full view on server
2271    {
2272      CContext* context = CContext::getCurrent();
2273      shared_ptr<CDistributedElement> elementFrom = make_shared<CDistributedElement>(event) ;
2274      elementFrom->addFullView() ;
2275      gathererConnector_ = make_shared<CGathererConnector>(elementFrom->getView(CElementView::FULL), localElement_->getView(CElementView::FULL)) ;
2276      gathererConnector_->computeConnector() ; 
2277    }
2278    else if (phasis==2)
2279    {
2280      elementFrom_ = make_shared<CDistributedElement>(event) ;
2281      elementFrom_->addFullView() ;
2282    }
2283    else if (phasis==3)
2284    {
2285      elementTo_ = make_shared<CDistributedElement>(event) ;
2286      elementTo_->addFullView() ;
2287    }
2288  }
2289  CATCH
2290
2291  void CDomain::setServerMask(CArray<bool,1>& serverMask, CContextClient* client)
2292  TRY
2293  {
2294    // nota : the client is needed to get the remote size for the scatterer connector. Maybe it is not the good place for this
2295    // Later, server to client connector can be computed on demand, with "client" as argument
2296    CContext* context = CContext::getCurrent();
2297    localElement_->addView(CElementView::WORKFLOW, serverMask) ;
2298    mask_1d.reference(serverMask.copy()) ;
2299 
2300    serverFromClientConnector_ = make_shared<CGathererConnector>(elementFrom_->getView(CElementView::FULL), localElement_->getView(CElementView::WORKFLOW)) ;
2301    serverFromClientConnector_->computeConnector() ;
2302    elementFrom_.reset() ;
2303     
2304    if (elementTo_)
2305    {
2306      serverToClientConnector_ = make_shared<CScattererConnector>(localElement_->getView(CElementView::WORKFLOW), elementTo_->getView(CElementView::FULL),
2307                                                                context->getIntraComm(), client->getRemoteSize()) ;
2308      serverToClientConnector_->computeConnector() ;
2309      elementTo_.reset() ;
2310    }
2311
2312  }
2313  CATCH_DUMP_ATTR
2314
2315
2316  void CDomain::sendDistributedAttributes(CContextClient* client, shared_ptr<CScattererConnector> scattererConnector,  const string& domainId)
2317  {
2318    string serverDomainId = domainId.empty() ? this->getId() : domainId ;
2319    CContext* context = CContext::getCurrent();
2320
2321    if (hasLonLat)
2322    {
2323      { // send longitude
2324        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2325        CMessage message ;
2326        message<<serverDomainId<<string("lon") ; 
2327        scattererConnector->transfer(lonvalue, client, event,message) ;
2328      }
2329     
2330      { // send latitude
2331        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2332        CMessage message ;
2333        message<<serverDomainId<<string("lat") ; 
2334        scattererConnector->transfer(latvalue, client, event, message) ;
2335      }
2336    }
2337
2338    if (hasBounds)
2339    { 
2340      { // send longitude boudaries
2341        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2342        CMessage message ;
2343        message<<serverDomainId<<string("boundslon") ; 
2344        scattererConnector->transfer(nvertex, bounds_lonvalue, client, event, message ) ;
2345      }
2346
2347      { // send latitude boudaries
2348        CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2349        CMessage message ;
2350        message<<serverDomainId<<string("boundslat") ; 
2351        scattererConnector->transfer(nvertex, bounds_latvalue, client, event, message ) ;
2352      }
2353    }
2354
2355    if (hasArea)
2356    {  // send area
2357      CEventClient event(getType(), EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE);
2358      CMessage message ;
2359      message<<serverDomainId<<string("area") ; 
2360      scattererConnector->transfer(areavalue, client, event,message) ;
2361    }
2362  }
2363
2364  void CDomain::recvDistributedAttributes(CEventServer& event)
2365  TRY
2366  {
2367    string domainId;
2368    string type ;
2369    for (auto& subEvent : event.subEvents) (*subEvent.buffer) >> domainId >> type ;
2370    get(domainId)->recvDistributedAttributes(event, type);
2371  }
2372  CATCH
2373
2374  void CDomain::recvDistributedAttributes(CEventServer& event, const string& type)
2375  TRY
2376  {
2377    if (type=="lon") 
2378    {
2379      CArray<double,1> value ;
2380      gathererConnector_->transfer(event, value, 0.); 
2381      lonvalue_2d.resize(ni,nj) ;
2382      if (lonvalue_2d.numElements()>0) lonvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2383    }
2384    else if (type=="lat")
2385    {
2386      CArray<double,1> value ;
2387      gathererConnector_->transfer(event, value, 0.); 
2388      latvalue_2d.resize(ni,nj) ;
2389      if (latvalue_2d.numElements()>0) latvalue_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2390    }
2391    else if (type=="boundslon")
2392    {
2393      CArray<double,1> value ;
2394      gathererConnector_->transfer(event, nvertex, value, 0.); 
2395      bounds_lon_2d.resize(nvertex,ni,nj) ;
2396      if (bounds_lon_2d.numElements()>0) bounds_lon_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2397    }
2398    else if (type=="boundslat")
2399    {
2400      CArray<double,1> value ;
2401      gathererConnector_->transfer(event, nvertex, value, 0.); 
2402      bounds_lat_2d.resize(nvertex,ni,nj) ;
2403      if (bounds_lat_2d.numElements()>0) bounds_lat_2d=CArray<double,3>(value.dataFirst(),shape(nvertex,ni,nj),neverDeleteData) ; 
2404    }
2405    else if (type=="area") 
2406    {
2407      CArray<double,1> value ;
2408      gathererConnector_->transfer(event, value, 0.); 
2409      area_2d.resize(ni,nj) ;
2410      if (area_2d.numElements()>0) area_2d=CArray<double,2>(value.dataFirst(),shape(ni,nj),neverDeleteData) ; 
2411    }
2412  }
2413  CATCH
2414   
2415  bool CDomain::dispatchEvent(CEventServer& event)
2416  TRY
2417  {
2418    if (SuperClass::dispatchEvent(event)) return true;
2419    else
2420    {
2421      switch(event.type)
2422      {
2423        case EVENT_ID_DOMAIN_DISTRIBUTION:
2424          recvDomainDistribution(event);
2425          return true;
2426          break;
2427        case EVENT_ID_SEND_DISTRIBUTED_ATTRIBUTE:
2428          recvDistributedAttributes(event);
2429          return true;
2430          break; 
2431        default:
2432          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
2433                << "Unknown Event");
2434          return false;
2435       }
2436    }
2437  }
2438  CATCH
2439
2440 
2441  /*!
2442    Compare two domain objects.
2443    They are equal if only if they have identical attributes as well as their values.
2444    Moreover, they must have the same transformations.
2445  \param [in] domain Compared domain
2446  \return result of the comparison
2447  */
2448  bool CDomain::isEqual(CDomain* obj)
2449  TRY
2450  {
2451    vector<StdString> excludedAttr;
2452    excludedAttr.push_back("domain_ref");
2453    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
2454    if (!objEqual) return objEqual;
2455
2456    TransMapTypes thisTrans = this->getAllTransformations();
2457    TransMapTypes objTrans  = obj->getAllTransformations();
2458
2459    TransMapTypes::const_iterator it, itb, ite;
2460    std::vector<ETranformationType> thisTransType, objTransType;
2461    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2462      thisTransType.push_back(it->first);
2463    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2464      objTransType.push_back(it->first);
2465
2466    if (thisTransType.size() != objTransType.size()) return false;
2467    for (int idx = 0; idx < thisTransType.size(); ++idx)
2468      objEqual &= (thisTransType[idx] == objTransType[idx]);
2469
2470    return objEqual;
2471  }
2472  CATCH_DUMP_ATTR
2473
2474/////////////////////////////////////////////////////////////////////////
2475///////////////             TRANSFORMATIONS                    //////////
2476/////////////////////////////////////////////////////////////////////////
2477
2478  std::map<StdString, ETranformationType> CDomain::transformationMapList_ = std::map<StdString, ETranformationType>();
2479  bool CDomain::dummyTransformationMapList_ = CDomain::initializeTransformationMap(CDomain::transformationMapList_);
2480
2481  bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
2482  TRY
2483  {
2484    m["zoom_domain"] = TRANS_EXTRACT_DOMAIN;
2485    m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
2486    m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
2487    m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
2488    m["expand_domain"] = TRANS_EXPAND_DOMAIN;
2489    m["reorder_domain"] = TRANS_REORDER_DOMAIN;
2490    m["extract_domain"] = TRANS_EXTRACT_DOMAIN;
2491    m["redistribute_domain"] = TRANS_REDISTRIBUTE_DOMAIN;
2492    return true;
2493  }
2494  CATCH
2495
2496
2497  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
2498  TRY
2499  {
2500    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
2501    return transformationMap_.back().second;
2502  }
2503  CATCH_DUMP_ATTR
2504
2505  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, CTransformation<CDomain>* transformation)
2506  TRY
2507  {
2508    transformationMap_.push_back(std::make_pair(transType, transformation));
2509    return transformationMap_.back().second;
2510  }
2511  CATCH_DUMP_ATTR
2512  /*!
2513    Check whether a domain has transformation
2514    \return true if domain has transformation
2515  */
2516  bool CDomain::hasTransformation()
2517  TRY
2518  {
2519    return (!transformationMap_.empty());
2520  }
2521  CATCH_DUMP_ATTR
2522
2523  /*!
2524    Set transformation for current domain. It's the method to move transformation in hierarchy
2525    \param [in] domTrans transformation on domain
2526  */
2527  void CDomain::setTransformations(const TransMapTypes& domTrans)
2528  TRY
2529  {
2530    transformationMap_ = domTrans;
2531  }
2532  CATCH_DUMP_ATTR
2533
2534  /*!
2535    Get all transformation current domain has
2536    \return all transformation
2537  */
2538  CDomain::TransMapTypes CDomain::getAllTransformations(void)
2539  TRY
2540  {
2541    return transformationMap_;
2542  }
2543  CATCH_DUMP_ATTR
2544
2545  void CDomain::duplicateTransformation(CDomain* src)
2546  TRY
2547  {
2548    if (src->hasTransformation())
2549    {
2550      this->setTransformations(src->getAllTransformations());
2551    }
2552  }
2553  CATCH_DUMP_ATTR
2554   
2555  /*!
2556   * Go through the hierarchy to find the domain from which the transformations must be inherited
2557   */
2558  void CDomain::solveInheritanceTransformation_old()
2559  TRY
2560  {
2561    if (hasTransformation() || !hasDirectDomainReference())
2562      return;
2563
2564    CDomain* domain = this;
2565    std::vector<CDomain*> refDomains;
2566    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
2567    {
2568      refDomains.push_back(domain);
2569      domain = domain->getDirectDomainReference();
2570    }
2571
2572    if (domain->hasTransformation())
2573      for (size_t i = 0; i < refDomains.size(); ++i)
2574        refDomains[i]->setTransformations(domain->getAllTransformations());
2575  }
2576  CATCH_DUMP_ATTR
2577
2578
2579  void CDomain::solveInheritanceTransformation()
2580  TRY
2581  {
2582    if (solveInheritanceTransformation_done_) return;
2583    else solveInheritanceTransformation_done_=true ;
2584
2585    CDomain* domain = this;
2586    CDomain* Lastdomain ;
2587    std::list<CDomain*> refDomains;
2588    bool out=false ;
2589    vector<StdString> excludedAttr;
2590    excludedAttr.push_back("domain_ref");
2591   
2592    refDomains.push_front(domain) ;
2593    while (domain->hasDirectDomainReference() && !out)
2594    {
2595      CDomain* lastDomain=domain ;
2596      domain = domain->getDirectDomainReference();
2597      domain->solveRefInheritance() ;
2598      if (!domain->SuperClass::isEqual(lastDomain,excludedAttr)) out=true ;
2599      refDomains.push_front(domain) ;
2600    }
2601
2602    CTransformationPaths::TPath path ;
2603    auto& pathList = std::get<2>(path) ;
2604    std::get<0>(path) = EElement::DOMAIN ;
2605    std::get<1>(path) = refDomains.front()->getId() ;
2606    for (auto& domain : refDomains)
2607    {
2608      CDomain::TransMapTypes transformations = domain->getAllTransformations();
2609      for(auto& transformation : transformations) pathList.push_back({transformation.second->getTransformationType(), 
2610                                                                      transformation.second->getId()}) ;
2611    }
2612    transformationPaths_.addPath(path) ;
2613
2614  }
2615  CATCH_DUMP_ATTR
2616 
2617
2618  bool CDomain::activateFieldWorkflow(CGarbageCollector& gc)
2619  TRY
2620  {
2621    if (!domain_ref.isEmpty())
2622    {
2623      CField* field=getFieldFromId(domain_ref) ;
2624      if (field!=nullptr)
2625      {
2626        bool ret = field->buildWorkflowGraph(gc) ;
2627        if (!ret) return false ; // cannot build workflow graph at this state
2628      }
2629      else 
2630      {
2631        CDomain* domain = get(domain_ref) ;
2632        bool ret = domain->activateFieldWorkflow(gc) ;
2633        if (!ret) return false ; // cannot build workflow graph at this state
2634        domain_ref=domain->getId() ; // replace domain_ref by solved reference
2635      }
2636    }
2637    activateFieldWorkflow_done_=true ;
2638    return true ;
2639  }
2640  CATCH_DUMP_ATTR
2641/////////////////////////////////////////////////////////////////////////////////////////////
2642/////////////////////////////////////////////////////////////////////////////////////////////
2643
2644  void CDomain::setContextClient(CContextClient* contextClient)
2645  TRY
2646  {
2647    if (clientsSet.find(contextClient)==clientsSet.end())
2648    {
2649      clients.push_back(contextClient) ;
2650      clientsSet.insert(contextClient);
2651    }
2652  }
2653  CATCH_DUMP_ATTR
2654
2655  /*!
2656    Parse children nodes of a domain in xml file.
2657    Whenver there is a new transformation, its type and name should be added into this function
2658    \param node child node to process
2659  */
2660  void CDomain::parse(xml::CXMLNode & node)
2661  TRY
2662  {
2663    SuperClass::parse(node);
2664
2665    if (node.goToChildElement())
2666    {
2667      StdString nodeElementName;
2668      do
2669      {
2670        StdString nodeId("");
2671        if (node.getAttributes().end() != node.getAttributes().find("id"))
2672        { nodeId = node.getAttributes()["id"]; }
2673
2674        nodeElementName = node.getElementName();
2675        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_.end(), it;
2676        it = transformationMapList_.find(nodeElementName);
2677        if (ite != it)
2678        {
2679          if (it->first == "zoom_domain")
2680          {
2681            info(0) << "WARNING : " << it->first << " is deprecated, replaced by extract_domain." << endl;
2682          }
2683          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
2684                                                                                                                nodeId,
2685                                                                                                                &node)));
2686        }
2687        else
2688        {
2689          ERROR("void CDomain::parse(xml::CXMLNode & node)",
2690                << "The transformation " << nodeElementName << " has not been supported yet.");
2691        }
2692      } while (node.goToNextElement()) ;
2693      node.goToParentElement();
2694    }
2695  }
2696  CATCH_DUMP_ATTR
2697   //----------------------------------------------------------------
2698
2699   DEFINE_REF_FUNC(Domain,domain)
2700
2701   ///---------------------------------------------------------------
2702
2703} // namespace xios
Note: See TracBrowser for help on using the repository browser.