RepairZCORN.hpp

Go to the documentation of this file.

/*
  Copyright 2017 SINTEF ICT, Applied Mathematics.
  Copyright 2017 Statoil ASA.
 
  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_REPAIRZCORN_HEADER_INCLUDED
#define OPM_REPAIRZCORN_HEADER_INCLUDED
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <exception>
#include <iterator>
#include <numeric>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
 
namespace Opm { namespace UgGridHelpers {
 
    class RepairZCORN
    {
    public:
        template <class CartDims>
        RepairZCORN(std::vector<double>&&   zcorn,
                    const std::vector<int>& actnum,
                    const CartDims&         cartDims)
            : active_   (actnum, cartDims)
            , zcorn_idx_(cartDims)
            , zcorn_    (std::move(zcorn))
        {
            this->ensureZCornIsDepth();
            this->ensureTopNotBelowBottom();
            this->ensureBottomNotBelowLowerTop(cartDims[2]);
        }
 
        struct ZCornChangeCount
        {
            std::size_t cells{0};
 
            std::size_t corners{0};
        };
 
        std::vector<double> destructivelyGrabSanitizedValues()
        {
            return std::move(this->zcorn_);
        }
 
        bool switchedToDepth() const
        {
            return this->switchedToDepth_;
        }
 
        const ZCornChangeCount& statTopBelowBottom() const
        {
            return this->topBelowBottom_;
        }
 
        const ZCornChangeCount& statBottomBelowLowerTop() const
        {
            return this->bottomBelowLowerTop_;
        }
 
    private:
        class ActiveCells
        {
        public:
            template <class CartDims>
            ActiveCells(const std::vector<int>& actnum,
                        const CartDims&         cartDims)
                : nx_(cartDims[0])
                , ny_(cartDims[1])
                , nz_(cartDims[2])
            {
                const auto nglob = nx_ * ny_ * nz_;
 
                if (actnum.empty()) {
                    this->is_active_.resize(nglob, true);
                    this->acells_   .resize(nglob, 0);
 
                    std::iota(std::begin(this->acells_),
                              std::end  (this->acells_), 0);
                }
                else if (actnum.size() == nglob) {
                    this->is_active_.resize(nglob, false);
                    this->acells_   .reserve(nglob);
 
                    for (auto i = 0*nglob; i < nglob; ++i) {
                        if (actnum[i] != 0) {
                            this->acells_.push_back(i);
                            this->is_active_[i] = true;
                        }
                    }
                }
                else {
                    throw std::invalid_argument {
                        "ACTNUM vector does not match global size"
                    };
                }
            }
 
            using IndexTuple = std::array<std::size_t, 3>;
 
            IndexTuple getCellIJK(const std::size_t globCell) const
            {
                auto c = globCell;
 
                auto i = c % nx_;  c /= nx_;
                auto j = c % ny_;  c /= ny_;
                auto k = c % nz_;
 
                return {{ i, j, k }};
            }
 
            const std::vector<std::size_t>& activeGlobal() const
            {
                return this->acells_;
            }
 
            std::size_t numGlobalCells() const
            {
                return this->nx_ * this->ny_ * this->nz_;
            }
 
            std::size_t neighbourBelow(IndexTuple ijk) const
            {
                if (ijk[2] >= nz_ - 1) {
                    return -1;
                }
 
                ijk[2] += 1;
 
                auto below = this->globalCellIdx(ijk);
                while ((below < this->numGlobalCells()) &&
                       (! this->is_active_[below]))
                {
                    ijk[2] += 1;
 
                    below = this->globalCellIdx(ijk);
                }
 
                return below;
            }
 
        private:
            const std::size_t nx_;
 
            const std::size_t ny_;
 
            const std::size_t nz_;
 
            std::vector<bool> is_active_;
 
            std::vector<std::size_t> acells_;
 
            std::size_t globalCellIdx(const IndexTuple& ijk) const
            {
                if (ijk[2] > nz_ - 1) { return -1; }
 
                return ijk[0] + nx_*(ijk[1] + ny_*ijk[2]);
            }
        };
 
        class ZCornIndex
        {
        public:
            template <class CartDims>
            ZCornIndex(const CartDims& cartDims)
                : nx_         (cartDims[0])
                , ny_         (cartDims[1])
                , nz_         (cartDims[2])
                , layerOffset_((2 * nx_) * (2 * ny_))
            {}
 
            struct PillarPointIDX
            {
                std::size_t top;
 
                std::size_t bottom;
            };
 
            template <class IndexTuple>
            std::array<PillarPointIDX, 4>
            pillarPoints(const IndexTuple& ijk) const
            {
                const auto start = this->getStartIDX(ijk);
 
                return {
                    this->p00(start),
                    this->p10(start),
                    this->p01(start),
                    this->p11(start)
                };
            }
 
        private:
            const std::size_t nx_;
 
            const std::size_t ny_;
 
            const std::size_t nz_;
 
            const std::size_t layerOffset_;
 
            template <class IndexTuple>
            std::size_t getStartIDX(const IndexTuple& ijk) const
            {
                return 2*ijk[0] + 2*nx_*(2*ijk[1] + 2*ny_*(2 * ijk[2]));
            }
 
            PillarPointIDX p00(const std::size_t start) const
            {
                return this->pillarpts(start, this->offset(0, 0));
            }
 
            PillarPointIDX p10(const std::size_t start) const
            {
                return this->pillarpts(start, this->offset(1, 0));
            }
 
            PillarPointIDX p01(const std::size_t start) const
            {
                return this->pillarpts(start, this->offset(0, 1));
            }
 
            PillarPointIDX p11(const std::size_t start) const
            {
                return this->pillarpts(start, this->offset(1, 1));
            }
 
            std::size_t offset(const std::size_t i, const std::size_t j) const
            {
                assert ((i == 0) || (i == 1));
                assert ((j == 0) || (j == 1));
 
                return i + j*2*nx_;
            }
 
            PillarPointIDX pillarpts(const std::size_t start,
                                     const std::size_t off) const
            {
                return {
                    start + off,
                    start + off + this->layerOffset_
                };
            }
        };
 
        const ActiveCells active_;
 
        const ZCornIndex zcorn_idx_;
 
        std::vector<double> zcorn_;
 
        bool switchedToDepth_{false};
 
        ZCornChangeCount topBelowBottom_;
 
        ZCornChangeCount bottomBelowLowerTop_;
 
        void ensureZCornIsDepth()
        {
            if (this->zcornIsElevation()) {
                for (auto& zc : this->zcorn_) {
                    zc = - zc;
                }
 
                this->switchedToDepth_ = true;
            }
        }
 
        void ensureTopNotBelowBottom()
        {
            for (const auto& globCell : this->active_.activeGlobal()) {
                this->ensureTopNotBelowBottom(globCell);
            }
        }
 
        void ensureBottomNotBelowLowerTop(const std::size_t nz)
        {
            if (nz == 0) { return; }
 
            auto bottomLayer = [nz](const std::size_t layerID)
            {
                return layerID == (nz - 1);
            };
 
            for (const auto& globCell : this->active_.activeGlobal()) {
                const auto& ijk = this->active_.getCellIJK(globCell);
 
                if (bottomLayer(ijk[2])) { continue; }
 
                this->ensureCellBottomNotBelowLowerTop(ijk);
            }
        }
 
        template <class CellIndex>
        void ensureTopNotBelowBottom(const CellIndex globCell)
        {
            const auto cornerCnt0 = this->topBelowBottom_.corners;
 
            const auto ijk = this->active_.getCellIJK(globCell);
 
            for (const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
                const auto zb = this->zcorn_[pt.bottom];
                auto&      zt = this->zcorn_[pt.top];
 
                if (zt > zb) {  // Top below bottom (ZCORN is depth)
                    zt = zb;
 
                    this->topBelowBottom_.corners += 1;
                }
            }
 
            this->topBelowBottom_.cells +=
                this->topBelowBottom_.corners > cornerCnt0;
        }
 
        template <class IndexTuple>
        void ensureCellBottomNotBelowLowerTop(const IndexTuple& ijk)
        {
            const auto below = this->active_.neighbourBelow(ijk);
 
            if (below >= this->active_.numGlobalCells()) {
                return;
            }
 
            const auto cornerCnt0 = this->bottomBelowLowerTop_.corners;
 
            const auto& up    = this->zcorn_idx_.pillarPoints(ijk);
            const auto  d_ijk = this->active_.getCellIJK(below);
            const auto& down  = this->zcorn_idx_.pillarPoints(d_ijk);
 
            for (auto n = up.size(), i = 0*n; i < n; ++i) {
                const auto zbu = this->zcorn_[up  [i].bottom];
                auto&      ztd = this->zcorn_[down[i].top];
 
                if (zbu > ztd) { // Bottom below lower top (ZCORN is depth)
                    ztd = zbu;
 
                    this->bottomBelowLowerTop_.corners += 1;
                }
            }
 
            this->bottomBelowLowerTop_.cells +=
                this->bottomBelowLowerTop_.corners > cornerCnt0;
        }
 
        bool zcornIsElevation() const
        {
            auto all_signs = std::vector<int>{};
            all_signs.reserve(this->active_.numGlobalCells());
 
            for (const auto& globCell : this->active_.activeGlobal()) {
                all_signs.push_back(this->getZCornSign(globCell));
            }
 
            // Ignore twisted cells (i.e., cells of indeterminate signs).
            const int ignore = 0;
 
            // Elevation implies that ZCORN values are decreasing which
            // means that the signs in all non-twisted cells equal -1.
            //
            // Note: This check is *NOT* equivalent to
            //
            //    allEqual(all_signs, ignore, -1)
            //
            // because that check returns \c true if ALL cells are ignored
            // whereas first()==-1 ensures that there is at least ONE cell
            // of determinate sign.
            return (first(all_signs, ignore) == -1)
                && allEqual(all_signs, ignore, -1);
        }
 
        template <typename CellIndex>
        int getZCornSign(const CellIndex globCell) const
        {
            auto sign = [](const double x) -> int
            {
                return (x > 0.0) - (x < 0.0);
            };
 
            const auto ijk = this->active_.getCellIJK(globCell);
 
            auto sgn = std::vector<int>{};  sgn.reserve(4);
 
            for (const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
                const auto dz =
                    this->zcorn_[pt.bottom] - this->zcorn_[pt.top];
 
                sgn.push_back(sign(dz));
            }
 
            const int ignore = 0;
 
            if (! allEqual(sgn, ignore)) {
                return 0;
            }
 
            return sgn.front();
        }
 
        bool allEqual(const std::vector<int>& coll,
                      const int               ignore) const
        {
            return this->allEqual(coll, ignore, ignore);
        }
 
        bool allEqual(const std::vector<int>& coll,
                      const int               ignore,
                      const int               lookfor) const
        {
            const auto x0 = (lookfor != ignore)
                ? lookfor : first(coll, ignore);
 
            return std::all_of(std::begin(coll), std::end(coll),
                               [x0, ignore](const int xi)
                   {
                       return (xi == x0) || (xi == ignore);
                   });
        }
 
        int first(const std::vector<int>& coll,
                  const int               ignore) const
        {
            auto e = std::end(coll);
 
            auto p = std::find_if(std::begin(coll), e,
                                  [ignore](const int xi)
                                  {
                                      return xi != ignore;
                                  });
 
            if (p == e) {
                return ignore;
            }
 
            return *p;
        }
    };
 
}} // namespace Opm::UgGridHelpers
 
#endif // OPM_REPAIRZCORN_HEADER_INCLUDED