#2393 Reservoir compaction derived result

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp

@@ -53,8 +53,31 @@
 
 #include <stdlib.h>
 #include <cmath>
+#include "RigHexIntersectionTools.h"
 
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+const std::string RigFemPartResultsCollection::FIELD_NAME_COMPACTION = "COMPACTION";
+
+//--------------------------------------------------------------------------------------------------
+/// Internal definitions
+//--------------------------------------------------------------------------------------------------
+class RefElement
+{
+public:
+    size_t              elementIdx;
+    float               intersectionPointToCurrentNodeDistance;
+    cvf::Vec3f          intersectionPoint;
+    std::vector<size_t> elementFaceNodeIdxs;
+};
+
+static std::vector<cvf::Vec3d> coordsFromNodeIndices(const RigFemPart& part, const std::vector<size_t>& nodeIdxs);
+static std::vector<size_t> nodesForElement(const RigFemPart& part, size_t elementIdx);
+static float horizontalDistance(const cvf::Vec3f& p1, const cvf::Vec3f& p2);
+static void findReferenceElementForNode(const RigFemPart& part, size_t nodeIdx, size_t kRefLayer, RefElement *refElement);
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -295,6 +318,7 @@ std::map<std::string, std::vector<std::string> > RigFemPartResultsCollection::sc
         {
             fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames();
             fieldCompNames["POR-Bar"];
+            fieldCompNames[FIELD_NAME_COMPACTION];
         }
         else if (resPos == RIG_ELEMENT_NODAL)
         {
@@ -559,7 +583,7 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateTimeLapseResult(
         frameCountProgress.setProgressDescription("Calculating " + QString::fromStdString(resVarAddr.fieldName + ": " + resVarAddr.componentName));
         frameCountProgress.setNextProgressIncrement(this->frameCount());
 
-        RigFemResultAddress resVarNative(resVarAddr.resultPosType, resVarAddr.fieldName, resVarAddr.componentName);
+        RigFemResultAddress resVarNative(resVarAddr.resultPosType, resVarAddr.fieldName, resVarAddr.componentName, RigFemResultAddress::NO_TIME_LAPSE, resVarAddr.refKLayerIndex);
         RigFemScalarResultFrames * srcDataFrames = this->findOrLoadScalarResult(partIndex, resVarNative);
         RigFemScalarResultFrames * dstDataFrames = m_femPartResults[partIndex]->createScalarResult(resVarAddr);
 
@@ -1427,6 +1451,63 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculatePrincipalStrainV
     return  requestedPrincipal;
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+RigFemScalarResultFrames* RigFemPartResultsCollection::calculateCompactionValues(int partIndex, const RigFemResultAddress &resVarAddr)
+{
+    CVF_ASSERT(resVarAddr.fieldName == FIELD_NAME_COMPACTION);
+
+    caf::ProgressInfo frameCountProgress(this->frameCount() + 1, "");
+    frameCountProgress.setProgressDescription("Calculating " + QString::fromStdString(resVarAddr.fieldName));
+
+    RigFemScalarResultFrames * u3Frames = this->findOrLoadScalarResult(partIndex, RigFemResultAddress(resVarAddr.resultPosType, "U", "U3"));
+    frameCountProgress.incrementProgress();
+
+    RigFemScalarResultFrames* compactionFrames = m_femPartResults[partIndex]->createScalarResult(resVarAddr);
+
+    const RigFemPart* part = m_femParts->part(partIndex);
+    for (int t = 0; t < u3Frames->frameCount(); t++)
+    {
+        std::vector<float>& compactionFrame = compactionFrames->frameData(t);
+        size_t nodeCount = part->nodes().nodeIds.size();
+
+        frameCountProgress.incrementProgress();
+
+        compactionFrame.resize(nodeCount);
+        for (int n = 0; n < nodeCount; n++)
+        {
+            RefElement refElement;
+            findReferenceElementForNode(*part, n, resVarAddr.refKLayerIndex, &refElement);
+
+            if (refElement.elementIdx != cvf::UNDEFINED_SIZE_T)
+            {
+                float shortestDist = HUGE_VALF;
+                size_t closestNodeIdx = cvf::UNDEFINED_SIZE_T;
+
+                for (size_t nodeIdx : refElement.elementFaceNodeIdxs)
+                {
+                    float dist = horizontalDistance(refElement.intersectionPoint, part->nodes().coordinates[nodeIdx]);
+                    if (dist < shortestDist)
+                    {
+                        shortestDist = dist;
+                        closestNodeIdx = nodeIdx;
+                    }
+                }
+
+                compactionFrame[n] = u3Frames->frameData(t)[n] - u3Frames->frameData(t)[closestNodeIdx];
+            }
+            else
+            {
+                compactionFrame[n] = HUGE_VAL;
+            }
+        }
+    }
+
+    RigFemScalarResultFrames* requestedPrincipal = this->findOrLoadScalarResult(partIndex, resVarAddr);
+    return  requestedPrincipal;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -1453,6 +1534,11 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult(in
         }
     }
 
+    if (resVarAddr.fieldName == FIELD_NAME_COMPACTION)
+    {
+        return calculateCompactionValues(partIndex, resVarAddr);
+    }
+
     if (resVarAddr.fieldName == "SE" && resVarAddr.componentName == "SFI")
     {
         return calculateSFI(partIndex, resVarAddr);    
@@ -2373,3 +2459,94 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar
         m_closestNodeId = femPart->nodes().nodeIds[closestNodeIdx];
     }
 }
+
+//--------------------------------------------------------------------------------------------------
+/// Internal functions
+//--------------------------------------------------------------------------------------------------
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+std::vector<cvf::Vec3d> coordsFromNodeIndices(const RigFemPart& part, const std::vector<size_t>& nodeIdxs)
+{
+    std::vector<cvf::Vec3d> out;
+    for (const auto& nodeIdx : nodeIdxs) out.push_back(cvf::Vec3d(part.nodes().coordinates[nodeIdx]));
+    return out;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+std::vector<size_t> nodesForElement(const RigFemPart& part, size_t elementIdx)
+{
+    std::vector<size_t> nodeIdxs;
+    const int* nodeConn = part.connectivities(elementIdx);
+    for (int n = 0; n < 8; n++) nodeIdxs.push_back(nodeConn[n]);
+    return nodeIdxs;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+float horizontalDistance(const cvf::Vec3f& p1, const cvf::Vec3f& p2)
+{
+    cvf::Vec3f p1_ = p1;
+    cvf::Vec3f p2_ = p2;
+    p1_.z() = p2_.z() = 0;
+    return p1_.pointDistance(p2_);
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void findReferenceElementForNode(const RigFemPart& part, size_t nodeIdx, size_t kRefLayer, RefElement *refElement)
+{
+    static const double zMin = -1e6, zMax = 1e6;
+
+    cvf::BoundingBox bb;
+    cvf::Vec3f currentNodeCoord = part.nodes().coordinates[nodeIdx];
+    cvf::Vec3f p1 = cvf::Vec3f(currentNodeCoord.x(), currentNodeCoord.y(), zMin);
+    cvf::Vec3f p2 = cvf::Vec3f(currentNodeCoord.x(), currentNodeCoord.y(), zMax);
+    bb.add(p1);
+    bb.add(p2);
+
+    std::vector<size_t> refElementCandidates;
+    part.findIntersectingCells(bb, &refElementCandidates);
+
+    const RigFemPartGrid* grid = part.structGrid();
+    const std::vector<cvf::Vec3f>& nodeCoords = part.nodes().coordinates;
+
+    refElement->elementIdx = cvf::UNDEFINED_SIZE_T;
+    refElement->intersectionPointToCurrentNodeDistance = HUGE_VALF;
+    size_t i, j, k;
+    for (const size_t elemIdx : refElementCandidates)
+    {
+        grid->ijkFromCellIndex(elemIdx, &i, &j, &k);
+        if (k == kRefLayer)
+        {
+            const std::vector<size_t> nodeIndices = nodesForElement(part, elemIdx);
+            CVF_ASSERT(nodeIndices.size() == 8);
+
+            std::vector<HexIntersectionInfo> intersections;
+            RigHexIntersectionTools::lineHexCellIntersection(cvf::Vec3d(p1), cvf::Vec3d(p2), coordsFromNodeIndices(part, nodeIndices).data(), elemIdx, &intersections);
+
+            for (const auto& intersection : intersections)
+            {
+                cvf::Vec3f intersectionPoint = cvf::Vec3f(intersection.m_intersectionPoint);
+
+                float nodeToIntersectionDistance = currentNodeCoord.pointDistance(intersectionPoint);
+                if (nodeToIntersectionDistance < refElement->intersectionPointToCurrentNodeDistance)
+                {
+                    cvf::ubyte faceNodes[4];
+                    grid->cellFaceVertexIndices(intersection.m_face, faceNodes);
+                    std::vector<size_t> topFaceCoords({ nodeIndices[faceNodes[0]], nodeIndices[faceNodes[1]], nodeIndices[faceNodes[2]], nodeIndices[faceNodes[3]] });
+
+                    refElement->elementIdx = elemIdx;
+                    refElement->intersectionPointToCurrentNodeDistance = nodeToIntersectionDistance;
+                    refElement->intersectionPoint = intersectionPoint;
+                    refElement->elementFaceNodeIdxs = topFaceCoords;
+                }
+            }
+        }
+    }
+}

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h

@@ -48,6 +48,8 @@ namespace caf
 class RigFemPartResultsCollection: public cvf::Object
 {
 public:
+    static const std::string FIELD_NAME_COMPACTION;
+
     RigFemPartResultsCollection(RifGeoMechReaderInterface* readerInterface, RifElementPropertyReader* elementPropertyReader, const RigFemPartCollection * femPartCollection);
     ~RigFemPartResultsCollection();
 
@@ -115,6 +117,7 @@ private:
     RigFemScalarResultFrames*                        calculateSurfaceAngles(int partIndex, const RigFemResultAddress& resVarAddr);
     RigFemScalarResultFrames*                        calculatePrincipalStressValues(int partIndex, const RigFemResultAddress &resVarAddr);
     RigFemScalarResultFrames*                        calculatePrincipalStrainValues(int partIndex, const RigFemResultAddress &resVarAddr);
+    RigFemScalarResultFrames*                        calculateCompactionValues(int partIndex, const RigFemResultAddress &resVarAddr);
 
     cvf::Collection<RigFemPartResults>               m_femPartResults;
     cvf::ref<RifGeoMechReaderInterface>              m_readerInterface;

ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h

@@ -36,24 +36,16 @@ public:
         componentName = "";
     }
 
-    RigFemResultAddress(RigFemResultPosEnum resPosType,
-                        const std::string& aFieldName,
-                        const std::string& aComponentName)
-        : resultPosType(resPosType),
-        fieldName(aFieldName),
-        componentName(aComponentName),
-        timeLapseBaseFrameIdx(-1)
-    {
-    }
-
     RigFemResultAddress(RigFemResultPosEnum resPosType,
                         const std::string& aFieldName,
                         const std::string& aComponentName,
-                        int aTimeLapseBaseFrame)
+                        int timeLapseBaseFrameIdx = NO_TIME_LAPSE,
+                        int refKLayerIndex = NO_COMPACTION)
         : resultPosType(resPosType),
         fieldName(aFieldName),
         componentName(aComponentName),
-        timeLapseBaseFrameIdx(aTimeLapseBaseFrame)
+        timeLapseBaseFrameIdx(timeLapseBaseFrameIdx),
+        refKLayerIndex(refKLayerIndex)
     {
     }
 
@@ -61,12 +53,17 @@ public:
     std::string         fieldName;
     std::string         componentName;
     int                 timeLapseBaseFrameIdx;
+    int                 refKLayerIndex;
 
     static const int ALL_TIME_LAPSES = -2;
+    static const int NO_TIME_LAPSE = -1;
+    static const int NO_COMPACTION = -1;
 
-    bool isTimeLapse() const { return timeLapseBaseFrameIdx >= 0;}
+    bool isTimeLapse() const { return timeLapseBaseFrameIdx > NO_TIME_LAPSE;}
     bool representsAllTimeLapses() const { return timeLapseBaseFrameIdx == ALL_TIME_LAPSES;}
 
+    bool isCompaction() const { return refKLayerIndex > NO_COMPACTION; }
+
     bool isValid() const
     {
         bool isTypeValid =     resultPosType == RIG_NODAL 
@@ -97,6 +94,11 @@ public:
             return (fieldName <  other.fieldName);
         }
  
+        if (refKLayerIndex != other.refKLayerIndex)
+        {
+            return refKLayerIndex < other.refKLayerIndex;
+        }
+
         return (componentName <  other.componentName);
   }
 

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp

@@ -24,6 +24,8 @@
 #include "RigFemPartResultsCollection.h"
 #include "RigFemResultAddress.h"
 #include "RigGeoMechCaseData.h"
+#include "RigFemPartCollection.h"
+#include "RigFemPartGrid.h"
 
 #include "RiaDefines.h"
 #include "RimGeoMechCase.h"
@@ -77,6 +79,12 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition(void)
     CAF_PDM_InitField(&m_timeLapseBaseTimestep, "TimeLapseBaseTimeStep", 0, "Base Time Step", "", "", "");
     m_timeLapseBaseTimestep.uiCapability()->setUiHidden(true);
 
+    CAF_PDM_InitField(&m_isCompaction, "IsCompaction", false, "Compaction", "", "", "");
+    m_isCompaction.uiCapability()->setUiHidden(true);
+
+    CAF_PDM_InitField(&m_compactionRefLayer, "CompactionRefLayer", 0, "Compaction Ref Layer", "", "", "");
+    m_compactionRefLayer.uiCapability()->setUiHidden(true);
+
     CAF_PDM_InitFieldNoDefault(&m_resultPositionTypeUiField, "ResultPositionTypeUi", "Result Position", "", "", "");
     m_resultPositionTypeUiField.xmlCapability()->setIOWritable(false);
     m_resultPositionTypeUiField.xmlCapability()->setIOReadable(false);
@@ -96,6 +104,14 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition(void)
     m_resultVariableUiField.uiCapability()->setUiEditorTypeName(caf::PdmUiListEditor::uiEditorTypeName());
     m_resultVariableUiField.uiCapability()->setUiLabelPosition(caf::PdmUiItemInfo::TOP);
 
+    CAF_PDM_InitField(&m_isCompactionUiField, "IsCompactionUi", false, "Enable Compaction", "", "", "");
+    m_isCompactionUiField.xmlCapability()->setIOWritable(false);
+    m_isCompactionUiField.xmlCapability()->setIOReadable(false);
+
+    CAF_PDM_InitField(&m_compactionRefLayerUiField, "CompactionRefLayerUi", 0, "Compaction Ref Layer", "", "", "");
+    m_compactionRefLayerUiField.xmlCapability()->setIOWritable(false);
+    m_compactionRefLayerUiField.xmlCapability()->setIOReadable(false);
+
     m_isChangedByField = false;
 }
 
@@ -123,6 +139,14 @@ void RimGeoMechResultDefinition::defineUiOrdering(QString uiConfigName, caf::Pdm
             timeLapseGr->add(&m_timeLapseBaseTimestepUiField);
     }
 
+    if (m_resultPositionTypeUiField() == RIG_NODAL)
+    {
+        caf::PdmUiGroup * compactionGroup = uiOrdering.addNewGroup("Compaction Options");
+        compactionGroup->add(&m_isCompactionUiField);
+        if (m_isCompactionUiField())
+            compactionGroup->add(&m_compactionRefLayerUiField);
+    }
+
     if (!m_isChangedByField)
     {
         m_resultPositionTypeUiField = m_resultPositionType;
@@ -156,6 +180,8 @@ QList<caf::PdmOptionItemInfo> RimGeoMechResultDefinition::calculateValueOptions(
 
             for (int oIdx = 0; oIdx < uiVarNames.size(); ++oIdx)
             {
+                if (uiVarNames[oIdx].startsWith(RigFemPartResultsCollection::FIELD_NAME_COMPACTION.data()) && !m_isCompactionUiField()) continue;
+
                 options.push_back(caf::PdmOptionItemInfo(uiVarNames[oIdx], varNames[oIdx]));
             }
         }
@@ -177,6 +203,14 @@ QList<caf::PdmOptionItemInfo> RimGeoMechResultDefinition::calculateValueOptions(
                 options.push_back(caf::PdmOptionItemInfo(QString::fromStdString(stepNames[stepIdx]), static_cast<int>(stepIdx)));
             }
         }
+        else if (&m_compactionRefLayerUiField == fieldNeedingOptions)
+        {
+            size_t kCount = m_geomCase->geoMechData()->femParts()->part(0)->structGrid()->gridPointCountK() - 1;
+            for (size_t layerIdx = 0; layerIdx < kCount; ++layerIdx)
+            {
+                options.push_back(caf::PdmOptionItemInfo(QString::number(layerIdx + 1), (int)layerIdx));
+            }
+        }
     }
 
     return options;
