Fix extraction egde cases.

OPM · Jan 17, 2024 · 761a12a · 761a12a
1 parent 70588df
commit 761a12a
Show file tree

Hide file tree

Showing 4 changed files with 68 additions and 45 deletions.
diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.cpp
@@ -162,6 +162,10 @@ double RimFaultReactivationDataAccessorStressEclipse::extractStressValue( Stress
             return RiaEclipseUnitTools::pascalToBar( RiaInterpolationTools::linear( xs, ys, position.z() ) );
         }
     }
+    else if ( position.z() <= intersections.back().z() )
+    {
+        return RiaEclipseUnitTools::pascalToBar( stressValues.back() );
+    }
 
     return std::numeric_limits<double>::infinity();
 }
@@ -219,35 +223,11 @@ std::vector<double>
     double              seaWaterLoad = gravity * std::abs( seabedDepth ) * waterDensity;
     std::vector<double> values       = { seaWaterLoad };
 
-    auto g           = model.grid( gridPart );
-    auto elementSets = g->elementSets();
-
-    auto getElementSet =
-        []( auto part, const cvf::Vec3d& point, const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets )
-        -> std::pair<bool, RimFaultReactivation::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 };
-                }
-            }
-        }
+    auto part = model.grid( gridPart );
+    CAF_ASSERT( part );
+    auto elementSets = part->elementSets();
 
-        return { false, RimFaultReactivation::ElementSets::OverBurden };
-    };
+    double previousDensity = densities.find( RimFaultReactivation::ElementSets::OverBurden )->second;
 
     for ( size_t i = 1; i < intersections.size(); i++ )
     {
@@ -256,26 +236,56 @@ std::vector<double>
         double currentDepth  = intersections[i].z();
 
         double deltaDepth = previousDepth - currentDepth;
+        double density    = previousDensity;
 
-        auto [isOk, elementSet] = getElementSet( g, intersections[i], elementSets );
+        auto [isOk, elementSet] = findElementSetForPoint( *part, intersections[i], elementSets );
         if ( isOk )
         {
             // Unit: kg/m^3
             CAF_ASSERT( densities.find( elementSet ) != densities.end() );
-            double density = densities.find( elementSet )->second;
-
-            double value = previousValue + density * gravity * deltaDepth;
-            values.push_back( value );
-        }
-        else
-        {
-            values.push_back( std::numeric_limits<double>::infinity() );
+            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 };
+};
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressEclipse.h
@@ -80,6 +80,11 @@ class RimFaultReactivationDataAccessorStressEclipse : public RimFaultReactivatio
 
     static void addOverburdenAndUnderburdenPoints( std::vector<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;

diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp
@@ -90,9 +90,11 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor()
     m_s22Frames       = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ), timeStepIndex );
 
     auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom();
-    auto   faultNormal                           = m_model->faultNormal();
-    double distanceFromFault                     = 1.0;
-    auto [topDepth, bottomDepth]                 = m_model->depthTopBottom();
+    auto faultNormal                             = m_model->faultNormal() ^ cvf::Vec3d::Z_AXIS;
+    faultNormal.normalize();
+
+    double distanceFromFault     = 1.0;
+    auto [topDepth, bottomDepth] = m_model->depthTopBottom();
 
     RigFemPartCollection* geoMechPartCollection = m_geoMechCaseData->femParts();
     std::string           errorName             = "fault reactivation data access";

diff --git a/...tionLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp b/...tionLibCode/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 };
 }
@@ -281,16 +285,18 @@ std::pair<std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>, std::
                                                                                         double                           seabedDepth )
 {
     auto [faultTopPosition, faultBottomPosition] = model.faultTopBottom();
-    auto   faultNormal                           = model.faultNormal();
-    double distanceFromFault                     = 1.0;
-    auto [topDepth, bottomDepth]                 = model.depthTopBottom();
+    auto faultNormal                             = model.faultNormal() ^ cvf::Vec3d::Z_AXIS;
+    faultNormal.normalize();
+
+    double distanceFromFault     = 1.0;
+    auto [topDepth, bottomDepth] = model.depthTopBottom();
 
     std::map<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,
@@ -350,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;
 }