Index: /XIOS/trunk/src/distribution_client.cpp
===================================================================
--- /XIOS/trunk/src/distribution_client.cpp	(revision 860)
+++ /XIOS/trunk/src/distribution_client.cpp	(revision 861)
@@ -16,9 +16,11 @@
    , nLocal_(), nGlob_(), nBeginLocal_(), nBeginGlobal_(),nZoomBegin_(), nZoomEnd_()
    , dataNIndex_(), dataDims_(), dataBegin_(), dataIndex_(), domainMasks_(), axisMasks_()
-   , gridMask_(), localDomainIndex_(), localAxisIndex_(), indexMap_(), indexDomainData_(), indexAxisData_()
+   , gridMask_(), indexMap_()
    , isDataDistributed_(true), axisNum_(0), domainNum_(0)
    , localDataIndex_(), localMaskIndex_()
    , globalLocalDataSendToServerMap_()
    , infoIndex_(), isComputed_(false)
+   , elementLocalIndex_(), elementGlobalIndex_(), elementIndexData_()
+   , elementZoomMask_(), elementNLocal_(), elementNGlobal_()
 {
 }
@@ -29,9 +31,11 @@
    , nLocal_(), nGlob_(), nBeginLocal_(), nBeginGlobal_(),nZoomBegin_(), nZoomEnd_()
    , dataNIndex_(), dataDims_(), dataBegin_(), dataIndex_(), domainMasks_(), axisMasks_()
-   , gridMask_(), localDomainIndex_(), localAxisIndex_(), indexMap_(), indexDomainData_(), indexAxisData_()
+   , gridMask_(), indexMap_()
    , isDataDistributed_(true), axisNum_(0), domainNum_(0)
    , localDataIndex_(), localMaskIndex_()
    , globalLocalDataSendToServerMap_()
    , infoIndex_(), isComputed_(false)
+   , elementLocalIndex_(), elementGlobalIndex_(), elementIndexData_()
+   , elementZoomMask_(), elementNLocal_(), elementNGlobal_()
 {
   readDistributionInfo(grid);
@@ -154,9 +158,21 @@
   idx = 0;
 
+  elementLocalIndex_.resize(numElement_);
+  elementGlobalIndex_.resize(numElement_);
+  elementIndexData_.resize(numElement_);
+  elementZoomMask_.resize(numElement_);
+  elementNLocal_.resize(numElement_);
+  elementNGlobal_.resize(numElement_);
+  elementNLocal_[0] = 1;
+  elementNGlobal_[0] = 1;
+  size_t localSize = 1, globalSize = 1;
+
   isDataDistributed_ = false;
   // Update all the vectors above
-  while (idx < numElement_)
+  for (idx = 0; idx < numElement_; ++idx)
   {
     bool isDomain = axisDomainOrder(idx);
+    elementNLocal_[idx] = localSize;
+    elementNGlobal_[idx] = globalSize;
 
     // If this is a domain
@@ -192,4 +208,6 @@
       isDataDistributed_ |= domList[domIndex]->isDistributed();
 
+      localSize *= nLocal_.at(indexMap_[idx]+1)* nLocal_.at(indexMap_[idx]);
+      globalSize *= nGlob_.at(indexMap_[idx]+1)* nGlob_.at(indexMap_[idx]);
       ++domIndex;
     }
@@ -211,7 +229,9 @@
       isDataDistributed_ |= axisList[axisIndex]->isDistributed();
 
+      localSize *= nLocal_.at(indexMap_[idx]);
+      globalSize *= nGlob_.at(indexMap_[idx]);
+
       ++axisIndex;
     }
-    ++idx;
   }
 }
