Changeset 1037 for XIOS/dev/branch_yushan/src/node

Timestamp:

01/25/17 16:25:17 (7 years ago)

Author:

yushan

Message:

initialize the branch

Location:

XIOS/dev/branch_yushan/src/node

Files:

• axis.cpp (modified) (1 diff)
• context.cpp (modified) (9 diffs)
• context.hpp (modified) (1 diff)
• domain.cpp (modified) (7 diffs)
• field.cpp (modified) (19 diffs)
• field.hpp (modified) (5 diffs)
• field_impl.hpp (modified) (1 diff)
• file.cpp (modified) (10 diffs)
• file.hpp (modified) (1 diff)
• grid.cpp (modified) (14 diffs)
• grid.hpp (modified) (1 diff)
• interpolate_domain.cpp (modified) (1 diff)
• interpolate_domain.hpp (modified) (1 diff)
• mesh.cpp (modified) (32 diffs)

XIOS/dev/branch_yushan/src/node/axis.cpp

@@ -264 +264 @@
       isDistributed_ = (!this->begin.isEmpty() && !this->n.isEmpty() && (this->begin + this->n < this->n_glo)) ||
                        (!this->n.isEmpty() && (this->n != this->n_glo));
-
-      // A same stupid condition to make sure that if there is only one client, axis
-      // should be considered to be distributed. This should be a temporary solution     
-      isDistributed_ |= (1 == CContext::getCurrent()->client->clientSize);
    }
 

XIOS/dev/branch_yushan/src/node/context.cpp

@@ -339 +339 @@
    }
 
+   //! Server side: Put server into a loop in order to listen message from client
+   bool CContext::eventLoop(void)
+   {
+     return server->eventLoop();
+   }
+
    //! Try to send the buffers and receive possible answers
    bool CContext::checkBuffersAndListen(void)
    {
      client->checkBuffers();
-
-     bool hasTmpBufferedEvent = client->hasTemporarilyBufferedEvent();
-     if (hasTmpBufferedEvent)
-       hasTmpBufferedEvent = !client->sendTemporarilyBufferedEvent();
-
-     // Don't process events if there is a temporarily buffered event
-     return server->eventLoop(!hasTmpBufferedEvent);
+     return server->eventLoop();
    }
 
@@ -389 +389 @@
    void CContext::closeDefinition(void)
    {
+     int myRank;
+     MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+     //printf("myRank = %d, hasClient = %d, hasServer = %d\n", myRank, hasClient, hasServer);
+
      // There is nothing client need to send to server
      if (hasClient)
      {
        // After xml is parsed, there are some more works with post processing
-       postProcessing();
-     }
-     setClientServerBuffer();
+       postProcessing(); 
+       //printf("myRank = %d,                postProcessing OK\n", myRank);
+     }
+     setClientServerBuffer(); //printf("myRank = %d, setClientServerBuffer OK\n", myRank);
 
      if (hasClient && !hasServer)
      {
       // Send all attributes of current context to server
-      this->sendAllAttributesToServer();
+      this->sendAllAttributesToServer(); //printf("myRank = %d, this->sendAllAttributesToServer OK\n", myRank);
 
       // Send all attributes of current calendar
       CCalendarWrapper::get(CCalendarWrapper::GetDefName())->sendAllAttributesToServer();
+      //printf("myRank = %d, CCalendarWrapper::get(CCalendarWrapper::GetDefName())->sendAllAttributesToServer OK\n", myRank);
 
       // We have enough information to send to server
       // First of all, send all enabled files
-       sendEnabledFiles();
+       sendEnabledFiles();  //printf("myRank = %d, sendEnabledFiles OK\n", myRank);
 
       // Then, send all enabled fields
-       sendEnabledFields();
+       sendEnabledFields();  //printf("myRank = %d, sendEnabledFields OK\n", myRank);
 
       // At last, we have all info of domain and axis, then send them
-       sendRefDomainsAxis();
+       sendRefDomainsAxis();  //printf("myRank = %d, sendRefDomainsAxis OK\n", myRank);
 
       // After that, send all grid (if any)
-       sendRefGrid();
+       sendRefGrid(); //printf("myRank = %d, sendRefGrid OK\n", myRank);
     }
 
     // We have a xml tree on the server side and now, it should be also processed
-    if (hasClient && !hasServer) sendPostProcessing();
+    if (hasClient && !hasServer) sendPostProcessing();  
 
     // There are some processings that should be done after all of above. For example: check mask or index
     if (hasClient)
     {
-      this->buildFilterGraphOfEnabledFields();
-      buildFilterGraphOfFieldsWithReadAccess();
-      this->solveAllRefOfEnabledFields(true);
+      this->buildFilterGraphOfEnabledFields();  //printf("myRank = %d, buildFilterGraphOfEnabledFields OK\n", myRank);
+      buildFilterGraphOfFieldsWithReadAccess();  //printf("myRank = %d, buildFilterGraphOfFieldsWithReadAccess OK\n", myRank);
+      this->solveAllRefOfEnabledFields(true);  //printf("myRank = %d, solveAllRefOfEnabledFields OK\n", myRank);
     }
 
@@ -434 +441 @@
 
     // Nettoyage de l'arborescence
-    if (hasClient && !hasServer) CleanTree(); // Only on client side??
+    if (hasClient && !hasServer) CleanTree();
 
     if (hasClient)
     {
-      sendCreateFileHeader();
-
-      startPrefetchingOfEnabledReadModeFiles();
+      sendCreateFileHeader();  //printf("myRank = %d, sendCreateFileHeader OK\n", myRank);
+
+      startPrefetchingOfEnabledReadModeFiles();  //printf("myRank = %d, startPrefetchingOfEnabledReadModeFiles OK\n", myRank);
     }
    }
@@ -479 +486 @@
    {
      int size = this->enabledFiles.size();
+     
      for (int i = 0; i < size; ++i)
      {
@@ -566 +574 @@
    {
       const std::vector<CFile*> allFiles = CFile::getAll();
-      const CDate& initDate = calendar->getInitDate();
 
       for (unsigned int i = 0; i < allFiles.size(); i++)
@@ -572 +579 @@
          {
             if (allFiles[i]->enabled.getValue()) // Si l'attribut 'enabled' est fixé à vrai.
-            {
-              if ((initDate + allFiles[i]->output_freq.getValue()) < (initDate + this->getCalendar()->getTimeStep()))
-              {
-                error(0)<<"WARNING: void CContext::findEnabledFiles()"<<endl
-                    << "Output frequency in file \""<<allFiles[i]->getFileOutputName()
-                    <<"\" is less than the time step. File will not be written."<<endl;
-              }
-              else
                enabledFiles.push_back(allFiles[i]);
-            }
          }
-         else
-         {
-           if ( (initDate + allFiles[i]->output_freq.getValue()) < (initDate + this->getCalendar()->getTimeStep()))
-           {
-             error(0)<<"WARNING: void CContext::findEnabledFiles()"<<endl
-                 << "Output frequency in file \""<<allFiles[i]->getFileOutputName()
-                 <<"\" is less than the time step. File will not be written."<<endl;
-           }
-           else
-             enabledFiles.push_back(allFiles[i]); // otherwise true by default
-         }
+         else enabledFiles.push_back(allFiles[i]); // otherwise true by default
+
 
       if (enabledFiles.size() == 0)
@@ -812 +801 @@
    void CContext::postProcessing()
    {
+     int myRank;
+     MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+     //printf("myRank = %d, in postProcessing, isPostProcessed = %d\n", myRank, isPostProcessed);
      if (isPostProcessed) return;
 
@@ -820 +813 @@
         ERROR("CContext::postProcessing()", << "A timestep must be defined for the context \"" << getId() << "!\"")
       // Calendar first update to set the current date equals to the start date
-      calendar->update(0);
+      calendar->update(0);  //printf("myRank = %d, calendar->update(0) OK\n", myRank);
 
       // Find all inheritance in xml structure
-      this->solveAllInheritance();
+      this->solveAllInheritance();  //printf("myRank = %d, this->solveAllInheritance OK\n", myRank);
 
       // Check if some axis, domains or grids are eligible to for compressed indexed output.
       // Warning: This must be done after solving the inheritance and before the rest of post-processing
-      checkAxisDomainsGridsEligibilityForCompressedOutput();
+      checkAxisDomainsGridsEligibilityForCompressedOutput();  //printf("myRank = %d, checkAxisDomainsGridsEligibilityForCompressedOutput OK\n", myRank);
 
       // Check if some automatic time series should be generated
       // Warning: This must be done after solving the inheritance and before the rest of post-processing
-      prepareTimeseries();
+      prepareTimeseries();  //printf("myRank = %d, prepareTimeseries OK\n", myRank);
 
       //Initialisation du vecteur 'enabledFiles' contenant la liste des fichiers à sortir.
-      this->findEnabledFiles();
-      this->findEnabledReadModeFiles();
+      this->findEnabledFiles();  //printf("myRank = %d, this->findEnabledFiles OK\n", myRank);
+      this->findEnabledReadModeFiles();  //printf("myRank = %d, this->findEnabledReadModeFiles OK\n", myRank);
 
       // Find all enabled fields of each file
-      this->findAllEnabledFields();
-      this->findAllEnabledFieldsInReadModeFiles();
+      this->findAllEnabledFields();  //printf("myRank = %d, this->findAllEnabledFields OK\n", myRank);
+      this->findAllEnabledFieldsInReadModeFiles();  //printf("myRank = %d, this->findAllEnabledFieldsInReadModeFiles OK\n", myRank);
 
      if (hasClient && !hasServer)
      {
       // Try to read attributes of fields in file then fill in corresponding grid (or domain, axis)
-      this->readAttributesOfEnabledFieldsInReadModeFiles();
+      this->readAttributesOfEnabledFieldsInReadModeFiles();  //printf("myRank = %d, this->readAttributesOfEnabledFieldsInReadModeFiles OK\n", myRank);
      }
 
       // Only search and rebuild all reference objects of enable fields, don't transform
-      this->solveOnlyRefOfEnabledFields(false);
+      this->solveOnlyRefOfEnabledFields(false);  //printf("myRank = %d, this->solveOnlyRefOfEnabledFields(false) OK\n", myRank);
 
       // Search and rebuild all reference object of enabled fields
-      this->solveAllRefOfEnabledFields(false);
+      this->solveAllRefOfEnabledFields(false);  //printf("myRank = %d, this->solveAllRefOfEnabledFields(false) OK\n", myRank);
 
       // Find all fields with read access from the public API
-      findFieldsWithReadAccess();
+      findFieldsWithReadAccess();  //printf("myRank = %d, findFieldsWithReadAccess OK\n", myRank);
       // and solve the all reference for them
-      solveAllRefOfFieldsWithReadAccess();
+      solveAllRefOfFieldsWithReadAccess();  //printf("myRank = %d, solveAllRefOfFieldsWithReadAccess OK\n", myRank);
 
       isPostProcessed = true;
@@ -1240 +1233 @@
          event.push(*itRank,1,msg);
        client->sendEvent(event);
-     }
-     else client->sendEvent(event);
+    }
+    else client->sendEvent(event);
   }
 

XIOS/dev/branch_yushan/src/node/context.hpp

@@ -93 +93 @@
 
          // Put sever or client into loop state
+         bool eventLoop(void);
+
          bool checkBuffersAndListen(void);
 

XIOS/dev/branch_yushan/src/node/domain.cpp

@@ -698 +698 @@
      computeNGlobDomain();
      checkZoom();
-     
+
      isDistributed_ = !((!ni.isEmpty() && (ni == ni_glo) && !nj.isEmpty() && (nj == nj_glo)) ||
                         (!i_index.isEmpty() && i_index.numElements() == ni_glo*nj_glo));
-
-     // A stupid condition to make sure that if there is only one client, domain
-     // should be considered to be distributed. This should be a temporary solution     
-     isDistributed_ |= (1 == CContext::getCurrent()->client->clientSize);
    }
 
