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 <boost/range.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;

        typedef boost::iterator_range<CellIter> Range;

        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