@@ -225,11 +245,4 @@
 void CDistributionClient::createLocalDomainDataIndex()
 {
-  int numDomain = 0;
-  for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
-    if (axisDomainOrder_(i)) ++numDomain;
-
-  localDomainIndex_.resize(numDomain*2);
-  indexDomainData_.resize(numDomain);
-
   int idxDomain = 0;
   for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
@@ -237,6 +250,7 @@
     if (axisDomainOrder_(i))
     {
-      int iIdx, jIdx = 0, count = 0;
-      indexDomainData_[idxDomain].resize(dataNIndex_[i], false);
+      elementIndexData_[i].resize(dataNIndex_[i]);
+      elementIndexData_[i] = false;
+      int iIdx, jIdx = 0, count = 0, localIndex;
       for (int j = 0; j < dataNIndex_[i]; ++j)
       {
@@ -249,7 +263,35 @@
            (domainMasks_[idxDomain](iIdx + jIdx*nLocal_[indexMap_[i]])))
         {
-          (localDomainIndex_[idxDomain]).push_back(iIdx);
-          (localDomainIndex_[idxDomain*2+1]).push_back(jIdx);
-          indexDomainData_[idxDomain][j] = true;
+          ++count;
+          elementIndexData_[i](j) = true;
+        }
+      }
+
+      elementLocalIndex_[i].resize(count);
+      elementGlobalIndex_[i].resize(count);
+      elementZoomMask_[i].resize(count);
+      elementZoomMask_[i] = false;
+      count = 0;
+      CArray<bool,1>& tmpIndexElementData = elementIndexData_[i];
+      CArray<bool,1>& tmpZoomMaskElement = elementZoomMask_[i];
+      CArray<int,1>& tmpLocalElementIndex = elementLocalIndex_[i];
+      CArray<size_t,1>& tmpGlobalElementIndex = elementGlobalIndex_[i];
+      for (int j = 0; j < dataNIndex_[i]; ++j)
+      {
+        if (tmpIndexElementData(j))
+        {
+          iIdx = getDomainIndex((*dataIndex_[indexMap_[i]])(j), (*dataIndex_[indexMap_[i]+1])(j),
+                                dataBegin_[indexMap_[i]], dataBegin_[indexMap_[i]+1],
+                                dataDims_[i], nLocal_[indexMap_[i]], jIdx);
+          localIndex = tmpLocalElementIndex(count) = iIdx + jIdx * nLocal_[indexMap_[i]];
+          tmpGlobalElementIndex(count) = (*infoIndex_[indexMap_[i]])(localIndex) + ((*infoIndex_[indexMap_[i]+1])(localIndex))*nGlob_[indexMap_[i]];
+          if ((((*infoIndex_[indexMap_[i]])(localIndex)) <= nZoomEnd_[indexMap_[i]])
+             && (nZoomBegin_[indexMap_[i]] <= ((*infoIndex_[indexMap_[i]])(localIndex)))
+             && (((*infoIndex_[indexMap_[i]+1])(localIndex)) <= nZoomEnd_[indexMap_[i]+1])
+             && (nZoomBegin_[indexMap_[i]+1] <= ((*infoIndex_[indexMap_[i]+1])(localIndex))))
+          {
+            tmpZoomMaskElement(count) = true;
+          }
+          ++count;
         }
       }
@@ -264,11 +306,4 @@
 void CDistributionClient::createLocalAxisDataIndex()
 {
-  int numAxis = 0;
-  for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
-    if (!axisDomainOrder_(i)) ++numAxis;
-
-  localAxisIndex_.resize(numAxis);
-  indexAxisData_.resize(numAxis);
-
   int idxAxis = 0;
   for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
@@ -276,6 +311,7 @@
     if (!axisDomainOrder_(i))
     {
-      int iIdx = 0;
-      indexAxisData_[idxAxis].resize(dataNIndex_[i], false);
+      elementIndexData_[i].resize(dataNIndex_[i]);
+      elementIndexData_[i] = false;
+      int iIdx = 0, count = 0, localIndex = 0;
       for (int j = 0; j < dataNIndex_[i]; ++j)
       {
@@ -284,12 +320,282 @@
            (iIdx < nLocal_[indexMap_[i]]) && (axisMasks_[idxAxis](iIdx)))
         {
-          localAxisIndex_[idxAxis].push_back(iIdx);
-          indexAxisData_[idxAxis][j] = true;
+          ++count;
+          elementIndexData_[i](j) = true;
         }
       }
+
+      elementLocalIndex_[i].resize(count);
+      elementGlobalIndex_[i].resize(count);
+      elementZoomMask_[i].resize(count);
+      elementZoomMask_[i] = false;
+      count = 0;
+      CArray<bool,1>& tmpIndexElementData = elementIndexData_[i];
+      CArray<bool,1>& tmpZoomMaskElement = elementZoomMask_[i];
+      CArray<int,1>& tmpLocalElementIndex = elementLocalIndex_[i];
+      CArray<size_t,1>& tmpGlobalElementIndex = elementGlobalIndex_[i];
+      for (int j = 0; j < dataNIndex_[i]; ++j)
+      {
+        if (tmpIndexElementData(j))
+        {
+          iIdx = tmpLocalElementIndex(count) = getAxisIndex((*dataIndex_[indexMap_[i]])(j), dataBegin_[indexMap_[i]], nLocal_[indexMap_[i]]);
+          tmpGlobalElementIndex(count) = (*infoIndex_[indexMap_[i]])(iIdx);
+          if ((((*infoIndex_[indexMap_[i]])(iIdx)) <= nZoomEnd_[indexMap_[i]])
+             && (nZoomBegin_[indexMap_[i]] <= ((*infoIndex_[indexMap_[i]])(iIdx))))
+          {
+            tmpZoomMaskElement(count) = true;
+          }
+          ++count;
+        }
+      }
       ++idxAxis;
     }
   }
 }
