#15088 Well Log Extractor: Use parametric tolerance for well path cell intersections

.github/workflows/ResInsightMac.yml

@@ -223,12 +223,11 @@ jobs:
       - name: Run pytest
         env:
           RESINSIGHT_EXECUTABLE: ${{ runner.workspace }}/ResInsight/cmakebuild/ResInsight.app/Contents/MacOS/ResInsight
-        # Deselected tests differ from the Linux/Windows reference values on macOS, see #15076
+        # The deselected test differs from the Linux/Windows reference values on macOS, see #15076
         run: |
           cd GrpcInterface/Python/rips
           ${{ steps.python-path.outputs.PYTHON_EXECUTABLE }} -m pytest --console \
-            --deselect rips/tests/test_completion_export.py::test_export_completion_files_unified \
-            --deselect rips/tests/test_simulation_wells.py::test_10k_acclength
+            --deselect rips/tests/test_completion_export.py::test_export_completion_files_unified
 
       - name: Write macOS launch instructions
         run: |

ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.cpp

@@ -33,7 +33,8 @@ int RigHexIntersectionTools::lineHexCellIntersection( const cvf::Vec3d&
                                                       const cvf::Vec3d&                 p2,
                                                       const std::array<cvf::Vec3d, 8>&  hexCorners,
                                                       const size_t                      hexIndex,
-                                                      std::vector<HexIntersectionInfo>* intersections )
+                                                      std::vector<HexIntersectionInfo>* intersections,
+                                                      double                            parametricTolerance )
 {
     CAF_ASSERT( intersections != nullptr );
 
@@ -61,7 +62,8 @@ int RigHexIntersectionTools::lineHexCellIntersection( const cvf::Vec3d&
                                                                                hexCorners[faceVertexIndices[next]],
                                                                                faceCenter,
                                                                                &intersection,
-                                                                               &isEntering );
+                                                                               &isEntering,
+                                                                               parametricTolerance );
             if ( intsStatus == 1 )
             {
                 uniqueIntersections.insert(

ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.h

@@ -56,7 +56,8 @@ struct RigHexIntersectionTools
                                         const cvf::Vec3d&                 p2,
                                         const std::array<cvf::Vec3d, 8>&  hexCorners,
                                         const size_t                      hexIndex,
-                                        std::vector<HexIntersectionInfo>* intersections );
+                                        std::vector<HexIntersectionInfo>* intersections,
+                                        double                            parametricTolerance = 0.0 );
 
     static bool lineIntersectsHexCell( const cvf::Vec3d& p1, const cvf::Vec3d& p2, const std::array<cvf::Vec3d, 8>& hexCorners );
 

ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp

@@ -89,7 +89,12 @@ void RigEclipseWellLogExtractor::calculateIntersection()
             if ( cell.isInvalid() || cell.subGrid() != nullptr ) continue;
 
             std::array<cvf::Vec3d, 8> hexCorners = m_caseData->mainGrid()->cellCornerVertices( globalCellIndex );
-            RigHexIntersectionTools::lineHexCellIntersection( p1, p2, hexCorners, globalCellIndex, &intersections );
+            RigHexIntersectionTools::lineHexCellIntersection( p1,
+                                                              p2,
+                                                              hexCorners,
+                                                              globalCellIndex,
+                                                              &intersections,
+                                                              RigWellLogExtractionTools::intersectionParametricTolerance );
         }
 
         if ( !isCellFaceNormalsOut )

ApplicationLibCode/ReservoirDataModel/Well/RigGeoMechWellLogExtractor.cpp

@@ -1161,9 +1161,12 @@ void RigGeoMechWellLogExtractor::calculateIntersection()
             hexCorners[6] = cvf::Vec3d( nodeCoords[cornerIndices[6]] );
             hexCorners[7] = cvf::Vec3d( nodeCoords[cornerIndices[7]] );
 
-            // int intersectionCount = RigHexIntersector::lineHexCellIntersection(p1, p2, hexCorners,
-            // closeCells[ccIdx], &intersections);
-            RigHexIntersectionTools::lineHexCellIntersection( p1, p2, hexCorners, closeCells[ccIdx], &intersections );
+            RigHexIntersectionTools::lineHexCellIntersection( p1,
+                                                              p2,
+                                                              hexCorners,
+                                                              closeCells[ccIdx],
+                                                              &intersections,
+                                                              RigWellLogExtractionTools::intersectionParametricTolerance );
         }
 
         // Now, with all the intersections of this piece of line, we need to