@@ -1307 +1303 @@
    void CDomain::checkAttributesOnClientAfterTransformation()
    {
-     CContext* context=CContext::getCurrent() ;
-
+     int myRank;
+     MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+     CContext* context=CContext::getCurrent() ;  //printf("myRank = %d, CContext::getCurrent OK\n", myRank);
+
+     //printf("myRank = %d, this->isClientAfterTransformationChecked = %d\n", myRank, this->isClientAfterTransformationChecked);
      if (this->isClientAfterTransformationChecked) return;
+     //printf("myRank = %d, context->hasClient = %d\n", myRank, context->hasClient);
      if (context->hasClient)
      {
-       this->checkMask();
-       if (hasLonLat || hasArea || isCompressible_) this->computeConnectedServer();
-       if (hasLonLat) this->completeLonLatClient();
+       this->checkMask();  //printf("myRank = %d, this->checkMask OK\n", myRank);
+       //printf("myRank = %d, hasLonLat = %d, hasArea = %d, isCompressible_ = %d\n", myRank, hasLonLat, hasArea, isCompressible_);
+       if (hasLonLat || hasArea || isCompressible_) 
+        {
+          //printf("myRank = %d, start this->computeConnectedServer\n", myRank);
+          this->computeConnectedServer();
+          //printf("myRank = %d, this->computeConnectedServer OK\n", myRank);
+        }
+       //printf("myRank = %d, hasLonLat = %d\n", myRank, hasLonLat);
+       if (hasLonLat) 
+        {
+          this->completeLonLatClient();
+          //printf("myRank = %d, this->completeLonLatClient OK\n", myRank);
+        }
      }
 
@@ -1448 +1460 @@
   void CDomain::computeConnectedServer(void)
   {
+    int myRank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
     CContext* context=CContext::getCurrent() ;
     CContextClient* client=context->client ;
@@ -1534 +1549 @@
       }
     }
+    
 
     size_t globalSizeIndex = 1, indexBegin, indexEnd;
@@ -1555 +1571 @@
       indexEnd = indexBegin + range - 1;
     }
-
+    
     CServerDistributionDescription serverDescription(nGlobDomain_, nbServer);
     if (isUnstructed_) serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 0);
     else serverDescription.computeServerGlobalIndexInRange(std::make_pair<size_t,size_t>(indexBegin, indexEnd), 1);
 
+    //printf("myRank = %d, check 7 OK\n", myRank);
+
     CClientServerMapping* clientServerMap = new CClientServerMappingDistributed(serverDescription.getGlobalIndexRange(),
                                                                                 client->intraComm);
-    clientServerMap->computeServerIndexMapping(globalIndexDomain);
+    //printf("myRank = %d new OK\n", myRank);
+
+    clientServerMap->computeServerIndexMapping(globalIndexDomain);  
+    //printf("myRank = %d, clientServerMap->computeServerIndexMapping(globalIndexDomain) OK\n", myRank);
+    
     const CClientServerMapping::GlobalIndexMap& globalIndexDomainOnServer = clientServerMap->getGlobalIndexOnServer();
+    //printf("myRank = %d, clientServerMap->getGlobalIndexOnServer OK\n", myRank);
+    
+    //printf("myRank = %d, check 8 OK\n", myRank);
 
     CClientServerMapping::GlobalIndexMap::const_iterator it  = globalIndexDomainOnServer.begin(),
@@ -1597 +1622 @@
     }
 
+    //printf("myRank = %d, check 9 OK\n", myRank);
+
     connectedServerRank_.clear();
     for (it = globalIndexDomainOnServer.begin(); it != ite; ++it) {
@@ -1605 +1632 @@
 
     delete clientServerMap;
+    //printf("myRank = %d, check 10 OK\n", myRank);
   }
 

XIOS/dev/branch_yushan/src/node/field.cpp

@@ -36 +36 @@
       , domAxisScalarIds_(vector<StdString>(3,"")), areAllReferenceSolved(false), isReferenceSolved(false)
       , useCompressedOutput(false)
-      , hasTimeInstant(false)
-      , hasTimeCentered(false)
-      , wasDataAlreadyReceivedFromServer(false)
-      , isEOF(false)
+      , isReadDataRequestPending(false)
    { setVirtualVariableGroup(CVariableGroup::create(getId() + "_virtual_variable_group")); }
 
@@ -50 +47 @@
       , domAxisScalarIds_(vector<StdString>(3,"")), areAllReferenceSolved(false), isReferenceSolved(false)
       , useCompressedOutput(false)
-      , hasTimeInstant(false)
-      , hasTimeCentered(false)
-      , wasDataAlreadyReceivedFromServer(false)
-      , isEOF(false)
+      , isReadDataRequestPending(false)
    { setVirtualVariableGroup(CVariableGroup::create(getId() + "_virtual_variable_group")); }
 
@@ -253 +247 @@
   }
 
-  bool CField::sendReadDataRequest(const CDate& tsDataRequested)
+  void CField::sendReadDataRequest(void)
   {
     CContext* context = CContext::getCurrent();
     CContextClient* client = context->client;
 
-    lastDataRequestedFromServer = tsDataRequested;
-
-    if (!isEOF) // No need to send the request if we already know we are at EOF
-    {
-      CEventClient event(getType(), EVENT_ID_READ_DATA);
-      if (client->isServerLeader())
-      {
-        CMessage msg;
-        msg << getId();
-        const std::list<int>& ranks = client->getRanksServerLeader();
-        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
-          event.push(*itRank, 1, msg);
-        client->sendEvent(event);
-      }
-      else client->sendEvent(event);
-    }
-    else
-      serverSourceFilter->signalEndOfStream(tsDataRequested);
-
-    return !isEOF;
+    lastDataRequestedFromServer = context->getCalendar()->getCurrentDate();
+    isReadDataRequestPending = true;
+
+    CEventClient event(getType(), EVENT_ID_READ_DATA);
+    if (client->isServerLeader())
+    {
+      CMessage msg;
+      msg << getId();
+      const std::list<int>& ranks = client->getRanksServerLeader();
+      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
+        event.push(*itRank, 1, msg);
+      client->sendEvent(event);
+    }
+    else client->sendEvent(event);
   }
 
@@ -288 +276 @@
     const CDate& currentDate = CContext::getCurrent()->getCalendar()->getCurrentDate();
 
-    bool dataRequested = false;
-
-    while (currentDate >= lastDataRequestedFromServer)
-    {
-      info(20) << "currentDate : " << currentDate << endl ;
-      info(20) << "lastDataRequestedFromServer : " << lastDataRequestedFromServer << endl ;
-      info(20) << "file->output_freq.getValue() : " << file->output_freq.getValue() << endl ;
-      info(20) << "lastDataRequestedFromServer + file->output_freq.getValue() : " << lastDataRequestedFromServer + file->output_freq << endl ;
-
-      dataRequested |= sendReadDataRequest(lastDataRequestedFromServer + file->output_freq);
-    }
-
-    return dataRequested;
+    bool requestData = (currentDate >= lastDataRequestedFromServer + file->output_freq.getValue());
+
+    if (requestData)
+    {
+      cout<<"currentDate : "<<currentDate<<endl ;
+      cout<<"lastDataRequestedFromServer : "<<lastDataRequestedFromServer<<endl ;
+      cout<<"file->output_freq.getValue() : "<<file->output_freq.getValue()<<endl ;
+      cout<<"lastDataRequestedFromServer + file->output_freq.getValue() : "<<lastDataRequestedFromServer + file->output_freq.getValue()<<endl ;
+
+      sendReadDataRequest();
+    }
+
+    return requestData;
   }
 
@@ -322 +310 @@
 
     map<int, CArray<double,1> >::iterator it;
