#11034 Fault Reactivation: integrate vertical stress from eclipse grid.

ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake

@@ -12,6 +12,7 @@ set(SOURCE_GROUP_HEADER_FILES
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorGeoMech.h
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStress.h
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStressGeoMech.h
+    ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStressEclipse.h
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorWellLogExtraction.h
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationEnums.h
 )
@@ -30,6 +31,7 @@ set(SOURCE_GROUP_SOURCE_FILES
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorGeoMech.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStress.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStressGeoMech.cpp
+    ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStressEclipse.cpp
     ${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorWellLogExtraction.cpp
 )
 

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp

@@ -35,6 +35,7 @@
 #include "RimFaultReactivationDataAccessor.h"
 #include "RimFaultReactivationDataAccessorGeoMech.h"
 #include "RimFaultReactivationDataAccessorPorePressure.h"
+#include "RimFaultReactivationDataAccessorStressEclipse.h"
 #include "RimFaultReactivationDataAccessorStressGeoMech.h"
 #include "RimFaultReactivationDataAccessorTemperature.h"
 #include "RimFaultReactivationDataAccessorVoidRatio.h"
@@ -56,6 +57,14 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( const RimFaultRe
     m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorPorePressure>( thecase, porePressureGradient, seabedDepth ) );
     m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorVoidRatio>( thecase, 0.0001 ) );
     m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorTemperature>( thecase, topTemperature, seabedDepth ) );
+
+    std::vector<RimFaultReactivation::Property> stressProperties = { RimFaultReactivation::Property::StressTop,
+                                                                     RimFaultReactivation::Property::DepthTop,
+                                                                     RimFaultReactivation::Property::StressBottom,
+                                                                     RimFaultReactivation::Property::DepthBottom,
+                                                                     RimFaultReactivation::Property::LateralStressComponentX,
+                                                                     RimFaultReactivation::Property::LateralStressComponentY };
+
     if ( geoMechCase )
     {
         std::vector<RimFaultReactivation::Property> properties = { RimFaultReactivation::Property::YoungsModulus,
@@ -67,18 +76,40 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( const RimFaultRe
             m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorGeoMech>( geoMechCase, property ) );
         }
 
-        std::vector<RimFaultReactivation::Property> stressProperties = { RimFaultReactivation::Property::StressTop,
-                                                                         RimFaultReactivation::Property::DepthTop,
-                                                                         RimFaultReactivation::Property::StressBottom,
-                                                                         RimFaultReactivation::Property::DepthBottom,
-                                                                         RimFaultReactivation::Property::LateralStressComponentX,
-                                                                         RimFaultReactivation::Property::LateralStressComponentY };
         for ( auto property : stressProperties )
         {
             m_accessors.push_back(
                 std::make_shared<RimFaultReactivationDataAccessorStressGeoMech>( geoMechCase, property, porePressureGradient, seabedDepth ) );
         }
     }
+    else
+    {
+        auto extractDensities = []( const RimFaultReactivationModel& model )
+        {
+            std::map<RimFaultReactivation::ElementSets, double> densities;
+            std::vector<RimFaultReactivation::ElementSets>      elementSets = { RimFaultReactivation::ElementSets::OverBurden,
+                                                                                RimFaultReactivation::ElementSets::UnderBurden,
+                                                                                RimFaultReactivation::ElementSets::Reservoir,
+                                                                                RimFaultReactivation::ElementSets::IntraReservoir };
+            for ( auto e : elementSets )
+            {
+                densities[e] = model.materialParameters( e )[2];
+            }
+            return densities;
+        };
+
+        std::map<RimFaultReactivation::ElementSets, double> densities = extractDensities( model );
+        for ( auto property : stressProperties )
+        {
+            m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorStressEclipse>( thecase,
+                                                                                                    property,
+                                                                                                    porePressureGradient,
+                                                                                                    seabedDepth,
+                                                                                                    model.lateralStressCoefficientX(),
+                                                                                                    model.lateralStressCoefficientY(),
+                                                                                                    densities ) );
+        }
+    }
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp

@@ -84,20 +84,20 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
     cvf::Vec3d topPosition( position.x(), position.y(), topDepth );
     cvf::Vec3d bottomPosition( position.x(), position.y(), bottomDepth );
 