@@ -229,7 +263,7 @@ void RimGeoMechResultDefinition::fieldChangedByUi(const caf::PdmFieldHandle* cha
     RimPlotCurve* curve = nullptr;
     this->firstAncestorOrThisOfType(curve);
 
-    if (&m_resultVariableUiField == changedField)
+    if (&m_resultVariableUiField == changedField || &m_compactionRefLayerUiField == changedField)
     {
         QStringList fieldComponentNames = m_resultVariableUiField().split(QRegExp("\\s+"));
         if (fieldComponentNames.size() > 0)
@@ -257,6 +291,8 @@ void RimGeoMechResultDefinition::fieldChangedByUi(const caf::PdmFieldHandle* cha
 
                 m_isTimeLapseResult     = m_isTimeLapseResultUiField();
                 m_timeLapseBaseTimestep = m_timeLapseBaseTimestepUiField();
+                m_isCompaction          = m_isCompactionUiField();
+                m_compactionRefLayer    = m_compactionRefLayerUiField();
             }
 
             if (m_geomCase->geoMechData()->femPartResults()->assertResultsLoaded(this->resultAddress()))
@@ -378,6 +414,8 @@ void RimGeoMechResultDefinition::initAfterRead()
     m_resultVariableUiField = composeFieldCompString(m_resultFieldName(), m_resultComponentName());
     m_isTimeLapseResultUiField = m_isTimeLapseResult;
     m_timeLapseBaseTimestepUiField = m_timeLapseBaseTimestep;
+    m_isCompactionUiField = m_isCompaction;
+    m_compactionRefLayerUiField = m_compactionRefLayer;
 }
 
 
@@ -392,16 +430,16 @@ void RimGeoMechResultDefinition::loadResult()
     }
 }
 