-    if (!grid->doGridHaveDataDistributed())
-    {
-       if (client->isServerLeader())
-       {
-          if (!data_srv.empty())
-          {
-            it = data_srv.begin();
-            const std::list<int>& ranks = client->getRanksServerLeader();
-            for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
-            {
-              msgs.push_back(CMessage());
-              CMessage& msg = msgs.back();
-              msg << getId();
-              if (hasData)
-                msg << getNStep() - 1 << it->second;
-              else
-                msg << int(-1);
-              event.push(*itRank, 1, msg);
-            }
-          }
-          client->sendEvent(event);
-       } 
-       else 
-       {
-          // if (!data_srv.empty())
-          // {
-          //   it = data_srv.begin();
-          //   const std::list<int>& ranks = client->getRanksServerNotLeader();
-          //   for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
-          //   {
-          //     msgs.push_back(CMessage());
-          //     CMessage& msg = msgs.back();
-          //     msg << getId();
-          //     if (hasData)
-          //       msg << getNStep() - 1 << it->second;
-          //     else
-          //       msg << int(-1);
-          //     event.push(*itRank, 1, msg);
-          //   }
-          // }
-          client->sendEvent(event);
-       }
-    }
-    else
-    {
-      for (it = data_srv.begin(); it != data_srv.end(); it++)
-      {
-        msgs.push_back(CMessage());
-        CMessage& msg = msgs.back();
-        msg << getId();
-        if (hasData)
-          msg << getNStep() - 1 << it->second;
-        else
-          msg << int(-1);
-        event.push(it->first, grid->nbSenders[it->first], msg);
-      }
-      client->sendEvent(event);
-    }
+    for (it = data_srv.begin(); it != data_srv.end(); it++)
+    {
+      msgs.push_back(CMessage());
+      CMessage& msg = msgs.back();
+      msg << getId();
+      if (hasData)
+        msg << getNStep() - 1 << it->second;
+      else
+        msg << int(-1);
+      event.push(it->first, grid->nbSenders[it->first], msg);
+    }
+    client->sendEvent(event);
   }
 
@@ -402 +344 @@
         this->incrementNStep();
 
+
+
         if (getNStep() > nstepMax && (getRelFile()->cyclic.isEmpty() || !getRelFile()->cyclic) )
           return false;
@@ -435 +379 @@
     std::map<int, CArray<double,1> > data;
 
+    bool isEOF = false;
+
     for (int i = 0; i < ranks.size(); i++)
     {
@@ -447 +393 @@
     }
 
-    if (wasDataAlreadyReceivedFromServer)
-      lastDataReceivedFromServer = lastDataReceivedFromServer + file->output_freq;
+    if (isEOF)
+      serverSourceFilter->signalEndOfStream(lastDataRequestedFromServer);
     else
-    {
-      lastDataReceivedFromServer = context->getCalendar()->getInitDate();
-      wasDataAlreadyReceivedFromServer = true;
-    }
-
-    if (isEOF)
-      serverSourceFilter->signalEndOfStream(lastDataReceivedFromServer);
-    else
-      serverSourceFilter->streamDataFromServer(lastDataReceivedFromServer, data);
+      serverSourceFilter->streamDataFromServer(lastDataRequestedFromServer, data);
+
+    isReadDataRequestPending = false;
   }
 
@@ -681 +621 @@
    void CField::solveAllReferenceEnabledField(bool doSending2Server)
    {
-     CContext* context = CContext::getCurrent();
-     solveOnlyReferenceEnabledField(doSending2Server);
+     int myRank;
+     MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+     CContext* context = CContext::getCurrent();  //printf("my_Rank = %d, CContext* context = CContext::getCurrent OK\n", myRank);
+     solveOnlyReferenceEnabledField(doSending2Server);  //printf("my_Rank = %d, solveOnlyReferenceEnabledField(doSending2Server) OK\n", myRank);
 
      if (!areAllReferenceSolved)
@@ -690 +633 @@
         if (context->hasClient)
         {
-          solveRefInheritance(true);
-          if (hasDirectFieldReference()) getDirectFieldReference()->solveAllReferenceEnabledField(false);
+          solveRefInheritance(true);  //printf("my_Rank = %d, solveRefInheritance(true) OK\n", myRank);
+          if (hasDirectFieldReference()) 
+          {
+            getDirectFieldReference()->solveAllReferenceEnabledField(false); 
+            //printf("my_Rank = %d, getDirectFieldReference()->solveAllReferenceEnabledField(false) OK\n", myRank);
+          }
         }
         else if (context->hasServer)
+        {
           solveServerOperation();
-
-        solveGridReference();
-     }
-
-     solveGridDomainAxisRef(doSending2Server);
+          //printf("my_Rank = %d, solveServerOperation OK\n", myRank);
+        }
+
+        solveGridReference();  //printf("my_Rank = %d, solveGridReference OK\n", myRank);
+     }
+
+     solveGridDomainAxisRef(doSending2Server);  //printf("my_Rank = %d, solveGridDomainAxisRef(doSending2Server) OK\n", myRank);
 
      if (context->hasClient)
      {
-       solveTransformedGrid();
-     }
-
-     solveCheckMaskIndex(doSending2Server);
+       solveTransformedGrid();  //printf("my_Rank = %d, solveTransformedGrid OK\n", myRank);
+     }
+
+     solveCheckMaskIndex(doSending2Server);  //printf("FIELD.CPP: my_Rank = %d, solveCheckMaskIndex(doSending2Server) OK\n", myRank);
    }
 
@@ -783 +733 @@
      {
        // Check if we have an expression to parse
-       if (hasExpression())
-       {
-         boost::scoped_ptr<IFilterExprNode> expr(parseExpr(getExpression() + '\0'));
-         boost::shared_ptr<COutputPin> filter = expr->reduce(gc, *this);
-
-         // Check if a spatial transformation is needed
-         if (!field_ref.isEmpty())
-         {
-           CGrid* gridRef = CField::get(field_ref)->grid;
-
-           if (grid && grid != gridRef && grid->hasTransform())
-           {
-             bool hasMissingValue = (!detect_missing_value.isEmpty() && !default_value.isEmpty() && detect_missing_value == true);
-             double defaultValue  = hasMissingValue ? default_value : (!default_value.isEmpty() ? default_value : 0.0);
-             std::pair<boost::shared_ptr<CFilter>, boost::shared_ptr<CFilter> > filters = CSpatialTransformFilter::buildFilterGraph(gc, gridRef, grid, hasMissingValue, defaultValue);
-
-             filter->connectOutput(filters.first, 0);
-             filter = filters.second;
-           }
-         }
-
-         instantDataFilter = filter;
+       if (!content.empty())
+       {
+         boost::scoped_ptr<IFilterExprNode> expr(parseExpr(content + '\0'));
+         instantDataFilter = expr->reduce(gc, *this);
        }
        // Check if we have a reference on another field
@@ -811 +743 @@
        // Check if the data is to be read from a file
        else if (file && !file->mode.isEmpty() && file->mode == CFile::mode_attr::read)
-         instantDataFilter = serverSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid,
-                                                                                                     freq_offset.isEmpty() ? NoneDu : freq_offset,
-                                                                                                     true));
+         instantDataFilter = serverSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(grid,
+                                                                                                     freq_offset.isEmpty() ? NoneDu : freq_offset));
        else // The data might be passed from the model
-       {
-          bool ignoreMissingValue = (!detect_missing_value.isEmpty() && !default_value.isEmpty() && detect_missing_value == true);
-          double defaultValue  = ignoreMissingValue ? default_value : (!default_value.isEmpty() ? default_value : 0.0);
-          instantDataFilter = clientSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, NoneDu, false,
-                                                                                                      ignoreMissingValue, defaultValue));       }
+         instantDataFilter = clientSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(grid));
      }
 
@@ -858 +785 @@
      // Check if a spatial transformation is needed
      if (grid && grid != fieldRef->grid && grid->hasTransform())
-     {       
-       bool hasMissingValue = (!detect_missing_value.isEmpty() && !default_value.isEmpty() && detect_missing_value == true);
-       double defaultValue  = hasMissingValue ? default_value : (!default_value.isEmpty() ? default_value : 0.0);                                
-       filters = CSpatialTransformFilter::buildFilterGraph(gc, fieldRef->grid, grid, hasMissingValue, defaultValue);
-     }
+     {
+       double defaultValue = 0.0;
+       if (!default_value.isEmpty()) defaultValue = this->default_value;
+       filters = CSpatialTransformFilter::buildFilterGraph(gc, fieldRef->grid, grid, defaultValue);
+     }
+
      else
        filters.first = filters.second = boost::shared_ptr<CFilter>(new CPassThroughFilter(gc));
@@ -882 +810 @@
    boost::shared_ptr<COutputPin> CField::getSelfReference(CGarbageCollector& gc)
    {
-     if (instantDataFilter || !hasExpression())
+     if (instantDataFilter || content.empty())
        ERROR("COutputPin* CField::getSelfReference(CGarbageCollector& gc)",
              "Impossible to add a self reference to a field which has already been parsed or which does not have an expression.");
@@ -891 +819 @@
        {
          if (!serverSourceFilter)
-           serverSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid,
-                                                                                   freq_offset.isEmpty() ? NoneDu : freq_offset,
-                                                                                   true));
+           serverSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(grid,
+                                                                                   freq_offset.isEmpty() ? NoneDu : freq_offset));
 
          selfReferenceFilter = serverSourceFilter;
        }
        else if (!field_ref.isEmpty())
-       {
-         CField* fieldRef = CField::get(field_ref);
-         fieldRef->buildFilterGraph(gc, false); 
-         selfReferenceFilter = fieldRef->getInstantDataFilter();
-       }
+         selfReferenceFilter = getFieldReference(gc);
        else
        {
          if (!clientSourceFilter)
-         {
-           bool ignoreMissingValue = (!detect_missing_value.isEmpty() && !default_value.isEmpty() && detect_missing_value == true);
-           double defaultValue  = ignoreMissingValue ? default_value : (!default_value.isEmpty() ? default_value : 0.0); 
-           clientSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, NoneDu, false,
-                                                                                   ignoreMissingValue, defaultValue));
-         }
+           clientSourceFilter = boost::shared_ptr<CSourceFilter>(new CSourceFilter(grid));
 
          selfReferenceFilter = clientSourceFilter;
@@ -945 +863 @@
 
        const bool ignoreMissingValue = (!detect_missing_value.isEmpty() && !default_value.isEmpty() && detect_missing_value == true);
