My Project
Loading...
Searching...
No Matches
RepairZCORN.hpp
Go to the documentation of this file.
1/*
2 Copyright 2017 SINTEF ICT, Applied Mathematics.
3 Copyright 2017 Statoil ASA.
4
5 This file is part of the Open Porous Media Project (OPM).
6
7 OPM is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 OPM is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with OPM. If not, see <http://www.gnu.org/licenses/>.
19*/
20
21#ifndef OPM_REPAIRZCORN_HEADER_INCLUDED
22#define OPM_REPAIRZCORN_HEADER_INCLUDED
23
24#include <algorithm>
25#include <array>
26#include <cassert>
27#include <cmath>
28#include <cstddef>
29#include <exception>
30#include <iterator>
31#include <numeric>
32#include <stdexcept>
33#include <type_traits>
34#include <utility>
35#include <vector>
36
57
58namespace Opm { namespace UgGridHelpers {
59
62 class RepairZCORN
63 {
64 public:
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))
90 {
91 this->ensureZCornIsDepth();
92 this->ensureTopNotBelowBottom();
93 this->ensureBottomNotBelowLowerTop(cartDims[2]);
94 }
95
97 struct ZCornChangeCount
98 {
100 std::size_t cells{0};
101
103 std::size_t corners{0};
104 };
105
110 std::vector<double> destructivelyGrabSanitizedValues()
111 {
112 return std::move(this->zcorn_);
113 }
114
121 bool switchedToDepth() const
122 {
123 return this->switchedToDepth_;
124 }
125
129 const ZCornChangeCount& statTopBelowBottom() const
130 {
131 return this->topBelowBottom_;
132 }
133
137 const ZCornChangeCount& statBottomBelowLowerTop() const
138 {
139 return this->bottomBelowLowerTop_;
140 }
141
142 private:
148 class ActiveCells
149 {
150 public:
168 template <class CartDims>
169 ActiveCells(const std::vector<int>& actnum,
170 const CartDims& cartDims)
171 : nx_(cartDims[0])
172 , ny_(cartDims[1])
173 , nz_(cartDims[2])
174 {
175 const auto nglob = nx_ * ny_ * nz_;
176
177 if (actnum.empty()) {
178 this->is_active_.resize(nglob, true);
179 this->acells_ .resize(nglob, 0);
180
181 std::iota(std::begin(this->acells_),
182 std::end (this->acells_), 0);
183 }
184 else if (actnum.size() == nglob) {
185 this->is_active_.resize(nglob, false);
186 this->acells_ .reserve(nglob);
187
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;
192 }
193 }
194 }
195 else {
196 throw std::invalid_argument {
197 "ACTNUM vector does not match global size"
198 };
199 }
200 }
201
203 using IndexTuple = std::array<std::size_t, 3>;
204
211 IndexTuple getCellIJK(const std::size_t globCell) const
212 {
213 auto c = globCell;
214
215 auto i = c % nx_; c /= nx_;
216 auto j = c % ny_; c /= ny_;
217 auto k = c % nz_;
218
219 return {{ i, j, k }};
220 }
221
223 const std::vector<std::size_t>& activeGlobal() const
224 {
225 return this->acells_;
226 }
227
229 std::size_t numGlobalCells() const
230 {
231 return this->nx_ * this->ny_ * this->nz_;
232 }
233
246 std::size_t neighbourBelow(IndexTuple ijk) const
247 {
248 if (ijk[2] >= nz_ - 1) {
249 return -1;
250 }
251
252 ijk[2] += 1;
253
254 auto below = this->globalCellIdx(ijk);
255 while ((below < this->numGlobalCells()) &&
256 (! this->is_active_[below]))
257 {
258 ijk[2] += 1;
259
260 below = this->globalCellIdx(ijk);
261 }
262
263 return below;
264 }
265
266 private:
268 const std::size_t nx_;
269
271 const std::size_t ny_;
272
274 const std::size_t nz_;
275
279 std::vector<bool> is_active_;
280
282 std::vector<std::size_t> acells_;
283
292 std::size_t globalCellIdx(const IndexTuple& ijk) const
293 {
294 if (ijk[2] > nz_ - 1) { return -1; }
295
296 return ijk[0] + nx_*(ijk[1] + ny_*ijk[2]);
297 }
298 };
299
302 class ZCornIndex
303 {
304 public:
306 template <class CartDims>
307 explicit ZCornIndex(const CartDims& cartDims)
308 : nx_ (cartDims[0])
309 , ny_ (cartDims[1])
310 , nz_ (cartDims[2])
311 , layerOffset_((2 * nx_) * (2 * ny_))
312 {}
313
316 struct PillarPointIDX
317 {
319 std::size_t top;
320
322 std::size_t bottom;
323 };
324
331 template <class IndexTuple>
332 std::array<PillarPointIDX, 4>
333 pillarPoints(const IndexTuple& ijk) const
334 {
335 const auto start = this->getStartIDX(ijk);
336
337 return {
338 this->p00(start),
339 this->p10(start),
340 this->p01(start),
341 this->p11(start)
342 };
343 }
344
345 private:
347 const std::size_t nx_;
348
350 const std::size_t ny_;
351
353 const std::size_t nz_;
354
356 const std::size_t layerOffset_;
357
363 template <class IndexTuple>
364 std::size_t getStartIDX(const IndexTuple& ijk) const
365 {
366 return 2*ijk[0] + 2*nx_*(2*ijk[1] + 2*ny_*(2 * ijk[2]));
367 }
368
376 PillarPointIDX p00(const std::size_t start) const
377 {
378 return this->pillarpts(start, this->offset(0, 0));
379 }
380
388 PillarPointIDX p10(const std::size_t start) const
389 {
390 return this->pillarpts(start, this->offset(1, 0));
391 }
392
400 PillarPointIDX p01(const std::size_t start) const
401 {
402 return this->pillarpts(start, this->offset(0, 1));
403 }
404
412 PillarPointIDX p11(const std::size_t start) const
413 {
414 return this->pillarpts(start, this->offset(1, 1));
415 }
416
426 std::size_t offset(const std::size_t i, const std::size_t j) const
427 {
428 assert ((i == 0) || (i == 1));
429 assert ((j == 0) || (j == 1));
430
431 return i + j*2*nx_;
432 }
433
441 PillarPointIDX pillarpts(const std::size_t start,
442 const std::size_t off) const
443 {
444 return {
445 start + off,
446 start + off + this->layerOffset_
447 };
448 }
449 };
450
452 const ActiveCells active_;
453
455 const ZCornIndex zcorn_idx_;
456
458 std::vector<double> zcorn_;
459
462 bool switchedToDepth_{false};
463
465 ZCornChangeCount topBelowBottom_;
466
468 ZCornChangeCount bottomBelowLowerTop_;
469
473 void ensureZCornIsDepth()
474 {
475 if (this->zcornIsElevation()) {
476 std::transform(this->zcorn_.begin(), this->zcorn_.end(),
477 this->zcorn_.begin(), std::negate<>{});
478
479 this->switchedToDepth_ = true;
480 }
481 }
482
487 void ensureTopNotBelowBottom()
488 {
489 for (const auto& globCell : this->active_.activeGlobal()) {
490 this->ensureTopNotBelowBottom(globCell);
491 }
492 }
493
501 void ensureBottomNotBelowLowerTop(const std::size_t nz)
502 {
503 if (nz == 0) { return; }
504
505 auto bottomLayer = [nz](const std::size_t layerID)
506 {
507 return layerID == (nz - 1);
508 };
509
510 for (const auto& globCell : this->active_.activeGlobal()) {
511 const auto& ijk = this->active_.getCellIJK(globCell);
512
513 if (bottomLayer(ijk[2])) { continue; }
514
515 this->ensureCellBottomNotBelowLowerTop(ijk);
516 }
517 }
518
527 template <class CellIndex>
528 void ensureTopNotBelowBottom(const CellIndex globCell)
529 {
530 const auto cornerCnt0 = this->topBelowBottom_.corners;
531
532 const auto ijk = this->active_.getCellIJK(globCell);
533
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];
537
538 if (zt > zb) { // Top below bottom (ZCORN is depth)
539 zt = zb;
540
541 this->topBelowBottom_.corners += 1;
542 }
543 }
544
545 this->topBelowBottom_.cells +=
546 this->topBelowBottom_.corners > cornerCnt0;
547 }
548
559 template <class IndexTuple>
560 void ensureCellBottomNotBelowLowerTop(const IndexTuple& ijk)
561 {
562 const auto below = this->active_.neighbourBelow(ijk);
563
564 if (below >= this->active_.numGlobalCells()) {
565 return;
566 }
567
568 const auto cornerCnt0 = this->bottomBelowLowerTop_.corners;
569
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);
573
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];
577
578 if (zbu > ztd) { // Bottom below lower top (ZCORN is depth)
579 ztd = zbu;
580
581 this->bottomBelowLowerTop_.corners += 1;
582 }
583 }
584
585 this->bottomBelowLowerTop_.cells +=
586 this->bottomBelowLowerTop_.corners > cornerCnt0;
587 }
588
592 bool zcornIsElevation() const
593 {
594 auto all_signs = std::vector<int>{};
595 all_signs.reserve(this->active_.numGlobalCells());
596
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); });
602
603 // Ignore twisted cells (i.e., cells of indeterminate signs).
604 const int ignore = 0;
605
606 // Elevation implies that ZCORN values are decreasing which
607 // means that the signs in all non-twisted cells equal -1.
608 //
609 // Note: This check is *NOT* equivalent to
610 //
611 // allEqual(all_signs, ignore, -1)
612 //
613 // because that check returns \c true if ALL cells are ignored
614 // whereas first()==-1 ensures that there is at least ONE cell
615 // of determinate sign.
616 return (first(all_signs, ignore) == -1)
617 && allEqual(all_signs, ignore, -1);
618 }
619
632 template <typename CellIndex>
633 int getZCornSign(const CellIndex globCell) const
634 {
635 auto sign = [](const double x) -> int
636 {
637 return (x > 0.0) - (x < 0.0);
638 };
639
640 const auto ijk = this->active_.getCellIJK(globCell);
641
642 auto sgn = std::vector<int>{}; sgn.reserve(4);
643
644 for (const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
645 const auto dz =
646 this->zcorn_[pt.bottom] - this->zcorn_[pt.top];
647
648 sgn.push_back(sign(dz));
649 }
650
651 const int ignore = 0;
652
653 if (! allEqual(sgn, ignore)) {
654 return 0;
655 }
656
657 return sgn.front();
658 }
659
670 bool allEqual(const std::vector<int>& coll,
671 const int ignore) const
672 {
673 return this->allEqual(coll, ignore, ignore);
674 }
675
697 bool allEqual(const std::vector<int>& coll,
698 const int ignore,
699 const int lookfor) const
700 {
701 const auto x0 = (lookfor != ignore)
702 ? lookfor : first(coll, ignore);
703
704 return std::all_of(std::begin(coll), std::end(coll),
705 [x0, ignore](const int xi)
706 {
707 return (xi == x0) || (xi == ignore);
708 });
709 }
710
721 int first(const std::vector<int>& coll,
722 const int ignore) const
723 {
724 auto e = std::end(coll);
725
726 auto p = std::find_if(std::begin(coll), e,
727 [ignore](const int xi)
728 {
729 return xi != ignore;
730 });
731
732 if (p == e) {
733 return ignore;
734 }
735
736 return *p;
737 }
738 };
739
740}} // namespace Opm::UgGridHelpers
741
742#endif // OPM_REPAIRZCORN_HEADER_INCLUDED
Holds the implementation of the CpGrid as a pimple.
Definition CellQuadrature.cpp:68