+
+/*!
+   Create global index on client
+   In order to do the mapping between client-server, each client creates its own
+global index of sending data. This global index is then used to calculate to which server
+the client needs to send it data as well as which part of data belongs to the server.
+So as to make clients and server coherent in order of index, global index is calculated by
+take into account of C-convention, the rightmost dimension varies faster.
+*/
+void CDistributionClient::createGlobalIndexSendToServer()
+{
+  if (isComputed_) return;
+  isComputed_ = true;
+  createLocalDomainDataIndex();
+  createLocalAxisDataIndex();
+
+  int idxDomain = 0, idxAxis = 0;
+  std::vector<int> eachElementSize(numElement_);
+
+  // Precompute size of the loop
+  for (int i = 0; i < numElement_; ++i)
+  {
+    eachElementSize[i] = elementLocalIndex_[i].numElements();
+  }
+
+  //   Compute size of the global index on client
+  std::vector<StdSize> idxLoop(numElement_,0);
+  std::vector<StdSize> currentIndex(numElement_,0);
+  std::vector<StdSize> currentGlobalIndex(numElement_,0);
+  int innerLoopSize = eachElementSize[0];
+  size_t idx = 0, indexLocalDataOnClientCount = 0, indexSend2ServerCount = 0;
+  size_t ssize = 1;
+  for (int i = 0; i < numElement_; ++i) ssize *= eachElementSize[i];
+  while (idx < ssize)
+  {
+    for (int i = 0; i < numElement_-1; ++i)
+    {
+      if (idxLoop[i] == eachElementSize[i])
+      {
+        idxLoop[i] = 0;
+        ++idxLoop[i+1];
+      }
+    }
+
+    // Find out outer index
+    // Depending the inner-most element is axis or domain,
+    // The outer loop index begins correspondingly at one (1) or zero (0)
+    for (int i = 1; i < numElement_; ++i)
+    {
+      currentIndex[i] = elementLocalIndex_[i](idxLoop[i]);
+    }
+
+    // Inner most index
+    for (int i = 0; i < innerLoopSize; ++i)
+    {
+      int gridMaskIndex = 0;
+      for (int k = 0; k < this->numElement_; ++k)
+      {
+        gridMaskIndex += (currentIndex[k])*elementNLocal_[k];
+      }
+
+      if (gridMask_(gridMaskIndex))
+      {
+        ++indexLocalDataOnClientCount;
+        bool isIndexOnServer = true;
+
+        for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
+        {
+          isIndexOnServer = isIndexOnServer && elementZoomMask_[idxElement](idxLoop[idxElement]);
+        }
+        if (isIndexOnServer) ++indexSend2ServerCount;
+      }
+    }
+    idxLoop[0] += innerLoopSize;
+    idx += innerLoopSize;
+  }
+
+  // Now allocate these arrays
+  localDataIndex_.resize(indexLocalDataOnClientCount);
+  localMaskIndex_.resize(indexSend2ServerCount);
+  globalLocalDataSendToServerMap_.rehash(std::ceil(indexSend2ServerCount/globalLocalDataSendToServerMap_.max_load_factor())); //globalLocalDataSendToServerMap_.reserve(indexSend2ServerCount);
+
+  // We need to loop with data index
+  idxLoop.assign(numElement_,0);
+  idx = indexLocalDataOnClientCount = indexSend2ServerCount = 0;
+  ssize = 1; for (int i = 0; i < numElement_; ++i) ssize *= dataNIndex_[i];
+  innerLoopSize = dataNIndex_[0];
+  int countLocalData = 0;
+  std::vector<int> correctOuterIndex(numElement_,0);
+  bool isOuterIndexCorrect = true;
+  while (idx < ssize)
+  {
+    for (int i = 0; i < numElement_-1; ++i)
+    {
+      if (idxLoop[i] == dataNIndex_[i])
+      {
+        idxLoop[i] = 0;
+        correctOuterIndex[i] = 0;
+        ++idxLoop[i+1];
+        if (isOuterIndexCorrect) ++correctOuterIndex[i+1];
+      }
+    }
+
+    // Depending the inner-most element axis or domain,
+    // The outer loop index begins correspondingly at one (1) or zero (0)
+    bool isIndexElementDataCorrect = true;
+    for (int i = 1; i < numElement_; ++i)
+    {
+      if (elementIndexData_[i](idxLoop[i]))
+      {
+        currentIndex[i] = elementLocalIndex_[i](correctOuterIndex[i]);
+        currentGlobalIndex[i] = elementGlobalIndex_[i](correctOuterIndex[i]);
+        isIndexElementDataCorrect &= true;
+      }
+      else isIndexElementDataCorrect = false;
+    }
+
+    isOuterIndexCorrect = isIndexElementDataCorrect;
+
+    // Inner most index
+    int correctIndexElement = 0;
+    for (int i = 0; i < innerLoopSize; ++i)
+    {
+      bool isCurrentIndexDataCorrect = isOuterIndexCorrect;
+      if (elementIndexData_[0](i))
+      {
+        currentIndex[0] = elementLocalIndex_[0](correctIndexElement);
+        currentGlobalIndex[0] = elementGlobalIndex_[0](correctIndexElement);
+        isCurrentIndexDataCorrect &= true;
+        ++correctIndexElement;
+      }
+      else isCurrentIndexDataCorrect = false;
+
+      if (isCurrentIndexDataCorrect)
+      {
+        int gridMaskIndex = 0; //currentIndex[0];
+        for (int k = 0; k < this->numElement_; ++k)
+        {
+          gridMaskIndex += (currentIndex[k])*elementNLocal_[k];
+        }
+
+        if (gridMask_(gridMaskIndex))
+        {
+          localDataIndex_[indexLocalDataOnClientCount] = countLocalData;
+          bool isIndexOnServer = true;
+          for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
+          {
+            isIndexOnServer = isIndexOnServer && elementZoomMask_[idxElement](idxLoop[idxElement]);
+          }
+
+          if (isIndexOnServer)
+          {
+            size_t globalIndex = 0; //currentGlobalIndex[0];
+            for (int k = 0; k < numElement_; ++k)
+            {
+              globalIndex += (currentGlobalIndex[k])*elementNGlobal_[k];
+            }
+            globalLocalDataSendToServerMap_[globalIndex] = indexLocalDataOnClientCount;
+            localMaskIndex_[indexSend2ServerCount] = gridMaskIndex;
+            ++indexSend2ServerCount;
+          }
+          ++indexLocalDataOnClientCount;
+        }
+      }
+      ++countLocalData;
+    }
+    idxLoop[0] += innerLoopSize;
+    idx += innerLoopSize;
+  }
+}
+
+///*!
+//  Create local index of domain(s).
+//  A domain can have data index which even contains the "ghost" points. Very often, these
+//data surround the true data. In order to send correct data to server,
+//a client need to know index of the true data.
+//*/
+//void CDistributionClient::createLocalDomainDataIndex()
+//{
+//  int numDomain = 0;
+//  for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
+//    if (axisDomainOrder_(i)) ++numDomain;
+//
+//  localDomainIndex_.resize(numDomain*2);
+//  indexDomainData_.resize(numDomain);
+//
+//  int idxDomain = 0;
+//  for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
+//  {
+//    if (axisDomainOrder_(i))
+//    {
+//      int iIdx, jIdx = 0, count = 0;
+//      indexDomainData_[idxDomain].resize(dataNIndex_[i], false);
+//      for (int j = 0; j < dataNIndex_[i]; ++j)
+//      {
+//        iIdx = getDomainIndex((*dataIndex_[indexMap_[i]])(j), (*dataIndex_[indexMap_[i]+1])(j),
+//                              dataBegin_[indexMap_[i]], dataBegin_[indexMap_[i]+1],
+//                              dataDims_[i], nLocal_[indexMap_[i]], jIdx);
+//
+//        if ((iIdx >= nBeginLocal_[indexMap_[i]]) && (iIdx < nLocal_[indexMap_[i]]) &&
+//           (jIdx >= nBeginLocal_[indexMap_[i]+1]) && (jIdx < nLocal_[indexMap_[i]+1]) &&
+//           (domainMasks_[idxDomain](iIdx + jIdx*nLocal_[indexMap_[i]])))
+//        {
+//          (localDomainIndex_[idxDomain]).push_back(iIdx);
+//          (localDomainIndex_[idxDomain*2+1]).push_back(jIdx);
+//          indexDomainData_[idxDomain][j] = true;
+//        }
+//      }
+//      ++idxDomain;
+//    }
+//  }
+//}
+
+///*!
+//  Create local index of axis.
+//*/
+//void CDistributionClient::createLocalAxisDataIndex()
+//{
+//  int numAxis = 0;
+//  for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
+//    if (!axisDomainOrder_(i)) ++numAxis;
+//
+//  localAxisIndex_.resize(numAxis);
+//  indexAxisData_.resize(numAxis);
+//
+//  int idxAxis = 0;
+//  for (int i = 0; i < axisDomainOrder_.numElements(); ++i)
+//  {
+//    if (!axisDomainOrder_(i))
+//    {
+//      int iIdx = 0;
+//      indexAxisData_[idxAxis].resize(dataNIndex_[i], false);
+//      for (int j = 0; j < dataNIndex_[i]; ++j)
+//      {
+//        iIdx = getAxisIndex((*dataIndex_[indexMap_[i]])(j), dataBegin_[indexMap_[i]], nLocal_[indexMap_[i]]);
+//        if ((iIdx >= nBeginLocal_[indexMap_[i]]) &&
+//           (iIdx < nLocal_[indexMap_[i]]) && (axisMasks_[idxAxis](iIdx)))
+//        {
+//          localAxisIndex_[idxAxis].push_back(iIdx);
+//          indexAxisData_[idxAxis][j] = true;
+//        }
+//      }
+//      ++idxAxis;
+//    }
+//  }
+//}
 
 void CDistributionClient::createGlobalIndex()