-       
+
        boost::shared_ptr<CTemporalFilter> temporalFilter(new CTemporalFilter(gc, operation,
                                                                              CContext::getCurrent()->getCalendar()->getInitDate(),
@@ -957 +875 @@
      return it->second;
    }
-
-  /*!
-    * Returns the temporal filter corresponding to the field's temporal operation
-    * for the specified operation frequency.
-    *
-    * \param gc the garbage collector to use
-    * \param outFreq the operation frequency, i.e. the frequency at which the output data will be computed
-    * \return the output pin corresponding to the requested temporal filter
-    */
-   
-   boost::shared_ptr<COutputPin> CField::getSelfTemporalDataFilter(CGarbageCollector& gc, CDuration outFreq)
-   {
-     if (instantDataFilter || !hasExpression())
-       ERROR("COutputPin* CField::getSelfTemporalDataFilter(CGarbageCollector& gc)",
-             "Impossible to add a self reference to a field which has already been parsed or which does not have an expression.");
-
-     if (!selfReferenceFilter) getSelfReference(gc) ;
-
-     if (serverSourceFilter || clientSourceFilter)
-     {
-       if (operation.isEmpty())
-         ERROR("void CField::getSelfTemporalDataFilter(CGarbageCollector& gc, CDuration outFreq)",
-               << "An operation must be defined for field \"" << getId() << "\".");
-
-       if (freq_op.isEmpty()) freq_op.setValue(TimeStep);
-       if (freq_offset.isEmpty()) freq_offset.setValue(NoneDu);
-
-       const bool ignoreMissingValue = (!detect_missing_value.isEmpty() && !default_value.isEmpty() && detect_missing_value == true);
-
-       boost::shared_ptr<CTemporalFilter> temporalFilter(new CTemporalFilter(gc, operation,
-                                                                             CContext::getCurrent()->getCalendar()->getInitDate(),
-                                                                             freq_op, freq_offset, outFreq,
-                                                                             ignoreMissingValue, ignoreMissingValue ? default_value : 0.0));
-       selfReferenceFilter->connectOutput(temporalFilter, 0);
-       return temporalFilter ;
-     }
-     else if (!field_ref.isEmpty())
-     {
-       CField* fieldRef = CField::get(field_ref);
-       fieldRef->buildFilterGraph(gc, false); 
-       return fieldRef->getTemporalDataFilter(gc, outFreq) ;
-     }
-  }
 
    //----------------------------------------------------------------
@@ -1153 +1028 @@
      else if (grid && grid->hasTransform() && !grid->isTransformed())
      {
-       // Temporarily deactivate the self-transformation of grid
-       //grid->transformGrid(grid);
+       grid->transformGrid(grid);
      }
    }
@@ -1430 +1304 @@
    }
 
-   /*!
-    * Returns string arithmetic expression associated to the field.
-    * \return if content is defined return content string, otherwise, if "expr" attribute is defined, return expr string.
-    */
-   const string& CField::getExpression(void)
-   {
-     if (!expr.isEmpty() && content.empty())
-     {
-       content = expr;
-       expr.reset();
-     }
-
-     return content;
-   }
-
-   bool CField::hasExpression(void) const
-   {
-     return (!expr.isEmpty() || !content.empty());
-   }
-
    DEFINE_REF_FUNC(Field,field)
 } // namespace xios

XIOS/dev/branch_yushan/src/node/field.hpp

@@ -124 +124 @@
          boost::shared_ptr<COutputPin> getSelfReference(CGarbageCollector& gc);
          boost::shared_ptr<COutputPin> getTemporalDataFilter(CGarbageCollector& gc, CDuration outFreq);
-         boost::shared_ptr<COutputPin> getSelfTemporalDataFilter(CGarbageCollector& gc, CDuration outFreq);
 
 //         virtual void fromBinary(StdIStream& is);
@@ -143 +142 @@
         void recvUpdateData(vector<int>& ranks, vector<CBufferIn*>& buffers);
         void writeField(void);
-        bool sendReadDataRequest(const CDate& tsDataRequested);
+        void sendReadDataRequest(void);
         bool sendReadDataRequestIfNeeded(void);
         static void recvReadDataRequest(CEventServer& event);
@@ -179 +178 @@
         const std::vector<StdString>& getRefDomainAxisIds();
 
-        const string& getExpression(void);
-        bool hasExpression(void) const;
-
       public:
          /// Propriétés privées ///
@@ -195 +191 @@
          bool isEOF;
          CDate lastlast_Write_srv, last_Write_srv, last_operation_srv;
-         CDate lastDataRequestedFromServer, lastDataReceivedFromServer;
-         bool wasDataAlreadyReceivedFromServer;
+         CDate lastDataRequestedFromServer;
 
          map<int,boost::shared_ptr<func::CFunctor> > foperation_srv;
@@ -206 +201 @@
          bool isReferenceSolved;
          std::vector<StdString> domAxisScalarIds_;
+         bool isReadDataRequestPending;
          bool useCompressedOutput;
-
-         // Two variable to identify the time_counter meta data written in file, which has no time_counter
-         bool hasTimeInstant;
-         bool hasTimeCentered;
 
          DECLARE_REF_FUNC(Field,field)

XIOS/dev/branch_yushan/src/node/field_impl.hpp

@@ -19 +19 @@
   {
     if (clientSourceFilter)
+    {
       clientSourceFilter->streamData(CContext::getCurrent()->getCalendar()->getCurrentDate(), _data);
-    else if (instantDataFilter)
+    }
+    else if (!field_ref.isEmpty() || !content.empty())
       ERROR("void CField::setData(const CArray<double, N>& _data)",
             << "Impossible to receive data from the model for a field [ id = " << getId() << " ] with a reference or an arithmetic operation.");

XIOS/dev/branch_yushan/src/node/file.cpp

@@ -17 +17 @@
 #include "mpi.hpp"
 
+
 namespace xios {
 
@@ -29 +30 @@
      setVirtualVariableGroup(CVariableGroup::create(getId() + "_virtual_variable_group"));
    }
+   
+
 
    CFile::CFile(const StdString & id)
@@ -38 +41 @@
       setVirtualVariableGroup(CVariableGroup::create(getId() + "_virtual_variable_group"));
     }
+   
 
    CFile::~CFile(void)
@@ -235 +239 @@
       const int recordOffset = record_offset.isEmpty() ? 0 : record_offset;
 
-      // set<CAxis*> setAxis;
-      // set<CDomain*> setDomains;
-      set<StdString> setAxis;
-      set<StdString> setDomains;
-      
+      set<CAxis*> setAxis;
+      set<CDomain*> setDomains;
+
       std::vector<CField*>::iterator it, end = this->enabledFields.end();
       for (it = this->enabledFields.begin(); it != end; it++)
@@ -247 +249 @@
          std::vector<CAxis*> vecAxis = field->grid->getAxis();
          for (size_t i = 0; i < vecAxis.size(); ++i)
-            setAxis.insert(vecAxis[i]->getAxisOutputName());
-            // setAxis.insert(vecAxis[i]);
+            setAxis.insert(vecAxis[i]);
          std::vector<CDomain*> vecDomains = field->grid->getDomains();
          for (size_t i = 0; i < vecDomains.size(); ++i)
-            setDomains.insert(vecDomains[i]->getDomainOutputName());
-            // setDomains.insert(vecDomains[i]);
+            setDomains.insert(vecDomains[i]);
 
          field->resetNStep(recordOffset);
@@ -264 +264 @@
       if (allDomainEmpty) MPI_Comm_free(&fileComm);
 
-      // if (time_counter.isEmpty()) time_counter.setValue(time_counter_attr::centered);
+      if (time_counter.isEmpty()) time_counter.setValue(time_counter_attr::centered);
       if (time_counter_name.isEmpty()) time_counter_name = "time_counter";
     }
@@ -375 +375 @@
          if (pos2!=std::string::npos)
          {
-           middlePart=filename.substr(pos1,pos2-pos1) ;           
+           middlePart=filename.substr(pos1,pos2-pos1) ;
+           cout<<pos2<<endl ;
            pos2+=strEndDate.size() ;
            lastPart=filename.substr(pos2,filename.size()-pos2) ;
@@ -580 +581 @@
           this->data_in->closeFile();
        }
+
       if (fileComm != MPI_COMM_NULL) MPI_Comm_free(&fileComm);
+
    }
    //----------------------------------------------------------------
