MswRollUp: Handle previous well result points that is inside current cell

This is an intermediate commit and does not compile
p4#: 22229

ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp

@@ -838,6 +838,11 @@ cvf::Vec3d interpolate3DPosition(const std::vector<SegmentPositionContribution>
         {
             return filteredPositions[i].m_connectionPosition;
         }
+        else if (distance < 1.0)
+        {
+            distance = 1.0;
+        }
+
 
         distance = 1.0 / distance;
         nominators[i] = distance;

ApplicationCode/ModelVisualization/RivWellPipesPartMgr.cpp

@@ -267,7 +267,7 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
                 pipeBranchesCellIds.back().push_back(*prevWellResPoint);
             }
 
-            // Loop over all the resultCells in the branch
+            // Loop over all the resultPoints in the branch
 
             const std::vector<RigWellResultPoint>& resBranchCells = resBranches[brIdx].m_branchResultPoints;
 
@@ -278,12 +278,17 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
 
                 const RigWellResultPoint& currentWellResPoint = resBranchCells[cIdx];
 
+                // Ignore invalid cells
+
                 if (!currentWellResPoint.isValid())
                 {
                     //CVF_ASSERT(false); // Some segments does not get anything yet.
                     continue;
                 }
 
+                // Add cl contribution for a geometrical resultPoint by adding exit point from previous cell, 
+                // and then the result point position 
+
                 if (!currentWellResPoint.isCell()) 
                 {
                     // Use the interpolated value of branch head
@@ -307,7 +312,7 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
 
                         cvf::Vec3d outOfPrevCell(centerPreviousCell);
 
-                        bool intersectionOk = prevCell.firstIntersectionPoint(rayToThisCell, &outOfPrevCell);
+                        int intersectionOk = prevCell.firstIntersectionPoint(rayToThisCell, &outOfPrevCell);
                         //CVF_ASSERT(intersectionOk);
                         //CVF_ASSERT(intersectionOk);
                         if ((currentPoint - outOfPrevCell).lengthSquared() > 1e-3)
@@ -318,18 +323,21 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
 
                     }
 
-                    pipeBranchesCLCoords.back().push_back(currentPoint);
+                    branchCLCoords.push_back(currentPoint);
-                    pipeBranchesCellIds.back().push_back(currentWellResPoint);
+                    branchCellIds.push_back(currentWellResPoint);
                     
                     prevWellResPoint = &currentWellResPoint;
 
                     continue;
                 }
 
+                //
+                // Handle currentWellResPoint as a real cell result points.
+                //
+
                 const RigCell& cell = rigReservoir->cellFromWellResultCell(currentWellResPoint);
 
                 // Check if this and the previous cells has shared faces
-              
 
                 cvf::StructGridInterface::FaceType sharedFace;
                 if (prevWellResPoint && prevWellResPoint->isCell() && rigReservoir->findSharedSourceFace(sharedFace, currentWellResPoint, *prevWellResPoint))