@@ -347,287 +653,287 @@
 
 
-/*!
-   Create global index on client
-   In order to do the mapping between client-server, each client creates its own
-global index of sending data. This global index is then used to calculate to which server
-the client needs to send it data as well as which part of data belongs to the server.
-So as to make clients and server coherent in order of index, global index is calculated by
-take into account of C-convention, the rightmost dimension varies faster.
-*/
-void CDistributionClient::createGlobalIndexSendToServer()
-{
-  if (isComputed_) return;
-  isComputed_ = true;
-  createLocalDomainDataIndex();
-  createLocalAxisDataIndex();
-
-  int idxDomain = 0, idxAxis = 0;
-  std::vector<int> eachElementSize(numElement_);
-
-  // Precompute size of the loop
-  for (int i = 0; i < numElement_; ++i)
-  {
-    if(axisDomainOrder_(i))
-    {
-      eachElementSize[i] = localDomainIndex_[idxDomain].size();
-      idxDomain += 2;
-    }
-    else
-    {
-      eachElementSize[i] = localAxisIndex_[idxAxis].size();
-      ++idxAxis;
-    }
-  }
-
-  //   Compute size of the global index on client
-  std::vector<StdSize> idxLoop(numElement_,0);
-  std::vector<StdSize> currentIndex(this->dims_,0);
-  int innerLoopSize = eachElementSize[0];
-  size_t idx = 0, indexLocalDataOnClientCount = 0, indexSend2ServerCount = 0;
-  size_t ssize = 1;
-  for (int i = 0; i < numElement_; ++i) ssize *= eachElementSize[i];
-  while (idx < ssize)
-  {
-    for (int i = 0; i < numElement_-1; ++i)
-    {
-      if (idxLoop[i] == eachElementSize[i])
-      {
-        idxLoop[i] = 0;
-        ++idxLoop[i+1];
-      }
-    }
-
-    // Find out outer index
-    // Depending the inner-most element is axis or domain,
-    // The outer loop index begins correspondingly at one (1) or zero (0)
-    idxDomain = idxAxis = 0;
-    if (axisDomainOrder_(0)) ++idxDomain;
-    else ++idxAxis;
-    for (int i = 1; i < numElement_; ++i)
-    {
-      if (axisDomainOrder_(i))
-      {
-        currentIndex[indexMap_[i]]   = localDomainIndex_[idxDomain][idxLoop[i]];
-        currentIndex[indexMap_[i]+1] = localDomainIndex_[idxDomain+1][idxLoop[i]];
-        idxDomain += 2;
-      }
-      else
-      {
-        currentIndex[indexMap_[i]]   = localAxisIndex_[idxAxis][idxLoop[i]];
-        ++idxAxis;
-      }
-    }
-
-    // Inner most index
-    idxDomain = idxAxis = 0;
-    for (int i = 0; i < innerLoopSize; ++i)
-    {
-      if (axisDomainOrder_(0))
-      {
-        currentIndex[0] = localDomainIndex_[idxDomain][i];
-        currentIndex[1] = localDomainIndex_[idxDomain+1][i];
-      }
-      else currentIndex[0]   = localAxisIndex_[idxAxis][i];
-
-      StdSize gridMaskIndex = currentIndex[0];
-      int mulDimMask = 1;
-      for (int k = 1; k < this->dims_; ++k)
-      {
-        mulDimMask *= nLocal_[k-1];
-        gridMaskIndex += (currentIndex[k])*mulDimMask;
-      }
-
-      if (gridMask_(gridMaskIndex))
-      {
-        ++indexLocalDataOnClientCount;
-        bool isIndexOnServer = true;
-
-        for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
-        {
-          int actualIdx = 0;
-          if (axisDomainOrder_(idxElement))
-          {
-            actualIdx = currentIndex[indexMap_[idxElement]]+currentIndex[indexMap_[idxElement]+1]*nLocal_[indexMap_[idxElement]];
-            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]])
-                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(actualIdx)))
-                                              && (((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]+1])
-                                              && (nZoomBegin_[indexMap_[idxElement]+1] <= ((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)));
-          }
-          else
-          {
-            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])) <= nZoomEnd_[indexMap_[idxElement]])
-                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])));
-          }
-        }
-        if (isIndexOnServer) ++indexSend2ServerCount;
-      }
-
-    }
-    idxLoop[0] += innerLoopSize;
-    idx += innerLoopSize;
-  }
-
-  // Now allocate these arrays
-  localDataIndex_.resize(indexLocalDataOnClientCount);
-  localMaskIndex_.resize(indexSend2ServerCount);
-  globalLocalDataSendToServerMap_.rehash(std::ceil(indexSend2ServerCount/globalLocalDataSendToServerMap_.max_load_factor())); //globalLocalDataSendToServerMap_.reserve(indexSend2ServerCount);
-
-  // We need to loop with data index
-  idxLoop.assign(numElement_,0);
-  idx = indexLocalDataOnClientCount = indexSend2ServerCount = 0;
-  ssize = 1; for (int i = 0; i < numElement_; ++i) ssize *= dataNIndex_[i];
-  innerLoopSize = dataNIndex_[0];
-  int countLocalData = 0;
-  std::vector<int> correctOuterIndex(numElement_,0);
-  bool isOuterIndexCorrect = true;
-  while (idx < ssize)
-  {
-    for (int i = 0; i < numElement_-1; ++i)
-    {
-      if (idxLoop[i] == dataNIndex_[i])
-      {
-        idxLoop[i] = 0;
-        correctOuterIndex[i] = 0;
-        ++idxLoop[i+1];
-        if (isOuterIndexCorrect) ++correctOuterIndex[i+1];
-      }
-    }
-
-    // Depending the inner-most element axis or domain,
-    // The outer loop index begins correspondingly at one (1) or zero (0)
-    idxDomain = idxAxis = 0;
-    if (axisDomainOrder_(0)) ++idxDomain;
-    else ++idxAxis;
-    bool isIndexDomainDataCorrect = true;
-    bool isIndexAxisDataCorrect = true;
-
-    for (int i = 1; i < numElement_; ++i)
-    {
-      if (axisDomainOrder_(i))
-      {
-        if (indexDomainData_[idxDomain][idxLoop[i]])
-        {
-          currentIndex[indexMap_[i]]   = localDomainIndex_[idxDomain][correctOuterIndex[i]];
-          currentIndex[indexMap_[i]+1] = localDomainIndex_[idxDomain*2+1][correctOuterIndex[i]];
-          isIndexDomainDataCorrect &= true;
-        }
-        else isIndexDomainDataCorrect = false;
-        ++idxDomain;
-      }
-      else
-      {
-        if (indexAxisData_[idxAxis][idxLoop[i]])
-        {
-          currentIndex[indexMap_[i]]   = localAxisIndex_[idxAxis][correctOuterIndex[i]];
-          isIndexAxisDataCorrect &= true;
-        }
-        else isIndexAxisDataCorrect = false;
-        ++idxAxis;
-      }
-    }
-
-    isOuterIndexCorrect = (isIndexAxisDataCorrect) && (isIndexDomainDataCorrect);
-
-    // Inner most index
-    idxDomain = idxAxis = 0;
-    int correctIndexDomain = 0, correctIndexAxis = 0;
-    for (int i = 0; i < innerLoopSize; ++i)
-    {
-      bool isCurrentIndexDomainDataCorrect = isIndexDomainDataCorrect;
-      bool isCurrentIndexAxisDataCorrect = isIndexAxisDataCorrect;
-
-      if (axisDomainOrder_(0))
-      {
-        if (indexDomainData_[idxDomain][i])
-        {
-          currentIndex[0] = localDomainIndex_[idxDomain][correctIndexDomain];
-          currentIndex[1] = localDomainIndex_[idxDomain+1][correctIndexDomain];
-          isCurrentIndexDomainDataCorrect &= true;
-          ++correctIndexDomain;
-        }
-        else isCurrentIndexDomainDataCorrect = false;
-      }
-      else
-      {
-        if (indexAxisData_[idxAxis][i])
-        {
-          currentIndex[0] = localAxisIndex_[idxAxis][correctIndexAxis];
-          isCurrentIndexAxisDataCorrect &= true;
-          ++correctIndexAxis;
-        }
-        else isCurrentIndexAxisDataCorrect = false;
-      }
-
-      int gridMaskIndex = currentIndex[0];
-      int mulDimMask = 1;
-      for (int k = 1; k < this->dims_; ++k)
-      {
-        mulDimMask *= nLocal_[k-1];
-        gridMaskIndex += (currentIndex[k])*mulDimMask;
-      }
-
-      if (isCurrentIndexDomainDataCorrect &&
-          isCurrentIndexAxisDataCorrect &&
-          gridMask_(gridMaskIndex))
-      {
-        localDataIndex_[indexLocalDataOnClientCount] = countLocalData;
-        bool isIndexOnServer = true;
-        for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
-        {
-          int actualIdx = 0;
-          if (axisDomainOrder_(idxElement))
-          {
-            actualIdx = currentIndex[indexMap_[idxElement]]+currentIndex[indexMap_[idxElement]+1]*nLocal_[indexMap_[idxElement]];
-            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]])
-                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(actualIdx)))
-                                              && (((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]+1])
-                                              && (nZoomBegin_[indexMap_[idxElement]+1] <= ((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)));
-          }
-          else
-          {
-            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])) <= nZoomEnd_[indexMap_[idxElement]])
-                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])));
-          }
-        }
-
-        if (isIndexOnServer)
-        {
-          int actualIdx = (axisDomainOrder_(0)) ? currentIndex[0]+currentIndex[1]*nLocal_[0]
-                                                : currentIndex[0];
-          size_t globalIndex = (*infoIndex_[0])(actualIdx); //idxLoop[0] + nBeginGlobal_[0];
-          size_t mulDim = 1;
-          for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
-          {
-            if (axisDomainOrder_(idxElement))
-            {
-              actualIdx = currentIndex[indexMap_[idxElement]]+currentIndex[indexMap_[idxElement]+1]*nLocal_[indexMap_[idxElement]];
-              int jb = (0 == idxElement) ? 1 : 0;
-              for (int j = jb; j <= 1; ++j)
-              {
-                mulDim *= nGlob_[indexMap_[idxElement]+j-1];
-                globalIndex += ((*infoIndex_[indexMap_[idxElement]+j])(actualIdx))*mulDim;
-              }
-            }
-            else
-            {
-              if (0 != idxElement)
-              {
-                mulDim *= nGlob_[indexMap_[idxElement]-1];
-                globalIndex += ((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]]))*mulDim;
-              }
-            }
-          }
-          globalLocalDataSendToServerMap_[globalIndex] = indexLocalDataOnClientCount;
-          localMaskIndex_[indexSend2ServerCount] = gridMaskIndex;
-          ++indexSend2ServerCount;
-        }
-        ++indexLocalDataOnClientCount;
-      }
-      ++countLocalData;
-    }
-    idxLoop[0] += innerLoopSize;
-    idx += innerLoopSize;
-  }
-}
+///*!
+//   Create global index on client
+//   In order to do the mapping between client-server, each client creates its own
+//global index of sending data. This global index is then used to calculate to which server
+//the client needs to send it data as well as which part of data belongs to the server.
+//So as to make clients and server coherent in order of index, global index is calculated by
+//take into account of C-convention, the rightmost dimension varies faster.
+//*/
+//void CDistributionClient::createGlobalIndexSendToServer()
+//{
+//  if (isComputed_) return;
+//  isComputed_ = true;
+//  createLocalDomainDataIndex();
+//  createLocalAxisDataIndex();
+//
+//  int idxDomain = 0, idxAxis = 0;
+//  std::vector<int> eachElementSize(numElement_);
+//
+//  // Precompute size of the loop
+//  for (int i = 0; i < numElement_; ++i)
+//  {
+//    if(axisDomainOrder_(i))
+//    {
+//      eachElementSize[i] = localDomainIndex_[idxDomain].size();
+//      idxDomain += 2;
+//    }
+//    else
+//    {
+//      eachElementSize[i] = localAxisIndex_[idxAxis].size();
+//      ++idxAxis;
+//    }
+//  }
+//
+//  //   Compute size of the global index on client
+//  std::vector<StdSize> idxLoop(numElement_,0);
+//  std::vector<StdSize> currentIndex(this->dims_,0);
+//  int innerLoopSize = eachElementSize[0];
+//  size_t idx = 0, indexLocalDataOnClientCount = 0, indexSend2ServerCount = 0;
+//  size_t ssize = 1;
+//  for (int i = 0; i < numElement_; ++i) ssize *= eachElementSize[i];
+//  while (idx < ssize)
+//  {
+//    for (int i = 0; i < numElement_-1; ++i)
+//    {
+//      if (idxLoop[i] == eachElementSize[i])
+//      {
+//        idxLoop[i] = 0;
+//        ++idxLoop[i+1];
+//      }
+//    }
+//
+//    // Find out outer index
+//    // Depending the inner-most element is axis or domain,
+//    // The outer loop index begins correspondingly at one (1) or zero (0)
+//    idxDomain = idxAxis = 0;
+//    if (axisDomainOrder_(0)) ++idxDomain;
+//    else ++idxAxis;
+//    for (int i = 1; i < numElement_; ++i)
+//    {
+//      if (axisDomainOrder_(i))
+//      {
+//        currentIndex[indexMap_[i]]   = localDomainIndex_[idxDomain][idxLoop[i]];
+//        currentIndex[indexMap_[i]+1] = localDomainIndex_[idxDomain+1][idxLoop[i]];
+//        idxDomain += 2;
+//      }
+//      else
+//      {
+//        currentIndex[indexMap_[i]]   = localAxisIndex_[idxAxis][idxLoop[i]];
+//        ++idxAxis;
+//      }
+//    }
+//
+//    // Inner most index
+//    idxDomain = idxAxis = 0;
+//    for (int i = 0; i < innerLoopSize; ++i)
+//    {
+//      if (axisDomainOrder_(0))
+//      {
+//        currentIndex[0] = localDomainIndex_[idxDomain][i];
+//        currentIndex[1] = localDomainIndex_[idxDomain+1][i];
+//      }
+//      else currentIndex[0]   = localAxisIndex_[idxAxis][i];
+//
+//      StdSize gridMaskIndex = currentIndex[0];
+//      int mulDimMask = 1;
+//      for (int k = 1; k < this->dims_; ++k)
+//      {
+//        mulDimMask *= nLocal_[k-1];
+//        gridMaskIndex += (currentIndex[k])*mulDimMask;
+//      }
+//
+//      if (gridMask_(gridMaskIndex))
+//      {
+//        ++indexLocalDataOnClientCount;
+//        bool isIndexOnServer = true;
+//
+//        for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
+//        {
+//          int actualIdx = 0;
+//          if (axisDomainOrder_(idxElement))
+//          {
+//            actualIdx = currentIndex[indexMap_[idxElement]]+currentIndex[indexMap_[idxElement]+1]*nLocal_[indexMap_[idxElement]];
+//            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]])
+//                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(actualIdx)))
+//                                              && (((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]+1])
+//                                              && (nZoomBegin_[indexMap_[idxElement]+1] <= ((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)));
+//          }
+//          else
+//          {
+//            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])) <= nZoomEnd_[indexMap_[idxElement]])
+//                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])));
+//          }
+//        }
+//        if (isIndexOnServer) ++indexSend2ServerCount;
+//      }
+//
+//    }
+//    idxLoop[0] += innerLoopSize;
+//    idx += innerLoopSize;
+//  }
+//
+//  // Now allocate these arrays
+//  localDataIndex_.resize(indexLocalDataOnClientCount);
+//  localMaskIndex_.resize(indexSend2ServerCount);
+//  globalLocalDataSendToServerMap_.rehash(std::ceil(indexSend2ServerCount/globalLocalDataSendToServerMap_.max_load_factor())); //globalLocalDataSendToServerMap_.reserve(indexSend2ServerCount);
+//
+//  // We need to loop with data index
+//  idxLoop.assign(numElement_,0);
+//  idx = indexLocalDataOnClientCount = indexSend2ServerCount = 0;
+//  ssize = 1; for (int i = 0; i < numElement_; ++i) ssize *= dataNIndex_[i];
+//  innerLoopSize = dataNIndex_[0];
+//  int countLocalData = 0;
+//  std::vector<int> correctOuterIndex(numElement_,0);
+//  bool isOuterIndexCorrect = true;
+//  while (idx < ssize)
+//  {
+//    for (int i = 0; i < numElement_-1; ++i)
+//    {
+//      if (idxLoop[i] == dataNIndex_[i])
+//      {
+//        idxLoop[i] = 0;
+//        correctOuterIndex[i] = 0;
+//        ++idxLoop[i+1];
+//        if (isOuterIndexCorrect) ++correctOuterIndex[i+1];
+//      }
+//    }
+//
+//    // Depending the inner-most element axis or domain,
+//    // The outer loop index begins correspondingly at one (1) or zero (0)
+//    idxDomain = idxAxis = 0;
+//    if (axisDomainOrder_(0)) ++idxDomain;
+//    else ++idxAxis;
+//    bool isIndexDomainDataCorrect = true;
+//    bool isIndexAxisDataCorrect = true;
+//
+//    for (int i = 1; i < numElement_; ++i)
+//    {
+//      if (axisDomainOrder_(i))
+//      {
+//        if (indexDomainData_[idxDomain][idxLoop[i]])
+//        {
+//          currentIndex[indexMap_[i]]   = localDomainIndex_[idxDomain][correctOuterIndex[i]];
+//          currentIndex[indexMap_[i]+1] = localDomainIndex_[idxDomain*2+1][correctOuterIndex[i]];
+//          isIndexDomainDataCorrect &= true;
+//        }
+//        else isIndexDomainDataCorrect = false;
+//        ++idxDomain;
+//      }
+//      else
+//      {
+//        if (indexAxisData_[idxAxis][idxLoop[i]])
+//        {
+//          currentIndex[indexMap_[i]]   = localAxisIndex_[idxAxis][correctOuterIndex[i]];
+//          isIndexAxisDataCorrect &= true;
+//        }
+//        else isIndexAxisDataCorrect = false;
+//        ++idxAxis;
+//      }
+//    }
+//
+//    isOuterIndexCorrect = (isIndexAxisDataCorrect) && (isIndexDomainDataCorrect);
+//
+//    // Inner most index
+//    idxDomain = idxAxis = 0;
+//    int correctIndexDomain = 0, correctIndexAxis = 0;
+//    for (int i = 0; i < innerLoopSize; ++i)
+//    {
+//      bool isCurrentIndexDomainDataCorrect = isIndexDomainDataCorrect;
+//      bool isCurrentIndexAxisDataCorrect = isIndexAxisDataCorrect;
+//
+//      if (axisDomainOrder_(0))
+//      {
+//        if (indexDomainData_[idxDomain][i])
+//        {
+//          currentIndex[0] = localDomainIndex_[idxDomain][correctIndexDomain];
+//          currentIndex[1] = localDomainIndex_[idxDomain+1][correctIndexDomain];
+//          isCurrentIndexDomainDataCorrect &= true;
+//          ++correctIndexDomain;
+//        }
+//        else isCurrentIndexDomainDataCorrect = false;
+//      }
+//      else
+//      {
+//        if (indexAxisData_[idxAxis][i])
+//        {
+//          currentIndex[0] = localAxisIndex_[idxAxis][correctIndexAxis];
+//          isCurrentIndexAxisDataCorrect &= true;
+//          ++correctIndexAxis;
+//        }
+//        else isCurrentIndexAxisDataCorrect = false;
+//      }
+//
+//      int gridMaskIndex = currentIndex[0];
+//      int mulDimMask = 1;
+//      for (int k = 1; k < this->dims_; ++k)
+//      {
+//        mulDimMask *= nLocal_[k-1];
+//        gridMaskIndex += (currentIndex[k])*mulDimMask;
+//      }
+//
+//      if (isCurrentIndexDomainDataCorrect &&
+//          isCurrentIndexAxisDataCorrect &&
+//          gridMask_(gridMaskIndex))
+//      {
+//        localDataIndex_[indexLocalDataOnClientCount] = countLocalData;
+//        bool isIndexOnServer = true;
+//        for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
+//        {
+//          int actualIdx = 0;
+//          if (axisDomainOrder_(idxElement))
+//          {
+//            actualIdx = currentIndex[indexMap_[idxElement]]+currentIndex[indexMap_[idxElement]+1]*nLocal_[indexMap_[idxElement]];
+//            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]])
+//                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(actualIdx)))
+//                                              && (((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)) <= nZoomEnd_[indexMap_[idxElement]+1])
+//                                              && (nZoomBegin_[indexMap_[idxElement]+1] <= ((*infoIndex_[indexMap_[idxElement]+1])(actualIdx)));
+//          }
+//          else
+//          {
+//            isIndexOnServer = isIndexOnServer && (((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])) <= nZoomEnd_[indexMap_[idxElement]])
+//                                              && (nZoomBegin_[indexMap_[idxElement]] <= ((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]])));
+//          }
+//        }
+//
+//        if (isIndexOnServer)
+//        {
+//          int actualIdx = (axisDomainOrder_(0)) ? currentIndex[0]+currentIndex[1]*nLocal_[0]
+//                                                : currentIndex[0];
+//          size_t globalIndex = (*infoIndex_[0])(actualIdx); //idxLoop[0] + nBeginGlobal_[0];
+//          size_t mulDim = 1;
+//          for (int idxElement = 0; idxElement < this->numElement_; ++idxElement)
+//          {
+//            if (axisDomainOrder_(idxElement))
+//            {
+//              actualIdx = currentIndex[indexMap_[idxElement]]+currentIndex[indexMap_[idxElement]+1]*nLocal_[indexMap_[idxElement]];
+//              int jb = (0 == idxElement) ? 1 : 0;
+//              for (int j = jb; j <= 1; ++j)
+//              {
+//                mulDim *= nGlob_[indexMap_[idxElement]+j-1];
+//                globalIndex += ((*infoIndex_[indexMap_[idxElement]+j])(actualIdx))*mulDim;
+//              }
+//            }
+//            else
+//            {
+//              if (0 != idxElement)
+//              {
+//                mulDim *= nGlob_[indexMap_[idxElement]-1];
+//                globalIndex += ((*infoIndex_[indexMap_[idxElement]])(currentIndex[indexMap_[idxElement]]))*mulDim;
+//              }
+//            }
+//          }
+//          globalLocalDataSendToServerMap_[globalIndex] = indexLocalDataOnClientCount;
+//          localMaskIndex_[indexSend2ServerCount] = gridMaskIndex;
+//          ++indexSend2ServerCount;
+//        }
+//        ++indexLocalDataOnClientCount;
+//      }
+//      ++countLocalData;
+//    }
+//    idxLoop[0] += innerLoopSize;
+//    idx += innerLoopSize;
+//  }
+//}
 
 /*!
@@ -675,5 +981,5 @@
   Return global local data mapping of client
 */
