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

Last change on this file since 2606 was 2606, checked in by jderouillat, 5 months ago

Compute the GridTransformConnector? (and the associated GridRemoteConnector?) only in a distributed context for the identification of templated elements

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