ResInsight/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp

/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) Statoil ASA
//  Copyright (C) Ceetron Solutions AS
// 
//  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 "RigGeoMechWellLogExtractor.h"

#include "RiaDefines.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigFemTypes.h"
#include "RigGeoMechBoreHoleStressCalculator.h"
#include "RigFemPart.h"
#include "RigFemPartCollection.h"
#include "RigGeoMechCaseData.h"
#include "RigFemPartResultsCollection.h"

#include "RigWellLogExtractionTools.h"
#include "RigWellPath.h"
#include "RigWellPathIntersectionTools.h"

#include "cafTensor3.h"
#include "cvfGeometryTools.h"
#include "cvfMath.h"

#include <type_traits>

const double RigGeoMechWellLogExtractor::UNIT_WEIGHT_OF_WATER = 9.81 * 1000.0; // N / m^3

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigGeoMechWellLogExtractor::RigGeoMechWellLogExtractor(RigGeoMechCaseData* aCase,
                                                       const RigWellPath*  wellpath,
                                                       const std::string&  wellCaseErrorMsgName)
    : RigWellLogExtractor(wellpath, wellCaseErrorMsgName)
    , m_caseData(aCase)
{
    calculateIntersection();
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::curveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values)
{   
    CVF_TIGHT_ASSERT(values);
    
    if (resAddr.resultPosType == RIG_WELLPATH_DERIVED)
    {
        if (resAddr.fieldName == RiaDefines::wellPathFGResultName().toStdString() || resAddr.fieldName == RiaDefines::wellPathSFGResultName().toStdString())
        {
            wellBoreWallCurveData(resAddr, frameIndex, values);
            return;
        }
        else if (resAddr.fieldName == "PP" || resAddr.fieldName == "OBG" || resAddr.fieldName == "SH")
        {
            wellPathScaledCurveData(resAddr, frameIndex, values);
            return;
        }
        else if (resAddr.fieldName == "Azimuth" || resAddr.fieldName == "Inclination")
        {
            wellPathAngles(resAddr, values);
            return;
        }

    }

    if (!resAddr.isValid()) return;

    RigFemResultAddress convResAddr = resAddr;

    // When showing POR results, always use the element nodal result, 
    // to get correct handling of elements without POR results
     
    if (convResAddr.fieldName == "POR-Bar") convResAddr.resultPosType = RIG_ELEMENT_NODAL;

    CVF_ASSERT(resAddr.resultPosType != RIG_WELLPATH_DERIVED);

    const std::vector<float>& resultValues    = m_caseData->femPartResults()->resultValues(convResAddr, 0, frameIndex);

    if (!resultValues.size()) return;

    values->resize(m_intersections.size());

    for (size_t intersectionIdx = 0; intersectionIdx < m_intersections.size(); ++intersectionIdx)
    {
        (*values)[intersectionIdx] = static_cast<double>(interpolateGridResultValue<float>(convResAddr.resultPosType, resultValues, intersectionIdx, false));
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
float RigGeoMechWellLogExtractor::calculatePorePressureInSegment(int64_t intersectionIdx, float averageSegmentPorePressureBar, double hydroStaticPorePressureBar, double effectiveDepthMeters, const std::vector<float>& poreElementPressuresPascal) const
{
    double porePressure = hydroStaticPorePressureBar;

    // 1: Try pore pressure from the grid
    if (porePressure == hydroStaticPorePressureBar && averageSegmentPorePressureBar != std::numeric_limits<double>::infinity() && averageSegmentPorePressureBar > 0.0)
    {
        porePressure = averageSegmentPorePressureBar;
    }

    // 2: Try mud weight from LAS-file to generate pore pressure
    if (porePressure == hydroStaticPorePressureBar && !m_wellLogMdAndMudWeightKgPerM3.empty())
    {
        double lasMudWeightKgPerM3 = getWellLogSegmentValue(intersectionIdx, m_wellLogMdAndMudWeightKgPerM3);
        if (lasMudWeightKgPerM3 != std::numeric_limits<double>::infinity())
        {
            double specificMudWeightNPerM3 = lasMudWeightKgPerM3 * 9.81;
            double porePressurePascal = specificMudWeightNPerM3 * effectiveDepthMeters;
            porePressure = pascalToBar(porePressurePascal);
        }
    }
    size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
    // 3: Try pore pressure from element property tables
    if (porePressure == hydroStaticPorePressureBar && elmIdx < poreElementPressuresPascal.size())
    {
        // Pore pressure from element property tables are in pascal.
        porePressure = pascalToBar(poreElementPressuresPascal[elmIdx]);
    }
    // 4: If no pore-pressure was found, the default value of hydrostatic pore pressure is used.

    CVF_ASSERT(porePressure >= 0.0);
    return porePressure;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
float RigGeoMechWellLogExtractor::calculatePoissonRatio(int64_t intersectionIdx, const std::vector<float>& poissonRatios) const
{
    const double defaultPoissonRatio = 0.25;

    double poissonRatio = defaultPoissonRatio;

    if (!m_wellLogMdAndPoissonRatios.empty())
    {
        double lasPoissionRatio = getWellLogSegmentValue(intersectionIdx, m_wellLogMdAndPoissonRatios);
        if (lasPoissionRatio != std::numeric_limits<double>::infinity())
        {
            poissonRatio = lasPoissionRatio;
        }
    }

    size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
    if (poissonRatio == defaultPoissonRatio && elmIdx < poissonRatios.size())
    {
        poissonRatio = poissonRatios[elmIdx];
    }
    return poissonRatio;
}


//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
float RigGeoMechWellLogExtractor::calculateUcs(int64_t intersectionIdx, const std::vector<float>& ucsValuesPascal) const
{
    // Typical UCS: http://ceae.colorado.edu/~amadei/CVEN5768/PDF/NOTES8.pdf
    // Typical UCS for Shale is 5 - 100 MPa -> 50 - 1000 bar.
    const double defaultUniaxialStrengthInBar = 100.0;

    double uniaxialStrengthInBar = defaultUniaxialStrengthInBar;
    if (!m_wellLogMdAndUcsBar.empty())
    {
        double lasUniaxialStrengthInBar = getWellLogSegmentValue(intersectionIdx, m_wellLogMdAndUcsBar);
        if (lasUniaxialStrengthInBar != std::numeric_limits<double>::infinity())
        {
            uniaxialStrengthInBar = lasUniaxialStrengthInBar;
        }
    }

    size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
    if (uniaxialStrengthInBar == defaultUniaxialStrengthInBar && elmIdx < ucsValuesPascal.size())
    {
        // Read UCS from element table in Pascal
        uniaxialStrengthInBar = pascalToBar(ucsValuesPascal[elmIdx]);
    }
    return uniaxialStrengthInBar;
}


//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::wellPathAngles(const RigFemResultAddress& resAddr, std::vector<double>* values)
{
    CVF_ASSERT(values);
    CVF_ASSERT(resAddr.fieldName == "Azimuth" || resAddr.fieldName == "Inclination");
    values->resize(m_intersections.size(), 0.0f);
    const double epsilon = 1.0e-6 * 360;
    const cvf::Vec3d trueNorth(0.0, 1.0, 0.0);
    const cvf::Vec3d up(0.0, 0.0, 1.0);
    for (int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx)
    {
        cvf::Vec3d wellPathTangent = calculateWellPathTangent(intersectionIdx, TangentFollowWellPathSegments);
   
        // Deviation from vertical. Since well path is tending downwards we compare with negative z.
        double inclination = cvf::Math::toDegrees(std::acos(cvf::Vec3d(0.0, 0.0, -1.0) * wellPathTangent.getNormalized()));       

        if (resAddr.fieldName == "Azimuth")
        {
            double azimuth = HUGE_VAL;
            
            // Azimuth is not defined when well path is vertical. We define it as infinite to avoid it showing up in the plot.
            if (cvf::Math::valueInRange(inclination, epsilon, 180.0 - epsilon))
            {
                cvf::Vec3d projectedTangentXY = wellPathTangent;
                projectedTangentXY.z() = 0.0;

                // Do tangentXY to true north for clockwise angles.
                double dotProduct = projectedTangentXY * trueNorth;
                double crossProduct = (projectedTangentXY ^ trueNorth) * up;
                // http://www.glossary.oilfield.slb.com/Terms/a/azimuth.aspx
                azimuth = cvf::Math::toDegrees(std::atan2(crossProduct, dotProduct));
                if (azimuth < 0.0)
                {
                    // Straight atan2 gives angle from -PI to PI yielding angles from -180 to 180
                    // where the negative angles are counter clockwise.
                    // To get all positive clockwise angles, we add 360 degrees to negative angles.
                    azimuth = azimuth + 360.0;
                }
            }

            (*values)[intersectionIdx] = azimuth;
        }
        else
        {
            (*values)[intersectionIdx] = inclination;
        }
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::wellPathScaledCurveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values)
{
    CVF_ASSERT(values);

    const RigFemPart* femPart = m_caseData->femParts()->part(0);
    RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults();

    std::string nativeFieldName;
    std::string nativeCompName;
    if (resAddr.fieldName == "PP")
    {
        nativeFieldName = "POR-Bar"; // More likely to be in memory than POR
    }
    else if (resAddr.fieldName == "OBG")
    {
        nativeFieldName = "ST";
        nativeCompName = "S33";
    }
    else if (resAddr.fieldName == "SH")
    {
        nativeFieldName = "ST";
        nativeCompName = "S3";
    }

    RigFemResultAddress nativeAddr(RIG_ELEMENT_NODAL, nativeFieldName, nativeCompName);
    RigFemResultAddress porElementResAddr(RIG_ELEMENT, "POR", "");

    std::vector<float> unscaledResultValues = resultCollection->resultValues(nativeAddr, 0, frameIndex);
    std::vector<float> poreElementPressuresPascal = resultCollection->resultValues(porElementResAddr, 0, frameIndex);

    std::vector<float> interpolatedInterfaceValues;
    interpolatedInterfaceValues.resize(m_intersections.size(), std::numeric_limits<double>::infinity());

#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;

        interpolatedInterfaceValues[intersectionIdx] = interpolateGridResultValue<float>(nativeAddr.resultPosType, unscaledResultValues, intersectionIdx, false);
    }

    values->resize(m_intersections.size(), 0.0f);

#pragma omp parallel for
    for (int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx)
    {
        // Set the value to invalid by default
        (*values)[intersectionIdx] = std::numeric_limits<double>::infinity();

        size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
        RigElementType elmType = femPart->elementType(elmIdx);

        if (!(elmType == HEX8 || elmType == HEX8P)) continue;

        cvf::Vec3f centroid = cellCentroid(intersectionIdx);

        double trueVerticalDepth = -centroid.z();
        
        double effectiveDepthMeters = trueVerticalDepth + m_rkbDiff;
        double hydroStaticPorePressureBar = pascalToBar(effectiveDepthMeters * UNIT_WEIGHT_OF_WATER);

        float averageUnscaledValue = std::numeric_limits<float>::infinity();
        bool validAverage = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfaceValues, std::numeric_limits<float>::infinity(), &averageUnscaledValue);

        if (resAddr.fieldName == "PP" && validAverage)
        {
            double segmentPorePressureFromGrid = averageUnscaledValue;
            averageUnscaledValue = calculatePorePressureInSegment(intersectionIdx, segmentPorePressureFromGrid, hydroStaticPorePressureBar, effectiveDepthMeters, poreElementPressuresPascal);

        }

        (*values)[intersectionIdx] = static_cast<double>(averageUnscaledValue) / hydroStaticPorePressureBar;            
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::wellBoreWallCurveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values)
{
    CVF_ASSERT(values);
    CVF_ASSERT(resAddr.fieldName == RiaDefines::wellPathFGResultName().toStdString() || resAddr.fieldName == RiaDefines::wellPathSFGResultName().toStdString());

    // The result addresses needed
    RigFemResultAddress stressResAddr(RIG_ELEMENT_NODAL, "ST", "");
    RigFemResultAddress porBarResAddr(RIG_ELEMENT_NODAL, "POR-Bar", "");
    // Allow POR as an element property value
    RigFemResultAddress porElementResAddr(RIG_ELEMENT, "POR", "");
    RigFemResultAddress poissonResAddr(RIG_ELEMENT, "RATIO", "");
    RigFemResultAddress ucsResAddr(RIG_ELEMENT, "UCS", "");
    
    const RigFemPart* femPart = m_caseData->femParts()->part(0);
    RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults();

    // Load results
    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));
    }
    std::vector<float> porePressures              = resultCollection->resultValues(porBarResAddr, 0, frameIndex);
    std::vector<float> poreElementPressuresPascal = resultCollection->resultValues(porElementResAddr, 0, frameIndex);
    std::vector<float> poissonRatios              = resultCollection->resultValues(poissonResAddr, 0, frameIndex);
    std::vector<float> ucsValuesPascal            = resultCollection->resultValues(ucsResAddr, 0, frameIndex);

    std::vector<float> interpolatedInterfacePorePressureBar;
    interpolatedInterfacePorePressureBar.resize(m_intersections.size(), std::numeric_limits<double>::infinity());

    std::vector<caf::Ten3d> interpolatedInterfaceStressBar;
    interpolatedInterfaceStressBar.resize(m_intersections.size());
#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;

        interpolatedInterfacePorePressureBar[intersectionIdx]     = interpolateGridResultValue(porBarResAddr.resultPosType, porePressures, intersectionIdx, false);
        interpolatedInterfaceStressBar[intersectionIdx] = interpolateGridResultValue(stressResAddr.resultPosType, vertexStresses, intersectionIdx, false);
    }

    values->resize(m_intersections.size(), 0.0f);
    
