Changeset 1568 for XIOS/dev/dev_olga/src/node

Timestamp:

08/28/18 16:38:30 (6 years ago)

Author:

oabramkina

Message:

Grid mask is now applied in the source filter of clients: unmasked values are replaced by NaN. It is not reconstructed any more by servers.

This needs to be tested more rigorously before commiting to trunk.

Location:

XIOS/dev/dev_olga/src/node

Files:

• axis.cpp (modified) (16 diffs)
• axis.hpp (modified) (1 diff)
• domain.cpp (modified) (2 diffs)
• field.cpp (modified) (6 diffs)
• grid.cpp (modified) (3 diffs)
• grid.hpp (modified) (5 diffs)

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

@@ -285 +285 @@
       }
 
-      // Remove this check because it doen't make sense in case of a hole or overlapping axes
       if (!this->value.isEmpty())
       {
-//        StdSize true_size = value.numElements();
-//        if (this->n.getValue() != true_size)
-//          ERROR("CAxis::checkAttributes(void)",
-//                << "[ id = '" << getId() << "' , context = '" << CObjectFactory::GetCurrentContextId() << "' ] "
-//                << "The axis is wrongly defined, attribute 'value' has a different size (" << true_size << ") than the one defined by the \'size\' attribute (" << n.getValue() << ").");
+        StdSize true_size = value.numElements();
+        if (this->n.getValue() != true_size)
+          ERROR("CAxis::checkAttributes(void)",
+                << "[ id = '" << getId() << "' , context = '" << CObjectFactory::GetCurrentContextId() << "' ] "
+                << "The axis is wrongly defined, attribute 'value' has a different size (" << true_size << ") than the one defined by the \'size\' attribute (" << n.getValue() << ").");
         this->hasValue = true;
       }
@@ -310 +309 @@
       Check the validity of data and fill in values if any.
    */
+
    void CAxis::checkData()
    {
@@ -328 +328 @@
       {
         data_index.resize(data_n);
-        for (int i = 0; i < data_n; ++i) data_index(i) = i;
-      }
+        for (int i = 0; i < data_n; ++i)
+        {
+          if ((i+data_begin) >= 0 && (i+data_begin<n))
+            data_index(i) = i+data_begin;
+          else
+            data_index(i) = -1;
+        }
+      }
+      else
+      {
+        if (data_index.numElements() != data_n)
+        {
+          ERROR("CAxis::checkData(void)",
+                << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
+                << "The size of data_index = "<< data_index.numElements() << "is not equal to the data size data_n = " << data_n.getValue() << ").");
+        }
+        data_index.resize(data_n);
+        for (int i = 0; i < data_n; ++i)
+        {
+          if ((i+data_begin) >= 0 && (i+data_begin<n))
+            data_index(i) = data_index(i);
+          else
+            data_index(i) = -1;
+        }
+      }
+
    }
 
@@ -344 +368 @@
       if (!mask.isEmpty())
       {
-         if (mask.extent(0) != n)
-           ERROR("CAxis::checkMask(void)",
-                 << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
-                 << "The mask does not have the same size as the local domain." << std::endl
-                 << "Local size is " << n.getValue() << "." << std::endl
-                 << "Mask size is " << mask.extent(0) << ".");
-      }
-      else // (mask.isEmpty())
-      { // If no mask was defined, we create a default one without any masked point.
-         mask.resize(n);
-         for (int i = 0; i < n; ++i)
-         {
-           mask(i) = true;
-         }
+        if (mask.extent(0) != n)
+        {
+          ERROR("CAxis::checkMask(void)",
+              << "[ id = " << this->getId() << " , context = '" << CObjectFactory::GetCurrentContextId() << " ] "
+              << "The mask does not have the same size as the local domain." << std::endl
+              << "Local size is " << n.getValue() << "." << std::endl
+              << "Mask size is " << mask.extent(0) << ".");
+        }
+      }
+      else
+      {
+        mask.resize(n);
+        mask = true;
       }
    }
@@ -547 +570 @@
 
         // Calculate the compressed index if any
