#1901 Use the more robust WellLogExtraction system for perforation interval intersection calculations.

ApplicationCode/Commands/CompletionCommands/RicFishbonesTransmissibilityCalculationFeatureImp.cpp

@@ -179,7 +179,8 @@ void RicFishbonesTransmissibilityCalculationFeatureImp::findFishboneImportedLate
     for (const RimFishboneWellPath* fishbonesPath : wellPath->fishbonesCollection()->wellPathCollection()->wellPaths())
     {
         std::vector<WellPathCellIntersectionInfo> intersectedCells = RigWellPathIntersectionTools::findCellsIntersectedByPath(settings.caseToApply->eclipseCaseData(), 
-                                                                                                                              fishbonesPath->coordinates());
+                                                                                                                              fishbonesPath->coordinates(),
+                                                                                                                              fishbonesPath->measuredDepths());
         for (auto& cell : intersectedCells)
         {
             if (std::find(wellPathCells.begin(), wellPathCells.end(), cell.cellIndex) != wellPathCells.end()) continue;
@@ -211,11 +212,12 @@ void RicFishbonesTransmissibilityCalculationFeatureImp::findMainWellBoreParts(st
     double wellPathEndMD = 0.0;
     if (wellPathMD.size() > 1) wellPathEndMD = wellPathMD.back();
 
-    std::vector<cvf::Vec3d> fishbonePerfWellPathCoords = wellPath->wellPathGeometry()->clippedPointSubset(wellPath->fishbonesCollection()->startMD(),
+    std::pair< std::vector<cvf::Vec3d>, std::vector<double> > fishbonePerfWellPathCoords = wellPath->wellPathGeometry()->clippedPointSubset(wellPath->fishbonesCollection()->startMD(),
-                                                                                                          wellPathEndMD);
+                                                                                                                                            wellPathEndMD);
 
     std::vector<WellPathCellIntersectionInfo> intersectedCellsIntersectionInfo = RigWellPathIntersectionTools::findCellsIntersectedByPath(settings.caseToApply->eclipseCaseData(),
-                                                                                                                                          fishbonePerfWellPathCoords);
+                                                                                                                                          fishbonePerfWellPathCoords.first,
+                                                                                                                                          fishbonePerfWellPathCoords.second);
 
     for (auto& cell : intersectedCellsIntersectionInfo)
     {

ApplicationCode/Commands/CompletionCommands/RicWellPathExportCompletionDataFeature.cpp

@@ -716,9 +716,11 @@ std::vector<RigCompletionData> RicWellPathExportCompletionDataFeature::generateP
     {
         if (!interval->isActiveOnDate(settings.caseToApply->timeStepDates()[settings.timeStep])) continue;
 
-        std::vector<cvf::Vec3d> perforationPoints = wellPath->wellPathGeometry()->clippedPointSubset(interval->startMD(), interval->endMD());
+        using namespace std;
+        pair<vector<cvf::Vec3d>, vector<double> > perforationPointsAndMD = wellPath->wellPathGeometry()->clippedPointSubset(interval->startMD(), interval->endMD());
         std::vector<WellPathCellIntersectionInfo> intersectedCells = RigWellPathIntersectionTools::findCellsIntersectedByPath(settings.caseToApply->eclipseCaseData(), 
-                                                                                                                              perforationPoints);
+                                                                                                                              perforationPointsAndMD.first,
+                                                                                                                              perforationPointsAndMD.second);
         for (auto& cell : intersectedCells)
         {
             bool cellIsActive = activeCellInfo->isActive(cell.cellIndex);
@@ -861,9 +863,24 @@ void RicWellPathExportCompletionDataFeature::calculateLateralIntersections(const
     for (WellSegmentLateral& lateral : location->laterals)
     {
         lateral.branchNumber = ++(*branchNum);
-        std::vector<cvf::Vec3d> lateralCoords = location->fishbonesSubs->coordsForLateral(location->subIndex, lateral.lateralIndex);
+
+        std::vector<std::pair<cvf::Vec3d, double> > lateralCoordMDPairs = location->fishbonesSubs->coordsAndMDForLateral(location->subIndex, lateral.lateralIndex);
+        
+        std::vector<cvf::Vec3d> lateralCoords;
+        std::vector<double> lateralMDs;
+
+        lateralCoords.reserve(lateralCoordMDPairs.size());
+        lateralMDs.reserve(lateralCoordMDPairs.size());
+
+        for (auto& coordMD : lateralCoordMDPairs)
+        {
+            lateralCoords.push_back(coordMD.first);
+            lateralMDs.push_back(coordMD.second);
+        }
+
         std::vector<WellPathCellIntersectionInfo> intersections = RigWellPathIntersectionTools::findCellsIntersectedByPath(caseToApply->eclipseCaseData(), 
-                                                                                                                           lateralCoords);
+                                                                                                                           lateralCoords,
+                                                                                                                           lateralMDs);
 
         auto intersection = intersections.cbegin();
         double length = 0;

ApplicationCode/ModelVisualization/RivWellPathPartMgr.cpp

@@ -200,16 +200,18 @@ void RivWellPathPartMgr::appendPerforationsToModel(const QDateTime& currentViewD
 
         if (currentViewDate.isValid() && !perforation->isActiveOnDate(currentViewDate)) continue;
 
-        std::vector<cvf::Vec3d> displayCoords = wellPathGeometry->clippedPointSubset(perforation->startMD(), perforation->endMD());
+        using namespace std;
+        pair<vector<cvf::Vec3d>, vector<double> >  displayCoordsAndMD = wellPathGeometry->clippedPointSubset(perforation->startMD(), 
+                                                                                                             perforation->endMD());
  
-        if (displayCoords.size() < 2) continue;
+        if (displayCoordsAndMD.first.size() < 2) continue;
 
-        for (cvf::Vec3d& point : displayCoords) point = displayCoordTransform->transformToDisplayCoord(point);
+        for (cvf::Vec3d& point : displayCoordsAndMD.first) point = displayCoordTransform->transformToDisplayCoord(point);
 
         cvf::ref<RivObjectSourceInfo> objectSourceInfo = new RivObjectSourceInfo(perforation);
 
         cvf::Collection<cvf::Part> parts;
-        geoGenerator.cylinderWithCenterLineParts(&parts, displayCoords, cvf::Color3f::GREEN, perforationRadius);
+        geoGenerator.cylinderWithCenterLineParts(&parts, displayCoordsAndMD.first, cvf::Color3f::GREEN, perforationRadius);
         for (auto part : parts)
         {
             part->setSourceInfo(objectSourceInfo.p());

ApplicationCode/ProjectDataModel/Completions/RimCompletionCellIntersectionCalc.cpp

@@ -148,9 +148,17 @@ void RimCompletionCellIntersectionCalc::calculateFishbonesIntersections(const Ri
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-void RimCompletionCellIntersectionCalc::calculatePerforationIntersections(const RimWellPath* wellPath, const RimPerforationInterval* perforationInterval, const RigMainGrid* grid, std::vector<double>& values)
+void RimCompletionCellIntersectionCalc::calculatePerforationIntersections(const RimWellPath* wellPath, 
+                                                                          const RimPerforationInterval* perforationInterval, 
+                                                                          const RigMainGrid* grid, 
+                                                                          std::vector<double>& values)
 {
-    std::vector<HexIntersectionInfo> intersections = RigWellPathIntersectionTools::getIntersectedCells(grid, wellPath->wellPathGeometry()->clippedPointSubset(perforationInterval->startMD(), perforationInterval->endMD()));
+    using namespace std;
+    pair<vector<cvf::Vec3d>, vector<double> > clippedWellPathData =  wellPath->wellPathGeometry()->clippedPointSubset(perforationInterval->startMD(), 
+                                                                                                                      perforationInterval->endMD());
+
+    std::vector<HexIntersectionInfo> intersections = RigWellPathIntersectionTools::getIntersectedCells(grid,  
+                                                                                                       clippedWellPathData.first);
     for (auto& intersection : intersections)
     {
         values[intersection.m_hexIndex] = RiaDefines::PERFORATION_INTERVAL;

ApplicationCode/ReservoirDataModel/RigEclipseWellLogExtractor.cpp

@@ -27,6 +27,7 @@
 
 #include "RigWellLogExtractionTools.h"
 #include "RigMainGrid.h"
+#include "RigWellPathIntersectionTools.h"
 
 //==================================================================================================
 /// 
@@ -131,6 +132,31 @@ void RigEclipseWellLogExtractor::curveData(const RigResultAccessor* resultAccess
    
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+std::vector<WellPathCellIntersectionInfo> RigEclipseWellLogExtractor::cellIntersectionInfo()
+{
+    std::vector<WellPathCellIntersectionInfo> cellIntersectionInfos;
+    cellIntersectionInfos.reserve(m_intersections.size()-1);
+   
+    for (size_t cpIdx = 0; cpIdx < m_intersections.size()-1; ++cpIdx)
+    {
+        size_t cellIdx1 = m_intersectedCells[cpIdx];
+        size_t cellIdx2 = m_intersectedCells[cpIdx+1];
+
+        if (cellIdx1 == cellIdx2)
+        {
+            cvf::Vec3d internalCellLengths;
+            internalCellLengths = RigWellPathIntersectionTools::calculateLengthInCell( m_caseData->mainGrid(), cellIdx1, m_intersections[cpIdx], m_intersections[cpIdx+1] );
+
+            cellIntersectionInfos.push_back(WellPathCellIntersectionInfo(cellIdx1, m_intersections[cpIdx], m_intersections[cpIdx+1], internalCellLengths));
+        }
+    }   
+
+    return  cellIntersectionInfos;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigEclipseWellLogExtractor.h

@@ -28,6 +28,7 @@
 
 #include <vector>
 #include "cvfStructGrid.h"
+#include "RigWellPathIntersectionTools.h"
 
 class RigEclipseCaseData;
 class RigWellPath;
@@ -48,6 +49,8 @@ public:
     void                        curveData(const RigResultAccessor* resultAccessor, std::vector<double>* values );
     const RigEclipseCaseData*          caseData()     { return m_caseData.p();}
 
+    std::vector<WellPathCellIntersectionInfo> cellIntersectionInfo();
+
 protected:
     void                        calculateIntersection();
     std::vector<size_t>         findCloseCells(const cvf::BoundingBox& bb);

ApplicationCode/ReservoirDataModel/RigWellPath.cpp

@@ -192,25 +192,29 @@ void RigWellPath::twoClosestPoints(const cvf::Vec3d& position, cvf::Vec3d* p1, c
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-std::vector<cvf::Vec3d> RigWellPath::clippedPointSubset(double startMD, double endMD) const
+std::pair<std::vector<cvf::Vec3d>, std::vector<double> > RigWellPath::clippedPointSubset(double startMD, double endMD) const
 {
-    std::vector<cvf::Vec3d> points;
+    std::pair<std::vector<cvf::Vec3d>, std::vector<double> >  pointsAndMDs;
-    if (m_measuredDepths.empty()) return points;
+    if (m_measuredDepths.empty()) return pointsAndMDs;
-    if (startMD > endMD) return points;
+    if (startMD > endMD) return pointsAndMDs;
+
+    pointsAndMDs.first.push_back(interpolatedPointAlongWellPath(startMD));
+    pointsAndMDs.second.push_back(startMD);
 
-    points.push_back(interpolatedPointAlongWellPath(startMD));
     for (size_t i = 0; i < m_measuredDepths.size(); ++i)
     {
         double measuredDepth = m_measuredDepths[i];
         if (measuredDepth > startMD && measuredDepth < endMD)
         {
-            points.push_back(m_wellPathPoints[i]);
+            pointsAndMDs.first.push_back(m_wellPathPoints[i]);
+            pointsAndMDs.second.push_back(measuredDepth);
         }
     }
-    points.push_back(interpolatedPointAlongWellPath(endMD));
+    pointsAndMDs.first.push_back(interpolatedPointAlongWellPath(endMD));
+    pointsAndMDs.second.push_back(endMD);
 
 
-    return points;
+    return pointsAndMDs;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigWellPath.h

@@ -43,7 +43,9 @@ public:
     cvf::Vec3d                  interpolatedPointAlongWellPath(double measuredDepth) const;
     double                      wellPathAzimuthAngle(const cvf::Vec3d& position) const;
     void                        twoClosestPoints(const cvf::Vec3d& position, cvf::Vec3d* p1, cvf::Vec3d* p2) const;
-    std::vector<cvf::Vec3d>     clippedPointSubset(double startMD, double endMD) const;
+
+    std::pair<std::vector<cvf::Vec3d>, std::vector<double> > 
+                                clippedPointSubset(double startMD, double endMD) const;
 
     std::vector<cvf::Vec3d>     wellPathPointsIncludingFractureIntersection(double fractureIntersectionMD) const;
 private:

ApplicationCode/ReservoirDataModel/RigWellPathIntersectionTools.cpp

@@ -28,68 +28,29 @@
 
 #include "cvfGeometryTools.h"
 #include "cvfMatrix3.h"
+#include "RigEclipseWellLogExtractor.h"
+#include "RimEclipseCase.h"
 
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
 std::vector<WellPathCellIntersectionInfo> RigWellPathIntersectionTools::findCellsIntersectedByPath(const RigEclipseCaseData* caseData, 
-                                                                                                   const std::vector<cvf::Vec3d>& pathCoords)
+                                                                                                   const std::vector<cvf::Vec3d>& pathCoords, 
+                                                                                                   const std::vector<double>& pathMds)
 {
     std::vector<WellPathCellIntersectionInfo> intersectionInfos;
     const RigMainGrid* grid = caseData->mainGrid();
 
     if (pathCoords.size() < 2) return intersectionInfos;
+    cvf::ref<RigWellPath> dummyWellPath = new RigWellPath;
+    dummyWellPath->m_wellPathPoints = pathCoords;
+    dummyWellPath->m_measuredDepths = pathMds;
 
-    std::vector<HexIntersectionInfo> intersections = getIntersectedCells(grid, pathCoords);
+    cvf::ref<RigEclipseWellLogExtractor> extractor = new RigEclipseWellLogExtractor(caseData, dummyWellPath.p(), caseData->ownerCase()->caseUserDescription().toStdString());
-    removeEnteringIntersections(&intersections);
 
-    if (intersections.empty()) return intersectionInfos;
+    return extractor->cellIntersectionInfo();
 
-    cvf::Vec3d            startPoint;
+ 
-    cvf::Vec3d            endPoint;
-    size_t                cellIndex;
-    cvf::Vec3d            internalCellLengths;
-
-    auto intersection = intersections.cbegin();
-
-    //start cell
-    bool foundCell;
-    startPoint = pathCoords[0];
-    cellIndex = findCellFromCoords(grid, startPoint, &foundCell);
-    if (foundCell)
-    {
-        endPoint = intersection->m_intersectionPoint;
-        internalCellLengths = calculateLengthInCell(grid, cellIndex, startPoint, endPoint);
-        intersectionInfos.push_back(WellPathCellIntersectionInfo(cellIndex, startPoint, endPoint, internalCellLengths));
-    }
-    else
-    {
-        RiaLogging::debug("Path starts outside valid cell");
-    }
-
-    //center cells
-    startPoint = intersection->m_intersectionPoint;
-    cellIndex = intersection->m_hexIndex;
-
-    ++intersection;
-
-    while (intersection != intersections.cend())
-    {
-        endPoint = intersection->m_intersectionPoint;
-        internalCellLengths = calculateLengthInCell(grid, cellIndex, startPoint, endPoint);
-        intersectionInfos.push_back(WellPathCellIntersectionInfo(cellIndex, startPoint, endPoint, internalCellLengths));
-
-        startPoint = endPoint;
-        cellIndex = intersection->m_hexIndex;
-        ++intersection;
-    }
-
-    //end cell
-    endPoint = pathCoords[pathCoords.size() - 1];
-    internalCellLengths = calculateLengthInCell(grid, cellIndex, startPoint, endPoint);
-    intersectionInfos.push_back(WellPathCellIntersectionInfo(cellIndex, startPoint, endPoint, internalCellLengths));
-
-    return intersectionInfos;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigWellPathIntersectionTools.h

@@ -53,7 +53,9 @@ struct WellPathCellIntersectionInfo {
 class RigWellPathIntersectionTools
 {
 public:
-    static std::vector<WellPathCellIntersectionInfo>   findCellsIntersectedByPath(const RigEclipseCaseData* caseData, const std::vector<cvf::Vec3d>& pathCoords);
+    static std::vector<WellPathCellIntersectionInfo>   findCellsIntersectedByPath(const RigEclipseCaseData* caseData, 
+                                                                                  const std::vector<cvf::Vec3d>& pathCoords,
+                                                                                  const std::vector<double>& pathMds);
 
     static std::vector<HexIntersectionInfo>            getIntersectedCells(const RigMainGrid* grid, const std::vector<cvf::Vec3d>& coords);