RegionMapping.hpp

Go to the documentation of this file.

/*
  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 <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED
#define OPM_REGIONMAPPING_HEADER_INCLUDED
 
#include <opm/grid/utility/IteratorRange.hpp>
 
#include <unordered_map>
#include <vector>
 
namespace Opm
{
 
    template < class Region = std::vector<int> >
    class RegionMapping {
    public:
        explicit
        RegionMapping(const Region& reg)
            : reg_(reg)
        {
            rev_.init(reg_);
        }
 
        typedef typename Region::value_type RegionId;
 
        typedef typename Region::size_type CellId;
 
        typedef typename std::vector<CellId>::const_iterator CellIter;
 
        using Range = iterator_range<CellIter>;
 
        RegionId
        region(const CellId c) const { return reg_[c]; }
 
        const std::vector<RegionId>&
        activeRegions() const
        {
            return rev_.active;
        }
 
        Range
        cells(const RegionId r) const {
            const auto id = rev_.binid.find(r);
 
            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]);
        }
 
    private:
        Region reg_;
 
        struct {
            typedef typename std::vector<CellId>::size_type Pos;
 
            std::unordered_map<RegionId, Pos> binid;
            std::vector<RegionId>             active;
 
            std::vector<Pos>    p;   
            std::vector<CellId> c;   
            void
            init(const Region& reg)
            {
                binid.clear();
                for (const auto& r : reg) {
                    ++binid[r];
                }
 
                p     .clear();  p.emplace_back(0);
                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] == static_cast<Pos>(reg.size()));
 
                c.resize(reg.size());
                {
                    CellId i = 0;
                    for (const auto& r : reg) {
                        auto& pos  = p[ binid[r] + 1 ];
                        c[ pos++ ] = i++;
                    }
                }
 
                p[0] = 0;
            }
        } rev_; 
    };
 
} // namespace Opm
 
#endif // OPM_REGIONMAPPING_HEADER_INCLUDED