diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2968d25f..010cb7f4 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -170,6 +170,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_wellsmanager.cpp tests/test_wellcontrols.cpp tests/test_equil.cpp + tests/test_regionmapping.cpp ) # originally generated with the command: @@ -364,6 +365,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/utility/NonuniformTableLinear.hpp opm/core/utility/NullStream.hpp opm/core/utility/PolynomialUtils.hpp + opm/core/utility/RegionMapping.hpp opm/core/utility/RootFinders.hpp opm/core/utility/SparseTable.hpp opm/core/utility/SparseVector.hpp diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 4ed2b619..fd1c6721 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -260,170 +261,6 @@ namespace Opm } // namespace miscibility - /** - * Forward and reverse mappings between cells and - * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', - * 'PVTNUM', or 'EQUILNUM' arrays). - * - * \tparam Region Type of a forward region mapping. Expected - * to provide indexed access through - * operator[]() as well as inner types - * 'value_type', 'size_type', and - * 'const_iterator'. - */ - template < class Region = std::vector > - class RegionMapping { - public: - /** - * Constructor. - * - * \param[in] reg Forward region mapping, restricted to - * active cells only. - */ - explicit - RegionMapping(const Region& reg) - : reg_(reg) - { - rev_.init(reg_); - } - - /** - * Type of forward (cell-to-region) mapping result. - * Expected to be an integer. - */ - typedef typename Region::value_type RegionId; - - /** - * Type of reverse (region-to-cell) mapping (element) - * result. - */ - typedef typename Region::size_type CellId; - - /** - * Type of reverse region-to-cell range bounds and - * iterators. - */ - typedef typename std::vector::const_iterator CellIter; - - /** - * 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; - - /** - * Beginning of cell range. - */ - const_iterator begin() const { return b_; } - - /** - * One past end of cell range. - */ - const_iterator end() const { return e_; } - - 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. - * - * \param[in] c Active cell - * \return Region to which @c c belongs. - */ - RegionId - region(const CellId c) const { return reg_[c]; } - - /** - * Extract active cells in particular region. - * - * \param[in] r Region number - * \returns Range of active cells in region @c r. - */ - CellRange - cells(const RegionId r) const { - const RegionId i = r - rev_.low; - return CellRange(rev_.c.begin() + rev_.p[i + 0], - rev_.c.begin() + rev_.p[i + 1]); - } - - private: - /** - * Copy of forward region mapping (cell-to-region). - */ - Region reg_; - - /** - * Reverse mapping (region-to-cell). - */ - struct { - typedef typename std::vector::size_type Pos; - std::vector p; /**< Region start pointers */ - std::vector c; /**< Region cells */ - RegionId low; /**< Smallest region number */ - - /** - * Compute reverse mapping. Standard linear insertion - * sort algorithm. - */ - void - init(const Region& reg) - { - typedef typename Region::const_iterator CI; - const std::pair - m = std::minmax_element(reg.begin(), reg.end()); - - 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::size_type - i = 1, sz = p.size(); i < sz; ++i) { - p[0] += p[i]; - p[i] = p[0] - p[i]; - } - - assert (p[0] == - static_cast(reg.size())); - - c.resize(reg.size()); - for (CellId i = 0, nc = reg.size(); i < nc; ++i) { - c[ p[ reg[i] - low + 1 ] ++ ] = i; - } - - p[0] = 0; - } - } rev_; /**< Reverse mapping instance */ - }; /** * Equilibration record. diff --git a/opm/core/utility/RegionMapping.hpp b/opm/core/utility/RegionMapping.hpp new file mode 100644 index 00000000..9845919f --- /dev/null +++ b/opm/core/utility/RegionMapping.hpp @@ -0,0 +1,195 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + + 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 . +*/ + +#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED +#define OPM_REGIONMAPPING_HEADER_INCLUDED + +#include + +namespace Opm +{ + + /** + * Forward and reverse mappings between cells and + * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', + * 'PVTNUM', or 'EQUILNUM' arrays). + * + * \tparam Region Type of a forward region mapping. Expected + * to provide indexed access through + * operator[]() as well as inner types + * 'value_type', 'size_type', and + * 'const_iterator'. + */ + template < class Region = std::vector > + class RegionMapping { + public: + /** + * Constructor. + * + * \param[in] reg Forward region mapping, restricted to + * active cells only. + */ + explicit + RegionMapping(const Region& reg) + : reg_(reg) + { + rev_.init(reg_); + } + + /** + * Type of forward (cell-to-region) mapping result. + * Expected to be an integer. + */ + typedef typename Region::value_type RegionId; + + /** + * Type of reverse (region-to-cell) mapping (element) + * result. + */ + typedef typename Region::size_type CellId; + + /** + * Type of reverse region-to-cell range bounds and + * iterators. + */ + typedef typename std::vector::const_iterator CellIter; + + /** + * 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; + + /** + * Beginning of cell range. + */ + const_iterator begin() const { return b_; } + + /** + * One past end of cell range. + */ + const_iterator end() const { return e_; } + + 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. + * + * \param[in] c Active cell + * \return Region to which @c c belongs. + */ + RegionId + region(const CellId c) const { return reg_[c]; } + + /** + * Extract active cells in particular region. + * + * \param[in] r Region number + * \returns Range of active cells in region @c r. + */ + CellRange + cells(const RegionId r) const { + const RegionId i = r - rev_.low; + return CellRange(rev_.c.begin() + rev_.p[i + 0], + rev_.c.begin() + rev_.p[i + 1]); + } + + private: + /** + * Copy of forward region mapping (cell-to-region). + */ + Region reg_; + + /** + * Reverse mapping (region-to-cell). + */ + struct { + typedef typename std::vector::size_type Pos; + std::vector p; /**< Region start pointers */ + std::vector c; /**< Region cells */ + RegionId low; /**< Smallest region number */ + + /** + * Compute reverse mapping. Standard linear insertion + * sort algorithm. + */ + void + init(const Region& reg) + { + typedef typename Region::const_iterator CI; + const std::pair + m = std::minmax_element(reg.begin(), reg.end()); + + 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::size_type + i = 1, sz = p.size(); i < sz; ++i) { + p[0] += p[i]; + p[i] = p[0] - p[i]; + } + + assert (p[0] == + static_cast(reg.size())); + + c.resize(reg.size()); + for (CellId i = 0, nc = reg.size(); i < nc; ++i) { + c[ p[ reg[i] - low + 1 ] ++ ] = i; + } + + p[0] = 0; + } + } rev_; /**< Reverse mapping instance */ + }; + +} // namespace Opm + +#endif // OPM_REGIONMAPPING_HEADER_INCLUDED diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index c4db1786..4519f5f5 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -287,12 +287,12 @@ BOOST_AUTO_TEST_CASE (RegMapping) G->cartdims, cdim, &cells[0], &eqlnum[0]); } - Opm::equil::RegionMapping<> eqlmap(eqlnum); + Opm::RegionMapping<> eqlmap(eqlnum); PPress ppress(2, PVal(G->number_of_cells, 0)); for (int r = 0, e = eqlmap.numRegions(); r != e; ++r) { - const Opm::equil::RegionMapping<>::CellRange& + const Opm::RegionMapping<>::CellRange& rng = eqlmap.cells(r); const int rno = r; @@ -301,7 +301,7 @@ BOOST_AUTO_TEST_CASE (RegMapping) Opm::equil::phasePressures(*G, region[rno], rng, grav); PVal::size_type i = 0; - for (Opm::equil::RegionMapping<>::CellRange::const_iterator + for (Opm::RegionMapping<>::CellRange::const_iterator c = rng.begin(), ce = rng.end(); c != ce; ++c, ++i) { diff --git a/tests/test_regionmapping.cpp b/tests/test_regionmapping.cpp new file mode 100644 index 00000000..ca2ad4ea --- /dev/null +++ b/tests/test_regionmapping.cpp @@ -0,0 +1,64 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + + 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 . +*/ + +#include "config.h" + +/* --- Boost.Test boilerplate --- */ +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE UnitsTest +#include +#include + +/* --- our own headers --- */ + +#include + + +BOOST_AUTO_TEST_SUITE () + + +BOOST_AUTO_TEST_CASE (RegionMapping) +{ + // 0 1 2 3 4 5 6 7 8 + std::vector 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 region_ids = { 2, 3, 4, 5, 6, 7 }; + std::vector< std::vector > region_cells = { { 0, 2, 4 }, { 7 }, { 3 }, { 1 }, { 6, 8 }, { 5 } }; + BOOST_REQUIRE_EQUAL(rm.numRegions(), region_ids.size()); + for (size_t i = 0; i < region_ids.size(); ++i) { + auto cells = rm.cells(region_ids[i]); + BOOST_REQUIRE_EQUAL(std::distance(cells.begin(), cells.end()), region_cells[i].size()); + size_t count = 0; + for (auto iter = cells.begin(); iter != cells.end(); ++iter, ++count) { + BOOST_CHECK_EQUAL(*iter, region_cells[i][count]); + } + } +} + + + +BOOST_AUTO_TEST_SUITE_END()