-    if ( isPositionValid( position, topPosition, bottomPosition ) )
+    if ( isPositionValid( position, topPosition, bottomPosition, gridPart ) )
     {
         if ( m_property == RimFaultReactivation::Property::StressTop )
         {
-            auto [porBar, extractionPos] = calculatePorBar( topPosition, m_gradient );
+            auto [porBar, extractionPos] = calculatePorBar( topPosition, m_gradient, gridPart );
             if ( std::isinf( porBar ) ) return porBar;
-            double s33 = extractStressValue( StressType::S33, extractionPos );
+            double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
             return RiaEclipseUnitTools::barToPascal( s33 - porBar );
         }
         else if ( m_property == RimFaultReactivation::Property::StressBottom )
         {
-            auto [porBar, extractionPos] = calculatePorBar( bottomPosition, m_gradient );
+            auto [porBar, extractionPos] = calculatePorBar( bottomPosition, m_gradient, gridPart );
             if ( std::isinf( porBar ) ) return porBar;
-            double s33 = extractStressValue( StressType::S33, extractionPos );
+            double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
             return RiaEclipseUnitTools::barToPascal( s33 - porBar );
         }
         else if ( m_property == RimFaultReactivation::Property::DepthTop )
@@ -110,21 +110,37 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
         }
         else if ( m_property == RimFaultReactivation::Property::LateralStressComponentX )
         {
-            auto [porBar, extractionPos] = calculatePorBar( position, m_gradient );
+            return lateralStressComponentX( position, gridPart );
-            if ( std::isinf( porBar ) ) return porBar;
-            double s11 = extractStressValue( StressType::S11, extractionPos );
-            double s33 = extractStressValue( StressType::S33, extractionPos );
-            return ( s11 - porBar ) / ( s33 - porBar );
         }
         else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY )
         {
-            auto [porBar, extractionPos] = calculatePorBar( position, m_gradient );
+            return lateralStressComponentY( position, gridPart );
-            if ( std::isinf( porBar ) ) return porBar;
-            double s22 = extractStressValue( StressType::S22, extractionPos );
-            double s33 = extractStressValue( StressType::S33, extractionPos );
-            return ( s22 - porBar ) / ( s33 - porBar );
         }
     }
 
     return std::numeric_limits<double>::infinity();
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorStress::lateralStressComponentX( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const
+{
+    auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, gridPart );
+    if ( std::isinf( porBar ) ) return porBar;
+    double s11 = extractStressValue( StressType::S11, extractionPos, gridPart );
+    double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
+    return ( s11 - porBar ) / ( s33 - porBar );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorStress::lateralStressComponentY( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const
+{
+    auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, gridPart );
+    if ( std::isinf( porBar ) ) return porBar;
+    double s22 = extractStressValue( StressType::S22, extractionPos, gridPart );
+    double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
+    return ( s22 - porBar ) / ( s33 - porBar );
+}

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h

@@ -61,11 +61,18 @@ public:
 protected:
     virtual bool isDataAvailable() const = 0;
 
-    virtual double extractStressValue( StressType stressType, const cvf::Vec3d& position ) const = 0;
+    virtual double extractStressValue( StressType stressType, const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const = 0;
 
-    virtual std::pair<double, cvf::Vec3d> calculatePorBar( const cvf::Vec3d& position, double gradient ) const = 0;
+    virtual std::pair<double, cvf::Vec3d>
+        calculatePorBar( const cvf::Vec3d& position, double gradient, RimFaultReactivation::GridPart gridPart ) const = 0;
 
-    virtual bool isPositionValid( const cvf::Vec3d& position, const cvf::Vec3d& topPosition, const cvf::Vec3d& bottomPosition ) const = 0;
+    virtual bool isPositionValid( const cvf::Vec3d&              position,
+                                  const cvf::Vec3d&              topPosition,
+                                  const cvf::Vec3d&              bottomPosition,
+                                  RimFaultReactivation::GridPart gridPart ) const = 0;
+
+    virtual double lateralStressComponentX( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const;
+    virtual double lateralStressComponentY( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const;
 
     RimFaultReactivation::Property m_property;
     double                         m_gradient;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp

@@ -0,0 +1,305 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  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 "RimFaultReactivationDataAccessorStressEclipse.h"
+
+#include "RiaEclipseUnitTools.h"
+#include "RiaInterpolationTools.h"
+#include "RiaLogging.h"
+#include "RiaWellLogUnitTools.h"
+
+#include "RigCaseCellResultsData.h"
+#include "RigEclipseResultAddress.h"
+#include "RigEclipseWellLogExtractor.h"
+#include "RigFaultReactivationModel.h"
+#include "RigGriddedPart3d.h"
+#include "RigResultAccessorFactory.h"
+#include "RigWellPath.h"
+
+#include "RimEclipseCase.h"
+#include "RimFaultReactivationDataAccessorStress.h"
+#include "RimFaultReactivationDataAccessorWellLogExtraction.h"
+#include "RimFaultReactivationEnums.h"
+
+#include "cvfObject.h"
+#include "cvfVector3.h"
+
+#include <cmath>
+#include <limits>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimFaultReactivationDataAccessorStressEclipse::RimFaultReactivationDataAccessorStressEclipse(
+    RimEclipseCase*                                            eclipseCase,
+    RimFaultReactivation::Property                             property,
+    double                                                     gradient,
+    double                                                     seabedDepth,
+    double                                                     lateralStressComponentX,
+    double                                                     lateralStressComponentY,
+    const std::map<RimFaultReactivation::ElementSets, double>& densities )
+    : RimFaultReactivationDataAccessorStress( property, gradient, seabedDepth )
+    , m_eclipseCase( eclipseCase )
+    , m_caseData( nullptr )
+    , m_mainGrid( nullptr )
+    , m_lateralStressComponentX( lateralStressComponentX )
+    , m_lateralStressComponentY( lateralStressComponentY )
+    , m_densities( densities )
+{
+    if ( m_eclipseCase )
+    {
+        m_caseData = m_eclipseCase->eclipseCaseData();
+        m_mainGrid = m_eclipseCase->mainGrid();
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimFaultReactivationDataAccessorStressEclipse::~RimFaultReactivationDataAccessorStressEclipse()
+{
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimFaultReactivationDataAccessorStressEclipse::updateResultAccessor()
+{
+    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, 0, RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress );
+
+    auto [wellPaths, extractors] =
+        RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData, m_seabedDepth );
+    m_wellPaths  = wellPaths;
+    m_extractors = extractors;
+
+    // TODO: get from model
+    double waterDensity = 1030.0;
+
+    for ( auto [gridPart, wellPath] : m_wellPaths )
+    {
+        auto                    extractor     = m_extractors[gridPart];
+        std::vector<cvf::Vec3d> intersections = extractor->intersections();
+
+        addOverburdenAndUnderburdenPoints( intersections, wellPath->wellPathPoints() );
+
+        m_stressValues[gridPart] =
+            integrateVerticalStress( *wellPath.p(), intersections, *m_model, gridPart, m_seabedDepth, waterDensity, m_densities );
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimFaultReactivationDataAccessorStressEclipse::addOverburdenAndUnderburdenPoints( std::vector<cvf::Vec3d>&       intersections,
+                                                                                       const std::vector<cvf::Vec3d>& wellPathPoints )
+{
+    // Insert points at top of overburden and under underburden
+    intersections.insert( intersections.begin(), wellPathPoints.front() );
+    intersections.push_back( wellPathPoints.back() );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RimFaultReactivationDataAccessorStressEclipse::isDataAvailable() const
+{
+    return m_mainGrid != nullptr && m_resultAccessor.notNull();
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorStressEclipse::extractStressValue( StressType                     stressType,
+                                                                          const cvf::Vec3d&              position,
+                                                                          RimFaultReactivation::GridPart gridPart ) const
+{
+    CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() );
+    auto extractor = m_extractors.find( gridPart )->second;
+
+    CAF_ASSERT( m_stressValues.find( gridPart ) != m_stressValues.end() );
+    auto stressValues = m_stressValues.find( gridPart )->second;
+
+    CAF_ASSERT( m_wellPaths.find( gridPart ) != m_wellPaths.end() );
+    auto wellPath = m_wellPaths.find( gridPart )->second;
+
+    auto intersections = extractor->intersections();
+    addOverburdenAndUnderburdenPoints( intersections, wellPath->wellPathPoints() );
+
+    CAF_ASSERT( stressValues.size() == intersections.size() );
+
+    auto [topIdx, bottomIdx] = RimFaultReactivationDataAccessorWellLogExtraction::findIntersectionsForTvd( intersections, position.z() );
+    if ( topIdx != -1 && bottomIdx != -1 )
+    {
+        double topValue    = stressValues[topIdx];
+        double bottomValue = stressValues[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 = { stressValues[bottomIdx], stressValues[topIdx] };
+            return RiaEclipseUnitTools::pascalToBar( RiaInterpolationTools::linear( xs, ys, position.z() ) );
+        }
+    }
+    else if ( position.z() <= intersections.back().z() )
+    {
+        return RiaEclipseUnitTools::pascalToBar( stressValues.back() );
+    }
+
+    return std::numeric_limits<double>::infinity();
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RimFaultReactivationDataAccessorStressEclipse::isPositionValid( const cvf::Vec3d&              position,
+                                                                     const cvf::Vec3d&              topPosition,
+                                                                     const cvf::Vec3d&              bottomPosition,
+                                                                     RimFaultReactivation::GridPart gridPart ) const
+{
+    auto [porBar, extractionPosition] = calculatePorBar( position, m_gradient, gridPart );
+    return !std::isinf( porBar ) && extractionPosition != cvf::Vec3d::UNDEFINED;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressEclipse::calculatePorBar( const cvf::Vec3d& position,
+                                                                                              double            gradient,
+                                                                                              RimFaultReactivation::GridPart gridPart ) const
+{
+    if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
+    {
+        CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() );
+        auto extractor = m_extractors.find( gridPart )->second;
+
+        CAF_ASSERT( m_wellPaths.find( gridPart ) != m_wellPaths.end() );
+        auto wellPath = m_wellPaths.find( gridPart )->second;
+
+        auto [values, intersections] =
+            RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );
+
+        auto [value, pos] = RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( intersections, values, position, m_gradient );
+        return { value, pos };
+    }
+
+    return { std::numeric_limits<double>::infinity(), cvf::Vec3d::UNDEFINED };
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<double>
+    RimFaultReactivationDataAccessorStressEclipse::integrateVerticalStress( const RigWellPath&               wellPath,
+                                                                            const std::vector<cvf::Vec3d>&   intersections,
+                                                                            const RigFaultReactivationModel& model,
+                                                                            RimFaultReactivation::GridPart   gridPart,
+                                                                            double                           seabedDepth,
+                                                                            double                           waterDensity,
+                                                                            const std::map<RimFaultReactivation::ElementSets, double>& densities )
+{
+    double              gravity      = RiaWellLogUnitTools<double>::gravityAcceleration();
+    double              seaWaterLoad = gravity * std::abs( seabedDepth ) * waterDensity;
+    std::vector<double> values       = { seaWaterLoad };
+
+    auto part = model.grid( gridPart );
+    CAF_ASSERT( part );
+    auto elementSets = part->elementSets();
+
+    double previousDensity = densities.find( RimFaultReactivation::ElementSets::OverBurden )->second;
+
+    for ( size_t i = 1; i < intersections.size(); i++ )
+    {
+        double previousValue = values[i - 1];
+        double previousDepth = intersections[i - 1].z();
+        double currentDepth  = intersections[i].z();
+
+        double deltaDepth = previousDepth - currentDepth;
+        double density    = previousDensity;
+
+        auto [isOk, elementSet] = findElementSetForPoint( *part, intersections[i], elementSets );
+        if ( isOk )
+        {
+            // Unit: kg/m^3
+            CAF_ASSERT( densities.find( elementSet ) != densities.end() );
+            density = densities.find( elementSet )->second;
+        }
+
+        double value = previousValue + density * gravity * deltaDepth;
+        values.push_back( value );
+
+        previousDensity = density;
+    }
+
+    return values;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<bool, RimFaultReactivation::ElementSets> RimFaultReactivationDataAccessorStressEclipse::findElementSetForPoint(
+    const RigGriddedPart3d&                                                       part,
+    const cvf::Vec3d&                                                             point,
+    const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets )
+{
+    for ( auto [elementSet, elements] : elementSets )
+    {
+        for ( unsigned int elementIndex : elements )
+        {
+            // Find nodes from element
+            auto positions = part.elementCorners( elementIndex );
+
+            std::array<cvf::Vec3d, 8> coordinates;
+            for ( size_t i = 0; i < positions.size(); ++i )
+            {
+                coordinates[i] = positions[i];
+            }
+
+            if ( RigHexIntersectionTools::isPointInCell( point, coordinates.data() ) )
+            {
+                return { true, elementSet };
+            }
+        }
+    }
+
+    return { false, RimFaultReactivation::ElementSets::OverBurden };
+};
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorStressEclipse::lateralStressComponentX( const cvf::Vec3d&              position,
+                                                                               RimFaultReactivation::GridPart gridPart ) const
+{
+    return m_lateralStressComponentX;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorStressEclipse::lateralStressComponentY( const cvf::Vec3d&              position,
+                                                                               RimFaultReactivation::GridPart gridPart ) const
+{
+    return m_lateralStressComponentY;
+}

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h

@@ -0,0 +1,99 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  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 "RimFaultReactivationDataAccessorStress.h"
+#include "RimFaultReactivationEnums.h"
+
+#include "cafPdmField.h"
+#include "cafPdmPtrField.h"
+
+#include "cvfObject.h"
+
+#include <vector>
+
+class RigGriddedPart3d;
+class RimModeledWellPath;
+class RigWellPath;
+class RimEclipseCase;
+class RigEclipseCaseData;
+class RigResultAccessor;
+class RigEclipseWellLogExtractor;
+
+//==================================================================================================
+///
+///
+//==================================================================================================
+class RimFaultReactivationDataAccessorStressEclipse : public RimFaultReactivationDataAccessorStress
+{
+public:
+    RimFaultReactivationDataAccessorStressEclipse( RimEclipseCase*                                            geoMechCase,
+                                                   RimFaultReactivation::Property                             property,
+                                                   double                                                     gradient,
+                                                   double                                                     seabedDepth,
+                                                   double                                                     lateralStressComponentX,
+                                                   double                                                     lateralStressComponentY,
+                                                   const std::map<RimFaultReactivation::ElementSets, double>& densities );
+    ~RimFaultReactivationDataAccessorStressEclipse() override;
+
+private:
+    void updateResultAccessor() override;
+
+    bool isDataAvailable() const override;
+
+    double extractStressValue( StressType stressType, const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const override;
+
+    std::pair<double, cvf::Vec3d>
+        calculatePorBar( const cvf::Vec3d& position, double gradient, RimFaultReactivation::GridPart gridPart ) const override;
+
+    bool isPositionValid( const cvf::Vec3d&              position,
+                          const cvf::Vec3d&              topPosition,
+                          const cvf::Vec3d&              bottomPosition,
+                          RimFaultReactivation::GridPart gridPart ) const override;
+
+    double lateralStressComponentX( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const override;
+    double lateralStressComponentY( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const override;
+
+    static std::vector<double> integrateVerticalStress( const RigWellPath&                                         wellPath,
+                                                        const std::vector<cvf::Vec3d>&                             intersections,
+                                                        const RigFaultReactivationModel&                           model,
+                                                        RimFaultReactivation::GridPart                             gridPart,
+                                                        double                                                     seabedDepth,
+                                                        double                                                     waterDensity,
+                                                        const std::map<RimFaultReactivation::ElementSets, double>& densities );
+
+    static void addOverburdenAndUnderburdenPoints( std::vector<cvf::Vec3d>& intersections, const std::vector<cvf::Vec3d>& wellPathPoints );
+
+    static std::pair<bool, RimFaultReactivation::ElementSets>
+        findElementSetForPoint( const RigGriddedPart3d&                                                       part,
+                                const cvf::Vec3d&                                                             point,
+                                const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets );
+
+    RimEclipseCase*             m_eclipseCase;
+    RigEclipseCaseData*         m_caseData;
+    const RigMainGrid*          m_mainGrid;
+    cvf::ref<RigResultAccessor> m_resultAccessor;
+    double                      m_lateralStressComponentX;
+    double                      m_lateralStressComponentY;
+
+    std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>                m_wellPaths;
+    std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>> m_extractors;
+    std::map<RimFaultReactivation::GridPart, std::vector<double>>                  m_stressValues;
+    std::map<RimFaultReactivation::ElementSets, double>                            m_densities;
+};

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp

@@ -90,8 +90,11 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor()
     m_s22Frames       = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ), timeStepIndex );
 
     auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom();
-    auto   faultNormal                           = m_model->faultNormal();
+    auto faultNormal                             = m_model->faultNormal() ^ cvf::Vec3d::Z_AXIS;
+    faultNormal.normalize();
+
     double distanceFromFault     = 1.0;
+    auto [topDepth, bottomDepth] = m_model->depthTopBottom();
 
     RigFemPartCollection* geoMechPartCollection = m_geoMechCaseData->femParts();
     std::string           errorName             = "fault reactivation data access";
@@ -101,6 +104,7 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor()
             RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
                                                                                    faultBottomPosition,
                                                                                    m_seabedDepth,
+                                                                                   bottomDepth,
                                                                                    faultNormal * distanceFromFault );
         m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
         m_partIndexA    = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] );
@@ -112,6 +116,7 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor()
             RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
                                                                                    faultBottomPosition,
                                                                                    m_seabedDepth,
+                                                                                   bottomDepth,
                                                                                    -faultNormal * distanceFromFault );
         m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
         m_partIndexB    = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] );
@@ -139,7 +144,9 @@ bool RimFaultReactivationDataAccessorStressGeoMech::isDataAvailable() const
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-double RimFaultReactivationDataAccessorStressGeoMech::extractStressValue( StressType stressType, const cvf::Vec3d& position ) const
+double RimFaultReactivationDataAccessorStressGeoMech::extractStressValue( StressType                     stressType,
+                                                                          const cvf::Vec3d&              position,
+                                                                          RimFaultReactivation::GridPart gridPart ) const
 {
     RimWellIADataAccess iaDataAccess( m_geoMechCase );
 
@@ -170,7 +177,9 @@ RigFemScalarResultFrames* RimFaultReactivationDataAccessorStressGeoMech::dataFra
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position, double gradient ) const
+std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position,
+                                                                                              double            gradient,
+                                                                                              RimFaultReactivation::GridPart gridPart ) const
 {
     int timeStepIndex = 0;
     int frameIndex    = m_s33Frames->frameCount( timeStepIndex ) - 1;
@@ -182,7 +191,8 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::cal
 //--------------------------------------------------------------------------------------------------
 bool RimFaultReactivationDataAccessorStressGeoMech::isPositionValid( const cvf::Vec3d&              position,
                                                                      const cvf::Vec3d&              topPosition,
-                                                                     const cvf::Vec3d& bottomPosition ) const
+                                                                     const cvf::Vec3d&              bottomPosition,
+                                                                     RimFaultReactivation::GridPart gridPart ) const
 {
     RimWellIADataAccess iaDataAccess( m_geoMechCase );
     int                 centerElementIdx = iaDataAccess.elementIndex( position );

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h

@@ -60,11 +60,15 @@ private:
 
     bool isDataAvailable() const override;
 
-    double extractStressValue( StressType stressType, const cvf::Vec3d& position ) const override;
+    double extractStressValue( StressType stressType, const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const override;
 
-    std::pair<double, cvf::Vec3d> calculatePorBar( const cvf::Vec3d& position, double gradient ) const override;
+    std::pair<double, cvf::Vec3d>
+        calculatePorBar( const cvf::Vec3d& position, double gradient, RimFaultReactivation::GridPart gridPart ) const override;
 
-    bool isPositionValid( const cvf::Vec3d& position, const cvf::Vec3d& topPosition, const cvf::Vec3d& bottomPosition ) const override;
+    bool isPositionValid( const cvf::Vec3d&              position,
+                          const cvf::Vec3d&              topPosition,
+                          const cvf::Vec3d&              bottomPosition,
+                          RimFaultReactivation::GridPart gridPart ) const override;
 
     static RigFemResultAddress getResultAddress( const std::string& fieldName, const std::string& componentName );
 

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp

@@ -115,6 +115,10 @@ std::pair<double, cvf::Vec3d>
             return { porBar, extractionPosition };
         }
     }
+    else if ( position.z() <= intersections.back().z() )
+    {
+        return { values.back(), intersections.back() };
+    }
 
     return { std::numeric_limits<double>::infinity(), cvf::Vec3d::UNDEFINED };
 }
@@ -244,12 +248,13 @@ void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithT
 std::vector<cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( const cvf::Vec3d& faultTopPosition,
                                                                                                const cvf::Vec3d& faultBottomPosition,
                                                                                                double            seabedDepth,
+                                                                                               double            bottomDepth,
                                                                                                const cvf::Vec3d& offset )
 {
     cvf::Vec3d faultTop = faultTopPosition + offset;
     cvf::Vec3d seabed( faultTop.x(), faultTop.y(), seabedDepth );
     cvf::Vec3d faultBottom = faultBottomPosition + offset;
-    cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), -10000.0 );
+    cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), bottomDepth );
     return { seabed, faultTop, faultBottom, underburdenBottom };
 }
 
