diff --git a/ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake b/ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake index c904cfb827..bf1a368e00 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}/RimFaultReactivationDataAccessorWellLogExtraction.h ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationEnums.h ) @@ -27,6 +28,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}/RimFaultReactivationDataAccessorWellLogExtraction.cpp ) list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES}) diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp index e9a386b1e9..7015dece50 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp @@ -17,21 +17,26 @@ ///////////////////////////////////////////////////////////////////////////////// #include "RimFaultReactivationDataAccessorPorePressure.h" +#include "RimEclipseCase.h" #include "RimFaciesProperties.h" +#include "RimFaultReactivationDataAccessorWellLogExtraction.h" #include "RimFaultReactivationEnums.h" #include "RiaDefines.h" +#include "RiaEclipseUnitTools.h" #include "RiaPorosityModel.h" #include "RigCaseCellResultsData.h" #include "RigEclipseCaseData.h" #include "RigEclipseResultAddress.h" +#include "RigEclipseWellLogExtractor.h" +#include "RigFaultReactivationModel.h" #include "RigMainGrid.h" #include "RigResultAccessorFactory.h" - -#include "RimEclipseCase.h" +#include "RigWellPath.h" #include +#include //-------------------------------------------------------------------------------------------------- /// @@ -72,6 +77,24 @@ void RimFaultReactivationDataAccessorPorePressure::updateResultAccessor() RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress ); + + auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom(); + auto faultNormal = m_model->faultNormal(); + double distanceFromFault = 1.0; + + std::string errorName = "fault reactivation data access"; + + { + std::vector wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, faultNormal * distanceFromFault ); + m_faceWellPathA = new RigWellPath( wellPoints, generateMds( wellPoints ) ); + m_extractorA = new RigEclipseWellLogExtractor( m_caseData, m_faceWellPathA.p(), errorName ); + } + + { + std::vector wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, -faultNormal * distanceFromFault ); + m_faceWellPathB = new RigWellPath( wellPoints, generateMds( wellPoints ) ); + m_extractorB = new RigEclipseWellLogExtractor( m_caseData, m_faceWellPathB.p(), errorName ); + } } } @@ -95,24 +118,23 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf: { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { - auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position ); - if ( cellIdx != cvf::UNDEFINED_SIZE_T ) - { - double value = m_resultAccessor->cellScalar( cellIdx ); - if ( !std::isinf( value ) ) - { - return 100000.0 * value; // return in pascal, not bar - } - } - } + // TODO: choose correct extractor + auto extractor = m_extractorA; + std::vector values; + extractor->curveData( m_resultAccessor.p(), &values ); - return calculatePorePressure( std::abs( position.z() ), m_defaultPorePressureGradient ); -} + auto intersections = extractor->intersections(); -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimFaultReactivationDataAccessorPorePressure::calculatePorePressure( double depth, double gradient ) -{ - return gradient * 9.81 * depth * 1000.0; + intersections.insert( intersections.begin(), m_faceWellPathA->wellPathPoints().front() ); + values.insert( values.begin(), std::numeric_limits::infinity() ); + + intersections.push_back( m_faceWellPathA->wellPathPoints().back() ); + values.push_back( std::numeric_limits::infinity() ); + + auto [value, pos] = + RimFaultReactivationDataAccessorWellLogExtraction::getPorBar( intersections, values, position, m_defaultPorePressureGradient ); + return RiaEclipseUnitTools::barToPascal( value ); + } + + return std::numeric_limits::infinity(); } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h index 647537a7cb..e6a6357cc3 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h @@ -28,6 +28,8 @@ class RimEclipseCase; class RigEclipseCaseData; class RigMainGrid; class RigResultAccessor; +class RigEclipseWellLogExtractor; +class RigWellPath; //================================================================================================== /// @@ -58,4 +60,9 @@ class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivation const RigMainGrid* m_mainGrid; double m_defaultPorePressureGradient; cvf::ref m_resultAccessor; + + cvf::ref m_faceWellPathA; + cvf::ref m_faceWellPathB; + cvf::ref m_extractorA; + cvf::ref m_extractorB; }; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp index 15f11388ce..12813b0674 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp @@ -34,6 +34,7 @@ #include "RigResultAccessorFactory.h" #include "RigWellPath.h" +#include "RimFaultReactivationDataAccessorWellLogExtraction.h" #include "RimFaultReactivationEnums.h" #include "RimFracture.h" #include "RimGeoMechCase.h" @@ -235,180 +236,5 @@ std::pair RimFaultReactivationDataAccessorStress::getPorBar( 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 calculatePorePressure = []( double depth, double gradient ) - { return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 ); }; - - auto computeGradient = []( double depth1, double value1, double depth2, double value2 ) - { return ( value2 - value1 ) / ( depth2 - depth1 ); }; - - auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); - - // Fill in overburden values using gradient - double topPorePressure = calculatePorePressure( std::abs( intersections[0].z() ), gradient ); - double overburdenGradient = - computeGradient( intersections[0].z(), topPorePressure, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] ); - - for ( int i = 0; i < lastOverburdenIndex; i++ ) - { - values[i] = computePorBarWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient ); - } - - // Fill in underburden values using gradient - int lastElementIndex = static_cast( values.size() ) - 1; - double bottomPorePressure = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient ); - double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(), - values[firstUnderburdenIndex], - intersections[lastElementIndex].z(), - bottomPorePressure ); - - for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- ) - { - values[i] = computePorBarWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient ); - } - - // Interpolate the missing values (should only be intra-reservoir by now) - 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 }; + return RimFaultReactivationDataAccessorWellLogExtraction::getPorBar( extractor->intersections(), values, position, gradient ); } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h index b50f128df2..cc31f6a655 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h @@ -78,23 +78,6 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc int timeStepIndex, int frameIndex ) const; - static std::pair - findElementSetContainingElement( const std::map>& elementSets, - unsigned int elmIdx ); - - 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; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp new file mode 100644 index 0000000000..07d1598c33 --- /dev/null +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp @@ -0,0 +1,245 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "RimFaultReactivationDataAccessorWellLogExtraction.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 + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RimFaultReactivationDataAccessorWellLogExtraction::RimFaultReactivationDataAccessorWellLogExtraction() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RimFaultReactivationDataAccessorWellLogExtraction::~RimFaultReactivationDataAccessorWellLogExtraction() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorWellLogExtraction::getPorBar( const std::vector& intersections, + std::vector& values, + const cvf::Vec3d& position, + double gradient ) +{ + // Fill in missing values + 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 RimFaultReactivationDataAccessorWellLogExtraction::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 RimFaultReactivationDataAccessorWellLogExtraction::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 RimFaultReactivationDataAccessorWellLogExtraction::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 RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValues( const std::vector& intersections, + std::vector& values, + double gradient ) +{ + CAF_ASSERT( intersections.size() == values.size() ); + + auto calculatePorePressure = []( double depth, double gradient ) + { return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 ); }; + + auto computeGradient = []( double depth1, double value1, double depth2, double value2 ) + { return ( value2 - value1 ) / ( depth2 - depth1 ); }; + + auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); + + // Fill in overburden values using gradient + double topPorePressure = calculatePorePressure( std::abs( intersections[0].z() ), gradient ); + double overburdenGradient = + computeGradient( intersections[0].z(), topPorePressure, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] ); + + for ( int i = 0; i < lastOverburdenIndex; i++ ) + { + values[i] = computePorBarWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient ); + } + + // Fill in underburden values using gradient + int lastElementIndex = static_cast( values.size() ) - 1; + double bottomPorePressure = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient ); + double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(), + values[firstUnderburdenIndex], + intersections[lastElementIndex].z(), + bottomPorePressure ); + + for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- ) + { + values[i] = computePorBarWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient ); + } + + // Interpolate the missing values (should only be intra-reservoir by now) + std::vector intersectionsZ; + for ( auto i : intersections ) + { + intersectionsZ.push_back( i.z() ); + } + + RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimFaultReactivationDataAccessorWellLogExtraction::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 }; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimFaultReactivationDataAccessorWellLogExtraction::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; +} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h new file mode 100644 index 0000000000..e856d9be17 --- /dev/null +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h @@ -0,0 +1,72 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "RimFaultReactivationDataAccessor.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 RimFaultReactivationDataAccessorWellLogExtraction +{ +public: + RimFaultReactivationDataAccessorWellLogExtraction(); + ~RimFaultReactivationDataAccessorWellLogExtraction(); + + static std::pair + getPorBar( const std::vector& intersections, std::vector& values, const cvf::Vec3d& position, double gradient ); + +protected: + 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 + generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset ); + + static std::vector generateMds( const std::vector& points ); +}; diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp index e922fa0ece..4ae750ad16 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp @@ -24,6 +24,7 @@ #include "RimFaultReactivationDataAccess.h" +#include "RimFaultReactivationEnums.h" #include "cafAssert.h" //-------------------------------------------------------------------------------------------------- @@ -145,8 +146,9 @@ void RigFaultReactivationModel::updateGeometry( size_t startCell, cvf::StructGri { reset(); - auto frontPart = m_3dparts[GridPart::FW]; - auto backPart = m_3dparts[GridPart::HW]; + auto frontPart = m_3dparts[GridPart::FW]; + auto backPart = m_3dparts[GridPart::HW]; + m_normalPointsAt = GridPart::FW; m_generator->generateGeometry( startCell, startFace, frontPart, backPart ); @@ -154,6 +156,7 @@ void RigFaultReactivationModel::updateGeometry( size_t startCell, cvf::StructGri { m_3dparts[GridPart::HW] = frontPart; m_3dparts[GridPart::FW] = backPart; + m_normalPointsAt = GridPart::HW; } auto& frontPoints = m_generator->frontPoints(); @@ -228,3 +231,11 @@ const std::pair RigFaultReactivationModel::faultTopBotto if ( m_generator.get() == nullptr ) return std::make_pair( cvf::Vec3d(), cvf::Vec3d() ); return m_generator->faultTopBottomPoints(); } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RimFaultReactivation::GridPart RigFaultReactivationModel::normalPointsAt() const +{ + return m_normalPointsAt; +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.h b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.h index afd7077a5c..029c15c1c9 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.h +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.h @@ -85,6 +85,8 @@ class RigFaultReactivationModel : public cvf::Object const cvf::Vec3d faultNormal() const; const std::pair faultTopBottom() const; + RimFaultReactivation::GridPart normalPointsAt() const; + private: std::shared_ptr m_generator; @@ -94,4 +96,5 @@ class RigFaultReactivationModel : public cvf::Object bool m_isValid; std::map m_3dparts; + RimFaultReactivation::GridPart m_normalPointsAt; };