-const CDistributionClient::GlobalLocalDataMap& CDistributionClient::getGlobalLocalDataSendToServer()
+CDistributionClient::GlobalLocalDataMap& CDistributionClient::getGlobalLocalDataSendToServer()
 {
   if (!isComputed_) createGlobalIndexSendToServer();
Index: /XIOS/trunk/src/distribution_client.hpp
===================================================================
--- /XIOS/trunk/src/distribution_client.hpp	(revision 860)
+++ /XIOS/trunk/src/distribution_client.hpp	(revision 861)
@@ -43,5 +43,5 @@
 
     virtual const std::vector<int>& getLocalDataIndexOnClient();
-    const GlobalLocalDataMap& getGlobalLocalDataSendToServer();
+    GlobalLocalDataMap& getGlobalLocalDataSendToServer();
     const std::vector<int>& getLocalMaskIndexOnClient();
 
@@ -101,11 +101,11 @@
     std::vector<CArray<bool,1> > axisMasks_; //!< Axis mask
 
-    std::vector<std::vector<int> > localDomainIndex_;
-    std::vector<std::vector<int> > localAxisIndex_;
     std::vector<int> indexMap_; //!< Mapping element index to dimension index
-
-    // The correct index of a domain has true value, the ghost one has false value
-    std::vector<std::vector<bool> > indexDomainData_;
-    std::vector<std::vector<bool> > indexAxisData_;
+    std::vector<CArray<int,1> > elementLocalIndex_;  //!< Local index of each element
+    std::vector<CArray<size_t,1> > elementGlobalIndex_; //!< Global index of each element
+    std::vector<CArray<bool,1> > elementIndexData_; //!< // The correct index of a domain has true value, the ghost one has false value
+    std::vector<CArray<bool,1> > elementZoomMask_; //!< Only zoomed region are true
+    std::vector<size_t> elementNLocal_;
+    std::vector<size_t> elementNGlobal_;
 
     //!< (Only for grid with one axis or scalar)Flag to determine whether data is distributed or not
