From 577b4e770e313ef4c8f8fc301240f9d3e199ca75 Mon Sep 17 00:00:00 2001 From: Stefan Crooijmans Date: Mon, 2 Nov 2026 21:21:46 +0100 Subject: [PATCH] #15088 Well Log Extractor: Use parametric tolerance for well path cell intersections --- .github/workflows/ResInsightMac.yml | 5 +- .../RigHexIntersectionTools.cpp | 6 +- .../RigHexIntersectionTools.h | 3 +- .../Well/RigEclipseWellLogExtractor.cpp | 7 +- .../Well/RigGeoMechWellLogExtractor.cpp | 9 +- .../Well/RigWellLogExtractionTools.h | 7 + .../ReservoirDataModel/cvfGeometryTools.cpp | 18 +- .../ReservoirDataModel/cvfGeometryTools.h | 3 +- .../RigHexIntersectionTools-Test.cpp | 170 ++++++++++++++++++ .../UnitTests/cvfGeometryTools-Test.cpp | 40 +++++ 10 files changed, 252 insertions(+), 16 deletions(-) diff --git a/.github/workflows/ResInsightMac.yml b/.github/workflows/ResInsightMac.yml index 8d01822ab2..56d2fa5421 100644 --- a/.github/workflows/ResInsightMac.yml +++ b/.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: | diff --git a/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.cpp b/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.cpp index dd8d9ec58d..50742b7a14 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.cpp @@ -33,7 +33,8 @@ int RigHexIntersectionTools::lineHexCellIntersection( const cvf::Vec3d& const cvf::Vec3d& p2, const std::array& hexCorners, const size_t hexIndex, - std::vector* intersections ) + std::vector* 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( diff --git a/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.h b/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.h index 67f2f4746c..c63a852bb4 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.h +++ b/ApplicationLibCode/ReservoirDataModel/RigHexIntersectionTools.h @@ -56,7 +56,8 @@ struct RigHexIntersectionTools const cvf::Vec3d& p2, const std::array& hexCorners, const size_t hexIndex, - std::vector* intersections ); + std::vector* intersections, + double parametricTolerance = 0.0 ); static bool lineIntersectsHexCell( const cvf::Vec3d& p1, const cvf::Vec3d& p2, const std::array& hexCorners ); diff --git a/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp b/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp index d278970963..190ade40af 100644 --- a/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp +++ b/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp @@ -89,7 +89,12 @@ void RigEclipseWellLogExtractor::calculateIntersection() if ( cell.isInvalid() || cell.subGrid() != nullptr ) continue; std::array hexCorners = m_caseData->mainGrid()->cellCornerVertices( globalCellIndex ); - RigHexIntersectionTools::lineHexCellIntersection( p1, p2, hexCorners, globalCellIndex, &intersections ); + RigHexIntersectionTools::lineHexCellIntersection( p1, + p2, + hexCorners, + globalCellIndex, + &intersections, + RigWellLogExtractionTools::intersectionParametricTolerance ); } if ( !isCellFaceNormalsOut ) diff --git a/ApplicationLibCode/ReservoirDataModel/Well/RigGeoMechWellLogExtractor.cpp b/ApplicationLibCode/ReservoirDataModel/Well/RigGeoMechWellLogExtractor.cpp index ca3eee54dc..8dd650c887 100644 --- a/ApplicationLibCode/ReservoirDataModel/Well/RigGeoMechWellLogExtractor.cpp +++ b/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 diff --git a/ApplicationLibCode/ReservoirDataModel/Well/RigWellLogExtractionTools.h b/ApplicationLibCode/ReservoirDataModel/Well/RigWellLogExtractionTools.h index 9e1922f0e0..2f2e665c68 100644 --- a/ApplicationLibCode/ReservoirDataModel/Well/RigWellLogExtractionTools.h +++ b/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; diff --git a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp index 83d92c18db..b82c490cfd 100644 --- a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp +++ b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp @@ -21,6 +21,8 @@ #include "cafAssert.h" +#include + 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 diff --git a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h b/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h index e539c0e534..c48b17f35e 100644 --- a/ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h +++ b/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 ); diff --git a/ApplicationLibCode/UnitTests/RigHexIntersectionTools-Test.cpp b/ApplicationLibCode/UnitTests/RigHexIntersectionTools-Test.cpp index 15e0cd7e2d..68fbce5125 100644 --- a/ApplicationLibCode/UnitTests/RigHexIntersectionTools-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigHexIntersectionTools-Test.cpp @@ -23,6 +23,13 @@ #include #include "RigHexIntersectionTools.h" +#include "Well/RigWellLogExtractionTools.h" + +#include "cvfGeometryTools.h" +#include "cvfStructGrid.h" + +#include +#include //-------------------------------------------------------------------------------------------------- /// @@ -97,3 +104,166 @@ TEST( RigHexIntersectionTools, DISABLED_planeHexCellIntersectionPerformanceTest qDebug() << "Time total elapsed: " << timeTotal.elapsed() << " intersected: " << static_cast( isCellIntersected ); } + +//-------------------------------------------------------------------------------------------------- +/// Hex corners in ResInsight corner order, slightly distorted, at UTM-like coordinates +//-------------------------------------------------------------------------------------------------- +static std::array 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 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 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 easting( 400000.0, 600000.0 ); + std::uniform_real_distribution northing( 6000000.0, 7000000.0 ); + std::uniform_real_distribution depth( -4000.0, -1000.0 ); + + int missed = 0; + + for ( int hexIdx = 0; hexIdx < hexCount; ++hexIdx ) + { + std::array 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( 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 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( 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 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 exact; + std::vector 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 exact; + std::vector 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 exact; + std::vector 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 ); + } +} diff --git a/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp b/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp index 58b064d974..9b6359e0b3 100644 --- a/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp +++ b/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 ) ); +}