#5591 Avoid saw tooth shape on normalized stresses

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.cpp

@@ -20,6 +20,8 @@
 #include "RigFemPartGrid.h"
 
 #include "RigFemPart.h"
+
+#include <array>
 #include <cmath>
 #include <limits.h>
 
@@ -481,8 +483,15 @@ void RigFemPartGrid::cellCornerVertices( size_t cellIndex, cvf::Vec3d vertices[8
 //--------------------------------------------------------------------------------------------------
 cvf::Vec3d RigFemPartGrid::cellCentroid( size_t cellIndex ) const
 {
-    CVF_ASSERT( false );
+    std::array<cvf::Vec3d, 8> cellVertices;
-    return cvf::Vec3d::ZERO;
+    this->cellCornerVertices( cellIndex, cellVertices.data() );
+
+    cvf::Vec3d centroid( 0.0, 0.0, 0.0 );
+    for ( int i = 0; i < 8; ++i )
+    {
+        centroid += cellVertices[i];
+    }
+    return centroid / 8.0;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.h

@@ -42,6 +42,7 @@ public:
 
     std::pair<cvf::Vec3st, cvf::Vec3st> reservoirIJKBoundingBox() const;
     void                                cellCornerVertices( size_t cellIndex, cvf::Vec3d vertices[8] ) const override;
+    cvf::Vec3d                          cellCentroid( size_t cellIndex ) const override;
 
 private:
     void generateStructGridData();
@@ -59,8 +60,7 @@ private: // Unused, Not implemented
     cvf::Vec3d maxCoordinate() const override;
     bool       cellIJKNeighbor( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const override;
 
-    bool       cellIJKFromCoordinate( const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k ) const override;
+    bool cellIJKFromCoordinate( const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k ) const override;
-    cvf::Vec3d cellCentroid( size_t cellIndex ) const override;
     void cellMinMaxCordinates( size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate ) const override;
     size_t     gridPointIndexFromIJK( size_t i, size_t j, size_t k ) const override;
     cvf::Vec3d gridPointCoordinate( size_t i, size_t j, size_t k ) const override;

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp

@@ -1213,14 +1213,25 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNormalizedResult
 {
     CVF_ASSERT( resVarAddr.normalizeByHydrostaticPressure() && isNormalizableResult( resVarAddr ) );
 
-    RigFemResultAddress unscaledResult             = resVarAddr;
+    RigFemResultAddress unscaledResult = resVarAddr;
+    if ( unscaledResult.resultPosType == RIG_NODAL && unscaledResult.fieldName == "POR-Bar" )
+        unscaledResult.resultPosType = RIG_ELEMENT_NODAL;
     unscaledResult.normalizedByHydrostaticPressure = false;
 
-    caf::ProgressInfo frameCountProgress( this->frameCount() * 3, "Calculating Normalized Result" );
+    CAF_ASSERT( unscaledResult.resultPosType == RIG_ELEMENT_NODAL );
 
+    caf::ProgressInfo frameCountProgress( this->frameCount() * 4, "Calculating Normalized Result" );
+
+    RigFemScalarResultFrames* porDataFrames = nullptr;
     RigFemScalarResultFrames* srcDataFrames = nullptr;
     RigFemScalarResultFrames* dstDataFrames = nullptr;
 
+    {
+        auto task     = frameCountProgress.task( "Loading POR Result", this->frameCount() );
+        porDataFrames = this->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_ELEMENT_NODAL, "POR-Bar", "" ) );
+        if ( !porDataFrames ) return nullptr;
+    }
+
     {
         auto task     = frameCountProgress.task( "Loading Unscaled Result", this->frameCount() );
         srcDataFrames = this->findOrLoadScalarResult( partIndex, unscaledResult );
@@ -1235,7 +1246,9 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNormalizedResult
     frameCountProgress.setProgressDescription( "Normalizing Result" );
     frameCountProgress.setNextProgressIncrement( 1u );
 
-    const RigFemPart*              femPart      = m_femParts->part( partIndex );
+    const RigFemPart*     femPart     = m_femParts->part( partIndex );
+    const RigFemPartGrid* femPartGrid = femPart->getOrCreateStructGrid();
+
     float                          inf          = std::numeric_limits<float>::infinity();
     int                            elmNodeCount = femPart->elementCount();
     const std::vector<cvf::Vec3f>& nodeCoords   = femPart->nodes().coordinates;
@@ -1243,36 +1256,60 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNormalizedResult
     int frameCount = srcDataFrames->frameCount();
     for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
     {
+        const std::vector<float>& porFrameData = porDataFrames->frameData( fIdx );
         const std::vector<float>& srcFrameData = srcDataFrames->frameData( fIdx );
         std::vector<float>&       dstFrameData = dstDataFrames->frameData( fIdx );
 
         size_t resultCount = srcFrameData.size();
         dstFrameData.resize( resultCount );
 
-        if ( resVarAddr.resultPosType == RIG_ELEMENT_NODAL )
+        if ( unscaledResult.resultPosType == RIG_ELEMENT_NODAL )
         {
 #pragma omp parallel for schedule( dynamic )
-            for ( int elmNodeResIdx = 0; elmNodeResIdx < resultCount; ++elmNodeResIdx )
+            for ( int elmIdx = 0; elmIdx < femPart->elementCount(); ++elmIdx )
             {
-                const int        nodeIdx             = femPart->nodeIdxFromElementNodeResultIdx( elmNodeResIdx );
+                RigElementType elmType = femPart->elementType( elmIdx );
-                const cvf::Vec3f node                = nodeCoords[nodeIdx];
+                if ( !( elmType == HEX8 || elmType == HEX8P ) ) continue;
-                double           tvdRKB              = std::abs( node.z() ) + m_normalizationAirGap;
+
-                double           hydrostaticPressure = RiaWellLogUnitTools::hydrostaticPorePressureBar( tvdRKB );
+                bool porRegion = false;
-                dstFrameData[elmNodeResIdx]          = srcFrameData[elmNodeResIdx] / hydrostaticPressure;
+                for ( int elmLocalNodeIdx = 0; elmLocalNodeIdx < 8; ++elmLocalNodeIdx )
+                {
+                    size_t    elmNodeResIdx     = femPart->elementNodeResultIdx( elmIdx, elmLocalNodeIdx );
+                    const int nodeIdx           = femPart->nodeIdxFromElementNodeResultIdx( elmNodeResIdx );
+                    dstFrameData[elmNodeResIdx] = srcFrameData[elmNodeResIdx];
+                    if ( porFrameData[elmNodeResIdx] != std::numeric_limits<float>::infinity() )
+                    {
+                        porRegion = true;
+                    }
+                }
+                if ( porRegion )
+                {
+                    // This is in the POR-region. Use hydrostatic pressure from the individual nodes
+                    for ( int elmLocalNodeIdx = 0; elmLocalNodeIdx < 8; ++elmLocalNodeIdx )
+                    {
+                        size_t    elmNodeResIdx       = femPart->elementNodeResultIdx( elmIdx, elmLocalNodeIdx );
+                        const int nodeIdx             = femPart->nodeIdxFromElementNodeResultIdx( elmNodeResIdx );
+                        double    tvdRKB              = std::abs( nodeCoords[nodeIdx].z() ) + m_normalizationAirGap;
+                        double    hydrostaticPressure = RiaWellLogUnitTools::hydrostaticPorePressureBar( tvdRKB );
+                        dstFrameData[elmNodeResIdx] /= hydrostaticPressure;
+                    }
+                }
+                else
+                {
+                    // Over/under/sideburden. Use hydrostatic pressure from cell centroid.
+                    cvf::Vec3d cellCentroid       = femPartGrid->cellCentroid( elmIdx );
+                    double     cellCentroidTvdRKB = std::abs( cellCentroid.z() ) + m_normalizationAirGap;
+                    double     cellCenterHydroStaticPressure =
+                        RiaWellLogUnitTools::hydrostaticPorePressureBar( cellCentroidTvdRKB );
+
+                    for ( int elmLocalNodeIdx = 0; elmLocalNodeIdx < 8; ++elmLocalNodeIdx )
+                    {
+                        size_t elmNodeResIdx = femPart->elementNodeResultIdx( elmIdx, elmLocalNodeIdx );
+                        dstFrameData[elmNodeResIdx] /= cellCenterHydroStaticPressure;
+                    }
+                }
             }
         }
-        else if ( resVarAddr.resultPosType == RIG_NODAL )
-        {
-#pragma omp parallel for schedule( dynamic )
-            for ( int nodeResIdx = 0; nodeResIdx < resultCount; ++nodeResIdx )
-            {
-                const cvf::Vec3f node                = nodeCoords[nodeResIdx];
-                double           tvdRKB              = std::abs( node.z() ) + m_normalizationAirGap;
-                double           hydrostaticPressure = RiaWellLogUnitTools::hydrostaticPorePressureBar( tvdRKB );
-                dstFrameData[nodeResIdx]             = srcFrameData[nodeResIdx] / hydrostaticPressure;
-            }
-        }
-        frameCountProgress.incrementProgress();
     }
     return dstDataFrames;
 }
@@ -1655,7 +1692,8 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAngles( i
     int                            frameCount      = this->frameCount();
 
     // HACK ! Todo : make it robust against other elements than Hex8
-    size_t valCount = femPart->elementCount() * 24; // Number of Elm Node Face results 24 = 4 * num faces = 3* numElmNodes
+    size_t valCount = femPart->elementCount() * 24; // Number of Elm Node Face results 24 = 4 * num faces = 3*
+                                                    // numElmNodes
 
     for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
     {
@@ -1693,8 +1731,9 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAngles( i
                     cvf::Mat3f rotMx = cvf::GeometryTools::computePlaneHorizontalRotationMx( quadVxs[2] - quadVxs[0],
                                                                                              quadVxs[3] - quadVxs[1] );
                     RiaOffshoreSphericalCoords sphCoord(
-                        cvf::Vec3f( rotMx.rowCol( 0, 2 ), rotMx.rowCol( 1, 2 ), rotMx.rowCol( 2, 2 ) ) ); // Use Ez from
+                        cvf::Vec3f( rotMx.rowCol( 0, 2 ), rotMx.rowCol( 1, 2 ), rotMx.rowCol( 2, 2 ) ) ); // Use Ez
-                                                                                                          // the matrix
+                                                                                                          // from the
+                                                                                                          // matrix
                                                                                                           // as plane
                                                                                                           // normal
 
@@ -2654,7 +2693,8 @@ std::vector<RigFemResultAddress>
     if ( fcIt != fieldAndComponentNames.end() )
     {
         std::vector<std::string> compNames = fcIt->second;
-        if ( !resVarAddr.componentName.empty() ) // If we did not request a particular component, do not add the components
+        if ( !resVarAddr.componentName.empty() ) // If we did not request a particular component, do not add the
+                                                 // components
         {
             for ( const auto& compName : compNames )
             {

ApplicationCode/ProjectDataModel/RimGeoMechCase.cpp

@@ -295,6 +295,7 @@ RimGeoMechCase::CaseOpenStatus RimGeoMechCase::openGeoMechCase( std::string* err
         fileNames.push_back( fileName.path() );
     }
     geoMechCaseData->femPartResults()->addElementPropertyFiles( fileNames );
+    geoMechCaseData->femPartResults()->setCalculationParameters( m_cohesion, cvf::Math::toRadians( m_frictionAngleDeg() ) );
 
     m_geoMechCaseData = geoMechCaseData;
 

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp

@@ -114,7 +114,7 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition( void )
                        "",
                        "",
                        "" );
-    CAF_PDM_InitField( &m_normalizationAirGap, "NormalizationAirGap", -1.0, "Air Gap", "", "", "" );
+    CAF_PDM_InitField( &m_normalizationAirGap, "NormalizationAirGap", 0.0, "Air Gap", "", "", "" );
     m_normalizationAirGap.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleValueEditor::uiEditorTypeName() );
 
     CAF_PDM_InitField( &m_compactionRefLayerUiField,
@@ -338,33 +338,13 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
             m_resultVariableUiField = "";
         }
     }
-    if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationAirGap < 0.0 )
+    if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationAirGap == 0.0 )
     {
-        RiaMedianCalculator<double> airGapCalc;
+        calculateNormalizationAirGapDefault();
-        for ( auto wellPath : RiaApplication::instance()->project()->allWellPaths() )
+        if ( m_normalizeByHydrostaticPressure )
         {
-            if ( wellPath->wellPathGeometry() )
+            m_geomCase->geoMechData()->femPartResults()->setNormalizationAirGap( m_normalizationAirGap );
-            {
-                double airGap = wellPath->wellPathGeometry()->rkbDiff();
-                if ( airGap > 0.0 )
-                {
-                    airGapCalc.add( airGap );
-                }
-            }
         }
-        if ( airGapCalc.valid() )
-        {
-            m_normalizationAirGap = airGapCalc.median();
-        }
-        else
-        {
-            m_normalizationAirGap = 0.0;
-        }
-    }
-
-    if ( m_normalizeByHydrostaticPressure && m_normalizationAirGap >= 0.0 )
-    {
-        m_geomCase->geoMechData()->femPartResults()->setNormalizationAirGap( m_normalizationAirGap );
     }
 
     // Get the possible property filter owner
@@ -456,6 +436,33 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimGeoMechResultDefinition::calculateNormalizationAirGapDefault()
+{
+    RiaMedianCalculator<double> airGapCalc;
+    for ( auto wellPath : RiaApplication::instance()->project()->allWellPaths() )
+    {
+        if ( wellPath->wellPathGeometry() )
+        {
+            double airGap = wellPath->wellPathGeometry()->rkbDiff();
+            if ( airGap > 0.0 )
+            {
+                airGapCalc.add( airGap );
+            }
+        }
+    }
+    if ( airGapCalc.valid() )
+    {
+        m_normalizationAirGap = airGapCalc.median();
+    }
+    else
+    {
+        m_normalizationAirGap = 0.0;
+    }
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -798,6 +805,22 @@ QString RimGeoMechResultDefinition::defaultLasUnits() const
     return units;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimGeoMechResultDefinition::normalizationAirGap() const
+{
+    return m_normalizationAirGap;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimGeoMechResultDefinition::setNormalizationAirGap( double airGap )
+{
+    m_normalizationAirGap = airGap;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h

@@ -71,6 +71,9 @@ public:
     QString currentResultUnits() const;
     QString defaultLasUnits() const;
 
+    double normalizationAirGap() const;
+    void   setNormalizationAirGap( double airGap );
+
     bool hasCategoryResult() const { return m_resultPositionType() == RIG_FORMATION_NAMES; }
 
     void updateLegendTextAndRanges( RimRegularLegendConfig* legendConfigToUpdate,
@@ -89,10 +92,13 @@ private:
     QList<caf::PdmOptionItemInfo> calculateValueOptions( const caf::PdmFieldHandle* fieldNeedingOptions,
                                                          bool*                      useOptionsOnly ) override;
     void                          fieldChangedByUi( const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue ) override;
-    void                          initAfterRead() override;
+
-    void                          defineEditorAttribute( const caf::PdmFieldHandle* field,
+    void calculateNormalizationAirGapDefault();
-                                                         QString                    uiConfigName,
+
-                                                         caf::PdmUiEditorAttribute* attribute ) override;
+    void initAfterRead() override;
+    void defineEditorAttribute( const caf::PdmFieldHandle* field,
+                                QString                    uiConfigName,
+                                caf::PdmUiEditorAttribute* attribute ) override;
     // Metadata and option build tools
 
     std::map<std::string, std::vector<std::string>> getResultMetaDataForUIFieldSetting();

ApplicationCode/ProjectDataModel/RimGeoMechView.cpp

@@ -169,6 +169,7 @@ void RimGeoMechView::onLoadDataAndUpdate()
     CVF_ASSERT( this->cellResult() != nullptr );
     if ( this->hasUserRequestedAnimation() )
     {
+        m_geomechCase->geoMechData()->femPartResults()->setNormalizationAirGap( this->cellResult()->normalizationAirGap() );
         m_geomechCase->geoMechData()->femPartResults()->assertResultsLoaded( this->cellResult()->resultAddress() );
     }
     progress.incrementProgress();

ApplicationCode/ProjectDataModel/RimWellLogExtractionCurve.cpp

@@ -222,6 +222,7 @@ void RimWellLogExtractionCurve::setPropertiesFromView( Rim3dView* view )
     if ( geoMechView )
     {
         m_geomResultDefinition->setResultAddress( geoMechView->cellResultResultDefinition()->resultAddress() );
+        m_geomResultDefinition->setNormalizationAirGap( geoMechView->cellResultResultDefinition()->normalizationAirGap() );
         m_timeStep = geoMechView->currentTimeStep();
     }
     else if ( geomCase )