RegionMapping<>: Support arbitrary region IDs

This commit introduces a new public method, activeRegions(), that
retrieves those region IDs that contain at least one active cell.
We furthermore extend the cells() method to support lookup of
arbitrary region IDs.  Non-active region IDs produce empty cell
ranges.

Intended use case is

    for (const auto& reg : rmap.activeRegions()) {
        const auto& c = rmap.cells(reg);

        // use c
    }

opm/core/simulator/initStateEquil.hpp

@@ -399,13 +399,10 @@ namespace Opm
                                  const Grid&                       G    ,
                                  const double grav)
                 {
-                    for (typename RMap::RegionId
+                    for (const auto& r : reg.activeRegions()) {
-                             r = 0, nr = reg.numRegions();
+                        const auto& cells = reg.cells(r);
-                         r < nr; ++r)
-                    {
-                        const typename RMap::CellRange cells = reg.cells(r);
-
                         const int repcell = *cells.begin();
+
                         const RhoCalc calc(props, repcell);
                         const EqReg eqreg(rec[r], calc,
                                           rs_func_[r], rv_func_[r],

opm/core/utility/RegionMapping.hpp

@@ -20,6 +20,9 @@
 #ifndef OPM_REGIONMAPPING_HEADER_INCLUDED
 #define OPM_REGIONMAPPING_HEADER_INCLUDED
 
+#include <boost/range.hpp>
+
+#include <unordered_map>
 #include <vector>
 
 namespace Opm
@@ -70,58 +73,7 @@ namespace Opm
          */
         typedef typename std::vector<CellId>::const_iterator CellIter;
 
-        /**
+        typedef boost::iterator_range<CellIter> Range;
-         * Range of cells.  Result from reverse (region-to-cell)
-         * mapping.
-         */
-        class CellRange {
-        public:
-            /**
-             * Constructor.
-             *
-             * \param[in] b Beginning of range.
-             * \param[in] e One past end of range.
-             */
-            CellRange(const CellIter b,
-                      const CellIter e)
-                : b_(b), e_(e)
-            {}
-
-            /**
-             * Read-only iterator on cell ranges.
-             */
-            typedef CellIter const_iterator;
-
-            /**
-             * Size type for this range.
-             */
-            typedef typename std::vector<CellId>::size_type size_type;
-
-            /**
-             * Beginning of cell range.
-             */
-            const_iterator begin() const { return b_; }
-
-            /**
-             * One past end of cell range.
-             */
-            const_iterator end()   const { return e_; }
-
-            /**
-             * Number of elements in the range.
-             */
-            size_type size() const { return e_ - b_; }
-
-        private:
-            const_iterator b_;
-            const_iterator e_;
-        };
-
-        /**
-         * Number of declared regions in cell-to-region mapping.
-         */
-        RegionId
-        numRegions() const { return RegionId(rev_.p.size()) - 1; }
 
         /**
          * Compute region number of given active cell.
@@ -132,16 +84,32 @@ namespace Opm
         RegionId
         region(const CellId c) const { return reg_[c]; }
 
+        const std::vector<RegionId>&
+        activeRegions() const
+        {
+            return rev_.active;
+        }
+
         /**
          * Extract active cells in particular region.
          *
          * \param[in] r Region number
-         * \returns Range of active cells in region @c r.
+         *
+         * \return Range of active cells in region @c r.  Empty if @c r is
+         * not an active region.
          */
-        CellRange
+        Range
         cells(const RegionId r) const {
-            const RegionId i = r - rev_.low;
+            const auto id = rev_.binid.find(r);
-            return CellRange(rev_.c.begin() + rev_.p[i + 0],
+
+            if (id == rev_.binid.end()) {
+                // Region 'r' not an active region.  Return empty.
+                return Range(rev_.c.end(), rev_.c.end());
+            }
+
+            const auto i = id->second;
+
+            return Range(rev_.c.begin() + rev_.p[i + 0],
                          rev_.c.begin() + rev_.p[i + 1]);
         }
 
@@ -156,9 +124,12 @@ namespace Opm
          */
         struct {
             typedef typename std::vector<CellId>::size_type Pos;
+
+            std::unordered_map<RegionId, Pos> binid;
+            std::vector<RegionId>             active;
+
             std::vector<Pos>    p;   /**< Region start pointers */
             std::vector<CellId> c;   /**< Region cells */
-            RegionId            low; /**< Smallest region number */
 
             /**
              * Compute reverse mapping.  Standard linear insertion
@@ -167,32 +138,37 @@ namespace Opm
             void
             init(const Region& reg)
             {
-                typedef typename Region::const_iterator CI;
+                binid.clear();
-                const std::pair<CI,CI>
+                for (const auto& r : reg) {
-                    m = std::minmax_element(reg.begin(), reg.end());
+                    ++binid[r];
-
-                low  = *m.first;
-
-                const typename Region::size_type
-                    n = *m.second - low + 1;
-
-                p.resize(n + 1);  std::fill(p.begin(), p.end(), Pos(0));
-                for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
-                    p[ reg[i] - low + 1 ] += 1;
                 }
 
-                for (typename std::vector<Pos>::size_type
+                p     .clear();  p.emplace_back(0);
-                         i = 1, sz = p.size(); i < sz; ++i) {
+                active.clear();
+                {
+                    Pos n = 0;
+                    for (auto& id : binid) {
+                        active.push_back(id.first);
+                        p     .push_back(id.second);
+
+                        id.second = n++;
+                    }
+                }
+
+                for (decltype(p.size()) i = 1, sz = p.size(); i < sz; ++i) {
                     p[0] += p[i];
                     p[i]  = p[0] - p[i];
                 }
 
-                assert (p[0] ==
+                assert (p[0] == static_cast<Pos>(reg.size()));
-                        static_cast<typename Region::size_type>(reg.size()));
 
                 c.resize(reg.size());
-                for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
+                {
-                    c[ p[ reg[i] - low + 1 ] ++ ] = i;
+                    CellId i = 0;
+                    for (const auto& r : reg) {
+                        auto& pos  = p[ binid[r] + 1 ];
+                        c[ pos++ ] = i++;
+                    }
                 }
 
                 p[0] = 0;

tests/test_equil.cpp

@@ -307,10 +307,8 @@ BOOST_AUTO_TEST_CASE (RegMapping)
     Opm::RegionMapping<> eqlmap(eqlnum);
 
     PPress ppress(2, PVal(G->number_of_cells, 0));
-    for (int r = 0, e = eqlmap.numRegions(); r != e; ++r)
+    for (const auto& r : eqlmap.activeRegions()) {
-    {
+        const auto& rng = eqlmap.cells(r);
-        const Opm::RegionMapping<>::CellRange&
-            rng = eqlmap.cells(r);
 
         const int    rno  = r;
         const double grav = 10;
@@ -318,14 +316,13 @@ BOOST_AUTO_TEST_CASE (RegMapping)
             Opm::Equil::phasePressures(*G, region[rno], rng, grav);
 
         PVal::size_type i = 0;
-        for (Opm::RegionMapping<>::CellRange::const_iterator
+        for (const auto& c : rng) {
-                 c = rng.begin(), ce = rng.end();
-             c != ce; ++c, ++i)
-        {
             assert (i < p[0].size());
 
-            ppress[0][*c] = p[0][i];
+            ppress[0][c] = p[0][i];
-            ppress[1][*c] = p[1][i];
+            ppress[1][c] = p[1][i];
+
+            ++i;
         }
     }
 

tests/test_regionmapping.cpp

@@ -26,42 +26,130 @@
 
 #define NVERBOSE  // Suppress own messages when throw()ing
 
-#define BOOST_TEST_MODULE UnitsTest
+#define BOOST_TEST_MODULE RegionMapping
 
 #include <opm/core/utility/platform_dependent/disable_warnings.h>
 #include <boost/test/unit_test.hpp>
-#include <boost/test/floating_point_comparison.hpp>
 #include <opm/core/utility/platform_dependent/reenable_warnings.h>
 
 /* --- our own headers --- */
 
 #include <opm/core/utility/RegionMapping.hpp>
 
+#include <algorithm>
+#include <map>
 
 BOOST_AUTO_TEST_SUITE ()
 
 
-BOOST_AUTO_TEST_CASE (RegionMapping)
+BOOST_AUTO_TEST_CASE (Forward)
+{
+    std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
+
+    Opm::RegionMapping<> rm(regions);
+
+    for (decltype(regions.size()) i = 0, n = regions.size(); i < n; ++i) {
+        BOOST_CHECK_EQUAL(rm.region(i), regions[i]);
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE (ActiveRegions)
 {
     //                           0  1  2  3  4  5  6  7  8
     std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
+
     Opm::RegionMapping<> rm(regions);
-    for (size_t i = 0; i < regions.size(); ++i) {
+
-        BOOST_CHECK_EQUAL(rm.region(i), regions[i]);
-    }
     std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
-    std::vector< std::vector<int> > region_cells = {  { 0, 2, 4 },  { 7 },  { 3 },  { 1 },  { 6, 8 },  { 5 }  };
+
-    BOOST_REQUIRE_EQUAL(rm.numRegions(), region_ids.size());
+    auto active = [&region_ids](const int reg)
-    for (size_t i = 0; i < region_ids.size(); ++i) {
+        {
-        auto cells = rm.cells(region_ids[i]);
+            auto b = region_ids.begin();
-        BOOST_REQUIRE_EQUAL(std::distance(cells.begin(), cells.end()), region_cells[i].size());
+            auto e = region_ids.end();
-        size_t count = 0;
+
-        for (auto iter = cells.begin(); iter != cells.end(); ++iter, ++count) {
+            return std::find(b, e, reg) != e;
-            BOOST_CHECK_EQUAL(*iter, region_cells[i][count]);
+        };
-        }
+
+    BOOST_CHECK_EQUAL(rm.activeRegions().size(), region_ids.size());
+
+    for (const auto& reg : rm.activeRegions()) {
+        BOOST_CHECK(active(reg));
     }
 }
 
 
+BOOST_AUTO_TEST_CASE (Consecutive)
+{
+    using RegionCells = std::map<int, std::vector<int>>;
+
+    //                           0  1  2  3  4  5  6  7  8
+    std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
+
+    Opm::RegionMapping<> rm(regions);
+
+    std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
+    RegionCells      region_cells;
+    {
+        using VT = RegionCells::value_type;
+
+        region_cells.insert(VT(2, { 0, 2, 4 }));
+        region_cells.insert(VT(3, { 7 }));
+        region_cells.insert(VT(4, { 3 }));
+        region_cells.insert(VT(5, { 1 }));
+        region_cells.insert(VT(6, { 6, 8 }));
+        region_cells.insert(VT(7, { 5 }));
+    }
+
+    for (const auto& reg : region_ids) {
+        const auto& cells  = rm.cells(reg);
+        const auto& expect = region_cells[reg];
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(),
+                                      expect.begin(), expect.end());
+    }
+
+    // Verify that there are no cells in unused regions 0 and 1.
+    for (const auto& r : { 0, 1 }) {
+        BOOST_CHECK(rm.cells(r).empty());
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE (NonConsecutive)
+{
+    using RegionCells = std::map<int, std::vector<int>>;
+
+    //                           0  1  2  3  4  5  6  7  8
+    std::vector<int> regions = { 2, 4, 2, 4, 2, 7, 6, 3, 6 };
+
+    Opm::RegionMapping<> rm(regions);
+
+    std::vector<int> region_ids = { 2, 3, 4, 6, 7 };
+    RegionCells      region_cells;
+    {
+        using VT = RegionCells::value_type;
+
+        region_cells.insert(VT(2, { 0, 2, 4 }));
+        region_cells.insert(VT(3, { 7 }));
+        region_cells.insert(VT(4, { 1, 3 }));
+        region_cells.insert(VT(6, { 6, 8 }));
+        region_cells.insert(VT(7, { 5 }));
+    }
+
+    for (const auto& reg : region_ids) {
+        const auto& cells  = rm.cells(reg);
+        const auto& expect = region_cells[reg];
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(),
+                                      expect.begin(), expect.end());
+    }
+
+    // Verify that there are no cells in unused regions 0, 1, and 5.
+    for (const auto& r : { 0, 1, 5 }) {
+        BOOST_CHECK(rm.cells(r).empty());
+    }
+}
+
 
 BOOST_AUTO_TEST_SUITE_END()