#10826 Fault Reactivation: improve pore pressure eclipse extraction.

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})

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 "RigWellPath.h"
-#include "RimEclipseCase.h"
 
 #include <cmath>
+#include <limits>
 
 //--------------------------------------------------------------------------------------------------
 ///
@@ -62,16 +67,32 @@ RimFaultReactivationDataAccessorPorePressure::~RimFaultReactivationDataAccessorP
 //--------------------------------------------------------------------------------------------------
 void RimFaultReactivationDataAccessorPorePressure::updateResultAccessor()
 {
-    if ( m_caseData )
+    if ( !m_caseData ) return;
-    {
+
     RigEclipseResultAddress resVarAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "PRESSURE" );
     m_eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( resVarAddress );
 
-        m_resultAccessor = RigResultAccessorFactory::createFromResultAddress( m_caseData,
+    m_resultAccessor =
-                                                                              0,
+        RigResultAccessorFactory::createFromResultAddress( m_caseData, 0, RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress );
-                                                                              RiaDefines::PorosityModelType::MATRIX_MODEL,
+
-                                                                              m_timeStep,
+    auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom();
-                                                                              resVarAddress );
+    auto   faultNormal                           = m_model->faultNormal();
+    double distanceFromFault                     = 1.0;
+
+    for ( auto gridPart : m_model->allGridParts() )
+    {
+        double                  sign = m_model->normalPointsAt() == gridPart ? 1.0 : -1.0;
+        std::vector<cvf::Vec3d> wellPoints =
+            RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
+                                                                                   faultBottomPosition,
+                                                                                   sign * faultNormal * distanceFromFault );
+        cvf::ref<RigWellPath> wellPath =
+            new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
+        m_wellPaths[gridPart] = wellPath;
+
+        std::string                          errorName = "fault reactivation data access";
+        cvf::ref<RigEclipseWellLogExtractor> extractor = new RigEclipseWellLogExtractor( m_caseData, wellPath.p(), errorName );
+        m_extractors[gridPart]                         = extractor;
     }
 }
 
@@ -95,24 +116,30 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
 {
     if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
     {
-        auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
+        CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() );
-        if ( cellIdx != cvf::UNDEFINED_SIZE_T )
+        auto extractor = m_extractors.find( gridPart )->second;
-        {
+
-            double value = m_resultAccessor->cellScalar( cellIdx );
+        // Extract values along well path
-            if ( !std::isinf( value ) )
+        std::vector<double> values;
-            {
+        extractor->curveData( m_resultAccessor.p(), &values );
-                return 100000.0 * value; // return in pascal, not bar
+
-            }
+        auto intersections = extractor->intersections();
-        }
+
+        CAF_ASSERT( m_wellPaths.find( gridPart ) != m_wellPaths.end() );
+        auto wellPath = m_wellPaths.find( gridPart )->second;
+
+        // Insert top of overburden point
+        intersections.insert( intersections.begin(), wellPath->wellPathPoints().front() );
+        values.insert( values.begin(), std::numeric_limits<double>::infinity() );
+
+        // Insert bottom of underburden point
+        intersections.push_back( wellPath->wellPathPoints().back() );
+        values.push_back( std::numeric_limits<double>::infinity() );
+
+        auto [value, pos] =
+            RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( intersections, values, position, m_defaultPorePressureGradient );
+        return RiaEclipseUnitTools::barToPascal( value );
     }
 
-    return calculatePorePressure( std::abs( position.z() ), m_defaultPorePressureGradient );
+    return std::numeric_limits<double>::infinity();
-}
-
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-double RimFaultReactivationDataAccessorPorePressure::calculatePorePressure( double depth, double gradient )
-{
-    return gradient * 9.81 * depth * 1000.0;
 }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h

@@ -18,16 +18,17 @@
 
 #pragma once
 
+#include "RimFaultReactivationDataAccessor.h"
+
 #include "cvfObject.h"
 #include "cvfVector3.h"
 
-#include "RimFaultReactivationDataAccessor.h"
-#include "RimFaultReactivationEnums.h"
-
 class RimEclipseCase;
 class RigEclipseCaseData;
 class RigMainGrid;
 class RigResultAccessor;
