source: XIOS/dev/branch_openmp/src/node/domain.cpp @ 1460

Last change on this file since 1460 was 1460, checked in by yushan, 6 years ago

branch_openmp merged with XIOS_DEV_CMIP6@1459

  • 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: 114.1 KB
RevLine 
[219]1#include "domain.hpp"
2
[352]3#include "attribute_template.hpp"
4#include "object_template.hpp"
5#include "group_template.hpp"
[219]6
[591]7#include "xios_spl.hpp"
[300]8#include "event_client.hpp"
9#include "event_server.hpp"
10#include "buffer_in.hpp"
[352]11#include "message.hpp"
12#include "type.hpp"
13#include "context.hpp"
14#include "context_client.hpp"
[676]15#include "context_server.hpp"
[369]16#include "array_new.hpp"
[676]17#include "distribution_client.hpp"
[553]18#include "server_distribution_description.hpp"
[569]19#include "client_server_mapping_distributed.hpp"
[274]20
[1460]21using namespace ep_lib;
22
[676]23#include <algorithm>
24
[335]25namespace xios {
[509]26
[1328]27   /// ////////////////////// Définitions ////////////////////// ///
[219]28
29   CDomain::CDomain(void)
30      : CObjectTemplate<CDomain>(), CDomainAttributes()
[1460]31      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_()
32      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
[665]33      , isClientAfterTransformationChecked(false), hasLonLat(false)
[1460]34      , isRedistributed_(false), hasPole(false), doZoomByIndex_(false)
35      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
36      , globalLocalIndexMap_(), computedWrittenIndex_(false)
37          , clients()
[878]38   {
[924]39   }
[219]40
41   CDomain::CDomain(const StdString & id)
42      : CObjectTemplate<CDomain>(id), CDomainAttributes()
[1460]43      , isChecked(false), relFiles(), isClientChecked(false), nbSenders(), indSrv_(), connectedServerRank_() 
44      , hasBounds(false), hasArea(false), isCompressible_(false), isUnstructed_(false)
[665]45      , isClientAfterTransformationChecked(false), hasLonLat(false)
[1460]46      , isRedistributed_(false), hasPole(false), doZoomByIndex_(false)
47      , lonvalue(), latvalue(), bounds_lonvalue(), bounds_latvalue()
48      , globalLocalIndexMap_(), computedWrittenIndex_(false)
49          , clients()
[878]50   {
[1460]51    }
[219]52
53   CDomain::~CDomain(void)
[509]54   {
[219]55   }
56
57   ///---------------------------------------------------------------
58
[924]59   void CDomain::assignMesh(const StdString meshName, const int nvertex)
[881]60   {
[924]61     mesh = CMesh::getMesh(meshName, nvertex);
[881]62   }
63
[622]64   CDomain* CDomain::createDomain()
65   {
66     CDomain* domain = CDomainGroup::get("domain_definition")->createChild();
67     return domain;
68   }
69
[1334]70
[1328]71   std::map<StdString, ETranformationType> *CDomain::transformationMapList_ptr = 0;
[836]72
73   bool CDomain::initializeTransformationMap(std::map<StdString, ETranformationType>& m)
74   {
75     m["zoom_domain"] = TRANS_ZOOM_DOMAIN;
76     m["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
77     m["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
[934]78     m["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
[935]79     m["expand_domain"] = TRANS_EXPAND_DOMAIN;
[1460]80     m["reorder_domain"] = TRANS_REORDER_DOMAIN;
[836]81   }
82
[1134]83   bool CDomain::initializeTransformationMap()
84   {
85     CDomain::transformationMapList_ptr = new std::map<StdString, ETranformationType>();
86     (*CDomain::transformationMapList_ptr)["zoom_domain"] = TRANS_ZOOM_DOMAIN;
87     (*CDomain::transformationMapList_ptr)["interpolate_domain"] = TRANS_INTERPOLATE_DOMAIN;
88     (*CDomain::transformationMapList_ptr)["generate_rectilinear_domain"] = TRANS_GENERATE_RECTILINEAR_DOMAIN;
89     (*CDomain::transformationMapList_ptr)["compute_connectivity_domain"] = TRANS_COMPUTE_CONNECTIVITY_DOMAIN;
90     (*CDomain::transformationMapList_ptr)["expand_domain"] = TRANS_EXPAND_DOMAIN;
[1460]91     (*CDomain::transformationMapList_ptr)["reorder_domain"] = TRANS_REORDER_DOMAIN;
[1134]92   }
93
[219]94   const std::set<StdString> & CDomain::getRelFiles(void) const
95   {
[509]96      return (this->relFiles);
[219]97   }
98
[895]99
[676]100   /*!
101     Returns the number of indexes written by each server.
102     \return the number of indexes written by each server
103   */
[1460]104   int CDomain::getNumberWrittenIndexes(ep_lib::MPI_Comm writtenCom)
[676]105   {
[1460]106     int writtenSize;
107     ep_lib::MPI_Comm_size(writtenCom, &writtenSize);
108     return numberWrittenIndexes_[writtenSize];
[676]109   }
110
111   /*!
112     Returns the total number of indexes written by the servers.
113     \return the total number of indexes written by the servers
114   */
[1460]115   int CDomain::getTotalNumberWrittenIndexes(ep_lib::MPI_Comm writtenCom)
[676]116   {
[1460]117     int writtenSize;
118     ep_lib::MPI_Comm_size(writtenCom, &writtenSize);
119     return totalNumberWrittenIndexes_[writtenSize];
[676]120   }
121
122   /*!
123     Returns the offset of indexes written by each server.
124     \return the offset of indexes written by each server
125   */
[1460]126   int CDomain::getOffsetWrittenIndexes(ep_lib::MPI_Comm writtenCom)
[676]127   {
[1460]128     int writtenSize;
129     ep_lib::MPI_Comm_size(writtenCom, &writtenSize);
130     return offsetWrittenIndexes_[writtenSize];
[676]131   }
132
[1460]133   CArray<int, 1>& CDomain::getCompressedIndexToWriteOnServer(ep_lib::MPI_Comm writtenCom)
134   {
135     int writtenSize;
136     ep_lib::MPI_Comm_size(writtenCom, &writtenSize);
137     return compressedIndexToWriteOnServer[writtenSize];
138   }
139
[676]140   //----------------------------------------------------------------
141
[731]142   /*!
143    * Compute the minimum buffer size required to send the attributes to the server(s).
144    *
145    * \return A map associating the server rank with its minimum buffer size.
146    */
[1460]147   std::map<int, StdSize> CDomain::getAttributesBufferSize(CContextClient* client, bool bufferForWriting /*= false*/)
[731]148   {
149
[1460]150     std::map<int, StdSize> attributesSizes = getMinimumBufferSizeForAttributes(client);
[731]151
152     if (client->isServerLeader())
153     {
[1460]154       // size estimation for sendDistributionAttribut
[731]155       size_t size = 11 * sizeof(size_t);
156
157       const std::list<int>& ranks = client->getRanksServerLeader();
158       for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
159       {
160         if (size > attributesSizes[*itRank])
161           attributesSizes[*itRank] = size;
162       }
163     }
164
[1460]165     boost::unordered_map<int, vector<size_t> >::const_iterator itIndexEnd = indSrv_[client->serverSize].end();
166     // std::map<int, std::vector<int> >::const_iterator itWrittenIndexEnd = indWrittenSrv_.end();
167     for (size_t k = 0; k < connectedServerRank_[client->serverSize].size(); ++k)
[731]168     {
[1460]169       int rank = connectedServerRank_[client->serverSize][k];
170       boost::unordered_map<int, std::vector<size_t> >::const_iterator it = indSrv_[client->serverSize].find(rank);
[764]171       size_t idxCount = (it != itIndexEnd) ? it->second.size() : 0;
172
[731]173       // size estimation for sendIndex (and sendArea which is always smaller or equal)
[764]174       size_t sizeIndexEvent = 2 * sizeof(size_t) + 2 * CArray<int,1>::size(idxCount);
[1460]175       // if (isCompressible_)
176       // {
177       //   std::map<int, std::vector<int> >::const_iterator itWritten = indWrittenSrv_.find(rank);
178       //   size_t writtenIdxCount = (itWritten != itWrittenIndexEnd) ? itWritten->second.size() : 0;
179       //   sizeIndexEvent += CArray<int,1>::size(writtenIdxCount);
180       // }
[731]181
182       // size estimation for sendLonLat
[764]183       size_t sizeLonLatEvent = CArray<double,1>::size(idxCount);
[731]184       if (hasBounds)
[764]185         sizeLonLatEvent += CArray<double,2>::size(nvertex * idxCount);
[731]186
187       size_t size = CEventClient::headerSize + getId().size() + sizeof(size_t) + std::max(sizeIndexEvent, sizeLonLatEvent);
[764]188       if (size > attributesSizes[rank])
189         attributesSizes[rank] = size;
[731]190     }
191
192     return attributesSizes;
193   }
194
195   //----------------------------------------------------------------
196
[219]197   bool CDomain::isEmpty(void) const
198   {
[1460]199      return ((this->zoom_i_index.isEmpty()) || (0 == this->zoom_i_index.numElements()));
200
[219]201   }
202
203   //----------------------------------------------------------------
[676]204
[219]205   bool CDomain::IsWritten(const StdString & filename) const
206   {
207      return (this->relFiles.find(filename) != this->relFiles.end());
208   }
209
[676]210   bool CDomain::isWrittenCompressed(const StdString& filename) const
211   {
212      return (this->relFilesCompressed.find(filename) != this->relFilesCompressed.end());
213   }
214
[219]215   //----------------------------------------------------------------
[676]216
[594]217   bool CDomain::isDistributed(void) const
218   {
[1460]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      distributed |= (1 == CContext::getCurrent()->client->clientSize);
222
223      return distributed;
[594]224   }
225
226   //----------------------------------------------------------------
[676]227
228   /*!
229    * Test whether the data defined on the domain can be outputted in a compressed way.
[679]230    *
[676]231    * \return true if and only if a mask was defined for this domain
232    */
233   bool CDomain::isCompressible(void) const
234   {
235      return isCompressible_;
236   }
237
[219]238   void CDomain::addRelFile(const StdString & filename)
239   {
240      this->relFiles.insert(filename);
241   }
[895]242
[676]243   void CDomain::addRelFileCompressed(const StdString& filename)
244   {
245      this->relFilesCompressed.insert(filename);
246   }
[219]247
248   StdString CDomain::GetName(void)   { return (StdString("domain")); }
249   StdString CDomain::GetDefName(void){ return (CDomain::GetName()); }
250   ENodeType CDomain::GetType(void)   { return (eDomain); }
251
252   //----------------------------------------------------------------
253
[687]254   /*!
[1064]255      Verify if all distribution information of a domain are available
256      This checking verifies the definition of distribution attributes (ni, nj, ibegin, jbegin)
257   */
258   bool CDomain::distributionAttributesHaveValue() const
259   {
260      bool hasValues = true;
261
262      if (ni.isEmpty() && ibegin.isEmpty() && i_index.isEmpty())
263      {
264        hasValues = false;
265        return hasValues;
266      }
267
268      return hasValues;
269   }
270
271   /*!
[1460]272     Redistribute RECTILINEAR or CURVILINEAR domain with a number of local domains.
[687]273   All attributes ni,nj,ibegin,jbegin (if defined) will be rewritten
274   The optional attributes lonvalue, latvalue will be added. Because this function only serves (for now)
275   for interpolation from unstructured domain to rectilinear one, range of latvalue is 0-360 and lonvalue is -90 - +90
276    \param [in] nbLocalDomain number of local domain on the domain destination
277   */
278   void CDomain::redistribute(int nbLocalDomain)
279   {
[715]280     if (this->isRedistributed_) return;
[727]281
[783]282     this->isRedistributed_ = true;
283     CContext* context = CContext::getCurrent();
[1460]284     // For now the assumption is that secondary server pools consist of the same number of procs.
285     // CHANGE the line below if the assumption changes.
286     CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[0] : context->client;
[783]287     int rankClient = client->clientRank;
288     int rankOnDomain = rankClient%nbLocalDomain;
[687]289
[783]290     if (ni_glo.isEmpty() || ni_glo <= 0 )
291     {
292        ERROR("CDomain::redistribute(int nbLocalDomain)",
293           << "[ Id = " << this->getId() << " ] "
294           << "The global domain is badly defined,"
295           << " check the \'ni_glo\'  value !")
296     }
[687]297
[783]298     if (nj_glo.isEmpty() || nj_glo <= 0 )
299     {
300        ERROR("CDomain::redistribute(int nbLocalDomain)",
301           << "[ Id = " << this->getId() << " ] "
302           << "The global domain is badly defined,"
303           << " check the \'nj_glo\'  value !")
304     }
[687]305
[783]306     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
307     {
[687]308        int globalDomainSize = ni_glo * nj_glo;
309        if (globalDomainSize <= nbLocalDomain)
310        {
311          for (int idx = 0; idx < nbLocalDomain; ++idx)
312          {
313            if (rankOnDomain < globalDomainSize)
314            {
315              int iIdx = rankOnDomain % ni_glo;
316              int jIdx = rankOnDomain / ni_glo;
317              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
318              ni.setValue(1); nj.setValue(1);
319            }
320            else
321            {
322              ibegin.setValue(0); jbegin.setValue(0);
323              ni.setValue(0); nj.setValue(0);
324            }
325          }
326        }
327        else
328        {
[688]329          float njGlo = nj_glo.getValue();
330          float niGlo = ni_glo.getValue();
331          int nbProcOnX, nbProcOnY, range;
332
[687]333          // Compute (approximately) number of segment on x and y axis
[688]334          float yOverXRatio = njGlo/niGlo;
335
[687]336          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio));
337          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX);
338
339          // Simple distribution: Sweep from top to bottom, left to right
340          // Calculate local begin on x
341          std::vector<int> ibeginVec(nbProcOnX,0), jbeginVec(nbProcOnY,0);
342          std::vector<int> niVec(nbProcOnX), njVec(nbProcOnY);
343          for (int i = 1; i < nbProcOnX; ++i)
344          {
345            range = ni_glo / nbProcOnX;
346            if (i < (ni_glo%nbProcOnX)) ++range;
347            niVec[i-1] = range;
348            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
349          }
350          niVec[nbProcOnX-1] = ni_glo - ibeginVec[nbProcOnX-1];
351
352          // Calculate local begin on y
353          for (int j = 1; j < nbProcOnY; ++j)
354          {
355            range = nj_glo / nbProcOnY;
356            if (j < (nj_glo%nbProcOnY)) ++range;
357            njVec[j-1] = range;
358            jbeginVec[j] = jbeginVec[j-1] + njVec[j-1];
359          }
360          njVec[nbProcOnY-1] = nj_glo - jbeginVec[nbProcOnY-1];
361
362          // Now assign value to ni, ibegin, nj, jbegin
363          int iIdx = rankOnDomain % nbProcOnX;
364          int jIdx = rankOnDomain / nbProcOnX;
365
366          if (rankOnDomain != (nbLocalDomain-1))
367          {
368            ibegin.setValue(ibeginVec[iIdx]);
369            jbegin.setValue(jbeginVec[jIdx]);
370            nj.setValue(njVec[jIdx]);
371            ni.setValue(niVec[iIdx]);
372          }
373          else // just merge all the remaining rectangle into the last one
374          {
375            ibegin.setValue(ibeginVec[iIdx]);
376            jbegin.setValue(jbeginVec[jIdx]);
377            nj.setValue(njVec[jIdx]);
378            ni.setValue(ni_glo - ibeginVec[iIdx]);
379          }
[1064]380        } 
[687]381     }
[783]382     else  // unstructured domain
383     {
[785]384       if (this->i_index.isEmpty())
385       {
386          int globalDomainSize = ni_glo * nj_glo;
387          if (globalDomainSize <= nbLocalDomain)
[783]388          {
[785]389            for (int idx = 0; idx < nbLocalDomain; ++idx)
[783]390            {
[785]391              if (rankOnDomain < globalDomainSize)
392              {
393                int iIdx = rankOnDomain % ni_glo;
394                int jIdx = rankOnDomain / ni_glo;
395                ibegin.setValue(iIdx); jbegin.setValue(jIdx);
396                ni.setValue(1); nj.setValue(1);
397              }
398              else
399              {
400                ibegin.setValue(0); jbegin.setValue(0);
401                ni.setValue(0); nj.setValue(0);
402              }
[783]403            }
[785]404          }
405          else
406          {
407            float njGlo = nj_glo.getValue();
408            float niGlo = ni_glo.getValue();
409            std::vector<int> ibeginVec(nbLocalDomain,0);
410            std::vector<int> niVec(nbLocalDomain);
411            for (int i = 1; i < nbLocalDomain; ++i)
[783]412            {
[785]413              int range = ni_glo / nbLocalDomain;
414              if (i < (ni_glo%nbLocalDomain)) ++range;
415              niVec[i-1] = range;
416              ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
[783]417            }
[785]418            niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
419
420            int iIdx = rankOnDomain % nbLocalDomain;
421            ibegin.setValue(ibeginVec[iIdx]);
422            jbegin.setValue(0);
423            ni.setValue(niVec[iIdx]);
424            nj.setValue(1);
[783]425          }
[1064]426
427          i_index.resize(ni);         
428          for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
[783]429        }
430        else
431        {
[785]432          ibegin.setValue(this->i_index(0));
[783]433          jbegin.setValue(0);
[785]434          ni.setValue(this->i_index.numElements());
[783]435          nj.setValue(1);
436        }
437     }
[785]438
439     checkDomain();
[687]440   }
441
442   /*!
[1064]443     Fill in longitude and latitude whose values are read from file
444   */
445   void CDomain::fillInLonLat()
446   {
447     switch (type)
448     {
449      case type_attr::rectilinear:
450        fillInRectilinearLonLat();
451        break;
452      case type_attr::curvilinear:
453        fillInCurvilinearLonLat();
454        break;
455      case type_attr::unstructured:
456        fillInUnstructuredLonLat();
457        break;
458
459      default:
460      break;
461     }
[1328]462     completeLonLatClient() ;
[1064]463
464   }
465
466   /*!
[687]467     Fill in the values for lonvalue_1d and latvalue_1d of rectilinear domain
468     Range of longitude value from 0 - 360
469     Range of latitude value from -90 - +90
470   */
471   void CDomain::fillInRectilinearLonLat()
472   {
[775]473     if (!lonvalue_rectilinear_read_from_file.isEmpty())
474     {
475       lonvalue_1d.resize(ni);
476       for (int idx = 0; idx < ni; ++idx)
477         lonvalue_1d(idx) = lonvalue_rectilinear_read_from_file(idx+ibegin);
478       lon_start.setValue(lonvalue_rectilinear_read_from_file(0));
479       lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1));
480     }
481     else
482     {
483       if (!lonvalue_2d.isEmpty()) lonvalue_2d.free();
484       lonvalue_1d.resize(ni);
485       double lonRange = lon_end - lon_start;
486       double lonStep = (1 == ni_glo.getValue()) ? lonRange : lonRange/double(ni_glo.getValue()-1);
[688]487
[775]488        // Assign lon value
489       for (int i = 0; i < ni; ++i)
490       {
491         if (0 == (ibegin + i))
492         {
493           lonvalue_1d(i) = lon_start;
494         }
495         else if (ni_glo == (ibegin + i + 1))
496         {
497           lonvalue_1d(i) = lon_end;
498         }
499         else
500         {
501           lonvalue_1d(i) = (ibegin + i) * lonStep  + lon_start;
502         }
503       }
504     }
[687]505
[727]506
[775]507     if (!latvalue_rectilinear_read_from_file.isEmpty())
[687]508     {
[775]509       latvalue_1d.resize(nj);
510       for (int idx = 0; idx < nj; ++idx)
511         latvalue_1d(idx) = latvalue_rectilinear_read_from_file(idx+jbegin);
512       lat_start.setValue(latvalue_rectilinear_read_from_file(0));
513       lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1));
[687]514     }
[775]515     else
516     {
517       if (!latvalue_2d.isEmpty()) latvalue_1d.free();
518       latvalue_1d.resize(nj);
[687]519
[775]520       double latRange = lat_end - lat_start;
521       double latStep = (1 == nj_glo.getValue()) ? latRange : latRange/double(nj_glo.getValue()-1);
522
523       for (int j = 0; j < nj; ++j)
[727]524       {
[775]525         if (0 == (jbegin + j))
526         {
527            latvalue_1d(j) = lat_start;
528         }
529         else if (nj_glo == (jbegin + j + 1))
530         {
531            latvalue_1d(j) = lat_end;
532         }
533         else
534         {
535           latvalue_1d(j) =  (jbegin + j) * latStep + lat_start;
536         }
[727]537       }
[687]538     }
539   }
540
[1064]541    /*
542      Fill in longitude and latitude of curvilinear domain read from a file
543      If there are already longitude and latitude defined by model. We just igonore reading value.
544    */
545   void CDomain::fillInCurvilinearLonLat()
546   {
547     if (!lonvalue_curvilinear_read_from_file.isEmpty() && lonvalue_2d.isEmpty())
548     {
549       lonvalue_2d.resize(ni,nj);
550       for (int jdx = 0; jdx < nj; ++jdx)
551        for (int idx = 0; idx < ni; ++idx)
[1460]552         lonvalue_2d(idx,jdx) = lonvalue_curvilinear_read_from_file(idx, jdx);
[809]553
[1064]554       lonvalue_curvilinear_read_from_file.free();
555     }
[809]556
[1064]557     if (!latvalue_curvilinear_read_from_file.isEmpty() && latvalue_2d.isEmpty())
558     {
559       latvalue_2d.resize(ni,nj);
560       for (int jdx = 0; jdx < nj; ++jdx)
561        for (int idx = 0; idx < ni; ++idx)
[1460]562           latvalue_2d(idx,jdx) = latvalue_curvilinear_read_from_file(idx, jdx);
[1064]563
564       latvalue_curvilinear_read_from_file.free();
565     }
566
567     if (!bounds_lonvalue_curvilinear_read_from_file.isEmpty() && bounds_lon_2d.isEmpty())
568     {
569       bounds_lon_2d.resize(nvertex,ni,nj);
570       for (int jdx = 0; jdx < nj; ++jdx)
571        for (int idx = 0; idx < ni; ++idx)
572          for (int ndx = 0; ndx < nvertex; ++ndx)
[1460]573            bounds_lon_2d(ndx,idx,jdx) = bounds_lonvalue_curvilinear_read_from_file(ndx,idx, jdx);
[1064]574
575       bounds_lonvalue_curvilinear_read_from_file.free();
576     }
577
578     if (!bounds_latvalue_curvilinear_read_from_file.isEmpty() && bounds_lat_2d.isEmpty())
579     {
580       bounds_lat_2d.resize(nvertex,ni,nj);
581       for (int jdx = 0; jdx < nj; ++jdx)
582        for (int idx = 0; idx < ni; ++idx)
583          for (int ndx = 0; ndx < nvertex; ++ndx)
[1460]584            bounds_lat_2d(ndx,idx,jdx) = bounds_latvalue_curvilinear_read_from_file(ndx,idx, jdx);
[1064]585
586       bounds_latvalue_curvilinear_read_from_file.free();
587     }
588
589   }
590
591    /*
592      Fill in longitude and latitude of unstructured domain read from a file
593      If there are already longitude and latitude defined by model. We just igonore reading value.
594    */
595   void CDomain::fillInUnstructuredLonLat()
596   {
597     if (i_index.isEmpty())
598     {
599       i_index.resize(ni);
600       for(int idx = 0; idx < ni; ++idx) i_index(idx)=ibegin+idx;
601     }
602
603     if (!lonvalue_unstructured_read_from_file.isEmpty() && lonvalue_1d.isEmpty())
604     {
605        lonvalue_1d.resize(ni);
606        for (int idx = 0; idx < ni; ++idx)
[1460]607          lonvalue_1d(idx) = lonvalue_unstructured_read_from_file(idx);
[1064]608
609        // We dont need these values anymore, so just delete them
610        lonvalue_unstructured_read_from_file.free();
611     } 
612
613     if (!latvalue_unstructured_read_from_file.isEmpty() && latvalue_1d.isEmpty())
614     {
615        latvalue_1d.resize(ni);
616        for (int idx = 0; idx < ni; ++idx)
[1460]617          latvalue_1d(idx) =  latvalue_unstructured_read_from_file(idx);
[1064]618
619        // We dont need these values anymore, so just delete them
620        latvalue_unstructured_read_from_file.free();
621     }
622
623     if (!bounds_lonvalue_unstructured_read_from_file.isEmpty() && bounds_lon_1d.isEmpty())
624     {
625        int nbVertex = nvertex;
626        bounds_lon_1d.resize(nbVertex,ni);
627        for (int idx = 0; idx < ni; ++idx)
628          for (int jdx = 0; jdx < nbVertex; ++jdx)
[1460]629            bounds_lon_1d(jdx,idx) = bounds_lonvalue_unstructured_read_from_file(jdx, idx);
[1064]630
631        // We dont need these values anymore, so just delete them
[1460]632        bounds_lonvalue_unstructured_read_from_file.free();
[1064]633     }
634
635     if (!bounds_latvalue_unstructured_read_from_file.isEmpty() && bounds_lat_1d.isEmpty())
636     {
637        int nbVertex = nvertex;
638        bounds_lat_1d.resize(nbVertex,ni);
639        for (int idx = 0; idx < ni; ++idx)
640          for (int jdx = 0; jdx < nbVertex; ++jdx)
[1460]641            bounds_lat_1d(jdx,idx) = bounds_latvalue_unstructured_read_from_file(jdx, idx);
[1064]642
643        // We dont need these values anymore, so just delete them
[1460]644        bounds_latvalue_unstructured_read_from_file.free();
[1064]645     }
646   }
647
648  /*
649    Get global longitude and latitude of rectilinear domain.
650  */
[809]651   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g)
[688]652   {
[1460]653     CContext* context = CContext::getCurrent();
654    // For now the assumption is that secondary server pools consist of the same number of procs.
655    // CHANGE the line below if the assumption changes.
656    CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[0] : context->client;
657     lon_g.resize(ni_glo) ;
658     lat_g.resize(nj_glo) ;
[815]659
660
[1460]661     int* ibegin_g = new int[client->clientSize] ;
662     int* jbegin_g = new int[client->clientSize] ;
663     int* ni_g = new int[client->clientSize] ;
664     int* nj_g = new int[client->clientSize] ;
665     int v ;
666     v=ibegin ;
667     ep_lib::MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,client->intraComm) ;
668     v=jbegin ;
669     ep_lib::MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,client->intraComm) ;
670     v=ni ;
671     ep_lib::MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,client->intraComm) ;
672     v=nj ;
673     ep_lib::MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,client->intraComm) ;
[815]674
[1460]675     ep_lib::MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,client->intraComm) ;
676     ep_lib::MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,client->intraComm) ;
[815]677
[809]678      delete[] ibegin_g ;
679      delete[] jbegin_g ;
680      delete[] ni_g ;
681      delete[] nj_g ;
682   }
683
684   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
685                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
686   {
[688]687     int i,j,k;
[815]688
[688]689     const int nvertexValue = 4;
[775]690     boundsLon.resize(nvertexValue,ni*nj);
[688]691
[899]692     if (ni_glo>1)
693     {
694       double lonStepStart = lon(1)-lon(0);
695       bounds_lon_start=lon(0) - lonStepStart/2;
696       double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2);
697       bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
698       double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
[810]699
[899]700       // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
701       if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
702       {
703         bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
704         bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
705       }
706     }
707     else
[810]708     {
[899]709       if (bounds_lon_start.isEmpty()) bounds_lon_start=-180. ;
710       if (bounds_lon_end.isEmpty()) bounds_lon_end=180.-1e-8 ;
[810]711     }
[906]712
[809]713     for(j=0;j<nj;++j)
714       for(i=0;i<ni;++i)
715       {
716         k=j*ni+i;
717         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start
718                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2;
719         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end
720                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2;
721       }
[727]722
[815]723
[809]724    boundsLat.resize(nvertexValue,nj*ni);
[810]725    bool isNorthPole=false ;
726    bool isSouthPole=false ;
[809]727    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
728    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
[808]729
[810]730    // lat boundaries beyond pole the assimilate it to pole
731    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
[899]732    if (nj_glo>1)
[810]733    {
[899]734      double latStepStart = lat(1)-lat(0);
735      if (isNorthPole) bounds_lat_start=lat(0);
736      else
[810]737      {
[899]738        bounds_lat_start=lat(0)-latStepStart/2;
739        if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
740        else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
741        else if (bounds_lat_start <= 90 && bounds_lat_start >= lat(0))
742        {
743          if ( std::abs(90-bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
744        }
745        else if (bounds_lat_start >= -90 && bounds_lat_start <= lat(0))
746        {
747          if ( std::abs(-90 - bounds_lat_start) <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
748        }
[810]749      }
[899]750
751      double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
752      if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
753      else
[810]754      {
[899]755        bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
756
757        if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
758        else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
759        else if (bounds_lat_end <= 90 && bounds_lat_end >= lat(nj_glo-1))
760        {
761          if ( std::abs(90-bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
762        }
763        else if (bounds_lat_end >= -90 && bounds_lat_end <= lat(nj_glo-1))
764        {
765          if ( std::abs(-90 - bounds_lat_end) <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
766        }
[810]767      }
768    }
769    else
770    {
[1062]771      if (bounds_lat_start.isEmpty()) bounds_lat_start=-90. ;
[899]772      if (bounds_lat_end.isEmpty()) bounds_lat_end=90 ;
[815]773    }
774
[809]775    for(j=0;j<nj;++j)
776      for(i=0;i<ni;++i)
777      {
778        k=j*ni+i;
779        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start
780                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2;
781        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end
782                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2;
783      }
[688]784   }
785
[1460]786   /*
787     General check of the domain to verify its mandatory attributes
788   */
[467]789   void CDomain::checkDomain(void)
[219]790   {
[664]791     if (type.isEmpty())
792     {
793       ERROR("CDomain::checkDomain(void)",
[679]794             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
795             << "The domain type is mandatory, "
796             << "please define the 'type' attribute.")
[664]797     }
[969]798
[953]799     if (type == type_attr::gaussian) 
800     {
[1460]801        hasPole=true ;
802        type.setValue(type_attr::unstructured) ;
803      }
804      else if (type == type_attr::rectilinear) hasPole=true ;
805
[679]806     if (type == type_attr::unstructured)
[664]807     {
[679]808        if (ni_glo.isEmpty())
[664]809        {
[679]810          ERROR("CDomain::checkDomain(void)",
811                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
812                << "The global domain is badly defined, "
813                << "the mandatory 'ni_glo' attribute is missing.")
[664]814        }
[679]815        else if (ni_glo <= 0)
816        {
817          ERROR("CDomain::checkDomain(void)",
818                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
819                << "The global domain is badly defined, "
820                << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
821        }
[664]822        isUnstructed_ = true;
823        nj_glo = 1;
824        nj = 1;
825        jbegin = 0;
[942]826        if (!i_index.isEmpty()) ni = i_index.numElements();
[664]827        j_index.resize(ni);
828        for(int i=0;i<ni;++i) j_index(i)=0;
[611]829
[664]830        if (!area.isEmpty())
831          area.transposeSelf(1, 0);
832     }
[679]833
834     if (ni_glo.isEmpty())
[664]835     {
[679]836       ERROR("CDomain::checkDomain(void)",
837             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
838             << "The global domain is badly defined, "
839             << "the mandatory 'ni_glo' attribute is missing.")
[664]840     }
[679]841     else if (ni_glo <= 0)
842     {
843       ERROR("CDomain::checkDomain(void)",
844             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
845             << "The global domain is badly defined, "
846             << "'ni_glo' attribute should be strictly positive so 'ni_glo = " << ni_glo.getValue() << "' is invalid.")
847     }
[509]848
[679]849     if (nj_glo.isEmpty())
850     {
851       ERROR("CDomain::checkDomain(void)",
852             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
853             << "The global domain is badly defined, "
854             << "the mandatory 'nj_glo' attribute is missing.")
855     }
856     else if (nj_glo <= 0)
857     {
858       ERROR("CDomain::checkDomain(void)",
859             << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
860             << "The global domain is badly defined, "
861             << "'nj_glo' attribute should be strictly positive so 'nj_glo = " << nj_glo.getValue() << "' is invalid.")
862     }
863
[664]864     checkLocalIDomain();
865     checkLocalJDomain();
[509]866
[664]867     if (i_index.isEmpty())
868     {
869       i_index.resize(ni*nj);
870       for (int j = 0; j < nj; ++j)
871         for (int i = 0; i < ni; ++i) i_index(i+j*ni) = i+ibegin;
872     }
[594]873
[664]874     if (j_index.isEmpty())
875     {
876       j_index.resize(ni*nj);
877       for (int j = 0; j < nj; ++j)
878         for (int i = 0; i < ni; ++i) j_index(i+j*ni) = j+jbegin;
879     }
[1078]880     
[821]881     checkZoom();
882   }
[663]883
[969]884   // Check global zoom of a domain
885   // If there is no zoom defined for the domain, zoom will have value of global doamin
[821]886   void CDomain::checkZoom(void)
887   {
888     if (global_zoom_ibegin.isEmpty())
889      global_zoom_ibegin.setValue(0);
890     if (global_zoom_ni.isEmpty())
891      global_zoom_ni.setValue(ni_glo);
892     if (global_zoom_jbegin.isEmpty())
893      global_zoom_jbegin.setValue(0);
894     if (global_zoom_nj.isEmpty())
895      global_zoom_nj.setValue(nj_glo);
[1460]896    if (zoom_i_index.isEmpty()) zoom_i_index.setValue(i_index.getValue());
897    if (zoom_j_index.isEmpty()) zoom_j_index.setValue(j_index.getValue());
898    if (zoom_ibegin.isEmpty()) zoom_ibegin.setValue(ibegin);
899    if (zoom_ni.isEmpty()) zoom_ni.setValue(ni);
900    if (zoom_jbegin.isEmpty()) zoom_jbegin.setValue(jbegin);
901    if (zoom_nj.isEmpty()) zoom_nj.setValue(nj);
[219]902   }
903
[1460]904   size_t CDomain::getGlobalWrittenSize(void)
905   {
906      return global_zoom_ni*global_zoom_nj ;
907   }
[219]908   //----------------------------------------------------------------
909
[969]910   // Check validity of local domain on using the combination of 3 parameters: ibegin, ni and i_index
[219]911   void CDomain::checkLocalIDomain(void)
912   {
[969]913      // If ibegin and ni are provided then we use them to check the validity of local domain
914      if (i_index.isEmpty() && !ibegin.isEmpty() && !ni.isEmpty())
[219]915      {
[969]916        if ((ni.getValue() < 0 || ibegin.getValue() < 0) || ((ibegin.getValue() + ni.getValue()) > ni_glo.getValue()))
917        {
918          ERROR("CDomain::checkLocalIDomain(void)",
919                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
920                << "The local domain is wrongly defined,"
921                << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
922        }
923      }
924
925      // i_index has higher priority than ibegin and ni
926      if (!i_index.isEmpty())
927      {
[975]928        int minIIndex = (0 < i_index.numElements()) ? i_index(0) : 0;
[969]929        if (ni.isEmpty()) 
930        {         
931         // No information about ni
932          int minIndex = ni_glo - 1;
933          int maxIndex = 0;
934          for (int idx = 0; idx < i_index.numElements(); ++idx)
935          {
936            if (i_index(idx) < minIndex) minIndex = i_index(idx);
937            if (i_index(idx) > maxIndex) maxIndex = i_index(idx);
938          }
939          ni = maxIndex - minIndex + 1; 
940          minIIndex = minIIndex;         
941        }
942
943        // It's not so correct but if ibegin is not the first value of i_index
944        // then data on local domain has user-defined distribution. In this case, ibegin, ni have no meaning.
945        if (ibegin.isEmpty()) ibegin = minIIndex;
946      }
947      else if (ibegin.isEmpty() && ni.isEmpty())
948      {
[594]949        ibegin = 0;
950        ni = ni_glo;
951      }
[971]952      else if ((!ibegin.isEmpty() && ni.isEmpty()) || (ibegin.isEmpty() && !ni.isEmpty()))
[219]953      {
[969]954        ERROR("CDomain::checkLocalIDomain(void)",
955              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
[971]956              << "The local domain is wrongly defined," << endl
957              << "i_index is empty and either 'ni' or 'ibegin' is not defined. " 
958              << "If 'ni' and 'ibegin' are used to define a domain, both of them must not be empty.");
[219]959      }
[969]960       
[219]961
[969]962      if ((ni.getValue() < 0 || ibegin.getValue() < 0))
[594]963      {
964        ERROR("CDomain::checkLocalIDomain(void)",
[679]965              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
[594]966              << "The local domain is wrongly defined,"
[679]967              << " check the attributes 'ni_glo' (" << ni_glo.getValue() << "), 'ni' (" << ni.getValue() << ") and 'ibegin' (" << ibegin.getValue() << ")");
[594]968      }
[219]969   }
970
[969]971   // Check validity of local domain on using the combination of 3 parameters: jbegin, nj and j_index
[219]972   void CDomain::checkLocalJDomain(void)
973   {
[969]974    // If jbegin and nj are provided then we use them to check the validity of local domain
975     if (j_index.isEmpty() && !jbegin.isEmpty() && !nj.isEmpty())
[657]976     {
[969]977       if ((nj.getValue() < 0 || jbegin.getValue() < 0) || (jbegin.getValue() + nj.getValue()) > nj_glo.getValue())
978       {
979         ERROR("CDomain::checkLocalJDomain(void)",
980                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
981                << "The local domain is wrongly defined,"
982                << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
983       }
[657]984     }
[969]985
986     if (!j_index.isEmpty())
[657]987     {
[975]988        int minJIndex = (0 < j_index.numElements()) ? j_index(0) : 0;
[969]989        if (nj.isEmpty()) 
990        {
991          // No information about nj
992          int minIndex = nj_glo - 1;
993          int maxIndex = 0;
994          for (int idx = 0; idx < j_index.numElements(); ++idx)
995          {
996            if (j_index(idx) < minIndex) minIndex = j_index(idx);
997            if (j_index(idx) > maxIndex) maxIndex = j_index(idx);
998          }
999          nj = maxIndex - minIndex + 1;
1000          minJIndex = minIndex; 
1001        } 
1002        // It's the same as checkLocalIDomain. It's not so correct but if jbegin is not the first value of j_index
1003        // then data on local domain has user-defined distribution. In this case, jbegin has no meaning.
1004       if (jbegin.isEmpty()) jbegin = minJIndex;       
[657]1005     }
[969]1006     else if (jbegin.isEmpty() && nj.isEmpty())
1007     {
1008       jbegin = 0;
1009       nj = nj_glo;
1010     }     
[219]1011
[969]1012
1013     if ((nj.getValue() < 0 || jbegin.getValue() < 0))
1014     {
1015       ERROR("CDomain::checkLocalJDomain(void)",
[679]1016              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
[594]1017              << "The local domain is wrongly defined,"
[679]1018              << " check the attributes 'nj_glo' (" << nj_glo.getValue() << "), 'nj' (" << nj.getValue() << ") and 'jbegin' (" << jbegin.getValue() << ")");
[969]1019     }
[219]1020   }
1021
1022   //----------------------------------------------------------------
[679]1023
[219]1024   void CDomain::checkMask(void)
1025   {
[664]1026      if (!mask_1d.isEmpty() && !mask_2d.isEmpty())
[657]1027        ERROR("CDomain::checkMask(void)",
[679]1028              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1029              << "Both mask_1d and mask_2d are defined but only one can be used at the same time." << std::endl
1030              << "Please define only one mask: 'mask_1d' or 'mask_2d'.");
[657]1031
[664]1032      if (!mask_1d.isEmpty() && mask_2d.isEmpty())
[219]1033      {
[664]1034        if (mask_1d.numElements() != i_index.numElements())
[657]1035          ERROR("CDomain::checkMask(void)",
[679]1036                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1037                << "'mask_1d' does not have the same size as the local domain." << std::endl
1038                << "Local size is " << i_index.numElements() << "." << std::endl
1039                << "Mask size is " << mask_1d.numElements() << ".");
[657]1040      }
[509]1041
[664]1042      if (mask_1d.isEmpty() && !mask_2d.isEmpty())
[657]1043      {
[679]1044        if (mask_2d.extent(0) != ni || mask_2d.extent(1) != nj)
1045          ERROR("CDomain::checkMask(void)",
1046                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1047                << "The mask does not have the same size as the local domain." << std::endl
1048                << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1049                << "Mask size is " << mask_2d.extent(0) << " x " << mask_2d.extent(1) << ".");
[657]1050      }
[509]1051
[664]1052      if (!mask_2d.isEmpty())
[657]1053      {
[1460]1054        domainMask.resize(mask_2d.extent(0) * mask_2d.extent(1));
[657]1055        for (int j = 0; j < nj; ++j)
[1460]1056          for (int i = 0; i < ni; ++i) domainMask(i+j*ni) = mask_2d(i,j);
1057//        mask_2d.reset();
[657]1058      }
[675]1059      else if (mask_1d.isEmpty())
[657]1060      {
[1460]1061        domainMask.resize(i_index.numElements());
1062        for (int i = 0; i < i_index.numElements(); ++i) domainMask(i) = true;
[657]1063      }
[1460]1064      else
1065      {
1066      domainMask.resize(mask_1d.numElements());
1067      domainMask=mask_1d ;
1068     }
[219]1069   }
1070
[676]1071   //----------------------------------------------------------------
[219]1072
1073   void CDomain::checkDomainData(void)
[509]1074   {
[679]1075      if (data_dim.isEmpty())
[219]1076      {
[679]1077        data_dim.setValue(1);
[219]1078      }
[679]1079      else if (!(data_dim.getValue() == 1 || data_dim.getValue() == 2))
[219]1080      {
[679]1081        ERROR("CDomain::checkDomainData(void)",
1082              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1083              << "The data dimension is invalid, 'data_dim' must be 1 or 2 not << " << data_dim.getValue() << ".");
[219]1084      }
1085
1086      if (data_ibegin.isEmpty())
[679]1087         data_ibegin.setValue(0);
[660]1088      if (data_jbegin.isEmpty())
[679]1089         data_jbegin.setValue(0);
[219]1090
[679]1091      if (data_ni.isEmpty())
[219]1092      {
[679]1093        data_ni.setValue((data_dim == 1) ? (ni.getValue() * nj.getValue()) : ni.getValue());
[219]1094      }
[679]1095      else if (data_ni.getValue() < 0)
[219]1096      {
[679]1097        ERROR("CDomain::checkDomainData(void)",
1098              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1099              << "The data size cannot be negative ('data_ni' = " << data_ni.getValue() << ").");
[219]1100      }
1101
[679]1102      if (data_nj.isEmpty())
[219]1103      {
[679]1104        data_nj.setValue((data_dim.getValue() == 1) ? (ni.getValue() * nj.getValue()) : nj.getValue());
[219]1105      }
[679]1106      else if (data_nj.getValue() < 0)
1107      {
1108        ERROR("CDomain::checkDomainData(void)",
1109              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1110              << "The data size cannot be negative ('data_nj' = " << data_nj.getValue() << ").");
1111      }
[219]1112   }
1113
1114   //----------------------------------------------------------------
[679]1115
[219]1116   void CDomain::checkCompression(void)
1117   {
1118      if (!data_i_index.isEmpty())
1119      {
[664]1120        if (!data_j_index.isEmpty() &&
[678]1121            data_j_index.numElements() != data_i_index.numElements())
[664]1122        {
[679]1123           ERROR("CDomain::checkCompression(void)",
1124                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1125                 << "'data_i_index' and 'data_j_index' arrays must have the same size." << std::endl
1126                 << "'data_i_index' size = " << data_i_index.numElements() << std::endl
1127                 << "'data_j_index' size = " << data_j_index.numElements());
[664]1128        }
1129
[679]1130        if (2 == data_dim)
1131        {
[664]1132          if (data_j_index.isEmpty())
[660]1133          {
[679]1134             ERROR("CDomain::checkCompression(void)",
1135                   << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1136                   << "'data_j_index' must be defined when 'data_i_index' is set and 'data_dim' is 2.");
[660]1137          }
[679]1138        }
1139        else // (1 == data_dim)
1140        {
[664]1141          if (data_j_index.isEmpty())
1142          {
[679]1143            data_j_index.resize(data_ni);
1144            for (int j = 0; j < data_ni; ++j) data_j_index(j) = 0;
[664]1145          }
[679]1146        }
[219]1147      }
1148      else
1149      {
[679]1150        if (data_dim == 2 && !data_j_index.isEmpty())
1151          ERROR("CDomain::checkCompression(void)",
1152                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1153                << "'data_i_index' must be defined when 'data_j_index' is set and 'data_dim' is 2.");
[219]1154
[679]1155        if (1 == data_dim)
1156        {
1157          data_i_index.resize(data_ni);
1158          data_j_index.resize(data_ni);
[666]1159
[679]1160          for (int i = 0; i < data_ni; ++i)
1161          {
1162            data_i_index(i) = i;
1163            data_j_index(i) = 0;
1164          }
1165        }
1166        else // (data_dim == 2)
1167        {
1168          const int dsize = data_ni * data_nj;
1169          data_i_index.resize(dsize);
1170          data_j_index.resize(dsize);
[509]1171
[679]1172          for(int count = 0, j = 0; j < data_nj; ++j)
1173          {
1174            for(int i = 0; i < data_ni; ++i, ++count)
[219]1175            {
[679]1176              data_i_index(count) = i;
1177              data_j_index(count) = j;
[219]1178            }
[679]1179          }
1180        }
[219]1181      }
1182   }
1183
[666]1184   //----------------------------------------------------------------
[893]1185   void CDomain::computeLocalMask(void)
1186   {
[1460]1187     localMask.resize(i_index.numElements()) ;
[893]1188     localMask=false ;
[676]1189
[893]1190     size_t dn=data_i_index.numElements() ;
1191     int i,j ;
1192     size_t k,ind ;
[906]1193
[893]1194     for(k=0;k<dn;k++)
1195     {
1196       if (data_dim==2)
1197       {
1198          i=data_i_index(k)+data_ibegin ;
1199          j=data_j_index(k)+data_jbegin ;
[1460]1200          if (i>=0 && i<ni && j>=0 && j<nj)
1201          {
1202            ind=j*ni+i ;
1203            localMask(ind)=domainMask(ind) ;
1204          }
[893]1205       }
1206       else
1207       {
[1460]1208          i=data_i_index(k)+data_ibegin ;
1209          if (i>=0 && i<i_index.numElements())
1210          {
1211            ind=i ;
1212            localMask(ind)=domainMask(ind) ;
1213          }
[893]1214       }
1215     }
1216   }
1217
[676]1218   void CDomain::checkEligibilityForCompressedOutput(void)
1219   {
1220     // We don't check if the mask or the indexes are valid here, just if they have been defined at this point.
1221     isCompressible_ = !mask_1d.isEmpty() || !mask_2d.isEmpty() || !data_i_index.isEmpty();
1222   }
1223
1224   //----------------------------------------------------------------
1225
[1460]1226   /*
1227     Fill in longitude and latitude value from clients (or models) into internal values lonvalue, latvalue which
1228     will be used by XIOS.
1229   */
[666]1230   void CDomain::completeLonLatClient(void)
[219]1231   {
[1460]1232     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1233     checkBounds() ;
1234     if (!lonvalue_2d.isEmpty() && !lonlatValueExisted)
[664]1235     {
[1460]1236       lonvalue.resize(ni * nj);
1237       latvalue.resize(ni * nj);
[761]1238       if (hasBounds)
1239       {
[1460]1240         bounds_lonvalue.resize(nvertex, ni * nj);
1241         bounds_latvalue.resize(nvertex, ni * nj);
[761]1242       }
1243
1244       for (int j = 0; j < nj; ++j)
1245       {
1246         for (int i = 0; i < ni; ++i)
1247         {
1248           int k = j * ni + i;
1249
[1460]1250           lonvalue(k) = lonvalue_2d(i,j);
1251           latvalue(k) = latvalue_2d(i,j);
[761]1252
1253           if (hasBounds)
1254           {
1255             for (int n = 0; n < nvertex; ++n)
1256             {
[1460]1257               bounds_lonvalue(n,k) = bounds_lon_2d(n,i,j);
1258               bounds_latvalue(n,k) = bounds_lat_2d(n,i,j);
[761]1259             }
1260           }
1261         }
1262       }
[664]1263     }
[1460]1264     else if (!lonvalue_1d.isEmpty()  && !lonlatValueExisted)
[664]1265     {
1266       if (type_attr::rectilinear == type)
1267       {
[679]1268         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
[664]1269         {
[1460]1270           lonvalue.resize(ni * nj);
1271           latvalue.resize(ni * nj);
[761]1272           if (hasBounds)
1273           {
[1460]1274             bounds_lonvalue.resize(nvertex, ni * nj);
1275             bounds_latvalue.resize(nvertex, ni * nj);
[761]1276           }
1277
1278           for (int j = 0; j < nj; ++j)
1279           {
1280             for (int i = 0; i < ni; ++i)
[664]1281             {
[761]1282               int k = j * ni + i;
1283
[1460]1284               lonvalue(k) = lonvalue_1d(i);
1285               latvalue(k) = latvalue_1d(j);
[761]1286
[664]1287               if (hasBounds)
1288               {
[679]1289                 for (int n = 0; n < nvertex; ++n)
[666]1290                 {
[1460]1291                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1292                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
[666]1293                 }
[664]1294               }
1295             }
[761]1296           }
1297         }
[1460]1298         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
[974]1299         {
[1460]1300           lonvalue.reference(lonvalue_1d);
1301           latvalue.reference(latvalue_1d);
[974]1302            if (hasBounds)
1303           {
[1460]1304             bounds_lonvalue.reference(bounds_lon_1d);
1305             bounds_latvalue.reference(bounds_lat_1d);
[974]1306           }
1307         }
[761]1308         else
1309           ERROR("CDomain::completeLonClient(void)",
1310                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1311                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
[974]1312                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1313                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1314                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1315                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
[664]1316       }
[1460]1317       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
[664]1318       {
[1460]1319         lonvalue.reference(lonvalue_1d);
1320         latvalue.reference(latvalue_1d);
[664]1321         if (hasBounds)
1322         {
[1460]1323           bounds_lonvalue.reference(bounds_lon_1d);
1324           bounds_latvalue.reference(bounds_lat_1d);
[664]1325         }
1326       }
1327     }
1328   }
1329
[1460]1330   /*
1331     Convert internal longitude latitude value used by XIOS to "lonvalue_*" which can be retrieved with Fortran interface
1332   */
1333   void CDomain::convertLonLatValue(void)
1334   {
1335     bool lonlatValueExisted = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1336     if (!lonvalue_2d.isEmpty() && lonlatValueExisted)
1337     {
1338       lonvalue_2d.resize(ni,nj);
1339       latvalue_2d.resize(ni,nj);
1340       if (hasBounds)
1341       {
1342         bounds_lon_2d.resize(nvertex, ni, nj);
1343         bounds_lat_2d.resize(nvertex, ni, nj);
1344       }
1345
1346       for (int j = 0; j < nj; ++j)
1347       {
1348         for (int i = 0; i < ni; ++i)
1349         {
1350           int k = j * ni + i;
1351
1352           lonvalue_2d(i,j) = lonvalue(k);
1353           latvalue_2d(i,j) = latvalue(k);
1354
1355           if (hasBounds)
1356           {
1357             for (int n = 0; n < nvertex; ++n)
1358             {
1359               bounds_lon_2d(n,i,j) = bounds_lonvalue(n,k);
1360               bounds_lat_2d(n,i,j) = bounds_latvalue(n,k);
1361             }
1362           }
1363         }
1364       }
1365     }
1366     else if (!lonvalue_1d.isEmpty()  && lonlatValueExisted)
1367     {
1368       if (type_attr::rectilinear == type)
1369       {
1370         if (ni == lonvalue_1d.numElements() && nj == latvalue_1d.numElements())
1371         {
1372           lonvalue.resize(ni * nj);
1373           latvalue.resize(ni * nj);
1374           if (hasBounds)
1375           {
1376             bounds_lonvalue.resize(nvertex, ni * nj);
1377             bounds_latvalue.resize(nvertex, ni * nj);
1378           }
1379
1380           for (int j = 0; j < nj; ++j)
1381           {
1382             for (int i = 0; i < ni; ++i)
1383             {
1384               int k = j * ni + i;
1385
1386               lonvalue(k) = lonvalue_1d(i);
1387               latvalue(k) = latvalue_1d(j);
1388
1389               if (hasBounds)
1390               {
1391                 for (int n = 0; n < nvertex; ++n)
1392                 {
1393                   bounds_lonvalue(n,k) = bounds_lon_1d(n,i);
1394                   bounds_latvalue(n,k) = bounds_lat_1d(n,j);
1395                 }
1396               }
1397             }
1398           }
1399         }
1400         else if (i_index.numElements() == lonvalue_1d.numElements() && j_index.numElements() == latvalue_1d.numElements()  && !lonlatValueExisted)
1401         {
1402           lonvalue.reference(lonvalue_1d);
1403           latvalue.reference(latvalue_1d);
1404            if (hasBounds)
1405           {
1406             bounds_lonvalue.reference(bounds_lon_1d);
1407             bounds_latvalue.reference(bounds_lat_1d);
1408           }
1409         }
1410         else
1411           ERROR("CDomain::completeLonClient(void)",
1412                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1413                 << "'lonvalue_1d' and 'latvalue_1d' does not have the same size as the local domain." << std::endl
1414                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() 
1415                 << " and 'latvalue_1d' size is " << latvalue_1d.numElements() << std::endl
1416                 << " They should be correspondingly " << ni.getValue() << " and "  << nj.getValue() << " or " << std::endl
1417                 << i_index.numElements() << " and "  << j_index.numElements() << ".");
1418       }
1419       else if (type == type_attr::curvilinear || type == type_attr::unstructured  && !lonlatValueExisted)
1420       {
1421         lonvalue.reference(lonvalue_1d);
1422         latvalue.reference(latvalue_1d);
1423         if (hasBounds)
1424         {
1425           bounds_lonvalue.reference(bounds_lon_1d);
1426           bounds_latvalue.reference(bounds_lat_1d);
1427         }
1428       }
1429     }
1430   }
1431
1432
[449]1433   void CDomain::checkBounds(void)
1434   {
[1460]1435     bool hasBoundValues = (0 != bounds_lonvalue.numElements()) || (0 != bounds_latvalue.numElements());
1436     if (!nvertex.isEmpty() && nvertex > 0 && !hasBoundValues)
[449]1437     {
[664]1438       if (!bounds_lon_1d.isEmpty() && !bounds_lon_2d.isEmpty())
1439         ERROR("CDomain::checkBounds(void)",
[679]1440               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1441               << "Only one longitude boundary attribute can be used but both 'bounds_lon_1d' and 'bounds_lon_2d' are defined." << std::endl
1442               << "Define only one longitude boundary attribute: 'bounds_lon_1d' or 'bounds_lon_2d'.");
[664]1443
1444       if (!bounds_lat_1d.isEmpty() && !bounds_lat_2d.isEmpty())
1445         ERROR("CDomain::checkBounds(void)",
[679]1446               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1447               << "Only one latitude boundary attribute can be used but both 'bounds_lat_1d' and 'bounds_lat_2d' are defined." << std::endl
1448               << "Define only one latitude boundary attribute: 'bounds_lat_1d' or 'bounds_lat_2d'.");
[664]1449
1450       if ((!bounds_lon_1d.isEmpty() && bounds_lat_1d.isEmpty()) || (bounds_lon_1d.isEmpty() && !bounds_lat_1d.isEmpty()))
1451       {
1452         ERROR("CDomain::checkBounds(void)",
[679]1453               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1454               << "Only 'bounds_lon_1d' or 'bounds_lat_1d' is defined." << std::endl
1455               << "Please define either both attributes or none.");
[664]1456       }
1457
1458       if ((!bounds_lon_2d.isEmpty() && bounds_lat_2d.isEmpty()) || (bounds_lon_2d.isEmpty() && !bounds_lat_2d.isEmpty()))
1459       {
1460         ERROR("CDomain::checkBounds(void)",
[679]1461               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1462               << "Only 'bounds_lon_2d' or 'bounds_lat_2d' is defined." << std::endl
1463               << "Please define either both attributes or none.");
[664]1464       }
1465
[691]1466       if (!bounds_lon_1d.isEmpty() && nvertex.getValue() != bounds_lon_1d.extent(0))
[679]1467         ERROR("CDomain::checkBounds(void)",
1468               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1469               << "'bounds_lon_1d' dimension is not compatible with 'nvertex'." << std::endl
1470               << "'bounds_lon_1d' dimension is " << bounds_lon_1d.extent(1)
1471               << " but nvertex is " << nvertex.getValue() << ".");
[664]1472
[691]1473       if (!bounds_lon_2d.isEmpty() && nvertex.getValue() != bounds_lon_2d.extent(0))
[679]1474         ERROR("CDomain::checkBounds(void)",
1475               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1476               << "'bounds_lon_2d' dimension is not compatible with 'nvertex'." << std::endl
1477               << "'bounds_lon_2d' dimension is " << bounds_lon_2d.extent(2)
1478               << " but nvertex is " << nvertex.getValue() << ".");
[664]1479
1480       if (!bounds_lon_1d.isEmpty() && lonvalue_1d.isEmpty())
[679]1481         ERROR("CDomain::checkBounds(void)",
1482               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1483               << "Since 'bounds_lon_1d' is defined, 'lonvalue_1d' must be defined too." << std::endl);
[664]1484
1485       if (!bounds_lon_2d.isEmpty() && lonvalue_2d.isEmpty())
[679]1486         ERROR("CDomain::checkBounds(void)",
1487               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1488               << "Since 'bounds_lon_2d' is defined, 'lonvalue_2d' must be defined too." << std::endl);
[664]1489
[691]1490       if (!bounds_lat_1d.isEmpty() && nvertex.getValue() != bounds_lat_1d.extent(0))
[679]1491         ERROR("CDomain::checkBounds(void)",
1492               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1493               << "'bounds_lat_1d' dimension is not compatible with 'nvertex'." << std::endl
1494               << "'bounds_lat_1d' dimension is " << bounds_lat_1d.extent(1)
1495               << " but nvertex is " << nvertex.getValue() << ".");
[664]1496
[691]1497       if (!bounds_lat_2d.isEmpty() && nvertex.getValue() != bounds_lat_2d.extent(0))
[679]1498         ERROR("CDomain::checkBounds(void)",
1499               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1500               << "'bounds_lat_2d' dimension is not compatible with 'nvertex'." << std::endl
1501               << "'bounds_lat_2d' dimension is " << bounds_lat_2d.extent(2)
1502               << " but nvertex is " << nvertex.getValue() << ".");
[664]1503
1504       if (!bounds_lat_1d.isEmpty() && latvalue_1d.isEmpty())
[679]1505         ERROR("CDomain::checkBounds(void)",
1506               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1507               << "Since 'bounds_lat_1d' is defined, 'latvalue_1d' must be defined too." << std::endl);
[664]1508
1509       if (!bounds_lat_2d.isEmpty() && latvalue_2d.isEmpty())
[679]1510         ERROR("CDomain::checkBounds(void)",
1511               << "Since 'bounds_lat_2d' is defined, 'latvalue_2d' must be defined too." << std::endl);
1512
1513       hasBounds = true;
[449]1514     }
[1460]1515     else if (hasBoundValues)
1516     {
1517       hasBounds = true;       
1518     }
[509]1519     else
[449]1520     {
[679]1521       hasBounds = false;
[449]1522     }
1523   }
[509]1524
[611]1525   void CDomain::checkArea(void)
1526   {
[1460]1527     bool hasAreaValue = (!areavalue.isEmpty() && 0 != areavalue.numElements());
[611]1528     hasArea = !area.isEmpty();
[1460]1529     if (hasArea && !hasAreaValue)
[611]1530     {
1531       if (area.extent(0) != ni || area.extent(1) != nj)
1532       {
[679]1533         ERROR("CDomain::checkArea(void)",
1534               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1535               << "The area does not have the same size as the local domain." << std::endl
1536               << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1537               << "Area size is " << area.extent(0) << " x " << area.extent(1) << ".");
[611]1538       }
[1460]1539       if (areavalue.isEmpty())
1540       {
1541          areavalue.resize(ni*nj);
1542         for (int j = 0; j < nj; ++j)
1543         {
1544           for (int i = 0; i < ni; ++i)
1545           {
1546             int k = j * ni + i;
1547             areavalue(k) = area(i,j);
1548           }
1549         }
1550       }
[611]1551     }
1552   }
1553
[665]1554   void CDomain::checkLonLat()
1555   {
[1460]1556     if (!hasLonLat) hasLonLat = (!latvalue_1d.isEmpty() && !lonvalue_1d.isEmpty()) ||
1557                                 (!latvalue_2d.isEmpty() && !lonvalue_2d.isEmpty());
1558     bool hasLonLatValue = (0 != lonvalue.numElements()) || (0 != latvalue.numElements());
1559     if (hasLonLat && !hasLonLatValue)
[666]1560     {
1561       if (!lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
[762]1562         ERROR("CDomain::checkLonLat()",
[679]1563               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1564               << "Only one longitude attribute can be used but both 'lonvalue_1d' and 'lonvalue_2d' are defined." << std::endl
1565               << "Define only one longitude attribute: 'lonvalue_1d' or 'lonvalue_2d'.");
[666]1566
1567       if (!lonvalue_1d.isEmpty() && lonvalue_2d.isEmpty())
1568       {
[687]1569         if ((type_attr::rectilinear != type) && (lonvalue_1d.numElements() != i_index.numElements()))
[762]1570           ERROR("CDomain::checkLonLat()",
[679]1571                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1572                 << "'lonvalue_1d' does not have the same size as the local domain." << std::endl
1573                 << "Local size is " << i_index.numElements() << "." << std::endl
1574                 << "'lonvalue_1d' size is " << lonvalue_1d.numElements() << ".");
[666]1575       }
1576
1577       if (lonvalue_1d.isEmpty() && !lonvalue_2d.isEmpty())
1578       {
[679]1579         if (lonvalue_2d.extent(0) != ni || lonvalue_2d.extent(1) != nj)
[762]1580           ERROR("CDomain::checkLonLat()",
[679]1581                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1582                 << "'lonvalue_2d' does not have the same size as the local domain." << std::endl
1583                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1584                 << "'lonvalue_2d' size is " << lonvalue_2d.extent(0) << " x " << lonvalue_2d.extent(1) << ".");
[666]1585       }
1586
1587       if (!latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
[762]1588         ERROR("CDomain::checkLonLat()",
[679]1589               << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1590               << "Only one latitude attribute can be used but both 'latvalue_1d' and 'latvalue_2d' are defined." << std::endl
1591               << "Define only one latitude attribute: 'latvalue_1d' or 'latvalue_2d'.");
[666]1592
1593       if (!latvalue_1d.isEmpty() && latvalue_2d.isEmpty())
1594       {
[687]1595         if ((type_attr::rectilinear != type) && (latvalue_1d.numElements() != i_index.numElements()))
[762]1596           ERROR("CDomain::checkLonLat()",
[679]1597                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1598                 << "'latvalue_1d' does not have the same size as the local domain." << std::endl
1599                 << "Local size is " << i_index.numElements() << "." << std::endl
1600                 << "'latvalue_1d' size is " << latvalue_1d.numElements() << ".");
[666]1601       }
1602
1603       if (latvalue_1d.isEmpty() && !latvalue_2d.isEmpty())
1604       {
[679]1605         if (latvalue_2d.extent(0) != ni || latvalue_2d.extent(1) != nj)
[762]1606           ERROR("CDomain::checkLonLat()",
[679]1607                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
1608                 << "'latvalue_2d' does not have the same size as the local domain." << std::endl
1609                 << "Local size is " << ni.getValue() << " x " << nj.getValue() << "." << std::endl
1610                 << "'latvalue_2d' size is " << latvalue_2d.extent(0) << " x " << latvalue_2d.extent(1) << ".");
[666]1611       }
1612     }
[665]1613   }
1614
[657]1615   void CDomain::checkAttributesOnClientAfterTransformation()
1616   {
1617     CContext* context=CContext::getCurrent() ;
1618
1619     if (this->isClientAfterTransformationChecked) return;
1620     if (context->hasClient)
1621     {
[1460]1622      this->computeConnectedClients();
1623       if (hasLonLat)
1624         if (!context->hasServer)
1625           this->completeLonLatClient();
[657]1626     }
1627
1628     this->isClientAfterTransformationChecked = true;
1629   }
1630
[219]1631   //----------------------------------------------------------------
[667]1632   // Divide function checkAttributes into 2 seperate ones
[509]1633   // This function only checks all attributes of current domain
1634   void CDomain::checkAttributesOnClient()
1635   {
1636     if (this->isClientChecked) return;
1637     CContext* context=CContext::getCurrent();
[219]1638
[1460]1639      if (context->hasClient && !context->hasServer)
1640      {
1641        this->checkDomain();
1642        this->checkBounds();
1643        this->checkArea();
1644        this->checkLonLat();
1645      }
[509]1646
[1460]1647      if (context->hasClient && !context->hasServer)
1648      { // Ct client uniquement
[509]1649         this->checkMask();
1650         this->checkDomainData();
1651         this->checkCompression();
[893]1652         this->computeLocalMask() ;
[509]1653      }
1654      else
[1460]1655      { // Ct serveur uniquement
[509]1656      }
1657
1658      this->isClientChecked = true;
1659   }
1660
1661   // Send all checked attributes to server
1662   void CDomain::sendCheckedAttributes()
1663   {
1664     if (!this->isClientChecked) checkAttributesOnClient();
[657]1665     if (!this->isClientAfterTransformationChecked) checkAttributesOnClientAfterTransformation();
[509]1666     CContext* context=CContext::getCurrent() ;
1667
1668     if (this->isChecked) return;
1669     if (context->hasClient)
1670     {
[1460]1671       sendAttributes();
[509]1672     }
1673     this->isChecked = true;
1674   }
1675
[219]1676   void CDomain::checkAttributes(void)
1677   {
1678      if (this->isChecked) return;
[347]1679      CContext* context=CContext::getCurrent() ;
[219]1680
[467]1681      this->checkDomain();
[665]1682      this->checkLonLat();
[449]1683      this->checkBounds();
[611]1684      this->checkArea();
[509]1685
[300]1686      if (context->hasClient)
[1460]1687      { // Ct client uniquement
[219]1688         this->checkMask();
1689         this->checkDomainData();
1690         this->checkCompression();
[893]1691         this->computeLocalMask() ;
[666]1692
[219]1693      }
1694      else
[1460]1695      { // Ct serveur uniquement
[219]1696      }
[509]1697
[300]1698      if (context->hasClient)
1699      {
[1460]1700        this->computeConnectedClients();
[666]1701        this->completeLonLatClient();
[300]1702      }
[509]1703
[219]1704      this->isChecked = true;
1705   }
[509]1706
[1460]1707  /*!
1708     Compute the connection of a client to other clients to determine which clients to send attributes to.
1709     The sending clients are supposed to already know the distribution of receiving clients (In simple cases, it's band)
1710     The connection among clients is calculated by using global index.
1711     A client connects to other clients which holds the same global index as it.     
1712  */
1713  void CDomain::computeConnectedClients()
[300]1714  {
[1460]1715    CContext* context=CContext::getCurrent() ;
1716   
1717    // This line should be changed soon.
1718    int nbSrvPools = (context->hasServer) ? (context->hasClient ? context->clientPrimServer.size() : 0) : 1;
[509]1719
[1460]1720    nbSenders.clear();
1721    connectedServerRank_.clear();
[657]1722
[1460]1723    for (int p = 0; p < nbSrvPools; ++p)
[595]1724    {
[1460]1725      CContextClient* client = (0 != context->clientPrimServer.size()) ? context->clientPrimServer[p] : context->client;
1726      int nbServer = client->serverSize;
1727      int nbClient = client->clientSize;
1728      int rank     = client->clientRank;
1729      bool doComputeGlobalIndexServer = true;
[509]1730
[1460]1731      if (connectedServerRank_.find(nbServer) == connectedServerRank_.end())
[595]1732      {
[509]1733
[1460]1734        if (indSrv_.find(nbServer) == indSrv_.end())
1735        {
1736          int i,j,i_ind,j_ind, nbIndex, nbIndexZoom;
1737          int global_zoom_iend=global_zoom_ibegin+global_zoom_ni-1;
1738          int global_zoom_jend=global_zoom_jbegin+global_zoom_nj-1;
[595]1739
[1460]1740           // Precompute number of index
1741           int globalIndexCountZoom = 0;
1742           nbIndex = i_index.numElements();
[467]1743
[1460]1744           if (doZoomByIndex_)
1745           {
1746             globalIndexCountZoom = zoom_i_index.numElements();
1747           }
1748           else
1749           {
1750             for (i = 0; i < nbIndex; ++i)
1751             {
1752               i_ind=i_index(i);
1753               j_ind=j_index(i);
[1078]1754
[1460]1755               if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1756               {
1757                  ++globalIndexCountZoom;
1758               }
1759             }
1760           }
[657]1761
[1460]1762           // Fill in index
1763           CArray<size_t,1> globalIndexDomainZoom(globalIndexCountZoom);
1764           CArray<size_t,1> localIndexDomainZoom(globalIndexCountZoom);
1765           CArray<size_t,1> globalIndexDomain(nbIndex);
1766           size_t globalIndex;
1767           int globalIndexCount = 0;
[449]1768
[1460]1769           for (i = 0; i < nbIndex; ++i)
1770           {
1771             i_ind=i_index(i);
1772             j_ind=j_index(i);
1773             globalIndex = i_ind + j_ind * ni_glo;
1774             globalIndexDomain(i) = globalIndex;
1775           }
[509]1776
[1460]1777           if (globalLocalIndexMap_.empty())
1778           {
1779             for (i = 0; i < nbIndex; ++i)
1780               globalLocalIndexMap_[globalIndexDomain(i)] = i;
1781           }
[657]1782
[1460]1783           globalIndexCountZoom = 0;
1784           if (doZoomByIndex_)
1785           {
1786             int nbIndexZoom = zoom_i_index.numElements();
[467]1787
[1460]1788             for (i = 0; i < nbIndexZoom; ++i)
1789             {
1790               i_ind=zoom_i_index(i);
1791               j_ind=zoom_j_index(i);
1792               globalIndex = i_ind + j_ind * ni_glo;
1793               globalIndexDomainZoom(globalIndexCountZoom) = globalIndex;
1794               ++globalIndexCountZoom;
1795             }
1796           }
1797           else
1798           {
1799             int global_zoom_iend=global_zoom_ibegin+global_zoom_ni-1;
1800             int global_zoom_jend=global_zoom_jbegin+global_zoom_nj-1;
1801             for (i = 0; i < nbIndex; ++i)
1802             {
1803               i_ind=i_index(i);
1804               j_ind=j_index(i);
1805               globalIndex = i_ind + j_ind * ni_glo;
1806               if (i_ind >= global_zoom_ibegin && i_ind <= global_zoom_iend && j_ind >= global_zoom_jbegin && j_ind <= global_zoom_jend)
1807               {
1808                  globalIndexDomainZoom(globalIndexCountZoom) = globalIndex;
1809                  ++globalIndexCountZoom;
1810               }
1811             }
1812
1813             int iend = ibegin + ni -1;
1814             int jend = jbegin + nj -1;
1815             zoom_ibegin = global_zoom_ibegin > ibegin ? global_zoom_ibegin : ibegin;
1816             int zoom_iend  = global_zoom_iend < iend ? zoom_iend : iend ;
1817             zoom_ni     = zoom_iend-zoom_ibegin+1 ;
1818
1819             zoom_jbegin = global_zoom_jbegin > jbegin ? global_zoom_jbegin : jbegin ;
1820             int zoom_jend   = global_zoom_jend < jend ? zoom_jend : jend;
1821             zoom_nj     = zoom_jend-zoom_jbegin+1;
1822           }
1823
1824           size_t globalSizeIndex = 1, indexBegin, indexEnd;
1825           int range, clientSize = client->clientSize;
1826           std::vector<int> nGlobDomain(2);
1827           nGlobDomain[0] = this->ni_glo;
1828           nGlobDomain[1] = this->nj_glo;
1829           for (int i = 0; i < nGlobDomain.size(); ++i) globalSizeIndex *= nGlobDomain[i];
1830           indexBegin = 0;
1831           if (globalSizeIndex <= clientSize)
1832           {
1833             indexBegin = rank%globalSizeIndex;
1834             indexEnd = indexBegin;
1835           }
1836           else
1837           {
1838             for (int i = 0; i < clientSize; ++i)
1839             {
1840               range = globalSizeIndex / clientSize;
1841               if (i < (globalSizeIndex%clientSize)) ++range;
1842               if (i == client->clientRank) break;
1843               indexBegin += range;
1844             }
1845             indexEnd = indexBegin + range - 1;
1846           }
1847
1848           // Even if servers have no index, they must received something from client
1849           // We only use several client to send "empty" message to these servers
1850           CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
1851           std::vector<int> serverZeroIndex;
1852           if (isUnstructed_) serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 0);
1853           else serverZeroIndex = serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 1);
1854
1855           std::list<int> serverZeroIndexLeader;
1856           std::list<int> serverZeroIndexNotLeader;
1857           CContextClient::computeLeader(client->clientRank, client->clientSize, serverZeroIndex.size(), serverZeroIndexLeader, serverZeroIndexNotLeader);
1858           for (std::list<int>::iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
1859              *it = serverZeroIndex[*it];
1860
1861           CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(),
1862                 client->intraComm);
1863           clientServerMap->computeServerIndexMapping(globalIndexDomain, nbServer);
1864           CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
1865
1866           CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
1867                 ite = globalIndexDomainOnServer.end();
1868           indSrv_[nbServer].swap(globalIndexDomainOnServer);
1869           connectedServerRank_[nbServer].clear();
1870           for (it = indSrv_[nbServer].begin(); it != ite; ++it)
1871             connectedServerRank_[nbServer].push_back(it->first);
1872
1873           for (std::list<int>::const_iterator it = serverZeroIndexLeader.begin(); it != serverZeroIndexLeader.end(); ++it)
1874              connectedServerRank_[nbServer].push_back(*it);
1875
1876           // Even if a client has no index, it must connect to at least one server and
1877           // send an "empty" data to this server
1878           if (connectedServerRank_[nbServer].empty())
1879              connectedServerRank_[nbServer].push_back(client->clientRank % client->serverSize);
1880
1881           nbSenders[nbServer] = clientServerMap->computeConnectedClients(client->serverSize, client->clientSize, client->intraComm, connectedServerRank_[nbServer]);
1882           delete clientServerMap;
[676]1883        }
1884      }
1885    }
[1460]1886  }
[676]1887
[1460]1888   /*!
1889     Compute index to write data. We only write data on the zoomed region, therefore, there should
1890     be a map between the complete grid and the reduced grid where we write data.
1891     By using global index we can easily create this kind of mapping.
1892   */
1893   void CDomain::computeWrittenIndex()
1894   { 
1895      if (computedWrittenIndex_) return;
1896      computedWrittenIndex_ = true;
[467]1897
[1460]1898      CContext* context=CContext::getCurrent();     
1899      CContextServer* server = context->server; 
[509]1900
[1460]1901      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
1902      nBegin[0]       = zoom_ibegin;  nBegin[1] = zoom_jbegin;
1903      nSize[0]        = zoom_ni;      nSize[1]  = zoom_nj;
1904      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
1905      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
1906      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
1907      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
1908
1909      size_t nbWritten = 0, indGlo;     
1910      boost::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
1911                                                          ite = globalLocalIndexMap_.end(), it;         
1912      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
1913                                       itSrve = writtenGlobalIndex.end(), itSrv;
1914
1915//      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1916//      {
1917//        indGlo = *itSrv;
1918//        if (ite != globalLocalIndexMap_.find(indGlo))
1919//        {
1920//          ++nbWritten;
1921//        }
1922//      }
1923
1924//      localIndexToWriteOnServer.resize(nbWritten);
1925      localIndexToWriteOnServer.resize(writtenGlobalIndex.numElements());
1926
1927      nbWritten = 0;
1928      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
[676]1929      {
[1460]1930        indGlo = *itSrv;
1931        if (ite != globalLocalIndexMap_.find(indGlo))
[676]1932        {
[1460]1933          localIndexToWriteOnServer(nbWritten) = globalLocalIndexMap_[indGlo];
1934          ++nbWritten;
[676]1935        }
[1460]1936        else
1937        {
1938          localIndexToWriteOnServer(nbWritten) = 0;
1939          ++nbWritten;
1940        }
[676]1941      }
[1460]1942     
1943      // if (isCompressible())
1944      // {
1945      //   nbWritten = 0;
1946      //   boost::unordered_map<size_t,size_t> localGlobalIndexMap;
1947      //   for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
1948      //   {
1949      //     indGlo = *itSrv;
1950      //     if (ite != globalLocalIndexMap_.find(indGlo))
1951      //     {
1952      //       localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
1953      //       ++nbWritten;
1954      //     }                 
1955      //   }
[569]1956
[1460]1957      //   nbWritten = 0;
1958      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1959      //   {
1960      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1961      //     {
1962      //       ++nbWritten;
1963      //     }
1964      //   }
[676]1965
[1460]1966      //   compressedIndexToWriteOnServer.resize(nbWritten);
1967      //   nbWritten = 0;
1968      //   for (int idx = 0; idx < data_i_index.numElements(); ++idx)
1969      //   {
1970      //     if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
1971      //     {
1972      //       compressedIndexToWriteOnServer(nbWritten) = localGlobalIndexMap[data_i_index(idx)];
1973      //       ++nbWritten;
1974      //     }
1975      //   }
[657]1976
[1460]1977      //   numberWrittenIndexes_ = nbWritten;
1978      //   if (isDistributed())
1979      //   {           
1980      //     ep_lib::MPI_Allreduce(&numberWrittenIndexes_, &totalNumberWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1981      //     ep_lib::MPI_Scan(&numberWrittenIndexes_, &offsetWrittenIndexes_, 1, MPI_INT, MPI_SUM, server->intraComm);
1982      //     offsetWrittenIndexes_ -= numberWrittenIndexes_;
1983      //   }
1984      //   else
1985      //     totalNumberWrittenIndexes_ = numberWrittenIndexes_;
1986      // }     
1987   }
[584]1988
[1460]1989  void CDomain::computeWrittenCompressedIndex(ep_lib::MPI_Comm writtenComm)
1990  {
1991    int writtenCommSize;
1992    ep_lib::MPI_Comm_size(writtenComm, &writtenCommSize);
1993    if (compressedIndexToWriteOnServer.find(writtenCommSize) != compressedIndexToWriteOnServer.end())
1994      return;
[668]1995
[1460]1996    if (isCompressible())
[584]1997    {
[1460]1998      size_t nbWritten = 0, indGlo;
1999      CContext* context=CContext::getCurrent();     
2000      CContextServer* server = context->server; 
2001
2002      std::vector<int> nBegin(2), nSize(2), nBeginGlobal(2), nGlob(2);
2003      nBegin[0]       = zoom_ibegin;  nBegin[1] = zoom_jbegin;
2004      nSize[0]        = zoom_ni;      nSize[1]  = zoom_nj;
2005      nBeginGlobal[0] = 0; nBeginGlobal[1] = 0;
2006      nGlob[0]        = ni_glo;   nGlob[1] = nj_glo;
2007      CDistributionServer srvDist(server->intraCommSize, nBegin, nSize, nBeginGlobal, nGlob); 
2008      const CArray<size_t,1>& writtenGlobalIndex  = srvDist.getGlobalIndex();
2009
2010      boost::unordered_map<size_t,size_t>::const_iterator itb = globalLocalIndexMap_.begin(),
2011                                                          ite = globalLocalIndexMap_.end(), it;   
2012      CArray<size_t,1>::const_iterator itSrvb = writtenGlobalIndex.begin(),
2013                                       itSrve = writtenGlobalIndex.end(), itSrv;
2014      boost::unordered_map<size_t,size_t> localGlobalIndexMap;
2015      for (itSrv = itSrvb; itSrv != itSrve; ++itSrv)
[584]2016      {
[1460]2017        indGlo = *itSrv;
2018        if (ite != globalLocalIndexMap_.find(indGlo))
[584]2019        {
[1460]2020          localGlobalIndexMap[localIndexToWriteOnServer(nbWritten)] = indGlo;
2021          ++nbWritten;
2022        }                 
2023      }
2024
2025      nbWritten = 0;
2026      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
2027      {
2028        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
2029        {
2030          ++nbWritten;
[584]2031        }
2032      }
[1460]2033
2034      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
2035      nbWritten = 0;
2036      for (int idx = 0; idx < data_i_index.numElements(); ++idx)
[676]2037      {
[1460]2038        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_i_index(idx)))
[676]2039        {
[1460]2040          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_i_index(idx)];
2041          ++nbWritten;
[676]2042        }
2043      }
[569]2044
[1460]2045      numberWrittenIndexes_[writtenCommSize] = nbWritten;
2046      if (isDistributed())
2047      {
2048             
2049        ep_lib::MPI_Allreduce(&numberWrittenIndexes_[writtenCommSize], &totalNumberWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2050        ep_lib::MPI_Scan(&numberWrittenIndexes_[writtenCommSize], &offsetWrittenIndexes_[writtenCommSize], 1, MPI_INT, MPI_SUM, writtenComm);
2051        offsetWrittenIndexes_[writtenCommSize] -= numberWrittenIndexes_[writtenCommSize];
2052      }
2053      else
2054        totalNumberWrittenIndexes_[writtenCommSize] = numberWrittenIndexes_[writtenCommSize];
2055      }
[300]2056  }
[467]2057
[1460]2058  /*!
2059    Send all attributes from client to connected clients
2060    The attributes will be rebuilt on receiving side
2061  */
2062  void CDomain::sendAttributes()
[657]2063  {
[1460]2064    sendDistributionAttributes();
2065    sendIndex();       
2066    sendMask();
2067    sendLonLat();
2068    sendArea();   
2069    sendDataIndex();
[657]2070  }
2071
[667]2072  /*!
[1460]2073    Send global index and zoom index from client to connected client(s)
2074    zoom index can be smaller than global index
[667]2075  */
[665]2076  void CDomain::sendIndex()
[300]2077  {
[610]2078    int ns, n, i, j, ind, nv, idx;
[1460]2079    std::list<CContextClient*>::iterator it;
2080    for (it=clients.begin(); it!=clients.end(); ++it)
2081    {
2082      CContextClient* client = *it;
[610]2083
[1460]2084      int serverSize = client->serverSize;
2085      CEventClient eventIndex(getType(), EVENT_ID_INDEX);
[665]2086
[1460]2087      list<CMessage> list_msgsIndex;
2088      list<CArray<int,1> > list_indZoom, list_writtenInd, list_indGlob;
[665]2089
[1460]2090      boost::unordered_map<int, vector<size_t> >::const_iterator itIndex, iteIndex;
2091      iteIndex = indSrv_[serverSize].end();
2092      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2093      {
2094        int nbIndGlob = 0;
2095        int rank = connectedServerRank_[serverSize][k];
2096        itIndex = indSrv_[serverSize].find(rank);
2097        if (iteIndex != itIndex)
2098          nbIndGlob = itIndex->second.size();
2099
2100        list_indGlob.push_back(CArray<int,1>(nbIndGlob));       
2101
2102        CArray<int,1>& indGlob = list_indGlob.back();
2103        for (n = 0; n < nbIndGlob; ++n)
2104        {
2105          indGlob(n) = static_cast<int>(itIndex->second[n]);
2106        }
2107
2108        list_msgsIndex.push_back(CMessage());
2109        list_msgsIndex.back() << this->getId() << (int)type; // enum ne fonctionne pour les message => ToFix
2110        list_msgsIndex.back() << isCurvilinear;
2111        list_msgsIndex.back() << list_indGlob.back(); //list_indi.back() << list_indj.back();
2112       
2113        eventIndex.push(rank, nbSenders[serverSize][rank], list_msgsIndex.back());
2114      }
2115
2116      client->sendEvent(eventIndex);
2117    }
2118  }
2119
2120  /*!
2121    Send distribution from client to other clients
2122    Because a client in a level knows correctly the grid distribution of client on the next level
2123    it calculates this distribution then sends it to the corresponding clients on the next level
2124  */
2125  void CDomain::sendDistributionAttributes(void)
2126  {
2127    std::list<CContextClient*>::iterator it;
2128    for (it=clients.begin(); it!=clients.end(); ++it)
[665]2129    {
[1460]2130      CContextClient* client = *it;
2131      int nbServer = client->serverSize;
2132      std::vector<int> nGlobDomain(2);
2133      nGlobDomain[0] = this->ni_glo;
2134      nGlobDomain[1] = this->nj_glo;
[665]2135
[1460]2136      CServerDistributionDescription serverDescription(nGlobDomain, nbServer);
2137      if (isUnstructed_) serverDescription.computeServerDistribution(false, 0);
2138      else serverDescription.computeServerDistribution(false, 1);
[665]2139
[1460]2140      std::vector<std::vector<int> > serverIndexBegin = serverDescription.getServerIndexBegin();
2141      std::vector<std::vector<int> > serverDimensionSizes = serverDescription.getServerDimensionSizes();
2142
2143      CEventClient event(getType(),EVENT_ID_SERVER_ATTRIBUT);
2144      if (client->isServerLeader())
[665]2145      {
[1460]2146        std::list<CMessage> msgs;
2147
2148        const std::list<int>& ranks = client->getRanksServerLeader();
2149        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
2150        {
2151          // Use const int to ensure CMessage holds a copy of the value instead of just a reference
2152          const int ibegin_srv = serverIndexBegin[*itRank][0];
2153          const int jbegin_srv = serverIndexBegin[*itRank][1];
2154          const int ni_srv = serverDimensionSizes[*itRank][0];
2155          const int nj_srv = serverDimensionSizes[*itRank][1];
2156
2157          msgs.push_back(CMessage());
2158          CMessage& msg = msgs.back();
2159          msg << this->getId() ;
2160          msg << isUnstructed_;
2161          msg << ni_srv << ibegin_srv << nj_srv << jbegin_srv;
2162          msg << global_zoom_ni.getValue() << global_zoom_ibegin.getValue() << global_zoom_nj.getValue() << global_zoom_jbegin.getValue();       
2163          msg << isCompressible_;
2164
2165          event.push(*itRank,1,msg);
2166        }
2167        client->sendEvent(event);
[665]2168      }
[1460]2169      else client->sendEvent(event);
2170    }
2171  }
[665]2172
[1460]2173  /*!
2174    Send mask index from client to connected(s) clients   
2175  */
2176  void CDomain::sendMask()
2177  {
2178    int ns, n, i, j, ind, nv, idx;
2179    std::list<CContextClient*>::iterator it;
2180    for (it=clients.begin(); it!=clients.end(); ++it)
2181    {
2182      CContextClient* client = *it;
2183      int serverSize = client->serverSize;
[665]2184
[1460]2185      // send area for each connected server
2186      CEventClient eventMask(getType(), EVENT_ID_MASK);
[665]2187
[1460]2188      list<CMessage> list_msgsMask;
2189      list<CArray<bool,1> > list_mask;
2190
2191      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2192      iteMap = indSrv_[serverSize].end();
2193      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
[676]2194      {
[1460]2195        int nbData = 0;
2196        int rank = connectedServerRank_[serverSize][k];
2197        it = indSrv_[serverSize].find(rank);
2198        if (iteMap != it)
2199          nbData = it->second.size();
2200        list_mask.push_back(CArray<bool,1>(nbData));
[676]2201
[1460]2202        const std::vector<size_t>& temp = it->second;
2203        for (n = 0; n < nbData; ++n)
2204        {
2205          idx = static_cast<int>(it->second[n]);
2206          list_mask.back()(n) = domainMask(globalLocalIndexMap_[idx]);
2207        }
[676]2208
[1460]2209        list_msgsMask.push_back(CMessage());
2210        list_msgsMask.back() << this->getId() << list_mask.back();
2211        eventMask.push(rank, nbSenders[serverSize][rank], list_msgsMask.back());
[676]2212      }
[1460]2213      client->sendEvent(eventMask);
[665]2214    }
2215  }
2216
[667]2217  /*!
[1460]2218    Send area from client to connected client(s)
[667]2219  */
[665]2220  void CDomain::sendArea()
2221  {
2222    if (!hasArea) return;
2223
2224    int ns, n, i, j, ind, nv, idx;
[1460]2225    std::list<CContextClient*>::iterator it;
[665]2226
[1460]2227    for (it=clients.begin(); it!=clients.end(); ++it)
2228    {
2229      CContextClient* client = *it;
2230      int serverSize = client->serverSize;
[665]2231
[1460]2232      // send area for each connected server
2233      CEventClient eventArea(getType(), EVENT_ID_AREA);
[665]2234
[1460]2235      list<CMessage> list_msgsArea;
2236      list<CArray<double,1> > list_area;
[665]2237
[1460]2238      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2239      iteMap = indSrv_[serverSize].end();
2240      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
[665]2241      {
[1460]2242        int nbData = 0;
2243        int rank = connectedServerRank_[serverSize][k];
2244        it = indSrv_[serverSize].find(rank);
2245        if (iteMap != it)
2246          nbData = it->second.size();
2247        list_area.push_back(CArray<double,1>(nbData));
2248
2249        const std::vector<size_t>& temp = it->second;
2250        for (n = 0; n < nbData; ++n)
2251        {
2252          idx = static_cast<int>(it->second[n]);
2253          list_area.back()(n) = areavalue(globalLocalIndexMap_[idx]);
2254        }
2255
2256        list_msgsArea.push_back(CMessage());
2257        list_msgsArea.back() << this->getId() << hasArea;
2258        list_msgsArea.back() << list_area.back();
2259        eventArea.push(rank, nbSenders[serverSize][rank], list_msgsArea.back());
[665]2260      }
[1460]2261      client->sendEvent(eventArea);
[665]2262    }
2263  }
2264
[667]2265  /*!
2266    Send longitude and latitude from client to servers
[1460]2267    Each client send long and lat information to corresponding connected clients(s).
[667]2268    Because longitude and latitude are optional, this function only called if latitude and longitude exist
2269  */
[665]2270  void CDomain::sendLonLat()
2271  {
2272    if (!hasLonLat) return;
2273
2274    int ns, n, i, j, ind, nv, idx;
[1460]2275    std::list<CContextClient*>::iterator it;
2276    for (it=clients.begin(); it!=clients.end(); ++it)
[300]2277    {
[1460]2278      CContextClient* client = *it;
2279      int serverSize = client->serverSize;
[584]2280
[1460]2281      // send lon lat for each connected server
2282      CEventClient eventLon(getType(), EVENT_ID_LON);
2283      CEventClient eventLat(getType(), EVENT_ID_LAT);
[509]2284
[1460]2285      list<CMessage> list_msgsLon, list_msgsLat;
2286      list<CArray<double,1> > list_lon, list_lat;
2287      list<CArray<double,2> > list_boundslon, list_boundslat;
[610]2288
[1460]2289      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2290      iteMap = indSrv_[serverSize].end();
2291      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
[467]2292      {
[1460]2293        int nbData = 0;
2294        int rank = connectedServerRank_[serverSize][k];
2295        it = indSrv_[serverSize].find(rank);
2296        if (iteMap != it)
2297          nbData = it->second.size();
[509]2298
[1460]2299        list_lon.push_back(CArray<double,1>(nbData));
2300        list_lat.push_back(CArray<double,1>(nbData));
2301
[467]2302        if (hasBounds)
[300]2303        {
[1460]2304          list_boundslon.push_back(CArray<double,2>(nvertex, nbData));
2305          list_boundslat.push_back(CArray<double,2>(nvertex, nbData));
2306        }
[610]2307
[1460]2308        CArray<double,1>& lon = list_lon.back();
2309        CArray<double,1>& lat = list_lat.back();
2310        const std::vector<size_t>& temp = it->second;
2311        for (n = 0; n < nbData; ++n)
2312        {
2313          idx = static_cast<int>(it->second[n]);
2314          int localInd = globalLocalIndexMap_[idx];
2315          lon(n) = lonvalue(localInd);
2316          lat(n) = latvalue(localInd);
2317
2318          if (hasBounds)
[449]2319          {
[1460]2320            CArray<double,2>& boundslon = list_boundslon.back();
2321            CArray<double,2>& boundslat = list_boundslat.back();
2322
2323            for (nv = 0; nv < nvertex; ++nv)
2324            {
2325              boundslon(nv, n) = bounds_lonvalue(nv, localInd);
2326              boundslat(nv, n) = bounds_latvalue(nv, localInd);
2327            }
[449]2328          }
[300]2329        }
[509]2330
[1460]2331        list_msgsLon.push_back(CMessage());
2332        list_msgsLat.push_back(CMessage());
[518]2333
[1460]2334        list_msgsLon.back() << this->getId() << hasLonLat;
2335        if (hasLonLat) 
2336          list_msgsLon.back() << list_lon.back();
2337        list_msgsLon.back()  << hasBounds;
2338        if (hasBounds)
2339        {
2340          list_msgsLon.back() << list_boundslon.back();
2341        }
[610]2342
[1460]2343        list_msgsLat.back() << this->getId() << hasLonLat;
2344        if (hasLonLat)
2345          list_msgsLat.back() << list_lat.back();
2346        list_msgsLat.back() << hasBounds;
2347        if (hasBounds)
2348        {         
2349          list_msgsLat.back() << list_boundslat.back();
2350        }
2351
2352        eventLon.push(rank, nbSenders[serverSize][rank], list_msgsLon.back());
2353        eventLat.push(rank, nbSenders[serverSize][rank], list_msgsLat.back());
[518]2354      }
[1460]2355      client->sendEvent(eventLon);
2356      client->sendEvent(eventLat);
[300]2357    }
2358  }
[509]2359
[667]2360  /*!
[1460]2361    Send data index to corresponding connected clients.
2362    Data index can be compressed however, we always send decompressed data index
2363    and they will be compressed on receiving.
2364    The compressed index are represented with 1 and others are represented with -1
[667]2365  */
[1460]2366  void CDomain::sendDataIndex()
[665]2367  {
[1460]2368    int ns, n, i, j, ind, nv, idx;
2369    std::list<CContextClient*>::iterator it;
2370    for (it=clients.begin(); it!=clients.end(); ++it)
2371    {
2372      CContextClient* client = *it;
2373
2374      int serverSize = client->serverSize;
2375
2376      // send area for each connected server
2377      CEventClient eventDataIndex(getType(), EVENT_ID_DATA_INDEX);
2378
2379      list<CMessage> list_msgsDataIndex;
2380      list<CArray<int,1> > list_data_i_index, list_data_j_index;
2381
2382      int nbIndex = i_index.numElements();
2383      int niByIndex = max(i_index) - min(i_index) + 1;
2384      int njByIndex = max(j_index) - min(j_index) + 1; 
2385      int dataIindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : niByIndex;
2386      int dataJindexBound = (1 == data_dim) ? (niByIndex * njByIndex) : njByIndex;
2387
2388     
2389      CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
2390      dataIIndex = -1; 
2391      dataJIndex = -1;
2392      ind = 0;
2393
2394      for (idx = 0; idx < data_i_index.numElements(); ++idx)
2395      {
2396        int dataIidx = data_i_index(idx) + data_ibegin;
2397        int dataJidx = data_j_index(idx) + data_jbegin;
2398        if ((0 <= dataIidx) && (dataIidx < dataIindexBound) &&
2399            (0 <= dataJidx) && (dataJidx < dataJindexBound))
2400        {
2401          dataIIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //i_index(dataIidx);//dataIidx;
2402          dataJIndex((1 == data_dim) ? dataIidx : dataJidx * ni + dataIidx) = 1; //j_index(dataJidx);//         
2403        }
2404      }
2405
2406      boost::unordered_map<int, vector<size_t> >::const_iterator it, iteMap;
2407      iteMap = indSrv_[serverSize].end();
2408      for (int k = 0; k < connectedServerRank_[serverSize].size(); ++k)
2409      {
2410        int nbData = 0;
2411        int rank = connectedServerRank_[serverSize][k];
2412        it = indSrv_[serverSize].find(rank);
2413        if (iteMap != it)
2414          nbData = it->second.size();
2415        list_data_i_index.push_back(CArray<int,1>(nbData));
2416        list_data_j_index.push_back(CArray<int,1>(nbData));
2417
2418        const std::vector<size_t>& temp = it->second;
2419        for (n = 0; n < nbData; ++n)
2420        {
2421          idx = static_cast<int>(it->second[n]);
2422          i = globalLocalIndexMap_[idx];
2423          list_data_i_index.back()(n) = dataIIndex(i);
2424          list_data_j_index.back()(n) = dataJIndex(i);
2425        }
2426
2427        list_msgsDataIndex.push_back(CMessage());
2428        list_msgsDataIndex.back() << this->getId();
2429        list_msgsDataIndex.back() << list_data_i_index.back() << list_data_j_index.back();
2430        eventDataIndex.push(rank, nbSenders[serverSize][rank], list_msgsDataIndex.back());
2431      }
2432      client->sendEvent(eventDataIndex);
2433    }
[665]2434  }
[1460]2435 
[300]2436  bool CDomain::dispatchEvent(CEventServer& event)
[610]2437  {
2438    if (SuperClass::dispatchEvent(event)) return true;
2439    else
2440    {
2441      switch(event.type)
[300]2442      {
[610]2443        case EVENT_ID_SERVER_ATTRIBUT:
[1460]2444          recvDistributionAttributes(event);
[610]2445          return true;
2446          break;
2447        case EVENT_ID_INDEX:
2448          recvIndex(event);
2449          return true;
2450          break;
[1460]2451        case EVENT_ID_MASK:
2452          recvMask(event);
2453          return true;
2454          break;
[610]2455        case EVENT_ID_LON:
2456          recvLon(event);
2457          return true;
2458          break;
2459        case EVENT_ID_LAT:
2460          recvLat(event);
2461          return true;
2462          break;
[611]2463        case EVENT_ID_AREA:
2464          recvArea(event);
2465          return true;
[1460]2466          break; 
2467        case EVENT_ID_DATA_INDEX:
2468          recvDataIndex(event);
2469          return true;
[611]2470          break;
[610]2471        default:
[762]2472          ERROR("bool CDomain::dispatchEvent(CEventServer& event)",
[610]2473                << "Unknown Event");
2474          return false;
2475       }
2476    }
2477  }
[509]2478
[667]2479  /*!
[1460]2480    Receive index event from clients(s)
2481    \param[in] event event contain info about rank and associated index
2482  */
2483  void CDomain::recvIndex(CEventServer& event)
2484  {
2485    string domainId;
2486    std::map<int, CBufferIn*> rankBuffers;
2487
2488    list<CEventServer::SSubEvent>::iterator it;
2489    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2490    {     
2491      CBufferIn* buffer = it->buffer;
2492      *buffer >> domainId;
2493      rankBuffers[it->rank] = buffer;       
2494    }
2495    get(domainId)->recvIndex(rankBuffers);
2496  }
2497
2498  /*!
2499    Receive index information from client(s). We use the global index for mapping index between
2500    sending clients and receiving clients.
2501    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
2502  */
2503  void CDomain::recvIndex(std::map<int, CBufferIn*>& rankBuffers)
2504  {
2505    int nbReceived = rankBuffers.size(), i, ind, index, type_int, iIndex, jIndex;
2506    recvClientRanks_.resize(nbReceived);       
2507
2508    std::map<int, CBufferIn*>::iterator it = rankBuffers.begin(), ite = rankBuffers.end();
2509    ind = 0;
2510    for (ind = 0; it != ite; ++it, ++ind)
2511    {       
2512       recvClientRanks_[ind] = it->first;
2513       CBufferIn& buffer = *(it->second);
2514       buffer >> type_int >> isCurvilinear >> indGlob_[it->first]; 
2515       type.setValue((type_attr::t_enum)type_int); // probleme des type enum avec les buffers : ToFix
2516    }
2517    int nbIndGlob = 0;
2518    for (i = 0; i < nbReceived; ++i)
2519    {
2520      nbIndGlob += indGlob_[recvClientRanks_[i]].numElements();
2521    }
2522   
2523    globalLocalIndexMap_.rehash(std::ceil(nbIndGlob/globalLocalIndexMap_.max_load_factor()));
2524    i_index.resize(nbIndGlob);
2525    j_index.resize(nbIndGlob);   
2526    int nbIndexGlobMax = nbIndGlob, nbIndLoc;
2527
2528    nbIndGlob = 0;
2529    for (i = 0; i < nbReceived; ++i)
2530    {
2531      CArray<int,1>& tmp = indGlob_[recvClientRanks_[i]];
2532      for (ind = 0; ind < tmp.numElements(); ++ind)
2533      {
2534         index = tmp(ind);
2535         if (0 == globalLocalIndexMap_.count(index))
2536         {
2537           iIndex = (index%ni_glo)-ibegin;
2538           iIndex = (iIndex < 0) ? 0 : iIndex;
2539           jIndex = (index/ni_glo)-jbegin;
2540           jIndex = (jIndex < 0) ? 0 : jIndex;
2541           nbIndLoc = iIndex + ni * jIndex;
2542           if (nbIndLoc < nbIndexGlobMax)
2543           {
2544             i_index(nbIndLoc) = index % ni_glo;
2545             j_index(nbIndLoc) = index / ni_glo;
2546             globalLocalIndexMap_[index] = nbIndLoc; 
2547             ++nbIndGlob;
2548           }
2549           // i_index(nbIndGlob) = index % ni_glo;
2550           // j_index(nbIndGlob) = index / ni_glo;
2551           // globalLocalIndexMap_[index] = nbIndGlob; 
2552           // ++nbIndGlob;
2553         } 
2554      } 
2555    } 
2556
2557    if (nbIndGlob==0)
2558    {
2559      i_index.resize(nbIndGlob);
2560      j_index.resize(nbIndGlob);
2561    }
2562    else
2563    {
2564      i_index.resizeAndPreserve(nbIndGlob);
2565      j_index.resizeAndPreserve(nbIndGlob);
2566    }
2567  }
2568
2569  /*!
[667]2570    Receive attributes event from clients(s)
2571    \param[in] event event contain info about rank and associated attributes
2572  */
[1460]2573  void CDomain::recvDistributionAttributes(CEventServer& event)
[300]2574  {
2575    CBufferIn* buffer=event.subEvents.begin()->buffer;
2576    string domainId ;
2577    *buffer>>domainId ;
[1460]2578    get(domainId)->recvDistributionAttributes(*buffer);
[300]2579  }
[509]2580
[667]2581  /*!
2582    Receive attributes from client(s): zoom info and begin and n of each server
2583    \param[in] rank rank of client source
2584    \param[in] buffer message containing attributes info
2585  */
[1460]2586  void CDomain::recvDistributionAttributes(CBufferIn& buffer)
[300]2587  {
[1460]2588    int ni_tmp, ibegin_tmp, nj_tmp, jbegin_tmp;
[821]2589    int global_zoom_ni_tmp, global_zoom_ibegin_tmp, global_zoom_nj_tmp, global_zoom_jbegin_tmp;
[1460]2590    buffer >> isUnstructed_ >> ni_tmp >> ibegin_tmp >> nj_tmp >> jbegin_tmp
2591           >> global_zoom_ni_tmp >> global_zoom_ibegin_tmp >> global_zoom_nj_tmp >> global_zoom_jbegin_tmp           
[676]2592           >> isCompressible_;
[1460]2593    ni.setValue(ni_tmp);
2594    ibegin.setValue(ibegin_tmp);
2595    nj.setValue(nj_tmp);
2596    jbegin.setValue(jbegin_tmp);
[300]2597
[821]2598    global_zoom_ni.setValue(global_zoom_ni_tmp);
2599    global_zoom_ibegin.setValue(global_zoom_ibegin_tmp);
2600    global_zoom_nj.setValue(global_zoom_nj_tmp);
2601    global_zoom_jbegin.setValue(global_zoom_jbegin_tmp);
2602
[1460]2603    int iend = ibegin + ni  - 1;
2604    int jend = jbegin + nj  - 1;
2605    int zoom_iend_glob = global_zoom_ibegin + global_zoom_ni - 1;
2606    int zoom_jend_glob = global_zoom_jbegin + global_zoom_nj - 1;
[509]2607
[1460]2608    zoom_ibegin.setValue(global_zoom_ibegin > ibegin ? global_zoom_ibegin : ibegin);
2609    int zoom_iend = zoom_iend_glob < iend ? zoom_iend_glob : iend ;
2610    zoom_ni.setValue(zoom_iend-zoom_ibegin+1);
[509]2611
[1460]2612    zoom_jbegin.setValue(global_zoom_jbegin > jbegin ? global_zoom_jbegin : jbegin);
2613    int zoom_jend = zoom_jend_glob < jend ? zoom_jend_glob : jend ;
2614    zoom_nj.setValue(zoom_jend-zoom_jbegin+1);
[300]2615
[1460]2616    if (zoom_ni<=0 || zoom_nj<=0)
[300]2617    {
[1460]2618      zoom_ni=0 ; zoom_ibegin=global_zoom_ibegin ; //=0; zoom_iend=0 ;
2619      zoom_nj=0 ; zoom_jbegin=global_zoom_jbegin ; //=0; zoom_jend=0 ;
[300]2620    }
[611]2621
[300]2622  }
[509]2623
[667]2624  /*!
[1460]2625    Receive area event from clients(s)
2626    \param[in] event event contain info about rank and associated area
[667]2627  */
[1460]2628  void CDomain::recvMask(CEventServer& event)
[610]2629  {
[1460]2630    string domainId;
2631    std::map<int, CBufferIn*> rankBuffers;
[676]2632
[610]2633    list<CEventServer::SSubEvent>::iterator it;
2634    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
[1460]2635    {     
[610]2636      CBufferIn* buffer = it->buffer;
2637      *buffer >> domainId;
[1460]2638      rankBuffers[it->rank] = buffer;     
[610]2639    }
[1460]2640    get(domainId)->recvMask(rankBuffers);
2641  }
[676]2642
2643
[667]2644  /*!
[1460]2645    Receive mask information from client(s)
2646    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
[667]2647  */
[1460]2648  void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)
[610]2649  {
[1460]2650    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2651    if (nbReceived != recvClientRanks_.size())
2652      ERROR("void CDomain::recvMask(std::map<int, CBufferIn*>& rankBuffers)",
2653           << "The number of sending clients is not correct."
2654           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
[676]2655
[1460]2656    vector<CArray<bool,1> > recvMaskValue(nbReceived);     
2657    for (i = 0; i < recvClientRanks_.size(); ++i)
[676]2658    {
[1460]2659      int rank = recvClientRanks_[i];
2660      CBufferIn& buffer = *(rankBuffers[rank]);     
2661      buffer >> recvMaskValue[i];
[676]2662    }
[1460]2663
2664    int nbMaskInd = 0;
2665    for (i = 0; i < nbReceived; ++i)
2666    {
2667      nbMaskInd += recvMaskValue[i].numElements();
2668    }
2669 
2670    if (nbMaskInd != globalLocalIndexMap_.size())
2671      info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2672               << "Something must be wrong with mask index "<< std::endl;
2673
2674    nbMaskInd = globalLocalIndexMap_.size();
2675    mask_1d.resize(nbMaskInd);
2676    domainMask.resize(nbMaskInd);
2677    mask_1d = false;
2678   
2679    for (i = 0; i < nbReceived; ++i)
2680    {
2681      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2682      CArray<bool,1>& tmp = recvMaskValue[i];
2683      for (ind = 0; ind < tmp.numElements(); ++ind)
2684      {
2685        lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2686        if (!mask_1d(lInd)) // Only rewrite mask_1d if it's not true
2687          mask_1d(lInd) = tmp(ind);
2688      }
2689    }
2690    domainMask=mask_1d ;
[610]2691  }
2692
[667]2693  /*!
2694    Receive longitude event from clients(s)
2695    \param[in] event event contain info about rank and associated longitude
2696  */
[518]2697  void CDomain::recvLon(CEventServer& event)
[300]2698  {
[1460]2699    string domainId;
2700    std::map<int, CBufferIn*> rankBuffers;
2701
[610]2702    list<CEventServer::SSubEvent>::iterator it;
2703    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
[1460]2704    {     
[610]2705      CBufferIn* buffer = it->buffer;
2706      *buffer >> domainId;
[1460]2707      rankBuffers[it->rank] = buffer;       
[300]2708    }
[1460]2709    get(domainId)->recvLon(rankBuffers);
[300]2710  }
[509]2711
[667]2712  /*!
2713    Receive longitude information from client(s)
[1460]2714    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
[667]2715  */
[1460]2716  void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)
[300]2717  {
[1460]2718    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2719    if (nbReceived != recvClientRanks_.size())
2720      ERROR("void CDomain::recvLon(std::map<int, CBufferIn*>& rankBuffers)",
2721           << "The number of sending clients is not correct."
2722           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
[518]2723
[1460]2724    vector<CArray<double,1> > recvLonValue(nbReceived);
2725    vector<CArray<double,2> > recvBoundsLonValue(nbReceived);   
2726    for (i = 0; i < recvClientRanks_.size(); ++i)
2727    {
2728      int rank = recvClientRanks_[i];
2729      CBufferIn& buffer = *(rankBuffers[rank]);
2730      buffer >> hasLonLat;
2731      if (hasLonLat)
2732        buffer >> recvLonValue[i];
2733      buffer >> hasBounds;
2734      if (hasBounds)
2735        buffer >> recvBoundsLonValue[i];
2736    }
[924]2737
[1460]2738    if (hasLonLat)
2739    {
2740      int nbLonInd = 0;
2741      for (i = 0; i < nbReceived; ++i)
2742      {
2743        nbLonInd += recvLonValue[i].numElements();
2744      }
2745   
2746      if (nbLonInd != globalLocalIndexMap_.size())
2747        info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2748                 << "Something must be wrong with longitude index "<< std::endl;
[518]2749
[1460]2750      nbLonInd = globalLocalIndexMap_.size();
2751      lonvalue.resize(nbLonInd);
[518]2752      if (hasBounds)
2753      {
[1460]2754        bounds_lonvalue.resize(nvertex,nbLonInd);
2755        bounds_lonvalue = 0.;
[518]2756      }
[1460]2757
2758      nbLonInd = 0;
2759      for (i = 0; i < nbReceived; ++i)
2760      {
2761        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2762        CArray<double,1>& tmp = recvLonValue[i];
2763        for (ind = 0; ind < tmp.numElements(); ++ind)
2764        {
2765          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2766          lonvalue(lInd) = tmp(ind); 
2767           if (hasBounds)
2768           {         
2769            for (int nv = 0; nv < nvertex; ++nv)
2770              bounds_lonvalue(nv, lInd) = recvBoundsLonValue[i](nv, ind);
2771           }                 
2772        }
2773      }       
[518]2774    }
2775  }
2776
[667]2777  /*!
2778    Receive latitude event from clients(s)
2779    \param[in] event event contain info about rank and associated latitude
2780  */
[518]2781  void CDomain::recvLat(CEventServer& event)
2782  {
[1460]2783    string domainId;
2784    std::map<int, CBufferIn*> rankBuffers;
2785
[610]2786    list<CEventServer::SSubEvent>::iterator it;
2787    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
[1460]2788    {     
[610]2789      CBufferIn* buffer = it->buffer;
2790      *buffer >> domainId;
[1460]2791      rankBuffers[it->rank] = buffer;   
[518]2792    }
[1460]2793    get(domainId)->recvLat(rankBuffers);
[518]2794  }
2795
[667]2796  /*!
2797    Receive latitude information from client(s)
[1460]2798    \param[in] rankBuffers rank of sending client and the corresponding receive buffer 
[667]2799  */
[1460]2800  void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)
[518]2801  {
[1460]2802    int nbReceived = rankBuffers.size(), i, ind, index, iindex, jindex, lInd;
2803    if (nbReceived != recvClientRanks_.size())
2804      ERROR("void CDomain::recvLat(std::map<int, CBufferIn*>& rankBuffers)",
2805           << "The number of sending clients is not correct."
2806           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
[509]2807
[1460]2808    vector<CArray<double,1> > recvLatValue(nbReceived);
2809    vector<CArray<double,2> > recvBoundsLatValue(nbReceived);   
2810    for (i = 0; i < recvClientRanks_.size(); ++i)
2811    {
2812      int rank = recvClientRanks_[i];
2813      CBufferIn& buffer = *(rankBuffers[rank]);
2814      buffer >> hasLonLat;
2815      if (hasLonLat)
2816        buffer >> recvLatValue[i];
2817      buffer >> hasBounds;
2818      if (hasBounds)
2819        buffer >> recvBoundsLatValue[i];
2820    }
[610]2821
[1460]2822    if (hasLonLat)
[300]2823    {
[1460]2824      int nbLatInd = 0;
2825      for (i = 0; i < nbReceived; ++i)
2826      {
2827        nbLatInd += recvLatValue[i].numElements();
2828      }
2829   
2830      if (nbLatInd != globalLocalIndexMap_.size())
2831        info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2832                << "Something must be wrong with latitude index "<< std::endl;
2833
2834      nbLatInd = globalLocalIndexMap_.size();
2835      latvalue.resize(nbLatInd);
[509]2836      if (hasBounds)
[449]2837      {
[1460]2838        bounds_latvalue.resize(nvertex,nbLatInd);
2839        bounds_latvalue = 0. ;
[449]2840      }
[1460]2841
2842      nbLatInd = 0;
2843      for (i = 0; i < nbReceived; ++i)
2844      {
2845        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2846        CArray<double,1>& tmp = recvLatValue[i];
2847        for (ind = 0; ind < tmp.numElements(); ++ind)
2848        {
2849          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2850          latvalue(lInd) = tmp(ind);   
2851           if (hasBounds)
2852           {
2853            CArray<double,2>& boundslat = recvBoundsLatValue[i];
2854            for (int nv = 0; nv < nvertex; ++nv)
2855              bounds_latvalue(nv, lInd) = boundslat(nv, ind);
2856           }   
2857          ++nbLatInd;
2858        }
2859      }       
[300]2860    }
2861  }
[553]2862
[667]2863  /*!
2864    Receive area event from clients(s)
2865    \param[in] event event contain info about rank and associated area
2866  */
[611]2867  void CDomain::recvArea(CEventServer& event)
2868  {
[1460]2869    string domainId;
2870    std::map<int, CBufferIn*> rankBuffers;
2871
[611]2872    list<CEventServer::SSubEvent>::iterator it;
2873    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
[1460]2874    {     
[611]2875      CBufferIn* buffer = it->buffer;
2876      *buffer >> domainId;
[1460]2877      rankBuffers[it->rank] = buffer;     
[611]2878    }
[1460]2879    get(domainId)->recvArea(rankBuffers);
[611]2880  }
2881
[667]2882  /*!
2883    Receive area information from client(s)
[1460]2884    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
[667]2885  */
[1460]2886  void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)
[611]2887  {
[1460]2888    int nbReceived = rankBuffers.size(), i, ind, index, lInd;
2889    if (nbReceived != recvClientRanks_.size())
2890      ERROR("void CDomain::recvArea(std::map<int, CBufferIn*>& rankBuffers)",
2891           << "The number of sending clients is not correct."
2892           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
[611]2893
[1460]2894    vector<CArray<double,1> > recvAreaValue(nbReceived);     
2895    for (i = 0; i < recvClientRanks_.size(); ++i)
2896    {
2897      int rank = recvClientRanks_[i];
2898      CBufferIn& buffer = *(rankBuffers[rank]);     
2899      buffer >> hasArea;
2900      if (hasArea)
2901        buffer >> recvAreaValue[i];
2902    }
[611]2903
[1460]2904    if (hasArea)
[611]2905    {
[1460]2906      int nbAreaInd = 0;
2907      for (i = 0; i < nbReceived; ++i)
2908      {     
2909        nbAreaInd += recvAreaValue[i].numElements();
2910      }
2911
2912      if (nbAreaInd != globalLocalIndexMap_.size())
2913        info (0) << "If the domain " << this->getDomainOutputName() <<" does not have overlapped region between processes."
2914                 << "Something must be wrong with area index "<< std::endl;
2915
2916      nbAreaInd = globalLocalIndexMap_.size();
2917      areavalue.resize(nbAreaInd);
2918      nbAreaInd = 0;     
2919      for (i = 0; i < nbReceived; ++i)
2920      {
2921        CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
2922        CArray<double,1>& tmp = recvAreaValue[i];
2923        for (ind = 0; ind < tmp.numElements(); ++ind)
2924        {
2925          lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
2926          areavalue(lInd) = tmp(ind);         
2927        }
2928      }
2929     
[611]2930    }
2931  }
2932
[1106]2933  /*!
2934    Compare two domain objects.
2935    They are equal if only if they have identical attributes as well as their values.
2936    Moreover, they must have the same transformations.
2937  \param [in] domain Compared domain
2938  \return result of the comparison
2939  */
2940  bool CDomain::isEqual(CDomain* obj)
2941  {
[1117]2942    vector<StdString> excludedAttr;
2943    excludedAttr.push_back("domain_ref");
2944    bool objEqual = SuperClass::isEqual(obj, excludedAttr);
[1106]2945    if (!objEqual) return objEqual;
2946
2947    TransMapTypes thisTrans = this->getAllTransformations();
2948    TransMapTypes objTrans  = obj->getAllTransformations();
2949
2950    TransMapTypes::const_iterator it, itb, ite;
2951    std::vector<ETranformationType> thisTransType, objTransType;
2952    for (it = thisTrans.begin(); it != thisTrans.end(); ++it)
2953      thisTransType.push_back(it->first);
2954    for (it = objTrans.begin(); it != objTrans.end(); ++it)
2955      objTransType.push_back(it->first);
2956
2957    if (thisTransType.size() != objTransType.size()) return false;
2958    for (int idx = 0; idx < thisTransType.size(); ++idx)
2959      objEqual &= (thisTransType[idx] == objTransType[idx]);
2960
2961    return objEqual;
2962  }
2963
[1460]2964  /*!
2965    Receive data index event from clients(s)
2966    \param[in] event event contain info about rank and associated index
2967  */
2968  void CDomain::recvDataIndex(CEventServer& event)
2969  {
2970    string domainId;
2971    std::map<int, CBufferIn*> rankBuffers;
2972
2973    list<CEventServer::SSubEvent>::iterator it;
2974    for (it = event.subEvents.begin(); it != event.subEvents.end(); ++it)
2975    {     
2976      CBufferIn* buffer = it->buffer;
2977      *buffer >> domainId;
2978      rankBuffers[it->rank] = buffer;       
2979    }
2980    get(domainId)->recvDataIndex(rankBuffers);
2981  }
2982
2983  /*!
2984    Receive data index information from client(s)
2985    A client receives data index from different clients to rebuild its own data index.
2986    Because we use global index + mask info to calculate the sending data to client(s),
2987    this data index must be updated with mask info (maybe it will change in the future)
2988    Because the data index is local, to rebuild data index of received client, we should use global index along with.
2989
2990    \param[in] rankBuffers rank of sending client and the corresponding receive buffer     
2991  */
2992  void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)
2993  {
2994    int nbReceived = rankBuffers.size(), i, ind, index, indexI, indexJ, type_int, lInd;   
2995    if (nbReceived != recvClientRanks_.size())
2996      ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
2997           << "The number of sending clients is not correct."
2998           << "Expected number: " << recvClientRanks_.size() << " but received " << nbReceived);
2999
3000    vector<CArray<int,1> > recvDataIIndex(nbReceived),recvDataJIndex(nbReceived);     
3001    for (i = 0; i < recvClientRanks_.size(); ++i)
3002    {
3003      int rank = recvClientRanks_[i];
3004      CBufferIn& buffer = *(rankBuffers[rank]);
3005      buffer >> recvDataIIndex[i];
3006      buffer >> recvDataJIndex[i];
3007    }
3008   
3009    int nbIndex = i_index.numElements();
3010    CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
3011    dataIIndex = -1; dataJIndex = -1;
3012     
3013    nbIndex = 0;
3014    for (i = 0; i < nbReceived; ++i)
3015    {     
3016      CArray<int,1>& tmpInd = indGlob_[recvClientRanks_[i]];
3017      CArray<int,1>& tmpI = recvDataIIndex[i];   
3018      CArray<int,1>& tmpJ = recvDataJIndex[i];     
3019      if ((tmpI.numElements() != tmpInd.numElements()) || (tmpJ.numElements() != tmpInd.numElements()))
3020          ERROR("void CDomain::recvDataIndex(std::map<int, CBufferIn*>& rankBuffers)",
3021             << "The number of global received index is not coherent with the number of received data index."
3022             << "Expected number of global index: " << tmpI.numElements() << " but received " << tmpInd.numElements());
3023
3024      for (ind = 0; ind < tmpI.numElements(); ++ind)
3025      {
3026         lInd = globalLocalIndexMap_[size_t(tmpInd(ind))];
3027         dataIIndex(lInd) = (-1 == dataIIndex(lInd)) ? tmpI(ind) : dataIIndex(lInd); // Only fill in dataIndex if there is no data
3028         dataJIndex(lInd) = (-1 == dataJIndex(lInd)) ? tmpJ(ind) : dataJIndex(lInd); 
3029
3030         if (!domainMask(lInd))   // Include mask info into data index on the RECEIVE getServerDimensionSizes   
3031         {
3032           dataIIndex(lInd) = dataJIndex(lInd) = -1;
3033         }
3034      } 
3035    }
3036
3037    int nbCompressedData = 0; 
3038    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3039    {
3040       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3041       if ((0 <= indexI) && (0 <= indexJ))
3042         ++nbCompressedData;
3043    }       
3044 
3045    data_i_index.resize(nbCompressedData);
3046    data_j_index.resize(nbCompressedData);
3047
3048    nbCompressedData = 0; 
3049    for (ind = 0; ind < dataIIndex.numElements(); ++ind)
3050    {
3051       indexI = dataIIndex(ind); indexJ = dataJIndex(ind);
3052       if ((0 <= indexI) && (0 <= indexJ))
3053       {
3054          data_i_index(nbCompressedData) = (1 == data_dim) ? ind : ind % ni;
3055          data_j_index(nbCompressedData) = (1 == data_dim) ? 0   : ind / ni;
3056         ++nbCompressedData;
3057       }
3058    }
3059
3060    // Reset data_ibegin, data_jbegin
3061    data_ibegin.setValue(0);
3062    data_jbegin.setValue(0);
3063  }
3064
[836]3065  CTransformation<CDomain>* CDomain::addTransformation(ETranformationType transType, const StdString& id)
3066  {
3067    transformationMap_.push_back(std::make_pair(transType, CTransformation<CDomain>::createTransformation(transType,id)));
3068    return transformationMap_.back().second;
3069  }
3070
[667]3071  /*!
3072    Check whether a domain has transformation
3073    \return true if domain has transformation
3074  */
[631]3075  bool CDomain::hasTransformation()
3076  {
3077    return (!transformationMap_.empty());
3078  }
3079
[667]3080  /*!
3081    Set transformation for current domain. It's the method to move transformation in hierarchy
3082    \param [in] domTrans transformation on domain
3083  */
[631]3084  void CDomain::setTransformations(const TransMapTypes& domTrans)
3085  {
3086    transformationMap_ = domTrans;
3087  }
3088
[667]3089  /*!
3090    Get all transformation current domain has
3091    \return all transformation
3092  */
[631]3093  CDomain::TransMapTypes CDomain::getAllTransformations(void)
3094  {
3095    return transformationMap_;
3096  }
3097
[823]3098  void CDomain::duplicateTransformation(CDomain* src)
3099  {
3100    if (src->hasTransformation())
3101    {
3102      this->setTransformations(src->getAllTransformations());
3103    }
3104  }
3105
[667]3106  /*!
[747]3107   * Go through the hierarchy to find the domain from which the transformations must be inherited
3108   */
[631]3109  void CDomain::solveInheritanceTransformation()
3110  {
[747]3111    if (hasTransformation() || !hasDirectDomainReference())
3112      return;
[631]3113
[747]3114    CDomain* domain = this;
3115    std::vector<CDomain*> refDomains;
3116    while (!domain->hasTransformation() && domain->hasDirectDomainReference())
[631]3117    {
[747]3118      refDomains.push_back(domain);
3119      domain = domain->getDirectDomainReference();
[631]3120    }
3121
[747]3122    if (domain->hasTransformation())
3123      for (size_t i = 0; i < refDomains.size(); ++i)
3124        refDomains[i]->setTransformations(domain->getAllTransformations());
[631]3125  }
3126
[1460]3127  void CDomain::setContextClient(CContextClient* contextClient)
3128  {
3129    if (clientsSet.find(contextClient)==clientsSet.end())
3130    {
3131      clients.push_back(contextClient) ;
3132      clientsSet.insert(contextClient);
3133    }
3134  }
3135
[667]3136  /*!
3137    Parse children nodes of a domain in xml file.
3138    Whenver there is a new transformation, its type and name should be added into this function
3139    \param node child node to process
3140  */
[631]3141  void CDomain::parse(xml::CXMLNode & node)
3142  {
3143    SuperClass::parse(node);
3144
3145    if (node.goToChildElement())
3146    {
[836]3147      StdString nodeElementName;
[631]3148      do
3149      {
[784]3150        StdString nodeId("");
3151        if (node.getAttributes().end() != node.getAttributes().find("id"))
3152        { nodeId = node.getAttributes()["id"]; }
3153
[836]3154        nodeElementName = node.getElementName();
[1134]3155        if(transformationMapList_ptr == 0) initializeTransformationMap();
[1334]3156
[1328]3157        std::map<StdString, ETranformationType>::const_iterator ite = transformationMapList_ptr->end(), it;
[1334]3158
[1328]3159        it = transformationMapList_ptr->find(nodeElementName);
[836]3160        if (ite != it)
[657]3161        {
[836]3162          transformationMap_.push_back(std::make_pair(it->second, CTransformation<CDomain>::createTransformation(it->second,
3163                                                                                                                nodeId,
3164                                                                                                                &node)));
[657]3165        }
[968]3166        else
3167        {
3168          ERROR("void CDomain::parse(xml::CXMLNode & node)",
3169                << "The transformation " << nodeElementName << " has not been supported yet.");
3170        }
[631]3171      } while (node.goToNextElement()) ;
3172      node.goToParentElement();
3173    }
3174  }
[219]3175   //----------------------------------------------------------------
[509]3176
[540]3177   DEFINE_REF_FUNC(Domain,domain)
[509]3178
[219]3179   ///---------------------------------------------------------------
3180
[335]3181} // namespace xios
Note: See TracBrowser for help on using the repository browser.