NNC processing first cut

ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp

@@ -392,15 +392,22 @@ bool RifReaderEclipseOutput::open(const QString& fileName, RigCaseData* eclipseC
 
     m_eclipseCase = eclipseCase;
 
-    progInfo.setProgressDescription("Reading Result index");
-    progInfo.setNextProgressIncrement(50);
-    
     // Build results meta data
+    progInfo.setProgressDescription("Reading Result index");
+    progInfo.setNextProgressIncrement(25);
     buildMetaData();
-
-    transferNNCData(mainEclGrid, m_ecl_init_file, eclipseCase->mainGrid());
-
     progInfo.incrementProgress();
+
+    progInfo.setProgressDescription("Reading NNC data");
+    progInfo.setNextProgressIncrement(5);
+    transferNNCData(mainEclGrid, m_ecl_init_file, eclipseCase->mainGrid());
+    progInfo.incrementProgress();
+
+    progInfo.setProgressDescription("Processing NNC's");
+    progInfo.setNextProgressIncrement(20);
+    eclipseCase->mainGrid()->nncData()->processConnections( *(eclipseCase->mainGrid()));
+    progInfo.incrementProgress();
+
     progInfo.setNextProgressIncrement(8);
     progInfo.setProgressDescription("Reading Well information");
     readWellCells(mainEclGrid);

ApplicationCode/ReservoirDataModel/RigNNCData.cpp

@@ -18,7 +18,7 @@
 
 #include "RigNNCData.h"
 #include "RigMainGrid.h"
-//#include "../ModelVisualization/cvfGeometryTools.h"
+#include "cvfGeometryTools.h"
 
 
 //--------------------------------------------------------------------------------------------------
@@ -35,14 +35,16 @@ RigNNCData::RigNNCData()
 //--------------------------------------------------------------------------------------------------
 void RigNNCData::processConnections(const RigMainGrid& mainGrid)
 {
-    /*
+
-    for (size_t cnIdx = 0; cnIdx < 0; ++cnIdx)
+    for (size_t cnIdx = 0; cnIdx < m_connections.size(); ++cnIdx)
     {
         const RigCell& c1 = mainGrid.cells()[m_connections[cnIdx].m_c1GlobIdx];
         const RigCell& c2 = mainGrid.cells()[m_connections[cnIdx].m_c2GlobIdx];
 
         // Try to find the shared face
+
         char hasNeighbourInAnyDirection = 0;
+        bool isPossibleNeighborInDirection[6]= {true, true, true, true, true, true};
 
         if (c1.hostGrid() == c2.hostGrid())
         {
@@ -51,7 +53,7 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid)
             size_t i2, j2, k2;
             c2.hostGrid()->ijkFromCellIndex(c2.cellIndex(), &i2, &j2, &k2);
 
-            bool isPossibleNeighborInDirection[6]= {false, false, false, false, false, false};
+          
             isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] = ((i1 + 1) == i2);
             isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_I] = ((i2 + 1) == i1);
             isPossibleNeighborInDirection[cvf::StructGridInterface::POS_J] = ((j1 + 1) == j2);
@@ -60,7 +62,7 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid)
             isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_K] = ((k2 + 1) == k1);
 
             hasNeighbourInAnyDirection = 
-              isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] 
+                isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] 
             + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_I]
             + isPossibleNeighborInDirection[cvf::StructGridInterface::POS_J]
             + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_J]
@@ -73,68 +75,95 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid)
 
             if (!hasNeighbourInAnyDirection)
             {
-                m_connections[cnIdx].m_hasNoSharedArea = true;
+                // Add to search map
+                m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c1GlobIdx][cvf::StructGridInterface::NO_FACE].push_back(cnIdx);
+                m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c2GlobIdx][cvf::StructGridInterface::NO_FACE].push_back(cnIdx);
+
                 continue; // to next connection
             }
+        }
 
-            if (hasNeighbourInAnyDirection == 1)
+        // Possibly do some testing to avoid unneccesary overlap calculations
+
+        cvf::Vec3d normal;
+        for (char fIdx = 0; fIdx < 6; ++fIdx)
+        {
+            if (isPossibleNeighborInDirection[fIdx])
             {
-                for (char fIdx = 0; fIdx < 6; ++fIdx)
+                cvf::Vec3d fc1 = c1.faceCenter((cvf::StructGridInterface::FaceType)(fIdx));
+                cvf::Vec3d fc2 = c2.faceCenter(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)));
+                cvf::Vec3d fc1ToFc2 = fc2 - fc1;
+                normal = c1.faceNormal((cvf::StructGridInterface::FaceType)(fIdx));
+                normal.normalize();
+                // Check that face centers are approx in the face plane
+                if (normal.dot(fc1ToFc2) < 0.01*fc1ToFc2.length()) 
                 {
-                    if (isPossibleNeighborInDirection[fIdx])
+
-                    {
-                        m_connections[cnIdx].m_c1Face = (cvf::StructGridInterface::FaceType)fIdx;
-                        break; // the face loop
-                    }
                 }
-
-                // calculate polygon for this face
-
-            }
-            else
-            {
-
-                cvf::Vec3d normal;
-                for (char fIdx = 0; fIdx < 6; ++fIdx)
-                {
-                    if (isPossibleNeighborInDirection[fIdx])
-                    {
-                        cvf::Vec3d fc1 = c1.faceCenter((cvf::StructGridInterface::FaceType)(fIdx));
-                        cvf::Vec3d fc2 = c2.faceCenter(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)));
-                        cvf::Vec3d fc1ToFc2 = fc2 - fc1;
-                        normal = c1.faceNormal((cvf::StructGridInterface::FaceType)(fIdx));
-                        normal.normalize();
-                        // Check that face centers are approx in the face plane
-                        if (normal.dot(fc1ToFc2) < 0.01*fc1ToFc2.length()) 
-                        {
-
-                        }
-
-                        // Calculate connection polygon
-
-                        std::vector<size_t> polygon;
-                        std::vector<cvf::Vec3d> intersections;
-                        bool isOk = false;
-                        caf::SizeTArray4 face1;
-                        caf::SizeTArray4 face2;
-                        c1.faceIndices((cvf::StructGridInterface::FaceType)(fIdx), &face1);
-                        c2.faceIndices(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)), &face2);
-
-                        isOk = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads(
-                            &polygon, 
-                            &intersections, 
-                            (cvf::EdgeIntersectStorage<size_t>*)NULL, 
-                            cvf::wrapArrayConst(&mainGrid.nodes()), 
-                            face1.data(), 
-                            face2.data(), 
-                            1e-6);
-
-                    }
-
-                    }
-                }
-
             }
         }
-        */
+
+        for (char fIdx = 0; fIdx < 6; ++fIdx)
+        {
+            if (!isPossibleNeighborInDirection[fIdx])
+            { 
+                continue;
+            }
+
+            // Calculate connection polygon
+
+            std::vector<size_t> polygon;
+            std::vector<cvf::Vec3d> intersections;
+            caf::SizeTArray4 face1;
+            caf::SizeTArray4 face2;
+            c1.faceIndices((cvf::StructGridInterface::FaceType)(fIdx), &face1);
+            c2.faceIndices(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)), &face2);
+
+            bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads(
+                &polygon, 
+                &intersections, 
+                (cvf::EdgeIntersectStorage<size_t>*)NULL, 
+                cvf::wrapArrayConst(&mainGrid.nodes()), 
+                face1.data(), 
+                face2.data(), 
+                1e-6);
+
+            if (foundOverlap)
+            {
+                // Found an overlap polygon. Store data about connection
+
+                m_connections[cnIdx].m_c1Face = (cvf::StructGridInterface::FaceType)fIdx;
+                for (size_t pIdx = 0; pIdx < polygon.size(); ++pIdx)
+                {
+                    if (polygon[pIdx] < mainGrid.nodes().size())
+                        m_connections[cnIdx].m_polygon.push_back(mainGrid.nodes()[polygon[pIdx]]);
+                    else
+                        m_connections[cnIdx].m_polygon.push_back(intersections[polygon[pIdx] - mainGrid.nodes().size()]);
+                }
+
+                // Add to search map
+                m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c1GlobIdx][fIdx].push_back(cnIdx);
+                m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c2GlobIdx][cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx))].push_back(cnIdx);
+
+                break; // The connection face is found. Stop looping over the cell faces. Jump to next connection
+            }
+        }
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+const std::vector<size_t>& RigNNCData::findConnectionIndices( size_t globalCellIndex, cvf::StructGridInterface::FaceType face) const
+{
+    ConnectionSearchMap::const_iterator it;
+    static std::vector<size_t> empty;
+
+    it = m_cellIdxToFaceToConnectionIdxMap.find(globalCellIndex);
+    if (it != m_cellIdxToFaceToConnectionIdxMap.end())
+    {
+        return it->second[face];
+    }
+
+    return empty;
 }

