Fault Reactivation: extract stress data using well log extraction.

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp

@@ -71,7 +71,7 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( RimEclipseCase*
                                                                          RimFaultReactivation::Property::LateralStressComponentY };
         for ( auto property : stressProperties )
         {
-            m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorStress>( geoMechCase, property ) );
+            m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorStress>( geoMechCase, property, 0.003 ) );
         }
     }
 }
@@ -126,12 +126,18 @@ std::vector<double> RimFaultReactivationDataAccess::extractModelData( const RigF
     };
 
     std::shared_ptr<RimFaultReactivationDataAccessor> accessor = getAccessor( property );
+
     if ( accessor )
     {
-        accessor->setTimeStep( timeStep );
+        accessor->setModelAndTimeStep( model, timeStep );
 
         auto grid = model.grid( gridPart );
 
+        const std::map<RimFaultReactivation::BorderSurface, std::vector<unsigned int>>& borderSurfaceElements = grid->borderSurfaceElements();
+        auto it = borderSurfaceElements.find( RimFaultReactivation::BorderSurface::Seabed );
+        CAF_ASSERT( it != borderSurfaceElements.end() && "Sea bed border surface does not exist" );
+        std::set<unsigned int> seabedElements( it->second.begin(), it->second.end() );
+
         std::vector<double> values;
 
         if ( nodeProperties.contains( property ) )
@@ -149,11 +155,15 @@ std::vector<double> RimFaultReactivationDataAccess::extractModelData( const RigF
             {
                 std::vector<cvf::Vec3d> corners = grid->elementCorners( elementIndex );
 
-                double topDepth    = computeAverageDepth( corners, { 0, 1, 2, 3 } );
+                // Move top of sea bed element down to end up inside top element
+                bool   isTopElement   = seabedElements.contains( static_cast<unsigned int>( elementIndex ) );
+                double topDepthAdjust = isTopElement ? 0.1 : 0.0;
+
+                double topDepth    = computeAverageDepth( corners, { 0, 1, 2, 3 } ) - topDepthAdjust;
                 double bottomDepth = computeAverageDepth( corners, { 4, 5, 6, 7 } );
 
                 cvf::Vec3d position = RigCaseToCaseCellMapperTools::calculateCellCenter( corners.data() );
-                double     value    = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth );
+                double     value    = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth, elementIndex );
                 values.push_back( value );
             }
         }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.cpp

@@ -18,6 +18,7 @@
 
 #include "RimFaultReactivationDataAccessor.h"
 
+#include "RigFaultReactivationModel.h"
 #include "RigMainGrid.h"
 
 #include "cafAssert.h"
@@ -28,6 +29,7 @@
 RimFaultReactivationDataAccessor::RimFaultReactivationDataAccessor()
 {
     m_timeStep = -1;
+    m_model    = nullptr;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -40,8 +42,9 @@ RimFaultReactivationDataAccessor::~RimFaultReactivationDataAccessor()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RimFaultReactivationDataAccessor::setTimeStep( size_t timeStep )
+void RimFaultReactivationDataAccessor::setModelAndTimeStep( const RigFaultReactivationModel& model, size_t timeStep )
 {
+    m_model    = &model;
     m_timeStep = timeStep;
     updateResultAccessor();
 }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h

@@ -38,7 +38,7 @@ public:
     RimFaultReactivationDataAccessor();
     ~RimFaultReactivationDataAccessor();
 
-    virtual void setTimeStep( size_t timeStep );
+    virtual void setModelAndTimeStep( const RigFaultReactivationModel& model, size_t timeStep );
 
     virtual bool isMatching( RimFaultReactivation::Property property ) const = 0;
 
@@ -46,10 +46,12 @@ public:
                                     const RigFaultReactivationModel& model,
                                     RimFaultReactivation::GridPart   gridPart,
                                     double                           topDepth     = std::numeric_limits<double>::infinity(),
-                                    double                           bottomDepth = std::numeric_limits<double>::infinity() ) const = 0;
+                                    double                           bottomDepth  = std::numeric_limits<double>::infinity(),
+                                    size_t                           elementIndex = std::numeric_limits<size_t>::max() ) const = 0;
 
 protected:
     virtual void updateResultAccessor() = 0;
 
+    const RigFaultReactivationModel* m_model;
     size_t                           m_timeStep;
 };

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp

@@ -101,8 +101,8 @@ double RimFaultReactivationDataAccessorGeoMech::valueAtPosition( const cvf::Vec3
                                                                  const RigFaultReactivationModel& model,
                                                                  RimFaultReactivation::GridPart   gridPart,
                                                                  double                           topDepth,
-                                                                 double                           bottomDepth ) const
+                                                                 double                           bottomDepth,
-
+                                                                 size_t                           elementIndex ) const
 {
     if ( !m_resultFrames ) return std::numeric_limits<double>::infinity();
 

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h

@@ -44,7 +44,8 @@ public:
                             const RigFaultReactivationModel& model,
                             RimFaultReactivation::GridPart   gridPart,
                             double                           topDepth     = std::numeric_limits<double>::infinity(),
-                            double                           bottomDepth = std::numeric_limits<double>::infinity() ) const override;
+                            double                           bottomDepth  = std::numeric_limits<double>::infinity(),
+                            size_t                           elementIndex = std::numeric_limits<size_t>::max() ) const override;
 
 private:
     void updateResultAccessor() override;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp

@@ -90,7 +90,8 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
                                                                       const RigFaultReactivationModel& model,
                                                                       RimFaultReactivation::GridPart   gridPart,
                                                                       double                           topDepth,
-                                                                      double                           bottomDepth ) const
+                                                                      double                           bottomDepth,
+                                                                      size_t                           elementIndex ) const
 {
     if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
     {

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h

@@ -45,7 +45,8 @@ public:
                             const RigFaultReactivationModel& model,
                             RimFaultReactivation::GridPart   gridPart,
                             double                           topDepth     = std::numeric_limits<double>::infinity(),
-                            double                           bottomDepth = std::numeric_limits<double>::infinity() ) const override;
+                            double                           bottomDepth  = std::numeric_limits<double>::infinity(),
+                            size_t                           elementIndex = std::numeric_limits<size_t>::max() ) const override;
 
 private:
     void updateResultAccessor() override;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp

@@ -19,19 +19,28 @@
 #include "RimFaultReactivationDataAccessorStress.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>
 
@@ -39,9 +48,11 @@
 ///
 //--------------------------------------------------------------------------------------------------
 RimFaultReactivationDataAccessorStress::RimFaultReactivationDataAccessorStress( RimGeoMechCase*                geoMechCase,
-                                                                                RimFaultReactivation::Property property )
+                                                                                RimFaultReactivation::Property property,
+                                                                                double                         gradient )
     : m_geoMechCase( geoMechCase )
     , m_property( property )
+    , m_gradient( gradient )
 {
     m_geoMechCaseData = geoMechCase->geoMechData();
 }
