#3114 compdat export. Improve performance

ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp

@@ -142,10 +142,14 @@ std::vector<RigCompletionData>
     // To handle several fractures in the same eclipse cell we need to keep track of the transmissibility
     // to the well from each fracture intersecting the cell and sum these transmissibilities at the end.
     // std::map <eclipseCellIndex ,map< fracture, trans> >
-    std::map<size_t, std::map<const RimFracture*, double>> eclCellIdxToTransPrFractureMap;
+    //std::map<size_t, std::map<const RimFracture*, double>> eclCellIdxToTransPrFractureMap;
 
-    for (const RimFracture* fracture : fractures)
+    std::vector<std::vector<RigCompletionData>> sharedComplForFracture(fractures.size());
+
+#pragma omp parallel for
+    for (int i = 0; i < (int)fractures.size(); i++)
     {
+        RimFracture* fracture = fractures[i];
         RimFractureTemplate* fracTemplate = fracture->fractureTemplate();
 
         if (!fracTemplate) continue;
@@ -320,7 +324,7 @@ std::vector<RigCompletionData>
                 double trans = transCondenser.condensedTransmissibility(
                     externalCell, {true, RigTransmissibilityCondenser::CellAddress::WELL, 1});
 
-                eclCellIdxToTransPrFractureMap[externalCell.m_globalCellIdx][fracture] = trans;
+                //eclCellIdxToTransPrFractureMap[externalCell.m_globalCellIdx][fracture] = trans;
 
                 RigCompletionData compDat(wellPathName,
                                           RigCompletionDataGridCell(externalCell.m_globalCellIdx, caseToApply->mainGrid()),
@@ -459,7 +463,8 @@ std::vector<RigCompletionData>
         }
 
         std::copy(
-            allCompletionsForOneFracture.begin(), allCompletionsForOneFracture.end(), std::back_inserter(fractureCompletions));
+            allCompletionsForOneFracture.begin(), allCompletionsForOneFracture.end(), std::back_inserter(sharedComplForFracture[i]));
+
 
         if (outputStreamForIntermediateResultsText)
         {
@@ -483,5 +488,9 @@ std::vector<RigCompletionData>
         }
     }
 
+    for (const auto& completions : sharedComplForFracture)
+    {
+        std::copy(completions.begin(), completions.end(), std::back_inserter(fractureCompletions));
+    }
     return fractureCompletions;
 }

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp

@@ -151,7 +151,7 @@ void RigCellGeometryTools::findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hex
 //--------------------------------------------------------------------------------------------------
 ///  
 //--------------------------------------------------------------------------------------------------
