21#ifndef OPM_REPAIRZCORN_HEADER_INCLUDED
22#define OPM_REPAIRZCORN_HEADER_INCLUDED
58namespace Opm {
namespace UgGridHelpers {
83 template <
class CartDims>
84 RepairZCORN(std::vector<double>&& zcorn,
85 const std::vector<int>& actnum,
86 const CartDims& cartDims)
87 : active_ (actnum, cartDims)
88 , zcorn_idx_(cartDims)
89 , zcorn_ (std::move(zcorn))
91 this->ensureZCornIsDepth();
92 this->ensureTopNotBelowBottom();
93 this->ensureBottomNotBelowLowerTop(cartDims[2]);
97 struct ZCornChangeCount
100 std::size_t cells{0};
103 std::size_t corners{0};
110 std::vector<double> destructivelyGrabSanitizedValues()
112 return std::move(this->zcorn_);
121 bool switchedToDepth()
const
123 return this->switchedToDepth_;
129 const ZCornChangeCount& statTopBelowBottom()
const
131 return this->topBelowBottom_;
137 const ZCornChangeCount& statBottomBelowLowerTop()
const
139 return this->bottomBelowLowerTop_;
168 template <
class CartDims>
169 ActiveCells(
const std::vector<int>& actnum,
170 const CartDims& cartDims)
175 const auto nglob = nx_ * ny_ * nz_;
177 if (actnum.empty()) {
178 this->is_active_.resize(nglob,
true);
179 this->acells_ .resize(nglob, 0);
181 std::iota(std::begin(this->acells_),
182 std::end (this->acells_), 0);
184 else if (actnum.size() == nglob) {
185 this->is_active_.resize(nglob,
false);
186 this->acells_ .reserve(nglob);
188 for (
auto i = 0*nglob; i < nglob; ++i) {
189 if (actnum[i] != 0) {
190 this->acells_.push_back(i);
191 this->is_active_[i] =
true;
196 throw std::invalid_argument {
197 "ACTNUM vector does not match global size"
203 using IndexTuple = std::array<std::size_t, 3>;
211 IndexTuple getCellIJK(
const std::size_t globCell)
const
215 auto i = c % nx_; c /= nx_;
216 auto j = c % ny_; c /= ny_;
219 return {{ i, j, k }};
223 const std::vector<std::size_t>& activeGlobal()
const
225 return this->acells_;
229 std::size_t numGlobalCells()
const
231 return this->nx_ * this->ny_ * this->nz_;
246 std::size_t neighbourBelow(IndexTuple ijk)
const
248 if (ijk[2] >= nz_ - 1) {
254 auto below = this->globalCellIdx(ijk);
255 while ((below < this->numGlobalCells()) &&
256 (! this->is_active_[below]))
260 below = this->globalCellIdx(ijk);
268 const std::size_t nx_;
271 const std::size_t ny_;
274 const std::size_t nz_;
279 std::vector<bool> is_active_;
282 std::vector<std::size_t> acells_;
292 std::size_t globalCellIdx(
const IndexTuple& ijk)
const
294 if (ijk[2] > nz_ - 1) {
return -1; }
296 return ijk[0] + nx_*(ijk[1] + ny_*ijk[2]);
306 template <
class CartDims>
307 explicit ZCornIndex(
const CartDims& cartDims)
311 , layerOffset_((2 * nx_) * (2 * ny_))
316 struct PillarPointIDX
331 template <
class IndexTuple>
332 std::array<PillarPointIDX, 4>
333 pillarPoints(
const IndexTuple& ijk)
const
335 const auto start = this->getStartIDX(ijk);
347 const std::size_t nx_;
350 const std::size_t ny_;
353 const std::size_t nz_;
356 const std::size_t layerOffset_;
363 template <
class IndexTuple>
364 std::size_t getStartIDX(
const IndexTuple& ijk)
const
366 return 2*ijk[0] + 2*nx_*(2*ijk[1] + 2*ny_*(2 * ijk[2]));
376 PillarPointIDX p00(
const std::size_t start)
const
378 return this->pillarpts(start, this->offset(0, 0));
388 PillarPointIDX p10(
const std::size_t start)
const
390 return this->pillarpts(start, this->offset(1, 0));
400 PillarPointIDX p01(
const std::size_t start)
const
402 return this->pillarpts(start, this->offset(0, 1));
412 PillarPointIDX p11(
const std::size_t start)
const
414 return this->pillarpts(start, this->offset(1, 1));
426 std::size_t offset(
const std::size_t i,
const std::size_t j)
const
428 assert ((i == 0) || (i == 1));
429 assert ((j == 0) || (j == 1));
441 PillarPointIDX pillarpts(
const std::size_t start,
442 const std::size_t off)
const
446 start + off + this->layerOffset_
452 const ActiveCells active_;
455 const ZCornIndex zcorn_idx_;
458 std::vector<double> zcorn_;
462 bool switchedToDepth_{
false};
465 ZCornChangeCount topBelowBottom_;
468 ZCornChangeCount bottomBelowLowerTop_;
473 void ensureZCornIsDepth()
475 if (this->zcornIsElevation()) {
476 std::transform(this->zcorn_.begin(), this->zcorn_.end(),
477 this->zcorn_.begin(), std::negate<>{});
479 this->switchedToDepth_ =
true;
487 void ensureTopNotBelowBottom()
489 for (
const auto& globCell : this->active_.activeGlobal()) {
490 this->ensureTopNotBelowBottom(globCell);
501 void ensureBottomNotBelowLowerTop(
const std::size_t nz)
503 if (nz == 0) {
return; }
505 auto bottomLayer = [nz](
const std::size_t layerID)
507 return layerID == (nz - 1);
510 for (
const auto& globCell : this->active_.activeGlobal()) {
511 const auto& ijk = this->active_.getCellIJK(globCell);
513 if (bottomLayer(ijk[2])) {
continue; }
515 this->ensureCellBottomNotBelowLowerTop(ijk);
527 template <
class CellIndex>
528 void ensureTopNotBelowBottom(
const CellIndex globCell)
530 const auto cornerCnt0 = this->topBelowBottom_.corners;
532 const auto ijk = this->active_.getCellIJK(globCell);
534 for (
const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
535 const auto zb = this->zcorn_[pt.bottom];
536 auto& zt = this->zcorn_[pt.top];
541 this->topBelowBottom_.corners += 1;
545 this->topBelowBottom_.cells +=
546 this->topBelowBottom_.corners > cornerCnt0;
559 template <
class IndexTuple>
560 void ensureCellBottomNotBelowLowerTop(
const IndexTuple& ijk)
562 const auto below = this->active_.neighbourBelow(ijk);
564 if (below >= this->active_.numGlobalCells()) {
568 const auto cornerCnt0 = this->bottomBelowLowerTop_.corners;
570 const auto& up = this->zcorn_idx_.pillarPoints(ijk);
571 const auto d_ijk = this->active_.getCellIJK(below);
572 const auto& down = this->zcorn_idx_.pillarPoints(d_ijk);
574 for (
auto n = up.size(), i = 0*n; i < n; ++i) {
575 const auto zbu = this->zcorn_[up [i].bottom];
576 auto& ztd = this->zcorn_[down[i].top];
581 this->bottomBelowLowerTop_.corners += 1;
585 this->bottomBelowLowerTop_.cells +=
586 this->bottomBelowLowerTop_.corners > cornerCnt0;
592 bool zcornIsElevation()
const
594 auto all_signs = std::vector<int>{};
595 all_signs.reserve(this->active_.numGlobalCells());
597 std::transform(this->active_.activeGlobal().begin(),
598 this->active_.activeGlobal().end(),
599 std::back_inserter(all_signs),
600 [
this](
const auto globCell)
601 { return this->getZCornSign(globCell); });
604 const int ignore = 0;
616 return (first(all_signs, ignore) == -1)
617 && allEqual(all_signs, ignore, -1);
632 template <
typename CellIndex>
633 int getZCornSign(
const CellIndex globCell)
const
635 auto sign = [](
const double x) ->
int
637 return (x > 0.0) - (x < 0.0);
640 const auto ijk = this->active_.getCellIJK(globCell);
642 auto sgn = std::vector<int>{}; sgn.reserve(4);
644 for (
const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
646 this->zcorn_[pt.bottom] - this->zcorn_[pt.top];
648 sgn.push_back(sign(dz));
651 const int ignore = 0;
653 if (! allEqual(sgn, ignore)) {
670 bool allEqual(
const std::vector<int>& coll,
671 const int ignore)
const
673 return this->allEqual(coll, ignore, ignore);
697 bool allEqual(
const std::vector<int>& coll,
699 const int lookfor)
const
701 const auto x0 = (lookfor != ignore)
702 ? lookfor : first(coll, ignore);
704 return std::all_of(std::begin(coll), std::end(coll),
705 [x0, ignore](
const int xi)
707 return (xi == x0) || (xi == ignore);
721 int first(
const std::vector<int>& coll,
722 const int ignore)
const
724 auto e = std::end(coll);
726 auto p = std::find_if(std::begin(coll), e,
727 [ignore](
const int xi)
Holds the implementation of the CpGrid as a pimple.
Definition CellQuadrature.cpp:68