-        std::set<int> writtenInd;
-        if (isCompressible_)
-        {
-          for (int idx = 0; idx < data_index.numElements(); ++idx)
-          {
-            int ind = CDistributionClient::getAxisIndex(data_index(idx), data_begin, ni);
-
-            if (ind >= 0 && ind < ni && mask(ind))
-            {
-              ind += ibegin;
-              writtenInd.insert(ind);
-            }
-          }
-        }
+//        std::set<int> writtenInd;
+//        if (isCompressible_)
+//        {
+//          for (int idx = 0; idx < data_index.numElements(); ++idx)
+//          {
+//            int ind = CDistributionClient::getAxisIndex(data_index(idx), data_begin, ni);
+//
+//            if (ind >= 0 && ind < ni && mask(ind))
+//            {
+//              ind += ibegin;
+//              writtenInd.insert(ind);
+//            }
+//          }
+//        }
 
         // Compute the global index of the current client (process) hold
@@ -725 +748 @@
         }                 
       }
+//
+//      nbWritten = 0;
+//      for (int idx = 0; idx < data_index.numElements(); ++idx)
+//      {
+//        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_index(idx)))
+//        {
+//          ++nbWritten;
+//        }
+//      }
+//
+//      compressedIndexToWriteOnServer[writtenCommSize].resize(nbWritten);
+//      nbWritten = 0;
+//      for (int idx = 0; idx < data_index.numElements(); ++idx)
+//      {
+//        if (localGlobalIndexMap.end() != localGlobalIndexMap.find(data_index(idx)))
+//        {
+//          compressedIndexToWriteOnServer[writtenCommSize](nbWritten) = localGlobalIndexMap[data_index(idx)];
+//          ++nbWritten;
+//        }
+//      }
 
       nbWritten = 0;
@@ -959 +1002 @@
 
   /*
-    Send attributes of axis from a group of client to other group of clients/servers 
-    on supposing that these attributes are distributed among the clients of the sending group
-    In the future, if new attributes are added, they should also be processed in this function
+    Send axis attributes from a group of clients to another group of clients/servers
+    supposing that these attributes are distributed among the clients of the sending group
+    In future, if new attributes are added, they should also be processed in this function
   */
   void CAxis::sendDistributedAttributes(void)
   {
-    int ns, n, i, j, ind, nv, idx;
+    int ind, idx;
     std::list<CContextClient*>::iterator it;
 
@@ -982 +1025 @@
       list<CArray<string,1> > list_label;
 
+      // Cut off the ghost points
       int nbIndex = index.numElements();
       CArray<int,1> dataIndex(nbIndex);
@@ -988 +1032 @@
       {
         if (0 <= data_index(idx) && data_index(idx) < nbIndex)
-          dataIndex(idx) = 1;
+          dataIndex(data_index(idx)) = 1;
       }
 
@@ -995 +1039 @@
       for (int k = 0; k < connectedServerRank_[nbServer].size(); ++k)
       {
-        int nbData = 0;
+        int nbData = 0, nbDataCount = 0;
         int rank = connectedServerRank_[nbServer][k];
         it = indSrv_[nbServer].find(rank);
@@ -1002 +1046 @@
 
         list_indi.push_back(CArray<int,1>(nbData));
-        list_dataInd.push_back(CArray<int,1>(nbData));        
+        list_dataInd.push_back(CArray<int,1>(nbData));
         list_mask.push_back(CArray<bool,1>(nbData));
 
@@ -1015 +1059 @@
 
         CArray<int,1>& indi = list_indi.back();
-        CArray<int,1>& dataIndi = list_dataInd.back();        
+        CArray<int,1>& dataIndi = list_dataInd.back();
+        dataIndi = -1;
         CArray<bool,1>& maskIndi = list_mask.back();
 
-        for (n = 0; n < nbData; ++n)
+        for (int n = 0; n < nbData; ++n)
         {
           idx = static_cast<int>(it->second[n]);
@@ -1162 +1207 @@
     nonCompressedData = -1;   
     mask.resize(nbData);
+    mask = true;
     if (hasValue)
       value.resize(nbData);
@@ -1199 +1245 @@
     }
     
-    int nbCompressedData = 0; 
+    int nbCompressedData = 0;
     for (idx = 0; idx < nonCompressedData.numElements(); ++idx)
     {
       if (0 <= nonCompressedData(idx))
-        ++nbCompressedData;        
+        ++nbCompressedData;
     }
 
@@ -1213 +1259 @@
       {
         data_index(nbCompressedData) = idx % n;
-        ++nbCompressedData;        
+        ++nbCompressedData;
       }
     }
 
     data_begin.setValue(0);
+    data_n.setValue(data_index.numElements());
   }
 