+class RigEclipseWellLogExtractor;
+class RigWellPath;
 
 //==================================================================================================
 ///
@@ -58,4 +59,7 @@ private:
     const RigMainGrid*          m_mainGrid;
     double                      m_defaultPorePressureGradient;
     cvf::ref<RigResultAccessor> m_resultAccessor;
+
+    std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>                m_wellPaths;
+    std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>> m_extractors;
 };

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp

@@ -19,7 +19,6 @@
 #include "RimFaultReactivationDataAccessorStress.h"
 
 #include "RiaEclipseUnitTools.h"
-#include "RiaInterpolationTools.h"
 #include "RiaLogging.h"
 
 #include "RigFaultReactivationModel.h"
@@ -34,8 +33,8 @@
 #include "RigResultAccessorFactory.h"
 #include "RigWellPath.h"
 
+#include "RimFaultReactivationDataAccessorWellLogExtraction.h"
 #include "RimFaultReactivationEnums.h"
-#include "RimFracture.h"
 #include "RimGeoMechCase.h"
 #include "RimWellIADataAccess.h"
 
@@ -97,15 +96,21 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor()
     std::string           errorName             = "fault reactivation data access";
 
     {
-        std::vector<cvf::Vec3d> wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, faultNormal * distanceFromFault );
+        std::vector<cvf::Vec3d> wellPoints =
-        m_faceAWellPath                    = new RigWellPath( wellPoints, generateMds( wellPoints ) );
+            RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
+                                                                                   faultBottomPosition,
+                                                                                   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 = generateWellPoints( faultTopPosition, faultBottomPosition, -faultNormal * distanceFromFault );
+        std::vector<cvf::Vec3d> wellPoints =
-        m_faceBWellPath                    = new RigWellPath( wellPoints, generateMds( wellPoints ) );
+            RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
+                                                                                   faultBottomPosition,
+                                                                                   -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 );
     }
@@ -156,14 +161,14 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
 
         if ( m_property == RimFaultReactivation::Property::StressTop )
         {
-            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, topPosition, m_gradient, timeStepIndex, frameIndex );
+            auto [porBar, extractionPos] = calculatePorBar( 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 )
         {
-            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, bottomPosition, m_gradient, timeStepIndex, frameIndex );
+            auto [porBar, extractionPos] = calculatePorBar( bottomPosition, m_gradient, timeStepIndex, frameIndex );
             if ( std::isinf( porBar ) ) return porBar;
             double s33 = interpolatedResultValue( iaDataAccess, m_femPart, extractionPos, s33Data );
             return RiaEclipseUnitTools::barToPascal( s33 - porBar );
@@ -178,7 +183,7 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
         }
         else if ( m_property == RimFaultReactivation::Property::LateralStressComponentX )
         {
-            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, position, m_gradient, timeStepIndex, frameIndex );
+            auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, timeStepIndex, frameIndex );
             if ( std::isinf( porBar ) ) return porBar;
 
             const std::vector<float>& s11Data = m_s11Frames->frameData( timeStepIndex, frameIndex );
@@ -188,7 +193,7 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
         }
         else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY )
         {
-            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, position, m_gradient, timeStepIndex, frameIndex );
+            auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, timeStepIndex, frameIndex );
             if ( std::isinf( porBar ) ) return porBar;
 
             const std::vector<float>& s22Data = m_s22Frames->frameData( timeStepIndex, frameIndex );