@@ -715 +718 @@
    {
      int size = this->enabledFields.size();
+     
      for (int i = 0; i < size; ++i)
      {
@@ -745 +749 @@
      int size = this->enabledFields.size();
      for (int i = 0; i < size; ++i)
-       this->enabledFields[i]->sendReadDataRequest(CContext::getCurrent()->getCalendar()->getCurrentDate());
+       this->enabledFields[i]->sendReadDataRequest();
    }
 

XIOS/dev/branch_yushan/src/node/file.hpp

@@ -12 +12 @@
 #include "attribute_enum_impl.hpp"
 #include "mpi.hpp"
+#ifdef _usingEP
+#include "ep_declaration.hpp"
+#endif
 
 namespace xios {

XIOS/dev/branch_yushan/src/node/grid.cpp

@@ -166 +166 @@
    void CGrid::checkAttributesAfterTransformation()
    {
-      setAxisList();
-      std::vector<CAxis*> axisListP = this->getAxis();
+      int myRank;
+      MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+      setAxisList();  //printf("myRank = %d, setAxisList OK\n", myRank);
+      std::vector<CAxis*> axisListP = this->getAxis();  //printf("myRank = %d, this->getAxis OK\n", myRank);
       if (!axisListP.empty())
       {
         int idx = 0;
         axisPositionInGrid_.resize(0);
+        
         for (int i = 0; i < axis_domain_order.numElements(); ++i)
         {
@@ -179 +183 @@
             axisPositionInGrid_.push_back(idx);
             ++idx;
+            //printf("myRank = %d, axisPositionInGrid_.push_back OK\n", myRank);
           }
           else if (2 == elementDimension) idx += 2;
@@ -186 +191 @@
         {
           axisListP[i]->checkAttributesOnClientAfterTransformation(globalDim_,axisPositionInGrid_[i]);
-        }
-      }
-
-      setDomainList();
-      std::vector<CDomain*> domListP = this->getDomains();
+          //printf("myRank = %d, axisListP[%d/%d]->checkAttributesOnClientAfterTransformation OK\n", myRank, i, axisListP.size());
+        }
+      }
+
+      setDomainList(); //printf("myRank = %d, setDomainList OK\n", myRank);
+      std::vector<CDomain*> domListP = this->getDomains(); //printf("myRank = %d, this->getDomains OK\n", myRank);
       if (!domListP.empty())
       {
         for (int i = 0; i < domListP.size(); ++i)
         {
+          //printf("myRank = %d, start domListP[%d]->checkAttributesOnClientAfterTransformation\n", myRank, i);
           domListP[i]->checkAttributesOnClientAfterTransformation();
+          //printf("myRank = %d, domListP[%d]->checkAttributesOnClientAfterTransformation OK\n", myRank, i);
         }
       }
@@ -272 +280 @@
    void CGrid::checkMaskIndex(bool doSendingIndex)
    {
-     CContext* context = CContext::getCurrent();
+     int myRank;
+     MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+     CContext* context = CContext::getCurrent(); 
      CContextClient* client=context->client;
 
@@ -278 +289 @@
      {
        if (context->hasClient)
-          if (this->isChecked && doSendingIndex && !isIndexSent) { sendIndexScalarGrid(); this->isIndexSent = true; }
+          if (this->isChecked && doSendingIndex && !isIndexSent) 
+            { 
+              sendIndexScalarGrid(); 
+              this->isIndexSent = true; 
+            }
 
        if (this->isChecked) return;
-
        if (context->hasClient)
        {
@@ -287 +301 @@
        }
 
-       if (!(this->hasTransform() && !this->isTransformed()))
-        this->isChecked = true;
+       this->isChecked = true;
        return;
      }
-
+     
      if (context->hasClient)
-      if (this->isChecked && doSendingIndex && !isIndexSent) { sendIndex(); this->isIndexSent = true; }
+      if (this->isChecked && doSendingIndex && !isIndexSent) 
+        {
+          sendIndex(); 
+          this->isIndexSent = true;
+        }
 
      if (this->isChecked) return;
@@ -299 +316 @@
      if (context->hasClient)
      {
-        this->checkAttributesAfterTransformation();
-        this->checkMask();
-        this->computeIndex();
+        this->checkAttributesAfterTransformation(); 
+        this->checkMask();  
+        this->computeIndex();  
      }
-
-     if (!(this->hasTransform() && !this->isTransformed())) 
-      this->isChecked = true;
-
-     if (!(this->hasTransform() && (!this->isGenerated())))
-      this->isChecked = true; 
+     this->isChecked = true;
    }
 
@@ -1051 +1063 @@
         int rank = *itRank;
         int nb = 1;
-        storeIndex_toSrv.insert(std::make_pair(rank, CArray<int,1>(nb)));        
+        storeIndex_toSrv.insert(std::make_pair(rank, CArray<int,1>(nb)));
         listOutIndex.push_back(CArray<size_t,1>(nb));
 
@@ -1063 +1075 @@
         }
 
-        storeIndex_fromSrv.insert(std::make_pair(rank, CArray<int,1>(outLocalIndexToServer)));
         listMsg.push_back(CMessage());
         listMsg.back() << getId( )<< isDataDistributed_ << isCompressible_ << listOutIndex.back();
@@ -1072 +1083 @@
     }
     else
-    {
-      const std::list<int>& ranks = client->getRanksServerNotLeader();
-      for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
-      {
-        int rank = *itRank;
-        int nb = 1;        
-        storeIndex_fromSrv.insert(std::make_pair(rank, CArray<int,1>(nb)));
-        CArray<int, 1>& outLocalIndexToServer = storeIndex_fromSrv[rank];
-        for (int k = 0; k < nb; ++k)
-        {          
-          outLocalIndexToServer(k)  = 0;
-        }
-      }
       client->sendEvent(event);
-    }
   }
 
@@ -1114 +1111 @@
           outLocalIndexToServer(idx) = itIndex->second;
         }
-        
+
+        //int nbClient = client->clientSize; // This stupid variable signals the servers the number of client connect to them
         const std::list<int>& ranks = client->getRanksServerLeader();
         for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
         {
           storeIndex_toSrv.insert(std::make_pair(*itRank, CArray<int,1>(outLocalIndexToServer)));
-          storeIndex_fromSrv.insert(std::make_pair(*itRank, CArray<int,1>(outLocalIndexToServer)));
           listOutIndex.push_back(CArray<size_t,1>(outGlobalIndexOnServer));
 
@@ -1129 +1126 @@
         client->sendEvent(event);
       }
-      else 
-      {
-        int indexSize = globalLocalIndexSendToServer.size();        
-        CArray<int,1> outLocalIndexToServer(indexSize);
-        for (int idx = 0; itIndex != iteIndex; ++itIndex, ++idx)
-        {          
-          outLocalIndexToServer(idx) = itIndex->second;
-        }
-        
-        const std::list<int>& ranks = client->getRanksServerNotLeader();
-        for (std::list<int>::const_iterator itRank = ranks.begin(), itRankEnd = ranks.end(); itRank != itRankEnd; ++itRank)
-        {          
-          storeIndex_fromSrv.insert(std::make_pair(*itRank, CArray<int,1>(outLocalIndexToServer)));
-        }
+      else
         client->sendEvent(event);
-      }
     }
     else
@@ -1177 +1160 @@
           nb = globalIndexTmp[rank].size();
 
-        storeIndex_toSrv.insert(make_pair(rank, CArray<int,1>(nb)));        
+        storeIndex_toSrv.insert(make_pair(rank, CArray<int,1>(nb)));
         listOutIndex.push_back(CArray<size_t,1>(nb));
 
@@ -1189 +1172 @@
         }
 
-        storeIndex_fromSrv.insert(make_pair(rank, CArray<int,1>(outLocalIndexToServer)));
         listMsg.push_back(CMessage());
         listMsg.back() << getId() << isDataDistributed_ << isCompressible_ << listOutIndex.back();

XIOS/dev/branch_yushan/src/node/grid.hpp

@@ -221 +221 @@
 
          map<int, CArray<int, 1> > storeIndex_toSrv;
-         map<int, CArray<int, 1> > storeIndex_fromSrv;
          map<int,int> nbSenders;
 

XIOS/dev/branch_yushan/src/node/interpolate_domain.cpp

@@ -48 +48 @@
     }
 
-    if (this->mode.isEmpty()) this->mode.setValue(mode_attr::compute);
-    if (this->write_weight.isEmpty()) this->write_weight.setValue(false);
+    StdString weightFile = "interpolation_weights_" + domainSrc->getDomainOutputName();
 
-    StdString weightFile;
-    switch (this->mode)
-    {
-      case mode_attr::read:
-        if (this->weight_filename.isEmpty())
-        {
-          if (!this->write_weight)
+    if (!this->mode.isEmpty())
+    { 
+      if (mode_attr::read == this->mode)
+      {
+         if (this->file.isEmpty())
+         {
             ERROR("void CInterpolateDomain::checkValid(CDomain* domainSrc)",
                  << "Read mode is activated but there is no file specified." << std::endl
                  << "Please define a correct file containing interpolation weights with option 'file'. ");
-        }
-        else
-        {
-          weightFile = this->weight_filename;
-          ifstream f(weightFile.c_str());
-          if (!f.good())
-            ERROR("void CInterpolateDomain::checkValid(CDomain* domainSrc)",
-                  << "Read mode is activated but file "  << weightFile << " doesn't exist." << std::endl
-                  << "Please check this file ");
-        }
-        break;
-      case mode_attr::compute:
-        break;
-      case mode_attr::read_or_compute:
-        if (!this->weight_filename.isEmpty() && !this->write_weight)
-        {
-          weightFile = this->weight_filename;
-          ifstream f(weightFile.c_str());
-          if (!f.good())
-            ERROR("void CInterpolateDomain::checkValid(CDomain* domainSrc)",
-                  << "read_or_compute mode is activated but file "  << weightFile << " doesn't exist." << std::endl
-                  << "Please check this file ");
-        }
-        break;
-      default:
-        break;
+         }
+         else
+         {
+           weightFile = this->file;
+           ifstream f(weightFile.c_str());
+           if (!f.good())
+             ERROR("void CInterpolateDomain::checkValid(CDomain* domainSrc)",
+                   << "Read mode is activated but file "  << weightFile << " doesn't exist." << std::endl
+                   << "Please check this file ");
+         }
+      }
+      else
+      {
+        if (file.isEmpty())
+          this->file.setValue(weightFile);
+      }   
+    }
+    else
+    {
+      if (!file.isEmpty()) // set mode read
+      {
+         weightFile = this->file;
+         ifstream f(weightFile.c_str());
+         if (!f.good()) // file doesn't exist
+           this->mode.setValue(mode_attr::write);
+         else
+           this->mode.setValue(mode_attr::read);
+      }
+      else // only calculate weights() Mode set write but there is no file)
+      {
+        this->mode.setValue(mode_attr::write);
+      }
     }
 

XIOS/dev/branch_yushan/src/node/interpolate_domain.hpp

@@ -50 +50 @@
       virtual ~CInterpolateDomain(void);
 
-      virtual void checkValid(CDomain* domainSource);
+      virtual void checkValid(CDomain* axisDest);
 
       /// Accesseurs statiques ///

XIOS/dev/branch_yushan/src/node/mesh.cpp

@@ -11 +11 @@
 /// ////////////////////// Définitions ////////////////////// ///
 
-  CMesh::CMesh(void) :  nbNodesGlo(0), nbEdgesGlo(0)
-            ,  node_start(0), node_count(0)
-            ,  edge_start(0), edge_count(0)
-            ,  nbFaces_(0), nbNodes_(0), nbEdges_(0)
-            ,  nodesAreWritten(false), edgesAreWritten(false), facesAreWritten(false)
+  CMesh::CMesh(void) :  nbNodesGlo{0}, nbEdgesGlo{0}
+            ,  node_start{0}, node_count{0}
+            ,  edge_start{0}, edge_count{0}
+            ,  nbFaces_{0}, nbNodes_{0}, nbEdges_{0}
+            ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false}
             ,  node_lon(), node_lat()
             ,  edge_lon(), edge_lat(), edge_nodes()
             ,  face_lon(), face_lat()
             ,  face_nodes()
-            ,  pNodeGlobalIndex(NULL), pEdgeGlobalIndex(NULL)
+            ,  pNodeGlobalIndex{NULL}, pEdgeGlobalIndex{NULL}
   {
   }
@@ -91 +91 @@
       seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
     }
