From 6f29e7469224614cc6554ae6dcbdc25787490e16 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Tue, 22 Mar 2022 13:12:58 +0100
Subject: [PATCH] Add Means of Detecting Valid Cell Geometries

We define a valid cell geometry as one in which all components of
all eight corner vertices of a cell have finite coordinate values
(absolute value less than 1.0E+20) and at least one pair of pillar
vertices are separated by a positive distance.
---
 .../eclipse/EclipseState/Grid/EclipseGrid.hpp | 18 ++++---
 .../eclipse/EclipseState/Grid/EclipseGrid.cpp | 54 +++++++++++++++----
 2 files changed, 56 insertions(+), 16 deletions(-)

diff --git a/opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp b/opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp
index f9ed932ee..67c38045b 100644
--- a/opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp
+++ b/opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp
@@ -17,21 +17,20 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-
 #ifndef OPM_PARSER_ECLIPSE_GRID_HPP
 #define OPM_PARSER_ECLIPSE_GRID_HPP
 
-
+#include <opm/input/eclipse/EclipseState/Grid/GridDims.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/MapAxes.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/MinpvMode.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/PinchMode.hpp>
-#include <opm/input/eclipse/EclipseState/Grid/GridDims.hpp>
 
 #include <array>
 #include <memory>
 #include <optional>
-#include <vector>
+#include <stdexcept>
 #include <unordered_set>
+#include <vector>
 
 namespace Opm {
 
@@ -177,6 +176,14 @@ namespace Opm {
             return cellActive(i, j, k);
         }
 
+        /// Whether or not given cell has a valid cell geometry
+        ///
+        /// Valid geometry is defined as all vertices have finite
+        /// coordinates and at least one pair of coordinates are separated
+        /// by a physical distance along a pillar.
+        bool isValidCellGeomtry(const std::size_t globalIndex,
+                                const UnitSystem& usys) const;
+
         double getCellDepth(size_t i,size_t j, size_t k) const;
         double getCellDepth(size_t globalIndex) const;
         ZcornMapper zcornMapper() const;
@@ -329,7 +336,4 @@ namespace Opm {
     };
 }
 
-
-
-
 #endif // OPM_PARSER_ECLIPSE_GRID_HPP
diff --git a/src/opm/input/eclipse/EclipseState/Grid/EclipseGrid.cpp b/src/opm/input/eclipse/EclipseState/Grid/EclipseGrid.cpp
index 0bdae226b..11c55820e 100644
--- a/src/opm/input/eclipse/EclipseState/Grid/EclipseGrid.cpp
+++ b/src/opm/input/eclipse/EclipseState/Grid/EclipseGrid.cpp
@@ -19,15 +19,8 @@
 */
 
 #define _USE_MATH_DEFINES
-#include <cmath>
-#include <cstring>
-#include <numeric>
 
-#include <iostream>
-#include <tuple>
-#include <functional>
-
-#include <fmt/format.h>
+#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
@@ -60,8 +53,16 @@
 #include <opm/input/eclipse/Parser/ParserKeywords/T.hpp>
 #include <opm/input/eclipse/Parser/ParserKeywords/Z.hpp>
 
-#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <algorithm>
+#include <cmath>
+#include <cstring>
+#include <functional>
+#include <initializer_list>
+#include <iostream>
+#include <numeric>
+#include <tuple>
 
+#include <fmt/format.h>
 
 namespace Opm {
 
@@ -1597,6 +1598,41 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
         }
     }
 
+    bool EclipseGrid::isValidCellGeomtry(const std::size_t globalIndex,
+                                         const UnitSystem& usys) const
+    {
+        const auto threshold = usys.to_si(UnitSystem::measure::length, 1.0e+20f);
+
+        auto is_finite = [threshold](const double coord_component)
+        {
+            return std::abs(coord_component) < threshold;
+        };
+
+        std::array<double,8> X, Y, Z;
+        this->getCellCorners(globalIndex, X, Y, Z);
+
+        const auto finite_coord = std::all_of(X.begin(), X.end(), is_finite)
+            && std::all_of(Y.begin(), Y.end(), is_finite)
+            && std::all_of(Z.begin(), Z.end(), is_finite);
+
+        if (! finite_coord) {
+            return false;
+        }
+
+        const auto max_pillar_point_distance = std::max({
+            Z[0 + 4] - Z[0 + 0],
+            Z[1 + 4] - Z[1 + 0],
+            Z[2 + 4] - Z[2 + 0],
+            Z[3 + 4] - Z[3 + 0],
+        });
+
+        // Define points as "well separated" if maximum distance exceeds
+        // 1e-4 length units (e.g., 0.1 mm in METRIC units).  May consider
+        // using a coarser tolerance/threshold here.
+        return max_pillar_point_distance
+            >  usys.to_si(UnitSystem::measure::length, 1.0e-4);
+    }
+
     double EclipseGrid::getCellDepth(size_t globalIndex) const {
         assertGlobalIndex( globalIndex );
         std::array<double,8> X;