@@ -280,19 +285,23 @@ std::pair<std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>, std::
                                                                                         double                           seabedDepth )
 {
     auto [faultTopPosition, faultBottomPosition] = model.faultTopBottom();
-    auto   faultNormal                           = model.faultNormal();
+    auto faultNormal                             = model.faultNormal() ^ cvf::Vec3d::Z_AXIS;
+    faultNormal.normalize();
+
     double distanceFromFault     = 1.0;
+    auto [topDepth, bottomDepth] = model.depthTopBottom();
 
     std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>                wellPaths;
     std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>> extractors;
 
     for ( auto gridPart : model.allGridParts() )
     {
-        double                  sign = model.normalPointsAt() == gridPart ? 1.0 : -1.0;
+        double                  sign = model.normalPointsAt() == gridPart ? -1.0 : 1.0;
         std::vector<cvf::Vec3d> wellPoints =
             RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
                                                                                    faultBottomPosition,
                                                                                    seabedDepth,
+                                                                                   bottomDepth,
                                                                                    sign * faultNormal * distanceFromFault );
         cvf::ref<RigWellPath> wellPath =
             new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
@@ -347,7 +356,7 @@ std::vector<double> RimFaultReactivationDataAccessorWellLogExtraction::extractDe
     std::vector<double> intersectionsZ;
     for ( auto i : intersections )
     {
-        intersectionsZ.push_back( i.z() );
+        intersectionsZ.push_back( -i.z() );
     }
     return intersectionsZ;
 }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h