-double RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(std::vector<cvf::Vec3d> polygonToCalcLengthOf)
+double RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(const std::vector<cvf::Vec3d>& polygonToCalcLengthOf)
 {
     //Find bounding box
     cvf::BoundingBox polygonBBox;
@@ -250,19 +250,20 @@ cvf::Vec3d fromClipperPoint(const ClipperLib::IntPoint& clipPoint)
 //--------------------------------------------------------------------------------------------------
 ///  
 //--------------------------------------------------------------------------------------------------
-std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::intersectPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2)
+std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::intersectPolygons(const std::vector<cvf::Vec3d>& polygon1,
+                                                                              const std::vector<cvf::Vec3d>& polygon2)
 {
     std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
     
     // Convert to int for clipper library and store as clipper "path"
     ClipperLib::Path polygon1path;
-    for (cvf::Vec3d& v : polygon1)
+    for (const cvf::Vec3d& v : polygon1)
     {
         polygon1path.push_back(toClipperPoint(v));
     }
 
     ClipperLib::Path polygon2path;
-    for (cvf::Vec3d& v : polygon2)
+    for (const cvf::Vec3d& v : polygon2)
     {
         polygon2path.push_back(toClipperPoint(v));
     }
@@ -458,28 +459,30 @@ std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygo
 //--------------------------------------------------------------------------------------------------
 ///  
 //--------------------------------------------------------------------------------------------------
-std::pair<cvf::Vec3d, cvf::Vec3d> RigCellGeometryTools::getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine)
+std::pair<cvf::Vec3d, cvf::Vec3d> RigCellGeometryTools::getLineThroughBoundingBox(const cvf::Vec3d& lineDirection,
+                                                                                  const cvf::BoundingBox& polygonBBox,
+                                                                                  const cvf::Vec3d& pointOnLine)
 {
     cvf::Vec3d bboxCorners[8];
     polygonBBox.cornerVertices(bboxCorners);
 
     cvf::Vec3d startPoint = pointOnLine;
     cvf::Vec3d endPoint = pointOnLine;
-    
+    cvf::Vec3d lineDir = lineDirection;
 
     //To avoid doing many iterations in loops below linedirection should be quite large. 
-    lineDirection.normalize();
+    lineDir.normalize();
-    lineDirection = lineDirection * polygonBBox.extent().length() / 5;
+    lineDir = lineDir * polygonBBox.extent().length() / 5;
 
     //Extend line in positive direction
     while (polygonBBox.contains(startPoint))
     {
-        startPoint = startPoint + lineDirection;
+        startPoint = startPoint + lineDir;
     }
     //Extend line in negative direction
     while (polygonBBox.contains(endPoint))
     {
-        endPoint = endPoint - lineDirection;
+        endPoint = endPoint - lineDir;
     }
 
     std::pair<cvf::Vec3d, cvf::Vec3d> line;
@@ -508,7 +511,7 @@ double RigCellGeometryTools::getLengthOfPolygonAlongLine(const std::pair<cvf::Ve
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-std::vector<cvf::Vec3d> RigCellGeometryTools::unionOfPolygons(std::vector<std::vector<cvf::Vec3d>> polygons)
+std::vector<cvf::Vec3d> RigCellGeometryTools::unionOfPolygons(const std::vector<std::vector<cvf::Vec3d>>& polygons)
 {
     // Convert to int for clipper library and store as clipper "path"
     std::vector<ClipperLib::Path> polygonPaths;

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h

@@ -38,10 +38,10 @@ public:
 
     static void findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hexCorners, cvf::Vec3d &localXdirection, cvf::Vec3d &localYdirection, cvf::Vec3d &localZdirection);
 
-    static double polygonLengthInLocalXdirWeightedByArea(std::vector<cvf::Vec3d> polygon2d);
+    static double polygonLengthInLocalXdirWeightedByArea(const std::vector<cvf::Vec3d>& polygon2d);
 
-    static std::vector<std::vector<cvf::Vec3d>> intersectPolygons(std::vector<cvf::Vec3d> polygon1,
+    static std::vector<std::vector<cvf::Vec3d>> intersectPolygons(const std::vector<cvf::Vec3d>& polygon1,
-                                                                  std::vector<cvf::Vec3d> polygon2);
+                                                                  const std::vector<cvf::Vec3d>& polygon2);
     
     static std::vector<std::vector<cvf::Vec3d>> subtractPolygons(const std::vector<cvf::Vec3d>&              sourcePolygon,
                                                                  const std::vector<std::vector<cvf::Vec3d>>& polygonToSubtract);
@@ -51,11 +51,11 @@ public:
                                                                         const std::vector<cvf::Vec3d>& polygon, 
                                                                         ZInterpolationType interpolType = USE_ZERO);
 
-    static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine);
+    static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(const cvf::Vec3d& lineDirection, const cvf::BoundingBox& polygonBBox, const cvf::Vec3d& pointOnLine);
 
     static double getLengthOfPolygonAlongLine(const std::pair<cvf::Vec3d, cvf::Vec3d>& line, const std::vector<cvf::Vec3d>& polygon);
 