+    return seed ;
+  }
+
+///----------------------------------------------------------------
+/*!
+ * \fn size_t generateEdgeIndex(size_t first, size_t second, int rank)
+ * Generates an edge index.
+ * If the same edge is generated by two processes, each process will have its own edge index.
+ * \param [in] first Edge node.
+ * \param [in] second Edge node.
+ * \param [in] rank MPI process rank.
+ */
+  size_t generateEdgeIndex(size_t first, size_t second, int rank)
+  {
+    size_t seed = rank ;
+    if (first < second)
+    {
+      seed = hashPair(seed, first);
+      seed = hashPair(seed, second);
+    }
+    else
+    {
+      seed = hashPair(seed, second);
+      seed = hashPair(seed, first);
+    }
+
     return seed ;
   }
@@ -115 +141 @@
     return seed ;
   }
-
-  ///----------------------------------------------------------------
-  /*!
-   * \fn size_t generateNodeIndex(vector<size_t>& valList)
-   * Generates a node index unique for all processes.
-   * \param [in] valList Vector storing four node hashes.
-   */
-    size_t generateNodeIndex(vector<size_t>& valList)
-    {
-      // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
-      vector<size_t> vec = valList;
-      sort (vec.begin(), vec.end());
-      size_t seed = vec[0] ;
-      int it = 1;
-      for(; it != vec.size(); ++it)
-      {
-         seed = hashPair(seed, vec[it]);
-      }
-      return seed ;
-    }
-
 
 ///----------------------------------------------------------------
@@ -595 +600 @@
         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
         CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4);
-        int nbHash = 0;
         for (int ne = 0; ne < nbEdges_; ++ne)
         {
@@ -605 +609 @@
               if (nodeHash2Idx[hashValues[nh]].size() == 0)
               {
-                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues));
+                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues, mpiRank));
                 nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
-                nodeHashList(nbHash) = hashValues[nh];
-                ++nbHash;
               }
-            }
-          }
-        }
-        nodeHashList.resizeAndPreserve(nbHash);
+              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
+            }
+          }
+        }
 
         // (2.2) Generating global node indexes
-        // The ownership criterion: priority of the process of smaller index
+        // The ownership criterion: priority of the process holding the smaller index
         // Maps generated in this step are:
-        // Maps generated in this step are:
-         // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
-         // nodeIdx2Idx = <idx, <rankOwner, idx>>
-
-         CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
-         dhtNodeHash.computeIndexInfoMapping(nodeHashList);
-         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
-
-
-         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
-         CArray<size_t,1> nodeIdxList(nbEdges_*nvertex*4);
-         size_t nIdx = 0;
-
-         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
-         {
-           size_t rankMin = (it->second)[1];
-           size_t idx = (it->second)[0];
-           for (int i = 2; i < (it->second).size();)
-           {
-             if ( (it->second)[i+1] < rankMin)
-             {
-               idx = (it->second)[i];
-               rankMin = (it->second)[i+1];
-               (it->second)[i+1] = (it->second)[i-1];
-             }
-             i += 2;
-           }
-           if (nodeIdx2Idx.count(idx) == 0)
-           {
-             if (mpiRank == rankMin)
-             {
-               nodeIdx2Idx[idx].push_back(rankMin);
-               nodeIdx2Idx[idx].push_back(idx);
-             }
-             nodeIdxList(nIdx) = idx;
-             ++nIdx;
-           }
-         }
-         nodeIdxList.resizeAndPreserve(nIdx);
-
-         // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
-         // Solution: global node indexing by hand.
-         // Maps modified in this step:
-         // nodeIdx2Idx = <idx, idxGlo>
-         int nodeCount = nodeIdx2Idx.size();
-         int nodeStart, nbNodes;
-         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
-         int nNodes = nodeStart;
-         MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
-         nbNodesGlo = nNodes;
-
-         nodeStart -= nodeCount;
-         node_start = nodeStart;
-         node_count = nodeCount;
-         CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
-         size_t count = 0;
-
-         for (int ne = 0; ne < nbEdges_; ++ne)
-         {
-           for (int nv = 0; nv < nvertex; ++nv)
-           {
-             vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
-             size_t nodeIdx = generateNodeIndex(hashValues);
-             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
-             if (it != nodeIdx2Idx.end())
-             {
-               if (dummyMap.count(nodeIdx) == 0)
-               {
-                 dummyMap[nodeIdx].push_back(nodeIdx);
-                 (it->second)[1] = node_start + count;
-                 ++count;
-               }
-             }
-           }
-         }
-
-         CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
-         dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
-         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
+        // nodeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
+        // nodeIdx2IdxMin = <idx, idxMin>
+        // nodeIdx2IdxGlo = <idxMin, idxMin>
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
+        CArray<size_t,1> nodeIdxMinList(nbEdges_*nvertex);
+
+        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
+        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
+        size_t iIdxMin = 0;
+
+        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
+        {
+          size_t idx = (it->second)[0];
+          size_t idxMin = (it->second)[0];
+          for (int i = 2; i < (it->second).size();)
+          {
+            if (mpiRank == (it->second)[i+1])
+            {
+              idx = (it->second)[i];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
+            }
+            i += 2;
+          }
+          (it->second)[0] = idxMin;
+          if (nodeIdx2IdxMin.count(idx) == 0)
+          {
+            nodeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
+            nodeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+        nodeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
+        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
+        // nodeIdx2IdxGlo holds global indexes only for nodes owned by a process
+        // nodeIdxMin2IdxGlo holds global indexes for all nodes generated by a process
 
         // (2.3) Saving variables: node_lon, node_lat, edge_nodes
         // Creating map nodeHash2IdxGlo <hash, idxGlo>
         // Creating map edgeHash2IdxGlo <hash, idxGlo>
-//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
-//        node_count = dhtNodeIdxGlo.getIndexCount();
-//        node_start = dhtNodeIdxGlo.getIndexStart();
+        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
+        node_count = dhtNodeIdxGlo.getIndexCount();
+        node_start = dhtNodeIdxGlo.getIndexStart();
         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
         node_lon.resize(node_count);
@@ -715 +685 @@
           {
             hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
-            size_t myIdx = generateNodeIndex(hashValues);
-            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxGlo.find(myIdx);
-            idxGlo = (itIdx->second)[1];
-
-            if (mpiRank == (itIdx->second)[0])
-            {
+            size_t myIdx = generateNodeIndex(hashValues, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdx2IdxGlo.find(myIdx);
+              idxGlo = (itIdxGlo->second)[0];
 //              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
               node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
               node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
             }
+            else
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdxMin2IdxGlo.find(ownerIdx);
+              idxGlo = (itIdxGlo->second)[0];
+            }
+
             edge_nodes(nv,ne) = idxGlo;
             for (int nh = 0; nh < 4; ++nh)
@@ -815 +793 @@
         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
-        CArray<size_t,1> edgeIdxList(nbFaces_*nvertex);
+        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
         size_t iIdx = 0;
 
@@ -846 +824 @@
               if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
               {
-                edgeIdxList(iIdx) = edgeIdxGlo;
+                edgeIdxGloList(iIdx) = edgeIdxGlo;
                 ++iIdx;
               }
               edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
+              edgeHash2Rank[edgeHash].push_back(itEdgeHash->first);
               edgeHash2Rank[edgeHash].push_back(mpiRank);
-              edgeHash2Rank[edgeHash].push_back(itEdgeHash->second[0]);
             }
             else
@@ -859 +837 @@
           }
         }
-        edgeIdxList.resizeAndPreserve(iIdx);
+        edgeIdxGloList.resizeAndPreserve(iIdx);
 
         // (1.3) Saving remaining variables edge_faces and face_faces
@@ -868 +846 @@
         dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
         CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
-
-        // edgeHash2Info = <edgeHash, < rank1, idxGlo, rank2, idxGlo>>
-        int edgeCount = 0;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo;
+
+        // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>>
+        // edgeHashMin2IdxGlo = <edgeHashMin, idxGlo>
         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
         {
           vector <size_t> edgeInfo = it->second;
-          if (edgeInfo[0] == mpiRank)
-          {
-            ++edgeCount;
-          }
-        }
-
-        int edgeStart, nbEdges;
-        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
-        int nEdges = edgeStart;
-        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
-        nbEdgesGlo = nEdges;
-
-        // edges to be splitted equally between procs
-        if ( (nbEdgesGlo % mpiSize) == 0)
-        {
-          edge_count = nbEdgesGlo/mpiSize;
-          edge_start = mpiRank*edge_count;
-        }
-        else
-        {
-          if (mpiRank == (mpiSize - 1) )
-          {
-            edge_count = nbEdgesGlo/mpiSize;
-            edge_start = mpiRank*(nbEdgesGlo/mpiSize + 1);
-          }
-          else
-          {
-            edge_count = nbEdgesGlo/mpiSize + 1;
-            edge_start = mpiRank*edge_count;
-          }
-        }
-        CArray<size_t,1> edgeIdxGloList(edge_count);
-        for (int i = 0; i < edge_count; ++i)
-        {
-          edgeIdxGloList(i) = i + edge_start;
-        }
-
-        CClientClientDHTSizet dhtEdgeIdxGlo2Face (edgeIdxGlo2Face, comm);
+          if (edgeInfo.size() == 4)                       // two processes generate the same edge
+            if (edgeInfo[1] > edgeInfo[3])
+              edgeInfo[1] = edgeInfo[3];
+          if (edgeInfo[1] == mpiRank)
+            edgeHashMin2IdxGlo[it->first].push_back(edgeInfo[0]);
+        }
+
+        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeHashMin2IdxGlo, comm);
+        edge_count = dhtEdgeIdxGlo.getIndexCount();
+        edge_start = dhtEdgeIdxGlo.getIndexStart();
+
+        edge_faces.resize(2, edge_count);
+        face_faces.resize(nvertex, nbFaces_);
+
         CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