@@ -1330 +1377 @@
       clientsSet.insert(contextClient);
     }
-}
+  }
 
   void CAxis::parse(xml::CXMLNode & node)

XIOS/dev/dev_olga/src/node/axis.hpp

@@ -126 +126 @@
         bool hasLabel;
 
-        CArray<size_t,1> localIndexToWriteOnServer;        
+        CArray<size_t,1> localIndexToWriteOnServer;
 
       private:

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

@@ -2447 +2447 @@
          {
            iIndex = (index%ni_glo)-ibegin;
-           iIndex = (iIndex < 0) ? 0 : iIndex;
+           iIndex = (iIndex < 0) ? 0 : iIndex;    // ??
            jIndex = (index/ni_glo)-jbegin;
-           jIndex = (jIndex < 0) ? 0 : jIndex;
+           jIndex = (jIndex < 0) ? 0 : jIndex;    // ??
            nbIndLoc = iIndex + ni * jIndex;
            if (nbIndLoc < nbIndexGlobMax)
@@ -2896 +2896 @@
       buffer >> recvDataJIndex[i];
     }
-   
+
     int nbIndex = i_index.numElements();
     CArray<int,1> dataIIndex(nbIndex), dataJIndex(nbIndex);
     dataIIndex = -1; dataJIndex = -1;
-     
+
     nbIndex = 0;
     for (i = 0; i < nbReceived; ++i)

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

@@ -1027 +1027 @@
      {
         if (!instantDataFilter)
-          instantDataFilter = clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid,true));
+          instantDataFilter = clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, false));
 
 
@@ -1043 +1043 @@
      {
        if (!instantDataFilter)
-         instantDataFilter = clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true));
+         instantDataFilter = clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, false, false));
 
              // If the field data is to be read by the client or/and written to a file
@@ -1089 +1089 @@
          {
            checkTimeAttributes();
-           instantDataFilter = serverSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, freq_offset, true,
+           instantDataFilter = serverSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, false, freq_offset, true,
                                                                                                        detectMissingValues, defaultValue));
          }
@@ -1095 +1095 @@
          {
             if (check_if_active.isEmpty()) check_if_active = false; 
-            instantDataFilter = clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, false, NoneDu, false,
+            instantDataFilter = clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, false, true, NoneDu, false,
                                                                                                         detectMissingValues, defaultValue));
          }
@@ -1176 +1176 @@
          {
            checkTimeAttributes();
-           serverSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, freq_offset, true,
+           serverSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, false, freq_offset, true,
                                                                                    detectMissingValues, defaultValue));
          }
@@ -1193 +1193 @@
          {
            if (check_if_active.isEmpty()) check_if_active = false;
-           clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, NoneDu, false,
+           clientSourceFilter = std::shared_ptr<CSourceFilter>(new CSourceFilter(gc, grid, true, true, NoneDu, false,
                                                                                    detectMissingValues, defaultValue));
          }

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

@@ -698 +698 @@
      clientDistribution_ = new CDistributionClient(rank, this);
      // Get local data index on client
-     storeIndex_client.resize(clientDistribution_->getLocalDataIndexOnClient().size());
-     int nbStoreIndex = storeIndex_client.numElements();
-     for (int idx = 0; idx < nbStoreIndex; ++idx) storeIndex_client(idx) = (clientDistribution_->getLocalDataIndexOnClient())[idx];
+     int nbStoreIndex = clientDistribution_->getLocalDataIndexOnClient().size();
+     storeIndex_client.resize(nbStoreIndex);
+     storeMask_client.resize(nbStoreIndex);
+     for (int idx = 0; idx < nbStoreIndex; ++idx)
+     {
+       storeIndex_client(idx) = (clientDistribution_->getLocalDataIndexOnClient())[idx];
+       storeMask_client(idx) = (clientDistribution_->getLocalMaskIndexOnClient())[idx];
+     }
 
      if (0 == serverDistribution_) isDataDistributed_= clientDistribution_->isDataDistributed();
@@ -1293 +1298 @@
    }
 
