Retrieve local ijk's of an LGR, given its name

CMakeLists_files.cmake

@@ -37,6 +37,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/grid/cpgrid/Iterators.cpp
   opm/grid/cpgrid/Indexsets.cpp
   opm/grid/cpgrid/PartitionTypeIndicator.cpp
+  opm/grid/cpgrid/CpGridUtilities.cpp
   opm/grid/cpgrid/processEclipseFormat.cpp
   opm/grid/common/GeometryHelpers.cpp
   opm/grid/common/GridPartitioning.cpp
@@ -115,6 +116,7 @@ if(Boost_VERSION_STRING VERSION_GREATER 1.53)
 	  tests/cpgrid/grid_lgr_test.cpp
 	  tests/cpgrid/inactiveCell_lgr_test.cpp
 	  tests/cpgrid/lgr_cartesian_idx_test.cpp
+	  tests/cpgrid/lgrIJK_test.cpp
 	  tests/cpgrid/lookUpCellCentroid_cpgrid_test.cpp
 	  tests/cpgrid/lookupdataCpGrid_test.cpp
 	  tests/cpgrid/replace_lgr1_corner_idx_by_lgr2_corner_idx_test.cpp
@@ -179,6 +181,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/grid/cpgrid/CartesianIndexMapper.hpp
   opm/grid/cpgrid/CpGridData.hpp
   opm/grid/cpgrid/CpGridDataTraits.hpp
+  opm/grid/cpgrid/CpGridUtilities.hpp
   opm/grid/cpgrid/DataHandleWrappers.hpp
   opm/grid/cpgrid/DefaultGeometryPolicy.hpp
   opm/grid/cpgrid/dgfparser.hh

opm/grid/cpgrid/CpGridUtilities.cpp

@@ -0,0 +1,73 @@
+/*
+  Copyright 2025
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/grid/cpgrid/CpGridUtilities.hpp>
+#include <opm/grid/cpgrid/LevelCartesianIndexMapper.hpp>
+
+#include <array>
+#include <string>
+#include <vector>
+
+namespace Opm {
+
+std::vector<std::array<int,3>> lgrIJK(const Dune::CpGrid& grid, const std::string& lgr_name)
+{
+    // Check 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 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
+
+    // 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::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).
+
+        const int lgr_cart_index = element.getLevelCartesianIdx();
+        lgrIJK[lgr_cart_index] = ijk;
+    }
+    return lgrIJK;
+}
+
+}

opm/grid/cpgrid/CpGridUtilities.hpp

@@ -0,0 +1,41 @@
+/*
+  Copyright 2025 ....
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_CPGRIDUTILITIES_HEADER_INCLUDED
+#define OPM_CPGRIDUTILITIES_HEADER_INCLUDED
+
+#include <opm/grid/CpGrid.hpp>
+
+namespace Opm
+{
+
+/// @brief Retrieves the local grid refinement (LGR) IJKs of cells.
+///
+/// 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}.
+///
+/// @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);
+
+}
+
+#endif // OPM_CPGRIDUTILITIES_HEADER_INCLUDED