-
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
 RigFemResultAddress RimGeoMechResultDefinition::resultAddress()
 {
-    if(m_isTimeLapseResult) 
+    return RigFemResultAddress(resultPositionType(),
-        return RigFemResultAddress(resultPositionType(), resultFieldName().toStdString(), resultComponentName().toStdString(), m_timeLapseBaseTimestep);
+                               resultFieldName().toStdString(), 
-    else                     
+                               resultComponentName().toStdString(),
-        return RigFemResultAddress(resultPositionType(), resultFieldName().toStdString(), resultComponentName().toStdString());
+                               m_isTimeLapseResult() ? m_timeLapseBaseTimestep() : RigFemResultAddress::NO_TIME_LAPSE,
+                               m_isCompaction() && resultFieldName().toStdString() == RigFemPartResultsCollection::FIELD_NAME_COMPACTION ? m_compactionRefLayer() : RigFemResultAddress::NO_COMPACTION);
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -479,11 +517,15 @@ void RimGeoMechResultDefinition::setResultAddress( const RigFemResultAddress& re
     m_resultFieldName       = QString::fromStdString(resultAddress.fieldName);
     m_resultComponentName   = QString::fromStdString(resultAddress.componentName);
     m_isTimeLapseResult     = resultAddress.isTimeLapse();
+    m_isCompaction          = resultAddress.isCompaction();
 
-    m_timeLapseBaseTimestep = m_isTimeLapseResult ? resultAddress.timeLapseBaseFrameIdx: 0;
+    m_timeLapseBaseTimestep = m_isTimeLapseResult ? resultAddress.timeLapseBaseFrameIdx: -1;
+    m_compactionRefLayer = m_isCompaction ? resultAddress.refKLayerIndex : -1;
 
     m_resultPositionTypeUiField = m_resultPositionType;
     m_resultVariableUiField = composeFieldCompString(m_resultFieldName(), m_resultComponentName());
     m_isTimeLapseResultUiField = m_isTimeLapseResult;
     m_timeLapseBaseTimestepUiField = m_timeLapseBaseTimestep;
+    m_isCompactionUiField = m_isCompaction;
+    m_compactionRefLayerUiField = m_compactionRefLayer;
 }

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h

@@ -101,6 +101,8 @@ private:
     caf::PdmField<QString>                            m_resultComponentName;
     caf::PdmField<bool>                               m_isTimeLapseResult;
     caf::PdmField<int>                                m_timeLapseBaseTimestep;
+    caf::PdmField<bool>                               m_isCompaction;
+    caf::PdmField<int>                                m_compactionRefLayer;
 
     // UI Fields only
 
@@ -112,7 +114,8 @@ private:
     caf::PdmField<QString>                            m_resultVariableUiField;
     caf::PdmField<bool>                               m_isTimeLapseResultUiField;
     caf::PdmField<int>                                m_timeLapseBaseTimestepUiField;
-
+    caf::PdmField<bool>                               m_isCompactionUiField;
+    caf::PdmField<int>                                m_compactionRefLayerUiField;
     caf::PdmPointer<RimGeoMechCase>                   m_geomCase;
 
     bool                                              m_isChangedByField;