+   void CGrid::maskField_arr(const double* const data, CArray<double, 1>& stored) const
+   {
+      const StdSize size = storeIndex_client.numElements();
+      stored.resize(size);
+
+      const double nanValue = std::numeric_limits<double>::quiet_NaN();
+      for(StdSize i = 0; i < size; i++) stored(i) = (storeMask_client(i)) ? data[storeIndex_client(i)] : nanValue;
+   }
+
    void CGrid::uncompressField_arr(const double* const data, CArray<double, 1>& out) const
    {
@@ -1749 +1763 @@
         }
 
-        modifyMaskSize(nSize, false);
-
-        // These below codes are reserved for future
-        CDistributionServer srvDist(server->intraCommRank, nBegin, nSize, nBeginGlobal, nGlob); 
-        map<int, CArray<size_t, 1> >::iterator itb = outGlobalIndexFromClient.begin(),
-                                               ite = outGlobalIndexFromClient.end(), it;  
-        const CDistributionServer::GlobalLocalMap&  globalLocalMask = srvDist.getGlobalLocalIndex();
-        CDistributionServer::GlobalLocalMap::const_iterator itSrv;
-        size_t nb = 0;
-        for (it = itb; it != ite; ++it)
-        {
-          CArray<size_t,1>& globalInd = it->second;
-          for (size_t idx = 0; idx < globalInd.numElements(); ++idx)
-          {
-            if (globalLocalMask.end() != globalLocalMask.find(globalInd(idx))) ++nb;
-          }
-        }
-        
-        CArray<int,1> indexToModify(nb);
-        nb = 0;    
-        for (it = itb; it != ite; ++it)
-        {
-          CArray<size_t,1>& globalInd = it->second;
-          for (size_t idx = 0; idx < globalInd.numElements(); ++idx)
-          {
-            itSrv = globalLocalMask.find(globalInd(idx));
-            if (globalLocalMask.end() != itSrv) 
-            {
-              indexToModify(nb) = itSrv->second;
-              ++nb;
-            }
-          }
-        }
-
-        modifyMask(indexToModify, true);
+//        modifyMaskSize(nSize, false);
+//        modifyMaskSize(nSize, true);   // Grid mask should be gone beyond client source filter
+//
+//        // These below codes are reserved for future
+//        CDistributionServer srvDist(server->intraCommRank, nBegin, nSize, nBeginGlobal, nGlob);
+//        map<int, CArray<size_t, 1> >::iterator itb = outGlobalIndexFromClient.begin(),
+//                                               ite = outGlobalIndexFromClient.end(), it;
+//        const CDistributionServer::GlobalLocalMap&  globalLocalMask = srvDist.getGlobalLocalIndex();
+//        CDistributionServer::GlobalLocalMap::const_iterator itSrv;
+//        size_t nb = 0;
+//        for (it = itb; it != ite; ++it)
+//        {
+//          CArray<size_t,1>& globalInd = it->second;
+//          for (size_t idx = 0; idx < globalInd.numElements(); ++idx)
+//          {
+//            if (globalLocalMask.end() != globalLocalMask.find(globalInd(idx))) ++nb;
+//          }
+//        }
+//
+//        CArray<int,1> indexToModify(nb);
+//        nb = 0;
+//        for (it = itb; it != ite; ++it)
+//        {
+//          CArray<size_t,1>& globalInd = it->second;
+//          for (size_t idx = 0; idx < globalInd.numElements(); ++idx)
+//          {
+//            itSrv = globalLocalMask.find(globalInd(idx));
+//            if (globalLocalMask.end() != itSrv)
+//            {
+//              indexToModify(nb) = itSrv->second;
+//              ++nb;
+//            }
+//          }
+//        }
+
+//        modifyMask(indexToModify, true);
       }
 

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

@@ -94 +94 @@
          template <int n>
          void inputField(const CArray<double,n>& field, CArray<double,1>& stored) const;
+         template <int n>
+         void maskField(const CArray<double,n>& field, CArray<double,1>& stored) const;
          template <int n>
          void outputField(const CArray<double,1>& stored, CArray<double,n>& field) const;  
@@ -208 +210 @@
       public:
          CArray<int, 1> storeIndex_client;
