ENH,CLN: add and migrate C++ methods for grid cell computing

src/lib/CMakeLists.txt

@@ -42,6 +42,7 @@ pybind11_add_module(
   "${SRC}/common/geometry/volumes.cpp"
   "${SRC}/grid3d/cell.cpp"
   "${SRC}/grid3d/grid.cpp"
+  "${SRC}/grid3d/grid_surf_oper.cpp"
   "${SRC}/cube/cube.cpp"
   "${SRC}/regsurf/sample_grid3d.cpp"
   "${SRC}/regsurf/utilities.cpp"

src/lib/include/xtgeo/grid3d.hpp

@@ -19,6 +19,16 @@ grid_cell_volumes(const size_t ncol,
                   const py::array_t<int> &actnumsv,
                   const int precision,
                   const bool asmasked);
+
+std::tuple<py::array_t<double>, py::array_t<double>, py::array_t<double>>
+grid_cell_centers(const size_t ncol,
+                  const size_t nrow,
+                  const size_t nlay,
+                  const py::array_t<double> &coordsv,
+                  const py::array_t<float> &zcornsv,
+                  const py::array_t<int> &actnumsv,
+                  const bool asmasked);
+
 std::tuple<py::array_t<double>, py::array_t<double>, py::array_t<double>>
 grid_height_above_ffl(const size_t ncol,
                       const size_t nrow,
@@ -53,6 +63,30 @@ get_depth_in_cell(const double x,
                   const std::array<double, 24> &corners,
                   int option);
 
+py::array_t<int>
+grid_assign_value_between_surfaces(const size_t ncol,
+                                   const size_t nrow,
+                                   const size_t nlay,
+                                   const py::array_t<double> &xmid,
+                                   const py::array_t<double> &ymid,
+                                   const py::array_t<double> &zmid,
+                                   const size_t top_ncol,
+                                   const size_t top_nrow,
+                                   const double top_xori,
+                                   const double top_yori,
+                                   const double top_xinc,
+                                   const double top_yinc,
+                                   const double top_rotation,
+                                   const py::array_t<double> &top_values,
+                                   const size_t bot_ncol,
+                                   const size_t bot_nrow,
+                                   const double bot_xori,
+                                   const double bot_yori,
+                                   const double bot_xinc,
+                                   const double bot_yinc,
+                                   const double bot_rotation,
+                                   const py::array_t<double> &bot_values);
+
 inline void
 init(py::module &m)
 {
@@ -61,6 +95,8 @@ init(py::module &m)
 
     m_grid3d.def("grid_cell_volumes", &grid_cell_volumes,
                  "Compute the bulk volume of cell.");
+    m_grid3d.def("grid_cell_centers", &grid_cell_centers,
+                 "Compute the cells centers coordinates as 3 arrays");
     m_grid3d.def("grid_height_above_ffl", &grid_height_above_ffl,
                  "Compute the height above a FFL (free fluid level).");
     m_grid3d.def("cell_corners", &cell_corners,
@@ -71,6 +107,9 @@ init(py::module &m)
                  "Determine if a XY point is inside a cell, top or base.");
     m_grid3d.def("get_depth_in_cell", &get_depth_in_cell,
                  "Determine the interpolated cell face Z from XY, top or base.");
+    m_grid3d.def("grid_assign_value_between_surfaces",
+                 &grid_assign_value_between_surfaces,
+                 "Make a property that is one if cell center is between two surfaces.");
 }
 
 }  // namespace xtgeo::grid3d

src/lib/include/xtgeo/regsurf.hpp

@@ -30,6 +30,18 @@ get_xy_from_ij(const size_t i,
                const size_t nrow,
                const double angle_deg);
 
+double
+get_z_from_xy(const double x,
+              const double y,
+              const double xori,
+              const double yori,
+              const double xinc,
+              const double yinc,
+              const size_t ncol,
+              const size_t nrow,
+              const double angle_deg,
+              const py::array_t<double> &values);
+
 std::tuple<int, int, int, int>
 find_cell_range(const double xmin,
                 const double xmax,
@@ -80,6 +92,8 @@ init(py::module &m)
                   "Get the outer corners of a regular surface.");
     m_regsurf.def("get_xy_from_ij", &get_xy_from_ij,
                   "Get the XY coordinates from the grid indices.");
+    m_regsurf.def("get_z_from_xy", &get_z_from_xy,
+                  "Get the Z value in a regsurf from a x y point.");
     m_regsurf.def("find_cell_range", &find_cell_range,
                   "Find the range of regular surface 2D nodes within a box.");
     m_regsurf.def("sample_grid3d_layer", &sample_grid3d_layer,

src/lib/include/xtgeo/xtgeo.h

@@ -1059,26 +1059,6 @@ extern "C"
                         long n_swig_np_dbl_inplace_v1,   // ntot,
                         metric m);
 
-    void grdcp3d_calc_xyz(long ncol,
-                          long nrow,
-                          long nlay,
-
-                          double *swig_np_dbl_in_v1,  // *coordsv,
-                          long n_swig_np_dbl_in_v1,   // ncoord,
-                          float *swig_np_flt_in_v1,   // *zcornsv,
-                          long n_swig_np_flt_in_v1,   // nzcorn,
-                          int *swig_np_int_in_v1,     // *actnumsv,
-                          long n_swig_np_int_in_v1,   // nactnum,
-
-                          int option,
-
-                          double *swig_np_dbl_aout_v1,  // xarr
-                          long n_swig_np_dbl_aout_v1,
-                          double *swig_np_dbl_aout_v2,  // yarr
-                          long n_swig_np_dbl_aout_v2,
-                          double *swig_np_dbl_aout_v3,  // zarr
-                          long n_swig_np_dbl_aout_v3);
-
     void grd3d_conv_roxapi_grid(int nx,
                                 int ny,
                                 int nz,
@@ -1394,34 +1374,6 @@ extern "C"
                               long n_swig_np_int_inplace_v1    // nactnum
     );
 
-    void grdcp3d_midpoint(long i,
-                          long j,
-                          long k,
-                          long ncol,
-                          long nrow,
-                          long nlay,
-                          double *coordsv,
-                          long ncoord,
-                          float *zcornsv,
-                          long nzcorn,
-                          double *x,
-                          double *y,
-                          double *z);
-
-    void grd3d_midpoint(int i,
-                        int j,
-                        int k,
-                        int nx,
-                        int ny,
-                        int nz,
-                        double *coordsv,
-                        long ncoord,
-                        double *zcornsv,
-                        long nzcorn,
-                        double *x,
-                        double *y,
-                        double *z);
-
     int grd3d_inact_outside_pol(double *swig_np_dbl_in_v1,  // polygons X
                                 long n_swig_np_dbl_in_v1,   // N
                                 double *swig_np_dbl_in_v2,  // polygons Y

src/lib/src/grid3d/grid.cpp

@@ -55,6 +55,60 @@ grid_cell_volumes(const size_t ncol,
     return cellvols;
 }
 
+/*
+ * Get cell centers for a grid.
+ *
+ * @param ncol Grid dimensions ncol/nx
+ * @param nrow Grid dimensions nrow/ny
+ * @param nlay Grid dimensions nlay/nz
+ * @param coordsv Grid Z coordinates vector
+ * @param zcornsv Grid Z corners vector
+ * @param actumsv Active cells vector
+ * @param asmasked Process grid cells as masked (return NaN for inactive cells)
+ * @return Arrays with the X, Y, Z coordinates of the cell centers
+ */
+std::tuple<py::array_t<double>, py::array_t<double>, py::array_t<double>>
+grid_cell_centers(const size_t ncol,
+                  const size_t nrow,
+                  const size_t nlay,
+                  const py::array_t<double> &coordsv,
+                  const py::array_t<float> &zcornsv,
+                  const py::array_t<int> &actnumsv,
+                  const bool asmasked = false)
+{
+    pybind11::array_t<double> xmid({ ncol, nrow, nlay });
+    pybind11::array_t<double> ymid({ ncol, nrow, nlay });
+    pybind11::array_t<double> zmid({ ncol, nrow, nlay });
+    auto xmid_ = xmid.mutable_unchecked<3>();
+    auto ymid_ = ymid.mutable_unchecked<3>();
+    auto zmid_ = zmid.mutable_unchecked<3>();
+    auto actnumsv_ = actnumsv.unchecked<3>();
+
+    for (size_t i = 0; i < ncol; i++) {
+        for (size_t j = 0; j < nrow; j++) {
+            for (size_t k = 0; k < nlay; k++) {
+                size_t idx = i * nrow * nlay + j * nlay + k;
+                if (asmasked && actnumsv_(i, j, k) == 0) {
+                    xmid_(i, j, k) = std::numeric_limits<double>::quiet_NaN();
+                    ymid_(i, j, k) = std::numeric_limits<double>::quiet_NaN();
+                    zmid_(i, j, k) = std::numeric_limits<double>::quiet_NaN();
+                    continue;
+                }
+                auto cr =
+                  grid3d::cell_corners(i, j, k, ncol, nrow, nlay, coordsv, zcornsv);
+
+                xmid_(i, j, k) = 0.125 * (cr[0] + cr[3] + cr[6] + cr[9] + cr[12] +
+                                          cr[15] + cr[18] + cr[21]);
+                ymid_(i, j, k) = 0.125 * (cr[1] + cr[4] + cr[7] + cr[10] + cr[13] +
+                                          cr[16] + cr[19] + cr[22]);
+                zmid_(i, j, k) = 0.125 * (cr[2] + cr[5] + cr[8] + cr[11] + cr[14] +
+                                          cr[17] + cr[20] + cr[23]);
+            }
+        }
+    }
+    return std::make_tuple(xmid, ymid, zmid);
+}
+
 /*
  * Compute cell height above ffl (free fluid level), as input to water saturation. Will
  * return hbot, htop, hmid (bottom of cell, top of cell, midpoint), but compute method