ApplicationLibCode/ReservoirDataModel/Well/RigWellLogExtractionTools.h

@@ -30,6 +30,13 @@
 
 struct RigWellLogExtractionTools
 {
+    // Parametric tolerance for the line segment/triangle test when a well path is intersected with the grid cells.
+    // Simulation well centre lines go through cell face centres, which are also the apex of the four triangles a
+    // face is split into, so the exact test (tolerance 0) depends on the rounding of the last bit and can miss the
+    // crossing on both segments that share the point (#15088). On a 100 m segment this is 1e-7 m, well above the
+    // rounding error at UTM coordinates and well below the measured depth tolerance.
+    static constexpr double intersectionParametricTolerance = 1e-9;
+
     static bool isEqualDepth( double d1, double d2, const double tolerance )
     {
         double depthDiff = d1 - d2;

ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -21,6 +21,8 @@
 
 #include "cafAssert.h"
 
+#include <algorithm>
+
 namespace cvf
 {
 //--------------------------------------------------------------------------------------------------
@@ -488,6 +490,9 @@ double GeometryTools::linePointSquareDist( const cvf::Vec3d& p1, const cvf::Vec3
 //
 /// Intersect a line segment with a 3D triangle
 ///    Input:  A line segment p0, p1. A triangle t0, t1, t2.
+///            parametricTolerance: the segment parameter r and the triangle coordinates s, t are accepted
+///            within [-parametricTolerance, 1 + parametricTolerance]. r is clamped to [0, 1] before the point is
+///            computed, so the reported point is always on the segment. The default 0.0 is the exact test.
 ///    Output: *intersectionPoint = intersection point (when it exists)
 ///    Return: -1 = triangle is degenerate (a segment or point)
 ///             0 = disjoint (no intersect)
@@ -505,7 +510,8 @@ int GeometryTools::intersectLineSegmentTriangle( const cvf::Vec3d& p0,
                                                  const cvf::Vec3d& t1,
                                                  const cvf::Vec3d& t2,
                                                  cvf::Vec3d*       intersectionPoint,
-                                                 bool*             isLineDirDotNormalNegative )
+                                                 bool*             isLineDirDotNormalNegative,
+                                                 double            parametricTolerance )
 {
     CAF_ASSERT( intersectionPoint != nullptr );
     CAF_ASSERT( isLineDirDotNormalNegative != nullptr );
@@ -538,12 +544,14 @@ int GeometryTools::intersectLineSegmentTriangle( const cvf::Vec3d& p0,
 
     // get intersect point of ray with triangle plane
     r = a / b;
-    if ( r < 0.0 ) // ray goes away from triangle
+    if ( r < -parametricTolerance ) // ray goes away from triangle
         return 0; // => no intersect
 
-    if ( r > 1.0 ) // Line segment does not reach triangle
+    if ( r > 1.0 + parametricTolerance ) // Line segment does not reach triangle
         return 0;
 
+    r = std::clamp( r, 0.0, 1.0 ); // keep the point on the segment when the tolerance accepted it
+
     *intersectionPoint = p0 + r * dir; // intersect point of ray and plane
 
     // is I inside T?
@@ -559,11 +567,11 @@ int GeometryTools::intersectLineSegmentTriangle( const cvf::Vec3d& p0,
     // get and test parametric coords
     double s, t;
     s = ( uv * wv - vv * wu ) / D;
-    if ( s < 0.0 || s > 1.0 ) // I is outside T
+    if ( s < -parametricTolerance || s > 1.0 + parametricTolerance ) // I is outside T
         return 0;
 
     t = ( uv * wu - uu * wv ) / D;
-    if ( t < 0.0 || ( s + t ) > 1.0 ) // I is outside T
+    if ( t < -parametricTolerance || ( s + t ) > 1.0 + parametricTolerance ) // I is outside T
         return 0;
 
     return 1; // I is in T

ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h

@@ -54,7 +54,8 @@ public:
                                                     const cvf::Vec3d& t1,
                                                     const cvf::Vec3d& t2,
                                                     cvf::Vec3d*       intersectionPoint,
-                                                    bool*             isLineDirDotNormalNegative );
+                                                    bool*             isLineDirDotNormalNegative,
+                                                    double            parametricTolerance = 0.0 );
     static cvf::Vec3d barycentricCoords( const cvf::Vec3d& t0, const cvf::Vec3d& t1, const cvf::Vec3d& t2, const cvf::Vec3d& p );
     static cvf::Vec4d
         barycentricCoords( const cvf::Vec3d& v0, const cvf::Vec3d& v1, const cvf::Vec3d& v2, const cvf::Vec3d& v3, const cvf::Vec3d& p );

ApplicationLibCode/UnitTests/RigHexIntersectionTools-Test.cpp

@@ -23,6 +23,13 @@
 #include <QElapsedTimer>
 
 #include "RigHexIntersectionTools.h"
+#include "Well/RigWellLogExtractionTools.h"
+
+#include "cvfGeometryTools.h"
+#include "cvfStructGrid.h"
+
+#include <iostream>
+#include <random>
 
 //--------------------------------------------------------------------------------------------------
 ///
@@ -97,3 +104,166 @@ TEST( RigHexIntersectionTools, DISABLED_planeHexCellIntersectionPerformanceTest
 
     qDebug() << "Time total elapsed: " << timeTotal.elapsed() << " intersected: " << static_cast<int>( isCellIntersected );
 }
+
+//--------------------------------------------------------------------------------------------------
+/// Hex corners in ResInsight corner order, slightly distorted, at UTM-like coordinates
+//--------------------------------------------------------------------------------------------------
+static std::array<cvf::Vec3d, 8> createDistortedHex( std::mt19937& rng, const cvf::Vec3d& origin )
+{
+    const double dx = 100.0;
+    const double dy = 100.0;
+    const double dz = 10.0;
+
+    std::array<cvf::Vec3d, 8> hexCorners = { origin + cvf::Vec3d( 0, 0, 0 ),
+                                             origin + cvf::Vec3d( dx, 0, 0 ),
+                                             origin + cvf::Vec3d( dx, dy, 0 ),
+                                             origin + cvf::Vec3d( 0, dy, 0 ),
+                                             origin + cvf::Vec3d( 0, 0, dz ),
+                                             origin + cvf::Vec3d( dx, 0, dz ),
+                                             origin + cvf::Vec3d( dx, dy, dz ),
+                                             origin + cvf::Vec3d( 0, dy, dz ) };
+
+    std::uniform_real_distribution<double> jitter( -5.0, 5.0 );
+    for ( auto& corner : hexCorners )
+    {
+        corner += cvf::Vec3d( jitter( rng ), jitter( rng ), 0.1 * jitter( rng ) );
+    }
+
+    return hexCorners;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Count the face centres where neither of the two well path segments sharing the point reports the crossing
+//--------------------------------------------------------------------------------------------------
+static int countMissedFaceCenterCrossings( double parametricTolerance, int hexCount )
+{
+    std::mt19937                           rng( 15088 );
+    std::uniform_real_distribution<double> easting( 400000.0, 600000.0 );
+    std::uniform_real_distribution<double> northing( 6000000.0, 7000000.0 );
+    std::uniform_real_distribution<double> depth( -4000.0, -1000.0 );
+
+    int missed = 0;
+
+    for ( int hexIdx = 0; hexIdx < hexCount; ++hexIdx )
+    {
+        std::array<cvf::Vec3d, 8> hexCorners = createDistortedHex( rng, cvf::Vec3d( easting( rng ), northing( rng ), depth( rng ) ) );
+
+        cvf::Vec3d cellCenter = cvf::Vec3d::ZERO;
+        for ( const auto& corner : hexCorners )
+            cellCenter += corner;
+        cellCenter *= 1.0 / 8.0;
+
+        for ( int face = 0; face < 6; ++face )
+        {
+            cvf::ubyte faceVertexIndices[4];
+            cvf::StructGridInterface::cellFaceVertexIndices( static_cast<cvf::StructGridInterface::FaceType>( face ), faceVertexIndices );
+
+            // Same expression as in RigHexIntersectionTools::lineHexCellIntersection, so the well path point is
+            // bit-identical to the triangle apex
+            cvf::Vec3d faceCenter = cvf::GeometryTools::computeFaceCenter( hexCorners[faceVertexIndices[0]],
+                                                                           hexCorners[faceVertexIndices[1]],
+                                                                           hexCorners[faceVertexIndices[2]],
+                                                                           hexCorners[faceVertexIndices[3]] );
+
+            cvf::Vec3d outsidePoint = faceCenter + ( faceCenter - cellCenter );
+
+            std::vector<HexIntersectionInfo> intersections;
+            RigHexIntersectionTools::lineHexCellIntersection( outsidePoint, faceCenter, hexCorners, 0, &intersections, parametricTolerance );
+            RigHexIntersectionTools::lineHexCellIntersection( faceCenter, cellCenter, hexCorners, 0, &intersections, parametricTolerance );
+
+            bool found = false;
+            for ( const auto& intersection : intersections )
+            {
+                if ( intersection.m_face == static_cast<cvf::StructGridInterface::FaceType>( face ) &&
+                     intersection.m_intersectionPoint.pointDistance( faceCenter ) < 1e-6 )
+                    found = true;
+            }
+
+            if ( !found ) ++missed;
+        }
+    }
+
+    return missed;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// #15088: a well path vertex on a cell face centre must be reported by one of the two segments sharing it
+//--------------------------------------------------------------------------------------------------
+TEST( RigHexIntersectionTools, lineHexCellIntersectionAtFaceCenters )
+{
+    const int hexCount = 500;
+
+    int missedExact = countMissedFaceCenterCrossings( 0.0, hexCount );
+    std::cout << "Face centre crossings missed with parametric tolerance 0: " << missedExact << " of " << 6 * hexCount << std::endl;
+
+    int missedWithTolerance = countMissedFaceCenterCrossings( RigWellLogExtractionTools::intersectionParametricTolerance, hexCount );
+    EXPECT_EQ( 0, missedWithTolerance );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// The parametric tolerance does not change the result for segments clearly inside or outside the cell, and a
+/// segment ending just short of a face is only accepted with a tolerance, with the point kept on the segment
+//--------------------------------------------------------------------------------------------------
+TEST( RigHexIntersectionTools, lineHexCellIntersectionParametricTolerance )
+{
+    std::array<cvf::Vec3d, 8> hexCorners = { cvf::Vec3d( 0, 0, 0 ),
+                                             cvf::Vec3d( 1, 0, 0 ),
+                                             cvf::Vec3d( 1, 1, 0 ),
+                                             cvf::Vec3d( 0, 1, 0 ),
+                                             cvf::Vec3d( 0, 0, 1 ),
+                                             cvf::Vec3d( 1, 0, 1 ),
+                                             cvf::Vec3d( 1, 1, 1 ),
+                                             cvf::Vec3d( 0, 1, 1 ) };
+
+    const double tolerance = RigWellLogExtractionTools::intersectionParametricTolerance;
+
+    {
+        std::vector<HexIntersectionInfo> exact;
+        std::vector<HexIntersectionInfo> withTolerance;
+        RigHexIntersectionTools::lineHexCellIntersection( cvf::Vec3d( 0.5, 0.5, -1 ), cvf::Vec3d( 0.5, 0.5, 2 ), hexCorners, 0, &exact );
+        RigHexIntersectionTools::lineHexCellIntersection( cvf::Vec3d( 0.5, 0.5, -1 ),
+                                                          cvf::Vec3d( 0.5, 0.5, 2 ),
+                                                          hexCorners,
+                                                          0,
+                                                          &withTolerance,
+                                                          tolerance );
+
+        ASSERT_EQ( 2u, exact.size() );
+        ASSERT_EQ( 2u, withTolerance.size() );
+        for ( size_t i = 0; i < exact.size(); ++i )
+        {
+            EXPECT_EQ( exact[i].m_face, withTolerance[i].m_face );
+            EXPECT_EQ( exact[i].m_isIntersectionEntering, withTolerance[i].m_isIntersectionEntering );
+            EXPECT_EQ( exact[i].m_intersectionPoint, withTolerance[i].m_intersectionPoint );
+        }
+    }
+
+    {
+        std::vector<HexIntersectionInfo> exact;
+        std::vector<HexIntersectionInfo> withTolerance;
+        RigHexIntersectionTools::lineHexCellIntersection( cvf::Vec3d( 0.5, 0.5, -1 ), cvf::Vec3d( 0.5, 0.5, -0.1 ), hexCorners, 0, &exact );
+        RigHexIntersectionTools::lineHexCellIntersection( cvf::Vec3d( 0.5, 0.5, -1 ),
+                                                          cvf::Vec3d( 0.5, 0.5, -0.1 ),
+                                                          hexCorners,
+                                                          0,
+                                                          &withTolerance,
+                                                          tolerance );
+
+        EXPECT_EQ( 0u, exact.size() );
+        EXPECT_EQ( 0u, withTolerance.size() );
+    }
+
+    {
+        const cvf::Vec3d p1( 0.5, 0.5, -1 );
+        const cvf::Vec3d p2( 0.5, 0.5, -1e-12 );
+
+        std::vector<HexIntersectionInfo> exact;
+        std::vector<HexIntersectionInfo> withTolerance;
+        RigHexIntersectionTools::lineHexCellIntersection( p1, p2, hexCorners, 0, &exact );
+        RigHexIntersectionTools::lineHexCellIntersection( p1, p2, hexCorners, 0, &withTolerance, tolerance );
+
+        EXPECT_EQ( 0u, exact.size() );
+        ASSERT_EQ( 1u, withTolerance.size() );
+        EXPECT_LT( withTolerance[0].m_intersectionPoint.pointDistance( p2 ), 1e-12 );
+    }
+}

ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp

@@ -766,3 +766,43 @@ TEST( GeometryToolsEdgeCases, PolygonVariations )
     double hexagonArea = GeometryTools::polygonArea( hexagonNodes );
     EXPECT_GT( hexagonArea, 0.0 );
 }
+
+//--------------------------------------------------------------------------------------------------
+/// intersectLineSegmentTriangle with a parametric tolerance accepts hits just outside r, s, t and s + t
+//--------------------------------------------------------------------------------------------------
+TEST( GeometryToolsEdgeCases, LineSegmentTriangleParametricTolerance )
+{
+    const cvf::Vec3d t0( 0, 0, 0 );
+    const cvf::Vec3d t1( 1, 0, 0 );
+    const cvf::Vec3d t2( 0, 1, 0 );
+
+    cvf::Vec3d intersection;
+    bool       isEntering = false;
+
+    // Through the apex t0: exact
+    EXPECT_EQ( 1,
+               GeometryTools::intersectLineSegmentTriangle( cvf::Vec3d( 0, 0, -1 ), cvf::Vec3d( 0, 0, 1 ), t0, t1, t2, &intersection, &isEntering ) );
+
+    // Just outside the edge s = 0
+    const cvf::Vec3d p0( -1e-12, 0.5, -1 );
+    const cvf::Vec3d p1( -1e-12, 0.5, 1 );
+    EXPECT_EQ( 0, GeometryTools::intersectLineSegmentTriangle( p0, p1, t0, t1, t2, &intersection, &isEntering ) );
+    EXPECT_EQ( 1, GeometryTools::intersectLineSegmentTriangle( p0, p1, t0, t1, t2, &intersection, &isEntering, 1e-9 ) );
+
+    // Just outside the edge s + t = 1
+    const cvf::Vec3d q0( 0.5 + 1e-12, 0.5, -1 );
+    const cvf::Vec3d q1( 0.5 + 1e-12, 0.5, 1 );
+    EXPECT_EQ( 0, GeometryTools::intersectLineSegmentTriangle( q0, q1, t0, t1, t2, &intersection, &isEntering ) );
+    EXPECT_EQ( 1, GeometryTools::intersectLineSegmentTriangle( q0, q1, t0, t1, t2, &intersection, &isEntering, 1e-9 ) );
+
+    // Segment ending just short of the plane: accepted with tolerance, point clamped to the segment end
+    const cvf::Vec3d r0( 0.25, 0.25, -1 );
+    const cvf::Vec3d r1( 0.25, 0.25, -1e-12 );
+    EXPECT_EQ( 0, GeometryTools::intersectLineSegmentTriangle( r0, r1, t0, t1, t2, &intersection, &isEntering ) );
+    EXPECT_EQ( 1, GeometryTools::intersectLineSegmentTriangle( r0, r1, t0, t1, t2, &intersection, &isEntering, 1e-9 ) );
+    EXPECT_LT( intersection.pointDistance( r1 ), 1e-12 );
+
+    // Clearly outside stays outside
+    EXPECT_EQ( 0,
+               GeometryTools::intersectLineSegmentTriangle( cvf::Vec3d( 2, 2, -1 ), cvf::Vec3d( 2, 2, 1 ), t0, t1, t2, &intersection, &isEntering, 1e-9 ) );
+}