@@ -60,10 +71,11 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor()
 {
     const int partIndex = 0;
 
-    auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr ) -> RigFemScalarResultFrames*
+    auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr, int timeStepIndex ) -> RigFemScalarResultFrames*
     {
         auto result     = femParts->findOrLoadScalarResult( partIndex, addr );
-        if ( result->frameData( 0, 0 ).empty() )
+        int  frameIndex = result->frameCount( timeStepIndex ) - 1;
+        if ( result->frameData( timeStepIndex, frameIndex ).empty() )
         {
             return nullptr;
         }
@@ -72,10 +84,63 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor()
 
     auto femParts     = m_geoMechCaseData->femPartResults();
     m_femPart         = femParts->parts()->part( partIndex );
-    m_s33Frames   = loadFrameLambda( femParts, getResultAddress( "ST", "S33" ) );
+    int timeStepIndex = 0;
-    m_s11Frames   = loadFrameLambda( femParts, getResultAddress( "ST", "S11" ) );
+    m_s33Frames       = loadFrameLambda( femParts, getResultAddress( "ST", "S33" ), timeStepIndex );
-    m_s22Frames   = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ) );
+    m_s11Frames       = loadFrameLambda( femParts, getResultAddress( "ST", "S11" ), timeStepIndex );
-    m_porFrames   = loadFrameLambda( femParts, RigFemAddressDefines::elementNodalPorBarAddress() );
+    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 = generateWellPoints( faultTopPosition, faultBottomPosition, faultNormal * distanceFromFault );
+        m_faceAWellPath                    = new RigWellPath( wellPoints, generateMds( wellPoints ) );
+        m_partIndexA                       = getPartIndexFromPoint( *geoMechPartCollection, wellPoints[1] );
+        m_extractorA                       = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceAWellPath.p(), errorName );
+    }
+
+    {
+        std::vector<cvf::Vec3d> wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, -faultNormal * distanceFromFault );
+        m_faceBWellPath                    = new RigWellPath( wellPoints, generateMds( wellPoints ) );
+        m_partIndexB                       = getPartIndexFromPoint( *geoMechPartCollection, wellPoints[1] );
+        m_extractorB                       = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceBWellPath.p(), errorName );
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+int RimFaultReactivationDataAccessorStress::getPartIndexFromPoint( const RigFemPartCollection& partCollection, const cvf::Vec3d& point )
+{
+    const int idx = 0;
+
+    // Find candidates for intersected global elements
+    const cvf::BoundingBox intersectingBb( point, point );
+    std::vector<size_t>    intersectedGlobalElementIndexCandidates;
+    partCollection.findIntersectingGlobalElementIndices( intersectingBb, &intersectedGlobalElementIndexCandidates );
+
+    if ( intersectedGlobalElementIndexCandidates.empty() ) return idx;
+
+    // Iterate through global element candidates and check if point is in hexCorners
+    for ( const auto& globalElementIndex : intersectedGlobalElementIndexCandidates )
+    {
+        const auto [part, elementIndex] = partCollection.partAndElementIndex( globalElementIndex );
+
+        // Find nodes from element
+        std::array<cvf::Vec3d, 8> coordinates;
+        const bool                isSuccess = part->fillElementCoordinates( elementIndex, coordinates );
+        if ( !isSuccess ) continue;
+
+        const bool isPointInCell = RigHexIntersectionTools::isPointInCell( point, coordinates.data() );
+        if ( isPointInCell ) return part->elementPartId();
+    }
+
+    // Utilize first part to have an id
+    return idx;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -101,9 +166,10 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
                                                                 const RigFaultReactivationModel& model,
                                                                 RimFaultReactivation::GridPart   gridPart,
                                                                 double                           topDepth,
-                                                                double                           bottomDepth ) const
+                                                                double                           bottomDepth,
+                                                                size_t                           elementIndex ) const
 {
-    if ( !m_porFrames || !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits<double>::infinity();
+    if ( !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits<double>::infinity();
 
     RimWellIADataAccess iaDataAccess( m_geoMechCase );
     int                 centerElementIdx = iaDataAccess.elementIndex( position );
@@ -116,23 +182,22 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
     if ( centerElementIdx != -1 && topElementIdx != -1 && bottomElementIdx != -1 )
     {
         int timeStepIndex = 0;
-        int frameIndex    = 0;
+        int frameIndex    = m_s33Frames->frameCount( timeStepIndex ) - 1;
 
-        const std::vector<float>& s11Data = m_s11Frames->frameData( timeStepIndex, frameIndex );
-        const std::vector<float>& s22Data = m_s22Frames->frameData( timeStepIndex, frameIndex );
         const std::vector<float>& s33Data = m_s33Frames->frameData( timeStepIndex, frameIndex );
-        const std::vector<float>& porData = m_porFrames->frameData( timeStepIndex, frameIndex );
 
         if ( m_property == RimFaultReactivation::Property::StressTop )
         {
-            double s33    = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, s33Data );
+            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, topPosition, m_gradient, timeStepIndex, frameIndex );
-            double porBar = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, porData );
+            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 )
         {
-            double s33    = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, s33Data );
+            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, bottomPosition, m_gradient, timeStepIndex, frameIndex );
-            double porBar = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, porData );
+            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::DepthTop )
@@ -145,16 +210,22 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
         }
         else if ( m_property == RimFaultReactivation::Property::LateralStressComponentX )
         {
-            double s11    = interpolatedResultValue( iaDataAccess, m_femPart, position, s11Data );
+            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, position, m_gradient, timeStepIndex, frameIndex );
-            double s33    = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data );
+            if ( std::isinf( porBar ) ) return porBar;
-            double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData );
+
+            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 );
             return ( s11 - porBar ) / ( s33 - porBar );
         }
         else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY )
         {
-            double s22    = interpolatedResultValue( iaDataAccess, m_femPart, position, s22Data );
+            auto [porBar, extractionPos] = getPorBar( iaDataAccess, m_femPart, position, m_gradient, timeStepIndex, frameIndex );
-            double s33    = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data );
+            if ( std::isinf( porBar ) ) return porBar;
-            double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData );
+
+            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 );
             return ( s22 - porBar ) / ( s33 - porBar );
         }
     }
@@ -172,3 +243,187 @@ double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellI
 {
     return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position );
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStress::getPorBar( RimWellIADataAccess& iaDataAccess,
+                                                                                 const RigFemPart*    femPart,
+                                                                                 const cvf::Vec3d&    position,
+                                                                                 double               gradient,
+                                                                                 int                  timeStepIndex,
+                                                                                 int                  frameIndex ) const
+{
+    RigFemPartCollection*                partCollection = m_geoMechCaseData->femParts();
+    cvf::ref<RigGeoMechWellLogExtractor> extractor      = m_partIndexA == getPartIndexFromPoint( *partCollection, 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 );
+
+    // 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<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 [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values );
+
+    // Fill in overburden values using gradient
+    for ( int i = 0; i < lastOverburdenIndex; i++ )
+    {
+        values[i] = computePorBarWithGradient( intersections, values, i, lastOverburdenIndex, gradient );
+    }
+
+    // Fill in underburden values using gradient
+    int lastElementIndex = static_cast<int>( values.size() ) - 1;
+    for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- )
+    {
+        values[i] = computePorBarWithGradient( intersections, values, i, firstUnderburdenIndex, gradient );
+    }
+
+    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

