All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
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 
58 namespace 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  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  for (auto& zc : this->zcorn_) {
477  zc = - zc;
478  }
479 
480  this->switchedToDepth_ = true;
481  }
482  }
483 
488  void ensureTopNotBelowBottom()
489  {
490  for (const auto& globCell : this->active_.activeGlobal()) {
491  this->ensureTopNotBelowBottom(globCell);
492  }
493  }
494 
502  void ensureBottomNotBelowLowerTop(const std::size_t nz)
503  {
504  if (nz == 0) { return; }
505 
506  auto bottomLayer = [nz](const std::size_t layerID)
507  {
508  return layerID == (nz - 1);
509  };
510 
511  for (const auto& globCell : this->active_.activeGlobal()) {
512  const auto& ijk = this->active_.getCellIJK(globCell);
513 
514  if (bottomLayer(ijk[2])) { continue; }
515 
516  this->ensureCellBottomNotBelowLowerTop(ijk);
517  }
518  }
519 
528  template <class CellIndex>
529  void ensureTopNotBelowBottom(const CellIndex globCell)
530  {
531  const auto cornerCnt0 = this->topBelowBottom_.corners;
532 
533  const auto ijk = this->active_.getCellIJK(globCell);
534 
535  for (const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
536  const auto zb = this->zcorn_[pt.bottom];
537  auto& zt = this->zcorn_[pt.top];
538 
539  if (zt > zb) { // Top below bottom (ZCORN is depth)
540  zt = zb;
541 
542  this->topBelowBottom_.corners += 1;
543  }
544  }
545 
546  this->topBelowBottom_.cells +=
547  this->topBelowBottom_.corners > cornerCnt0;
548  }
549 
560  template <class IndexTuple>
561  void ensureCellBottomNotBelowLowerTop(const IndexTuple& ijk)
562  {
563  const auto below = this->active_.neighbourBelow(ijk);
564 
565  if (below >= this->active_.numGlobalCells()) {
566  return;
567  }
568 
569  const auto cornerCnt0 = this->bottomBelowLowerTop_.corners;
570 
571  const auto& up = this->zcorn_idx_.pillarPoints(ijk);
572  const auto d_ijk = this->active_.getCellIJK(below);
573  const auto& down = this->zcorn_idx_.pillarPoints(d_ijk);
574 
575  for (auto n = up.size(), i = 0*n; i < n; ++i) {
576  const auto zbu = this->zcorn_[up [i].bottom];
577  auto& ztd = this->zcorn_[down[i].top];
578 
579  if (zbu > ztd) { // Bottom below lower top (ZCORN is depth)
580  ztd = zbu;
581 
582  this->bottomBelowLowerTop_.corners += 1;
583  }
584  }
585 
586  this->bottomBelowLowerTop_.cells +=
587  this->bottomBelowLowerTop_.corners > cornerCnt0;
588  }
589 
593  bool zcornIsElevation() const
594  {
595  auto all_signs = std::vector<int>{};
596  all_signs.reserve(this->active_.numGlobalCells());
597 
598  for (const auto& globCell : this->active_.activeGlobal()) {
599  all_signs.push_back(this->getZCornSign(globCell));
600  }
601 
602  // Ignore twisted cells (i.e., cells of indeterminate signs).
603  const int ignore = 0;
604 
605  // Elevation implies that ZCORN values are decreasing which
606  // means that the signs in all non-twisted cells equal -1.
607  //
608  // Note: This check is *NOT* equivalent to
609  //
610  // allEqual(all_signs, ignore, -1)
611  //
612  // because that check returns \c true if ALL cells are ignored
613  // whereas first()==-1 ensures that there is at least ONE cell
614  // of determinate sign.
615  return (first(all_signs, ignore) == -1)
616  && allEqual(all_signs, ignore, -1);
617  }
618 
631  template <typename CellIndex>
632  int getZCornSign(const CellIndex globCell) const
633  {
634  auto sign = [](const double x) -> int
635  {
636  return (x > 0.0) - (x < 0.0);
637  };
638 
639  const auto ijk = this->active_.getCellIJK(globCell);
640 
641  auto sgn = std::vector<int>{}; sgn.reserve(4);
642 
643  for (const auto& pt : this->zcorn_idx_.pillarPoints(ijk)) {
644  const auto dz =
645  this->zcorn_[pt.bottom] - this->zcorn_[pt.top];
646 
647  sgn.push_back(sign(dz));
648  }
649 
650  const int ignore = 0;
651 
652  if (! allEqual(sgn, ignore)) {
653  return 0;
654  }
655 
656  return sgn.front();
657  }
658 
669  bool allEqual(const std::vector<int>& coll,
670  const int ignore) const
671  {
672  return this->allEqual(coll, ignore, ignore);
673  }
674 
696  bool allEqual(const std::vector<int>& coll,
697  const int ignore,
698  const int lookfor) const
699  {
700  const auto x0 = (lookfor != ignore)
701  ? lookfor : first(coll, ignore);
702 
703  return std::all_of(std::begin(coll), std::end(coll),
704  [x0, ignore](const int xi)
705  {
706  return (xi == x0) || (xi == ignore);
707  });
708  }
709 
720  int first(const std::vector<int>& coll,
721  const int ignore) const
722  {
723  auto e = std::end(coll);
724 
725  auto p = std::find_if(std::begin(coll), e,
726  [ignore](const int xi)
727  {
728  return xi != ignore;
729  });
730 
731  if (p == e) {
732  return ignore;
733  }
734 
735  return *p;
736  }
737  };
738 
739 }} // namespace Opm::UgGridHelpers
740 
741 #endif // OPM_REPAIRZCORN_HEADER_INCLUDED