@@ -343,6 +351,7 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
                 else
                 {
                     // This and the previous cell does not share a face.
+                    // Then we need to calculate the exit of the previous cell, and the entry point into this cell
 
                     cvf::Vec3d centerPreviousCell(cvf::Vec3d::ZERO);
                     cvf::Vec3d centerThisCell = cell.center();
@@ -372,12 +381,12 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
                     if (   wellResults->isMultiSegmentWell()
                             || !isAutoDetectBranches
                             || (prevWellResPoint == whResCell) 
-                            || distanceToWellHeadIsLonger
+                            || distanceToWellHeadIsLonger)
-                        )
                     {
                         // Not starting a "display" branch for normal wells
-                        
+                        // Calculate the exit of the previous cell, and the entry point into this cell
-                        cvf::Vec3d intoThisCell(centerThisCell);
+
+                        cvf::Vec3d intoThisCell(centerThisCell); // Use cell center as default for "into" point.
 
                         if (prevWellResPoint && prevWellResPoint->isValid())
                         {
@@ -390,7 +399,8 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
 
                             // Intersect with the current cell to find a better entry point than the cell center
 
-                            cell.firstIntersectionPoint(rayToThisCell, &intoThisCell);
+                            int intersectionCount = cell.firstIntersectionPoint(rayToThisCell, &intoThisCell);
+                            bool isPreviousResPointInsideCurrentCell = (intersectionCount % 2); // Must intersect uneven times to be inside. (1 % 2 = 1)
 
                             // If we have a real previous cell, we need to go out of it, before entering this.
                             // That is: add a CL-point describing where it leaves the previous cell.
@@ -409,6 +419,13 @@ void RivWellPipesPartMgr::calculateWellPipeCenterline(  std::vector< std::vector
                                     branchCellIds.push_back(RigWellResultPoint());
                                 }
                             }
+                            else if (isPreviousResPointInsideCurrentCell)
+                            {
+                                // Since the previous point actually is inside this cell, 
+                                /// use that as the entry point into this cell
+                                intoThisCell = centerPreviousCell;
+                            }
+
                         }
 
                         branchCLCoords.push_back(intoThisCell);

ApplicationCode/ReservoirDataModel/RigCell.cpp

@@ -242,9 +242,11 @@ cvf::Vec3d RigCell::faceCenter(cvf::StructGridInterface::FaceType face) const
 
 //--------------------------------------------------------------------------------------------------
 /// Find the intersection between the cell and the ray. The point closest to the ray origin is returned
-/// if no intersection is found, the intersection point is untouched.
+/// in \a intersectionPoint, while the return value is the total number of intersections with the 24 triangles
+/// the cell is interpreted as.
+/// If no intersection is found, the intersection point is untouched.
 //--------------------------------------------------------------------------------------------------
-bool RigCell::firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersectionPoint) const
+int RigCell::firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersectionPoint) const
 {
     CVF_ASSERT(intersectionPoint != NULL);
 
@@ -254,6 +256,7 @@ bool RigCell::firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersecti
 
     cvf::Vec3d firstIntersection(cvf::Vec3d::ZERO);
     double minLsq = HUGE_VAL;
+    int intersectionCount = 0; 
 
     for (face = 0; face < 6 ; ++face)
     {
@@ -261,9 +264,6 @@ bool RigCell::firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersecti
         cvf::Vec3d intersection;
         cvf::Vec3d faceCenter = this->faceCenter(static_cast<cvf::StructGridInterface::FaceType>(face));
 
-        ray.triangleIntersect(nodes[m_cornerIndices[faceVertexIndices[0]]], 
-            nodes[m_cornerIndices[faceVertexIndices[1]]], faceCenter, &intersection);
-
         for (size_t i = 0; i < 4; ++i)
         {
             size_t next = i < 3 ? i+1 : 0;
@@ -272,6 +272,7 @@ bool RigCell::firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersecti
                                         faceCenter, 
                                         &intersection))
             {
+                intersectionCount++;
                 double lsq = (intersection - ray.origin() ).lengthSquared();
                 if (lsq < minLsq)
                 {
@@ -282,12 +283,11 @@ bool RigCell::firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersecti
         }
     }
 
-    if (minLsq != HUGE_VAL)
+    if (intersectionCount > 0)
     {
         *intersectionPoint = firstIntersection;
-        return true;
     }
 
-    return false;
+    return intersectionCount;
 }
 

ApplicationCode/ReservoirDataModel/RigCell.h

@@ -62,7 +62,7 @@ public:
 
     cvf::Vec3d              center() const;
     cvf::Vec3d              faceCenter(cvf::StructGridInterface::FaceType face) const;
-    bool                    firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersectionPoint) const;
+    int                     firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersectionPoint) const;
     bool                    isLongPyramidCell(double maxHeightFactor = 5, double nodeNearTolerance = 1e-3 ) const;
 private:
     caf::SizeTArray8        m_cornerIndices;