Local Cartesian Index for level grids

CMakeLists_files.cmake

@@ -122,6 +122,7 @@ if(Boost_VERSION_STRING VERSION_GREATER 1.53)
 	  tests/cpgrid/geometry_test.cpp
 	  tests/cpgrid/grid_lgr_test.cpp
 	  tests/cpgrid/inactiveCell_lgr_test.cpp
+	  tests/cpgrid/lgr_cartesian_idx_test.cpp
 	  tests/cpgrid/lookUpCellCentroid_cpgrid_test.cpp
 	  tests/cpgrid/lookupdataCpGrid_test.cpp
 	  tests/cpgrid/shifted_cart_test.cpp

opm/grid/CpGrid.hpp

@@ -1088,7 +1088,14 @@ namespace Dune
 
     public:
 
-
+        /// @brief Compute for each level grid, a map from the global_cell_[ cell index in level grid ] to the leaf index of the equivalent cell
+        ///        on the leaf grid view.
+        ///        Notice that cells that vanished and do not appear on the leaf grid view will not be considered.
+        ///        global_cell_[ cell index in level grid ] coincide with (local) Cartesian Index.
+        std::vector<std::unordered_map<std::size_t, std::size_t>> mapLocalCartesianIndexSetsToLeafIndexSet() const;
+
+        /// @brief Reverse map: from leaf index cell to { level, local/level Cartesian index of the cell }
+        std::vector<std::array<int,2>> mapLeafIndexSetToLocalCartesianIndexSets() const;
 
         /// \brief Size of the overlap on the leaf level
         unsigned int overlapSize(int) const;

opm/grid/common/CartesianIndexMapper.hpp

@@ -43,12 +43,6 @@ namespace Dune
             return 0;
         }
 
-        /** \brief return number of cells in the active level zero grid. Only relevant for CpGrid specialization. */
-        int compressedLevelZeroSize() const
-        {
-            return 0;
-        }
-
         /** \brief return index of the cells in the logical Cartesian grid */
         int cartesianIndex( const int /* compressedElementIndex */) const
         {
@@ -60,11 +54,25 @@ namespace Dune
         {
         }
 
+        /// Only relevant for CpGrid specialization.
+        /** \brief return number of cells in the active level zero grid. Only relevant for CpGrid specialization. */
+        int compressedLevelZeroSize() const
+        {
+            return 0;
+        }
+
         /** \brief return Cartesian coordinate, i.e. IJK, for a given cell. Only relevant for CpGrid specialization.*/
         void cartesianCoordinateLevel(const int /* compressedElementIndexOnLevel */,
                                       std::array<int,dimension>& /* coordsOnLevel */, int /*level*/) const
         {
         }
+
+        /** \brief return index of the cells in the logical Cartesian grid, for refined level grids. Only relevant for CpGrid specialization. */
+        int cartesianIndexLevel( const int /* compressedElementIndex */ , const int level) const
+        {
+            return 0;
+        }
+        /// Only relevant for CpGrid specialization. END
     };
 
 } // end namespace Opm

opm/grid/cpgrid/CartesianIndexMapper.hpp

@@ -50,11 +50,6 @@ namespace Dune
             return grid_.globalCell().size();
         }
 
-        int compressedLevelZeroSize() const
-        {
-            return (*grid_.currentData()[0]).size(0);
-        }
-
         int cartesianIndex( const int compressedElementIndex ) const
         {
             assert(  compressedElementIndex >= 0 && compressedElementIndex < compressedSize() );
@@ -66,13 +61,29 @@ namespace Dune
             grid_.getIJK( compressedElementIndex, coords );
         }
 
+        /// Additional methods related to LGRs
+        int compressedLevelZeroSize() const
+        {
+            return (*grid_.currentData()[0]).size(0);
+        }
+
+        int cartesianIndexLevel( const int compressedElementIndex, const int level) const
+        {
+            if ((level < 0) || (level > grid_.maxLevel())) {
+                throw std::invalid_argument("Invalid level.\n");
+            }
+            assert(  compressedElementIndex >= 0 && compressedElementIndex <  grid_.currentData()[level]->size(0) );
+            return grid_.currentData()[level]->globalCell()[compressedElementIndex];
+        }
+
         void cartesianCoordinateLevel(const int compressedElementIndexOnLevel, std::array<int,dimension>& coordsOnLevel, int level) const
         {
             if ((level < 0) || (level > grid_.maxLevel())) {
                 throw std::invalid_argument("Invalid level.\n");
             }
-            (*grid_.currentData()[level]).getIJK( compressedElementIndexOnLevel, coordsOnLevel);
+            grid_.currentData()[level]->getIJK( compressedElementIndexOnLevel, coordsOnLevel);
         }