+         CArray<bool, 1> storeMask_client;
 
 /** Map containing indexes that will be sent in sendIndex(). */
@@ -273 +276 @@
         void restoreField_arr(const CArray<double, 1>& stored, double* const data) const;
         void uncompressField_arr(const double* const data, CArray<double, 1>& outData) const;
+        void maskField_arr(const double* const data, CArray<double, 1>& stored) const;
 
         void setVirtualDomainGroup(CDomainGroup* newVDomainGroup);
@@ -385 +389 @@
 
    template <int n>
+   void CGrid::maskField(const CArray<double,n>& field, CArray<double,1>& stored) const
+   {
+//#ifdef __XIOS_DEBUG
+      if (this->getDataSize() != field.numElements())
+         ERROR("void CGrid::inputField(const  CArray<double,n>& field, CArray<double,1>& stored) const",
+                << "[ Awaiting data of size = " << this->getDataSize() << ", "
+                << "Received data size = "      << field.numElements() << " ] "
+                << "The data array does not have the right size! "
+                << "Grid = " << this->getId())
+//#endif
+      this->maskField_arr(field.dataFirst(), stored);
+   }
+
+   template <int n>
    void CGrid::outputField(const CArray<double,1>& stored, CArray<double,n>& field) const
    {
@@ -417 +435 @@
                              bool createMask)
    {
-     if (!gridMask.isEmpty() || createMask)
-     {
-       int idx = 0;
-       int numElement = axisDomainOrder.numElements();
-       int dim = domainMasks.size() * 2 + axisMasks.size();
-       std::vector<CDomain*> domainP = this->getDomains();
-       std::vector<CAxis*> axisP = this->getAxis();
-
-       std::vector<int> idxLoop(dim,0), indexMap(numElement), eachDimSize(dim);
-       std::vector<int> currentIndex(dim);
-       int idxDomain = 0, idxAxis = 0;
+     int idx = 0;
+     int numElement = axisDomainOrder.numElements();
+     int dim = domainMasks.size() * 2 + axisMasks.size();
+     std::vector<CDomain*> domainP = this->getDomains();
+     std::vector<CAxis*> axisP = this->getAxis();
+
+     std::vector<int> idxLoop(dim,0), indexMap(numElement), eachDimSize(dim);
+     std::vector<int> currentIndex(dim);
+     int idxDomain = 0, idxAxis = 0;
+    for (int i = 0; i < numElement; ++i)
+    {
+      indexMap[i] = idx;
+      if (2 == axisDomainOrder(i)) {
+          eachDimSize[indexMap[i]]   = domainP[idxDomain]->ni;
+          eachDimSize[indexMap[i]+1] = domainP[idxDomain]->nj;
+          idx += 2; ++idxDomain;
+      }
+      else if (1 == axisDomainOrder(i)) {
+//        eachDimSize[indexMap[i]] = axisMasks[idxAxis]->numElements();
+        eachDimSize[indexMap[i]] = axisP[idxAxis]->n;
+        ++idx; ++idxAxis;
+      }
+      else {};
+    }
+
+    if (!gridMask.isEmpty() && !createMask)
+    {
+      for (int i = 0; i < dim; ++i)
+      {
+        if (gridMask.extent(i) != eachDimSize[i])
+          ERROR("CGrid::checkMask(void)",
+                << "The mask has one dimension whose size is different from the one of the local grid." << std::endl
+                << "Local size of dimension " << i << " is " << eachDimSize[i] << "." << std::endl
+                << "Mask size for dimension " << i << " is " << gridMask.extent(i) << "." << std::endl
+                << "Grid = " << this->getId())
+      }
+    }
+    else {
+        CArrayBoolTraits<CArray<bool,N> >::resizeArray(gridMask,eachDimSize);
+        gridMask = true;
+    }
+
+    int ssize = gridMask.numElements();
+    idx = 0;
+    while (idx < ssize)
+    {
+      for (int i = 0; i < dim-1; ++i)
+      {
+        if (idxLoop[i] == eachDimSize[i])
+        {
+          idxLoop[i] = 0;
+          ++idxLoop[i+1];
+        }
+      }
+
+      // Find out outer index
+      idxDomain = idxAxis = 0;
+      bool maskValue = true;
       for (int i = 0; i < numElement; ++i)
       {
-        indexMap[i] = idx;
-        if (2 == axisDomainOrder(i)) {
-            eachDimSize[indexMap[i]]   = domainP[idxDomain]->ni;
-            eachDimSize[indexMap[i]+1] = domainP[idxDomain]->nj;
-            idx += 2; ++idxDomain;
+        if (2 == axisDomainOrder(i))
+        {
+          int idxTmp = idxLoop[indexMap[i]] + idxLoop[indexMap[i]+1] * eachDimSize[indexMap[i]];
+          if (idxTmp < (*domainMasks[idxDomain]).numElements())
+            maskValue = maskValue && (*domainMasks[idxDomain])(idxTmp);
+          else
+            maskValue = false;
+          ++idxDomain;
         }
-        else if (1 == axisDomainOrder(i)) {
-  //        eachDimSize[indexMap[i]] = axisMasks[idxAxis]->numElements();
-          eachDimSize[indexMap[i]] = axisP[idxAxis]->n;
-          ++idx; ++idxAxis;
+        else if (1 == axisDomainOrder(i))
+        {
+          int idxTmp = idxLoop[indexMap[i]];
+          if (idxTmp < (*axisMasks[idxAxis]).numElements())
+            maskValue = maskValue && (*axisMasks[idxAxis])(idxTmp);
+          else
+            maskValue = false;
+
+          ++idxAxis;
         }
-        else {};
-      }
-
-//      if (!gridMask.isEmpty() && !createMask)
-      if (!createMask)
+      }
+
+      int maskIndex = idxLoop[0];
+      int mulDim = 1;
+      for (int k = 1; k < dim; ++k)
       {
-        for (int i = 0; i < dim; ++i)
-        {
-          if (gridMask.extent(i) != eachDimSize[i])
-            ERROR("CGrid::checkMask(void)",
-                  << "The mask has one dimension whose size is different from the one of the local grid." << std::endl
-                  << "Local size of dimension " << i << " is " << eachDimSize[i] << "." << std::endl
-                  << "Mask size for dimension " << i << " is " << gridMask.extent(i) << "." << std::endl
-                  << "Grid = " << this->getId())
-        }
-      }
-      else {
-          CArrayBoolTraits<CArray<bool,N> >::resizeArray(gridMask,eachDimSize);
-          gridMask = true;
-      }
-
-      int ssize = gridMask.numElements();
-      idx = 0;
-      while (idx < ssize)
-      {
-        for (int i = 0; i < dim-1; ++i)
-        {
-          if (idxLoop[i] == eachDimSize[i])
-          {
-            idxLoop[i] = 0;
-            ++idxLoop[i+1];
-          }
-        }
-
-        // Find out outer index
-        idxDomain = idxAxis = 0;
-        bool maskValue = true;
-        for (int i = 0; i < numElement; ++i)
-        {
-          if (2 == axisDomainOrder(i))
-          {
-            int idxTmp = idxLoop[indexMap[i]] + idxLoop[indexMap[i]+1] * eachDimSize[indexMap[i]];
-            if (idxTmp < (*domainMasks[idxDomain]).numElements())
-              maskValue = maskValue && (*domainMasks[idxDomain])(idxTmp);
-            else
-              maskValue = false;
-            ++idxDomain;
-          }
-          else if (1 == axisDomainOrder(i))
-          {
-            int idxTmp = idxLoop[indexMap[i]];
-            if (idxTmp < (*axisMasks[idxAxis]).numElements())
-              maskValue = maskValue && (*axisMasks[idxAxis])(idxTmp);
-            else
-              maskValue = false;
-
-            ++idxAxis;
-          }
-        }
-
-        int maskIndex = idxLoop[0];
-        int mulDim = 1;
-        for (int k = 1; k < dim; ++k)
-        {
-          mulDim *= eachDimSize[k-1];
-          maskIndex += idxLoop[k]*mulDim;
-        }
-        *(gridMask.dataFirst()+maskIndex) &= maskValue;
-
-        ++idxLoop[0];
-        ++idx;
-      }
-     }
+        mulDim *= eachDimSize[k-1];
+        maskIndex += idxLoop[k]*mulDim;
+      }
+      *(gridMask.dataFirst()+maskIndex) &= maskValue;
+
+      ++idxLoop[0];
+      ++idx;
+    }
    }