@@ -215,12 +220,8 @@ double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellI
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStress::getPorBar( RimWellIADataAccess& iaDataAccess,
+std::pair<double, cvf::Vec3d>
-                                                                                 const RigFemPart*    femPart,
+    RimFaultReactivationDataAccessorStress::calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const
-                                                                                 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
@@ -235,180 +236,5 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStress::getPorBar(
     std::vector<double> values;
     extractor->curveData( resAddr, timeStepIndex, frameIndex, &values );
 
-    // Fill in missing values
+    return RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient );
-    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<double> xs     = { intersections[bottomIdx].z(), intersections[topIdx].z() };
-            std::vector<double> 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<double>::infinity(), cvf::Vec3d::UNDEFINED };
-}
-
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-std::pair<int, int> RimFaultReactivationDataAccessorStress::findIntersectionsForTvd( const std::vector<cvf::Vec3d>& 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<int>( i ) - 1;
-                bottomIdx = static_cast<int>( i );
-                break;
-            }
-        }
-    }
-    return std::make_pair( topIdx, bottomIdx );
-}
-
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-std::pair<int, int> RimFaultReactivationDataAccessorStress::findOverburdenAndUnderburdenIndex( const std::vector<double>& values )
-{
-    auto findLastOverburdenIndex = []( const std::vector<double>& values )
-    {
-        for ( size_t i = 0; i < values.size(); i++ )
-        {
-            if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
-        }
-
-        return -1;
-    };
-
-    auto findFirstUnderburdenIndex = []( const std::vector<double>& values )
-    {
-        for ( size_t i = values.size() - 1; i > 0; i-- )
-        {
-            if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
-        }
-
-        return -1;
-    };
-
-    int lastOverburdenIndex   = findLastOverburdenIndex( values );
-    int firstUnderburdenIndex = findFirstUnderburdenIndex( values );
-    return { lastOverburdenIndex, firstUnderburdenIndex };
-}
-
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-double RimFaultReactivationDataAccessorStress::computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
-                                                                          const std::vector<double>&     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<cvf::Vec3d>& intersections,
-                                                                  std::vector<double>&           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<int>( 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<double> intersectionsZ;
-    for ( auto i : intersections )
-    {
-        intersectionsZ.push_back( i.z() );
-    }
-
-    RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values );
-}
-
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-std::vector<double> RimFaultReactivationDataAccessorStress::generateMds( const std::vector<cvf::Vec3d>& points )
-{
-    CAF_ASSERT( points.size() >= 2 );
-
-    // Assume first at zero, all other points relative to that.
-    std::vector<double> 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<cvf::Vec3d> 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 };
 }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h

@@ -71,29 +71,7 @@ private:
                                     const cvf::Vec3d&         position,
                                     const std::vector<float>& scalarResults ) const;
 