+        /// Additional methods related to LGRs. END
     };
 
 } // end namespace Opm

opm/grid/cpgrid/CpGrid.cpp

@@ -679,6 +679,26 @@ void CpGrid::computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell
     }
 }
 
+std::vector<std::unordered_map<std::size_t, std::size_t>> CpGrid::mapLocalCartesianIndexSetsToLeafIndexSet() const
+{
+    std::vector<std::unordered_map<std::size_t, std::size_t>> localCartesianIdxSets_to_leafIdx(maxLevel()+1); // Plus level 0
+    for (const auto& element : elements(leafGridView())) {
+        const auto& global_cell_level = currentData()[element.level()]->globalCell()[element.getEquivLevelElem().index()];
+        localCartesianIdxSets_to_leafIdx[element.level()][global_cell_level] = element.index();
+    }
+    return localCartesianIdxSets_to_leafIdx;
+}
+
+std::vector<std::array<int,2>> CpGrid::mapLeafIndexSetToLocalCartesianIndexSets() const
+{
+    std::vector<std::array<int,2>> leafIdx_to_localCartesianIdxSets(currentData().back()->size(0));
+    for (const auto& element : elements(leafGridView())) {
+        const auto& global_cell_level = currentData()[element.level()]->globalCell()[element.getEquivLevelElem().index()];
+        leafIdx_to_localCartesianIdxSets[element.index()] = {element.level(), global_cell_level};
+    }
+    return leafIdx_to_localCartesianIdxSets;
+}
+
 void CpGrid::getIJK(const int c, std::array<int,3>& ijk) const
 {
     current_view_data_->getIJK(c, ijk);
@@ -1996,6 +2016,8 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
     computeGlobalCellLeafGridViewWithLgrs(global_cell_leaf);
     (*data[levels + preAdaptMaxLevel +1]).global_cell_.swap(global_cell_leaf);
 
+    mapLocalCartesianIndexSetsToLeafIndexSet();
+
     updateCornerHistoryLevels(cornerInMarkedElemWithEquivRefinedCorner,
                               elemLgrAndElemLgrCorner_to_refinedLevelAndRefinedCorner,
                               adaptedCorner_to_elemLgrAndElemLgrCorner,

opm/grid/cpgrid/CpGridData.hpp

@@ -306,8 +306,13 @@ class CpGridData
         ijk = getIJK(global_cell_[c], logical_cartesian_size_);
     }
 
+    const std::vector<int>& globalCell() const
+    {
+        return  global_cell_;
+    }
+
     /// @brief Extract Cartesian index triplet (i,j,k) given an index between 0 and NXxNYxNZ -1
-    ///        where NX, NY, and NZ is the total amoung of cells in each direction x-,y-,and z- respectively.
+    ///    where NX, NY, and NZ is the total amoung of cells in each direction x-,y-,and z- respectively.
     ///
     /// @param [in] idx      Integer between 0 and cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2]-1
     /// @param [in] cells_per_dim

opm/grid/polyhedralgrid/cartesianindexmapper.hh

@@ -31,7 +31,7 @@ namespace Dune
 
         const std::array<int, dimension>& cartesianDimensions() const
         {
-          return grid_.logicalCartesianSize();
+            return grid_.logicalCartesianSize();
         }
 
         int cartesianSize() const
@@ -44,12 +44,6 @@ namespace Dune
             return grid_.size( 0 );
         }
 
-        // Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
-        int compressedLevelZeroSize() const
-        {
-            return grid_.size( 0 );
-        }
-
         int cartesianIndex( const int compressedElementIndex ) const
         {
             assert( compressedElementIndex >= 0 && compressedElementIndex < compressedSize() );
@@ -73,6 +67,19 @@ namespace Dune
               coords[ 0 ] = gc ;
         }
 
+        // Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
+        int compressedLevelZeroSize() const
+        {
+            return grid_.size( 0 );
+        }
+        // Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
+        int cartesianIndexLevel( const int compressedElementIndex, const int level ) const
+        {
+            if (level) {
+                throw std::invalid_argument("Invalid level.\n");
+            }
+            return cartesianIndex(compressedElementIndex);
+        }
         // Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
         void cartesianCoordinateLevel(const int compressedElementIndexOnLevel, std::array<int,dimension>& coordsOnLevel, int level) const
         {