CLN: move towards structs in C++/pybind11

.clang-format

@@ -41,7 +41,7 @@ BreakBeforeBraces: Mozilla
 BreakBeforeInheritanceComma: true
 BreakBeforeTernaryOperators: true
 BreakConstructorInitializersBeforeComma: false
-BreakConstructorInitializers: BeforeComma
+BreakConstructorInitializers: AfterColon
 BreakAfterJavaFieldAnnotations: false
 BreakStringLiterals: true
 ColumnLimit:     88

src/lib/include/xtgeo/cube.hpp

@@ -4,25 +4,41 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
 #include <cstddef>
+#include <stdexcept>
+#include <xtgeo/types.hpp>
 
 namespace py = pybind11;
 
 namespace xtgeo::cube {
 
 std::unordered_map<std::string, py::array_t<double>>
-cube_stats_along_z(const size_t ncol,
-                   const size_t nrow,
-                   const size_t nlay,
-                   const py::array_t<float> &cubev);
+cube_stats_along_z(const Cube &cube_cpp);
 
 inline void
 init(py::module &m)
 {
     auto m_cube =
       m.def_submodule("cube", "Internal functions for operations on 3d cubes.");
 
-    m_cube.def("cube_stats_along_z", &cube_stats_along_z,
-               "Compute various statistics for cube along the Z axis, returning maps.");
+    py::class_<Cube>(m_cube, "Cube")
+      .def(py::init<>())
+      .def(py::init<const py::object &>(), py::arg("cube"))
+      .def_readonly("ncol", &Cube::ncol)
+      .def_readonly("nrow", &Cube::nrow)
+      .def_readonly("nlay", &Cube::nlay)
+      .def_readonly("xori", &Cube::xori)
+      .def_readonly("yori", &Cube::yori)
+      .def_readonly("zori", &Cube::zori)
+      .def_readonly("xinc", &Cube::xinc)
+      .def_readonly("yinc", &Cube::yinc)
+      .def_readonly("zinc", &Cube::zinc)
+      .def_readonly("rotation", &Cube::rotation)
+      .def_readonly("values", &Cube::values)
+
+      .def("cube_stats_along_z", &cube_stats_along_z,
+           "Compute various statistics for cube along the Z axis, returning maps.")
+
+      ;
 }
 
 }  // namespace xtgeo::cube

src/lib/include/xtgeo/geometry.hpp

@@ -6,10 +6,7 @@
 #include <array>
 #include <cmath>
 #include <vector>
-
-#ifndef M_PI
-#define M_PI 3.14159265358979323846  // seems like Windows does not define M_PI i cmath
-#endif
+#include <xtgeo/types.hpp>
 
 namespace py = pybind11;
 
@@ -65,57 +62,53 @@ constexpr int TETRAHEDRON_VERTICES[4][6][4] = {
 };
 
 inline double
-hexahedron_dz(const std::array<double, 24> &corners)
+hexahedron_dz(const grid3d::CellCorners &corners)
 {
     // TODO: This does not account for overall zflip ala Petrel or cells that
     // are malformed
     double dzsum = 0.0;
-    for (auto i = 0; i < 4; i++) {
-        dzsum += std::abs(corners[3 * i + 2] - corners[3 * i + 2 + 12]);
-    }
+    dzsum += std::abs(corners.upper_sw.z - corners.lower_sw.z);
+    dzsum += std::abs(corners.upper_se.z - corners.lower_se.z);
+    dzsum += std::abs(corners.upper_nw.z - corners.lower_nw.z);
+    dzsum += std::abs(corners.upper_ne.z - corners.lower_ne.z);
     return dzsum / 4.0;
 }
 
 inline double
-triangle_area(const std::array<double, 2> &p1,
-              const std::array<double, 2> &p2,
-              const std::array<double, 2> &p3)
+triangle_area(const xyz::Point &p1, const xyz::Point &p2, const xyz::Point &p3)
 {
-    return 0.5 * std::abs(p1[0] * (p2[1] - p3[1]) + p2[0] * (p3[1] - p1[1]) +
-                          p3[0] * (p1[1] - p2[1]));
+    return 0.5 *
+           std::abs(p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y));
 }
 
 double
-hexahedron_volume(const std::array<double, 24> &corners, const int precision);
+hexahedron_volume(const grid3d::CellCorners &corners, const int precision);
 
 bool