-    std::pair<double, cvf::Vec3d> getPorBar( RimWellIADataAccess& iaDataAccess,
+    std::pair<double, cvf::Vec3d> calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const;
-                                             const RigFemPart*    femPart,
-                                             const cvf::Vec3d&    position,
-                                             double               gradient,
-                                             int                  timeStepIndex,
-                                             int                  frameIndex ) const;
-
-    static std::pair<bool, RimFaultReactivation::ElementSets>
-        findElementSetContainingElement( const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets,
-                                         unsigned int                                                                  elmIdx );
-
-    static std::pair<int, int> findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd );
-    static std::pair<int, int> findOverburdenAndUnderburdenIndex( const std::vector<double>& values );
-    static double              computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
-                                                          const std::vector<double>&     values,
-                                                          int                            i1,
-                                                          int                            i2,
-                                                          double                         gradient );
-    static void fillInMissingValues( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, double gradient );
-
-    static std::vector<double> generateMds( const std::vector<cvf::Vec3d>& points );
-    static std::vector<cvf::Vec3d>
-        generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset );
 
     RimGeoMechCase*                m_geoMechCase;
     RimFaultReactivation::Property m_property;

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 <http://www.gnu.org/licenses/gpl.html>
+//  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 <cmath>
+#include <limits>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimFaultReactivationDataAccessorWellLogExtraction::RimFaultReactivationDataAccessorWellLogExtraction()
+{
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimFaultReactivationDataAccessorWellLogExtraction::~RimFaultReactivationDataAccessorWellLogExtraction()
+{
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( const std::vector<cvf::Vec3d>& intersections,
+                                                                                                  std::vector<double>& 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<double> xs     = { intersections[bottomIdx].z(), intersections[topIdx].z() };
+            std::vector<double> 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<double>::infinity(), cvf::Vec3d::UNDEFINED };
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findIntersectionsForTvd( const std::vector<cvf::Vec3d>& 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<int>( i ) - 1;
+                bottomIdx = static_cast<int>( i );
+                break;
+            }
+        }
+    }
+    return std::make_pair( topIdx, bottomIdx );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findOverburdenAndUnderburdenIndex( const std::vector<double>& values )
+{
+    auto findLastOverburdenIndex = []( const std::vector<double>& values )
+    {
+        for ( size_t i = 0; i < values.size(); i++ )
+        {
+            if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
+        }
+
+        return -1;
+    };
+
+    auto findFirstUnderburdenIndex = []( const std::vector<double>& values )
+    {
+        for ( size_t i = values.size() - 1; i > 0; i-- )
+        {
+            if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
+        }
+
+        return -1;
+    };
+
+    int lastOverburdenIndex   = findLastOverburdenIndex( values );
+    int firstUnderburdenIndex = findFirstUnderburdenIndex( values );
+    return { lastOverburdenIndex, firstUnderburdenIndex };
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorWellLogExtraction::computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
+                                                                                     const std::vector<double>&     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<cvf::Vec3d>& intersections,
+                                                                             std::vector<double>&           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<int>( 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<double> intersectionsZ;
+    for ( auto i : intersections )
+    {
+        intersectionsZ.push_back( i.z() );
+    }
+
+    RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<cvf::Vec3d> 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<double> RimFaultReactivationDataAccessorWellLogExtraction::generateMds( const std::vector<cvf::Vec3d>& points )
+{
+    CAF_ASSERT( points.size() >= 2 );
+
+    // Assume first at zero, all other points relative to that.
+    std::vector<double> 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;
+}

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h

@@ -0,0 +1,52 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  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 <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#pragma once
+
+#include "RimFaultReactivationDataAccessor.h"
+#include "RimFaultReactivationEnums.h"
+
+#include <vector>
+
+//==================================================================================================
+///
+///
+//==================================================================================================
+class RimFaultReactivationDataAccessorWellLogExtraction
+{
+public:
+    RimFaultReactivationDataAccessorWellLogExtraction();
+    ~RimFaultReactivationDataAccessorWellLogExtraction();
+
+    static std::pair<double, cvf::Vec3d>
+        calculatePorBar( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, const cvf::Vec3d& position, double gradient );
+    static std::vector<cvf::Vec3d>
+        generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset );
+
+    static std::vector<double> generateMds( const std::vector<cvf::Vec3d>& points );
+
+protected:
+    static std::pair<int, int> findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd );
+    static std::pair<int, int> findOverburdenAndUnderburdenIndex( const std::vector<double>& values );
+    static double              computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
+                                                          const std::vector<double>&     values,
+                                                          int                            i1,
+                                                          int                            i2,
+                                                          double                         gradient );
+    static void fillInMissingValues( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, double gradient );
+};

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp

@@ -20,11 +20,6 @@
 
 #include "RigFaultReactivationModelGenerator.h"
 #include "RigGriddedPart3d.h"
-#include "RigPolyLinesData.h"
-
-#include "RimFaultReactivationDataAccess.h"
-
-#include "cafAssert.h"
 
 //--------------------------------------------------------------------------------------------------
 ///
@@ -147,6 +142,7 @@ void RigFaultReactivationModel::updateGeometry( size_t startCell, cvf::StructGri
 
     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 +150,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 +225,11 @@ const std::pair<cvf::Vec3d, cvf::Vec3d> 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;
+}

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.h

@@ -85,6 +85,8 @@ public:
     const cvf::Vec3d                        faultNormal() const;
     const std::pair<cvf::Vec3d, cvf::Vec3d> faultTopBottom() const;
 
+    RimFaultReactivation::GridPart normalPointsAt() const;
+
 private:
     std::shared_ptr<RigFaultReactivationModelGenerator> m_generator;
 
@@ -94,4 +96,5 @@ private:
     bool m_isValid;
 
     std::map<GridPart, RigGriddedPart3d*> m_3dparts;
+    RimFaultReactivation::GridPart        m_normalPointsAt;
 };