From a908b0c9cfb463d78d475796bf33766b97826123 Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Thu, 30 Jun 2022 15:29:42 +0200 Subject: [PATCH] #9098 Thermal Fracture: interpolate missing values inside convex hull of points. --- .../ReservoirDataModel/CMakeLists_files.cmake | 2 + .../ReservoirDataModel/RigConvexHull.cpp | 89 +++++++++++++++++++ .../ReservoirDataModel/RigConvexHull.h | 37 ++++++++ .../RigThermalFractureResultUtil.cpp | 84 +++++++++++------ .../RigThermalFractureResultUtil.h | 6 ++ .../UnitTests/CMakeLists_files.cmake | 1 + .../UnitTests/RigConvexHull-Test.cpp | 56 ++++++++++++ 7 files changed, 246 insertions(+), 29 deletions(-) create mode 100644 ApplicationLibCode/ReservoirDataModel/RigConvexHull.cpp create mode 100644 ApplicationLibCode/ReservoirDataModel/RigConvexHull.h create mode 100644 ApplicationLibCode/UnitTests/RigConvexHull-Test.cpp diff --git a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake index f66ae9f408..d732b87996 100644 --- a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake @@ -64,6 +64,7 @@ set(SOURCE_GROUP_HEADER_FILES ${CMAKE_CURRENT_LIST_DIR}/RigThermalFractureResultUtil.h ${CMAKE_CURRENT_LIST_DIR}/RigFractureGrid.h ${CMAKE_CURRENT_LIST_DIR}/RigFractureCell.h + ${CMAKE_CURRENT_LIST_DIR}/RigConvexHull.h ${CMAKE_CURRENT_LIST_DIR}/RigWellResultPoint.h ${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryTools.h ${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryExporter.h @@ -155,6 +156,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RigThermalFractureResultUtil.cpp ${CMAKE_CURRENT_LIST_DIR}/RigFractureGrid.cpp ${CMAKE_CURRENT_LIST_DIR}/RigFractureCell.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigConvexHull.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWellResultPoint.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryTools.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryExporter.cpp diff --git a/ApplicationLibCode/ReservoirDataModel/RigConvexHull.cpp b/ApplicationLibCode/ReservoirDataModel/RigConvexHull.cpp new file mode 100644 index 0000000000..88e21c7da2 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigConvexHull.cpp @@ -0,0 +1,89 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2022 - Equinor ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigConvexHull.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigConvexHull::compute2d( const std::vector& points ) +{ + // Sort our points from left to right + std::vector sortedPoints = sortPoints( points ); + + // Find the lower half of the convex hull. + std::vector lower; + for ( auto it = sortedPoints.begin(); it != sortedPoints.end(); ++it ) + { + // Remove any points that does not make a convex angle with the current point + removePointsWithoutConvexAngle( lower, *it ); + lower.push_back( *it ); + } + + // Find the upper half of the convex hull. + std::vector upper; + for ( auto it = sortedPoints.rbegin(); it != sortedPoints.rend(); ++it ) + { + // Remove any points that does not make a convex angle with the current point + removePointsWithoutConvexAngle( upper, *it ); + upper.push_back( *it ); + } + + // Concatenation of the lower and upper hulls gives the convex hull. + // Last point of each list is omitted because it is repeated at the beginning of the other list. + std::vector hull; + hull.insert( hull.end(), lower.begin(), lower.end() - 1 ); + hull.insert( hull.end(), upper.begin(), upper.end() - 1 ); + return hull; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigConvexHull::sortPoints( const std::vector& unsorted ) +{ + // Returns true if a is left of b + auto isLeftOf = []( const cvf::Vec3d& a, const cvf::Vec3d& b ) { + return ( a.x() < b.x() || ( a.x() == b.x() && a.y() < b.y() ) ); + }; + + std::vector sorted = unsorted; + std::sort( sorted.begin(), sorted.end(), isLeftOf ); + return sorted; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigConvexHull::removePointsWithoutConvexAngle( std::vector& points, const cvf::Vec3d& current ) +{ + // 2D cross product of (a, b) and (a, c) vectors, i.e. z-component of their 3D cross product. + // Returns a positive value, if c->(a,b) makes a counter-clockwise turn, + // negative for clockwise turn, and zero if the points are collinear. + auto counterClockWise = []( const cvf::Vec3d& a, const cvf::Vec3d& b, const cvf::Vec3d& c ) { + return ( b.x() - a.x() ) * ( c.y() - a.y() ) - ( b.y() - a.y() ) * ( c.x() - a.x() ); + }; + + // Remove all points that is not a counter-clockwise turn from current point + while ( points.size() >= 2 && counterClockWise( *( points.rbegin() + 1 ), *( points.rbegin() ), current ) >= 0 ) + { + points.pop_back(); + } +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigConvexHull.h b/ApplicationLibCode/ReservoirDataModel/RigConvexHull.h new file mode 100644 index 0000000000..c3804b5733 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigConvexHull.h @@ -0,0 +1,37 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2022 - Equinor ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "cvfVector3.h" + +#include + +//================================================================================================== +/// +/// +//================================================================================================== +class RigConvexHull +{ +public: + static std::vector compute2d( const std::vector& points ); + +private: + static std::vector sortPoints( const std::vector& unsorted ); + static void removePointsWithoutConvexAngle( std::vector& points, const cvf::Vec3d& current ); +}; diff --git a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp index 3543fe2d41..87df13fa9d 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp @@ -20,7 +20,9 @@ #include "RiaLogging.h" +#include "RiaWeightedMeanCalculator.h" #include "RigCellGeometryTools.h" +#include "RigConvexHull.h" #include "RigFractureGrid.h" #include "RigStatisticsMath.h" #include "RigThermalFractureDefinition.h" @@ -87,40 +89,27 @@ std::vector> vec.push_back( junk ); } - // Loop through each polygon and check if point is inside. - for ( size_t nodeIndex = 0; nodeIndex < fractureDefinition->numNodes(); nodeIndex++ ) + // Find the boundary of the fracture (i.e. convex hull around the points) + std::vector fractureBoundary = RigConvexHull::compute2d( relativePos ); + fractureBoundary.push_back( fractureBoundary.front() ); + + for ( int i = 0; i < static_cast( xCoords.size() ) - 1; i++ ) { - bool placed = false; - - for ( int i = 0; i < static_cast( xCoords.size() ) - 1; i++ ) + for ( int j = 0; j < static_cast( depthCoords.size() ) - 1; j++ ) { - for ( int j = 0; j < static_cast( depthCoords.size() ) - 1; j++ ) + cvf::BoundingBox bb; + bb.add( cvf::Vec3d( xCoords[i], depthCoords[j], 0.0 ) ); + bb.add( cvf::Vec3d( xCoords[i + 1], depthCoords[j], 0.0 ) ); + bb.add( cvf::Vec3d( xCoords[i + 1], depthCoords[j + 1], 0.0 ) ); + bb.add( cvf::Vec3d( xCoords[i], depthCoords[j + 1], 0.0 ) ); + + if ( RigCellGeometryTools::pointInsidePolygon2D( bb.center(), fractureBoundary ) ) { - std::vector cellPolygon; - cellPolygon.push_back( cvf::Vec3d( xCoords[i], depthCoords[j], 0.0 ) ); - cellPolygon.push_back( cvf::Vec3d( xCoords[i + 1], depthCoords[j], 0.0 ) ); - cellPolygon.push_back( cvf::Vec3d( xCoords[i + 1], depthCoords[j + 1], 0.0 ) ); - cellPolygon.push_back( cvf::Vec3d( xCoords[i], depthCoords[j + 1], 0.0 ) ); - // First and last polygon must match - cellPolygon.push_back( cvf::Vec3d( xCoords[i], depthCoords[j], 0.0 ) ); - - if ( RigCellGeometryTools::pointInsidePolygon2D( relativePos[nodeIndex], cellPolygon ) ) - { - double value = fractureDefinition->getPropertyValue( propertyIndex, nodeIndex, timeStepIndex ); - vec[j][i] = value; - placed = true; - } + double value = + interpolateProperty( bb.center(), relativePos, fractureDefinition, propertyIndex, timeStepIndex ); + vec[j][i] = value; } } - - if ( !placed ) - { - RiaLogging::error( QString( "Unable to place point: %1 [%2 %3 %4]" ) - .arg( nodeIndex ) - .arg( relativePos[nodeIndex].x() ) - .arg( relativePos[nodeIndex].y() ) - .arg( relativePos[nodeIndex].z() ) ); - } } return vec; @@ -443,7 +432,44 @@ std::vector for ( auto& r : relativePos ) { r.transformVector( rotMat ); + r.z() = 0.0; } return relativePos; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigThermalFractureResultUtil::interpolateProperty( const cvf::Vec3d& position, + const std::vector& points, + std::shared_ptr fractureDefinition, + int propertyIndex, + size_t timeStepIndex ) +{ + // Compute the distance to the other points + std::vector> distances; + + int index = 0; + for ( auto p : points ) + { + double distance = position.pointDistance( p ); + distances.push_back( std::make_pair( distance, index ) ); + index++; + } + + // Sort by distance + std::sort( distances.begin(), distances.end() ); + + // Create distance-weighthed mean of first few points + size_t numPoints = 3; + RiaWeightedMeanCalculator calc; + for ( size_t i = 0; i < numPoints; i++ ) + { + auto [distance, nodeIndex] = distances[i]; + double value = fractureDefinition->getPropertyValue( propertyIndex, nodeIndex, timeStepIndex ); + calc.addValueAndWeight( value, distance ); + } + + return calc.weightedMean(); +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h index 7333fc6db7..cc71bf9ca5 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h +++ b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h @@ -94,4 +94,10 @@ private: static std::vector getRelativeCoordinates( std::shared_ptr fractureDefinition, size_t timeStepIndex ); + + static double interpolateProperty( const cvf::Vec3d& position, + const std::vector& points, + std::shared_ptr fractureDefinition, + int propertyIndex, + size_t timeStepIndex ); }; diff --git a/ApplicationLibCode/UnitTests/CMakeLists_files.cmake b/ApplicationLibCode/UnitTests/CMakeLists_files.cmake index b6cb8bb1b6..c08b76929a 100644 --- a/ApplicationLibCode/UnitTests/CMakeLists_files.cmake +++ b/ApplicationLibCode/UnitTests/CMakeLists_files.cmake @@ -77,6 +77,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RifStimPlanModelDeviationFrkExporter-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RifSummaryDataReader-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigSlice2D-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigConvexHull-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigSurfaceResampler-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigSurfaceStatisticsCalculator-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/StructGridInterface-Test.cpp diff --git a/ApplicationLibCode/UnitTests/RigConvexHull-Test.cpp b/ApplicationLibCode/UnitTests/RigConvexHull-Test.cpp new file mode 100644 index 0000000000..f78989703f --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigConvexHull-Test.cpp @@ -0,0 +1,56 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2022 Equinor ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RigConvexHull.h" + +#include "cvfVector3.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigConvexHullTests, SimpleExample ) +{ + std::vector points = { cvf::Vec3d( 0.0, 0.0, 0.0 ), + cvf::Vec3d( 0.0, 1.0, 0.0 ), + cvf::Vec3d( 1.0, 1.0, 0.0 ), + cvf::Vec3d( 1.0, 0.0, 0.0 ), + cvf::Vec3d( 1.0, 1.0, 0.0 ), + cvf::Vec3d( 0.5, 0.5, 0.0 ), + cvf::Vec3d( 0.25, 0.25, 0.0 ), + cvf::Vec3d( 0.5, 1.25, 0.0 ), + cvf::Vec3d( 0.5, 0.75, 0.0 ) }; + + std::vector convexHull = RigConvexHull::compute2d( points ); + + std::vector expectedPoints = { cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 0, 1, 0 ), + cvf::Vec3d( 0.5, 1.25, 0 ), + cvf::Vec3d( 1, 1, 0 ), + cvf::Vec3d( 1, 0, 0 ) }; + + EXPECT_EQ( 5u, convexHull.size() ); + for ( size_t i = 0; i < convexHull.size(); i++ ) + { + EXPECT_EQ( expectedPoints[i].x(), convexHull[i].x() ); + EXPECT_EQ( expectedPoints[i].y(), convexHull[i].y() ); + } +}