@@ -54,6 +54,7 @@ public:
     static std::vector<cvf::Vec3d> generateWellPoints( const cvf::Vec3d& faultTopPosition,
                                                        const cvf::Vec3d& faultBottomPosition,
                                                        double            seabedDepth,
+                                                       double            bottomDepth,
                                                        const cvf::Vec3d& offset );
 
     static std::vector<double> generateMds( const std::vector<cvf::Vec3d>& points );
@@ -61,9 +62,9 @@ public:
     static std::pair<std::vector<double>, std::vector<cvf::Vec3d>> extractValuesAndIntersections( const RigResultAccessor& resultAccessor,
                                                                                                   RigEclipseWellLogExtractor& extractor,
                                                                                                   const RigWellPath&          wellPath );
+    static std::pair<int, int> findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd );
 
 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              computeValueWithGradient( const std::vector<cvf::Vec3d>& intersections,
                                                          const std::vector<double>&     values,

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp

@@ -130,6 +130,9 @@ RimFaultReactivationModel::RimFaultReactivationModel()
     CAF_PDM_InitField( &m_frictionAngleDeg, "FrictionAngle", 20.0, "Friction Angle [degree]" );
     CAF_PDM_InitField( &m_seabedTemperature, "SeabedTemperature", 5.0, "Seabed Temperature [C]" );
 
+    CAF_PDM_InitField( &m_lateralStressCoefficientX, "LateralStressCoefficientX", 0.5, "Lateral Stress Coeff. X" );
+    CAF_PDM_InitField( &m_lateralStressCoefficientY, "LateralStressCoefficientY", 0.5, "Lateral Stress Coeff. Y" );
+
     CAF_PDM_InitFieldNoDefault( &m_targets, "Targets", "Targets" );
     m_targets.uiCapability()->setUiEditorTypeName( caf::PdmUiTableViewEditor::uiEditorTypeName() );
     m_targets.uiCapability()->setUiTreeChildrenHidden( true );
@@ -486,15 +489,20 @@ void RimFaultReactivationModel::defineUiOrdering( QString uiConfigName, caf::Pdm
     {
         propertiesGrp->add( &m_useGridDensity );
         propertiesGrp->add( &m_useGridElasticProperties );
-        propertiesGrp->add( &m_useGridStress );
         propertiesGrp->add( &m_waterDensity );
-
+        propertiesGrp->add( &m_useGridStress );
         propertiesGrp->add( &m_seabedTemperature );
 
         bool useTemperatureFromGrid = m_useGridTemperature();
         m_seabedTemperature.uiCapability()->setUiReadOnly( !useTemperatureFromGrid );
     }
 
+    if ( !hasGeoMechCase || !m_useGridStress() )
+    {
+        propertiesGrp->add( &m_lateralStressCoefficientX );
+        propertiesGrp->add( &m_lateralStressCoefficientY );
+    }
+
     propertiesGrp->add( &m_frictionAngleDeg );
 
     uiOrdering.skipRemainingFields();
@@ -811,3 +819,19 @@ double RimFaultReactivationModel::seabedTemperature() const
 {
     return m_seabedTemperature;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationModel::lateralStressCoefficientX() const
+{
+    return m_lateralStressCoefficientX;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationModel::lateralStressCoefficientY() const
+{
+    return m_lateralStressCoefficientY;
+}

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.h

@@ -129,6 +129,8 @@ public:
     double waterDensity() const;
     double frictionAngleDeg() const;
     double seabedTemperature() const;
+    double lateralStressCoefficientX() const;
+    double lateralStressCoefficientY() const;
 
     RimEclipseCase* eclipseCase() const;
     RimGeoMechCase* geoMechCase() const;
@@ -185,6 +187,8 @@ private:
     caf::PdmField<double> m_waterDensity;
     caf::PdmField<double> m_frictionAngleDeg;
     caf::PdmField<double> m_seabedTemperature;
+    caf::PdmField<double> m_lateralStressCoefficientX;
+    caf::PdmField<double> m_lateralStressCoefficientY;
 
     caf::PdmField<size_t> m_startCellIndex;
     caf::PdmField<int>    m_startCellFace;

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp

@@ -21,6 +21,8 @@
 #include "RigFaultReactivationModelGenerator.h"
 #include "RigGriddedPart3d.h"
 
+#include <limits>
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -226,6 +228,16 @@ const std::pair<cvf::Vec3d, cvf::Vec3d> RigFaultReactivationModel::faultTopBotto
     return m_generator->faultTopBottomPoints();
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, double> RigFaultReactivationModel::depthTopBottom() const
+{
+    if ( m_generator.get() == nullptr )
+        return std::make_pair( std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity() );
+    return m_generator->depthTopBottom();
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.h

@@ -84,6 +84,7 @@ public:
 
     const cvf::Vec3d                        faultNormal() const;
     const std::pair<cvf::Vec3d, cvf::Vec3d> faultTopBottom() const;
+    std::pair<double, double>               depthTopBottom() const;
 
     RimFaultReactivation::GridPart normalPointsAt() const;
 

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp

@@ -42,6 +42,7 @@ RigFaultReactivationModelGenerator::RigFaultReactivationModelGenerator( cvf::Vec
     , m_bufferAboveFault( 0.0 )
     , m_bufferBelowFault( 0.0 )
     , m_startDepth( 0.0 )
+    , m_bottomDepth( 0.0 )
     , m_depthBelowFault( 100.0 )
     , m_horzExtentFromFault( 1000.0 )
     , m_modelThickness( 100.0 )
@@ -295,13 +296,13 @@ void RigFaultReactivationModelGenerator::generatePointsFrontBack()
     alongModel.normalize();
 
     double top_depth = -m_startDepth;
-    double bottom_depth = m_bottomFault.z() - m_depthBelowFault;
+    m_bottomDepth    = m_bottomFault.z() - m_depthBelowFault;
 
     cvf::Vec3d edge_front = m_startPosition - m_horzExtentFromFault * alongModel;
     cvf::Vec3d edge_back  = m_startPosition + m_horzExtentFromFault * alongModel;
 
     points[8]     = m_bottomFault;
-    points[8].z() = bottom_depth;
+    points[8].z() = m_bottomDepth;
 
     points[9]  = m_bottomFault;
     points[10] = m_bottomReservoirBack;
@@ -724,3 +725,11 @@ const std::pair<cvf::Vec3d, cvf::Vec3d> RigFaultReactivationModelGenerator::faul
 {
     return std::make_pair( m_topFault, m_bottomFault );
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, double> RigFaultReactivationModelGenerator::depthTopBottom() const
+{
+    return { -m_startDepth, m_bottomDepth };
+}

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.h

@@ -66,6 +66,7 @@ public:
     const std::array<cvf::Vec3d, 12>&       backPoints() const;
     const cvf::Vec3d                        normal() const;
     const std::pair<cvf::Vec3d, cvf::Vec3d> faultTopBottomPoints() const;
+    std::pair<double, double>               depthTopBottom() const;
 
 protected:
     static const std::array<int, 4>      faceIJCornerIndexes( cvf::StructGridInterface::FaceType face );
@@ -103,6 +104,7 @@ private:
     double m_bufferBelowFault;
 
     double m_startDepth;
+    double m_bottomDepth;
     double m_depthBelowFault;
     double m_horzExtentFromFault;
     double m_modelThickness;