#2680. First reasonable Fracture Gradient and Shear Failure Gradient estimation.

Caveats:
* Hard coded poissonRatio = 0.25
* Hard coded UCS to 100 bar (fairly close to average value for shale in literature).

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.cpp

@@ -98,14 +98,16 @@ void RigFemPart::assertNodeToElmIndicesIsCalculated()
 void RigFemPart::calculateNodeToElmRefs()
 {
     m_nodeToElmRefs.resize(nodes().nodeIds.size());
+    m_nodeGlobalToLocalIndices.resize(nodes().nodeIds.size());
 
     for (int eIdx = 0; eIdx < static_cast<int>(m_elementId.size()); ++eIdx)
     {
         int elmNodeCount = RigFemTypes::elmentNodeCount(elementType(eIdx));
         const int* elmNodes = connectivities(eIdx);
-        for (int enIdx = 0; enIdx < elmNodeCount; ++enIdx)
+        for (int localIdx = 0; localIdx < elmNodeCount; ++localIdx)
         {
-            m_nodeToElmRefs[elmNodes[enIdx]].push_back(eIdx);
+            m_nodeToElmRefs[elmNodes[localIdx]].push_back(eIdx);
+            m_nodeGlobalToLocalIndices[elmNodes[localIdx]].push_back(localIdx);
         }
     }
 }
@@ -113,17 +115,17 @@ void RigFemPart::calculateNodeToElmRefs()
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-const int* RigFemPart::elementsUsingNode(int nodeIndex)
+const std::vector<int>& RigFemPart::elementsUsingNode(int nodeIndex) const
 {
-   return &(m_nodeToElmRefs[nodeIndex][0]);
+   return m_nodeToElmRefs[nodeIndex];
 }
 
 //--------------------------------------------------------------------------------------------------
-/// 
+///
 //--------------------------------------------------------------------------------------------------
-int RigFemPart::numElementsUsingNode(int nodeIndex)
+const std::vector<unsigned char>& RigFemPart::elementLocalIndicesForNode(int nodeIndex) const
 {
-    return static_cast<int>(m_nodeToElmRefs[nodeIndex].size());
+    return m_nodeGlobalToLocalIndices[nodeIndex];
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -168,16 +170,15 @@ void RigFemPart::calculateElmNeighbors()
             candidates.clear();
             {
                 int firstNodeIdxOfFace = elmNodes[localFaceIndices[0]];
-                int candidateCount1 = this->numElementsUsingNode(firstNodeIdxOfFace);
+                const std::vector<int>& candidates1 = this->elementsUsingNode(firstNodeIdxOfFace);
-                const int* candidates1 = this->elementsUsingNode(firstNodeIdxOfFace);
 
-                if (candidateCount1)
+                if (!candidates1.empty())
                 {
                     // Get neighbor candidates from the diagonal node
 
                     int thirdNodeIdxOfFace = elmNodes[localFaceIndices[3]];
-                    int candidateCount2 = this->numElementsUsingNode(thirdNodeIdxOfFace);
+                    
-                    const int* candidates2 = this->elementsUsingNode(thirdNodeIdxOfFace);
+                    const std::vector<int>& candidates2 = this->elementsUsingNode(thirdNodeIdxOfFace);
 
                     // The candidates are sorted from smallest to largest, so we do a linear search to find the 
                     // (two) common cells in the two arrays, and leaving this element out, we have one candidate left
@@ -185,7 +186,7 @@ void RigFemPart::calculateElmNeighbors()
                     int idx1 = 0;
                     int idx2 = 0;
 
-                    while (idx1 < candidateCount1 && idx2 < candidateCount2)
+                    while (idx1 < candidates1.size() && idx2 < candidates2.size())
                     {
                         if (candidates1[idx1] < candidates2[idx2]){ ++idx1; continue; }
                         if (candidates1[idx1] > candidates2[idx2]){ ++idx2; continue; }

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.h

@@ -67,9 +67,9 @@ public:
     RigFemPartNodes&            nodes()                                    {return m_nodes;}
     const RigFemPartNodes&      nodes() const                              {return m_nodes;}
 
-    void                        assertNodeToElmIndicesIsCalculated();
+    void                              assertNodeToElmIndicesIsCalculated();
-    const int*                  elementsUsingNode(int nodeIndex);
+    const std::vector<int>&           elementsUsingNode(int nodeIndex) const;
-    int                         numElementsUsingNode(int nodeIndex);
+    const std::vector<unsigned char>& elementLocalIndicesForNode(int nodeIndex) const;
     
     void                        assertElmNeighborsIsCalculated();
     int                         elementNeighbor(int elementIndex, int faceIndex) const
@@ -100,12 +100,13 @@ private:
     mutable cvf::ref<RigFemPartGrid>    m_structGrid;
 
     void calculateNodeToElmRefs();
-    std::vector<std::vector<int> > m_nodeToElmRefs; // Needs a more memory friendly structure
+    std::vector<std::vector<int>>           m_nodeToElmRefs; // Needs a more memory friendly structure
+    std::vector<std::vector<unsigned char>> m_nodeGlobalToLocalIndices;
   
     void calculateElmNeighbors();
     struct Neighbors { int indicesToNeighborElms[6]; char faceInNeighborElm[6];};
-    std::vector<  Neighbors > m_elmNeighbors;
+    std::vector<  Neighbors >   m_elmNeighbors;
-    std::vector<int> m_possibleGridCornerElements;
+    std::vector<int>            m_possibleGridCornerElements;
 
     mutable float               m_characteristicElementSize;
     mutable cvf::BoundingBox    m_boundingBox;

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp

@@ -482,6 +482,11 @@ std::map<std::string, std::vector<std::string> > RigFemPartResultsCollection::sc
                 fieldCompNames[field];
             }
         }
+        else if (resPos == RIG_WELLPATH_DERIVED)
+        {
+            fieldCompNames["FractureGradient"];
+            fieldCompNames["ShearFailureGradient"];
+        }
     }
 
     return fieldCompNames;

ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultPosEnum.h

@@ -25,5 +25,6 @@ enum RigFemResultPosEnum {
     RIG_INTEGRATION_POINT,
     RIG_ELEMENT_NODAL_FACE,
     RIG_FORMATION_NAMES,
-    RIG_ELEMENT
+    RIG_ELEMENT,
+    RIG_WELLPATH_DERIVED
 };

ApplicationCode/ProjectDataModel/Rim3dWellLogExtractionCurve.cpp

@@ -28,6 +28,7 @@
 #include "RigResultAccessorFactory.h"
 #include "RigCaseCellResultsData.h"
 #include "RigFemPartResultsCollection.h"
+#include "RigWellPath.h"
 #include "RimEclipseCase.h"
 #include "RimGeoMechCase.h"
 #include "Rim3dView.h"
@@ -79,6 +80,7 @@ Rim3dWellLogExtractionCurve::Rim3dWellLogExtractionCurve()
     m_geomResultDefinition.uiCapability()->setUiHidden(true);
     m_geomResultDefinition.uiCapability()->setUiTreeChildrenHidden(true);
     m_geomResultDefinition = new RimGeoMechResultDefinition;
+    m_geomResultDefinition->setAddWellPathDerivedResults(true);
 
     CAF_PDM_InitFieldNoDefault(&m_nameConfig, "NameConfig", "", "", "", "");
     m_nameConfig = new RimWellLogExtractionCurveNameConfig(this);    
@@ -199,10 +201,9 @@ void Rim3dWellLogExtractionCurve::curveValuesAndMds(std::vector<double>* values,
     else if (geomExtractor.notNull())
     {
         *measuredDepthValues = geomExtractor->measuredDepth();
-
         m_geomResultDefinition->loadResult();
-
+        geomExtractor->setRkbDiff(rkbDiff());
-        geomExtractor->curveData(m_geomResultDefinition->resultAddress(), m_timeStep, values);
+        geomExtractor->curveData(m_geomResultDefinition->resultAddress(), m_timeStep, values);        
     }
 }
 
@@ -280,6 +281,41 @@ QString Rim3dWellLogExtractionCurve::createCurveAutoName() const
     return generatedCurveName.join(", ");
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+double Rim3dWellLogExtractionCurve::rkbDiff() const
+{
+    RimWellPath* wellPath;
+    firstAncestorOrThisOfType(wellPath);
+
+    if (wellPath && wellPath->wellPathGeometry())
+    {
+        RigWellPath* geo = wellPath->wellPathGeometry();
+
+        if (geo->hasDatumElevation())
+        {
+            return geo->datumElevation();
+        }
+
+        // If measured depth is zero, use the z-value of the well path points
+        if (geo->m_wellPathPoints.size() > 0 && geo->m_measuredDepths.size() > 0)
+        {
+            double epsilon = 1e-3;
+
+            if (cvf::Math::abs(geo->m_measuredDepths[0]) < epsilon)
+            {
+                double diff = geo->m_measuredDepths[0] - (-geo->wellPathPoints()[0].z());
+
+                return diff;
+            }
+        }
+    }
+
+    return HUGE_VAL;
+}
+
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -293,7 +329,7 @@ caf::PdmFieldHandle* Rim3dWellLogExtractionCurve::userDescriptionField()
 //--------------------------------------------------------------------------------------------------
 void Rim3dWellLogExtractionCurve::fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue)
 {   
-    if (changedField == &m_timeStep || changedField == &m_case)
+    if (changedField == &m_case)
     {
         this->resetMinMaxValuesAndUpdateUI();
     }

ApplicationCode/ProjectDataModel/Rim3dWellLogExtractionCurve.h

@@ -48,6 +48,7 @@ public:
 
     virtual QString name() const override;
     virtual QString createCurveAutoName() const override;
+    double          rkbDiff() const;
 protected:
     virtual caf::PdmFieldHandle*            userDescriptionField() override;
     virtual void                            fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue) override;

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp

@@ -51,6 +51,7 @@ void caf::AppEnum< RigFemResultPosEnum >::setUp()
     addItem(RIG_ELEMENT_NODAL_FACE, "ELEMENT_NODAL_FACE", "Element Nodal On Face");
     addItem(RIG_FORMATION_NAMES, "FORMATION_NAMES", "Formation Names");
     addItem(RIG_ELEMENT, "ELEMENT", "Element");
+    addItem(RIG_WELLPATH_DERIVED, "WELLPATH_DERIVED", "Well Path Derived");
     setDefault(RIG_NODAL);
 }
 }
@@ -108,6 +109,7 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition(void)
     m_compactionRefLayerUiField.xmlCapability()->setIOReadable(false);
 
     m_isChangedByField = false;
+    m_addWellPathDerivedResults = false;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -169,7 +171,19 @@ QList<caf::PdmOptionItemInfo> RimGeoMechResultDefinition::calculateValueOptions(
 
     if (m_geomCase)
     {
-        if (&m_resultVariableUiField == fieldNeedingOptions)
+        if (&m_resultPositionTypeUiField == fieldNeedingOptions)
+        {
+            std::vector<RigFemResultPosEnum> optionItems = { RIG_NODAL, RIG_ELEMENT_NODAL, RIG_INTEGRATION_POINT, RIG_ELEMENT_NODAL_FACE, RIG_FORMATION_NAMES, RIG_ELEMENT };
+            if (true || m_addWellPathDerivedResults)
+            {
+                optionItems.push_back(RIG_WELLPATH_DERIVED);
+            }
+            for (RigFemResultPosEnum value : optionItems)
+            {
+                options.push_back(caf::PdmOptionItemInfo(caf::AppEnum<RigFemResultPosEnum>::uiText(value), QVariant(value)));
+            }
+        }
+        else if (&m_resultVariableUiField == fieldNeedingOptions)
         {
             std::map<std::string, std::vector<std::string> >  fieldCompNames = getResultMetaDataForUIFieldSetting();
 
@@ -439,10 +453,29 @@ void RimGeoMechResultDefinition::loadResult()
 {
     if (m_geomCase && m_geomCase->geoMechData())
     {
-        m_geomCase->geoMechData()->femPartResults()->assertResultsLoaded(this->resultAddress());
+        if (this->resultAddress().fieldName == "FractureGradient" ||
+            this->resultAddress().fieldName == "ShearFailureGradient")
+        {
+            RigFemResultAddress stressResAddr(RIG_ELEMENT_NODAL, std::string("ST"), "");
+            RigFemResultAddress porBarResAddr(RIG_ELEMENT_NODAL, std::string("POR-Bar"), "");
+            m_geomCase->geoMechData()->femPartResults()->assertResultsLoaded(stressResAddr);
+            m_geomCase->geoMechData()->femPartResults()->assertResultsLoaded(porBarResAddr);
+        }
+        else
+        {
+            m_geomCase->geoMechData()->femPartResults()->assertResultsLoaded(this->resultAddress());
+        }
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimGeoMechResultDefinition::setAddWellPathDerivedResults(bool addWellPathDerivedResults)
+{
+    m_addWellPathDerivedResults = addWellPathDerivedResults;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h

@@ -48,6 +48,7 @@ public:
     RigGeoMechCaseData*                               ownerCaseData();
     bool                                              hasResult(); 
     void                                              loadResult();
+    void                                              setAddWellPathDerivedResults(bool addWellPathDerivedResults);
 
     RigFemResultAddress                               resultAddress();
 
@@ -117,4 +118,5 @@ private:
     caf::PdmPointer<RimGeoMechCase>                   m_geomCase;
 
     bool                                              m_isChangedByField;
+    bool                                              m_addWellPathDerivedResults;
 };

ApplicationCode/ProjectDataModel/RimWellLogExtractionCurve.cpp

@@ -111,6 +111,7 @@ RimWellLogExtractionCurve::RimWellLogExtractionCurve()
     m_geomResultDefinition.uiCapability()->setUiHidden(true);
     m_geomResultDefinition.uiCapability()->setUiTreeChildrenHidden(true);
     m_geomResultDefinition = new RimGeoMechResultDefinition;
+    m_geomResultDefinition->setAddWellPathDerivedResults(true);
 
     CAF_PDM_InitField(&m_timeStep, "CurveTimeStep", 0,"Time Step", "", "", "");
 
@@ -388,7 +389,7 @@ void RimWellLogExtractionCurve::onLoadDataAndUpdate(bool updateParentPlot)
             tvDepthValues = geomExtractor->trueVerticalDepth();
 
             m_geomResultDefinition->loadResult();
-
+            geomExtractor->setRkbDiff(rkbDiff());
             geomExtractor->curveData(m_geomResultDefinition->resultAddress(), m_timeStep, &values);
         }
 

ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake

@@ -66,6 +66,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RigFractureCell.h
 ${CMAKE_CURRENT_LIST_DIR}/RigWellResultPoint.h
 ${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryTools.h
 ${CMAKE_CURRENT_LIST_DIR}/RigCaseRealizationParameters.h
+${CMAKE_CURRENT_LIST_DIR}/RigGeoMechBoreHoleStressCalculator.h
 )
 
 
@@ -131,6 +132,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RigFractureCell.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigWellResultPoint.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryTools.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigCaseRealizationParameters.cpp
+${CMAKE_CURRENT_LIST_DIR}/RigGeoMechBoreHoleStressCalculator.cpp
 )
 
 list(APPEND CODE_HEADER_FILES

ApplicationCode/ReservoirDataModel/RigGeoMechBoreHoleStressCalculator.cpp

@@ -0,0 +1,265 @@
+#include "RigGeoMechBoreHoleStressCalculator.h"
+
+//==================================================================================================
+/// Internal root finding class to find a Well Pressure that gives:
+/// a) a zero SigmaT for estimating the fracture gradient.
+/// b) a solution to the Stassi-d'Alia failure criterion for estimating the shear failure gradient.
+//==================================================================================================
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RigGeoMechBoreHoleStressCalculator::RigGeoMechBoreHoleStressCalculator(const caf::Ten3d& tensor,
+                                                                       double            porePressure,
+                                                                       double            poissonRatio,
+                                                                       double            uniaxialCompressiveStrength,
+                                                                       int               nThetaSubSamples)
+    : m_tensor(tensor)
+    , m_porePressure(porePressure)
+    , m_poissonRatio(poissonRatio)
+    , m_uniaxialCompressiveStrength(uniaxialCompressiveStrength)
+    , m_nThetaSubSamples(nThetaSubSamples)
+{
+    calculateStressComponents();
+}
+
+
+//--------------------------------------------------------------------------------------------------
+/// Simple bisection method for now
+//--------------------------------------------------------------------------------------------------
+double RigGeoMechBoreHoleStressCalculator::solveFractureGradient(double* thetaOut)
+{
+    MemberFunc fn = &RigGeoMechBoreHoleStressCalculator::sigmaTMinOfMin;
+    return solveSecant(fn, thetaOut);
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RigGeoMechBoreHoleStressCalculator::solveStassiDalia(double* thetaOut)
+{
+    MemberFunc fn = &RigGeoMechBoreHoleStressCalculator::stassiDalia;
+    return solveSecant(fn, thetaOut);
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::Vec3d RigGeoMechBoreHoleStressCalculator::principleStressesAtWall(double pw, double theta) const
+{
+    cvf::Vec4d stressComponentsForAngle = calculateStressComponentsForSegmentAngle(theta);
+    double sigma_theta = stressComponentsForAngle[1] - pw;
+    const double& sigma_z = stressComponentsForAngle[2];
+    double tauSqrx4 = std::pow(stressComponentsForAngle[3], 2) * 4.0;
+
+    double sigmaComponent1 = sigma_z + sigma_theta;
+    double sigmaComponent2 = std::sqrt(std::pow(sigma_z - sigma_theta, 2) + tauSqrx4);    
+    return cvf::Vec3d(pw - m_porePressure, 0.5 * (sigmaComponent1 + sigmaComponent2) - m_porePressure, 0.5 * (sigmaComponent1 - sigmaComponent2) - m_porePressure);
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Bi-section root finding method: https://en.wikipedia.org/wiki/Bisection_method
+/// Used as fall-back in case the secant method doesn't converge.
+//--------------------------------------------------------------------------------------------------
+double RigGeoMechBoreHoleStressCalculator::solveBisection(double minPw, double maxPw, MemberFunc fn, double* thetaOut)
+{
+    const int N = 50;
+    const double epsilon = 1.0e-10;
+
+    double theta = 0.0;
+
+    std::pair<double, double> largestNegativeValue(0.0, -std::numeric_limits<double>::infinity());
+    std::pair<double, double> smallestPositiveValue (0.0, std::numeric_limits<double>::infinity());
+        
+    for (int i = 0; i <= N; ++i)
+    {
+        double pw = minPw + (maxPw - minPw) * i / static_cast<double>(N);
+        double f_pw = std::invoke(fn, this, pw, &theta);
+        if (f_pw >= 0.0 && f_pw < smallestPositiveValue.second)
+        {
+            smallestPositiveValue = std::make_pair(pw, f_pw);
+        }
+        if (f_pw < 0.0 && f_pw > largestNegativeValue.second)
+        {
+            largestNegativeValue = std::make_pair(pw, f_pw);
+        }        
+    }
+
+    // TODO: Provide a warning if there was no solution to the equation
+    if (largestNegativeValue.second == -std::numeric_limits<double>::infinity())
+    {
+        // No solution. Function is always positive. Pick smallest value.
+        return smallestPositiveValue.first;
+    }
+    if (smallestPositiveValue.second == std::numeric_limits<double>::infinity())
+    {
+        // No solution. Function is always negative. Pick largest value.
+        return largestNegativeValue.first;        
+    }
+
+    double minPwFuncVal = std::invoke(fn, this, largestNegativeValue.first, &theta);
+    double maxPwFuncVal = std::invoke(fn, this, smallestPositiveValue.first, &theta);
+    double range = maxPw - minPw;
+    
+    int i = 0;
+    for (; i <= N && range > m_porePressure * epsilon; ++i)
+    {
+        double midPw = (minPw + maxPw) * 0.5;
+        double midPwFuncVal = std::invoke(fn, this, midPw, &theta);
+        if (midPwFuncVal * minPwFuncVal < 0.0)
+        {
+            maxPw = midPw;
+            maxPwFuncVal = midPwFuncVal;
+        }
+        else
+        {
+            minPw = midPw;
+            minPwFuncVal = midPwFuncVal;
+        }
+        range = maxPw - minPw;
+    }
+    CVF_ASSERT(i < N); // Otherwise it hasn't converged                       
+
+    if (thetaOut)
+    {
+        *thetaOut = theta;
+    }
+    
+    // Return average of minPw and maxPw.
+    return 0.5 * (maxPw + minPw);
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Secant root finding method: https://en.wikipedia.org/wiki/Secant_method
+/// Basically a Newton's method using finite differences for the derivative.
+//--------------------------------------------------------------------------------------------------
+double RigGeoMechBoreHoleStressCalculator::solveSecant(MemberFunc fn, double* thetaOut)
+{
+    const double epsilon = 1.0e-10;
+    const int N = 50;
+    double theta = 0.0;
+
+    double x_0 = 0.0;    
+    double f_x0 = std::invoke(fn, this, x_0, &theta);
+    double x_1 = m_porePressure;
+    double f_x1 = std::invoke(fn, this, x_1, &theta);
+    double x = 0.0;
+    double f_x = 0.0;
+    int i = 0;
+    for (; i < N && std::abs(f_x1 - f_x0) > epsilon; ++i)
+    {
+        x = x_1 - f_x1 * (x_1 - x_0) / (f_x1 - f_x0);
+        f_x = std::invoke(fn, this, x, &theta);
+        if (std::abs(f_x) < epsilon * m_porePressure) break;
+
+        // Update iteration variables
+        x_0 = x_1;
+        f_x0 = f_x1;
+        x_1 = x;
+        f_x1 = f_x;        
+    }
+
+    if (i == N)
+    {
+        // Fallback to bisection if secant doesn't converge.
+        return solveBisection(0.0, m_porePressure * 2.0, fn, thetaOut);
+    }
+
+    if (thetaOut)
+    {
+        *thetaOut = theta;
+    }
+    return x;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RigGeoMechBoreHoleStressCalculator::sigmaTMinOfMin(double wellPressure, double* thetaAtMin) const
+{
+    CVF_ASSERT(thetaAtMin);
+    double sigma_t_min_min = std::numeric_limits<double>::max();
+    for (const cvf::Vec4d& stressComponentsForAngle : m_stressComponents)
+    {
+        // Perform all these internal calculations in double to reduce significance errors
+        double sigma_theta = stressComponentsForAngle[1] - wellPressure;
+        const double& sigma_z = stressComponentsForAngle[2];
+        double tauSqrx4 = std::pow(stressComponentsForAngle[3], 2) * 4.0;
+        double sigma_t_min = 0.5 * ((sigma_z + sigma_theta) - std::sqrt(std::pow(sigma_z - sigma_theta, 2) + tauSqrx4)) - m_porePressure;
+        if (sigma_t_min < sigma_t_min_min)
+        {
+            sigma_t_min_min = sigma_t_min;
+            *thetaAtMin = stressComponentsForAngle[0];
+        }
+    }
+    return sigma_t_min_min;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RigGeoMechBoreHoleStressCalculator::stassiDalia(double wellPressure, double* thetaAtMin) const
+{
+    CVF_ASSERT(thetaAtMin);
+    double minStassiDalia = std::numeric_limits<double>::max();
+    for (const cvf::Vec4d& stressComponentsForAngle : m_stressComponents)
+    {
+        double sigma_theta = stressComponentsForAngle[1] - wellPressure;
+        const double& sigma_z = stressComponentsForAngle[2];
+        double tauSqrx4 = std::pow(stressComponentsForAngle[3], 2) * 4.0;
+
+        double sigma_1 = wellPressure - m_porePressure;
+        double sigma_2 = 0.5 * ((sigma_z + sigma_theta) + std::sqrt(std::pow(sigma_z - sigma_theta, 2) + tauSqrx4)) - m_porePressure;
+        double sigma_3 = 0.5 * ((sigma_z + sigma_theta) - std::sqrt(std::pow(sigma_z - sigma_theta, 2) + tauSqrx4)) - m_porePressure;
+
+        double stassiDalia = std::pow(sigma_1 - sigma_2, 2) + std::pow(sigma_2 - sigma_3, 2) + std::pow(sigma_1 - sigma_3, 2)
+                          - 2 * m_uniaxialCompressiveStrength * (sigma_1 + sigma_2 + sigma_3);
+
+        if (stassiDalia < minStassiDalia)
+        {
+            minStassiDalia = stassiDalia;
+            *thetaAtMin = stressComponentsForAngle[0];
+        }
+    }
+    return minStassiDalia;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigGeoMechBoreHoleStressCalculator::calculateStressComponents()
+{
+    m_stressComponents.reserve(m_nThetaSubSamples);
+
+    for (int i = 0; i < m_nThetaSubSamples; ++i)
+    {
+        double theta = (i *cvf::PI_F) / (m_nThetaSubSamples - 1.0);
+        cvf::Vec4d stressComponentsForAngle = calculateStressComponentsForSegmentAngle(theta);
+        m_stressComponents.push_back(stressComponentsForAngle);
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::Vec4d RigGeoMechBoreHoleStressCalculator::calculateStressComponentsForSegmentAngle(double theta) const
+{
+    cvf::Vec4d stressComponents;
+
+    const double& sx = m_tensor[caf::Ten3d::SXX];
+    const double& sy = m_tensor[caf::Ten3d::SYY];
+    const double& sz = m_tensor[caf::Ten3d::SZZ];
+    const double& txy = m_tensor[caf::Ten3d::SXY];
+    const double& txz = m_tensor[caf::Ten3d::SZX];
+    const double& tyz = m_tensor[caf::Ten3d::SYZ];
+
+    stressComponents[0] = theta;
+    stressComponents[1] = sx + sy - 2 * (sx - sy) * cos(2 * theta) - 4 * txy * sin(2 * theta);
+    stressComponents[2] = sz - m_poissonRatio * (2 * (sx - sy) * cos(2 * theta) + 4 * txy * sin(2 * theta));
+    stressComponents[3] = 2 * (tyz * cos(theta) - txz * sin(theta));
+
+    return stressComponents;
+}
+
+
+

ApplicationCode/ReservoirDataModel/RigGeoMechBoreHoleStressCalculator.h

@@ -0,0 +1,50 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) Statoil 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 "cafTensor3.h"
+#include "cvfVector3.h"
+#include "cvfVector4.h"
+
+#include <vector>
+
+class RigGeoMechBoreHoleStressCalculator
+{
+public:
+    RigGeoMechBoreHoleStressCalculator(const caf::Ten3d& tensor, double porePressure, double poissonRatio, double uniaxialCompressiveStrength, int nThetaSubSamples);
+    double solveFractureGradient(double* thetaOut = nullptr);
+    double solveStassiDalia(double* thetaOut = nullptr);
+    cvf::Vec3d principleStressesAtWall(double pw, double theta) const;
+private:
+    typedef double (RigGeoMechBoreHoleStressCalculator::*MemberFunc)(double pw, double* thetaOut) const;
+    double solveBisection(double minPw, double maxPw, MemberFunc fn, double* thetaOut);
+    double solveSecant(MemberFunc fn, double* thetaOut);
+    double sigmaTMinOfMin(double wellPressure, double* thetaAtMin) const;
+    double stassiDalia(double wellPressure, double* thetaAtMin) const;
+    void calculateStressComponents();
+    cvf::Vec4d calculateStressComponentsForSegmentAngle(double theta) const;
+
+    caf::Ten3d m_tensor;
+    double m_porePressure;
+    double m_poissonRatio;
+    double m_uniaxialCompressiveStrength;
+    int m_nThetaSubSamples;
+    std::vector<cvf::Vec4d> m_stressComponents;
+};
+

ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp

@@ -21,6 +21,8 @@
 /// 
 //==================================================================================================
 #include "RigGeoMechWellLogExtractor.h"
+
+#include "RigGeoMechBoreHoleStressCalculator.h"
 #include "RigFemPart.h"
 #include "RigFemPartCollection.h"
 #include "RigGeoMechCaseData.h"
@@ -28,9 +30,13 @@
 
 #include "RigWellLogExtractionTools.h"
 #include "RigWellPath.h"
-#include "cvfGeometryTools.h"
 #include "RigWellPathIntersectionTools.h"
 
+#include "cafTensor3.h"
+#include "cvfGeometryTools.h"
+
+
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -50,7 +56,13 @@ void RigGeoMechWellLogExtractor::curveData(const RigFemResultAddress& resAddr, i
 {   
     CVF_TIGHT_ASSERT(values);
     
-    if (!resAddr.isValid()) return ;
+    if (resAddr.fieldName == "FractureGradient" || resAddr.fieldName == "ShearFailureGradient")
+    {
+        wellPathDerivedCurveData(resAddr, frameIndex, values);
+        return;
+    }
+
+    if (!resAddr.isValid()) return;
 
     RigFemResultAddress convResAddr = resAddr;
 
@@ -59,69 +71,225 @@ void RigGeoMechWellLogExtractor::curveData(const RigFemResultAddress& resAddr, i
      
     if (convResAddr.fieldName == "POR-Bar") convResAddr.resultPosType = RIG_ELEMENT_NODAL;
 
+    CVF_ASSERT(resAddr.resultPosType != RIG_WELLPATH_DERIVED);
+
     const RigFemPart* femPart                 = m_caseData->femParts()->part(0);
-    const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
     const std::vector<float>& resultValues    = m_caseData->femPartResults()->resultValues(convResAddr, 0, frameIndex);
 
     if (!resultValues.size()) return;
 
     values->resize(m_intersections.size());
 
-    for (size_t cpIdx = 0; cpIdx < m_intersections.size(); ++cpIdx)
+    for (size_t intersectionIdx = 0; intersectionIdx < m_intersections.size(); ++intersectionIdx)
     {
-        size_t elmIdx = m_intersectedCellsGlobIdx[cpIdx];
+        (*values)[intersectionIdx] = static_cast<double>(interpolateGridResultValue<float>(convResAddr.resultPosType, resultValues, intersectionIdx, false));
+    }
+  
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigGeoMechWellLogExtractor::wellPathDerivedCurveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values)
+{
+    // TODO: Read in these values:
+    const double poissonRatio = 0.25; // TODO: Read this in.
+    // Typical UCS: http://ceae.colorado.edu/~amadei/CVEN5768/PDF/NOTES8.pdf
+    // Typical UCS for Shale is 5 - 100 MPa -> 50 - 1000 bar.
+    const double uniaxialStrengthInBars = 100.0;
+
+    CVF_TIGHT_ASSERT(values);
+    CVF_ASSERT(resAddr.fieldName == "FractureGradient" || resAddr.fieldName == "ShearFailureGradient");
+
+    const RigFemPart* femPart = m_caseData->femParts()->part(0);
+    const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
+    RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults();
+
+    RigFemResultAddress stressResAddr(RIG_ELEMENT_NODAL, std::string("ST"), "");
+    stressResAddr.fieldName = std::string("ST");
+
+    RigFemResultAddress porBarResAddr(RIG_ELEMENT_NODAL, std::string("POR-Bar"), "");
+
+    std::vector<caf::Ten3f> vertexStressesFloat = resultCollection->tensors(stressResAddr, 0, frameIndex);
+    
+    if (!vertexStressesFloat.size()) return;
+
+    std::vector<caf::Ten3d> vertexStresses; vertexStresses.reserve(vertexStressesFloat.size());
+    for (const caf::Ten3f& floatTensor : vertexStressesFloat)
+    {
+        vertexStresses.push_back(caf::Ten3d(floatTensor));
+    }
+
+    values->resize(m_intersections.size(), 0.0f);
+
+    std::vector<float> porePressures = resultCollection->resultValues(porBarResAddr, 0, frameIndex);
+
+#pragma omp parallel for
+    for (int64_t intersectionIdx = 0; intersectionIdx < (int64_t) m_intersections.size(); ++intersectionIdx)
+    {        
+        size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
         RigElementType elmType = femPart->elementType(elmIdx);
 
-        if (!(elmType == HEX8  || elmType == HEX8P)) continue;
+        if (!(elmType == HEX8 || elmType == HEX8P)) continue;
 
-        if (convResAddr.resultPosType == RIG_ELEMENT)
+        double trueVerticalDepth = -m_intersections[intersectionIdx].z();
+        double porePressure = trueVerticalDepth * 9.81 / 100.0;
+        if (false && !porePressures.empty())
         {
-            (*values)[cpIdx] = resultValues[elmIdx];
+            float interpolatedPorePressure = interpolateGridResultValue(porBarResAddr.resultPosType, porePressures, intersectionIdx, true);
-            continue;
+            if (interpolatedPorePressure != std::numeric_limits<float>::infinity() &&
+                interpolatedPorePressure != -std::numeric_limits<float>::infinity())
+            {
+                porePressure = static_cast<double>(interpolatedPorePressure);
+            }
         }
 
-        cvf::StructGridInterface::FaceType cellFace = m_intersectedCellFaces[cpIdx];
+        caf::Ten3d interpolatedStress = interpolateGridResultValue(stressResAddr.resultPosType, vertexStresses, intersectionIdx, true);
+        cvf::Vec3d wellPathTangent = calculateWellPathTangent(intersectionIdx);
+        caf::Ten3d wellPathStressFloat = transformTensorToWellPathOrientation(wellPathTangent, interpolatedStress);
+        caf::Ten3d wellPathStressDouble(wellPathStressFloat);
 
-        int faceNodeCount = 0;
+        RigGeoMechBoreHoleStressCalculator sigmaCalculator(wellPathStressDouble, porePressure, poissonRatio, uniaxialStrengthInBars, 32);
-        const int* faceLocalIndices = RigFemTypes::localElmNodeIndicesForFace(elmType, cellFace, &faceNodeCount);
+        double resultValue = 0.0;
-        const int* elmNodeIndices = femPart->connectivities(elmIdx);
+        if (resAddr.fieldName == "FractureGradient")
-
-        cvf::Vec3d v0(nodeCoords[elmNodeIndices[faceLocalIndices[0]]]);
-        cvf::Vec3d v1(nodeCoords[elmNodeIndices[faceLocalIndices[1]]]);
-        cvf::Vec3d v2(nodeCoords[elmNodeIndices[faceLocalIndices[2]]]);
-        cvf::Vec3d v3(nodeCoords[elmNodeIndices[faceLocalIndices[3]]]);
-
-        size_t resIdx0 = cvf::UNDEFINED_SIZE_T;
-        size_t resIdx1 = cvf::UNDEFINED_SIZE_T;
-        size_t resIdx2 = cvf::UNDEFINED_SIZE_T;
-        size_t resIdx3 = cvf::UNDEFINED_SIZE_T;
-
-        if (convResAddr.resultPosType ==  RIG_NODAL)
         {
-            resIdx0 = elmNodeIndices[faceLocalIndices[0]];
+            resultValue = sigmaCalculator.solveFractureGradient();
-            resIdx1 = elmNodeIndices[faceLocalIndices[1]];
-            resIdx2 = elmNodeIndices[faceLocalIndices[2]];
-            resIdx3 = elmNodeIndices[faceLocalIndices[3]];
         }
         else
         {
-            resIdx0 = (size_t)femPart->elementNodeResultIdx((int)elmIdx, faceLocalIndices[0]);
+            CVF_ASSERT(resAddr.fieldName == "ShearFailureGradient");
-            resIdx1 = (size_t)femPart->elementNodeResultIdx((int)elmIdx, faceLocalIndices[1]);
+            resultValue = sigmaCalculator.solveStassiDalia();
-            resIdx2 = (size_t)femPart->elementNodeResultIdx((int)elmIdx, faceLocalIndices[2]);
+        }
-            resIdx3 = (size_t)femPart->elementNodeResultIdx((int)elmIdx, faceLocalIndices[3]);
+        double effectiveDepth = -m_intersections[intersectionIdx].z() + m_rkbDiff;
+        if (effectiveDepth > 1.0e-8)
+        {
+            resultValue *= 100.0 / (effectiveDepth * 9.81);
+        }
+        else
+        {
+            resultValue = std::numeric_limits<double>::infinity();
         }
 
-        double interpolatedValue = cvf::GeometryTools::interpolateQuad( 
+        (*values)[intersectionIdx] = resultValue;
-            v0, resultValues[resIdx0],
-            v1, resultValues[resIdx1],
-            v2, resultValues[resIdx2],
-            v3, resultValues[resIdx3],
-            m_intersections[cpIdx]
-        );  
-
-        (*values)[cpIdx] = interpolatedValue;
     }
-  
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+const RigGeoMechCaseData* RigGeoMechWellLogExtractor::caseData()
+{
+    return m_caseData.p();
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigGeoMechWellLogExtractor::setRkbDiff(double rkbDiff)
+{
+    m_rkbDiff = rkbDiff;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+template<typename T>
+T RigGeoMechWellLogExtractor::interpolateGridResultValue(RigFemResultPosEnum   resultPosType,
+                                                      const std::vector<T>& gridResultValues,
+                                                      int64_t               intersectionIdx,
+                                                      bool                  averageNodeElementResults) const
+{
+    const RigFemPart* femPart = m_caseData->femParts()->part(0);
+    const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
+
+    size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
+    RigElementType elmType = femPart->elementType(elmIdx);
+
+    if (!(elmType == HEX8 || elmType == HEX8P)) return T();
+
+
+    if (resultPosType == RIG_ELEMENT)
+    {
+        return gridResultValues[elmIdx];        
+    }
+
+    cvf::StructGridInterface::FaceType cellFace = m_intersectedCellFaces[intersectionIdx];
+
+    if (cellFace == cvf::StructGridInterface::NO_FACE)
+    {
+        // TODO: Should interpolate within the whole hexahedron. This requires converting to locals coordinates.
+        // For now just pick the average value for the cell.
+        T sumOfVertexValues = gridResultValues[femPart->elementNodeResultIdx(static_cast<int>(elmIdx), 0)];
+        for (int i = 1; i < 8; ++i)
+        {
+            sumOfVertexValues = sumOfVertexValues + gridResultValues[femPart->elementNodeResultIdx(static_cast<int>(elmIdx), i)];
+        }
+        return sumOfVertexValues * (1.0 / 8.0);
+    }
+
+    int faceNodeCount = 0;
+    const int* faceLocalIndices = RigFemTypes::localElmNodeIndicesForFace(elmType, cellFace, &faceNodeCount);
+    const int* elmNodeIndices = femPart->connectivities(elmIdx);
+
+    cvf::Vec3d v0(nodeCoords[elmNodeIndices[faceLocalIndices[0]]]);
+    cvf::Vec3d v1(nodeCoords[elmNodeIndices[faceLocalIndices[1]]]);
+    cvf::Vec3d v2(nodeCoords[elmNodeIndices[faceLocalIndices[2]]]);
+    cvf::Vec3d v3(nodeCoords[elmNodeIndices[faceLocalIndices[3]]]);
+
+    std::vector<size_t> nodeResIdx(4, cvf::UNDEFINED_SIZE_T);
+
+    if (resultPosType == RIG_NODAL)
+    {
+        for (size_t i = 0; i < nodeResIdx.size(); ++i)
+        {
+            nodeResIdx[i] = elmNodeIndices[faceLocalIndices[i]];
+        }
+    }
+    else
+    {
+        for (size_t i = 0; i < nodeResIdx.size(); ++i)
+        {
+            nodeResIdx[i] = (size_t)femPart->elementNodeResultIdx((int)elmIdx, faceLocalIndices[i]);
+        }
+    }
+
+    std::vector<T> nodeResultValues;
+    nodeResultValues.reserve(4);
+    if (resultPosType == RIG_ELEMENT_NODAL && averageNodeElementResults)
+    {
+        // Estimate nodal values as the average of the node values from each connected element.
+        for (size_t i = 0; i < nodeResIdx.size(); ++i)
+        {
+            int nodeIndex = femPart->nodeIdxFromElementNodeResultIdx(nodeResIdx[i]);
+            const std::vector<int>& elements = femPart->elementsUsingNode(nodeIndex);
+            const std::vector<unsigned char>& localIndices = femPart->elementLocalIndicesForNode(nodeIndex);
+            size_t otherNodeResIdx = femPart->elementNodeResultIdx(elements[0], static_cast<int>(localIndices[0]));
+            T nodeResultValue = gridResultValues[otherNodeResIdx];
+            for (size_t j = 1; j < elements.size(); ++j)
+            {
+                otherNodeResIdx = femPart->elementNodeResultIdx(elements[j], static_cast<int>(localIndices[j]));
+                nodeResultValue = nodeResultValue + gridResultValues[otherNodeResIdx];
+            }
+            nodeResultValue = nodeResultValue * (1.0 / elements.size());
+            nodeResultValues.push_back(nodeResultValue);
+        }
+    }
+    else {
+        for (size_t i = 0; i < nodeResIdx.size(); ++i)
+        {
+            nodeResultValues.push_back(gridResultValues[nodeResIdx[i]]);
+        }
+    }
+
+    T interpolatedValue = cvf::GeometryTools::interpolateQuad<T>(
+        v0, nodeResultValues[0],
+        v1, nodeResultValues[1],
+        v2, nodeResultValues[2],
+        v3, nodeResultValues[3],
+        m_intersections[intersectionIdx]
+    );
+
+    return interpolatedValue;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -224,3 +392,37 @@ cvf::Vec3d RigGeoMechWellLogExtractor::calculateLengthInCell(size_t cellIndex, c
     return RigWellPathIntersectionTools::calculateLengthInCell(hexCorners, startPoint, endPoint); 
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::Vec3d RigGeoMechWellLogExtractor::calculateWellPathTangent(int64_t intersectionIdx) const
+{
+    cvf::Vec3d wellPathTangent;
+    if (intersectionIdx % 2 == 0)
+    {
+        wellPathTangent = m_intersections[intersectionIdx + 1] - m_intersections[intersectionIdx];
+    }
+    else
+    {
+        wellPathTangent = m_intersections[intersectionIdx] - m_intersections[intersectionIdx - 1];
+    }
+    CVF_ASSERT(wellPathTangent.length() > 1.0e-7);
+    return wellPathTangent.getNormalized();
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+caf::Ten3d RigGeoMechWellLogExtractor::transformTensorToWellPathOrientation(const cvf::Vec3d& wellPathTangent,
+                                                                            const caf::Ten3d& tensor)
+{
+    // Create local coordinate system for well path segment
+    cvf::Vec3d local_z = wellPathTangent;
+    cvf::Vec3d local_x = local_z.perpendicularVector().getNormalized();
+    cvf::Vec3d local_y = (local_z ^ local_x).getNormalized();
+    // Calculate the rotation matrix from global i, j, k to local x, y, z.
+    cvf::Mat4d rotationMatrix = cvf::Mat4d::fromCoordSystemAxes(&local_x, &local_y, &local_z);
+
+    return tensor.rotated(rotationMatrix.toMatrix3());
+}
+

ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h

@@ -21,17 +21,21 @@
 
 #include "RigWellLogExtractor.h"
 
+#include "cafTensor3.h"
+
 #include "cvfBase.h"
 #include "cvfObject.h"
 #include "cvfMath.h"
+#include "cvfStructGrid.h"
 #include "cvfVector3.h"
 
 #include <vector>
-#include "cvfStructGrid.h"
 
-class RigWellPath;
+enum RigElementType;
-class RigGeoMechCaseData;
+enum RigFemResultPosEnum;
 class RigFemResultAddress;
+class RigGeoMechCaseData;
+class RigWellPath;
 
 namespace cvf {
     class BoundingBox;
@@ -46,16 +50,25 @@ public:
     RigGeoMechWellLogExtractor(RigGeoMechCaseData* aCase, const RigWellPath* wellpath, const std::string& wellCaseErrorMsgName);
 
     void                         curveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values );
-    const RigGeoMechCaseData*    caseData()     { return m_caseData.p();}
+    void                         wellPathDerivedCurveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values);
-
+    const RigGeoMechCaseData*    caseData();
+    void                         setRkbDiff(double rkbDiff);
 private:
+  
+   
+    template<typename T>
+    T                            interpolateGridResultValue(RigFemResultPosEnum resultPosType, const std::vector<T>& gridResultValues, int64_t intersectionIdx, bool averageNodeElementResults) const;
     void                         calculateIntersection();
     std::vector<size_t>          findCloseCells(const cvf::BoundingBox& bb);
     virtual cvf::Vec3d           calculateLengthInCell(size_t cellIndex, 
                                                        const cvf::Vec3d& startPoint, 
                                                        const cvf::Vec3d& endPoint) const override;
+    cvf::Vec3d                   calculateWellPathTangent(int64_t intersectionIdx) const;
+    static caf::Ten3d            transformTensorToWellPathOrientation(const cvf::Vec3d& wellPathTangent,
+                                                                      const caf::Ten3d& wellPathTensor);
 
     cvf::ref<RigGeoMechCaseData> m_caseData;
+    double                       m_rkbDiff;
 };
 
 

ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -589,20 +589,6 @@ cvf::Vec4d GeometryTools::barycentricCoords(const cvf::Vec3d& v0,
     return w;
 }
 
-//--------------------------------------------------------------------------------------------------
-/// 
-//--------------------------------------------------------------------------------------------------
-double GeometryTools::interpolateQuad(const cvf::Vec3d& v1, double s1, 
-                                      const cvf::Vec3d& v2, double s2, 
-                                      const cvf::Vec3d& v3, double s3, 
-                                      const cvf::Vec3d& v4, double s4, 
-                                      const cvf::Vec3d& point)
-{
-    cvf::Vec4d bc = barycentricCoords(v1, v2, v3, v4, point);
-
-    return s1*bc[0] + s2*bc[1] + s3*bc[2] + s4*bc[3];
-}
-
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/cvfGeometryTools.h

@@ -47,11 +47,12 @@ public:
                                                     bool * isLineDirDotNormalNegative);
     static cvf::Vec3d barycentricCoords(const cvf::Vec3d&  t0, const cvf::Vec3d&  t1, const cvf::Vec3d&  t2, const cvf::Vec3d&  p);
     static cvf::Vec4d barycentricCoords(const cvf::Vec3d&  v0, const cvf::Vec3d&  v1, const cvf::Vec3d&  v2, const cvf::Vec3d& v3, const cvf::Vec3d&  p);
-    static double interpolateQuad(const cvf::Vec3d& v1, double s1,
+    template<typename DataType>
-                                  const cvf::Vec3d& v2, double s2,
+    static DataType interpolateQuad(const cvf::Vec3d& v1, DataType s1,
-                                  const cvf::Vec3d& v3, double s3,
+                                    const cvf::Vec3d& v2, DataType s2,
-                                  const cvf::Vec3d& v4, double s4,
+                                    const cvf::Vec3d& v3, DataType s3,
-                                  const cvf::Vec3d& point);
+                                    const cvf::Vec3d& v4, DataType s4,
+                                    const cvf::Vec3d& point);
 
     static int        findClosestAxis(const cvf::Vec3d& vec );
     static double     getAngle(const cvf::Vec3d& positiveNormalAxis, const cvf::Vec3d& v1, const cvf::Vec3d& v2);

ApplicationCode/ReservoirDataModel/cvfGeometryTools.inl

@@ -24,6 +24,22 @@
 namespace cvf
 {
 
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+template<typename DataType>
+DataType GeometryTools::interpolateQuad(const cvf::Vec3d& v1, DataType s1,
+                                        const cvf::Vec3d& v2, DataType s2,
+                                        const cvf::Vec3d& v3, DataType s3,
+                                        const cvf::Vec3d& v4, DataType s4,
+                                        const cvf::Vec3d& point)
+{
+    cvf::Vec4d bc = barycentricCoords(v1, v2, v3, v4, point);
+
+    return s1 * bc[0] + s2 * bc[1] + s3 * bc[2] + s4 * bc[3];
+}
+
 //--------------------------------------------------------------------------------------------------
 /// Inserts the vertex into the polygon if it fits along one of the edges within the tolerance.
 /// The method returns true if it was inserted, or if it was already in the polygon, or if it was 

Fwk/AppFwk/cafTensor/cafTensor3.h

@@ -32,6 +32,8 @@ public:
     Tensor3() {}
     Tensor3(S sxx, S syy, S szz, S sxy, S syz, S szx); 
     Tensor3(const Tensor3& other);
+    template<typename T>
+    explicit Tensor3(const Tensor3<T>& other);
 
     inline Tensor3&     operator=(const Tensor3& rhs);
     inline Tensor3      operator+(const Tensor3& rhs) const;
@@ -55,6 +57,7 @@ public:
 };
 
 typedef Tensor3<float> Ten3f;
+typedef Tensor3<double> Ten3d;
 
 }
 

Fwk/AppFwk/cafTensor/cafTensor3.inl

@@ -37,6 +37,21 @@ inline Tensor3<S>::Tensor3(const Tensor3& other)
     cvf::System::memcpy(m_tensor, sizeof(m_tensor), other.m_tensor, sizeof(other.m_tensor));
 }
 
+//----------------------------------------------------------------------------------------------------
+/// Explicit Cast constructor
+//----------------------------------------------------------------------------------------------------
+template <typename S>
+template <typename T>
+Tensor3<S>::Tensor3(const Tensor3<T>& other)
+{
+    m_tensor[SXX] = other[Tensor3<T>::SXX];
+    m_tensor[SYY] = other[Tensor3<T>::SYY];
+    m_tensor[SZZ] = other[Tensor3<T>::SZZ];
+    m_tensor[SXY] = other[Tensor3<T>::SXY];
+    m_tensor[SYZ] = other[Tensor3<T>::SYZ];
+    m_tensor[SZX] = other[Tensor3<T>::SZX];    
+}
+
 //----------------------------------------------------------------------------------------------------
 /// Constructor with explicit initialization of all tensor elements.
 ///