-        dhtEdgeIdxGlo2Face.computeIndexInfoMapping(edgeIdxGloList);
-        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdgeIdxGlo2Face.getInfoIndexMap();
-        dhtEdge2Face.computeIndexInfoMapping(edgeIdxList);
-        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2FaceIdx = dhtEdge2Face.getInfoIndexMap();
-
-
-        edge_faces.resize(2, edge_count);
-        for (int i = 0; i < edge_count; ++i)
-        {
-          CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdxGlo2FaceIdx.find(i + edge_start);
-          int indexGlo = it->first;
-          vector<size_t> faces = it->second;
-          int face1 = faces[0];
-          edge_faces(0, indexGlo - edge_start) = face1;
-          if (faces.size() == 2)
-          {
-            int face2 = faces[1];
-            edge_faces(1, indexGlo - edge_start) = face2;
-          }
-          else
-          {
-            edge_faces(1, indexGlo - edge_start) = -999;
-          }
-        }
-
-        size_t tmp;
-        vector<size_t> tmpVec;
-        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2FaceIdx.begin(); it != edgeIdx2FaceIdx.end(); it++)
-        {
-          tmp = it->first;
-          tmpVec = it->second;
-          tmp++;
-        }
-
-        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2, itIndex;
-        face_faces.resize(nvertex, nbFaces_);
+        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2;
+
         for (int nf = 0; nf < nbFaces_; ++nf)
         {
@@ -973 +896 @@
               size_t faceIdxGlo = nbFacesAccum + nf;
               size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
-              itEdgeHash = edgeHash2Info.find(edgeHash);
-              int edgeIdxGlo = (itEdgeHash->second)[1];
-
-              if ( (itEdgeHash->second)[0] == mpiRank)
+              itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
+              if (itEdgeHash != edgeHashMin2IdxGlo.end())
               {
-                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
+                int edgeIdxGlo = itEdgeHash->second[0];
+                itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                 int face1 = itFace1->second[0];
                 if (itFace1->second.size() == 1)
                 {
+                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                  edge_faces(1, edgeIdxGlo - edge_start) = -999;
                   face_faces(nv1, nf) = 999999;
                 }
@@ -987 +911 @@
                 {
                   int face2 = itFace1->second[1];
+                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                  edge_faces(1, edgeIdxGlo - edge_start) = face2;
                   face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
                 }
@@ -992 +918 @@
               else
               {
-                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
+                itEdgeHash = edgeHashMin2IdxGlo.find(edgeHash);
+                size_t edgeIdxGlo = itEdgeHash->second[0];
+                itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                 int face1 = itFace1->second[0];
                 int face2 = itFace1->second[1];
@@ -1029 +957 @@
         CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
         CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
-        int nEdgeHash = 0;
         for (int nf = 0; nf < nbFaces_; ++nf)
         {
@@ -1051 +978 @@
             face_nodes(nv1,nf) = it1->second[0];
             size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
-            if (edgeHash2Idx.count(edgeHash) == 0)
-            {
-              edgeHash2Idx[edgeHash].push_back(edgeHash);
-              edgeHash2Idx[edgeHash].push_back(mpiRank);
-              edgeHashList(nEdgeHash) = edgeHash;
-              ++nEdgeHash;
-            }
-          }
-        }
-        edgeHashList.resizeAndPreserve(nEdgeHash);
-
-        // (2.3) Generating global edge indexes
-        // The ownership criterion: priority of the process with smaller rank
+            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(it1->second[0], it2->second[0], mpiRank));
+            edgeHash2Idx[edgeHash].push_back(mpiRank);
+            edgeHashList(nf*nvertex + nv1) = edgeHash;
+          }
+        }
+
+        // (2.2) Generating global edge indexes
+        // The ownership criterion: priority of the process holding the smaller index
         // Maps generated in this step are:
-        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
-        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
+        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
+        // edgeIdx2IdxMin = = <idx, idxMin>
+        // edgeIdx2IdxGlo = <idxMin, idxGlo>
 
         CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
         dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
         CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