-    static std::vector<cvf::Vec3d> unionOfPolygons(std::vector<std::vector<cvf::Vec3d>> polygons);
+    static std::vector<cvf::Vec3d> unionOfPolygons(const std::vector<std::vector<cvf::Vec3d>>& polygons);
 
 private:
     static std::vector<cvf::Vec3d> ajustPolygonToAvoidIntersectionsAtVertex(const std::vector<cvf::Vec3d>& polyLine, 

ApplicationCode/ReservoirDataModel/RigHexIntersectionTools.cpp

@@ -114,12 +114,16 @@ bool RigHexIntersectionTools::isPointInCell(const cvf::Vec3d point,
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-bool RigHexIntersectionTools::planeHexCellIntersection(cvf::Vec3d* hexCorners, cvf::Plane fracturePlane, std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >& intersectionLineSegments)
+bool RigHexIntersectionTools::planeHexCellIntersection(cvf::Vec3d* hexCorners, const cvf::Plane& fracturePlane, std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >& intersectionLineSegments)
 {
     bool                                    isCellIntersected = false;
+    cvf::ubyte                              faceVertexIndices[4];
+    caf::HexGridIntersectionTools::ClipVx   triangleIntersectionPoint1;
+    caf::HexGridIntersectionTools::ClipVx   triangleIntersectionPoint2;
+    bool                                    isMostVxesOnPositiveSideOfP1 = false;
+
     for (int face = 0; face < 6; ++face)
     {
-        cvf::ubyte faceVertexIndices[4];
         cvf::StructGridInterface::cellFaceVertexIndices(static_cast<cvf::StructGridInterface::FaceType>(face), faceVertexIndices);
 
         cvf::Vec3d faceCenter = cvf::GeometryTools::computeFaceCenter(hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]]);
@@ -127,10 +131,6 @@ bool RigHexIntersectionTools::planeHexCellIntersection(cvf::Vec3d* hexCorners, c
         for (int i = 0; i < 4; i++)
         {
             int next = i < 3 ? i + 1 : 0;
-            caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint1;
-            caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint2;
-
-            bool isMostVxesOnPositiveSideOfP1 = false;
 
             bool isIntersectingPlane = caf::HexGridIntersectionTools::planeTriangleIntersection(fracturePlane,
                 hexCorners[faceVertexIndices[i]], 0,
@@ -141,7 +141,7 @@ bool RigHexIntersectionTools::planeHexCellIntersection(cvf::Vec3d* hexCorners, c
             if (isIntersectingPlane)
             {
                 isCellIntersected = true;
-                intersectionLineSegments.push_back({ triangleIntersectionPoint1.vx, triangleIntersectionPoint2.vx });
+                intersectionLineSegments.emplace_back( triangleIntersectionPoint1.vx, triangleIntersectionPoint2.vx );
             }
         }
     }    return isCellIntersected;

ApplicationCode/ReservoirDataModel/RigHexIntersectionTools.h

@@ -64,7 +64,7 @@ struct RigHexIntersectionTools
     static bool isPointInCell(const cvf::Vec3d point, const cvf::Vec3d hexCorners[8]);
 
     static bool planeHexCellIntersection(cvf::Vec3d* hexCorners,
-                                         cvf::Plane fracturePlane,
+                                         const cvf::Plane& fracturePlane,
                                          std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >& intersectionLineSegments);
 
     static bool planeHexIntersectionPolygons(std::array<cvf::Vec3d, 8> hexCorners,

ApplicationCode/UnitTests/RigHexIntersectionTools-Test.cpp

@@ -19,6 +19,8 @@
 #include "gtest/gtest.h"
 
 #include "RigHexIntersectionTools.h"
+#include "QDateTime"
+#include "QDebug"
 
 //--------------------------------------------------------------------------------------------------
 /// 
@@ -48,3 +50,48 @@ TEST(RigHexIntersectionTools, planeHexCellIntersectionTest)
     EXPECT_FALSE(isCellIntersected);
 
 }
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+TEST(RigHexIntersectionTools, DISABLED_planeHexCellIntersectionPerformanceTest)
+{
+    cvf::Vec3d hexCorners[8];
+    hexCorners[0] = cvf::Vec3d(0, 0, 0);
+    hexCorners[1] = cvf::Vec3d(1, 0, 0);
+    hexCorners[2] = cvf::Vec3d(0, 1, 0);
+    hexCorners[3] = cvf::Vec3d(0, 0, 1);
+    hexCorners[4] = cvf::Vec3d(0, 1, 1);
+    hexCorners[5] = cvf::Vec3d(1, 1, 0);
+    hexCorners[6] = cvf::Vec3d(1, 0, 1);
+    hexCorners[7] = cvf::Vec3d(1, 1, 1);
+
+    bool isCellIntersected = false;
+    cvf::Plane fracturePlaneIntersecting;
+    cvf::Plane fracturePlaneNotIntersecting;
+
+    fracturePlaneNotIntersecting.setFromPointAndNormal(cvf::Vec3d(1.5, 1.5, 1.5), cvf::Vec3d(1, 0, 0));
+    fracturePlaneIntersecting.setFromPointAndNormal(cvf::Vec3d(0.5, 0.5, 0.5), cvf::Vec3d(1, 0, 0));
+
+    QTime timeTotal;
+    timeTotal.start();
+
+    for (int run = 0; run < 5; run++)
+    {
+        QTime timeLocal;
+        timeLocal.start();
+
+        for (int i = 0; i < 2000000; i++)
+        {
+            std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;
+
+            isCellIntersected = RigHexIntersectionTools::planeHexCellIntersection(hexCorners, fracturePlaneIntersecting, intersectionLineSegments);
+            isCellIntersected = RigHexIntersectionTools::planeHexCellIntersection(hexCorners, fracturePlaneNotIntersecting, intersectionLineSegments);
+        }
+
+        qDebug() << "Time rim elapsed: " << timeLocal.elapsed();
+    }
+
+    qDebug() << "Time total elapsed: " << timeTotal.elapsed();
+}
+