Loading...
Searching...
No Matches
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>
explicit 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()) {
std::transform(this->zcorn_.begin(), this->zcorn_.end(),
this->zcorn_.begin(), std::negate<>{});
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());
std::transform(this->active_.activeGlobal().begin(),
this->active_.activeGlobal().end(),
std::back_inserter(all_signs),
[this](const auto globCell)
{ return 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
Opm
Holds the implementation of the CpGrid as a pimple.
Definition CellQuadrature.cpp:68