#4669 Compute stress gradient.

ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt

@@ -33,9 +33,11 @@ add_library( ${PROJECT_NAME}
 	RigFemResultAddress.h
 	RigFemResultPosEnum.h
 	RimFemResultObserver.h
-    RimFemResultObserver.cpp
-    RimGeoMechGeometrySelectionItem.h
-    RimGeoMechGeometrySelectionItem.cpp
+	RimFemResultObserver.cpp
+	RigHexGradientTools.h
+	RigHexGradientTools.cpp
+	RimGeoMechGeometrySelectionItem.h
+	RimGeoMechGeometrySelectionItem.cpp
 )
 
 target_include_directories(${PROJECT_NAME}

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp

@@ -36,6 +36,7 @@
 #include "RigFemPartResults.h"
 #include "RigFemScalarResultFrames.h"
 #include "RigFormationNames.h"
+#include "RigHexGradientTools.h"
 #include "RigHexIntersectionTools.h"
 #include "RigStatisticsDataCache.h"
 #include "RigWbsParameter.h"
@@ -45,11 +46,10 @@
 #include "RimWellLogPlot.h"
 #include "RimWellLogPlotCollection.h"
 
-#include "cafProgressInfo.h"
-#include "cvfBoundingBox.h"
-
 #include "cafProgressInfo.h"
 #include "cafTensor3.h"
+
+#include "cvfBoundingBox.h"
 #include "cvfGeometryTools.h"
 #include "cvfMath.h"
 
@@ -1008,6 +1008,10 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateStressGradient(
 
     frameCountProgress.incrementProgress();
 
+    const RigFemPart*              femPart      = m_femParts->part( partIndex );
+    int                            elementCount = femPart->elementCount();
+    const std::vector<cvf::Vec3f>& nodeCoords   = femPart->nodes().coordinates;
+
     int frameCount = st11->frameCount();
     for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
     {
@@ -1017,10 +1021,40 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateStressGradient(
         size_t              valCount     = st11Data.size();
         dstFrameData.resize( valCount );
 
-#pragma omp parallel for
-        for ( long vIdx = 0; vIdx < static_cast<long>( valCount ); ++vIdx )
+#pragma omp parallel for schedule( dynamic )
+        for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx )
         {
-            dstFrameData[vIdx] = st11Data[vIdx];
+            const int*     cornerIndices = femPart->connectivities( elmIdx );
+            RigElementType elmType       = femPart->elementType( elmIdx );
+
+            if ( !( elmType == HEX8P || elmType == HEX8 ) ) continue;
+
+            // Find the corner coordinates for element
+            std::array<cvf::Vec3d, 8> hexCorners;
+            for ( int i = 0; i < 8; i++ )
+            {
+                hexCorners[i] = cvf::Vec3d( nodeCoords[cornerIndices[i]] );
+            }
+
+            // Find the corresponding corner values for the element
+            std::array<double, 8> cornerValues;
+
+            int elmNodeCount = RigFemTypes::elmentNodeCount( elmType );
+            for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx )
+            {
+                size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
+                int    nodeIdx      = femPart->nodeIdxFromElementNodeResultIdx( elmNodResIdx );
+
+                cornerValues[elmNodIdx] = st11Data[nodeIdx];
+            }
+
+            std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
+
+            for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx )
+            {
+                size_t elmNodResIdx        = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
+                dstFrameData[elmNodResIdx] = gradients[elmNodIdx].x();
+            }
         }
 
         frameCountProgress.incrementProgress();

ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.cpp

@@ -0,0 +1,77 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2020-     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 "RigHexGradientTools.h"
+
+#include "../cafHexInterpolator/cafHexInterpolator.h"
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::array<cvf::Vec3d, 8> RigHexGradientTools::gradients( const std::array<cvf::Vec3d, 8>& hexCorners,
+                                                          const std::array<double, 8>&     values )
+{
+    //
+    std::array<cvf::Vec3d, 8> gradientsUVW;
+    gradientsUVW[0] = cvf::Vec3d( forwardFD( values, 0, 1 ), forwardFD( values, 0, 3 ), forwardFD( values, 0, 4 ) );
+    gradientsUVW[1] = cvf::Vec3d( backwardFD( values, 1, 0 ), forwardFD( values, 1, 2 ), forwardFD( values, 1, 5 ) );
+    gradientsUVW[2] = cvf::Vec3d( backwardFD( values, 2, 3 ), backwardFD( values, 2, 1 ), forwardFD( values, 2, 6 ) );
+    gradientsUVW[3] = cvf::Vec3d( forwardFD( values, 3, 2 ), backwardFD( values, 3, 0 ), forwardFD( values, 3, 7 ) );
+    gradientsUVW[4] = cvf::Vec3d( forwardFD( values, 4, 5 ), forwardFD( values, 4, 7 ), backwardFD( values, 4, 0 ) );
+    gradientsUVW[5] = cvf::Vec3d( backwardFD( values, 5, 4 ), forwardFD( values, 5, 6 ), backwardFD( values, 5, 1 ) );
+    gradientsUVW[6] = cvf::Vec3d( backwardFD( values, 6, 7 ), backwardFD( values, 6, 5 ), backwardFD( values, 6, 2 ) );
+    gradientsUVW[7] = cvf::Vec3d( forwardFD( values, 7, 6 ), backwardFD( values, 7, 4 ), backwardFD( values, 7, 3 ) );
+
+    std::array<cvf::Vec3d, 8> NC;
+    NC[0] = {-1, -1, -1};
+    NC[1] = {1, -1, -1};
+    NC[2] = {1, 1, -1};
+    NC[3] = {-1, 1, -1};
+    NC[4] = {-1, -1, 1};
+    NC[5] = {1, -1, 1};
+    NC[6] = {1, 1, 1};
+    NC[7] = {-1, 1, 1};
+
+    std::array<cvf::Vec3d, 8> gradientsXYZ;
+
+    for ( int i = 0; i < 8; i++ )
+    {
+        bool       isInvertPossible = false;
+        cvf::Mat3d jacobian         = caf::HexInterpolator::jacobi( hexCorners, NC[i] );
+        gradientsXYZ[i]             = gradientsUVW[i].getTransformedVector( jacobian.getInverted( &isInvertPossible ) );
+    }
+
+    return gradientsXYZ;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RigHexGradientTools::forwardFD( const std::array<double, 8>& values, int from, int to )
+{
+    double h = 2.0;
+    return ( values[to] - values[from] ) / h;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RigHexGradientTools::backwardFD( const std::array<double, 8>& values, int from, int to )
+{
+    return forwardFD( values, to, from );
+}

ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.h

@@ -0,0 +1,43 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2020-     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 "RigHexGradientTools.h"
+
+#include "cvfVector3.h"
+
+#include <array>
+
+//==================================================================================================
+///
+//==================================================================================================
+
+class RigHexGradientTools
+{
+public:
+    static std::array<cvf::Vec3d, 8> gradients( const std::array<cvf::Vec3d, 8>& hexCorners,
+                                                const std::array<double, 8>&     values );
+
+private:
+    // Private to avoid instantiation
+    RigHexGradientTools();
+
+    static double forwardFD( const std::array<double, 8>& values, int from, int to );
+    static double backwardFD( const std::array<double, 8>& values, int from, int to );
+};