-is_xy_point_in_polygon(const double x,
-                       const double y,
-                       const std::vector<std::array<double, 2>> &polygon);
+is_xy_point_in_polygon(const double x, const double y, const xyz::Polygon &polygon);
 
 bool
 is_xy_point_in_quadrilateral(const double x,
                              const double y,
-                             const std::array<double, 3> &p1,
-                             const std::array<double, 3> &p2,
-                             const std::array<double, 3> &p3,
-                             const std::array<double, 3> &p4);
-
+                             const xyz::Point &p1,
+                             const xyz::Point &p2,
+                             const xyz::Point &p3,
+                             const xyz::Point &p4);
 double
 interpolate_z_4p_regular(const double x,
                          const double y,
-                         const std::array<double, 3> &p1,
-                         const std::array<double, 3> &p2,
-                         const std::array<double, 3> &p3,
-                         const std::array<double, 3> &p4);
+                         const xyz::Point &p1,
+                         const xyz::Point &p2,
+                         const xyz::Point &p3,
+                         const xyz::Point &p4);
 
 double
 interpolate_z_4p(const double x,
                  const double y,
-                 const std::array<double, 3> &p1,
-                 const std::array<double, 3> &p2,
-                 const std::array<double, 3> &p3,
-                 const std::array<double, 3> &p4);
+                 const xyz::Point &p1,
+                 const xyz::Point &p2,
+                 const xyz::Point &p3,
+                 const xyz::Point &p4);
 
 // functions exposed to Python:
 inline void

src/lib/include/xtgeo/grid3d.hpp

@@ -5,112 +5,102 @@
 #include <cstddef>
 #include <cstdint>
 #include <optional>
+#include <stdexcept>
 #include <tuple>
