diff --git a/ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake b/ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake index bf1a368e00..e819a4deab 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake +++ b/ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake @@ -11,6 +11,7 @@ set(SOURCE_GROUP_HEADER_FILES ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorTemperature.h ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorGeoMech.h ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStress.h + ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStressGeoMech.h ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorWellLogExtraction.h ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationEnums.h ) @@ -28,6 +29,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorTemperature.cpp ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorGeoMech.cpp ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStress.cpp + ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStressGeoMech.cpp ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorWellLogExtraction.cpp ) diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp index 44275676aa..9fd5d690a6 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp @@ -35,7 +35,7 @@ #include "RimFaultReactivationDataAccessor.h" #include "RimFaultReactivationDataAccessorGeoMech.h" #include "RimFaultReactivationDataAccessorPorePressure.h" -#include "RimFaultReactivationDataAccessorStress.h" +#include "RimFaultReactivationDataAccessorStressGeoMech.h" #include "RimFaultReactivationDataAccessorTemperature.h" #include "RimFaultReactivationDataAccessorVoidRatio.h" #include "RimFaultReactivationEnums.h" @@ -76,7 +76,7 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( const RimFaultRe for ( auto property : stressProperties ) { m_accessors.push_back( - std::make_shared( geoMechCase, property, porePressureGradient, seabedDepth ) ); + std::make_shared( geoMechCase, property, porePressureGradient, seabedDepth ) ); } } } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp index 9dc4b5afc2..4d4baeae42 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp @@ -30,7 +30,6 @@ #include "RigGeoMechCaseData.h" #include "RigGeoMechWellLogExtractor.h" #include "RigGriddedPart3d.h" -#include "RigResultAccessorFactory.h" #include "RigWellPath.h" #include "RimFaultReactivationDataAccessorWellLogExtraction.h" @@ -46,16 +45,13 @@ //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -RimFaultReactivationDataAccessorStress::RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, - RimFaultReactivation::Property property, +RimFaultReactivationDataAccessorStress::RimFaultReactivationDataAccessorStress( RimFaultReactivation::Property property, double gradient, double seabedDepth ) - : m_geoMechCase( geoMechCase ) - , m_property( property ) + : m_property( property ) , m_gradient( gradient ) , m_seabedDepth( seabedDepth ) { - m_geoMechCaseData = geoMechCase->geoMechData(); } //-------------------------------------------------------------------------------------------------- @@ -65,69 +61,6 @@ RimFaultReactivationDataAccessorStress::~RimFaultReactivationDataAccessorStress( { } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimFaultReactivationDataAccessorStress::updateResultAccessor() -{ - const int partIndex = 0; - - auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr, int timeStepIndex ) -> RigFemScalarResultFrames* - { - 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 ); - 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 = - RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, - faultBottomPosition, - m_seabedDepth, - faultNormal * distanceFromFault ); - m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); - m_partIndexA = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); - m_extractorA = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceAWellPath.p(), errorName ); - } - - { - std::vector wellPoints = - RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, - faultBottomPosition, - m_seabedDepth, - -faultNormal * distanceFromFault ); - m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); - m_partIndexB = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); - m_extractorB = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceBWellPath.p(), errorName ); - } -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -RigFemResultAddress RimFaultReactivationDataAccessorStress::getResultAddress( const std::string& fieldName, const std::string& componentName ) -{ - return RigFemResultAddress( RIG_ELEMENT_NODAL, fieldName, componentName ); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -146,35 +79,25 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d double bottomDepth, size_t elementIndex ) const { - if ( !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits::infinity(); - - RimWellIADataAccess iaDataAccess( m_geoMechCase ); - int centerElementIdx = iaDataAccess.elementIndex( position ); + if ( !isDataAvailable() ) return std::numeric_limits::infinity(); cvf::Vec3d topPosition( position.x(), position.y(), topDepth ); - int topElementIdx = iaDataAccess.elementIndex( topPosition ); - cvf::Vec3d bottomPosition( position.x(), position.y(), bottomDepth ); - int bottomElementIdx = iaDataAccess.elementIndex( bottomPosition ); - if ( centerElementIdx != -1 && topElementIdx != -1 && bottomElementIdx != -1 ) + + if ( isPositionValid( position, topPosition, bottomPosition ) ) { - int timeStepIndex = 0; - int frameIndex = m_s33Frames->frameCount( timeStepIndex ) - 1; - - const std::vector& s33Data = m_s33Frames->frameData( timeStepIndex, frameIndex ); - if ( m_property == RimFaultReactivation::Property::StressTop ) { - auto [porBar, extractionPos] = calculatePorBar( topPosition, m_gradient, timeStepIndex, frameIndex ); + auto [porBar, extractionPos] = calculatePorBar( topPosition, m_gradient ); if ( std::isinf( porBar ) ) return porBar; - double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data ); + double s33 = extractStressValue( StressType::S33, extractionPos ); return RiaEclipseUnitTools::barToPascal( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::StressBottom ) { - auto [porBar, extractionPos] = calculatePorBar( bottomPosition, m_gradient, timeStepIndex, frameIndex ); + auto [porBar, extractionPos] = calculatePorBar( bottomPosition, m_gradient ); if ( std::isinf( porBar ) ) return porBar; - double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data ); + double s33 = extractStressValue( StressType::S33, extractionPos ); return RiaEclipseUnitTools::barToPascal( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::DepthTop ) @@ -187,58 +110,21 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d } else if ( m_property == RimFaultReactivation::Property::LateralStressComponentX ) { - auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, timeStepIndex, frameIndex ); + auto [porBar, extractionPos] = calculatePorBar( position, m_gradient ); 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 ); + double s11 = extractStressValue( StressType::S11, extractionPos ); + double s33 = extractStressValue( StressType::S33, extractionPos ); return ( s11 - porBar ) / ( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY ) { - auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, timeStepIndex, frameIndex ); + auto [porBar, extractionPos] = calculatePorBar( position, m_gradient ); 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 ); + double s22 = extractStressValue( StressType::S22, extractionPos ); + double s33 = extractStressValue( StressType::S33, extractionPos ); return ( s22 - porBar ) / ( s33 - porBar ); } } return std::numeric_limits::infinity(); } - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellIADataAccess& iaDataAccess, - const RigFemPart* femPart, - const cvf::Vec3d& position, - const std::vector& scalarResults ) const -{ - return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::pair - RimFaultReactivationDataAccessorStress::calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const -{ - RigFemPartCollection* partCollection = m_geoMechCaseData->femParts(); - cvf::ref extractor = m_partIndexA == partCollection->getPartIndexFromPoint( 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 ); - - return RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient ); -} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h index d14c8147a6..4d43169f8b 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h @@ -21,8 +21,6 @@ #include "RimFaultReactivationDataAccessor.h" #include "RimFaultReactivationEnums.h" -#include "RigFemResultAddress.h" - #include #include "cafPdmField.h" @@ -30,17 +28,9 @@ #include "cvfObject.h" -class RigFemPart; -class RigFemPartResultsCollection; -class RigFemScalarResultFrames; -class RigGeoMechCaseData; class RigGriddedPart3d; -class RimGeoMechCase; -class RimWellIADataAccess; class RimModeledWellPath; class RigWellPath; -class RigFemPartCollection; -class RigGeoMechWellLogExtractor; //================================================================================================== /// @@ -49,11 +39,15 @@ class RigGeoMechWellLogExtractor; class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAccessor { public: - RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, - RimFaultReactivation::Property property, - double gradient, - double seabedDepth ); - ~RimFaultReactivationDataAccessorStress(); + enum class StressType + { + S11, + S22, + S33 + }; + + RimFaultReactivationDataAccessorStress( RimFaultReactivation::Property property, double gradient, double seabedDepth ); + virtual ~RimFaultReactivationDataAccessorStress(); bool isMatching( RimFaultReactivation::Property property ) const override; @@ -64,32 +58,16 @@ public: double bottomDepth = std::numeric_limits::infinity(), size_t elementIndex = std::numeric_limits::max() ) const override; -private: - void updateResultAccessor() override; +protected: + virtual bool isDataAvailable() const = 0; - static RigFemResultAddress getResultAddress( const std::string& fieldName, const std::string& componentName ); + virtual double extractStressValue( StressType stressType, const cvf::Vec3d& position ) const = 0; - double interpolatedResultValue( RimWellIADataAccess& iaDataAccess, - const RigFemPart* femPart, - const cvf::Vec3d& position, - const std::vector& scalarResults ) const; + virtual std::pair calculatePorBar( const cvf::Vec3d& position, double gradient ) const = 0; - std::pair calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const; + virtual bool isPositionValid( const cvf::Vec3d& position, const cvf::Vec3d& topPosition, const cvf::Vec3d& bottomPosition ) const = 0; - RimGeoMechCase* m_geoMechCase; RimFaultReactivation::Property m_property; double m_gradient; double m_seabedDepth; - RigGeoMechCaseData* m_geoMechCaseData; - RigFemScalarResultFrames* m_s11Frames; - RigFemScalarResultFrames* m_s22Frames; - RigFemScalarResultFrames* m_s33Frames; - 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/RimFaultReactivationDataAccessorStressGeoMech.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp new file mode 100644 index 0000000000..dbd39426e1 --- /dev/null +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp @@ -0,0 +1,227 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2023 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 "RimFaultReactivationDataAccessorStressGeoMech.h" + +#include "RiaEclipseUnitTools.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 "RimFaultReactivationDataAccessorStress.h" +#include "RimFaultReactivationDataAccessorWellLogExtraction.h" +#include "RimFaultReactivationEnums.h" +#include "RimGeoMechCase.h" +#include "RimWellIADataAccess.h" + +#include "cvfVector3.h" + +#include +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RimFaultReactivationDataAccessorStressGeoMech::RimFaultReactivationDataAccessorStressGeoMech( RimGeoMechCase* geoMechCase, + RimFaultReactivation::Property property, + double gradient, + double seabedDepth ) + : RimFaultReactivationDataAccessorStress( property, gradient, seabedDepth ) + , m_geoMechCase( geoMechCase ) +{ + m_geoMechCaseData = geoMechCase->geoMechData(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RimFaultReactivationDataAccessorStressGeoMech::~RimFaultReactivationDataAccessorStressGeoMech() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor() +{ + const int partIndex = 0; + + auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr, int timeStepIndex ) -> RigFemScalarResultFrames* + { + 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 ); + 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 = + RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, + faultBottomPosition, + m_seabedDepth, + faultNormal * distanceFromFault ); + m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); + m_partIndexA = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); + m_extractorA = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceAWellPath.p(), errorName ); + } + + { + std::vector wellPoints = + RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, + faultBottomPosition, + m_seabedDepth, + -faultNormal * distanceFromFault ); + m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); + m_partIndexB = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); + m_extractorB = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceBWellPath.p(), errorName ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemResultAddress RimFaultReactivationDataAccessorStressGeoMech::getResultAddress( const std::string& fieldName, + const std::string& componentName ) +{ + return RigFemResultAddress( RIG_ELEMENT_NODAL, fieldName, componentName ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RimFaultReactivationDataAccessorStressGeoMech::isDataAvailable() const +{ + return m_s11Frames && m_s22Frames && m_s33Frames && m_femPart; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RimFaultReactivationDataAccessorStressGeoMech::extractStressValue( StressType stressType, const cvf::Vec3d& position ) const +{ + RimWellIADataAccess iaDataAccess( m_geoMechCase ); + + int timeStepIndex = 0; + + RigFemScalarResultFrames* frames = dataFrames( stressType ); + int frameIndex = frames->frameCount( timeStepIndex ) - 1; + const std::vector& s11Data = frames->frameData( timeStepIndex, frameIndex ); + return interpolatedResultValue( iaDataAccess, m_femPart, position, s11Data ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* RimFaultReactivationDataAccessorStressGeoMech::dataFrames( StressType stressType ) const +{ + if ( stressType == StressType::S11 ) + return m_s11Frames; + else if ( stressType == StressType::S22 ) + return m_s22Frames; + else + { + CAF_ASSERT( stressType == StressType::S33 ); + return m_s33Frames; + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position, double gradient ) const +{ + int timeStepIndex = 0; + int frameIndex = m_s33Frames->frameCount( timeStepIndex ) - 1; + return calculatePorBar( position, m_gradient, timeStepIndex, frameIndex ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RimFaultReactivationDataAccessorStressGeoMech::isPositionValid( const cvf::Vec3d& position, + const cvf::Vec3d& topPosition, + const cvf::Vec3d& bottomPosition ) const +{ + RimWellIADataAccess iaDataAccess( m_geoMechCase ); + int centerElementIdx = iaDataAccess.elementIndex( position ); + int bottomElementIdx = iaDataAccess.elementIndex( bottomPosition ); + int topElementIdx = iaDataAccess.elementIndex( topPosition ); + return ( centerElementIdx != -1 && topElementIdx != -1 && bottomElementIdx != -1 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RimFaultReactivationDataAccessorStressGeoMech::interpolatedResultValue( RimWellIADataAccess& iaDataAccess, + const RigFemPart* femPart, + const cvf::Vec3d& position, + const std::vector& scalarResults ) const +{ + return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position, + double gradient, + int timeStepIndex, + int frameIndex ) const +{ + RigFemPartCollection* partCollection = m_geoMechCaseData->femParts(); + cvf::ref extractor = m_partIndexA == partCollection->getPartIndexFromPoint( 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 ); + + return RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient ); +} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h new file mode 100644 index 0000000000..dd3a9a6dbb --- /dev/null +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h @@ -0,0 +1,93 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2023 - 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 "RimFaultReactivationDataAccessorStress.h" +#include "RimFaultReactivationEnums.h" + +#include "RigFemResultAddress.h" + +#include + +#include "cafPdmField.h" +#include "cafPdmPtrField.h" + +#include "cvfObject.h" + +class RigFemPart; +class RigFemPartResultsCollection; +class RigFemScalarResultFrames; +class RigGeoMechCaseData; +class RigGriddedPart3d; +class RimGeoMechCase; +class RimWellIADataAccess; +class RimModeledWellPath; +class RigWellPath; +class RigFemPartCollection; +class RigGeoMechWellLogExtractor; + +//================================================================================================== +/// +/// +//================================================================================================== +class RimFaultReactivationDataAccessorStressGeoMech : public RimFaultReactivationDataAccessorStress +{ +public: + RimFaultReactivationDataAccessorStressGeoMech( RimGeoMechCase* geoMechCase, + RimFaultReactivation::Property property, + double gradient, + double seabedDepth ); + ~RimFaultReactivationDataAccessorStressGeoMech() override; + +private: + void updateResultAccessor() override; + + bool isDataAvailable() const override; + + double extractStressValue( StressType stressType, const cvf::Vec3d& position ) const override; + + std::pair calculatePorBar( const cvf::Vec3d& position, double gradient ) const override; + + bool isPositionValid( const cvf::Vec3d& position, const cvf::Vec3d& topPosition, const cvf::Vec3d& bottomPosition ) const override; + + static RigFemResultAddress getResultAddress( const std::string& fieldName, const std::string& componentName ); + + double interpolatedResultValue( RimWellIADataAccess& iaDataAccess, + const RigFemPart* femPart, + const cvf::Vec3d& position, + const std::vector& scalarResults ) const; + + std::pair calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const; + + RigFemScalarResultFrames* dataFrames( StressType stressType ) const; + + RimGeoMechCase* m_geoMechCase; + RigGeoMechCaseData* m_geoMechCaseData; + RigFemScalarResultFrames* m_s11Frames; + RigFemScalarResultFrames* m_s22Frames; + RigFemScalarResultFrames* m_s33Frames; + 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; +};