Modify lgrIJK and simplify its test

opm/grid/cpgrid/CpGridUtilities.cpp

@@ -22,52 +22,47 @@
 
 #include <array>
 #include <string>
+#include <tuple>
+#include <unordered_map>
 #include <vector>
 
-namespace Opm {
+namespace Opm
+{
 
-std::vector<std::array<int,3>> lgrIJK(const Dune::CpGrid& grid, const std::string& lgr_name)
+std::tuple<std::vector<int>, std::unordered_map<int, int>, std::vector<std::array<int, 3>>>
+lgrIJK(const Dune::CpGrid& grid, const std::string& lgr_name)
 {
-    // Check lgr_name exists in lgr_names_
+    // Check if lgr_name exists in lgr_names_
     const auto& lgr_names = grid.getLgrNameToLevel();
     auto it = lgr_names.find(lgr_name);
     if (it == lgr_names.end()) {
         OPM_THROW(std::runtime_error, "LGR name not found: " + lgr_name);
     }
 
+    const auto level = it->second;
     const Opm::LevelCartesianIndexMapper<Dune::CpGrid> levelCartMapper(grid);
-    const auto& level = it->second;
-
-    // Determine the logical Cartesian size of the LGR and total cells (including inactive ones)
-    const auto lgr_dim = grid.currentData()[level]->logicalCartesianSize();
-    const int lgr_cells = lgr_dim[0] * lgr_dim[1] * lgr_dim[2]; // (including inactive ones)
-    // Actual size is given by grid.levelGridView(level).size(0)
-
-    // Initialize level IJKs with inactive values
-    static constexpr std::array<int, 3> inactiveIJK = { -1, -1, -1 };
-    std::vector<std::array<int, 3>> lgrIJK(lgr_cells);
-    lgrIJK.assign(lgr_cells, inactiveIJK);  // Ensures all elements are set to inactive
+    const auto& levelView = grid.levelGridView(level);
+    const auto numCells = levelView.size(0);
 
-    // Iterate over active elements in the grid. Rewrite ijks at the specified refined level active cells.
-    //
-    // Note: to avoid go over all the leaf elements and replace this approach with the
-    // corresponding grid.levelGridView(level), a map relating level_to_lgr_idx is needed.
-    for (const auto& element : Dune::elements(grid.leafGridView())) {
-        if (element.level() != level) {
-            continue;
-        }
+    std::vector<std::array<int, 3>> lgrIJK(numCells);
+    std::vector<int> cellIdxToLgrCartesianIdx(numCells);
+    std::unordered_map<int, int> lgrCartesianIdxToCellIdx;
+    lgrCartesianIdxToCellIdx.reserve(numCells);
 
+    // Iterate over (active) elements in the grid and populate the structures
+    for (const auto& element : Dune::elements(grid.levelGridView(level))) {
         std::array<int, 3> ijk;
-        levelCartMapper.cartesianCoordinate(element.getLevelElem().index(), ijk, level);
-        // In general, element.getLevelElem().index() != element.getLevelCartesianIdx().
-        // - element.getLevelElem().index() : integer between 0 and "total active cells - 1" on the level grid.
-        // - element.getLevelCartesianIdx() : integer between 0 and "lgr_cells - 1", representing an index of
-        //                                    the underlying Cartesian grid for the level (with/without INACTIVE cells).
+        levelCartMapper.cartesianCoordinate(element.index(), ijk, level);
 
-        const int lgr_cart_index = element.getLevelCartesianIdx();
-        lgrIJK[lgr_cart_index] = ijk;
+        const int cellIndex = element.index();
+        const int cartesianIdx = element.getLevelCartesianIdx();
+
+        lgrIJK[cellIndex] = ijk;
+        cellIdxToLgrCartesianIdx[cellIndex] = cartesianIdx;
+        lgrCartesianIdxToCellIdx[cartesianIdx] = cellIndex;
     }
-    return lgrIJK;
-}
 
+    return std::make_tuple(cellIdxToLgrCartesianIdx, lgrCartesianIdxToCellIdx, lgrIJK);
 }
+
+} // namespace Opm

opm/grid/cpgrid/CpGridUtilities.hpp

@@ -25,17 +25,20 @@
 namespace Opm
 {
 
-/// @brief Retrieves the local grid refinement (LGR) IJKs of cells.
+/// @brief Retrieves Cartesian indices for a specified Local Grid Refinement (LGR) level in a Dune::CpGrid.
 ///
-/// This method returns a vector of IJK indices corresponding to underlying Cartesian grid of the local
-/// grid refinement (LGR) specified by its name. If the specified name is not found, an exception is thrown.
-/// For inactive cells—i.e., non-existing child cells of inactive coarse cells defined by the CARFIN keyword—
-/// the corresponding IJK index is set to {-1, -1, -1}.
+/// This function extracts the mapping between active cell indices and Cartesian indices for a given LGR name.
+/// If the specified name is not found, an exception is thrown.
 ///
-/// @param lgr_name The name of the LGR whose local IJK coordinates are requested.
-/// @return std::vector<std::array<int, 3>> A list of (i, j, k)'s for each cell in the LGR.
-std::vector<std::array<int,3>> lgrIJK(const Dune::CpGrid& grid, const std::string& lgr_name);
-
-}
+/// @param grid The Dune::CpGrid
+/// @param lgr_name The name of the LGR whose indices are to be retrieved.
+/// @return A tuple containing:
+///   - A std::vector<int> mapping cell indices to their corresponding Cartesian indices in the LGR.
+///   - A std::unordered_map<int, int> mapping Cartesian indices back to cell indices (handles inactive parent cells).
+///   - A std::vector<std::array<int, 3>> storing the (i, j, k) Cartesian coordinates for active cells.
+std::tuple<std::vector<int>, std::unordered_map<int, int>, std::vector<std::array<int, 3>>>
+lgrIJK(const Dune::CpGrid& grid, const std::string& lgr_name);
+
+} // namespace Opm
 
 #endif // OPM_CPGRIDUTILITIES_HEADER_INCLUDED

opm/grid/cpgrid/Entity.hpp

@@ -607,8 +607,7 @@ template<int codim>
 int Dune::cpgrid::Entity<codim>::getLevelCartesianIdx() const
 {
     const auto& level_data = (*(pgrid_ -> level_data_ptr_))[level()].get();
-    // getLevelElem() throws when the entity does not belong to the leaf grid view.
-    return level_data -> global_cell_[getLevelElem().index()];
+    return level_data -> global_cell_[getEquivLevelElem().index()];
 }
 
 } // namespace cpgrid