#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;

        cvf::Vec3f centroid = cellCentroid(intersectionIdx);

        double trueVerticalDepth = -centroid.z();
        double effectiveDepthMeters = trueVerticalDepth + m_rkbDiff;
        double hydroStaticPorePressureBar = pascalToBar(effectiveDepthMeters * UNIT_WEIGHT_OF_WATER);

        float averagePorePressureBar = std::numeric_limits<float>::infinity();
        bool validGridPorePressure = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfacePorePressureBar, std::numeric_limits<float>::infinity(), &averagePorePressureBar);
        bool isFGregion = validGridPorePressure; // FG is for sands, SFG for shale. Sands has PP, shale does not.

        double porePressureBar = calculatePorePressureInSegment(intersectionIdx, averagePorePressureBar, hydroStaticPorePressureBar, effectiveDepthMeters, poreElementPressuresPascal);
        double poissonRatio    = calculatePoissonRatio(intersectionIdx, poissonRatios);
        double ucsBar          = calculateUcs(intersectionIdx, ucsValuesPascal);

        caf::Ten3d segmentStress;
        bool validSegmentStress = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfaceStressBar, caf::Ten3d::invalid(), &segmentStress);

        cvf::Vec3d wellPathTangent = calculateWellPathTangent(intersectionIdx, TangentConstantWithinCell);
        caf::Ten3d wellPathStressFloat = transformTensorToWellPathOrientation(wellPathTangent, segmentStress);
        caf::Ten3d wellPathStressDouble(wellPathStressFloat);

        RigGeoMechBoreHoleStressCalculator sigmaCalculator(wellPathStressDouble, porePressureBar, poissonRatio, ucsBar, 32);
        double resultValue = std::numeric_limits<double>::infinity();
        if (resAddr.fieldName == RiaDefines::wellPathFGResultName().toStdString())
        {
            if (isFGregion && validSegmentStress)
            {
                resultValue = sigmaCalculator.solveFractureGradient();
            }
        }
        else
        {
            CVF_ASSERT(resAddr.fieldName == RiaDefines::wellPathSFGResultName().toStdString());
            if (!isFGregion && validSegmentStress)
            {
                resultValue = sigmaCalculator.solveStassiDalia();
            }
        }
        if (resultValue != std::numeric_limits<double>::infinity())
        {
            if (hydroStaticPorePressureBar > 1.0e-8)
            {
                resultValue /= hydroStaticPorePressureBar;
            }
        }
        (*values)[intersectionIdx] = resultValue;
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const RigGeoMechCaseData* RigGeoMechWellLogExtractor::caseData()
{
    return m_caseData.p();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::setRkbDiff(double rkbDiff)
{
    m_rkbDiff = rkbDiff;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::setWellLogMdAndMudWeightKgPerM3(const std::vector<std::pair<double, double>>& porePressures)
{
    m_wellLogMdAndMudWeightKgPerM3 = porePressures;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::setWellLogMdAndUcsBar(const std::vector<std::pair<double, double>>& ucsValues)
{
    m_wellLogMdAndUcsBar = ucsValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::setWellLogMdAndPoissonRatio(const std::vector<std::pair<double, double>>& poissonRatios)
{
    m_wellLogMdAndPoissonRatios = poissonRatios;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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_FORMATION_NAMES)
    {
        resultPosType = RIG_ELEMENT_NODAL; // formation indices are stored per element node result.
    }

    if (resultPosType == RIG_ELEMENT)
    {
        return gridResultValues[elmIdx];        
    }

    cvf::StructGridInterface::FaceType cellFace = m_intersectedCellFaces[intersectionIdx];

    if (cellFace == cvf::StructGridInterface::NO_FACE)
    {
        if (resultPosType == RIG_ELEMENT_NODAL_FACE)
        {
            return std::numeric_limits<T>::infinity(); // undefined value. ELEMENT_NODAL_FACE values are only defined on a face.
        }
        // TODO: Should interpolate within the whole hexahedron. This requires converting to locals coordinates.
        // For now just pick the average value for the cell.
        size_t gridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(elmIdx), 0);
        T sumOfVertexValues = gridResultValues[gridResultValueIdx];
        for (int i = 1; i < 8; ++i)
        {
            gridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(elmIdx), i);
            sumOfVertexValues = sumOfVertexValues + gridResultValues[gridResultValueIdx];
        }
        return sumOfVertexValues * (1.0 / 8.0);
    }

    int faceNodeCount = 0;
    const int* elementLocalIndicesForFace = RigFemTypes::localElmNodeIndicesForFace(elmType, cellFace, &faceNodeCount);
    const int* elmNodeIndices = femPart->connectivities(elmIdx);

    cvf::Vec3d v0(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[0]]]);
    cvf::Vec3d v1(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[1]]]);
    cvf::Vec3d v2(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[2]]]);
    cvf::Vec3d v3(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[3]]]);

    std::vector<size_t> nodeResIdx(4, cvf::UNDEFINED_SIZE_T);

    for (size_t i = 0; i < nodeResIdx.size(); ++i)
    {
        if (resultPosType == RIG_ELEMENT_NODAL_FACE)
        {
            nodeResIdx[i] = gridResultIndexFace(elmIdx, cellFace, static_cast<int>(i));
        }
        else
        {
            nodeResIdx[i] = femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(elmIdx), elementLocalIndicesForFace[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 otherGridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, elements[0], static_cast<int>(localIndices[0]));
            T nodeResultValue = gridResultValues[otherGridResultValueIdx];
            for (size_t j = 1; j < elements.size(); ++j)
            {
                otherGridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, elements[j], static_cast<int>(localIndices[j]));
                nodeResultValue = nodeResultValue + gridResultValues[otherGridResultValueIdx];
            }
            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;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigGeoMechWellLogExtractor::gridResultIndexFace(size_t elementIdx, cvf::StructGridInterface::FaceType cellFace, int faceLocalNodeIdx) const
{
    CVF_ASSERT(cellFace != cvf::StructGridInterface::NO_FACE && faceLocalNodeIdx < 4);
    return elementIdx * 24 + static_cast<int>(cellFace) * 4 + faceLocalNodeIdx;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::calculateIntersection()
{
    CVF_ASSERT(m_caseData->femParts()->partCount() == 1);

    std::map<RigMDCellIdxEnterLeaveKey, HexIntersectionInfo > uniqueIntersections;

    const RigFemPart* femPart = m_caseData->femParts()->part(0);
    const std::vector<cvf::Vec3f>& nodeCoords =  femPart->nodes().coordinates;

    for (size_t wpp = 0; wpp < m_wellPath->m_wellPathPoints.size() - 1; ++wpp)
    {
        std::vector<HexIntersectionInfo> intersections;
        cvf::Vec3d p1 = m_wellPath->m_wellPathPoints[wpp];
        cvf::Vec3d p2 = m_wellPath->m_wellPathPoints[wpp+1];

        cvf::BoundingBox bb;

        bb.add(p1);
        bb.add(p2);

        std::vector<size_t> closeCells = findCloseCells(bb);

        cvf::Vec3d hexCorners[8];
        for (size_t ccIdx = 0; ccIdx < closeCells.size(); ++ccIdx)
        {
            RigElementType elmType = femPart->elementType(closeCells[ccIdx]);
            if (!(elmType == HEX8 || elmType == HEX8P)) continue;

            const int* cornerIndices = femPart->connectivities(closeCells[ccIdx]);

            hexCorners[0] = cvf::Vec3d(nodeCoords[cornerIndices[0]]);
            hexCorners[1] = cvf::Vec3d(nodeCoords[cornerIndices[1]]);
            hexCorners[2] = cvf::Vec3d(nodeCoords[cornerIndices[2]]);
            hexCorners[3] = cvf::Vec3d(nodeCoords[cornerIndices[3]]);
            hexCorners[4] = cvf::Vec3d(nodeCoords[cornerIndices[4]]);
            hexCorners[5] = cvf::Vec3d(nodeCoords[cornerIndices[5]]);
            hexCorners[6] = cvf::Vec3d(nodeCoords[cornerIndices[6]]);
            hexCorners[7] = cvf::Vec3d(nodeCoords[cornerIndices[7]]);

            //int intersectionCount = RigHexIntersector::lineHexCellIntersection(p1, p2, hexCorners, closeCells[ccIdx], &intersections);
            RigHexIntersectionTools::lineHexCellIntersection(p1, p2, hexCorners, closeCells[ccIdx], &intersections);
        }

        // Now, with all the intersections of this piece of line, we need to 
        // sort them in order, and set the measured depth and corresponding cell index

        // Inserting the intersections in this map will remove identical intersections
        // and sort them according to MD, CellIdx, Leave/enter

        double md1 = m_wellPath->m_measuredDepths[wpp];
        double md2 = m_wellPath->m_measuredDepths[wpp+1];

        insertIntersectionsInMap(intersections,
                                 p1, md1, p2, md2,
                                 &uniqueIntersections);
    }

    this->populateReturnArrays(uniqueIntersections);
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
std::vector<size_t> RigGeoMechWellLogExtractor::findCloseCells(const cvf::BoundingBox& bb)
{
    std::vector<size_t> closeCells;

    if (m_caseData->femParts()->partCount())
    {
        m_caseData->femParts()->part(0)->findIntersectingCells(bb, &closeCells);
    }
    return closeCells;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigGeoMechWellLogExtractor::calculateLengthInCell(size_t cellIndex, const cvf::Vec3d& startPoint, const cvf::Vec3d& endPoint) const
{
    std::array<cvf::Vec3d, 8> hexCorners;

    const RigFemPart* femPart = m_caseData->femParts()->part(0);
    const std::vector<cvf::Vec3f>& nodeCoords =  femPart->nodes().coordinates;
    const int* cornerIndices = femPart->connectivities(cellIndex);

    hexCorners[0] = cvf::Vec3d(nodeCoords[cornerIndices[0]]);
    hexCorners[1] = cvf::Vec3d(nodeCoords[cornerIndices[1]]);
    hexCorners[2] = cvf::Vec3d(nodeCoords[cornerIndices[2]]);
    hexCorners[3] = cvf::Vec3d(nodeCoords[cornerIndices[3]]);
    hexCorners[4] = cvf::Vec3d(nodeCoords[cornerIndices[4]]);
    hexCorners[5] = cvf::Vec3d(nodeCoords[cornerIndices[5]]);
    hexCorners[6] = cvf::Vec3d(nodeCoords[cornerIndices[6]]);
    hexCorners[7] = cvf::Vec3d(nodeCoords[cornerIndices[7]]);

    return RigWellPathIntersectionTools::calculateLengthInCell(hexCorners, startPoint, endPoint); 
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigGeoMechWellLogExtractor::calculateWellPathTangent(int64_t                    intersectionIdx,
                                                                WellPathTangentCalculation calculationType) const
{
    if (calculationType == TangentFollowWellPathSegments)
    {
        cvf::Vec3d segmentStart, segmentEnd;
        m_wellPath->twoClosestPoints(m_intersections[intersectionIdx], &segmentStart, &segmentEnd);
        return (segmentEnd - segmentStart).getNormalized();
    }
    else
    {
        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());
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3f RigGeoMechWellLogExtractor::cellCentroid(size_t intersectionIdx) 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);
    int elementNodeCount = RigFemTypes::elmentNodeCount(elmType);

    const int* elmNodeIndices = femPart->connectivities(elmIdx);

    cvf::Vec3f centroid(0.0, 0.0, 0.0);
    for (int i = 0; i < elementNodeCount; ++i)
    {
        centroid += nodeCoords[elmNodeIndices[i]];
    }
    return centroid / elementNodeCount;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigGeoMechWellLogExtractor::getWellLogSegmentValue(size_t intersectionIdx, const std::vector<std::pair<double, double>>& wellLogValues) const
{
    double startMD, endMD;
    if (intersectionIdx % 2 == 0)
    {
        startMD = m_intersectionMeasuredDepths[intersectionIdx];
        endMD   = m_intersectionMeasuredDepths[intersectionIdx + 1];
    }
    else
    {
        startMD = m_intersectionMeasuredDepths [intersectionIdx - 1];
        endMD = m_intersectionMeasuredDepths[intersectionIdx];
    }

    RiaWeightedMeanCalculator<double> averageCalc;
    for (auto& depthAndValue : wellLogValues)
    {
        if (cvf::Math::valueInRange(depthAndValue.first, startMD, endMD))
        {
            cvf::Vec3d position = m_wellPath->interpolatedPointAlongWellPath(depthAndValue.first);
            cvf::Vec3d centroid(cellCentroid(intersectionIdx));
            double weight = 1.0;
            double dist = (position - centroid).length();
            if (dist > 1.0)
            {
                weight = 1.0 / dist;
            }
            averageCalc.addValueAndWeight(depthAndValue.second, weight);
        }
    }
    if (averageCalc.validAggregatedWeight())
    {
        return averageCalc.weightedMean();
    }

    return std::numeric_limits<double>::infinity();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigGeoMechWellLogExtractor::pascalToBar(double pascalValue)
{
    return pascalValue * 1.0e-5;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
template<typename T>
bool RigGeoMechWellLogExtractor::averageIntersectionValuesToSegmentValue(size_t intersectionIdx, const std::vector<T>& values, const T& invalidValue, T* averagedCellValue) const
{
    CVF_ASSERT(values.size() >= 2);
    
    *averagedCellValue = invalidValue;

    T value1, value2;
    cvf::Vec3d centroid(cellCentroid(intersectionIdx));
    double dist1 = 0.0, dist2 = 0.0;
    if (intersectionIdx % 2 == 0)
    {
        value1 = values[intersectionIdx];
        value2 = values[intersectionIdx + 1];

        dist1 = (centroid - m_intersections[intersectionIdx]).length();
        dist2 = (centroid - m_intersections[intersectionIdx + 1]).length();
    }
    else {
        value1 = values[intersectionIdx - 1];
        value2 = values[intersectionIdx];

        dist1 = (centroid - m_intersections[intersectionIdx - 1]).length();
        dist2 = (centroid - m_intersections[intersectionIdx]).length();
    }

    if (invalidValue == value1 || invalidValue == value2)
    {
        return false;
    }

    RiaWeightedMeanCalculator<T> averageCalc;
    averageCalc.addValueAndWeight(value1, dist2);
    averageCalc.addValueAndWeight(value2, dist1);
    if (averageCalc.validAggregatedWeight())
    {
        *averagedCellValue = averageCalc.weightedMean();
    }
    return true;
}