+#include <xtgeo/types.hpp>
 
 namespace py = pybind11;
 
 namespace xtgeo::grid3d {
 
 py::array_t<double>
-grid_cell_volumes(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 int precision,
-                  const bool asmasked);
+get_cell_volumes(const Grid &grid_cpp, 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);
+get_cell_centers(const Grid &grid_cpp, 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,
-                      const size_t nlay,
-                      const py::array_t<double> &coordsv,
-                      const py::array_t<float> &zcornsv,
-                      const py::array_t<int> &actnumsv,
-                      const py::array_t<float> &ffl,
-                      const size_t option);
-std::array<double, 24>
-cell_corners(const size_t i,
-             const size_t j,
-             const size_t k,
-             const size_t ncol,
-             const size_t nrow,
-             const size_t nlay,
-             const py::array_t<double> &coordsv,
-             const py::array_t<float> &zcornsv);
+get_height_above_ffl(const Grid &grid_cpp,
+                     const py::array_t<float> &ffl,
+                     const size_t option);
+
+CellCorners
+get_cell_corners_from_ijk(const Grid &grid_cpp,
+                          const size_t i,
+                          const size_t j,
+                          const size_t k);
 
 std::vector<double>
-get_corners_minmax(std::array<double, 24> &corners);
+get_corners_minmax(CellCorners &get_cell_corners_from_ijk);
 
 bool
 is_xy_point_in_cell(const double x,
                     const double y,
-                    const std::array<double, 24> &corners,
+                    const CellCorners &corners,
                     int option);
 
 double
 get_depth_in_cell(const double x,
                   const double y,
-                  const std::array<double, 24> &corners,
+                  const CellCorners &corners,
                   int option);
 
 py::array_t<int8_t>
-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);
+get_gridprop_value_between_surfaces(const Grid &grd,
+                                    const regsurf::RegularSurface &top,
+                                    const regsurf::RegularSurface &bot);
 
 inline void
 init(py::module &m)
 {
     auto m_grid3d =
       m.def_submodule("grid3d", "Internal functions for operations on 3d grids.");
 
-    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,
-                 "Get a vector containing the corners of a cell.");
+    py::class_<Grid>(m_grid3d, "Grid")
+      .def(py::init<>())
+      .def(py::init<const py::object &>(), py::arg("grid"))
+
+      .def_readonly("ncol", &Grid::ncol)
+      .def_readonly("nrow", &Grid::nrow)
+      .def_readonly("nlay", &Grid::nlay)
+      .def_readonly("coordsv", &Grid::coordsv)
+      .def_readonly("zcornsv", &Grid::zcornsv)
+      .def_readonly("actnumsv", &Grid::actnumsv)
+
+      .def("get_cell_volumes", &get_cell_volumes, "Compute the bulk volume of cell.")
+
+      .def("get_cell_centers", &get_cell_centers,
+           "Compute the cells centers coordinates as 3 arrays")
+      .def("get_gridprop_value_between_surfaces", &get_gridprop_value_between_surfaces,
+           "Make a property that is one if cell center is between two surfaces.")
+      .def("get_height_above_ffl", &get_height_above_ffl,
+           "Compute the height above a FFL (free fluid level).")
+      .def("get_cell_corners_from_ijk", &get_cell_corners_from_ijk,
+           "Get a vector containing the corners of a specified IJK cell.")
+
+      ;
+
+    py::class_<CellCorners>(m_grid3d, "CellCorners")
+      .def(py::init<>())
+      .def_readonly("upper_sw", &CellCorners::upper_sw)
+      .def_readonly("upper_se", &CellCorners::upper_se)
+      .def_readonly("upper_nw", &CellCorners::upper_nw)
+      .def_readonly("upper_ne", &CellCorners::upper_ne)
+      .def_readonly("lower_sw", &CellCorners::lower_sw)
+      .def_readonly("lower_se", &CellCorners::lower_se)
+      .def_readonly("lower_nw", &CellCorners::lower_nw)
+      .def_readonly("lower_ne", &CellCorners::lower_ne)
+      .def("to_numpy", &CellCorners::to_numpy);
+
+    m_grid3d.def("arrange_corners", &CellCorners::arrange_corners,
+                 "Arrange the corners in a single array for easier access.");
+
     m_grid3d.def("get_corners_minmax", &get_corners_minmax,
                  "Get a vector containing the minmax of a single corner set");
     m_grid3d.def("is_xy_point_in_cell", &is_xy_point_in_cell,
                  "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/logging.hpp

@@ -11,8 +11,7 @@ namespace py = pybind11;
 class Logger
 {
 public:
-    explicit Logger(const std::string &name)
-      : logger_name(name)
+    explicit Logger(const std::string &name) : logger_name(name)
     {
         py::gil_scoped_acquire acquire;  // Acquire GIL (Global Interpreter Lock)
         py::object logging = py::module::import("logging");

src/lib/include/xtgeo/numerics.hpp

@@ -2,36 +2,32 @@
 #define XTGEO_NUMERICS_HPP_
 #include <pybind11/pybind11.h>
 #include <limits>
+#include <xtgeo/types.hpp>
 
 namespace py = pybind11;
 
 namespace xtgeo::numerics {
 
-constexpr double UNDEF_DOUBLE = std::numeric_limits<double>::max();
-constexpr double EPSILON = std::numeric_limits<double>::epsilon();
-constexpr double TOLERANCE = 1e-6;
-constexpr double QUIET_NAN = std::numeric_limits<double>::quiet_NaN();
-
 template<typename T>
 struct Vec3
 {
     T x, y, z;
 };
 
-inline Vec3<double>
+inline xyz::Point
 lerp3d(double x1, double y1, double z1, double x2, double y2, double z2, double t)
 {
-    return Vec3<double>{ x1 + t * (x2 - x1), y1 + t * (y2 - y1), z1 + t * (z2 - z1) };
+    return xyz::Point{ x1 + t * (x2 - x1), y1 + t * (y2 - y1), z1 + t * (z2 - z1) };
 }
 
 inline void
 init(py::module &m)
 {
     auto m_numerics = m.def_submodule("numerics", "Internal functions for numerics.");
 
-    m_numerics.attr("UNDEF_DOUBLE") = UNDEF_DOUBLE;
-    m_numerics.attr("EPSILON") = EPSILON;
-    m_numerics.attr("TOLERANCE") = TOLERANCE;
+    m_numerics.attr("UNDEF_DOUBLE") = numerics::UNDEF_DOUBLE;
+    m_numerics.attr("EPSILON") = numerics::EPSILON;
+    m_numerics.attr("TOLERANCE") = numerics::TOLERANCE;
 }
 
 }  // namespace xtgeo::numerics