@@ -25,12 +25,22 @@
 
 #include <vector>
 
-class RigFemPartResultsCollection;
+#include "cafPdmField.h"
-class RimGeoMechCase;
+#include "cafPdmPtrField.h"
-class RigGeoMechCaseData;
+
-class RigFemScalarResultFrames;
+#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;
 
 //==================================================================================================
 ///
@@ -39,7 +49,7 @@ class RimWellIADataAccess;
 class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAccessor
 {
 public:
-    RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property );
+    RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property, double gradient );
     ~RimFaultReactivationDataAccessorStress();
 
     bool isMatching( RimFaultReactivation::Property property ) const override;
@@ -48,7 +58,8 @@ public:
                             const RigFaultReactivationModel& model,
                             RimFaultReactivation::GridPart   gridPart,
                             double                           topDepth     = std::numeric_limits<double>::infinity(),
-                            double                           bottomDepth = std::numeric_limits<double>::infinity() ) const override;
+                            double                           bottomDepth  = std::numeric_limits<double>::infinity(),
+                            size_t                           elementIndex = std::numeric_limits<size_t>::max() ) const override;
 
 private:
     void updateResultAccessor() override;
@@ -60,12 +71,44 @@ private:
                                     const cvf::Vec3d&         position,
                                     const std::vector<float>& scalarResults ) const;
 
+    std::pair<double, cvf::Vec3d> getPorBar( RimWellIADataAccess& iaDataAccess,
+                                             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 int                 getPartIndexFromPoint( const RigFemPartCollection& partCollection, const cvf::Vec3d& point );
+    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;
+    double                         m_gradient;
     RigGeoMechCaseData*            m_geoMechCaseData;
     RigFemScalarResultFrames*      m_s11Frames;
     RigFemScalarResultFrames*      m_s22Frames;
     RigFemScalarResultFrames*      m_s33Frames;
-    RigFemScalarResultFrames*      m_porFrames;
     const RigFemPart*              m_femPart;
+
+    cvf::ref<RigWellPath>                m_faceAWellPath;
+    cvf::ref<RigWellPath>                m_faceBWellPath;
+    int                                  m_partIndexA;
+    int                                  m_partIndexB;
+    cvf::ref<RigGeoMechWellLogExtractor> m_extractorA;
+    cvf::ref<RigGeoMechWellLogExtractor> m_extractorB;
 };

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp

@@ -87,7 +87,8 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf::
                                                                      const RigFaultReactivationModel& model,
                                                                      RimFaultReactivation::GridPart   gridPart,
                                                                      double                           topDepth,
-                                                                     double                           bottomDepth ) const
+                                                                     double                           bottomDepth,
+                                                                     size_t                           elementIndex ) const
 {
     if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
     {

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h

@@ -45,7 +45,8 @@ public:
                             const RigFaultReactivationModel& model,
                             RimFaultReactivation::GridPart   gridPart,
                             double                           topDepth     = std::numeric_limits<double>::infinity(),
-                            double                           bottomDepth = std::numeric_limits<double>::infinity() ) const override;
+                            double                           bottomDepth  = std::numeric_limits<double>::infinity(),
+                            size_t                           elementIndex = std::numeric_limits<size_t>::max() ) const override;
 
 private:
     void updateResultAccessor() override;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp

@@ -88,7 +88,8 @@ double RimFaultReactivationDataAccessorVoidRatio::valueAtPosition( const cvf::Ve
                                                                    const RigFaultReactivationModel& model,
                                                                    RimFaultReactivation::GridPart   gridPart,
                                                                    double                           topDepth,
-                                                                   double                           bottomDepth ) const
+                                                                   double                           bottomDepth,
+                                                                   size_t                           elementIndex ) const
 {
     if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
     {

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h

@@ -45,7 +45,8 @@ public:
                             const RigFaultReactivationModel& model,
                             RimFaultReactivation::GridPart   gridPart,
                             double                           topDepth     = std::numeric_limits<double>::infinity(),
-                            double                           bottomDepth = std::numeric_limits<double>::infinity() ) const override;
+                            double                           bottomDepth  = std::numeric_limits<double>::infinity(),
+                            size_t                           elementIndex = std::numeric_limits<size_t>::max() ) const override;
 
 private:
     void updateResultAccessor() override;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp

@@ -709,6 +709,7 @@ bool RimFaultReactivationModel::extractAndExportModelData()
     // extract data for each timestep
     m_dataAccess = std::make_shared<RimFaultReactivationDataAccess>( eCase, geoMechCase(), selectedTimeStepIndexes );
     m_dataAccess->extractModelData( *model() );
+    m_dataAccess.reset();
 
     return true;
 }