ApplicationCode/ReservoirDataModel/RigNNCData.h

@@ -35,22 +35,19 @@ class RigConnection
 public:
     RigConnection( ) 
         : m_c1GlobIdx(cvf::UNDEFINED_SIZE_T),
-          m_c1Face(cvf::StructGridInterface::POS_I),
+          m_c1Face(cvf::StructGridInterface::NO_FACE),
           m_c2GlobIdx(cvf::UNDEFINED_SIZE_T),
-          m_c2Face(cvf::StructGridInterface::NEG_I),
-          m_hasNoSharedArea(false),
           m_transmissibility(0.0)
     {}
 
-
-
-    bool                                m_hasNoSharedArea;
     size_t                              m_c1GlobIdx;
     cvf::StructGridInterface::FaceType  m_c1Face;
     size_t                              m_c2GlobIdx;
-    cvf::StructGridInterface::FaceType  m_c2Face; //7 Probably Unused. Remove
 
     double                              m_transmissibility;
+
+    std::vector<cvf::Vec3d>             m_polygon;
+
  /*   enum NNCType
     {
         
@@ -64,12 +61,14 @@ class RigNNCData : public cvf::Object
 public:
     RigNNCData();
 
+    const std::vector<size_t>& findConnectionIndices(size_t globalCellIndex, cvf::StructGridInterface::FaceType face) const;
     void processConnections(const RigMainGrid& mainGrid);
   
     std::vector<RigConnection>&         connections()        { return m_connections; }
     const std::vector<RigConnection>&   connections() const  { return m_connections; };
 
 private:
-
+    typedef std::map<size_t, caf::FixedArray<std::vector<size_t>, 7 > > ConnectionSearchMap;
+    ConnectionSearchMap m_cellIdxToFaceToConnectionIdxMap;
     std::vector<RigConnection> m_connections; 
 };