-        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
-
-        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
+        CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex);
+        size_t iIdxMin = 0;
 
         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
         {
-          size_t rankMin = (it->second)[1];
+          size_t idxMin = (it->second)[0];
           size_t idx = (it->second)[0];
-
           for (int i = 2; i < (it->second).size();)
           {
-            if ((it->second)[i+1] < rankMin)
-            {
-              rankMin = (it->second)[i+1];
+            if (mpiRank == (it->second)[i+1])
+            {
               idx = (it->second)[i];
-              (it->second)[i+1] = (it->second)[i-1];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
             }
             i += 2;
           }
-          if (edgeIdx2Idx.count(idx) == 0)
-          {
-            if (mpiRank == rankMin)
-            {
-              edgeIdx2Idx[idx].push_back(rankMin);
-              edgeIdx2Idx[idx].push_back(idx);
-            }
-          }
-        }
-
-        int edgeCount = edgeIdx2Idx.size();
-        int edgeStart, nbEdges;
-        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
-        int nEdges = edgeStart;
-        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
-        nbEdgesGlo = nEdges;
-
-        edgeStart -= edgeCount;
-        edge_start = edgeStart;
-        edge_count = edgeCount;
-        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
-        int count = 0;
+          (it->second)[0] = idxMin;
+          if (edgeIdx2IdxMin.count(idx) == 0)
+          {
+            edgeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
+            edgeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+        edgeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
+        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
+        // edgeIdx2IdxGlo holds global indexes only for edges owned by a process
+        // edgeIdxMin2IdxGlo holds global indexes for all edges generated by a process
+
+        // (2.3) Saving variables: edge_lon, edge_lat, face_edges
+        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
+        edge_count = dhtEdgeIdxGlo.getIndexCount();
+        edge_start = dhtEdgeIdxGlo.getIndexStart();
+        edge_lon.resize(edge_count);
+        edge_lat.resize(edge_count);
+        edge_nodes.resize(2, edge_count);
+        face_edges.resize(nvertex, nbFaces_);
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
+        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
+        size_t iIdx = 0;
 
         for (int nf = 0; nf < nbFaces_; ++nf)
@@ -1133 +1068 @@
               it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
             }
-            size_t nodeIdxGlo1 = it1->second[0];
-            size_t nodeIdxGlo2 = it2->second[0];
-
-            if (nodeIdxGlo1 != nodeIdxGlo2)
-            {
-              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
-              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
-              if (it != edgeIdx2Idx.end())
-              {
-                if (dummyEdgeMap.count(edgeIdx) == 0)
-                {
-                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
-                  (it->second)[1] = edge_start + count;
-                  ++count;
-                }
-              }
-            }
-          }
-        }
-
-        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
-        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
-        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
-
-        // (2.4) Saving variables: edge_lon, edge_lat, face_edges
-        edge_lon.resize(edge_count);
-        edge_lat.resize(edge_count);
-        edge_nodes.resize(2, edge_count);
-        face_edges.resize(nvertex, nbFaces_);
-
-        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
-        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
-        size_t iIdx = 0;
-
-        for (int nf = 0; nf < nbFaces_; ++nf)
-        {
-          for (int nv1 = 0; nv1 < nvertex; ++nv1)
-          {
-            // Getting global indexes of edge's nodes
-            int nh1 = 0;
-            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
-            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
-            while (it1 == nodeHash2IdxGlo.end())
-            {
-              ++nh1;
-              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
-            }
-            int nh2 = 0;
-            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
-            while (it2 == nodeHash2IdxGlo.end())
-            {
-              ++nh2;
-              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
-            }
             // Getting edge global index
             size_t nodeIdxGlo1 = it1->second[0];
             size_t nodeIdxGlo2 = it2->second[0];
-            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
+            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
             if (nodeIdxGlo1 != nodeIdxGlo2)
             {
-              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
-              int edgeIdxGlo = (itIdx->second)[1];
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+              size_t ownerIdx = (itIdx->second)[0];
+              int edgeIdxGlo = 0;
               size_t faceIdxGlo = nbFacesAccum + nf;
 
-              if (mpiRank == (itIdx->second)[0])
+              if (myIdx == ownerIdx)
               {
+                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+                edgeIdxGlo = (itIdxGlo->second)[0];
                 double edgeLon;
                 double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
@@ -1212 +1096 @@
                 edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
               }
+              else
+              {
+                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+                edgeIdxGlo = (itIdxGlo->second)[0];
+              }
               face_edges(nv1,nf) = edgeIdxGlo;
               if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
@@ -1228 +1117 @@
         edgeIdxGloList.resizeAndPreserve(iIdx);
 
-        // (2.5) Saving remaining variables edge_faces and face_faces
+        // (2.4) Saving remaining variables edge_faces and face_faces
         edge_faces.resize(2, edge_count);
         face_faces.resize(nvertex, nbFaces_);
@@ -1236 +1125 @@
         CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
         CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
-        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx;
 
         for (int nf = 0; nf < nbFaces_; ++nf)
@@ -1261 +1149 @@
             size_t nodeIdxGlo2 = it2->second[0];
 
-            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
-            itIdx = edgeIdx2IdxGlo.find(myIdx);
+            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
             size_t faceIdxGlo = nbFacesAccum + nf;
-            int edgeIdxGlo = (itIdx->second)[1];
-
-            if (mpiRank == (itIdx->second)[0])
-            {
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+              int edgeIdxGlo = (itIdxGlo->second)[0];
               it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
               int face1 = it1->second[0];
@@ -1286 +1176 @@
             else
             {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+              size_t edgeIdxGlo = (itIdxGlo->second)[0];
               it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
               int face1 = it1->second[0];
@@ -1308 +1200 @@
           {
             hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
-//            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
-            size_t nodeIndex = generateNodeIndex(hashValues);
+            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
             for (int nh = 0; nh < 4; ++nh)
             {
@@ -1325 +1216 @@
 
         // (3.2) Generating global node indexes
-        // The ownership criterion: priority of the process with smaller rank.
-        // With any other criterion it is not possible to have consistent node indexing for different number of procs.
+        // The ownership criterion: priority of the process holding the smaller index
         // Maps generated in this step are:
-        // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
-        // nodeIdx2Idx = <idx, <rankOwner, idx>>
+        // nodeHash2Info = <hash, [[idxMin, rankMin], [idx1, rank1], [idx3, rank3]..]>
+        // nodeIdx2IdxMin = = <idx, idxMin>
+        // nodeIdx2IdxGlo = <idxMin, idxGlo>
 
         CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
@@ -1335 +1226 @@
         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
 
-        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
-        CArray<size_t,1> nodeIdxList(nbFaces_*nvertex*4);
-        size_t nIdx = 0;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
+        CArray<size_t,1> nodeIdxMinList(nbFaces_*nvertex*4);
+        size_t iIdxMin = 0;
 
         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
         {
-          size_t rankMin = (it->second)[1];
+          size_t idxMin = (it->second)[0];
           size_t idx = (it->second)[0];
           for (int i = 2; i < (it->second).size();)
           {
-            if ( (it->second)[i+1] < rankMin)
+            if (mpiRank == (it->second)[i+1])
             {
               idx = (it->second)[i];
-              rankMin = (it->second)[i+1];
-              (it->second)[i+1] = (it->second)[i-1];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
             }
             i += 2;
           }
-          if (nodeIdx2Idx.count(idx) == 0)
-          {
-            if (mpiRank == rankMin)
-            {
-              nodeIdx2Idx[idx].push_back(rankMin);
-              nodeIdx2Idx[idx].push_back(idx);
-            }
-            nodeIdxList(nIdx) = idx;
-            ++nIdx;
-          }
-        }
-
-//        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2Idx, comm);
-        // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
-        // Solution: global node indexing by hand.
-        // Maps modified in this step:
-        // nodeIdx2Idx = <idx, idxGlo>
-        int nodeCount = nodeIdx2Idx.size();
-        int nodeStart, nbNodes;
-        MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
-        int nNodes = nodeStart;
-        MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
-        nbNodesGlo = nNodes;
-
-        nodeStart -= nodeCount;
-        node_start = nodeStart;
-        node_count = nodeCount;
-        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
-        size_t count = 0;
-
-        for (int nf = 0; nf < nbFaces_; ++nf)
-        {
-          for (int nv = 0; nv < nvertex; ++nv)
-          {
-            vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
-            size_t nodeIdx = generateNodeIndex(hashValues);
-            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
-            if (it != nodeIdx2Idx.end())
-            {
-              if (dummyMap.count(nodeIdx) == 0)
-              {
-                dummyMap[nodeIdx].push_back(nodeIdx);
-                (it->second)[1] = node_start + count;
-                ++count;
-              }
-            }
-          }
-        }
-        nodeIdxList.resizeAndPreserve(nIdx);
-        CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
-        dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
-        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
+          (it->second)[0] = idxMin;
+          if (nodeIdx2IdxMin.count(idx) == 0)
+          {
+            nodeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
+            nodeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+
+        nodeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
+        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
 
         // (3.3) Saving node data: node_lon, node_lat, and face_nodes
         // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
-//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
-//        node_count = dhtNodeIdxGlo.getIndexCount();
-//        node_start = dhtNodeIdxGlo.getIndexStart();
+        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
+        node_count = dhtNodeIdxGlo.getIndexCount();
+        node_start = dhtNodeIdxGlo.getIndexStart();
         node_lon.resize(node_count);
         node_lat.resize(node_count);
@@ -1427 +1286 @@
             vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
             vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
-            size_t nodeIdx1 = generateNodeIndex(hashValues1);
-            size_t nodeIdx2 = generateNodeIndex(hashValues2);
-            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
-            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
-            size_t ownerRank = (itNodeIdx1->second)[0];
-            nodeIdxGlo1 = (itNodeIdx1->second)[1];
-            nodeIdxGlo2 = (itNodeIdx2->second)[1];
-
-            if (mpiRank == ownerRank)
-            {
+            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
+            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxMin.find(myNodeIdx1);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxMin.find(myNodeIdx2);
+            size_t ownerNodeIdx = (itNodeIdx1->second)[0];
+
+            if (myNodeIdx1 == ownerNodeIdx)
+            {
+              itNodeIdxGlo1 = nodeIdx2IdxGlo.find(myNodeIdx1);
+              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
               node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
               node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
             }
+            else
+            {
+              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find(ownerNodeIdx);
+              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+            }
+            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]);
+            nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
             if (nodeIdxGlo1 != nodeIdxGlo2)
             {
               size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
-              edgeHash2Idx[edgeHash].push_back(edgeHash);
+              edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank));
               edgeHash2Idx[edgeHash].push_back(mpiRank);
               edgeHashList(nEdgeHash) = edgeHash;
@@ -1455 +1321 @@
         // (3.4) Generating global edge indexes
         // Maps generated in this step are:
-        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
-        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
+        // edgeIdx2IdxMin = = <idx, idxMin>
+        // edgeIdx2IdxGlo = <idxMin, idxGlo>
 
         CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
@@ -1463 +1329 @@
         // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
 
-        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
+        CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex);
+        iIdxMin = 0;
 
         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
         {
-          size_t rankMin = (it->second)[1];
+          size_t idxMin = (it->second)[0];
           size_t idx = (it->second)[0];
 
           for (int i = 2; i < (it->second).size();)
           {
-            if ((it->second)[i+1] < rankMin)
-            {
-              rankMin = (it->second)[i+1];
+            if (mpiRank == (it->second)[i+1])
+            {
               idx = (it->second)[i];
-              (it->second)[i+1] = (it->second)[i-1];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
             }
             i += 2;
           }
-          if (edgeIdx2Idx.count(idx) == 0)
-          {
-            if (mpiRank == rankMin)
-            {
-              edgeIdx2Idx[idx].push_back(rankMin);
-              edgeIdx2Idx[idx].push_back(idx);
-            }
-          }
-        }
-
-        int edgeCount = edgeIdx2Idx.size();
-        int edgeStart, nbEdges;
-        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
-        int nEdges = edgeStart;
-        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
-        nbEdgesGlo = nEdges;
-
-        edgeStart -= edgeCount;
-        edge_start = edgeStart;
-        edge_count = edgeCount;
-        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
-        count = 0;
-
-        for (int nf = 0; nf < nbFaces_; ++nf)
-        {
-          for (int nv1 = 0; nv1 < nvertex; ++nv1)
-          {
-            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
-            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
-            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
-            size_t nodeIdx1 = generateNodeIndex(hashValues1);
-            size_t nodeIdx2 = generateNodeIndex(hashValues2);
-            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
-            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
-            nodeIdxGlo1 = (itNodeIdx1->second)[1];
-            nodeIdxGlo2 = (itNodeIdx2->second)[1];
-
-            if (nodeIdxGlo1 != nodeIdxGlo2)
-            {
-              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
-              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
-              if (it != edgeIdx2Idx.end())
-              {
-                if (dummyEdgeMap.count(edgeIdx) == 0)
-                {
-                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
-                  (it->second)[1] = edge_start + count;
-                  ++count;
-                }
-              }
-            }
-          }
-        }
-        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
-        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
-        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
+          (it->second)[0] = idxMin;
+          if (edgeIdx2IdxMin.count(idx) == 0)
+          {
+            edgeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
+            edgeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+        edgeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
+        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
 
         // (3.5) Saving variables: edge_lon, edge_lat, face_edges
         // Creating map edgeIdxGlo2Face <idxGlo, face>
-//        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
-//        edge_count = dhtEdgeIdxGlo.getIndexCount();
-//        edge_start = dhtEdgeIdxGlo.getIndexStart();
-
+
+        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
+        edge_count = dhtEdgeIdxGlo.getIndexCount();
+        edge_start = dhtEdgeIdxGlo.getIndexStart();
         edge_lon.resize(edge_count);
         edge_lat.resize(edge_count);
@@ -1562 +1393 @@
             vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
 
-            size_t nodeIdx1 = generateNodeIndex(hashValues1);
-            size_t nodeIdx2 = generateNodeIndex(hashValues2);
-            it1 = nodeIdx2IdxGlo.find(nodeIdx1);
-            it2 = nodeIdx2IdxGlo.find(nodeIdx2);
-            size_t nodeIdxGlo1 = (it1->second)[1];
-            size_t nodeIdxGlo2 = (it2->second)[1];
+            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
+            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
+            it1 = nodeIdx2IdxMin.find(myNodeIdx1);
+            it2 = nodeIdx2IdxMin.find(myNodeIdx2);
+            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
+            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
+            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
 
             if (nodeIdxGlo1 != nodeIdxGlo2)
             {
-              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
-              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
-              int edgeIdxGlo = (itIdx->second)[1];
+              size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+              size_t ownerIdx = (itIdx->second)[0];
+              int edgeIdxGlo;
               size_t faceIdxGlo = nbFacesAccum + nf;
 
-              if (mpiRank == (itIdx->second)[0])
+              if (myIdx == ownerIdx)
               {
+                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+                edgeIdxGlo = (itIdxGlo->second)[0];
                 double edgeLon;
                 double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
@@ -1591 +1427 @@
                 edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
               }
+              else
+              {
+                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+                edgeIdxGlo = (itIdxGlo->second)[0];
+              }
               face_edges(nv1,nf) = edgeIdxGlo;
               if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
@@ -1624 +1465 @@
             vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
 
-            size_t myNodeIdx1 = generateNodeIndex(hashValues1);
-            size_t myNodeIdx2 = generateNodeIndex(hashValues2);
+            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
+            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
             if (myNodeIdx1 != myNodeIdx2)
             {
-              it1 = nodeIdx2IdxGlo.find(myNodeIdx1);
-              it2 = nodeIdx2IdxGlo.find(myNodeIdx2);
-              size_t nodeIdxGlo1 = (it1->second)[1];
-              size_t nodeIdxGlo2 = (it2->second)[1];
-              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
-              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
-              int edgeIdxGlo = (itIdx->second)[1];
-
+              it1 = nodeIdx2IdxMin.find(myNodeIdx1);
+              it2 = nodeIdx2IdxMin.find(myNodeIdx2);
+              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
+              itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
+              size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+              size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
+
+              size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+              size_t ownerIdx = (itIdx->second)[0];
               size_t faceIdxGlo = nbFacesAccum + nf;
 
-              if (mpiRank == (itIdx->second)[0])
+              if (myIdx == ownerIdx)
               {
+                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+                int edgeIdxGlo = (itIdxGlo->second)[0];
                 it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                 int face1 = it1->second[0];
@@ -1655 +1500 @@
                   face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
                 }
-              }
+              } // myIdx == ownerIdx
               else
               {
+                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+                size_t edgeIdxGlo = (itIdxGlo->second)[0];
                 it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                 int face1 = it1->second[0];
                 int face2 = it1->second[1];
                 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
-              }
+              } // myIdx != ownerIdx
             } // myNodeIdx1 != myNodeIdx2
             else
@@ -1668 +1515 @@
           }
         }
-
       }
       facesAreWritten = true;