From 4d12d0f406e1f57d47ade70615a3a3afe4970ddf Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Fri, 8 Dec 2023 16:00:20 +0100 Subject: [PATCH] Fault Reactivation: extract stress data using well log extraction. --- .../Faults/RimFaultReactivationDataAccess.cpp | 18 +- .../RimFaultReactivationDataAccessor.cpp | 5 +- .../Faults/RimFaultReactivationDataAccessor.h | 10 +- ...imFaultReactivationDataAccessorGeoMech.cpp | 4 +- .../RimFaultReactivationDataAccessorGeoMech.h | 5 +- ...ltReactivationDataAccessorPorePressure.cpp | 3 +- ...aultReactivationDataAccessorPorePressure.h | 5 +- ...RimFaultReactivationDataAccessorStress.cpp | 307 ++++++++++++++++-- .../RimFaultReactivationDataAccessorStress.h | 57 +++- ...ultReactivationDataAccessorTemperature.cpp | 3 +- ...FaultReactivationDataAccessorTemperature.h | 5 +- ...FaultReactivationDataAccessorVoidRatio.cpp | 3 +- ...imFaultReactivationDataAccessorVoidRatio.h | 5 +- .../Faults/RimFaultReactivationModel.cpp | 1 + 14 files changed, 376 insertions(+), 55 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp index 725eab421f..f2e7691c0f 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp @@ -71,7 +71,7 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( RimEclipseCase* RimFaultReactivation::Property::LateralStressComponentY }; for ( auto property : stressProperties ) { - m_accessors.push_back( std::make_shared( geoMechCase, property ) ); + m_accessors.push_back( std::make_shared( geoMechCase, property, 0.003 ) ); } } } @@ -126,12 +126,18 @@ std::vector RimFaultReactivationDataAccess::extractModelData( const RigF }; std::shared_ptr accessor = getAccessor( property ); + if ( accessor ) { - accessor->setTimeStep( timeStep ); + accessor->setModelAndTimeStep( model, timeStep ); auto grid = model.grid( gridPart ); + const std::map>& borderSurfaceElements = grid->borderSurfaceElements(); + auto it = borderSurfaceElements.find( RimFaultReactivation::BorderSurface::Seabed ); + CAF_ASSERT( it != borderSurfaceElements.end() && "Sea bed border surface does not exist" ); + std::set seabedElements( it->second.begin(), it->second.end() ); + std::vector values; if ( nodeProperties.contains( property ) ) @@ -149,11 +155,15 @@ std::vector RimFaultReactivationDataAccess::extractModelData( const RigF { std::vector corners = grid->elementCorners( elementIndex ); - double topDepth = computeAverageDepth( corners, { 0, 1, 2, 3 } ); + // Move top of sea bed element down to end up inside top element + bool isTopElement = seabedElements.contains( static_cast( elementIndex ) ); + double topDepthAdjust = isTopElement ? 0.1 : 0.0; + + double topDepth = computeAverageDepth( corners, { 0, 1, 2, 3 } ) - topDepthAdjust; double bottomDepth = computeAverageDepth( corners, { 4, 5, 6, 7 } ); cvf::Vec3d position = RigCaseToCaseCellMapperTools::calculateCellCenter( corners.data() ); - double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth ); + double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth, elementIndex ); values.push_back( value ); } } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.cpp index c12f04b509..209b3b2bff 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.cpp @@ -18,6 +18,7 @@ #include "RimFaultReactivationDataAccessor.h" +#include "RigFaultReactivationModel.h" #include "RigMainGrid.h" #include "cafAssert.h" @@ -28,6 +29,7 @@ RimFaultReactivationDataAccessor::RimFaultReactivationDataAccessor() { m_timeStep = -1; + m_model = nullptr; } //-------------------------------------------------------------------------------------------------- @@ -40,8 +42,9 @@ RimFaultReactivationDataAccessor::~RimFaultReactivationDataAccessor() //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RimFaultReactivationDataAccessor::setTimeStep( size_t timeStep ) +void RimFaultReactivationDataAccessor::setModelAndTimeStep( const RigFaultReactivationModel& model, size_t timeStep ) { + m_model = &model; m_timeStep = timeStep; updateResultAccessor(); } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h index 3d3d3504fc..27fe4c14a6 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h @@ -38,18 +38,20 @@ class RimFaultReactivationDataAccessor RimFaultReactivationDataAccessor(); ~RimFaultReactivationDataAccessor(); - virtual void setTimeStep( size_t timeStep ); + virtual void setModelAndTimeStep( const RigFaultReactivationModel& model, size_t timeStep ); virtual bool isMatching( RimFaultReactivation::Property property ) const = 0; virtual double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const = 0; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const = 0; protected: virtual void updateResultAccessor() = 0; - size_t m_timeStep; + const RigFaultReactivationModel* m_model; + size_t m_timeStep; }; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp index 7da47f84c5..4c81dfd05b 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp @@ -101,8 +101,8 @@ double RimFaultReactivationDataAccessorGeoMech::valueAtPosition( const cvf::Vec3 const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const - + double bottomDepth, + size_t elementIndex ) const { if ( !m_resultFrames ) return std::numeric_limits::infinity(); diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h index e89241d175..5e7f9bdb06 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h @@ -43,8 +43,9 @@ class RimFaultReactivationDataAccessorGeoMech : public RimFaultReactivationDataA double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp index 9d32921eb2..e9a386b1e9 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp @@ -90,7 +90,8 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf: const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h index 6ae024a7a0..647537a7cb 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivation double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp index e14c1d7bee..b5d843e36b 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp @@ -19,19 +19,28 @@ #include "RimFaultReactivationDataAccessorStress.h" #include "RiaEclipseUnitTools.h" +#include "RiaInterpolationTools.h" +#include "RiaLogging.h" +#include "RigFaultReactivationModel.h" #include "RigFemAddressDefines.h" #include "RigFemPartCollection.h" #include "RigFemPartResultsCollection.h" #include "RigFemResultAddress.h" #include "RigFemScalarResultFrames.h" #include "RigGeoMechCaseData.h" +#include "RigGeoMechWellLogExtractor.h" +#include "RigGriddedPart3d.h" #include "RigResultAccessorFactory.h" +#include "RigWellPath.h" #include "RimFaultReactivationEnums.h" +#include "RimFracture.h" #include "RimGeoMechCase.h" #include "RimWellIADataAccess.h" +#include "cvfVector3.h" + #include #include @@ -39,9 +48,11 @@ /// //-------------------------------------------------------------------------------------------------- RimFaultReactivationDataAccessorStress::RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, - RimFaultReactivation::Property property ) + RimFaultReactivation::Property property, + double gradient ) : m_geoMechCase( geoMechCase ) , m_property( property ) + , m_gradient( gradient ) { m_geoMechCaseData = geoMechCase->geoMechData(); } @@ -60,22 +71,76 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor() { const int partIndex = 0; - auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr ) -> RigFemScalarResultFrames* + auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr, int timeStepIndex ) -> RigFemScalarResultFrames* { - auto result = femParts->findOrLoadScalarResult( partIndex, addr ); - if ( result->frameData( 0, 0 ).empty() ) + auto result = femParts->findOrLoadScalarResult( partIndex, addr ); + int frameIndex = result->frameCount( timeStepIndex ) - 1; + if ( result->frameData( timeStepIndex, frameIndex ).empty() ) { return nullptr; } return result; }; - auto femParts = m_geoMechCaseData->femPartResults(); - m_femPart = femParts->parts()->part( partIndex ); - m_s33Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S33" ) ); - m_s11Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S11" ) ); - m_s22Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ) ); - m_porFrames = loadFrameLambda( femParts, RigFemAddressDefines::elementNodalPorBarAddress() ); + auto femParts = m_geoMechCaseData->femPartResults(); + m_femPart = femParts->parts()->part( partIndex ); + int timeStepIndex = 0; + m_s33Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S33" ), timeStepIndex ); + m_s11Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S11" ), timeStepIndex ); + m_s22Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ), timeStepIndex ); + + auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom(); + auto faultNormal = m_model->faultNormal(); + double distanceFromFault = 1.0; + + RigFemPartCollection* geoMechPartCollection = m_geoMechCaseData->femParts(); + std::string errorName = "fault reactivation data access"; + + { + std::vector wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, faultNormal * distanceFromFault ); + m_faceAWellPath = new RigWellPath( wellPoints, generateMds( wellPoints ) ); + m_partIndexA = getPartIndexFromPoint( *geoMechPartCollection, wellPoints[1] ); + m_extractorA = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceAWellPath.p(), errorName ); + } + + { + std::vector wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, -faultNormal * distanceFromFault ); + m_faceBWellPath = new RigWellPath( wellPoints, generateMds( wellPoints ) ); + m_partIndexB = getPartIndexFromPoint( *geoMechPartCollection, wellPoints[1] ); + m_extractorB = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceBWellPath.p(), errorName ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +int RimFaultReactivationDataAccessorStress::getPartIndexFromPoint( const RigFemPartCollection& partCollection, const cvf::Vec3d& point ) +{ + const int idx = 0; + + // Find candidates for intersected global elements + const cvf::BoundingBox intersectingBb( point, point ); + std::vector intersectedGlobalElementIndexCandidates; + partCollection.findIntersectingGlobalElementIndices( intersectingBb, &intersectedGlobalElementIndexCandidates ); + + if ( intersectedGlobalElementIndexCandidates.empty() ) return idx; + + // Iterate through global element candidates and check if point is in hexCorners + for ( const auto& globalElementIndex : intersectedGlobalElementIndexCandidates ) + { + const auto [part, elementIndex] = partCollection.partAndElementIndex( globalElementIndex ); + + // Find nodes from element + std::array coordinates; + const bool isSuccess = part->fillElementCoordinates( elementIndex, coordinates ); + if ( !isSuccess ) continue; + + const bool isPointInCell = RigHexIntersectionTools::isPointInCell( point, coordinates.data() ); + if ( isPointInCell ) return part->elementPartId(); + } + + // Utilize first part to have an id + return idx; } //-------------------------------------------------------------------------------------------------- @@ -101,9 +166,10 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { - if ( !m_porFrames || !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits::infinity(); + if ( !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits::infinity(); RimWellIADataAccess iaDataAccess( m_geoMechCase ); int centerElementIdx = iaDataAccess.elementIndex( position ); @@ -116,23 +182,22 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d if ( centerElementIdx != -1 && topElementIdx != -1 && bottomElementIdx != -1 ) { int timeStepIndex = 0; - int frameIndex = 0; + int frameIndex = m_s33Frames->frameCount( timeStepIndex ) - 1; - const std::vector& s11Data = m_s11Frames->frameData( timeStepIndex, frameIndex ); - const std::vector& s22Data = m_s22Frames->frameData( timeStepIndex, frameIndex ); const std::vector& s33Data = m_s33Frames->frameData( timeStepIndex, frameIndex ); - const std::vector& porData = m_porFrames->frameData( timeStepIndex, frameIndex ); if ( m_property == RimFaultReactivation::Property::StressTop ) { - double s33 = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, porData ); + auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, topPosition, m_gradient, timeStepIndex, frameIndex ); + if ( std::isinf( porBar ) ) return porBar; + double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data ); return RiaEclipseUnitTools::barToPascal( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::StressBottom ) { - double s33 = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, porData ); + auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, bottomPosition, m_gradient, timeStepIndex, frameIndex ); + if ( std::isinf( porBar ) ) return porBar; + double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data ); return RiaEclipseUnitTools::barToPascal( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::DepthTop ) @@ -145,16 +210,22 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d } else if ( m_property == RimFaultReactivation::Property::LateralStressComponentX ) { - double s11 = interpolatedResultValue( iaDataAccess, m_femPart, position, s11Data ); - double s33 = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData ); + auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, position, m_gradient, timeStepIndex, frameIndex ); + if ( std::isinf( porBar ) ) return porBar; + + const std::vector& s11Data = m_s11Frames->frameData( timeStepIndex, frameIndex ); + double s11 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s11Data ); + double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data ); return ( s11 - porBar ) / ( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY ) { - double s22 = interpolatedResultValue( iaDataAccess, m_femPart, position, s22Data ); - double s33 = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData ); + auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, position, m_gradient, timeStepIndex, frameIndex ); + if ( std::isinf( porBar ) ) return porBar; + + const std::vector& s22Data = m_s22Frames->frameData( timeStepIndex, frameIndex ); + double s22 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s22Data ); + double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data ); return ( s22 - porBar ) / ( s33 - porBar ); } } @@ -172,3 +243,187 @@ double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellI { return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position ); } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStress::getPorBar( RimWellIADataAccess& iaDataAccess, + const RigFemPart* femPart, + const cvf::Vec3d& position, + double gradient, + int timeStepIndex, + int frameIndex ) const +{ + RigFemPartCollection* partCollection = m_geoMechCaseData->femParts(); + cvf::ref extractor = m_partIndexA == getPartIndexFromPoint( *partCollection, position ) ? m_extractorA + : m_extractorB; + if ( !extractor->valid() ) + { + RiaLogging::error( "Invalid extractor when extracting PorBar" ); + return { std::numeric_limits::infinity(), cvf::Vec3d::UNDEFINED }; + } + + RigFemResultAddress resAddr = RigFemAddressDefines::nodalPorBarAddress(); + std::vector values; + extractor->curveData( resAddr, timeStepIndex, frameIndex, &values ); + + // Fill in missing values + auto intersections = extractor->intersections(); + fillInMissingValues( intersections, values, gradient ); + + // Linear interpolation between two points + auto lerp = []( const cvf::Vec3d& start, const cvf::Vec3d& end, double t ) { return start + t * ( end - start ); }; + + auto [topIdx, bottomIdx] = findIntersectionsForTvd( intersections, position.z() ); + if ( topIdx != -1 && bottomIdx != -1 ) + { + double topValue = values[topIdx]; + double bottomValue = values[bottomIdx]; + if ( !std::isinf( topValue ) && !std::isinf( bottomValue ) ) + { + // Interpolate value from the two closest points. + std::vector xs = { intersections[bottomIdx].z(), intersections[topIdx].z() }; + std::vector ys = { values[bottomIdx], values[topIdx] }; + double porBar = RiaInterpolationTools::linear( xs, ys, position.z() ); + + // Interpolate position from depth + double fraction = RiaInterpolationTools::linear( xs, { 0.0, 1.0 }, position.z() ); + cvf::Vec3d extractionPosition = lerp( intersections[bottomIdx], intersections[topIdx], fraction ); + return { porBar, extractionPosition }; + } + } + + return { std::numeric_limits::infinity(), cvf::Vec3d::UNDEFINED }; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStress::findIntersectionsForTvd( const std::vector& intersections, double tvd ) +{ + int topIdx = -1; + int bottomIdx = -1; + if ( intersections.size() >= 2 ) + { + for ( size_t i = 1; i < intersections.size(); i++ ) + { + auto top = intersections[i - 1]; + auto bottom = intersections[i]; + if ( top.z() > tvd && bottom.z() < tvd ) + { + topIdx = static_cast( i ) - 1; + bottomIdx = static_cast( i ); + break; + } + } + } + return std::make_pair( topIdx, bottomIdx ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStress::findOverburdenAndUnderburdenIndex( const std::vector& values ) +{ + auto findLastOverburdenIndex = []( const std::vector& values ) + { + for ( size_t i = 0; i < values.size(); i++ ) + { + if ( !std::isinf( values[i] ) ) return static_cast( i ); + } + + return -1; + }; + + auto findFirstUnderburdenIndex = []( const std::vector& values ) + { + for ( size_t i = values.size() - 1; i > 0; i-- ) + { + if ( !std::isinf( values[i] ) ) return static_cast( i ); + } + + return -1; + }; + + int lastOverburdenIndex = findLastOverburdenIndex( values ); + int firstUnderburdenIndex = findFirstUnderburdenIndex( values ); + return { lastOverburdenIndex, firstUnderburdenIndex }; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RimFaultReactivationDataAccessorStress::computePorBarWithGradient( const std::vector& intersections, + const std::vector& values, + int i1, + int i2, + double gradient ) +{ + double tvdDiff = intersections[i2].z() - intersections[i1].z(); + return tvdDiff * gradient + values[i2]; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFaultReactivationDataAccessorStress::fillInMissingValues( const std::vector& intersections, + std::vector& values, + double gradient ) +{ + CAF_ASSERT( intersections.size() == values.size() ); + + auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); + + // Fill in overburden values using gradient + for ( int i = 0; i < lastOverburdenIndex; i++ ) + { + values[i] = computePorBarWithGradient( intersections, values, i, lastOverburdenIndex, gradient ); + } + + // Fill in underburden values using gradient + int lastElementIndex = static_cast( values.size() ) - 1; + for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- ) + { + values[i] = computePorBarWithGradient( intersections, values, i, firstUnderburdenIndex, gradient ); + } + + std::vector intersectionsZ; + for ( auto i : intersections ) + { + intersectionsZ.push_back( i.z() ); + } + + RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimFaultReactivationDataAccessorStress::generateMds( const std::vector& points ) +{ + CAF_ASSERT( points.size() >= 2 ); + + // Assume first at zero, all other points relative to that. + std::vector mds = { 0.0 }; + double sum = 0.0; + for ( size_t i = 1; i < points.size(); i++ ) + { + sum += points[i - 1].pointDistance( points[i] ); + mds.push_back( sum ); + } + return mds; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimFaultReactivationDataAccessorStress::generateWellPoints( const cvf::Vec3d& faultTopPosition, + const cvf::Vec3d& faultBottomPosition, + const cvf::Vec3d& offset ) +{ + cvf::Vec3d faultTop = faultTopPosition + offset; + cvf::Vec3d seabed( faultTop.x(), faultTop.y(), 0.0 ); + cvf::Vec3d faultBottom = faultBottomPosition + offset; + cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), -10000.0 ); + return { seabed, faultTop, faultBottom, underburdenBottom }; +} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h index 4fdcfa4bd2..cf667f14e2 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h @@ -25,12 +25,22 @@ #include +#include "cafPdmField.h" +#include "cafPdmPtrField.h" + +#include "cvfObject.h" + +class RigFemPart; class RigFemPartResultsCollection; -class RimGeoMechCase; -class RigGeoMechCaseData; class RigFemScalarResultFrames; -class RigFemPart; +class RigGeoMechCaseData; +class RigGriddedPart3d; +class RimGeoMechCase; class RimWellIADataAccess; +class RimModeledWellPath; +class RigWellPath; +class RigFemPartCollection; +class RigGeoMechWellLogExtractor; //================================================================================================== /// @@ -39,7 +49,7 @@ class RimWellIADataAccess; class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAccessor { public: - RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property ); + RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property, double gradient ); ~RimFaultReactivationDataAccessorStress(); bool isMatching( RimFaultReactivation::Property property ) const override; @@ -47,8 +57,9 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; @@ -60,12 +71,44 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc const cvf::Vec3d& position, const std::vector& scalarResults ) const; + std::pair getPorBar( RimWellIADataAccess& iaDataAccess, + const RigFemPart* femPart, + const cvf::Vec3d& position, + double gradient, + int timeStepIndex, + int frameIndex ) const; + + static std::pair + findElementSetContainingElement( const std::map>& elementSets, + unsigned int elmIdx ); + + static int getPartIndexFromPoint( const RigFemPartCollection& partCollection, const cvf::Vec3d& point ); + static std::pair findIntersectionsForTvd( const std::vector& intersections, double tvd ); + static std::pair findOverburdenAndUnderburdenIndex( const std::vector& values ); + static double computePorBarWithGradient( const std::vector& intersections, + const std::vector& values, + int i1, + int i2, + double gradient ); + static void fillInMissingValues( const std::vector& intersections, std::vector& values, double gradient ); + + static std::vector generateMds( const std::vector& points ); + static std::vector + generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset ); + RimGeoMechCase* m_geoMechCase; RimFaultReactivation::Property m_property; + double m_gradient; RigGeoMechCaseData* m_geoMechCaseData; RigFemScalarResultFrames* m_s11Frames; RigFemScalarResultFrames* m_s22Frames; RigFemScalarResultFrames* m_s33Frames; - RigFemScalarResultFrames* m_porFrames; const RigFemPart* m_femPart; + + cvf::ref m_faceAWellPath; + cvf::ref m_faceBWellPath; + int m_partIndexA; + int m_partIndexB; + cvf::ref m_extractorA; + cvf::ref m_extractorB; }; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp index 10607338c2..52c5eef8bf 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp @@ -87,7 +87,8 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf:: const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h index 8cb50c065c..3d310152d6 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorTemperature : public RimFaultReactivationD double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp index c1dc6ccfb3..96fe0c6f3b 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp @@ -88,7 +88,8 @@ double RimFaultReactivationDataAccessorVoidRatio::valueAtPosition( const cvf::Ve const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h index 069b3b344e..519a5727da 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorVoidRatio : public RimFaultReactivationDat double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp index 0ded3a19f1..46bcd52f72 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp @@ -709,6 +709,7 @@ bool RimFaultReactivationModel::extractAndExportModelData() // extract data for each timestep m_dataAccess = std::make_shared( eCase, geoMechCase(), selectedTimeStepIndexes ); m_dataAccess->extractModelData( *model() ); + m_dataAccess.reset(); return true; }