ENH: optimize speed for surface_from_grid method

src/lib/CMakeLists.txt

@@ -12,6 +12,10 @@ find_package(
 include(UseSWIG)
 find_package(pybind11 REQUIRED)
 
+if(UNIX AND NOT APPLE)
+  find_package(OpenMP REQUIRED)
+endif()
+
 message(STATUS "Python executable     : ${Python_EXECUTABLE}")
 message(STATUS "Python include dirs   : ${Python_INCLUDE_DIRS}")
 message(STATUS "numpy include path    : ${Python_NumPy_INCLUDE_DIRS}")
@@ -43,6 +47,11 @@ pybind11_add_module(
   "${SRC}/regsurf/utilities.cpp"
   ${SOURCES})
 target_include_directories(_internal PRIVATE ${CMAKE_CURRENT_LIST_DIR}/include)
+# Link OpenMP
+
+if(UNIX AND NOT APPLE AND OpenMP_CXX_FOUND)
+  target_link_libraries(_internal PRIVATE OpenMP::OpenMP_CXX)
+endif()
 
 message(STATUS "Compiling swig bindings")
 

src/lib/include/xtgeo/geometry.hpp

@@ -65,7 +65,7 @@ constexpr int TETRAHEDRON_VERTICES[4][6][4] = {
 };
 
 inline double
-hexahedron_dz(const std::vector<double> &corners)
+hexahedron_dz(const std::array<double, 24> &corners)
 {
     // TODO: This does not account for overall zflip ala Petrel or cells that
     // are malformed
@@ -86,7 +86,7 @@ triangle_area(const std::array<double, 2> &p1,
 }
 
 double
-hexahedron_volume(const std::vector<double> &corners, const int precision);
+hexahedron_volume(const std::array<double, 24> &corners, const int precision);
 
 bool
 is_xy_point_in_polygon(const double x,

src/lib/include/xtgeo/grid3d.hpp

@@ -28,7 +28,7 @@ grid_height_above_ffl(const size_t ncol,
                       const py::array_t<int> &actnumsv,
                       const py::array_t<float> &ffl,
                       const size_t option);
-std::vector<double>
+std::array<double, 24>
 cell_corners(const size_t i,
              const size_t j,
              const size_t k,
@@ -39,33 +39,19 @@ cell_corners(const size_t i,
              const py::array_t<float> &zcornsv);
 
 std::vector<double>
-get_corners_minmax(std::vector<double> &corners);
+get_corners_minmax(std::array<double, 24> &corners);
 
 bool
 is_xy_point_in_cell(const double x,
                     const double y,
-                    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,
-                    const int option);
+                    const std::array<double, 24> &corners,
+                    int option);
 
 double
 get_depth_in_cell(const double x,
                   const double y,
-                  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,
-                  const int option);
+                  const std::array<double, 24> &corners,
+                  int option);
 
 inline void
 init(py::module &m)

src/lib/include/xtgeo/numerics.hpp

@@ -1,8 +1,14 @@
 #ifndef XTGEO_NUMERICS_HPP_
 #define XTGEO_NUMERICS_HPP_
+#include <limits>
 
 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
 {

src/lib/include/xtgeo/regsurf.hpp

@@ -44,7 +44,11 @@ find_cell_range(const double xmin,
                 const size_t nrow,
                 const int expand);
 
-std::tuple<py::array_t<int>, py::array_t<int>, py::array_t<double>>
+std::tuple<py::array_t<int>,
+           py::array_t<int>,
+           py::array_t<double>,
+           py::array_t<double>,
+           py::array_t<bool>>
 sample_grid3d_layer(const size_t ncol,
                     const size_t nrow,
                     const double xori,
@@ -58,9 +62,7 @@ sample_grid3d_layer(const size_t ncol,
                     const py::array_t<double> &coordsv,
                     const py::array_t<float> &zcornsv,
                     const py::array_t<int> &actnumsv,
-                    const size_t klayer,
-                    const int option,
-                    const int activeonly);
+                    const size_t klayer);
 inline void
 init(py::module &m)
 {

src/lib/src/common/geometry/interpolate.cpp

@@ -1,5 +1,6 @@
 #include <xtgeo/constants.hpp>
 #include <xtgeo/geometry.hpp>
+#include <xtgeo/numerics.hpp>
 
 namespace xtgeo::geometry {
 
@@ -28,7 +29,7 @@ interpolate_z_triangles(const double x,
     double z_estimated = std::numeric_limits<double>::quiet_NaN();
 
     // Check if the point is in the triangle
-    if (std::abs(total_area - (area1 + area2 + area3)) < xtgeo::constants::EPSILON) {
+    if (std::abs(total_area - (area1 + area2 + area3)) < numerics::TOLERANCE) {
         double w1 = area1 / total_area;
         double w2 = area2 / total_area;
         double w3 = area3 / total_area;
@@ -45,7 +46,7 @@ interpolate_z_triangles(const double x,
     total_area = triangle_area(p1_2d, p3_2d, p4_2d);
 
     // Check if the point is in the triangle
-    if (std::abs(total_area - (area1 + area2 + area3)) < xtgeo::constants::EPSILON) {
+    if (std::abs(total_area - (area1 + area2 + area3)) < numerics::TOLERANCE) {
         double w1 = area1 / total_area;
         double w2 = area2 / total_area;
         double w3 = area3 / total_area;
@@ -86,7 +87,7 @@ interpolate_z_4p_regular(const double x,
 
     // Quick check if the point is inside the quadrilateral; note ordering of points
     if (!is_xy_point_in_quadrilateral(x, y, p1, p2, p4, p3)) {
-        return std::numeric_limits<double>::quiet_NaN();
+        return numerics::QUIET_NAN;
     }
 
     // Extract coordinates
@@ -127,7 +128,7 @@ interpolate_z_4p(const double x,
 {
     // Quick check if the point is inside the quadrilateral
     if (!is_xy_point_in_quadrilateral(x, y, p1, p2, p4, p3)) {
-        return std::numeric_limits<double>::quiet_NaN();
+        return numerics::QUIET_NAN;
     }
 
     double z1 = interpolate_z_triangles(x, y, p1, p2, p3, p4);

src/lib/src/common/geometry/volumes.cpp

@@ -1,6 +1,7 @@
 #include <pybind11/stl.h>
 #include <vector>
 #include <xtgeo/geometry.hpp>
+#include <xtgeo/numerics.hpp>
 #include <xtgeo/xtgeo.h>
 
 namespace xtgeo::geometry {
@@ -28,10 +29,10 @@ namespace xtgeo::geometry {
  * @return The volume of the hexahadron
  */
 double
-hexahedron_volume(const std::vector<double> &corners, const int precision)
+hexahedron_volume(const std::array<double, 24> &corners, const int precision)
 {
     // Avoid cells that collapsed in some way
-    if (hexahedron_dz(corners) < std::numeric_limits<double>::epsilon()) {
+    if (hexahedron_dz(corners) < numerics::EPSILON) {
         return 0.0;
     }
 

src/lib/src/cube/cube.cpp

@@ -8,6 +8,7 @@
 #include <string>
 #include <unordered_map>
 #include <vector>
+#include <xtgeo/numerics.hpp>
 
 namespace py = pybind11;
 
@@ -131,14 +132,14 @@ cube_stats_along_z(const size_t ncol,
                 meanabsv_(i, j) = sum_abs / n_items;
                 sumabsv_(i, j) = sum_abs;
             } else {
-                minv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                maxv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                meanv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                varv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                rmsv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                maxabsv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                meanabsv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                sumabsv_(i, j) = std::numeric_limits<double>::quiet_NaN();
+                minv_(i, j) = numerics::QUIET_NAN;
+                maxv_(i, j) = numerics::QUIET_NAN;
+                meanv_(i, j) = numerics::QUIET_NAN;
+                varv_(i, j) = numerics::QUIET_NAN;
+                rmsv_(i, j) = numerics::QUIET_NAN;
+                maxabsv_(i, j) = numerics::QUIET_NAN;
+                meanabsv_(i, j) = numerics::QUIET_NAN;
+                sumabsv_(i, j) = numerics::QUIET_NAN;
             }
 
             // Compute statistics for positive values
@@ -147,9 +148,9 @@ cube_stats_along_z(const size_t ncol,
                 meanposv_(i, j) = sum_pos / n_items_pos;
                 sumposv_(i, j) = sum_pos;
             } else {
-                maxposv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                meanposv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                sumposv_(i, j) = std::numeric_limits<double>::quiet_NaN();
+                maxposv_(i, j) = numerics::QUIET_NAN;
+                meanposv_(i, j) = numerics::QUIET_NAN;
+                sumposv_(i, j) = numerics::QUIET_NAN;
             }
 
             // Compute statistics for negative values
@@ -158,9 +159,9 @@ cube_stats_along_z(const size_t ncol,
                 meannegv_(i, j) = sum_neg / n_items_neg;
                 sumnegv_(i, j) = sum_neg;
             } else {
-                maxnegv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                meannegv_(i, j) = std::numeric_limits<double>::quiet_NaN();
-                sumnegv_(i, j) = std::numeric_limits<double>::quiet_NaN();
+                maxnegv_(i, j) = numerics::QUIET_NAN;
+                meannegv_(i, j) = numerics::QUIET_NAN;
+                sumnegv_(i, j) = numerics::QUIET_NAN;
             }
         }
     }