Grid Utility Functions for Unstructured Grid Conversion into CPG format

CMakeLists_files.cmake

@@ -46,6 +46,7 @@ list (APPEND MAIN_SOURCE_FILES
       opm/common/utility/numeric/calculateCellVol.cpp
       opm/common/utility/numeric/GeometryUtil.cpp
       opm/common/utility/numeric/VectorUtil.cpp
+      opm/common/utility/numeric/GridUtil.cpp
       opm/common/utility/numeric/MonotCubicInterpolator.cpp
       opm/common/utility/numeric/RootFinders.cpp
       opm/material/common/Spline.cpp
@@ -477,6 +478,7 @@ list (APPEND TEST_SOURCE_FILES
       tests/test_sparsevector.cpp
       tests/test_ThreadSafeMapBuilder.cpp
       tests/test_uniformtablelinear.cpp
+      tests/test_Uns2CPG.cpp
       tests/material/test_2dtables.cpp
       tests/material/test_blackoilfluidstate.cpp
       tests/material/test_components.cpp
@@ -851,6 +853,7 @@ list( APPEND PUBLIC_HEADER_FILES
       opm/common/utility/numeric/calculateCellVol.hpp
       opm/common/utility/numeric/GeometryUtil.hpp
       opm/common/utility/numeric/VectorUtil.hpp
+      opm/common/utility/numeric/GridUtil.hpp
       opm/common/utility/numeric/buildUniformMonotoneTable.hpp
       opm/common/utility/numeric/linearInterpolation.hpp
       opm/common/utility/numeric/MonotCubicInterpolator.hpp

opm/common/utility/numeric/GridUtil.cpp

@@ -0,0 +1,96 @@
+#include <opm/common/utility/numeric/GridUtil.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/GridDims.hpp>
+#include <tuple>
+namespace Opm::GridUtil {
+
+
+static std::vector<double> pillar_to_flat_array(const std::vector<std::array<std::array<double, 3>, 2>>& pillar)
+{
+    std::vector<double> flat_array;
+    // Pre-allocate space for efficiency: 2 points per pillar entry, 3 components each
+    flat_array.reserve(pillar.size() * 2 * 3);
+    // Iterate through each pair of points in the pillar
+    for (const auto& pair : pillar) {
+        // Insert all components of the first point (x, y, z)
+        flat_array.insert(flat_array.end(), pair[0].begin(), pair[0].end());
+        // Insert all components of the second point (x, y, z)
+        flat_array.insert(flat_array.end(), pair[1].begin(), pair[1].end());
+    }
+    return flat_array;
+}
+
+
+std::tuple<std::vector<double>,std::vector<double>> convertUnsToCPG(
+                     const std::vector<std::array<double, 3>>& coord_uns, 
+                     const std::vector<std::array<std::size_t, 8>>& element, 
+                     std::size_t nx, std::size_t ny, std::size_t nz) 
+{                        
+    // nx, ny, nz are the number of cells in the x, y and z directions
+    // converts unstructured mesh grid described by the coord_uns and element arrays
+    // element contains a referece to the nodes described by coords_uns
+    const std::size_t element_size = element.size();
+    const std::size_t num_pillars = (nx+1)*(ny+1);
+    GridDims cpg_grid = GridDims(nx, ny, nz);
+    GridDims pillar_grid = GridDims(nx+1, ny+1, 0);
+
+
+
+    auto ij_pillars = [ &pillar_grid](std::size_t i, std::size_t j) {
+        return pillar_grid.getGlobalIndex(i, j, 0);
+    };
+
+    auto compute_zcornind = [&nx, &ny]
+            (std::size_t  i, std::size_t j, std::size_t  k) 
+    {
+        std::array<std::size_t, 8> zind;
+        const std::size_t num = k * nx * ny * 8 + j * nx * 4 + i * 2;
+        const std::size_t offset_layer = 4 * nx * ny;
+        zind[0] = num;
+        zind[1] = num + 1;
+        zind[2] = num + 2 * nx;
+        zind[3] = zind[2] + 1;
+        zind[4] = num + offset_layer;
+        zind[5] = zind[1] + offset_layer;
+        zind[6] = zind[2] + offset_layer;
+        zind[7] = zind[3] + offset_layer;
+        return zind;
+    };    
+
+
+    std::vector<std::array<std::array<double, 3>,2>> pillars(num_pillars);
+
+    // looping through the first layer of elements
+    for (std::size_t index = 0; index < nx*ny; index++) {
+       auto [ii,jj,kk] = cpg_grid.getIJK((index));
+       const std::array<std::size_t,8>& element_nodes = element[index];
+       pillars[ij_pillars(ii,jj)][0] = coord_uns[element_nodes[0]];
+       pillars[ij_pillars(ii+1,jj)][0] = coord_uns[element_nodes[1]];
+       pillars[ij_pillars(ii,jj+1)][0] = coord_uns[element_nodes[2]];
+       pillars[ij_pillars(ii+1,jj+1)][0] = coord_uns[element_nodes[3]];
+    }
+    // // looping through the last layer of elements
+    for (std::size_t index = (element_size - nx*ny); index < element_size; index++) {
+       auto [ii,jj,kk] = cpg_grid.getIJK((index));
+       const std::array<std::size_t,8>& element_nodes = element[index];
+       pillars[ij_pillars(ii,jj)][1] = coord_uns[element_nodes[4]];
+       pillars[ij_pillars(ii+1,jj)][1] = coord_uns[element_nodes[5]];
+       pillars[ij_pillars(ii,jj+1)][1] = coord_uns[element_nodes[6]];
+       pillars[ij_pillars(ii+1,jj+1)][1] = coord_uns[element_nodes[7]];
+    }
+    std::vector<double> coord_cpg = pillar_to_flat_array(pillars);
+    std::vector<double> zcorn_cpg(element_size*8);
+    for (std::size_t index = 0; index < element_size; index++) {
+       auto [ii,jj,kk] = cpg_grid.getIJK((index));
+       std::array<std::size_t,8> z_ref = compute_zcornind(ii,jj,kk);       
+       const std::array<std::size_t,8>& element_nodes = element[index];
+       for (std::size_t index_el = 0; index_el < 8; index_el++) {
+           std::size_t local_z = z_ref[index_el];
+           std::size_t local_node = element_nodes[index_el];
+           zcorn_cpg[local_z] = coord_uns[local_node][2];
+       }
+    }
+    return std::make_tuple(coord_cpg,zcorn_cpg);
+}
+
+
+} // namespace GeometryUtils

opm/common/utility/numeric/GridUtil.hpp

@@ -0,0 +1,22 @@
+#ifndef GRIDUTIL_H
+#define GRIDUTIL_H
+
+#include <cmath> 
+#include <tuple>
+#include <vector>
+namespace Opm::GridUtil {
+
+
+std::tuple<std::vector<double>,std::vector<double>> convertUnsToCPG(
+                     const std::vector<std::array<double, 3>>&, 
+                     const std::vector<std::array<std::size_t, 8>>&, 
+                     std::size_t, std::size_t, std::size_t);
+
+
+
+
+// Example for other utilities...
+
+} // namespace GeometryUtils
+
+#endif // GEOMETRY_UTILS_H

tests/test_Uns2CPG.cpp

@@ -0,0 +1,98 @@
+/*
+ Copyright (C) 2023 Equinor
+  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/>.
+ */
+
+#define BOOST_TEST_MODULE TestUnstoCPG
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/test_tools.hpp>
+#include <opm/input/eclipse/Deck/Deck.hpp>
+#include <opm/input/eclipse/Parser/Parser.hpp>
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/LgrCollection.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/Carfin.hpp>
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/common/utility/numeric/GridUtil.hpp>
+
+
+
+using namespace Opm;
+
+bool are_equal(const std::vector<double>& a, const std::vector<double>& b, double epsilon = 1e-8) {
+    return a.size() == b.size() &&
+           std::equal(a.begin(), a.end(), b.begin(),
+                     [epsilon](double x, double y) {
+                         return std::abs(x - y) < epsilon;
+                     });
+}
+
+BOOST_AUTO_TEST_CASE(TestGrid2CPG1) { 
+std::size_t nx, ny, nz;
+nx = 2;
+ny=  1;
+nz = 2;   
+std::vector<std::array<double, 3>> coord = {
+{1.00000000,1.00000000,-0.05000000},
+{1.00000000,0.00000000,-0.05000000},
+{0.45000000,0.00000000,-0.01535455},
+{0.45000000,1.00000000,0.03099800},
+{0.00000000,1.00000000,0.05000000},
+{0.00000000,0.00000000,0.05000000},
+{1.00000000,1.00000000,0.20000000},
+{1.00000000,0.00000000,0.20000000},
+{0.47498716,0.00000000,0.24219037},
+{0.47501669,1.00000000,0.26565554},
+{0.00000000,0.00000000,0.30000000},
+{0.00000000,1.00000000,0.30000000},
+{1.00000000,1.00000000,0.45000000},
+{1.00000000,0.00000000,0.45000000},
+{0.50000000,0.00000000,0.50000000},
+{0.50000000,1.00000000,0.50000000},
+{0.00000000,0.00000000,0.55000000},
+{0.00000000,1.00000000,0.55000000},
+};
+
+std::vector<std::array<std::size_t, 8>> element = {
+{5,2,4,3,10,8,11,9},
+{2,1,3,0,8,7,9,6},
+{10,8,11,9,16,14,17,15},
+{8,7,9,6,14,13,15,12},
+};
+
+std::vector<double> pCOORDS ={
+    0.00000000,0.00000000,0.05000000,0.00000000,0.00000000,
+    0.55000000,0.45000000,0.00000000,-0.01535455,0.50000000,
+    0.00000000,0.50000000,1.00000000,0.00000000,-0.05000000,
+    1.00000000,0.00000000,0.45000000,0.00000000,1.00000000,
+    0.05000000,0.00000000,1.00000000,0.55000000,0.45000000,
+    1.00000000,0.03099800,0.50000000,1.00000000,0.50000000,
+    1.00000000,1.00000000,-0.05000000,1.00000000,1.00000000,
+    0.45000000};
+
+std::vector<double> pZCORN = {
+0.05000000,-0.01535455,-0.01535455,-0.05000000,0.05000000,
+0.03099800,0.03099800,-0.05000000,0.30000000,0.24219037,
+0.24219037,0.20000000,0.30000000,0.26565554,0.26565554,
+0.20000000,0.30000000,0.24219037,0.24219037,0.20000000,
+0.30000000,0.26565554,0.26565554,0.20000000,0.55000000,
+0.50000000,0.50000000,0.45000000,0.55000000,0.50000000,
+0.50000000,0.45000000
+};
+
+auto [COORDS, ZCORN] = 
+                        Opm::GridUtil::convertUnsToCPG(coord,element, nx, ny,nz);
+BOOST_CHECK(are_equal(pCOORDS, COORDS));
+BOOST_CHECK(are_equal(pZCORN, ZCORN));
+
+}