(#399) WellLog extraction from GeoMech cases works.

Needs more testing, and performance must be fixed see (#415)

ApplicationCode/ProjectDataModel/RimWellLogExtractionCurve.cpp

@@ -133,7 +133,6 @@ void RimWellLogExtractionCurve::updatePlotData()
 
         std::vector<double> filteredValues;
         std::vector<double> filteredDepths;
-        bool hasData = false;
 
         if (eclExtractor.notNull())
         {
@@ -153,32 +152,28 @@ void RimWellLogExtractionCurve::updatePlotData()
             if (resAcc.notNull())
             {
                 eclExtractor->curveData(resAcc.p(), &values);
-                hasData = true;
             }
 
-            // Remove values that are too difficult for Qwt to handle
+            filterPlotValues( depthValues,  filteredDepths, 
-            
+                              values,       filteredValues );
-            filteredValues.reserve(values.size());
-            filteredDepths.reserve(values.size());
-            for (size_t vIdx = 0; vIdx < values.size(); ++vIdx)
-            {
-                if (values[vIdx] == HUGE_VAL || values[vIdx] == -HUGE_VAL || (values[vIdx] != values[vIdx]))
-                {
-                    continue;
-                }
 
-                filteredDepths.push_back(depthValues[vIdx]);
-                filteredValues.push_back(values[vIdx]);
-            }
         }
-        else if (false) // geomExtractor
+        else if (geomExtractor.notNull()) // geomExtractor
         {
-            // Todo: do geomech log extraction
+            const std::vector<double>& depthValues = geomExtractor->measuredDepth();
+            m_geomResultDefinition->loadResult();
+            std::vector<double> values;
+            geomExtractor->curveData(m_geomResultDefinition->resultAddress() , m_timeStep, &values);
+            
+
+            filterPlotValues( depthValues,  filteredDepths, 
+                              values,       filteredValues );
+
         }
 
         m_plotCurve->setSamples(filteredValues.data(), filteredDepths.data(), (int)filteredValues.size());
 
-        if (hasData)
+        if (filteredValues.size())
         {
             RimWellLogPlot* wellLogPlot;
             firstAnchestorOrThisOfType(wellLogPlot);
@@ -319,4 +314,26 @@ void RimWellLogExtractionCurve::defineUiTreeOrdering(caf::PdmUiTreeOrdering& uiT
     uiTreeOrdering.setForgetRemainingFields(true);
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void RimWellLogExtractionCurve::filterPlotValues(const std::vector<double>& depthValues, std::vector<double> &filteredDepths, 
+                                                 const std::vector<double> &values,      std::vector<double> &filteredValues  )
+{
+    // Remove values that are too difficult for Qwt to handle
+
+    filteredValues.reserve(values.size());
+    filteredDepths.reserve(values.size());
+    for (size_t vIdx = 0; vIdx < values.size(); ++vIdx)
+    {
+        if (values[vIdx] == HUGE_VAL || values[vIdx] == -HUGE_VAL || (values[vIdx] != values[vIdx]))
+        {
+            continue;
+        }
+
+        filteredDepths.push_back(depthValues[vIdx]);
+        filteredValues.push_back(values[vIdx]);
+    }
+}
+
 

ApplicationCode/ProjectDataModel/RimWellLogExtractionCurve.h

@@ -52,6 +52,9 @@ protected:
 
     virtual void defineUiTreeOrdering(caf::PdmUiTreeOrdering& uiTreeOrdering, QString uiConfigName = "");
 
+    void filterPlotValues(const std::vector<double>& depthValues, std::vector<double> &filteredDepths, 
+                          const std::vector<double> &values,      std::vector<double> &filteredValues );
+
     caf::PdmPtrField<RimWellPath*>                  m_wellPath;
     caf::PdmPtrField<RimCase*>                      m_case;
     caf::PdmChildField<RimEclipseResultDefinition*> m_eclipseResultDefinition;

ApplicationCode/ProjectDataModel/RimWellLogPlotCollection.cpp

@@ -80,24 +80,24 @@ RigEclipseWellLogExtractor* RimWellLogPlotCollection::findOrCreateExtractor(RimW
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-RigGeoMechWellLogExtractor* RimWellLogPlotCollection::findOrCreateExtractor(RimWellPath* wellPath, RimGeoMechCase* eclCase)
+RigGeoMechWellLogExtractor* RimWellLogPlotCollection::findOrCreateExtractor(RimWellPath* wellPath, RimGeoMechCase* geomCase)
 {
-    if (!(wellPath && eclCase && wellPath->wellPathGeometry() && eclCase->geoMechData()))
+    if (!(wellPath && geomCase && wellPath->wellPathGeometry() && geomCase->geoMechData()))
     {
         return NULL;
     }
 
-    RigGeoMechCaseData* eclCaseData = eclCase->geoMechData();
+    RigGeoMechCaseData* geomCaseData = geomCase->geoMechData();
     RigWellPath* wellPathGeom = wellPath->wellPathGeometry();
     for (size_t exIdx = 0; exIdx < m_geomExtractors.size(); ++exIdx)
     {
-         if (m_geomExtractors[exIdx]->caseData() == eclCaseData && m_geomExtractors[exIdx]->wellPathData() == wellPathGeom)
+         if (m_geomExtractors[exIdx]->caseData() == geomCaseData && m_geomExtractors[exIdx]->wellPathData() == wellPathGeom)
          {
             return m_geomExtractors[exIdx].p();
          }
     }
 
-    cvf::ref<RigGeoMechWellLogExtractor> extractor = new RigGeoMechWellLogExtractor(eclCaseData, wellPathGeom);
+    cvf::ref<RigGeoMechWellLogExtractor> extractor = new RigGeoMechWellLogExtractor(geomCaseData, wellPathGeom);
     m_geomExtractors.push_back(extractor.p());
 
     return extractor.p();

ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp

@@ -43,12 +43,17 @@ RigGeoMechWellLogExtractor::RigGeoMechWellLogExtractor(RigGeoMechCaseData* aCase
 void RigGeoMechWellLogExtractor::curveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values)
 {   
     CVF_TIGHT_ASSERT(values);
-    values->resize(m_intersections.size());// + 1); // Plus one for the end of the wellpath stopping inside a cell
+    
+    if (!resAddr.isValid()) return ;
 
     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(resAddr, 0, frameIndex);
 
+    if (!resultValues.size()) return;
+
+    values->resize(m_intersections.size());// + 1); // Plus one for the end of the wellpath stopping inside a cell
+
     for (size_t cpIdx = 0; cpIdx < m_intersections.size(); ++cpIdx)
     {
         size_t elmIdx = m_intersectedCells[cpIdx];
@@ -127,11 +132,11 @@ void RigGeoMechWellLogExtractor::calculateIntersection()
         std::vector<HexIntersectionInfo> intersections;
 
         cvf::Vec3d hexCorners[8];
-        for (size_t cIdx = 0; cIdx < closeCells.size(); ++cIdx)
+        for (size_t ccIdx = 0; ccIdx < closeCells.size(); ++ccIdx)
         {
-            if (femPart->elementType(cIdx) != HEX8) continue;
+            if (femPart->elementType(closeCells[ccIdx]) != HEX8) continue;
 
-            const int* cornerIndices = femPart->connectivities(cIdx);
+            const int* cornerIndices = femPart->connectivities(closeCells[ccIdx]);
 
             hexCorners[0] = cvf::Vec3d(nodeCoords[cornerIndices[0]]);
             hexCorners[1] = cvf::Vec3d(nodeCoords[cornerIndices[1]]);
@@ -142,7 +147,7 @@ void RigGeoMechWellLogExtractor::calculateIntersection()
             hexCorners[6] = cvf::Vec3d(nodeCoords[cornerIndices[6]]);
             hexCorners[7] = cvf::Vec3d(nodeCoords[cornerIndices[7]]);
 
-            int intersectionCount = RigHexIntersector::lineHexCellIntersection(p1, p2, hexCorners, closeCells[cIdx], &intersections);
+            int intersectionCount = RigHexIntersector::lineHexCellIntersection(p1, p2, hexCorners, closeCells[ccIdx], &intersections);
         }
 
         // Now, with all the intersections of this piece of line, we need to 
@@ -193,14 +198,17 @@ std::vector<size_t> RigGeoMechWellLogExtractor::findCloseCells(const cvf::Boundi
         for (size_t elmIdx = 0; elmIdx < elmCount; ++elmIdx)
         {
             const int* elmNodeIndices = femPart->connectivities(elmIdx);
-            int elmNodeCount = RigFemTypes::elmentNodeCount( femPart->elementType(elmIdx));
+            int elmNodeCount = RigFemTypes::elmentNodeCount(femPart->elementType(elmIdx));
+            cvf::BoundingBox cellBB;
+
             for (int enIdx = 0; enIdx < elmNodeCount; ++enIdx)
             {
-                if (bb.contains(cvf::Vec3d(nodeCoords[elmNodeIndices[enIdx]])))
+                cellBB.add(cvf::Vec3d(nodeCoords[elmNodeIndices[enIdx]]));
-                {
+            }
-                    closeCells.push_back(elmIdx);
+
-                    break;
+            if (bb.intersects(cellBB))
-                }
+            {
+                closeCells.push_back(elmIdx);
             }
         }
     }

ApplicationCode/ReservoirDataModel/RigMainGrid.cpp

@@ -23,6 +23,7 @@
 #include "cvfAssert.h"
 #include "RimDefines.h"
 #include "RigFault.h"
+#include "cvfBoundingBoxTree.h"
 
 
 RigMainGrid::RigMainGrid(void)
@@ -441,3 +442,38 @@ const RigFault* RigMainGrid::findFaultFromCellIndexAndCellFace(size_t reservoirC
 #endif
     return NULL;
 }
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void RigMainGrid::findIntersectingCells(const cvf::BoundingBox& inputBB, std::vector<size_t>* cellIndices)
+{
+    if (m_cellSearchTree.isNull())
+    {
+        // build tree
+
+        size_t cellCount = m_cells.size();
+
+        std::vector<cvf::BoundingBox> cellBoundingBoxes;
+        cellBoundingBoxes.resize(cellCount);
+
+        for (size_t cIdx = 0; cIdx < cellCount; ++cIdx)
+        {
+            const caf::SizeTArray8& cellIndices = m_cells[cIdx].cornerIndices();
+            cvf::BoundingBox& cellBB = cellBoundingBoxes[cIdx];
+            cellBB.add(m_nodes[cellIndices[0]]);
+            cellBB.add(m_nodes[cellIndices[1]]);
+            cellBB.add(m_nodes[cellIndices[2]]);
+            cellBB.add(m_nodes[cellIndices[3]]);
+            cellBB.add(m_nodes[cellIndices[4]]);
+            cellBB.add(m_nodes[cellIndices[5]]);
+            cellBB.add(m_nodes[cellIndices[6]]);
+            cellBB.add(m_nodes[cellIndices[7]]);
+        }
+
+        m_cellSearchTree = new cvf::BoundingBoxTree;
+        m_cellSearchTree->buildTreeFromBoundingBoxes(cellBoundingBoxes, NULL);
+    }
+
+    m_cellSearchTree->findIntersections(inputBB, cellIndices);
+}

ApplicationCode/ReservoirDataModel/RigMainGrid.h

@@ -30,6 +30,11 @@
 #include <QtGlobal>
 #include "RigNNCData.h"
 
+namespace cvf
+{
+    class BoundingBoxTree;
+}
+
 class RigMainGrid : public RigGridBase
 {
 public:
@@ -63,7 +68,8 @@ public:
     void                                    setDisplayModelOffset(cvf::Vec3d offset);
 
     void                                    setFlipAxis(bool flipXAxis, bool flipYAxis);
-
+    void                                    findIntersectingCells(const cvf::BoundingBox& inputBB, std::vector<size_t>* cellIndices);
+                 
 private:
     void                                    initAllSubGridsParentGridPointer();
     void                                    initAllSubCellsMainGridCellIndex();
@@ -81,6 +87,7 @@ private:
     cvf::ref<RigFaultsPrCellAccumulator>    m_faultsPrCellAcc;
 
     cvf::Vec3d                              m_displayModelOffset;
+    cvf::ref<cvf::BoundingBoxTree>          m_cellSearchTree;    
 
     bool                                    m_flipXAxis;
     bool                                    m_flipYAxis;