Skip to content

Commit 6651530

Browse files
author
Lisa Julia Nebel
authored
Merge pull request #836 from aritorto/loopElem
[small refactor] Lgr ZCORN computation
2 parents c0d2beb + c030941 commit 6651530

File tree

1 file changed

+42
-54
lines changed

1 file changed

+42
-54
lines changed

‎opm/grid/cpgrid/CpGridUtilities.cpp‎

Lines changed: 42 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -96,9 +96,6 @@ lgrCOORDandZCORN(const Dune::CpGrid& grid,
9696
// Initialize all ZCORN as inactive (setting values to std::numeric_limits<double>::max()).
9797
std::vector<double> lgrZCORN(8*nx*ny*nz, std::numeric_limits<double>::max());
9898

99-
// Store cells per layer (via their ijk's, i.e. cells_per_layer[ k ] = { {i0, j0}, {i1, j1}, ... }).
100-
std::vector<std::vector<std::pair<int,int>>> cells_per_layer(nz);
101-
10299
// Map to determine min and max k per cell column (i, j) (min/max_k = 0, ..., nz-1).
103100
// Initialized as {nz, -1} to detect inactive cell columns.
104101
std::vector<std::array<int,2>> minMaxPerCellPillar(nx*ny, {nz, -1});
@@ -112,61 +109,52 @@ lgrCOORDandZCORN(const Dune::CpGrid& grid,
112109

113110
minMax[0] = std::min(ijk[2], minMax[0]);
114111
minMax[1] = std::max(ijk[2], minMax[1]);
115-
116-
// Store cells per layer
117-
cells_per_layer[ijk[2]].emplace_back(std::make_pair(ijk[0],ijk[1]));
118112
}
119113

120-
for (int layer = 0; layer < nz; ++layer) {
121-
for (const auto& ij : cells_per_layer[layer]) {
122-
int cell_lgr_cartesian_idx = (layer*nx*ny) + (ij.second *nx) + ij.first;
123-
124-
const auto& elemIdx = lgrCartesianIdxToCellIdx.at(cell_lgr_cartesian_idx);
125-
const auto& elem = Dune::cpgrid::Entity<0>(levelGrid, elemIdx, true);
126-
127-
// For a grid with nz layers, ZCORN values are ordered:
128-
//
129-
// top layer nz-1
130-
// bottom layer nz-1
131-
// top layer nz-2
132-
// bottom layer nz-2
133-
// ...
134-
// top layer 1
135-
// bottom layer 1
136-
// top layer 0
137-
// bottom layer 0
138-
139-
int zcorn_top_00_idx = ((nz-1-layer)*8*nx*ny) + (ij.second*4*nx) + (2*ij.first); // assoc. w. elem corner 4
140-
141-
// Bottom indices
142-
int zcorn_top_10_idx = zcorn_top_00_idx + 1; // assoc. w. elem corner 5
143-
int zcorn_top_01_idx = zcorn_top_00_idx + (2*nx); // assoc. w. elem corner 6
144-
int zcorn_top_11_idx = zcorn_top_01_idx + 1; // assoc. w. elem corner 7
145-
146-
// Top indices
147-
int zcorn_bottom_00_idx = zcorn_top_00_idx + (4*nx*ny); // assoc. w. elem corner 0
148-
int zcorn_bottom_10_idx = zcorn_bottom_00_idx + 1; // assoc. w. elem corner 1
149-
int zcorn_bottom_01_idx = zcorn_bottom_00_idx + (2*nx); // assoc. w. elem corner 2
150-
int zcorn_bottom_11_idx = zcorn_bottom_01_idx + 1; // assoc. w. elem corner
151-
152-
// Note: zcorn_idx + 1 moves to the next position along the x-axis (i+1, j, k)
153-
// zcorn_idx + (2*nx) moves to the next position along the y-axis (i, j+1, k)
154-
// zcorn_idx + (4*nx*ny) moves to the next position along the z-axis (i,j, k+1)
155-
156-
// Assign ZCORN values
157-
lgrZCORN[zcorn_top_00_idx] = elem.subEntity<3>(4).geometry().center()[2];
158-
lgrZCORN[zcorn_top_10_idx] = elem.subEntity<3>(5).geometry().center()[2];
159-
lgrZCORN[zcorn_top_01_idx] = elem.subEntity<3>(6).geometry().center()[2];
160-
lgrZCORN[zcorn_top_11_idx] = elem.subEntity<3>(7).geometry().center()[2];
161-
162-
lgrZCORN[zcorn_bottom_00_idx] = elem.subEntity<3>(0).geometry().center()[2];
163-
lgrZCORN[zcorn_bottom_10_idx] = elem.subEntity<3>(1).geometry().center()[2];
164-
lgrZCORN[zcorn_bottom_01_idx] = elem.subEntity<3>(2).geometry().center()[2];
165-
lgrZCORN[zcorn_bottom_11_idx] = elem.subEntity<3>(3).geometry().center()[2];
166-
}
114+
for (const auto& elem : elements(grid.levelGridView(level))) {
115+
const auto& elemIJK = lgrIJK[elem.index()];
116+
117+
// For a grid with nz layers, ZCORN values are ordered:
118+
//
119+
// top layer nz-1
120+
// bottom layer nz-1
121+
// top layer nz-2
122+
// bottom layer nz-2
123+
// ...
124+
// top layer 1
125+
// bottom layer 1
126+
// top layer 0
127+
// bottom layer 0
128+
129+
int zcorn_top_00_idx = ((nz-1-elemIJK[2])*8*nx*ny) + (elemIJK[1]*4*nx) + (2*elemIJK[0]); // assoc. w. elem corner 4
130+
131+
// Bottom indices
132+
int zcorn_top_10_idx = zcorn_top_00_idx + 1; // assoc. w. elem corner 5
133+
int zcorn_top_01_idx = zcorn_top_00_idx + (2*nx); // assoc. w. elem corner 6
134+
int zcorn_top_11_idx = zcorn_top_01_idx + 1; // assoc. w. elem corner 7
135+
136+
// Top indices
137+
int zcorn_bottom_00_idx = zcorn_top_00_idx + (4*nx*ny); // assoc. w. elem corner 0
138+
int zcorn_bottom_10_idx = zcorn_bottom_00_idx + 1; // assoc. w. elem corner 1
139+
int zcorn_bottom_01_idx = zcorn_bottom_00_idx + (2*nx); // assoc. w. elem corner 2
140+
int zcorn_bottom_11_idx = zcorn_bottom_01_idx + 1; // assoc. w. elem corner
141+
142+
// Note: zcorn_idx + 1 moves to the next position along the x-axis (i+1, j, k)
143+
// zcorn_idx + (2*nx) moves to the next position along the y-axis (i, j+1, k)
144+
// zcorn_idx + (4*nx*ny) moves to the next position along the z-axis (i,j, k+1)
145+
146+
// Assign ZCORN values
147+
lgrZCORN[zcorn_top_00_idx] = elem.subEntity<3>(4).geometry().center()[2];
148+
lgrZCORN[zcorn_top_10_idx] = elem.subEntity<3>(5).geometry().center()[2];
149+
lgrZCORN[zcorn_top_01_idx] = elem.subEntity<3>(6).geometry().center()[2];
150+
lgrZCORN[zcorn_top_11_idx] = elem.subEntity<3>(7).geometry().center()[2];
151+
152+
lgrZCORN[zcorn_bottom_00_idx] = elem.subEntity<3>(0).geometry().center()[2];
153+
lgrZCORN[zcorn_bottom_10_idx] = elem.subEntity<3>(1).geometry().center()[2];
154+
lgrZCORN[zcorn_bottom_01_idx] = elem.subEntity<3>(2).geometry().center()[2];
155+
lgrZCORN[zcorn_bottom_11_idx] = elem.subEntity<3>(3).geometry().center()[2];
167156
}
168157

169-
170158
// Rewrite values for active pillars
171159
for (int j = 0; j < ny; ++j) {
172160
for (int i = 0; i < nx; ++i) {

0 commit comments

Comments
 (0)