From 761a12aa3aa821357c90937f989d5d6ed031a147 Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Wed, 17 Jan 2024 14:08:08 +0100 Subject: [PATCH] Fix extraction egde cases. --- ...tReactivationDataAccessorStressEclipse.cpp | 84 +++++++++++-------- ...ultReactivationDataAccessorStressEclipse.h | 5 ++ ...tReactivationDataAccessorStressGeoMech.cpp | 8 +- ...ctivationDataAccessorWellLogExtraction.cpp | 16 ++-- 4 files changed, 68 insertions(+), 45 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp index cb46d62a3d..ea11ea6836 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp @@ -162,6 +162,10 @@ double RimFaultReactivationDataAccessorStressEclipse::extractStressValue( Stress return RiaEclipseUnitTools::pascalToBar( RiaInterpolationTools::linear( xs, ys, position.z() ) ); } } + else if ( position.z() <= intersections.back().z() ) + { + return RiaEclipseUnitTools::pascalToBar( stressValues.back() ); + } return std::numeric_limits::infinity(); } @@ -219,35 +223,11 @@ std::vector double seaWaterLoad = gravity * std::abs( seabedDepth ) * waterDensity; std::vector values = { seaWaterLoad }; - auto g = model.grid( gridPart ); - auto elementSets = g->elementSets(); - - auto getElementSet = - []( auto part, const cvf::Vec3d& point, const std::map>& elementSets ) - -> std::pair - { - for ( auto [elementSet, elements] : elementSets ) - { - for ( unsigned int elementIndex : elements ) - { - // Find nodes from element - auto positions = part->elementCorners( elementIndex ); - - std::array coordinates; - for ( size_t i = 0; i < positions.size(); ++i ) - { - coordinates[i] = positions[i]; - } - - if ( RigHexIntersectionTools::isPointInCell( point, coordinates.data() ) ) - { - return { true, elementSet }; - } - } - } + auto part = model.grid( gridPart ); + CAF_ASSERT( part ); + auto elementSets = part->elementSets(); - return { false, RimFaultReactivation::ElementSets::OverBurden }; - }; + double previousDensity = densities.find( RimFaultReactivation::ElementSets::OverBurden )->second; for ( size_t i = 1; i < intersections.size(); i++ ) { @@ -256,26 +236,56 @@ std::vector double currentDepth = intersections[i].z(); double deltaDepth = previousDepth - currentDepth; + double density = previousDensity; - auto [isOk, elementSet] = getElementSet( g, intersections[i], elementSets ); + auto [isOk, elementSet] = findElementSetForPoint( *part, intersections[i], elementSets ); if ( isOk ) { // Unit: kg/m^3 CAF_ASSERT( densities.find( elementSet ) != densities.end() ); - double density = densities.find( elementSet )->second; - - double value = previousValue + density * gravity * deltaDepth; - values.push_back( value ); - } - else - { - values.push_back( std::numeric_limits::infinity() ); + density = densities.find( elementSet )->second; } + + double value = previousValue + density * gravity * deltaDepth; + values.push_back( value ); + + previousDensity = density; } return values; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStressEclipse::findElementSetForPoint( + const RigGriddedPart3d& part, + const cvf::Vec3d& point, + const std::map>& elementSets ) +{ + for ( auto [elementSet, elements] : elementSets ) + { + for ( unsigned int elementIndex : elements ) + { + // Find nodes from element + auto positions = part.elementCorners( elementIndex ); + + std::array coordinates; + for ( size_t i = 0; i < positions.size(); ++i ) + { + coordinates[i] = positions[i]; + } + + if ( RigHexIntersectionTools::isPointInCell( point, coordinates.data() ) ) + { + return { true, elementSet }; + } + } + } + + return { false, RimFaultReactivation::ElementSets::OverBurden }; +}; + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h index 31397f951f..86050d2f02 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h @@ -80,6 +80,11 @@ class RimFaultReactivationDataAccessorStressEclipse : public RimFaultReactivatio static void addOverburdenAndUnderburdenPoints( std::vector& intersections, const std::vector& wellPathPoints ); + static std::pair + findElementSetForPoint( const RigGriddedPart3d& part, + const cvf::Vec3d& point, + const std::map>& elementSets ); + RimEclipseCase* m_eclipseCase; RigEclipseCaseData* m_caseData; const RigMainGrid* m_mainGrid; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp index e36d6c0775..8ce971af81 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp @@ -90,9 +90,11 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor() m_s22Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ), timeStepIndex ); auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom(); - auto faultNormal = m_model->faultNormal(); - double distanceFromFault = 1.0; - auto [topDepth, bottomDepth] = m_model->depthTopBottom(); + auto faultNormal = m_model->faultNormal() ^ cvf::Vec3d::Z_AXIS; + faultNormal.normalize(); + + double distanceFromFault = 1.0; + auto [topDepth, bottomDepth] = m_model->depthTopBottom(); RigFemPartCollection* geoMechPartCollection = m_geoMechCaseData->femParts(); std::string errorName = "fault reactivation data access"; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp index adf6789209..d6e4b6bb8c 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp @@ -115,6 +115,10 @@ std::pair return { porBar, extractionPosition }; } } + else if ( position.z() <= intersections.back().z() ) + { + return { values.back(), intersections.back() }; + } return { std::numeric_limits::infinity(), cvf::Vec3d::UNDEFINED }; } @@ -281,16 +285,18 @@ std::pair>, std:: double seabedDepth ) { auto [faultTopPosition, faultBottomPosition] = model.faultTopBottom(); - auto faultNormal = model.faultNormal(); - double distanceFromFault = 1.0; - auto [topDepth, bottomDepth] = model.depthTopBottom(); + auto faultNormal = model.faultNormal() ^ cvf::Vec3d::Z_AXIS; + faultNormal.normalize(); + + double distanceFromFault = 1.0; + auto [topDepth, bottomDepth] = model.depthTopBottom(); std::map> wellPaths; std::map> extractors; for ( auto gridPart : model.allGridParts() ) { - double sign = model.normalPointsAt() == gridPart ? 1.0 : -1.0; + double sign = model.normalPointsAt() == gridPart ? -1.0 : 1.0; std::vector wellPoints = RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, faultBottomPosition, @@ -350,7 +356,7 @@ std::vector RimFaultReactivationDataAccessorWellLogExtraction::extractDe std::vector intersectionsZ; for ( auto i : intersections ) { - intersectionsZ.push_back( i.z() ); + intersectionsZ.push_back( -i.z() ); } return intersectionsZ; }