From 2b17c9d9afbe2d6768b3adc98e3b24c288d72f1f Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Wed, 29 Jun 2022 14:47:28 +0200 Subject: [PATCH] Rotate fracture plane to align with x-y plane --- .../RigThermalFractureResultUtil.cpp | 98 +++++++++++++++---- .../RigThermalFractureResultUtil.h | 9 +- 2 files changed, 88 insertions(+), 19 deletions(-) diff --git a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp index 87c01a73d5..3543fe2d41 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp @@ -30,8 +30,8 @@ #include -int numSamplesX = 40; -int numSamplesY = 40; +int numSamplesX = 60; +int numSamplesY = 60; //-------------------------------------------------------------------------------------------------- /// @@ -58,14 +58,15 @@ std::vector> { int propertyIndex = fractureDefinition->getPropertyIndex( resultName ); - std::vector relativePos = fractureDefinition->relativeCoordinates( timeStepIndex ); - CAF_ASSERT( relativePos.size() == fractureDefinition->numNodes() ); + std::vector relativePos = getRelativeCoordinates( fractureDefinition, timeStepIndex ); // Create bounding box cvf::BoundingBox boundingBox; for ( auto p : relativePos ) boundingBox.add( p ); + boundingBox.expand( 1.0 ); + // Generate a uniform mesh auto [xCoordsAtNodes, yCoordsAtNodes] = generateUniformMesh( boundingBox, numSamplesX, numSamplesY ); @@ -82,31 +83,44 @@ std::vector> // Fill with invalid value for ( int i = 0; i < numSamplesY; i++ ) { - std::vector junk( numSamplesX, HUGE_VAL ); + std::vector junk( numSamplesX, std::numeric_limits::infinity() ); vec.push_back( junk ); } // Loop through each polygon and check if point is inside. - for ( int i = 0; i < static_cast( xCoords.size() ) - 1; i++ ) + for ( size_t nodeIndex = 0; nodeIndex < fractureDefinition->numNodes(); nodeIndex++ ) { - for ( int j = 0; j < static_cast( depthCoords.size() ) - 1; j++ ) - { - 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 ) ); + bool placed = false; - for ( size_t nodeIndex = 0; nodeIndex < fractureDefinition->numNodes(); nodeIndex++ ) + for ( int i = 0; i < static_cast( xCoords.size() ) - 1; i++ ) + { + for ( int j = 0; j < static_cast( depthCoords.size() ) - 1; j++ ) { + 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; + vec[j][i] = value; + placed = true; } } } + + 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; @@ -125,13 +139,15 @@ void RigThermalFractureResultUtil::createFractureTriangleGeometry( std::vector* triangleIndices ) { // Convert to coordinates relative to center node - std::vector points = fractureDefinition->relativeCoordinates( activeTimeStepIndex ); + std::vector points = getRelativeCoordinates( fractureDefinition, activeTimeStepIndex ); // Create bounding box cvf::BoundingBox boundingBox; for ( auto p : points ) boundingBox.add( p ); + boundingBox.expand( 1.0 ); + // Generate a uniform mesh auto [xCoords, depthCoords] = generateUniformMesh( boundingBox, numSamplesX, numSamplesY ); @@ -222,7 +238,7 @@ cvf::cref RiaDefines::EclipseUnitSystem requiredUnitSet ) { // Convert to coordinates relative to center node - std::vector points = fractureDefinition->relativeCoordinates( activeTimeStepIndex ); + std::vector points = getRelativeCoordinates( fractureDefinition, activeTimeStepIndex ); QString conductivityUnitTextOnFile; @@ -255,6 +271,8 @@ cvf::cref for ( auto p : points ) boundingBox.add( p ); + boundingBox.expand( 1.0 ); + // Generate a uniform mesh auto [xCoordsAtNodes, yCoordsAtNodes] = generateUniformMesh( boundingBox, numSamplesX, numSamplesY ); @@ -385,3 +403,47 @@ void RigThermalFractureResultUtil::appendDataToResultStatistics( std::shared_ptr } } } + +//-------------------------------------------------------------------------------------------------- +/// Taken from OverlayNavigationCube::computeNewUpVector +/// Consider move to geometry util class +//-------------------------------------------------------------------------------------------------- +cvf::Mat4d RigThermalFractureResultUtil::rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 ) +{ + using namespace cvf; + + Vec3d rotAxis = v1 ^ v2; + rotAxis.normalize(); + + // Guard acos against out-of-domain input + const double dotProduct = Math::clamp( v1 * v2, -1.0, 1.0 ); + const double angle = Math::acos( dotProduct ); + Mat4d rotMat = Mat4d::fromRotation( rotAxis, angle ); + return rotMat; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector + RigThermalFractureResultUtil::getRelativeCoordinates( std::shared_ptr fractureDefinition, + size_t timeStepIndex ) +{ + std::vector relativePos = fractureDefinition->relativeCoordinates( timeStepIndex ); + CAF_ASSERT( relativePos.size() == fractureDefinition->numNodes() ); + + cvf::Plane plane; + plane.setFromPoints( relativePos[0], relativePos[1], relativePos[2] ); + + cvf::Vec3d planeNormal = plane.normal().getNormalized(); + RiaLogging::info( + QString( "Plane normal: [%1 %2 %3]" ).arg( planeNormal.x() ).arg( planeNormal.y() ).arg( planeNormal.z() ) ); + auto rotMat = rotationMatrixBetweenVectors( planeNormal, cvf::Vec3d::Z_AXIS.getNormalized() ); + + for ( auto& r : relativePos ) + { + r.transformVector( rotMat ); + } + + return relativePos; +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h index 51c17eb53b..7333fc6db7 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h +++ b/ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h @@ -20,6 +20,7 @@ #include "RiaDefines.h" +#include "cvfMatrix4.h" #include "cvfObject.h" #include "cvfVector3.h" @@ -37,7 +38,7 @@ class PosNegAccumulator; namespace cvf { class BoundingBox; -}; +}; // namespace cvf //-------------------------------------------------------------------------------------------------- /// @@ -87,4 +88,10 @@ private: generateUniformMesh( const cvf::BoundingBox& bb, int numSamplesX, int numSamplesY ); static double linearSampling( double minValue, double maxValue, int numSamples, std::vector& samples ); + + static cvf::Mat4d rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 ); + + static std::vector + getRelativeCoordinates( std::shared_ptr fractureDefinition, + size_t timeStepIndex ); };