OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Refactor: create stress accessor base class.

Browse Source
This commit is contained in:
Kristian Bendiksen
2024-01-10 10:16:59 +01:00
parent 14faf607f4
commit cb74bf2f1e
6 changed files with 353 additions and 167 deletions
Download Patch File Download Diff File Expand all files Collapse all files

144
ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp
View File

@@ -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<cvf::Vec3d> 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<cvf::Vec3d> 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<double>::infinity();
RimWellIADataAccess iaDataAccess( m_geoMechCase );
int centerElementIdx = iaDataAccess.elementIndex( position );
if ( !isDataAvailable() ) return std::numeric_limits<double>::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<float>& 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<float>& 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<float>& 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<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellIADataAccess& iaDataAccess,
const RigFemPart* femPart,
const cvf::Vec3d& position,
const std::vector<float>& scalarResults ) const
{
return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<double, cvf::Vec3d>
RimFaultReactivationDataAccessorStress::calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const
{
RigFemPartCollection* partCollection = m_geoMechCaseData->femParts();
cvf::ref<RigGeoMechWellLogExtractor> extractor = m_partIndexA == partCollection->getPartIndexFromPoint( position ) ? m_extractorA
: m_extractorB;
if ( !extractor->valid() )
{
RiaLogging::error( "Invalid extractor when extracting PorBar" );
return { std::numeric_limits<double>::infinity(), cvf::Vec3d::UNDEFINED };
}
RigFemResultAddress resAddr = RigFemAddressDefines::nodalPorBarAddress();
std::vector<double> values;
extractor->curveData( resAddr, timeStepIndex, frameIndex, &values );
return RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient );
}